Anomalous Electrodynamics and Quantum Geometry in a Graphene Heterobilayer
AAnomalous Electrodynamics and Quantum Geometry in a Graphene Heterobilayer
Abigail Timmel and E. J. Mele ∗ Department of Physics and AstronomyUniversity of PennsylvaniaPhiladelphia, PA 19104 (Dated: February 24, 2021)Graphene heterobilayers with layer antisymmetric strains are studied using the Dirac-Harpermodel which describes a pair of single layer Dirac Hamiltonians coupled by a matrix-valued one-dimensional moir´e-periodic interlayer potential. We find that this model hosts low energy, nearly dis-persionless bands near charge neutrality that support anomalous polarizations of multipole distribu-tions of the charge density. These are studied by formulating symplectic two-forms that encode field-induced multipole dynamics allowed in a chiral medium with time reversal symmetry. The modelidentifies reciprocity relations between the responses to layer-symmetric and layer-antisymmetricin-plane electric fields and reveals momentum-space quantum oscillations that can quantify thestructure of its laterally modulated band inversions.
Response functions describing the electrodynamicproperties of electrons in crystals can contain anoma-lous terms arising from momentum dependence of inter-nal spin and orbital attributes in Bloch bands. Therehas been interest in interpreting these as geometric phe-nomena associated with a quantum metric [1, 2]. Thisapproach generalizes to driving fields that vary in space.When a smooth spatial variation lowers the symmetry ofthe system, it can activate responses that are symmetry-forbidden for strictly spatially uniform fields. These spa-tially dispersive responses can sometimes be expressedas responses of a system to the field gradients . Arti-ficial structures fabricated by stacking atomically thinlayers are of great current interest in this regard, andin these systems the field gradient is replaced by a dis-cretized version that represents the response to drivingfields that can be distinguished on the individual lay-ers [3–5]. In twisted graphene bilayers these discretizedcoupled layer responses can greatly exceed spatially dis-persive responses found in conventional chiral materials[6].In this Letter we formulate electrodynamic responsesof the Dirac Harper (DH) model for a graphene moir´e het-erobilayer in terms of a symplectic two-form on the pro-jective Hilbert space [2, 7]. This generalizes the anoma-lous velocity associated with a Berry curvature to a dif-ferent family of anomalous charge displacements that canoccur in a time-reversal invariant medium. We illustratethis effect in the DH model which presents a one dimen-sional version of the moir´e heterostructures encounteredin twisted bilayer graphenes and in other layered two di-mensional van der Waals heterostructures. Instead of atwist, in the DH model the interlayer atomic registry ismodulated by introducing a layer-antisymmetric latticestrains (Figure 1a). Because of the lattice periodicityof the individual layers, the interlayer registry varies pe-riodically on a long moir´e length scale while retaininga short range periodicity in the orthogonal coordinate.In this way the problem is separated into a family of one dimensional problems, with each sector indexed bya conserved crystal momentum k y in Figure 1b and sup-porting a one-dimensional spinorial variant of the Harperproblem in the modulation ( x ) direction. The DH Hamil-tonian introduced in [8] gives a spinorial Harper equation t ( k ) ψ n +1 + t † ( k ) ψ n − + ( t + ν ( x )) ψ n = (cid:15) ψ n (1)where ψ are four component fields with the ordering ofamplitudes on two sublattices and two layers: ψ T =( a , b , a , b ). Here t ( k ) are block diagonal 4 × ν ( x ) = ν ( x + L ) is the spatially modulated interlayer cou-pling with moir´e period L [9]. Defining the Pauli matrices σ i and τ i acting on sublattice index and layer index re-spectively, the products σ x τ = γ and σ y τ = γ forma basis for t . The symmetries of this model depend on theform of the lattice strain and the orientation of the modu-lation direction. An important limiting case studied hereoccurs for an antisymmetric shear strain along the arm-chair (AC) direction which produces a chiral structurewith C x , C y , and C z symmetries shown in Figure 1a.The k y -projected Hamiltonians individually break timereversal symmetry T though this is restored by the k y -space integrated response with a T -invariant population.As discussed in [8] the modulated interlayer couplingmatrices can be represented by a set of Dirac matri-ces γ = σ τ x , γ = σ x τ x , and γ = σ y τ y weightedby moir´e periodic (real) amplitudes. Sign changes ofthese amplitudes can spatially confine electronic statesto selected regions of the modulated structure. Modestrapped in the AA regions are confined to a fixed re-gion along the modulated coordinate for all q r within theband-inverted region, and the density of another suchconfined state for the AC moir´e pattern is plotted as afunction of q r = ( t (cid:48) a/ (cid:126) v F ) k y , (where t (cid:48) is the interlayercoupling strength and v F = √ ta/ (cid:126) is the Fermi veloc-ity for lattice coupling strength t on a honeycomb latticewith lattice constant a ) in Figure 2a. These modes coun- a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b terpropagate from the AB and BA regions towards the AA region as k y approaches a gap closure at q r = 1. Inthe semiclassical theory this counterpropagation could bedriven by an electric field (cid:126) ˙ k y = − eE y , but it transportsno charge along x . The spectral flow illustrated in Fig-ure 2a can nonetheless be associated with transport ofhigher order multipole densities. FIG. 1. (a) Shear-strained moir´e structure. (b) Low energyband structure for the shear-strained moir´e at k x = 0. Allbands are two-fold degenerate. Inside the gap, a set of eightweakly dispersing bands appear near zero energy. The outer-most four bands have AA -peaked modes, the innermost fourhave AB/BA -peaked bands.
Multipole distributions that are antisymmetric in x canassign a net polarization to the counterpropagating peaksby tagging the two peaks with opposite signs. The low-est allowed multipole moment that respects C x , C y ,and C z symmetry of the shear-strained AC moir´e is anoctopole. Thus a spatially antisymmetric density can beproduced by charge distributions that locally have the ro-tational symmetry of a quadrupole in the yz plane. Suchan octopole is explicitly present in the DH Hamiltonian:the mass term ˆ γ sin( Gx ) combines the quadrupole op-erator ˆ γ with an antisymmetric function that distin-guishes the AB and BA stacked regions of the struc-ture. The sublattice correlations measured by ˆ γ encodea property of the wavefunction which switches sign inthese two regions and tags the charge peaks with oppositesign. Similarly, the quadrupole density locally measuredby ˆ γ = σ z τ z describes a local layer and sublattice polar-ization that is also spatially antisymmetric (Figure 2).The spectral flow of a multipole density can be associ-ated with a generalization of the symplectic 2-form that FIG. 2. (a) Anomalous drift of the charge peaks in
AB/BA -peaked bands as a function of q r . The gap closure occursat q r = 1. (b) Multipole densities ρ and ρ at q r = 0 . defines the Berry curvature: ω m = (cid:88) i,α,β m αβ d u ∗ iα ∧ d u iβ Here u iα denotes the wavefunction component for the α th sublattice of the i th cell, and m αβ is a local mul-tipole weighting the sublattice degrees of freedom. For m αβ = I , we recover the usual Berry curvature responsi-ble for the anomalous Hall conductance. More generally,as long as m αβ has a trivial nullspace, then ω m is alsoa nondegenerate (symplectic) 2-form. By the DarbouxTheorem, this can be re-expressed in the familiar formof the Berry curvature (cid:80) i,α d u (cid:48)∗ iα ∧ d u (cid:48) iα via a coordinatetransformation, but for ease of interpretation we shallcontinue to work in the ( a , b , a , b ) basis. When ω m in this basis is applied to tangent vectors ∂ λ ν and ∂ λ µ at | u (cid:105) , the resulting curvature is weighted by the multipole M ≡ m αβ ⊗ I N . ω m ( ∂∂λ µ , ∂∂λ ν ) = (cid:104) ∂u∂λ µ | M | ∂u∂λ ν (cid:105) − (cid:104) ∂u∂λ ν | M | ∂u∂λ µ (cid:105) These generalized curvatures encode local multipoledensities via m αβ . To choose parameters λ , we observethat the one-dimensional moir´e structure is naturally de-scribed using two extended coordinates: a momentumcoordinate in the y direction and a spatial coordinate inthe x direction. The two dual coordinates k x and y pro-vide little additional information about the system due toits large aspect ratio. For definiteness we choose k x = 0promoting the model with symmetries that can be onlyslightly broken for this system. For the extended spatialcoordinate, we introduce the spatially varying mass fieldin the Hamiltonian. The masses appear in a sum overthe first star of G -vectors, ν ( x ) = (cid:88) G,α γ α ν G,α e iG ( x + x α ) where ν G,α = ν ∗− G,α and α indexes three symmetry-allowed Dirac mass matrices γ , γ , and γ . The x α are spatial parameters that allow us to vary the originfor each mass field individually.With our Hamiltonian defined by parameters k y and x α , we study the curvature: K ( M ) α = (cid:104) ∂u∂k y | M | ∂u∂x α (cid:105) − c.c. (2)Note that k y paired with any x α parametrizes a two-dimensional subspace of the projective Hilbert space,but the ground state lies in the one-dimensional sec-tor at x α = 0 where the mass fields are pinned bythe moir´e structure. The two-dimensional tangent space parametrized by ∂ x α , ∂ k y nevertheless encodes an observ-able of these x α = 0 states which will be described later.This curvature can be calculated in a two-fold degenerateband subspace using a modified Wilson loop [9].Now we study the symmetry of the curvature under thethree C rotations and time-reversal T . We assume that M is Hermitian so that ω m is imaginary. The tangentvectors ∂ x α and ∂ k y define a quadrupole in the xy plane. K ( M ) must have the symmetry of an octopole to respect C symmetry, and this is satisfied if M has the symmetryof a z -directed dipole. Here we use m = τ z as a physicallyintuitive candidate. Under T , a sign change of ∂ k y alongwith complex conjugation brings the imaginary curvatureback to itself, so K ( M ) does not require T -breaking todescribe an net response.Among the symmetries of the system, C x T is of par-ticularly interest because it is local in both extended co-ordinates x and k y . The symmetry of m αβ under this op-eration therefore determines whether the curvature K ( M ) or the multipole density ρ m vanish locally in momentumspace and position space respectively. For the curvature,one finds [9] (cid:104) ∂u∂k y | M | ∂u∂x α (cid:105) = (cid:104) ∂u∂k y | M (cid:48) | ∂u∂x α (cid:105) ∗ for M (cid:48) = ( C † x M C x ) ∗ . When M (cid:48) = M , the curvature isequal to its complex conjugate and therefore it vanisheseverywhere. Indeed this is the case for the usual Berrycurvature with M = I , but the generalized extensionscan incorporate multipoles such as τ z where M (cid:48) = − M ,allowing this response to be nonvanishing.For the multipole density ρ m , the same manipulationgives ρ m = ρ ∗ m (cid:48) where m (cid:48) is the x -projected version of M (cid:48) .Since ρ m is real, this quantity vanishes when m (cid:48) = − m .Thus we find that the symmetry of m under C x T di-vides multipole densities into two classes: those that canhave nonvanishing density but vanishing curvature and those that have vanishing density but can have a non-vanishing curvature. The anomalous response measuredby Equation 2 therefore is physically distinct from theconstruction using the antisymmetrically-tagged driftingpeaks constructed in Figure 2. Dynamical effects occurin quantities with no accumulation of a multipole den-sity within the bulk, and quantities with accumulationcan only display a static polarization.The triad of operators ∂ x α , ∂ k y , τ z can now be iden-tified with two response functions. The tangent vector ∂ k y couples the system to a perturbing electric field inthe y -direction, and the tangent vector ∂ x α arises fromthe observable (cid:126) ˙ k ( α ) x = [ ∂ x α , H ] which describes an in-duced force in the x -direction for the α -th mass channel.To introduce the third member of the triad, τ z , we canconsider either a layer-antisymmetric electric field or alayer-antisymmetric response. For the former, τ z entersin the driving term Eτ z ˆ y , and the force is given by f yα = ieE (cid:88) m (cid:54) = n (cid:104) n | [ ∂∂x α , H ] | m (cid:105)(cid:104) m | τ z | ∂n∂k y (cid:105) (cid:15) m − (cid:15) n − c.c. (3)where the first subscript labels the force by an inducedlattice effect, to be described later. For the latter, τ z weights the observable i∂ x α with opposite sign in eachlayer: f xα = ieE (cid:88) m (cid:54) = n (cid:104) n | [ τ z ∂∂x α , H ] | m (cid:105)(cid:104) m | ∂n∂k y (cid:105) (cid:15) m − (cid:15) n − c.c. (4)Both of these can be written as curvatures K ( τ z ) α in Equa-tion 2, and thus they are equivalent [9]. Numerical cal-culation for each mass channel α is plotted as a functionof q r in Figure 3.Equations 2, 3 and 4 are the main results of this workand generalize the constitutive relations for a mediumwith natural optical activity [3, 10–12]. The differentmass channels each target the region of charge densityconfined by the associated mass inversion. This is shownconcretely in Figure 3c, where the AA peaks definedby γ inversion are shifted by x (i) and the AB/BA peaks defined by γ inversion are shifted with x (ii). x displaces little charge (iii) because γ does not con-fine charge via a sign change [8]. The equivalence ofEquations 3 and 4 identifies a reciprocity relation statingthat a force induced by a multipole-weighted field is equalin magnitude to the multipole-weighted force induced bya uniform field. For our example of a z -directed dipole τ z ,these responses are driven by a layer-antisymmetric field (cid:126)E (cid:28) and a layer symmetric field (cid:126)E ⇒ respectively. For theformer (Equation 3), summing over mass channels givesa force depinning charge in the negative x direction asshown in Figure 3c iv. This can produce a counterforceon the lattice that acts to bring the moir´e back in phasewith the charge by sliding the layers relative to each other FIG. 3. (a) Reciprocal systems described by i. Equation 3and ii. Equation 4. (b) Curvatures K ( τ z ) α for each of the masschannels. (c) Change of the charge density from i. x , pri-marily affecting the AA charge peak, ii. x , primarily affect-ing the AB/BA peaks, iii. x , causing little displacement,and iv. all the mass terms together, inducing a constant shift. along y , transporting the moir´e pattern in x [13]. Theantisymmetrically weighted response in Equation 4 cor-responds to a similar depinning of charge, but weightedin opposite directions in the two layers. The analogouslattice counterforce would act to bring the charge backin phase across layers by sliding the layers in oppositedirections along x . The subscripts on Equations 3 and4 should now be clear: f y induces a y -directed shear be-tween layers and f x induces an x -directed shear. For ageneral y -directed electric field (cid:126)E = (cid:126)E (cid:28) + (cid:126)E ⇒ , the in-duced layer shear then acts at a general angle φ in the xy plane, which by reciprocity of the curvature is givenby φ = tan − ( f y /f x ) = tan − ( E (cid:28) /E ⇒ ).All three of the curvatures shown in Figure 3 oscillateas a function of k y . This behavior manifests oscillationsin the energies of the AB/BA bands, partially shownin Figure 1b, which are encoded in the response func-tion via the velocity-terms in the Kubo formula. The ob-served oscillations have a smoothly evolving period whichis nearly constant over the width of the band-invertedregion | q r | < γ potential, shownin Figure 4. A constant γ potential displaces the two hy- bridized Dirac cones in energy so that they intersect on aline degeneracy to form a low energy critical ring. Underzone folding from the moir´e structure, this would pro-duce a series of zero-energy bands crossing whose spacingevolves independently for bands originating on the ± q x sides of the ring and rapidly decreases toward an accu-mulation point at the edge of the critical ring q r = 1 [9].By contrast, a spatially varying γ potential produces in-stead a regular, slowly evolving pattern of band crossingsas seen in Figure 5a. This can be understood by squaringthe low-energy Dirac equation, (cid:0) i (cid:126) v F γ ∂∂x + (cid:126) v F k y γ + ν ( x ) γ (cid:1) ψ = (cid:15)ψ to form the effective potential problem ∂ ∂x ψ = (cid:16) − ν ( x ) (cid:126) v F + k y (cid:17) ψ ≡ V eff ψ (5)The tight-binding solutions shown in Figure 4 are foundto oscillate in the region where V eff is negative and de-cay in regions where it is positive. Wavelengths of theoscillating part of the wavefunctions match the value of √ V eff averaged over each period. Consecutive zero en-ergy band crossings differ by a single oscillation in thewavefunctions confined within the AA basin. Thus theregular band spacing of a varying γ potential arises froma resonance condition for the zero energy solutions withinthis effective potential. FIG. 4. (a) Zero energy crossings of the tight-binding spec-trum for the ν ( x ) potential. Tight-binding wavefunctionsplotted against the effective potential for (b) k y = .
08 and (c) k y = . These observations carry over to the full DH Hamilto-nian. Incorporating a term proportional to γ sin( Gx )opens a gap around the isolating the lowest energy bandmanifold without changing the spacing of these zero en-ergy crossings [9]. Away from the SP region, γ de-creases to zero. Therefore, oscillatory behavior in theenergies is inherited from the k y spacing of zero en-ergy crossings intrinsic to the γ part of the potential.Although the full moir´e AB/BA wavefunctions decaywithin the positive effective potential region, the sensi-tivity of the energies to these special values of k y sug-gests that the wavefunctions retain a small oscillatorypart stretching over the positive potential region flankedby the AB/BA charge peaks. Interestingly these oscil-lations would be absent from chiral models for twistedbilayer graphene which exclude this coupling to enforcea chiral symmetry [14]. Here we see that they controla physically measurable response and measurement ofthese oscillations can be used to determine the scale ofthe breaking of the chiral symmetry.The sympletic forms developed here can also be ap-plied to two dimensional twisted bilayers. In this casethe DH theory generalizes to a two dimensional theoryin a moir´e supercell and the forces driving the multipoledistributions analogous to Equations 3 and 4 are simi-larly promoted to two dimensional vectors obtained byintegrating over a population in the 2D folded Brillouinzone. The simplest case couples the driving fields to thetotal charge density recovering the charge pumping phe-nomena proposed for twisted graphene bilayers [13, 15].The generalized Berry curvature provides a unified for-mulation that associates this electromechanical responsewith the more general problem of manipulating multipoledistributions using applied fields. The generalized con-stitutive relations for electric field-driven multipole re-sponses also provides a geometric formulation of the nat-ural optical activity of twisted graphene bilayers studiedexperimentally [6] and theoretically [3–5] and to chiralWeyl semimetals [16].This work was supported by the Department of Energyunder grant DE-FG02-ER45118. ∗ [email protected][1] R. Resta “The Insulating State of Matter: A GeometricalTheory.” Eur. Phys. J. B , (2011).[2] T. Neupert, C. Chamon, and C. Mudry “Measuring theQuantum Geometry of Bloch Bands with Current Noise.”Phys. Rev. B , 245103 (2013).[3] E. Morell, L. Chico, and L. Brey, “Twisting Diracfermions: circular dicroism in bilayer graphene.” 2DMater. , 035015 (2017).[4] T. Stauber, T. Low, and G. G´omez-Santos, “Chiral Re-sponse of Twisted Bilayer Graphene.” Phys. Rev. Lett. , 046801 (2018).[5] D. Nguyen and D. Son, “Electrodynamics of Thin Sheetsof Twisted Material” arXiv:2008.02812v1 (2020).[6] C. Kim, A. S´anchez-Castillo, Z. Ziegler, Y. Ogawa, C.Noguez, and J. Park, “Chiral atomically thin films” Nat.Nanotechol. , 520 (2016).[7] J. Provost and G. Vallee, “Riemannian Structure onManifolds of Quantum States.” Commun. Math. Phys. ,166803 (2020).[9] Supplementary Information[10] F. Hidalgo, A. S´anchez-Castillo, and C. Noguez, “Effi-cient first-principles method for calculating the circular dichroism of nanostructures.” Phys. Rev. B , 075438(2009)[11] C. Noguez, and F. Hidalgo, “Ab initioelectronic circulardichroism of fullerenes, single-walled carbon nanotubes,and ligand-protected metal nanoparticles.” Chirality ,553–562 (2014).[12] L. D. Barron, “Molecular Light Scattering and OpticalActivity” (Cambridge University Press. 2009).[13] Y. Zhang, Y. Gao, and D. Xiao “Topological ChargePumping in Twisted Bilayer Graphene.” Phys. Rev. B , 041410 (2020).[14] G. Tarnopolsky, A. Kruchov, and A. Vishwanath “Originof Magic Angles in Twisted Bilayer Graphene.” Phys.Rev. Lett. , 106405 (2019).[15] M. Fujimoto, H. Koschke, and M. Koshino “TopologicalCharge Pumping by a Sliding Moir´e Pattern.” Phys. Rev.B , 041112 (2020).[16] M. Kargarian, M. Randeria, and N. Trivedi, “Theory ofKerr and Faraday rotations and linear dichroism in Topo-logical Weyl Semimetals.” Scientific Reports , 12683(2015). Supplementary Information: Anomalous Electrodynamics and Quantum Geometry inGraphene Heterobilayers
STRUCTURES
The kinetic terms t and t in main text derive from the tight-binding model for single-layer graphene t = t , t = t e i ( kx + √ ky ) e i ( kx − √ ky ) e i ( kx + √ ky ) e i ( kx − √ ky ) The interlayer coupling ν ( x ) is given by a series in the lowest G -vectors fixing the three high-symmetry stackingconfigurations AA , AB , and BA . The result is ν ( x ) = t (cid:48) (cid:16) (1 + cos(2 πx/L )) γ + (1 − πx/L )) γ + √ sin(2 πx/L ) γ (cid:17) LOCAL ANTIUNITARY SYMMETRY
Here we show in detail the conditions for vanishing curvature and multipole densities under a local antiunitarysymmetry U K , where K denotes complex conjugation and U is a unitary operator. First, we choose an orthonormalbasis of U K eigenstates with eigenvalues +1. In the case of shear-strained graphene, the energy bands are two-folddegenerate, so we can construct these explicitly with the following algorithm. Readers who are satisfied that such abasis can be constructed may skip this part.First we pick an arbitrary state | v (cid:105) in the band subspace. Since U K is a symmetry of the Hamiltonian, U K| v (cid:105) isan eigenstate with the same energy. Then we construct within the band subspace | u (cid:105) = √ ( | v (cid:105) + U K| v (cid:105) ) , | u (cid:105) = √ e iπ/ ( | v (cid:105) − U K| v (cid:105) )These are eigenstates of U K both with eigenvalue +1. To make them orthogonal, we observe that (cid:104) u | u (cid:105) = e iπ/ ( (cid:104) v | v (cid:105) − (cid:104) v ∗ | v ∗ (cid:105) + (cid:104) v ∗ | U † | v (cid:105) − (cid:104) v | U | v ∗ (cid:105) )= e iπ/ ( (cid:104) v ∗ | U † | v (cid:105) − (cid:104) v | U | v ∗ (cid:105) ) = e iπ/ Im( (cid:104) v ∗ | U † | v (cid:105) )This means that | u (cid:105) and | u (cid:105) are orthogonal if and only if (cid:104) v ∗ | U † | v (cid:105) is real. We can enforce this by replacing | v (cid:105) → e − iθ/ | v (cid:105) where θ = arg( (cid:104) v ∗ | U † | v (cid:105) ) in the above construction of the | u (cid:105) .In this basis | u (cid:105) , we now have that U K| u (cid:105) = | u (cid:105) ⇒ | u ∗ (cid:105) = U † | u (cid:105) The unitary part U can then be inserted into the curvature form as follows (cid:104) ∂u∂k y | M | ∂u∂x α (cid:105) = (cid:104) ∂u∂k y | U U † M U U † | ∂u∂x α (cid:105) = (cid:104) ∂u ∗ ∂k y | U † M U | ∂u ∗ ∂x α (cid:105) = ( (cid:104) ∂u∂k y | ( U † M U ) ∗ | ∂u∂x α (cid:105) ) ∗ ≡ ( (cid:104) ∂u∂k y | M (cid:48) | ∂u∂x α (cid:105) ) ∗ which is the result from the main text.For the multipole density, we follow a similar procedure on the x -projected wavefunctions u ( x ). Since we assume U K is local in x , we have U K u ( x ) = u ( x ) ⇒ u ∗ ( x ) = U † u ( x )Then the multipole density satisfies ρ m = u ∗ ( x ) mu ( x ) = u ∗ ( x ) U U † mU U † u ( x ) = u ( x ) U † mU u ∗ ( x )= ( u ∗ ( x )( U † mU ) ∗ u ( x )) ∗ ≡ ( u ∗ ( x ) m (cid:48) u ( x )) ∗ = ρ ∗ m (cid:48) as desired. KUBO FORMULAE
Here we derive the curvature form (Equation 2 of main text) from Equations 3 and 4 of the main text. Startingwith Equation 3, we have f yα = ieE (cid:88) m (cid:54) = n (cid:104) n | [ ∂∂x α , H ] | m (cid:105)(cid:104) m | τ z | ∂n∂k y (cid:105) (cid:15) m − (cid:15) n − (cid:104) ∂n∂k y | τ z | m (cid:105)(cid:104) m | [ ∂∂x α , H ] | n (cid:105) (cid:15) m − (cid:15) n = ieE (cid:88) m (cid:54) = n (cid:104) n | (cid:16) ( (cid:15) m − (cid:15) n ) ∂∂x α + ∂(cid:15) m ∂x α (cid:17) | m (cid:105)(cid:104) m | τ z | ∂n∂k y (cid:105) (cid:15) m − (cid:15) n − (cid:104) ∂n∂k y | τ z | m (cid:105)(cid:104) m | (cid:16) ( (cid:15) n − (cid:15) m ) ∂∂x α + ∂(cid:15) n ∂x α (cid:17) | n (cid:105) (cid:15) m − (cid:15) n Since | n (cid:105) and | m (cid:105) are orthogonal, the ∂(cid:15) n ∂x α terms drop out. Then we can cancel the energy denominator, giving= ieE (cid:88) m (cid:54) = n (cid:16) (cid:104) n | ∂∂x α | m (cid:105)(cid:104) m | τ z | ∂n∂k y (cid:105) + (cid:104) ∂n∂k y | τ z | m (cid:105)(cid:104) m | ∂∂x α | n (cid:105) (cid:17) Integrating the first term by parts, we have= ieE (cid:88) m (cid:54) = n (cid:16) −(cid:104) ∂n∂x α | m (cid:105)(cid:104) m | τ z | ∂n∂k y (cid:105) + (cid:104) ∂n∂k y | τ z | m (cid:105)(cid:104) m | ∂n∂x α (cid:105) (cid:17) Now we can add and subtract (cid:104) ∂n∂x α | n (cid:105)(cid:104) n | τ z | ∂n∂k y (cid:105) = (cid:104) ∂n∂k y | τ z | n (cid:105)(cid:104) n | ∂n∂x α (cid:105) to complete the sum.= ieE (cid:88) m (cid:16) −(cid:104) ∂n∂x α | m (cid:105)(cid:104) m | τ z | ∂n∂k y (cid:105) + (cid:104) ∂n∂k y | τ z | m (cid:105)(cid:104) m | ∂n∂x α (cid:105) (cid:17) = ieE (cid:16) −(cid:104) ∂n∂x α | τ z | ∂n∂k y (cid:105) + (cid:104) ∂n∂k y | τ z | ∂n∂x α (cid:105) (cid:17) as desired. Equation 4 follows a similar derivation, f xα = ieE (cid:88) m (cid:54) = n (cid:104) n | [ τ z ∂∂x α , H ] | m (cid:105)(cid:104) m | ∂n∂k y (cid:105) (cid:15) m − (cid:15) n − (cid:104) ∂n∂k y | m (cid:105)(cid:104) m | [ τ z ∂∂x α , H ] | n (cid:105) (cid:15) m − (cid:15) n = ieE (cid:88) m (cid:54) = n (cid:104) n | (cid:16) ( (cid:15) m − (cid:15) n ) τ z ∂∂x α + τ z ∂(cid:15) m ∂x α (cid:17) | m (cid:105)(cid:104) m | ∂n∂k y (cid:105) (cid:15) m − (cid:15) n − (cid:104) ∂n∂k y | m (cid:105)(cid:104) m | (cid:16) ( (cid:15) n − (cid:15) m ) τ z ∂∂x α + τ z ∂(cid:15) n ∂x α (cid:17) | n (cid:105) (cid:15) m − (cid:15) n = ieE (cid:88) m (cid:54) = n (cid:16) (cid:104) n | τ z ∂∂x α | m (cid:105)(cid:104) m | ∂n∂k y (cid:105) + (cid:104) ∂n∂k y | m (cid:105)(cid:104) m | τ z ∂∂x α | n (cid:105) (cid:17) = ieE (cid:16) −(cid:104) ∂n∂x α | τ z | ∂n∂k y (cid:105) + (cid:104) ∂n∂k y | τ z | ∂n∂x α (cid:105) (cid:17) again as desired. WILSON LOOP
To calculate the generalized Berry curvature using a Wilson loop, we encounter a difficulty from the vanishing ofthe zeroth order term in each link. The usual matrix-valued Wilson loop isIm log det N (cid:89) i =1 U † i U i +1 where U i is a matrix of n states in the band at the i -th point on the loop, and N + 1 ≡
1. The determinant is invariantunder unitary transformations of the U i , so we can pick a bases for the U i in which the zeroth order part of each linkis the identity. Taking the first order approximation U i +1 ≈ U i + δU i , the Wilson loop becomesIm log det N (cid:89) i =1 (cid:16) I + U † i δU i (cid:17) Taking the first order terms from each subsequent operation, we haveIm log det N (cid:89) i =1 (cid:16) I + U † i δU i (cid:17) ≈ Im log det (cid:32) I + N (cid:88) i =1 U † i δU i (cid:33) ≈ Im log N (cid:88) i =1 n (cid:88) j =1 (cid:104) u ij | δu ij (cid:105) where u ij is the j -th diagonal element of U † i δU i . Finally, the first order part of the log gives us the desired resultIm N (cid:88) i =1 n (cid:88) j =1 (cid:104) u ij | δu ij (cid:105) We would like to modify this derivation with a multipole moment M weighting each of the links. The previoussection showed that a local antiunitary symmetry C x T implies that nonvanishing curvatures correspond to vanishingdensities, so the crucial zeroth order diagonal term in the link U † i M U i +1 vanishes. However, we also have the relationthat nonvanishing densities correspond to vanishing curvatures, suggesting that if we weight the links by ( I + M ),the contribution of I to the curvature will vanish. This way, we can acquire a zeroth order diagonal term in the linkswithout changing the first-order result.We must be careful about zeroth order off-diagonal terms (cid:104) u i | M | u j (cid:105) , i (cid:54) = j which are allowed to be nonzero,imaginary. Due to the complication of tracking these terms, we only prove the following result within a single two-fold degenerate band subspace where each C y symmetry sector contains only a single state. This restriction canbe implemented for low-energy states in Figure 1b of the main paper where the bands do not cross over a regionwell-within the gap.For the following derivation, we utilize a basis with simultaneous C x T and C y symmetry. Starting with a C x T eigenbasis | u (cid:105) , | u (cid:105) constructed using the algorithm in the previous section, we observe that any linear combinationof the basis states with real coefficients is again a C x T eigenstate. We would like to diagonalize C y in the bandsubspace ( C y ) ij = (cid:104) u i | C y | u j (cid:105) , i, j ∈ { , } Using C x T symmetry, we have (cid:104) u i | C y | u j (cid:105) = (cid:104) u i | C x C † x C y C x C † x | u j (cid:105) = (cid:104) u ∗ i | C † x C y C x | u ∗ j (cid:105) = (cid:104) u ∗ i | C y | u ∗ j (cid:105) This shows that C y is real in the C x T basis. It follows that eigenstates of C y will be real linear combinations of | u (cid:105) , | u (cid:105) , which thus remain simultaneous C x T eigenstates.Now suppose that | u i (cid:105) , | v i (cid:105) is a basis of simultaneous C x T , C y eigenstates on the i th point of the Wilson loop,where we adopt the convention that | u i (cid:105) is the +1 C y eigenstate and | v i (cid:105) is the − C y eigenstate. Let U i be thematrix with columns | u i (cid:105) , | v i (cid:105) . We introduce M with a perturbative parameter λ .Im log det N (cid:89) i =1 U † i ( I + λM ) U i +1 = Im log det N (cid:89) i =1 (cid:18) (cid:104) u i | ( I + λM ) | u i +1 (cid:105) (cid:104) u i | ( I + λM ) | v i +1 (cid:105)(cid:104) v i | ( I + λM ) | u i +1 (cid:105) (cid:104) v i | ( I + λM ) | v i +1 (cid:105) (cid:19) = Im log det N (cid:89) i =1 (cid:18) (cid:104) u i | δu i (cid:105) + λ (cid:104) u i | M | δu i (cid:105) λ (cid:104) u i | M | v i (cid:105) + λ (cid:104) u i | M | δv i (cid:105) λ (cid:104) v i | M | u i (cid:105) + λ (cid:104) v i | M | δu i (cid:105) (cid:104) v i | δv i (cid:105) + λ (cid:104) v i | M | δv i (cid:105) (cid:19) where the (cid:104) u i | v i +1 (cid:105) , (cid:104) v i | u i +1 (cid:105) terms vanished because u and v are in different C y symmetry sectors, and the (cid:104) u i | M | u i (cid:105) , (cid:104) v i | M | v i (cid:105) terms vanished because M is assumed to be odd under C x T . The term (cid:104) u i | M | v i (cid:105) is pure imaginary by C x T symmetry, so let us write it as iQ M .= Im log det N (cid:89) i =1 (cid:18) (cid:104) u i | δu i (cid:105) + λ (cid:104) u i | M | δu i (cid:105) iλQ M + λ (cid:104) u i | M | δv i (cid:105)− iλQ M + λ (cid:104) v i | M | δu i (cid:105) (cid:104) v i | δv i (cid:105) + λ (cid:104) v i | M | δv i (cid:105) (cid:19) If we carry out the product keeping only terms to first order in λ and/or δ , we get= Im log det (cid:18) a b ∗ b c (cid:19) for a = 1 + (cid:88) i Here we offer further evidence that the oscillations found in the curvatures and energies originate from an effectivepotential constructed out of the γ part of the potential. To clarify the different character of the low-energy spectrumof a constant versus varying γ potential, these are plotted in Figures S1 and S2 respectively. One can observe that forthe constant potential, not only do the spacings of the bands decrease rapidly to the critical ring edge, but there aretwo independently-evolving spacing periods. This is because the nodal ring, formed by vertical displacement of theDirac cones, is not perfectly symmetric. Under zone-folding, bands displaced from ± G do not perfectly coincide butinstead evolve out of phase. In contrast, the varying γ potential gives rise to a single, slowly evolving spacing period.The qualitative difference suggests different mechanisms determining the energy spacing in these two situations. Theclose agreement between the numerical wavefunctions and the effective potential formulated in Equation 5 of themain paper for differing values of k y (Figure S3) gives strong evidence that the spacing for the varying potential isdetermined by a resonance condition within the effective potential.The preservation of the γ zero energy crossings under incorporation of the γ part of the potential (Figure S4)shows that this mechanism carries over to the full moir´e. Similary between this spectrum and the low bands of thefull moir´e shown in Figure 1 of the main paper support our claim that the γ part of the potential only exertsperturbative influence on this part of the spectrum. FIG. S1. Low energy bands for a zone-folded critical ring in a constant γ potential, showing the two independently-evolvingband spacings.FIG. S2. Low energy bands for a varying γ potential, displaying a single slowly-evolving band spacing in contrast with theconstant γ potential.FIG. S3. Effective potentials and numerical charge densities for k y = 0 . 06 and k y = . 086 indicating close agreement as therange of the AA basin changes.FIG. S4. Low energy bands with the γ and γ components of the moir´e potential. The zero-energy crossings line up withthose of the varying γ4