Superpersistent currents and whispering gallery modes in relativistic quantum chaotic systems

Persistent currents (PCs), one of the most intriguing manifestations of the Aharonov-Bohm (AB) effect, are known to vanish for Schrödinger particles in the presence of random scatterings, e.g., due to classical chaos. But would this still be the case for Dirac fermions? Addressing this question is of significant value due to the tremendous recent interest in two-dimensional Dirac materials. We investigate relativistic quantum AB rings threaded by a magnetic flux and find that PCs are extremely robust. Even for highly asymmetric rings that host fully developed classical chaos, the amplitudes of PCs are of the same order of magnitude as those for integrable rings, henceforth the term superpersistent currents (SPCs). A striking finding is that the SPCs can be attributed to a robust type of relativistic quantum states, i.e., Dirac whispering gallery modes (WGMs) that carry large angular momenta and travel along the boundaries. We propose an experimental scheme using topological insulators to observe and characterize Dirac WGMs and SPCs, and speculate that these features can potentially be the base for a new class of relativistic qubit systems. Our discovery of WGMs in relativistic quantum systems is remarkable because, although WGMs are common in photonic systems, they are relatively rare in electronic systems.

A remarkable phenomenon in the quantum world is persistent currents (PCs), permanent currents without any external source 1 , which are generated by the Aharonov-Bohm (AB) effect 2 in non-superconducting systems. PCs have been observed experimentally in metallic [3][4][5][6][7] and semiconductor [8][9][10] rings in the mesoscopic regime. Theoretical efforts have been focused on the effects of bulk disorders [11][12][13] , electron-electron interactions [14][15][16] , spin-orbital interactions 17,18 , and electromagnetic radiation 19,20 on PCs, typically studied in the diffusive regime using idealized circular-symmetric rings and cylinders. Rapid advances in nanotechnology have made it feasible to fabricate mesoscopic devices with mean free path larger than their sizes at sufficiently low temperatures (the ballistic transport regime) 21 . The AB system can thus be modeled as a quantum ballistic billiard in which the particles are scattered at the boundaries of the domain. As a result, the boundary shape becomes highly relevant. In experiments, uncontrollable boundary imperfections are inevitable 22,23 even when there are no bulk disorders. It is thus of interest to study the effects of boundaries, e.g., those that generate chaos in the classical limit, on PCs. In general, an asymmetric boundary destroys angular momentum conservation and introduces irregular scattering. Theoretical [24][25][26][27][28][29] and experimental 22 studies have shown that, similar to the effects of bulk disorder, symmetry breaking of the boundary can result in opening of gaps at the degeneracy points of the energy levels, leading to level repulsion, a typical manifestation of classical chaos. Energy gap opening can diminish AB oscillations through pinning of the corresponding states, leading to vanishing PCs. Since fully chaotic domains are associated with a strong degree of symmetry breaking, PCs are not expected to arise 30,31 . In nonrelativistic quantum systems governed by the Schrödinger equation, PCs are thus fragile.
In this paper, we address the question of whether, in relativistic quantum systems, PCs can arise and sustain in the presence of symmetry-breaking perturbations. Besides the importance of this question to fundamental physics, our work was motivated by the tremendous recent research on two-dimensional Dirac materials 32 such as graphene [33][34][35][36][37][38][39] , topological insulators 40 , molybdenum disulfide (MoS 2 ) 41,42 , HITP [Ni 3 (HITP) 2 ] 43 , and topological Dirac semimetals 44,45 . The physics of these materials is governed by the Dirac equation. This is thus interest in investigating relativistic PCs in Dirac fermion systems [46][47][48][49][50][51][52][53][54][55][56][57] . Existing theoretical works on relativistic PCs, however, assumed idealized circular-symmetric rings in the ballistic limit. Whether AB oscillations and consequently relativistic PCs can exist in asymmetric rings that exhibit chaos in the classical limit is a fundamental question. Our finding is that, even in the presence of significant boundary deformations so that the classical dynamics becomes fully chaotic, robust PCs can occur in relativistic quantum, Dirac fermion systems, henceforth the term superpersistent currents (SPCs). A more striking finding is that SPCs are generated by localized states at the domain boundaries, which are effectively chaotic Dirac whispering gallery modes (WGMs) that carry larger angular momenta. WhileWGMs are common in photonic systems [58][59][60][61] , its emergence in electronic systems 62 , especially in relativistic quantum systems, is rare and surprising. We develop a physical understanding of the counterintuitive phenomenon of SPCs by analytically exploiting the properties of the spinor wavefunctions in an idealized relativistic quantum system.
The significance of our results lies in the perspective of observing SPCs in the presence of strong random scattering, implying that they may occur in systems of size beyond the mesoscopic limit. This can potentially be a relativistic quantum phenomenon occurring at relatively large scales. There can be significant applications of SPCs in quantum information processing, for which we speculate on a scheme of Dirac WGM-based qubit. We note that, previous experimental studies of AB oscillations in graphene 63 and topological insulators 64,65 make it possible to experimentally test the phenomena of DiracWGMs and SPCs that we predict in this paper. To motivate experimental verification, we propose a feasible scheme using threedimensional topological insulators.
Theoretical model of the relativistic AB system. Consider a single massless spin-half particle of charge 2q confined by hard walls in a domain B with a ring topology in the plane r 5 (x, y), as shown in Fig. 1. Applying a single line of magnetic flux (AB flux) W at its origin and utilizing an infinite mass term outside the domain to model the confinement, we obtain the following Hamiltonian in the position representation:Ĥ~vŝ where v is the Fermi velocity,ŝ~ŝ x ,ŝ y À Á andŝ z are the Pauli matrices, and M 5 0 in the ring domain but M 5 ' otherwise (for hard-wall confinement -used previously, e.g., in the study of graphene rings 66 , graphene quantum dots 67 , and topological insulator quantum dots 68 ). The vector potential A(r) can be chosen as any vector field satisfying =|A r ð Þ~êaW 0 d r ð Þ, whereê is the unit vector perpendicular to B and a ; W/W 0 is the dimensionless quantum flux parameter, with W 0~2 p h=q being the flux quantum.
The HamiltonianĤ acting on the two-component spinor wavefunction y(r) 5 [y 1 , y 2 ] T has eigenvalue E: where D 5 = 1 ia is a compact notation for the covariant derivative with a ; 2pA/W 0 . In Eq. (2), Ms z represents the mass potential that takes into account the Klein tunneling effect, which can confine massless Dirac fermions in a finite domain. Such a confinement itself will break the time-reversal symmetry (T-breaking) in absence of any external magnetic field. However, the T-breaking due to mass potential confinement is different from that due to magnetic field in a classical picture in which no Lorenz force acts on the particle so that the geodesic motions are still characterized by straight lines within the confined domain. In relativistic quantum theory of electrons, magnetic field is introduced through the form of minimal coupling in the corresponding vector potential that is different from the mass term. As a result, in the conventional (3 1 1)-dimensional spacetime, the mass itself cannot break the time-reversal symmetry. However, in (2 1 1)-dimensional systems studied in this paper (as for twodimensional Dirac materials such as graphene and topological insulators), the mass term will induce chiral anomaly and hence T-breaking. Thus the T-breaking caused by mass confinement has different features from that due to the magnetic field 69,70 .
To better understand the physical origin of the mass confinement term Ms z , we consider a single particle in absence of magnetic field and compare the following two situations: (a) the particle is in a Dirac ring with hard-wall confinement of mass potential and (b) the particle is in a Schrödinger ring with the conventional, electrical potential (hard-wall) confinement. We assume that the classical orbits are identical for both cases. Apparently, in case (a), the T-breaking is intrinsic but there is no T-breaking in case (b). In fact, as indicated by Sir Berry 71 , the semiclassical origin of T-breaking induced by mass confinement is quite intriguing. One aim of our work is to uncover new and interesting phenomena caused by the mass confinement Ms z in (2 1 1)-dimensional Dirac systems. As we argue and demonstrate later, this type of confinement can be experimentally realized by exploiting the surface states of three-dimensional topological insulators, where the mass term in the Dirac equation is originated from a Zeeman term induced by local exchange coupling.
Some basic properties of Eq. (2) are the following. Firstly, the confinement condition of imposing infinite mass outside B naturally takes into account the Klein paradox for relativistic quantum particles and thus guarantees that our study is conducted in the singleparticle framework, which is relevant to the intrinsic physics of a single Dirac cone in graphene or topological insulator. Secondly, both reduced spatial dimension together with mass confinement and applied magnetic flux can break the time-reversal symmetry of H:T,Ĥ Â Ã =0 if M ? 0 or a ? 0, whereT~is yK andK denotes complex conjugate. Thirdly, for M 5 0 and A 5 0 in Eq. (2) (i.e. free massless particle), there exist plane wave solutions whose positive energy part has the following form: where k is a wave vector that makes an angle h with the x axis. Fourthly, by using the Dirac equation iL t y~Ĥy and defining r 5 y { y as the local probability density, we have the following continuity . It is therefore natural to define vŝ as the local current operator, so the local current density in state y(r) is given by j r ð Þ: To obtain solutions of Eq. (2), a proper treatment of the boundary condition is necessary. As done in previous works, we use the infinite where sgn(?) stands for the signum function and h(s) denotes the angle made by the outward unit normal n with the x axis at an arbitrary boundary point s, as shown in Fig. 1. Substituting Eq. (4) into the current density formula, one can show that the boundary current j(s) 5 2vjy 1 j 2 (2 sin h, cos h) is polarized along the boundary: clockwise and counterclockwise for the inner and outer edges, respectively. It is remarkable that this polarized property is independent of the shape of the confinement potential M(r) and is thus topologically protected from irregular boundary scatterings, even though the magnitude of the edge current can be affected. An analysis of the general properties of the a (magnetic flux) dependent relativistic quantum spectrum {E j (a)}, as determined by Eq. (2) under the boundary condition Eq. (4), reveals that the first ''Brillouin zone'' is given by 21/2 # a # 1/2 (Supplementary Note 1). To calculate a large number of relativistic eigenvalues and eigenstates with high accuracy, we use the conformal-mapping method [75][76][77] (Supplementary Note 2).
Whispering gallery modes and superpersistent currents. To demonstrate our findings, we deform a circular ring domain j 5 0.5 # r # 1, using the mapping w(z) 5 h[z 1 0.05az 2 1 0.18a exp(iv)z 5 ], where v 5 p/2, a g [0, 1] is the deformation parameter that controls the classical dynamics. Specifically, when increasing a from zero to unity the deformed ring will undergo a transition from being regular to mixed and finally to being fully chaotic. The normalization coefficient tees that the domain area is invariant for arbitrary values of the deformation parameters {a, v, j}. Four representative domains are shown in the top row in Fig. 2 where, classically, the left most domain is integrable, the right most domain is fully chaotic, and the two middle domains have mixed phase space. The middle and bottom rows of Fig. 2 show the lowest 10 energy levels as functions of the quantum flux parameter a, i.e., energy-flux dispersions, for Schrödinger and Dirac particles, respectively. We see that E j (a) 5 E j (2a) holds for the Schrödinger particle, but for the Dirac fermion, the symmetry is broken: E j (a) ? E(2a). However, for both nonrelativistic and relativistic spectra, we have E j (a) 5 E j (a 1 1). For the circular-symmetric ring (a 5 0), AB oscillations in the energy levels have the period W 0 (i.e., a 5 1) and there are level crossings. Making the domain less symmetric by tuning up the value of the deformation parameter a leads to classical mixed phase space (regular and chaotic), and eventually to full chaos (a 5 1). We see that, for the Schrödinger particle, emergence of a chaotic component in the classical space leads to opening of energy gaps, generating level repulsion and flattening the AB oscillations associated with the corresponding energy levels. Surprisingly, for the Dirac fermion, the AB oscillations are much more robust against asymmetric deformations. In particular, for the fully chaotic case, AB oscillations for the Schrödinger particle disappear almost entirely while those for the Dirac fermion persist with amplitudes of the same order of magnitudes as the integrable case.
We now present evidence of Dirac WGMs for the case of fully chaotic AB ring domain. By examining the eigenstates, we note that, for low energy levels, the Schrödinger particle is strongly localized throughout the domain, as shown in Figs. 3(a-c), leading to a flat energy-flux dispersion. However, the Dirac fermion typically travels around the ring's boundaries, forming relativistic WGMs that persist under irregular boundary scattering due to chaos and are magneticflux dependent, as shown in Fig. 3(d-f). Conventional wisdom for Schrödinger particle stipulates that asymmetry in the domain geometry can mix/couple well-defined angular momentum states, opening energy gaps and leading to localization of lower states in the entire domain region, so AB oscillations would vanish, as demonstrated both theoretically and experimentally 22,24 . However, for Dirac fermion, this picture breaks down -there are robust AB oscillations even in the fully chaotic domain and the particle tends to execute motions corresponding to WGMs.
The robust AB oscillations in chaotic Dirac rings lead to SPCs. The total persistent current can be calculated at zero temperature through 1,66 I a ð Þ~{ X j LE j La À Á , where the sum runs over all occupied states with E j . 0. Due to periodicity in the energy: E j (a) 5 E j (a 1 1), the current is also periodic in a with the fundamental period a 5 1. Figure 4 shows, for the nonrelativistic (top panels) and relativistic (bottom panels) cases, PCs resulted from the lowest three states (including spin) in regular (left column) and chaotic (right column) rings. We see that, for classical integrable dynamics, PC oscillations display a common sawtooth form. However, at zero flux, PC is zero for the nonrelativistic case [ Fig. 4 (a)], while it has a finite value for the relativistic case due to breaking of the time-reversal symmetry. In the chaotic case, the oscillations become smooth due to level repulsion in the corresponding energyflux dispersion pattern. As a result, PCs carried by the Schrödinger particle practically vanish as compared with the integrable case but, strikingly, the Dirac fermion still carries a persistent current with amplitude of the same order of magnitude as that for the integrable case -SPCs. Intuitively, SPCs carried by the Dirac fermion as an ''exceptional'' magnetic response are associated with the chaotic Dirac WGMs exemplified in Fig. 3.
Origin of WGMs and SPCs. The origin of the ''exceptional'' magnetic response of the chaotic Dirac fermion can be understood through the behavior of the current carried by the particle at the boundary interface. We have developed an analytic understanding to predict the occurrence of Dirac WGMs and, consequently, SPCs. In particular, we consider the following problem: a plane wave incident obliquely on a straight potential jump M(x, y) given by as shown in Fig. 5. Without loss of generality, we let the incident wave y i be described by the wave vector k 0 5 (k cos h 0 , k sin h 0 ), the reflected wave y r by k 1 5 (k cos h 1 , k sin h 1 , and the transmitted wave y t by u 5 (iq, K), where h 1 5 p 2 h 0 and K ; k sin h 0 . We focus on the situation where the energy of the incident wave satisfies E , V 0 , which corresponds to the total reflection case. Referring to Fig. 5, we have that the wave in region I (x , 0) can be written as and the wave in region II is where the coefficients R and T are to be determined by matching the waves at x 5 0. In the following, we treat the Schrödinger scalar wave and Dirac spinor wave separately.   Schrödinger scalar plane wave. For the nonrelativistic quantum case as shown in Fig. 5(b), we have, in region I (x , 0), and in region II (x . 0), where q and K are related to each other through with m denoting the mass of the particle. Matching the waves and their derivatives at x 5 0, we obtain where the parameter b is defined through Given the wave function Y S (r), the associated probability current density is where c.c. denotes complex conjugate. We therefore obtain, in region I, In region I, the y component of the probability current is the sum of two terms: (1) the term resulting from the sum of the currents associated with the incident and reflected waves, and (2) the term containing the factor cos(2kx cos h 0 1 2b) that accounts for the interference between the incident and reflected waves. In region II, the probability current is also parallel to the y-axis, and it decays exponentially as an evanescent wave. Note that ± J II S À Á y ?0 as q R ' (i.e., V 0 R ').
Planar Dirac spinor wave. The relativistic case is shown in Fig. 5(c). We proceed in the same manner as for the nonrelativistic case. Expressing the wave in terms of massless spinor planar waves that are solutions of the Dirac equation, we have, in region I (x , 0), where l 1~ffi and q~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with v being the Fermi velocity.
Matching boundary conditions at x 5 0, we obtain R~i where the parameters c and l are defined through and For the spinor wave Y D (r) describing a Dirac fermion in two dimensions, the corresponding probability current density is whereŝ~s x ,s y À Á and Y D 5 [y 1 , y 2 ] T . We have, in region I,  Fig. 6. We see that the nonrelativistic transverse current density (J S ) y is antisymmetric with respect to h 0 , leading to zero contribution to J y S , while the relativistic transverse current density (J D ) y is a nonnegative monotonic function of h 0 so that there exists a finite transverse current even for a fully chaotic ring when all  possible incident angles are taken into account. Such a finite averaged transverse current density J y D with magnitude one half the maximum makes it possible to form chaotic Dirac whispering gallery modes that carry considerable directional currents. This is the fundamental mechanism for the phenomenon of SPCs.
Note that, although Eqs. (25) and (26) appear different in form, at the interfact (x 5 0) they give exactly the same current density (see Supplementary Note 3). For V 0 R ', the hard-wall boundary condition is restored (the infinity-mass boundary condition). Physically the counterintuitive phenomenon can be understood, as follows. The incoming wave from the metal region (mass term M 5 0) is spin polarized along its momentum (current) direction. After entering the insulator region, a finite mass term acting on s z will change the direction of the spin and hence affect the current via the spin-momentum locking term k ? s.
The relevant Hall-like phenomenon associated with the T-breaking mass potential has been uncovered in a very recent work 78 . In addition, we note that the results do not depend on the special form of the wave function in region II (although such a symmetric form is convenient for analysis). In fact, we obtain the same results 38 when choosing the wavefunction to have the form , [1, C(E)] T .
Experimental scheme. A possible experimental scheme to observe and characterize Dirac WGMs and SPCs is, as follows. A 3D topological insulator supports a (2D) gapless state on its surface, with low-energy excitations described by the massless Dirac Hamiltonian 40,79Ĥ surface~{ i hv Fŝ : =, whereŝ characterizes the spin. The surface electronic structure is similar to that of graphene, except that there is only a single Dirac point. Different from graphene, the Dirac surface states of a topological insulator are associated with strong spin-orbit interactions. In spintronics applications of topological insulators, it is desirable to introduce a gap into the surface states. This can be done by breaking the timereversal symmetry using a ferromagnet insulator (FMI) deposited on the top of a topological insulator 68 . The exchange coupling induced due to proximity to the ferromagnet insulator will give rise to a local exchange field that lifts the Kramers degeneracy at the surface Dirac point and introduces a mass term into the Dirac Hamiltonian. Thus, generally, we haveĤ~vŝ :p zeA zM r ð Þŝ z zc z Bŝ z , where the vector potential A accounts for the effect of the external magnetic field B~=|A~Bẑ, with an additional Zeeman splitting correction in the last term. The controllable mass term Mŝ z , responsible for the local exchange coupling with a FMI cap layer, makes 3D topological insulators a potential experiment platform for observing and characterizing Dirac WGMs and SPCs.
Robust relativistic qubit based on WGMs. We speculate on a potential application of DiracWGMs in quantum information technology. Similar to the proposal of qubit in a two-dimensional topological insulator based on the two different sets of helical edge states localized at the boundaries 56 and the idea of charge qubit in a double quantum dot 80 , we present our qubit based on the chaotic Dirac WGMs guided by the inner and outer surfaces in opposite directions. Such an edge degree of freedom can be used to form a two-state system, denoted by the states jonae and joffae for the outer and inner chaotic DiracWGMs, respectively. The set {jonae,joffae} thus constitutes a complete diabatic basis. Generally, these two states correspond to two energy levels {E L (a), E U (a)} with a rather larger difference (mismatch) and can be coupled and superpositioned for different values of the magnetic flux. As a result, two different levels of the system arise, say {E L (a9), E U (a9)}. The corresponding instantaneous eigenstates {jLae,jUae} constitute an adiabatic basis. An effective a-dependent Hamiltonian describing the flux-tunable qubit can then be written in the adiabatic basis asĤ a qubi a ð Þ~E L a ð Þ L j i L h jz E U a ð Þ U j i U h j, providing a base for exploiting the Dirac WGMs as a qubit system.

Conclusions.
We formulate a relativistic version of AB chaotic billiards to study PCs in Dirac rings. We find that, in contrast to the nonrelativistic quantum counterpart where PCs vanish for chaotic rings, the currents continue to exist in the relativistic chaotic AB rings and, in this sense, they are superpersistent. We demonstrate that SPCs are a consequence of Dirac WGMs, and we develop an analytic understanding of their emergence in relativistic quantumsystems. We also propose that, experimentally, chaotic rings patterned by magnetic domain heterostructures deposited on the surface of a 3D topological insulator can be a feasible scheme to observe and characterize chaotic Dirac WGMs and SPCs. The coexistence of inner and outer chaotic Dirac WGMs naturally forms a flux-tunable two-level system. To investigate the magnetic response of chaotic Dirac fermions is not only fundamental to the emerging field of relativistic quantum chaos, but also relevant to device applications based on Dirac materials.