Quantum chaotic tunneling in graphene systems with electron-electron interactions

An outstanding and fundamental problem in contemporary physics is to include and probe the many-body effect in the study of relativistic quantum manifestations of classical chaos. We address this problem using graphene systems described by the Hubbard Hamiltonian in the setting of resonant tunneling. Such a system consists of two symmetric potential wells separated by a potential barrier, and the geometric shape of the whole domain can be chosen to generate integrable or chaotic dynamics in the classical limit. Employing a standard mean-field approach to calculating a large number of eigenenergies and eigenstates, we uncover a class of localized states with near-zero tunneling in the integrable systems. These states are not the edge states typically seen in graphene systems, and as such they are the consequence of many-body interactions. The physical origin of the non-edge-state type of localized states can be understood by the one-dimensional relativistic quantum tunneling dynamics through the solutions of the Dirac equation with appropriate boundary conditions. We demonstrate that, when the geometry of the system is modified to one with chaos, the localized states are effectively removed, implying that in realistic situations where many-body interactions are present, classical chaos is capable of facilitating greatly quantum tunneling. This result, besides its fundamental importance, can be useful for the development of nanoscale devices such as graphene-based resonant-tunneling diodes.

Figure 1

Schematic illustration of four classes of geometrical domains for graphene billiards studied in this paper: (a) rectangle, (b) stadium, (c) bowtie, and (d) mushroom. The respective classical dynamics are integrable (a), chaotic with neutral periodic orbits (b), hyperbolic with all periodic orbits being unstable (c), and nonhyperbolic with mixed phase space (d). The thin gray region along a symmetric line represents the potential barrier.

Figure 2

Representative confined eigenstates associated with spin-up (red/solid curves) and spin-down (blue/dashed curves) electrons: (a) Probability density profile of the $5849\mathrm{th}$ eigenstate at $y=D/2$ without a potential barrier, where it extends in both potential wells. (b) Probability density profile of the $5860\mathrm{th}$ eigenstate at $y=D/2$ with a potential barrier of height $V_0=0.766t$ and width $w=2.5a$ at $x=0$ (represented by the gray rectangle). In (b) there is spin polarization, i.e., spin-up electrons reside in the right well and spin-down electrons reside in the left well. The corresponding eigenenergies are $E_{5849}\approx 0.6819t$ for (a) and $E_{5860}\approx 0.6924t$ for (b). The insets in both panels show the corresponding LDS patterns in the entire domain. Note that the wave functions have a zigzag appearance because they are plotted for both graphene sublattices, denoted by $A$ and $B$. (A plot of the wave function on one sublattice would appear more smooth.)

Figure 3

(a) Eigenenergy levels versus the eigenstate index. The red and blue bars with up and down arrows, respectively, represent a pair of polarized states ($n=5760$ and $m=5765$), where the corresponding eigenenergies are $E_n=0.6203t$, $E_{n+1}=0.6204t$, $E_m=0.6235t$, and $E_{m+1}=0.6237t$. The black bar without any arrow corresponds to a nonpolarized state. (b) and (c) Probability densities for the polarized state pairs ($n$, $n+1$) and ($m$, $m+1$), respectively. Note that the energy difference for each pair is quite small.

Figure 4

Tunneling rate $1/\Delta T$ versus $\Delta P$ for the rectangular graphene billiard for two cases where (a) there are electron-electron interactions as described by the Hubbard model and (b) there is only a single electron in the system.

Figure 5

Tunneling rates versus $\Delta P$ for chaotic geometries: (a) Stadium and (b) bowtie graphene billiards. The green dashed horizontal line indicates the separation of the two classes of spin-polarized (lower) and unpolarized (upper) states. The black solid up- and down-triangles correspond to the typical LDS patterns in each class for spin-up and spin-down states, respectively.

Figure 6

Tunneling rates versus $\Delta P$ for the nonhyperbolic mushroom billiard. The green dashed horizontal lines indicate the separation of the data points into three different classes. The middle and lower classes correspond to the chaotic components. The upper class is originated from the stem of the mushroom billiard. The black solid up- and down-triangles correspond to the LDS patterns for spin-up and spin-down states in the right-side panels, respectively.

Figure 7

(a) Illustration of a rectangular graphene domain of length $L$ and width $D$ with a potential barrier in the middle. The zigzag and armchair edges are along the $x$ and $y$ axes, respectively. The filled red (blue) circles at the right (left) zigzag boundaries denote the positive (negative) magnetic moments, the radii of which represent the strength of the spin density. (b) The effective potential profile ${\tilde{V}}_{\uparrow }\left(x\right)$ at the position $y=D/2$ for spin-up electrons, where the positive and negative potentials near the boundaries represent the strength of the respective spin density, which approach the value of $M_0$ ($-M_0$) at the left (right) zigzag edges in (a). The effective potential for spin-down electrons is reversed at the boundaries as compared with that for spin-up electrons.

Figure 8

(a) Illustration of the shift in the wave function caused by the effective potential $M\left(x\right)$. The black solid vertical lines are the actual domain boundaries. The gray dashed lines represent the effective boundaries. The black solid curves show the actual probability of an evenly confined spin-up state without the potential barrier, and the red dashed curves represent the effective wave function after the “shift.” (b) and (c) Antiphase and in-phase wave functions for sublattice $A$ and $B$, respectively, where red circles represent the simulation results from the mean-field Hubbard Hamiltonian and the gray solid and dashed curves are predictions from theory [Eqs. (9), (11), and (13)]. The boundary shift is $\delta =0.032a$. The wave function mode in the $y$ direction is chosen to be $n_y=16$ for (b) and $n_y=112$ for (c). The numerical results are from the (arbitrarily chosen) states $n=5547$ and 5548 for (b) and $n=5568$ and 5569 for (c).

Figure 9

Time evolutions of $P^L\left(t\right)$ for spin-up states corresponding to $n=5858$ (a), 5826 (b), and 5692 (c), and their projected states with the corresponding weighting coefficients. The green dashed-dotted curves are fitted from Eq. (15) (a) and (b) and Eq. (16) (c).