{'fraction_non_alphanumeric': 0.027004166536136087, 'fraction_numerical': 0.00754988878794524, 'mean_word_length': 4.551652173913044, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anomalous Electrodynamics and Quantum Geometry in the Dirac-Harper model for a Graphene Bilayer

Graphene bilayers with layer antisymmetric strains are studied using the Dirac-Harper model for a pair of single layer Dirac Hamiltonians coupled by a one-dimensional moir\'e-periodic interlayer tunneling amplitude. This model hosts low energy, nearly dispersionless bands near charge neutrality that support anomalous polarizations of its charge multipole distributions. These are analyzed introducing a generalized Berry curvature that encodes the field-induced dynamics of multipole fields allowed in a chiral medium with time reversal symmetry. The formulation identifies a reciprocity relation between responses to layer-symmetric and layer-antisymmetric in-plane electric fields and reveals momentum-space quantum oscillations produced by a spatial pattern of band inversions on the moir\'e scale.


I. INTRODUCTION
The electrodynamic properties of electrons in crystals can include anomalous responses due to the momentum dependence of internal spin and orbital attributes in Bloch bands. There has been interest in interpreting these as geometric phenomena associated with a quantum metric. 1,2 This concept generalizes to driving fields that vary in space. When a smooth spatial variation lowers the symmetry of the system, it can activate responses that would be forbidden for spatially uniform fields. These spatially dispersive responses are sometimes expressed as responses of a system to the field gradients. Artificial structures fabricated by stacking atomically thin layers are of great current interest in this regard, and in these systems the field gradient is replaced by a discretized set of driving fields on separate layers. [3][4][5] In twisted graphene bilayers discretized coupled layer responses can greatly exceed spatially dispersive responses found in conventional chiral materials. 6 Here we formulate coupled-layer electrodynamic responses in Bloch bands as quantum geometric quantities. We find that charge multipole distributions in the bilayer have dynamics that are naturally represented by a generalization of the Berry curvature to a class of curvature forms 7 in a parameter space which combines position and momentum. This approach generalizes the concept of an anomalous charge velocity induced by Berry curvature to a wider class of multipole electrodynamic responses. We illustrate this for the specific case of a graphene bilayer with a one-dimensional moiré periodicity induced by layer-antisymmetric strain.
The paper is organized as follows. In Section II for completeness we briefly review the construction of the Dirac-Harper model for our model system, a shearstrained one-dimensional moiré structure. 8 In Section III * mele@physics.upenn.edu we describe features of the low energy spectrum for this model which reveals a set of weakly-dispersing modes with charge peaks that counterpropogate in the moiré as a function of the moiré-transverse crystal momentum. In Section IV we characterize multipole densities which associate this counterpropagation with a net polarization of a multipole distribution. In Section V we formulate a generalized Berry curvature that describes the anomalous response of these multipole densities to the electric field. In Section VI we examine the symmetry constraints of both the multipole susceptibilities (Section V) and the static multipole polarizations (Section IV) which turn out to be distinct phenomena. In Section VII we explain the origin of oscillatory behavior found numerically in calculations of the multipole susceptibilities.

II. DIRAC-HARPER MODEL
The system we study is a strained graphene bilayer in which the stacking is modulated in a single crystallographic direction x (Figure 1a) while retaining short range lattice periodicity in the orthogonal coordinate y.
Here for completeness we summarize the formulation of this problem as a spinorial Harper model. 8 The Hamiltonian for this system includes a spatially varying term ν(x) which couples the two layers and a set of kinetic terms t which couple adjacent cells in the same layer. The conserved crystal momentum k y allows this system to be separated into a family of one-dimensional problems along the x coordinate where the effects of the y-directed couplings are absorbed into a k y dependence of the kinetic energy terms. The result is a generalization the scalar Harper equation to a spinorial Dirac-Harper problem 8 t(k)ψ n+1 + t † (k)ψ n−1 + t 0 ψ n + ν(x)ψ n = ψ n (1) Here ψ n are four component fields with amplitudes on two sublattices (a and b) and two layers (1 and 2) with the ordering ψ T = (a 1 , b 1 , a 2 , b 2 ). The integers n in arXiv:2102.11850v2 [cond-mat.mes-hall] 12 Aug 2021 Equation 1 index the primitive cells of the microscopic graphene lattice along the x direction. The t(k) are block diagonal 4×4 matrices that can be derived from a nearest neighbor tight-binding kinetic energy for the individual graphene layers. 9 Defining the Pauli matrices σ i and τ i acting on sublattice index and layer index respectively, the products σ x τ 0 = γ 15 and σ y τ 0 = γ 25 form a basis for t.
The interlayer coupling ν(x) = ν(x + L) is a smooth spatial variation through three high-symmetry stacking configurations: AA where equivalent sublattices in the two layers align vertically, AB where a 1 sublattices eclipse b 2 sublattices, and BA where b 1 sublattices eclipse a 2 sublattices. The smoothest interpolation between these three configurations is found using a Fourier series in the lowest one dimensional G-vectors with matrix-valued coefficients ν G . The high-symmetry stacking configurations utilize three γ matrices: γ 4 = σ 0 τ x , γ 14 = σ y τ y , and γ 32 = σ x τ x , so the series can be re-expressed as Figure 1b.  1. (a) A section of the shear-strained (AC) moiré structure studied in this paper, which is taken to repeat indefinitely in the x and y directions. Typical moiré lengths are around 900 lattice constants in the x direction. (b) Scalar amplitudes of the three mass terms γ4, γ14, and γ32 normalized by the interlayer coupling constant t .
The symmetries of this model depend on the orientation of the modulation direction with respect to the underlying graphene lattice, and this crucial feature is encoded in the exact form of the kinetic terms t. An important limiting case studied here occurs for modulation along the armchair (AC) direction, which is achieved by imposing a layer-antisymmetric shear strain (Figure 1a). This produces a chiral structure with the two-fold rotational symmetries around three orthogonal axes, C 2x , C 2y , and C 2z , but no mirror symmetries. In contrast, modulation along the zigzag (ZZ) direction requires uniaxial strain and the resulting structure respects mirror symmetries. The presence of these mirror symmetries is a nongeneric feature of the DH model which occurs only for special orientations of the modulation direction. The kinetic terms are expressed in a low-energy effective theory by taking ψ n±1 ≈ ψ n in Equation 1 and linearizing the combination t(k) + t † (k) + t 0 in small momentum q about a Dirac point. For AC, this linearization gives a kinetic energy operator v F (γ 15 q x − γ 25 q y ) where v F = √ 3at/2, t is the intralayer coupling constant, and a is the graphene lattice constant. To compare, ZZ is a π/2 rotation of AC which swaps the roles of γ 15 and γ 25 . Figure 2d shows the band structure of the AC structure as a function of k y . Features of the spectrum can be understood in the long moiré limit directly from the commutation relations of the kinetic terms with the mass terms ν(x) introduced by the interlayer tunneling. 8 In the absence of interlayer coupling, the spectra consist of a pair of layer-degenerate Dirac cones shown in Figure 2a. In a semi-classical picture, modification of this spectrum by local mass terms describes behavior of Bloch electrons localized to specific regions within the moiré. The γ 4 mass term describes same-sublattice hopping between the two layers, characteristic of an AA-registered region. This mass term commutes with both kinetic terms, splitting two hybridized Dirac cones in energy so that they intersect on a ring in momentum space ( Figure 2b) which we refer to as the critical ring. The mass term γ 14 describes an asymmetry in the coupling between opposite sublattices on the two layers which occurs at a generic position in the superlattice cell and is strongest in the AB and BA regions. This term anticommutes with γ 15 , gapping out the critical ring everywhere a except at two nodal points where the critical ring intersects the q y axis (Figure 2c).

III. LOW ENERGY SPECTRA
We can use these local four-band features to analyze the full moiré spectrum by noticing that in the long-moiré limit, the slowly-varying stacking allows modes to localize in regions of nearly constant stacking. Referring to Figure 1b, the AA region is a location where two out of three mass terms have extrema. Indeed, we can readily identify the width of the region supporting zero-energy states in the full moiré spectrum (Figure 2d) as the peak width of the critical ring in AA. This critical ring occurs at |q| = v F /t a ≡ q c where t is the interlayer coupling strength. The gap opened around zero energy with two point closures at k y = ±q c arises from γ 14 changing its sign in this region.
As shown in the inset of Figure 2d, a set of nearly dispersionless states appears within this gap. These modes can be understood as produced by sign changes of the in- terlayer couplings as a function of lateral position x. 8 One set of modes associated with a sign change of the γ 14 mass confines charge density to the AA region of the moiré. A second set of modes arises from a sign change of the γ 4 mass term giving a k y -dependent potential that produces a charge density ( Figure 3a) with two peaks that counterpropagate as a function of k y . These features migrate from the AB and BA regions towards the AA region as k y approaches the gap closure at k y /q c ≡ q r = 1. In a semiclassical description, this counterpropagation could be driven by an electric field k y = −eE y . This counterpropagation transports no charge along x but can be associated with transport of higher order multipole distribution.

IV. MULTIPOLE POLARIZATION
We begin by inspecting multipole distributions that are antisymmetric in x because these quantities tag the counterpropogating peaks (Figure 3a) with opposite signs and associate them with net polarization of a charge-neutral multipole field. The lowest nontrivial multipole distribution that respects C 2x , C 2y , and C 2z symmetries of the shear-strained AC moiré is an octopole. A multipole distribution that is odd in x can therefore be produced by a quantity that has the local rotational symmetry of a quadrupole in the yz plane. Such an octopolar dis- tribution is already explicitly present in the DH model; the mass termγ 14 sin(Gx) is a product of a quadrupolar operatorγ 14 with an antisymmetric function of x that distinguishes AB and BA registered regions. The multipole distribution ρ 14 (x) = ψ † (x)γ 14 ψ(x) manifests this property in an observable by coupling the quadrupolar γ 14 to an antisymmetric spatial distribution (Figure 3b). Similarly, a quadrupole density locally measured by the sublattice-and layer-odd operatorγ 5 = σ z τ z describes a local layer-and sublattice-polarization that is also spatially antisymmetric (Figure 3b). We refer to transverse polarizations intrinsic to these antisymmetric multipole distributions as anomalous polarizations.

V. ELECTRODYNAMIC RESPONSES
The response of an anomalous multipole polarization to an applied electric field can be expressed as a quantum geometric quantity. In a Bloch band with broken time reversal or inversion symmetry, an anomalous velocity for charge is produced by the Berry curvature ∂ λµ u|∂ λν u − c.c.. To extend this to multipole dynamics, we examine the generalizated curvature tensor where M ≡ m αβ ⊗ I N and m αβ is a local multipole operator weighting the sublattice and layer degrees of freedom. For m αβ = I, we recover the usual Berry curvature responsible for the anomalous Hall conductance. More generally, ω m can describe the anomalous dynamics of a charge-neutral multipole distribution. A mathematical discussion of this curvature is given in the Supplementary Information. 9 The parameters λ specify a tangent space for a state in the Bloch band which will determine the association of this curvature with a physical response function. To identify appropriate choices for these parameters for the Dirac-Harper model, we observe that the onedimensional moiré structure is most naturally described using its two extended coordinates: a momentum coordinate in the y direction and a spatial coordinate in the x direction. The two dual coordinates k x and y are insignificant: zone-folding reduces the width the Brillouin zone G x by a factor of the inverse moiré length 1/L, and the y coordinate does not appear in the DH equation. For definiteness we choose k x = 0, promoting the model with symmetries that can be only weakly violated at order 1/L.
We decompose the extended spatial coordinate by introducing three parameters, each of which controls the origin of its respective mass term in the Hamiltonian. These spatial parameters x α appear in ν(x) as This parameterization allows us to vary the origin for each mass field individually. Each x α paired with k y parameterizes a two-dimensional subspace of the projective Hilbert space. Note that the ground state lies in the one-dimensional sector at x α = 0 where the mass fields are pinned by the moiré structure. Nonetheless the two-dimensional tangent space generated by differential changes ∂ xα , ∂ ky will encode an observable property of these states which we denote Notice that K (M ) α must have the symmetry of an octopole to respect C 2 symmetry. This is satisfied if m has the symmetry of a z-directed dipole, and we use τ z as a simple and useful example.
The triad of operators ∂ xα , ∂ ky , τ z can now be identified with two different response functions. The operator ∂ ky couples the band to a perturbing electric field applied in the y-direction, and the operator ∂ xα arises from the observable k (α) x = [∂ xα , H] which describes an x-directed force on a component of the wavefunction coupled to the α-th mass term. The coupling to individual mass terms will be explained in more detail in the next paragraph, but for now it is adequate to think of this as a force acting on component of the charge density. To introduce the third member of the triad, τ z , we can consider either a layer-antisymmetric electric field or a layer-antisymmetric response. For the former, τ z enters in the driving term Eτ zŷ , and the force is given by where the first subscript labels the force by an induced lattice effect, to be described later. For the latter, τ z weights the observable i∂ xα with opposite sign in each layer: Both of these can be formally expressed as manifestation of same generalized curvature K (τz) α in Equation 3, and thus they are equivalent. 9 Numerical calculation for each mass channel α is plotted as a function of q r in Figure 4b. These results can be regarded as variants of the constitutive relations that describe natural optical activity in chiral media. 3,10-12 Natural optical activity arises from a correlated electric dipole-quadrupole response driven by an optical field, while the multipole susceptibilities obtained here describe an analogous coupling of two different charge multipoles in the static limit. Each mass channel describes a force acting on the component of the charge density confined by the associated mass term. This is demonstrated by examining the differential changes of the charge density produced by shifts of each x α , shown explicitly in Figure 4c. The AA peaks defined by the γ 14 mass inversion are shifted by x 14 (i) and the AB/BA peaks defined by γ 4 inversion are shifted with x 4 (ii). Coupling to x 32 produces almost no effect (iii) because γ 32 does not confine charge in the moiré cell via a sign change. 8 The combination of all three x α produce uniform charge shift (iv), so the sum of K α over mass channels a gives the anomalous susceptibility for polarization of a multipole-weighted density (τ z in this example).
The equivalence of Equations 4 and 5 identifies a reciprocity relation stating that a force induced by a multipole-weighted field is equal in magnitude to the multipole-weighted force induced by a uniform field. For our example of a z-directed dipole τ z , these responses are driven by a layer-antisymmetric field E and a layer symmetric field E ⇒ respectively. The layerantisymmetric field can be thought of the low-frequency limit of a plane wave polarized in the y direction. For case of E (Equation 4), the response summed over mass channels produces a force that acts to depin charge in the negative x direction as shown in Figure 4c iv. The E ⇒ case (Equation 5) corresponds to a similar displacement of charge, but weighted in opposite directions in the two layers.
Displacements of various components of the charge density can produce a counterforce acting on the lattice to restore the moiré to an equilibrium stacking arrangement with respect to the charge. The coupling of lattice strain fields to the applied field must respect the C 2 symmetries, so we can deduce the allowed combinations. Focusing our discussion on uniform layer-antisymmetric displacements, the field E (odd in y and z) couples to displacement along y (also odd in y in z), and the field E ⇒ (odd only in y) couples to displacement along x (odd in x and z). The subscripts in Equations 4 and 5 are now clear: f y (f x ) induces a layer shear along y (x). Intuitively, the y-directed shear transports the full moiré pattern along the x direction 13 to bring the lattice back in phase with the charge, and the x-directed shear counters the layer antisymmetric distortion of the charge along x induced by the initial electrodynamic response. We can write a general y-directed electric field E y (z) as a sum of E + E ⇒ , where E = [E y (+z) − E y (−z)]/2 and E ⇒ = [E y (+z) + E y (−z)]/2. Then the induced layer shear acts at a general angle φ in the xy plane, which using reciprocity of the curvature is given by φ = tan −1 (f y /f x ) = tan −1 (E /E ⇒ ). φ = 0 corresponds to opposing layer displacements in the x direction driven by a uniform electric field in the y direction.

VI. SYMMETRY REQUIREMENTS
The response functions described by K α are symmetric under time-reversal T : since K α is imaginary and contains only one momentum coordinate, it accumulates two cancelling sign changes under T . This means that the anomalous response described by K α does not require T -breaking to be present. In contrast to the usual Berry curvature which vanishes everywhere in a Bloch band that respects both time-reversal and inversion symmetry, the generalized susceptibility studied here vanishes everywhere for systems with time-reversal and mirror symmetry. We can gain additional insight using C 2x T symmetry, which is a local symmetry in the parameter space (x, k y ). The symmetry of m αβ under this operation determines whether the curvature K (M ) or its associated multipole density ρ m must vanish locally in the mixed momentum-position space. For the curvature, one finds 9 ∂u ∂ky |M | ∂u ∂xα = ∂u ∂ky |M | ∂u ∂xα * for M = (C † 2x M C 2x ) * . When M = M , the curvature is equal to its complex conjugate and therefore it vanishes everywhere. Indeed this is the case for the charge susceptibility M = I, but the generalized extensions can incorporate multipoles such as τ z where M = −M , allowing this response to be nonzero. For the multipole density ρ m , the same manipulation gives ρ m = ρ * m where m is the x-projected version of M . Since ρ m is real, this quantity vanishes when m = −m. Thus we find that the symmetry of m under C 2x T divides multipole densities into two distinct classes: those that can have nonvanishing density but vanishing curvature and those that have vanishing density but a nonvanishing curvature. The anomalous response measured by Equation 3 therefore is physically distinct from the construction obtained by antisymmetrically-tagging the counterpropagating peaks shown in Figure 3a. Dynamical effects occur in quantities with no bulk accumulation of a multipole density, and quantities with bulk accumulation can only display a static polarization.
For unitary multipole operators, 9 the generalized Berry curvature is related to the standard Berry curvature by diffeomorphisms on the Hilbert space. This is a consequence of the Darboux theorem which states that any symplectic two-form can be re-expressed in a standard form on an open set via a coordinate change. Thus, we arrive at an alternative interpretation of the generalized Berry curvature as the standard Berry curvature computed on a different sector of the Hilbert space. As an example, we can put the τ z curvature into standard form by complex-conjugation of the second-layer amplitudes of the wavefunctions, which has the effect of negating the second-layer component of the complex curvature. The new wavefunctions do not have obvious physical interpretation, but the standard form of the Berry curvature lends computational and conceptual power. Calculations of the generalized Berry curvature can be very difficult in non-standard form, but after applying the diffeomorphism, a standard Wilson loop can be utilized. 9 Figure 4b oscillate as a function of k y . This behavior manifests oscillations in the energies of the AB/BA bands, partially shown in Figure 2d, which are encoded in the response function via the velocity terms in the Kubo formula. The observed oscillations have a smoothly evolving period which is nearly constant over the width of the band-inverted region |q r | < 1 in momentum space. We can understand this in a minimal model by examining the low-energy behavior of spatially varying γ 4 potential, shown in Figure 5. One might be tempted to interpret the spacing of zero energy crossings as a result of zone-folding of the critical ring (Figure 2b) into the first moiré Brillouin zone. Such an effect arises from a constant γ 4 potential, which we can inspect numerically using a tightbinding model. This calculation reveals that a nonnegligible ±q x anisotropy of the critical ring produces two independently-evolving series of zero-energy crossings in the band structure that rapidly converge near the edge of the critical ring q r = ±1. 9 Such behavior is quite different from the slowly modulated period of the oscillations we observe in the τ z curvature. By contrast, a spatially varying γ 4 potential produces a nearly uniform pattern of band crossings (Figure 6a) which closely matches the period of the curvature oscillations. This behavior can be understood by squaring the low-energy Dirac equation,

All three of the curvatures shown in
to form the effective potential problem The spacing of zero-crossings in this situation arises from a resonance condition within a basin where V eff is negative. The wavefunctions produced by the ν 4 (x)γ 4 potential in a tight-binding calculation (Figure 5b) oscillate in the region where V eff is negative and decay in regions where it is positive. Wavelengths of the oscillating part of the wavefunctions match the value of √ V eff averaged over each period. Consecutive zero energy band crossings differ by a single wavefunction oscillation within the V eff basin, which is indicative that the spacing of the zero-energy crossings is a consequence of wavefunction resonance within this effective potential.
These observations about the effective potential produced by a pure ν 4 potential carry over to the full DH Hamiltonian. The mass term proportional to γ 14 sin(Gx) has the effect of opening a gap around the lowest energy band manifold without changing the spacing of these zero energy crossings. 9 Away from the saddle point configuration in the stacking texture (SP ), γ 32 decreases to zero and has only a perturbative effect. Therefore, oscillatory behavior in the energies is inherited from the k y spacing of zero energy crossings intrinsic to the γ 4 part of the potential. Although the full moiré AB/BA wavefunctions decay within the positive effective potential region, the sensitivity of the energies to these special values of k y suggests that the wavefunctions retain a small oscillatory part stretching over the positive potential region flanked by the AB/BA charge peaks. Interestingly these oscillations would be absent from chiral models for twisted bilayer graphene which exclude this coupling to enforce a chiral symmetry. 14 Here we see that they control a physically measurable response, and measurement of these oscillations can be used to determine the scale of the breaking of the chiral symmetry.

VIII. CONCLUSION
The curvature forms developed here can also be applied to two dimensional twisted bilayers. In this case the DH theory generalizes to a two dimensional theory in a moiré supercell and the forces driving the multipole distributions analogous to Equations 4 and 5 are similarly promoted to two dimensional vectors obtained by integrating over a population in the 2D folded Brillouin zone. The simplest case couples the driving fields to the total charge density recovering the charge pumping phenomena proposed for twisted graphene bilayers. 13,15 The generalized Berry curvature provides a unified formulation that associates this electromechanical response with the more general problem of manipulating multipole distributions using applied fields. The generalized constitutive relations for electric field-driven multipole responses also provide a geometric formulation of the natural optical activity of twisted graphene bilayers studied experimentally 6 and theoretically [3][4][5] and to chiral Weyl semimetals. 16 This work was supported by the Department of Energy under grant DE-FG02-84ER45118.

I. STRUCTURES
The kinetic terms t 0 and t in main text derive from the tight-binding model for single-layer graphene The interlayer coupling ν(x) is given by a series in the lowest G-vectors fixing the three high-symmetry stacking configurations AA, AB, and BA. The result is

II. GENERALIZED BERRY CURVATURE
The usual Berry curvature on the projective Hilbert space CP 4N takes the form where u iα denotes the wavefunction component for the αth sublattice of the ith cell. We insert the multipole by contracting a tensor m αβ on the sublattice degrees of freedom When ω m is applied to tangent vectors ∂ λν and ∂ λµ at |u (where λ µ , λ ν parametrize the band subspace) we obtain the familiar form of the Berry curvature weighted by the multipole M ≡ m αβ ⊗ I N in Equation 2 of the main text. As long as m αβ is invertible, the generalized curvature is related to the usual Berry curvature by a diffeomorphism. This can be appreciated by noting that the generalized Berry curvature is a symplectic form on CP 4N , and by the Darboux Theorem for every neighborhood U there is a choice of coordinates in which ω m takes the standard form of the Berry curvature ω. The diffeomorphism generally changes the sector of the Hilbert space associated with a Bloch band, implying that integration is performed over different submanifolds of the Hilbert space and yields independent results. Furthermore, unless the multipole operator satisfies certain conditions, the diffeomorphism may map to unnormalized wavefunctions. Assuming that the multipole operator is Hermitian, the condition for preserving normalization is that the multipole operator must have eigenvalues of unit modulus, i.e. no scaling of coordinates occurs, only a unitary transformation. When this condition is satisfied, we may use the diffeomorphism to obtain a set of wavefunctions on which the usual calculational techniques and topological analysis related to the Berry curvature can be applied to yield results for the multipole-weighted curvature. In particular, a Wilson loop calculated over the new wavefunctions provides a simple method for calculating the generalized curvature.
To give a concrete example, when m = τ z , the coordinate change is complex conjugation of the second-layer amplitudes: We see that We can thus view the ω τz curvature integrated over a Bloch band as equivalent to the standard Berry curvature integrated over a different sector of the Hilbert space, obtained by complex conjugating all second-layer amplitudes in the wavefunctions.

III. LOCAL ANTIUNITARY SYMMETRY
Here we show in detail the conditions for vanishing curvature and multipole densities under a local antiunitary symmetry U K, where K denotes complex conjugation and U is a unitary operator. First, we choose an orthonormal basis of U K eigenstates with eigenvalues +1. In the case of shear-strained graphene, the energy bands are two-fold degenerate, so we can construct these explicitly with the following algorithm. Readers who are satisfied that such a basis can be constructed may skip this part.
First we pick an arbitrary state |v in the band subspace. Since U K is a symmetry of the Hamiltonian, U K|v is an eigenstate with the same energy. Then we construct within the band subspace These are eigenstates of U K both with eigenvalue +1. To make them orthogonal, we observe that This means that |u 1 and |u 2 are orthogonal if and only if v * |U † |v is real. We can enforce this by replacing |v → e −iθ/2 |v where θ = arg( v * |U † |v ) in the above construction of the |u .
In this basis |u , we now have that The unitary part U can then be inserted into the curvature form as follows ∂u ∂ky |M | ∂u ∂xα = ∂u ∂ky |U U † M U U † | ∂u ∂xα = ∂u * ∂ky |U † M U | ∂u * ∂xα = ( ∂u ∂ky |(U † M U ) * | ∂u ∂xα ) * ≡ ( ∂u ∂ky |M | ∂u ∂xα ) * 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 Then the multipole density satisfies  = ieE − ∂n ∂xα |τ z | ∂n ∂ky + ∂n ∂ky |τ z | ∂n ∂xα again as desired.

V. ENERGY OSCILLATIONS
Here we offer further evidence that the oscillations found in the curvatures and energies originate from an effective potential constructed out of the γ 4 part of the potential. To clarify the different character of the low-energy spectrum of a constant versus varying γ 4 potential, these are plotted in Figures S1 and S2 respectively. One can observe that for the constant potential, not only do the spacings of the bands decrease rapidly to the critical ring edge, but there are two independently-evolving spacing periods. This is because the nodal ring, formed by vertical displacement of the Dirac cones, is not perfectly symmetric. Under zone-folding, bands displaced from ±G do not perfectly coincide but instead evolve out of phase. In contrast, the varying γ 4 potential gives rise to a single, slowly evolving spacing period. The qualitative difference suggests different mechanisms determining the energy spacing in these two situations. The close agreement between the numerical wavefunctions and the effective potential formulated in Equation 6 of the main paper for differing values of k y ( Figure S3) gives strong evidence that the spacing for the varying potential is determined by a resonance condition within the effective potential.
The preservation of the γ 4 zero energy crossings under incorporation of the γ 14 part of the potential ( Figure S4) shows that this mechanism carries over to the full moiré. Similary between this spectrum and the low bands of the full moiré shown in Figure 1 of the main paper support our claim that the γ 32 part of the potential only exerts perturbative influence on this part of the spectrum.