Chiral Approximation to Twisted Bilayer Graphene: Exact Intra-Valley Inversion Symmetry, Nodal Structure and Implications for Higher Magic Angles

Jie Wang, ∗ Yunqin Zheng, 3 Andrew J. Millis, 4, † and Jennifer Cano 5, ‡ Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan Kavli Institute for the Physics and Mathematics of the Universe, University of Tokyo, Kashiwa, Chiba 277-8583, Japan Department of Physics, Columbia University, 538 W 120th Street, New York, New York 10027, USA Department of Physics and Astronomy, Stony Brook University, Stony Brook, New York 11974, USA


I. INTRODUCTION
When one graphene layer is stacked on top of another layer with small relative twist angle, a moiré super-lattice pattern is created. At particular twist angles, referred to by Bistritzer and MacDonald as "magic angles" 1 , the bands near the chemical potential are dramatically flattened and separated from other bands [1][2][3][4] . Experiments on "magic angle" bilayers report interesting phenomena including superconductivity, interaction-driven insulating states and anomalous Hall effects  .
There are eight flat bands arising from the combinations of degrees of freedom in the conduction bands of the component graphene layers [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44] . The states that comprise these bands may be labeled by a spin degree of freedom and two additional indices labeling the layer and sublattice of the component graphene sheets. Much of the novel physics of twisted bilayer graphene is believed to arise when the symmetries corresponding to these quantum numbers are spontaneously or explicitly broken. Interestingly, many of the broken symmetry states appear to have a topological character, revealed for example by anomalous Hall effects 7,8,[45][46][47][48] , and at least at integer filling the topological character is believed to be inherent in the single-particle wavefunctions. An improved understanding of the single-particle wavefunctions is therefore important both for improved understandings of the observed and potentially observable topological phases and as a basis for theories of interaction effects in magic angle bilayer graphene.
Recently, Tarnopolsky, Kruchkov and Vishwanath 49 drew attention to a particular "chiral" model in which the interlayer tunneling Hamiltonian contains no terms in which an electron hops from one layer to the same sublattice on the other layer. They showed that in this case the eight weakly dispersing bands become exactly flat (dispersionless) at certain twist angles. They further constructed explicit expressions for the zero mode wavefunctions, and noticed that their solutions exhibited a holomorphic character reminiscent of the lowest Landau level quantum Hall physics 49,50 . This holomorphic character can give rise to a nontrivial topology of the flatbands, explaining the anomalous Hall effects.
In this paper, we study the zero mode wavefunctions of the chiral model 49 of twisted bilayer graphene in more detail. We identify an exact intra-valley inversion symmetry of the chiral model and show how this symmetry implies that the flatband wavefunctions found by Tarpolsky, Kruchkov and Vishwanath can be written (up to a normalization factor) as: where Φ k is a quantum Hall wavefunction of the lowest Landau level, the function G(r) can be interpreted as a quantum Hall wavefunction in a magnetic field oppositely directed to that of Φ k , and η = ±1 is the inversion eigenvalue. The entire dependence on the crystal momentum k is carried by the quantum Hall wavefunction Φ k , which exhibits one node at a k−dependent position, while G, which is independent of k, has a number of nodes that increases as the magic angle index increases, indicating a similarity between higher magic angles and higher Landau levels. This structure is revealed in FIG. 1, which for the first three magic angles presents the norm of each component of φ k for the case where k is fixed at the moiré Dirac point K and implies a charge variation that can be detected by scanning tunneling spectroscopy.
We show that Eqn. (1) explains how the wavefunction φ k can simultaneously have the Abelian translation symmetry of the usual Bloch wavefunction and give rise to the anomalous Hall effect. Further, the quantum Hall anti-quantum Hall structure implies that the wavefunction nodes have a chirality and we find that nodes of a given chirality are concentrated in particular regions of the unit cell, implying intra-cell circulating currents that grow in magnitude as the magic angle increases. The increased density of nodes at higher magic angles will also affect the project of electron-electron interactions onto the flatbands in a manner similar to that occurring at higher Landau levels in the quantum Hall problem. (1) at the moiré Dirac point K, plotted at the first three magic angles (three rows). The upper left and lower right corners of the unit cells are the BA (r0) and AB (−r0) stacking points, as marked. The two columns correspond to the bottom and top components of the wavefunction. Each of them has clear symmetry and zero-structures. The wavefunctions at other Bloch momentum have a similar zero-structure, as explained in the text. The zeros are classified by their chirality i.e. whether the wavefunction's phase advances by ±2π when encircling the zero once. Remarkably, the wavefunction associated with the n th magic angle has 3(n − 1) zeros located at the unit cell center, all of which have the same chirality. We discuss the mathematical structure in Section IV and Section V, and implications for experimental observables in Section VI.
The paper is organized as follows. Section II reviews the continuum model of twisted bilayer graphene and the chiral model defined from it, to establish the notation and approximations used here. Section III introduces our intra-valley inversion symmetry and derives some of the FIG. 2: Left: moiré Brillouin zone. Right: real space moiré unit cell. Shown are the reciprocal lattice vectors b1,2, the moiré Dirac points K, K , the real space lattice vectors a1,2 and the wavevectors q0,1,2, and the BA stacking point r0. The AB and AA stacking points are located respectively at −r0 and the origin of the unit cell.
properties that follow from it. Then in Section IV, we reexamine the derivation of the flatband wavefunctions and derive their spinor-structure. We then discuss the nodal structure of the flatband wavefunctions in Section V. In the last part of this work, Section VI, we discuss how our findings can impact experimental observables. Section VII is a summary and conclusion.

II. MODEL HAMILTONIANS
We start this section by reviewing the continuum model 1-3 and the chiral model 49 of twisted bilayer graphene to establish the notation.
When two parallel graphene sheets (top, bottom) are stacked with any one of an infinite set of relative commensurate angles θ, a moiré pattern forms, in which the combined system retains the basic hexagonal lattice structure of graphene, but with a much larger unit cell containing a number of carbon atoms ∼ θ −2 . The corresponding reciprocal space unit cell, which we refer to as the moiré Brillouin zone, is illustrated in FIG. 2.
As shown in FIG. 2, a i=1,2 indicate the two dimensional basis vectors of the moiré unit cell. The area of the moiré unit cell is 2πS=|a 1 × a 2 |. We denote the reciprocal lattice vectors as b i=1,2 . Throughout this paper, we define the unit length by setting √ S=1. The fundamental single-particle Hamiltonian for twisted bilayer graphene consists of a standard singlelayer graphene Hamiltonian for the top/bottom layer, h G (r, r ), and an interlayer coupling T (r, r ) whose periodicity defines the moiré superlattice. Schematically the Hamiltonian operator is where t/b stands for the top/bottom graphene sheet. It is generally agreed that, as proposed by Bistritzer and MacDonald 1 , the low energy properties of twisted bilayer graphene can be adequately described by a model with three key features. The first is a continuum description of the physics in each graphene sheet, obtained by linearizing the graphene Hamiltonian h G near the Dirac points (we denote the linearized Dirac Hamiltonian as h D ). The second is that the interlayer hopping only couples states near one Dirac point in one layer with states near the same graphene Dirac point in the other layer. This means that the relevant Hamiltonian is the product of two copies, one for each valley. A third simplification proposed by Bistritzer and MacDonald is that the interlayer hopping, in principle a function of r in one layer and r in the other becomes a function only of r with r =r. This is a coarse-graining approximation based on the notion that T (r) has a range of the order of the carboncarbon distance so if the wavefunctions vary slowly on this scale we can ignore the detailed local structure.
Following Bistritzer and MacDonald 1 , the effective continuum Hamiltonian of a single valley is, (3) A related Hamiltonian can be found for the opposite valley by acting with time reversal symmetry. The wavefunction Ψ BM (r) is a four-component spinor, with the lower two components the two sublattices of the top layer, and the upper two the two sublattices of the bottom layer: We have suppressed the spin index because the global SU (2) spin invariance implies that the single-particle Hamiltonian is spin-diagonal. The continuum approximation to the Dirac Hamiltonian of a layer λ = t, b is: where K t/b + is the graphene Dirac point K + rotated by ±θ/2. As shown in FIG. 2, we define the moiré Dirac points as is the moiré Gamma point labeled in graphene's reciprocal lattice coordinates. The interlayer tunneling potential T (r) is constrained by the symmetries of a single valley: C 3 , M y and C 2 T , as discussed in Section III A. In the Bistritzer-MacDonald model, the interlayer hopping is with φ=2π/3, the T j is: The chiral model 49 is obtained by setting ω 0 = 0 in Eqn. (7). The chiral model for a single valley is written in a different basis as H BM in Eqn. (3): where Ψ c (r) is a four-component spinor whose upper two components (φ) correspond to the A sublattice of the bottom and top layer, and the lower two components (χ) the B sublattice of the bottom and top layer: where we have suppressed the Bloch momentum k.
The unitary transformation between the non-chiral basis Eqn. (4) and the chiral basis Eqn. (9) is: where in Eqn. (10), we have used σ and τ for Pauli matrices acting on the sublattice and layer degrees of freedom respectively: σ : sublattices; τ : layers.
In Eqn. (10), we have also shifted the center of the Bloch momentum of the chiral basis to the moiré Gamma point. The Bloch boundary condition of the chiral basis is: where the details of Eqn. (10) and Eqn. (11) can be found in Appendix A. The operators D † (r) and D(r) in Eqn. (8) are: where U φ (r) is: As usual, we have defined z=x+iy, ∂=∂ x −i∂ y . The parameter α is determined by the twisted angle: α=(3w 1 a 0 )/(8πv 0 sin θ 2 ) where v 0 is the graphene's Fermi velocity and a 0 is the graphene's lattice constant. The vectors q 0,1,2 are specified in FIG. 2.
The chiral Hamiltonian anti-commutes with the chiral matrix σ z . As a consequence, the single-particle spectrum is particle-hole symmetric. In the next section, we review symmetries of twisted bilayer graphene, and introduce the intra-valley inversion symmetry.

III. INTRA-VALLEY INVERSION SYMMETRY
In this section, we start by discussing the symmetries of twisted bilayer graphene with an emphasis on how C 2 T symmetry constrains the tunneling terms. In Section III B, we introduce the exact intra-valley inversion symmetry of the chiral model, and derive some properties that follow from it.
The "crystal symmetries" of twisted bilayer graphene are generated by the moiré translation symmetry, C 6 rotational rotation, and mirror symmetry M y . Time reversal symmetry, T , is also present. In addition, in the continuum model the charge conservation of each valley, i.e. U (1) valley symmetry, is assumed. The symmetries that keep each valley invariant (C 2 T , C 3 and M y ) constrain the single valley Hamiltonian in Eqn. (3) and Eqn. (8). Here, the important constraint for us is that C 2 T symmetry requires the tunneling term in Eqn. (3) satisfy (proof in Appendix B): where, as in the previous section, σ x acts on sublattice space. In the chiral basis, this means that the off-diagonal elements of D (and D † ) are related by r ↔ −r, as we shown in Eqn. (12).

B. Exact intra-valley inversion symmetry of the chiral model
Here we show that the chiral model enjoys an exact intra-valley inversion symmetry as constrained by C 2 T symmetry and the linearized Dirac fermion. We then discuss properties that follow from it, including the symmetries of the spectrum and single-particle states. In the end, we show a numerically observed alternating pattern of magic angle inversion parities. Lemma 1. The zero-mode operator satisfies, The calculation follows from the C 2 T constraint in Eqn. (14) and the definition of D † in Eqn. (12). For Lemma 1 to hold, we need the off-diagonal elements of the D(r) operator to be related by r ↔ −r. In other words, in the chiral basis (Eqn. (9)), the interlayer tunneling potential from top to bottom layer is identical to that from bottom to top layer with spatial inversion. As we discussed in Section III A, this is guaranteed by the C 2 T symmetry.
We now define the intra-valley inversion symmetry.
Theorem 1. The chiral model of twisted bilayer graphene has an exact intra-valley inversion symmetry, whose operator is, such that, Again, σ and τ are Pauli matrices acting on the sublattice and layer degrees of freedom respectively.
Proof. It is straightforward to prove by using Lemma 1: We call Eqn. (17) the intra-valley inversion symmetry in order to distinguish it from the crystalline 2D inversion symmetry, C 2 . Since the C 2 symmetry mixes valleys of twisted bilayer graphene, it is not a symmetry of the single-valley continuum models in Eqn. (3) and Eqn. (8). In contrast, the intra-valley inversion maps k to −k within the moiré Brillouin zone and thus does not mix valleys. As shown in Eqn. (17), the intra-valley inversion symmetry is an exact symmetry for the single valley chiral model Eqn. (8).
We emphasize that the intra-valley inversion symmetry requires no extra assumptions beyond the chiral model in Eqn. (8). The only requirement is the crystal symmetry C 2 T and the linearized Dirac fermion, which are already present in the chiral model in Eqn. (8).
The τ y operator of Eqn. (15) has appeared in recent literature. In Ref. (51) it is referred to as the involution operator. In Ref. (35) and a very recent paper Ref. (40), by the same authors, a similar operator iτ y plus r ↔ −r is termed the unitary particle-hole operator. This is different than our intra-valley inversion symmetry: our τ y operates on the chiral basis in Eqn. (9), while the "unitary particle-hole" acts on the non-chiral basis in Eqn. (4). Since the unitary transformation between these two bases in Eqn. (10) does not commute with τ y , these two symmetries are distinct. It is also important to emphasize that our intra-valley inversion is an exact symmetry of the chiral model, while the unitary particle-hole symmetry is an approximate symmetry for both the continuum model in Eqn. (3) and the chiral model in Eqn. (8), according to Ref. (35) and Ref. (40).
Many interesting facts follow from the intra-valley inversion symmetry, as we describe here and in the next section.
Corollary 1. At all twist angles, the single particle spectrum of the chiral model is not only particle-hole symmetric, but also inversion symmetric.
This follows directly from Theorem 1. Denote the sublattice A/B wavefunctions as where each of φ k and χ k is a two-component spinor representing the bottom and top layer's degrees of freedom.
If we know Ψ k as an eigenstate of energy E at Bloch momentum k, then IΨ k (−r) is the eigenstate of the same energy but with an opposite Bloch momentum: We have thus proved the spectrum inversion symmetry by explicitly constructing eigenstates of the same energy and opposite Bloch momentum. This construction in fact also illustrates a spinor structure of the eigenstates.
Theorem 2. At all twist angles for any Bloch momentum k, there exists a phase ζ k , such that, Proof. Below Corollary 1, we explicitly constructed the eigenstate of opposite Bloch momentum. At nondegenerate k, our constructed wavefunction must be proportional to the wavefunction at −k up to a U (1) phase, from which Eqn. (19) follows immediately. The fact that the phase ζ k is anti-symmetric is seen by applying Eqn. (19) twice. For degenerate zero modes, one can label them by the chiral eigenvalue and find the same conclusion.
Theorem 2 can be regarded as a gauge fixing condition. One can perform gauge transformations, to tune the ζ k field: The only obstruction of such tuning is at inversion symmetric points where k inv = −k inv modulo reciprocal lattice vectors: there ζ k and ζ −k cancel, and ζ kinv is either 0 or π. In practice, the intra-valley inversion eigenvalue can be read off from the transformation property of the moiré Gamma point wavefunction (or from wavefunctions at other k inv ): We numerically observed (for the first three magic angles) that there is a coincidence between the zero mode's intra-valley inversion eigenvalue and the parity of magic angle: we found η=+1 for the 1 st , 3 rd magic angles, while for the 2 nd magic angle η=−1. In FIG. 3, we monitored the evolution of the low lying eigenstates at the Gamma point from the 1 st to the 2 nd magic angle (we plot non-negative energies only since the full spectrum has particle-hole symmetry). Eigenstates at the Gamma point are either singlet or doublet, as they are one and two dimensional irreducible representation of the symmetry group (generated by C 3 and M y ) 33  We hypothesize that the alternating parity of magic angle zero mode wavefunctions is a generic feature and will hold for all magic angles: that is, the inversion eigenvalue of n th magic angle flatband wavefunction is −(−1) n .

IV. ZERO MODE WAVEFUNCTIONS AND THE SPINOR STRUCTURE
In this section, we reexamine the zero mode solution of Ref. (49), derive the spinor structure of the zero mode wavefunction as shown in Eqn. (1), and show its intricate relation to quantum Hall physics.

A. Chiral model and the zero modes
Following Ref. (49), we show that the chiral Hamiltonian in Eqn. (8) has two zero modes. The eigenvectors of these two zero modes must satisfy: In Ref. (49), the authors found and proved that, for a discrete series of values of α (corresponding to magic angles), the chiral model admits exact flatbands. At the crux of their analysis is the fact that at magic angles, both components of the moiré Dirac K point wavefunction φ K = (φ b,K , φ t,K ) T vanish at a common point: the BA stacking point, which permits an explicit construction of the zero mode wavefunctions. Here we review several key steps (Theorem 3 to Theorem 5) of Ref. (49) in deriving the zero mode wavefunctions. We refer the readers to Ref. (49) for more details. A crucial step in deriving the zero mode solutions in Ref. (49) is Theorem 3, which follows from the translation and C 3 rotation symmetry: is independent of r.
Proof. It is straightforward to find v F (α) is holomorphic, i.e.∂v F (α) = 0, by using the zero mode equations that φ l,K satisfy. Then, v F (α) must be a constant since it is also cell-periodic.
At magic angles (where the low lying two bands become dispersionless), the Fermi velocity goes to zero. Since the top component of φ K vanishes for all twist angles at ±r 0 , it follows from the vanishing Fermi velocity that the bottom component must at least have one common zero with the top component, at either +r 0 or −r 0 . In fact, the exact flatband condition coincides with the condition that two components of the wavefunction have a common zero, as pointed out in Ref. (49), where the authors proved it by explicitly constructing the zero mode wavefunctions.
Theorem 5. The magic angle zero mode wavefunctions take the following form 49,52 (up to a normalization factor), where z 0 and z k are the complex coordinates of r 0 , the BA stacking point defined in Eqn. (24), and: The complex coordinate for a vector r is defined as usual: Note here we have written the zero mode wavefunction in terms of the "modified Weierstrass sigma" function σ(z), which is slightly different from Ref. (49), where the authors used Jacobi theta functions. It has been shown 53-56 that both the sigma function and theta function can be used to define the quantum Hall states, and the advantage of the former is modular invariance. The Weierstrass sigma function satisfies a similar quasiperiodic boundary condition as the Jacobi theta function: where a i=1,2 are the complex coordinates of the primitive lattice vectors a 1,2 shown in FIG. (2). The quantum Hall wavefunction and the modified Weierstrass sigma function σ(z) are reviewed in detail in Appendix C. Note that the factor exp(− 1 2 |z k | 2 ) in Eqn. (26) is needed to ensure that the normalization is periodic in k.
The presence of the quasi-periodic elliptic function in the zero mode solution is reminiscent of the lowest Landau level physics on torus 50,57 . We find it conceptually and practically advantageous to rewrite Eqn. (26) in the following form, as a product of a quantum Hall wavefunction and a quasi-periodic spinor wavefunction: where G 1/2 (r) ≡ φ K,b/t (r)/ σ(z − z 0 )e − 1 2 |z| 2 and the quantum Hall wavefunction Φ k is, whose boundary condition can be found in Eqn. (C14) in Appendix C.
Reformulating the zero mode wavefunction in this way makes the subsequent discussions in Section V more clear.

B. Spinor structure of zero mode wavefunctions
The intra-valley inversion implies that the two components of the (magic angle) zero mode wavefunctions are not independent. Theorem 6. The zero mode wavefunction can be written as Eqn. (1), which we copy below, where η = ±1 is the intra-valley inversion eigenvalue from Eqn. (22) and Φ k (r) is the quantum Hall wavefunction Eqn. (31).
The boundary condition of G(r) is derived in Eqn. (C17) of Appendix C.
To conclude, following the intra-valley inversion symmetry, we have derived the spinor structure of the zero mode wavefunction as shown in Eqn. (1) and have demonstrated explicitly its connection to the lowest Landau level wavefunctions. The η in Eqn. (1) is the intravalley inversion eigenvalue, which can be read off from Eqn. (22).

V. NODAL STRUCTURE
In the previous sections, we described an intra-valley inversion symmetry of the chiral model, which led to the discovery of the spinor structure of the zero mode wavefunctions. There we factorized the wavefunction into a quantum Hall wavefunction and a pre-factor G(r).
However, so far the physical interpretation of the function G(r) remains mysterious, as does the structure of zeros in FIG. 1. One hint is that the zero modes must be Bloch functions that transform under the usual translation group, while quantum Hall states transform under the magnetic translation group. Hence G(r) must also be quasi-periodic to "cancel" the magnetic translation effects of the quantum Hall wavefunction.
In this section, we resolve this puzzle by demonstrating mathematically and numerically that G(r) can be regarded as an anti-quantum Hall wavefunction at a certain Landau level, i.e. a quantum Hall state in a magnetic field oppositely directed to that of Φ k , with the order of the magic angle serving the role of the Landau level index. In this way, the zero mode wavefunction is a product of a quantum Hall and an anti-quantum Hall state, whose net magnetic fluxes passing through the moiré unit cell cancel, allowing the whole wavefunction to be a usual Bloch function. We then discuss the zeros in more detail. In the next section, we discuss its experimental implications.

A. Analytical expansion of G(r)
To demonstrate the anti-quantum Hall nature of G(r), we will start by showing that the leading order expansion near r 0 is anti-holomorphic: Hence we can peel off an anti-quantum Hall wavefunction from the zero mode wavefunction and rewrite its components as Eqn. (37) and Eqn. (39).
Since G(r) is independent of Bloch momentum, without loss of generality we can consider the moiré K point sublattice A wavefunction φ K =(φ b K , φ t K ) T to analyze. Its zero mode equation D † (r)φ(r)=0 implies a relation between the top and bottom components φ t K =i∂φ b K /(αU φ ). Theorem 3 tells us that φ t K must have zeros 49 at ±r 0 . From the form of the zero mode wavefunction Eqn. (1), we know that the +r 0 and −r 0 zeros of φ t K come, respectively, from the quantum Hall part Φ k and G(−r). Therefore, near r 0 , φ t K must vanish holomorphically: Then, by using U φ (r 0 ) = 3 and the C 3 symmetry, one can see that φ b K must have a second order zero at r 0 , vanishing as: φ b K (r 0 + r)∼zz. Again according to Eqn. (1), z andz of the bottom component φ b K come, respectively, from the quantum Hall wavefunction and G(r). We hence justified Eqn. (35).

B. Zero mode wavefunction revisited
The vanishing behavior of G(r) near r 0 shows it is possible to factorize out an anti-quantum Hall wavefunction (a quantum Hall state in a magnetic field oppositely directed to that of Φ k , which we denote asΦ k ≡(Φ k ) * ) from it without encountering singularities. The Bloch momentum k ofΦ k is determined by the Bloch translation symmetry of the whole wavefunction. After some algebra, we end up with the final expression: where we introduced a function ρ(r) which must be cellperiodic due to the cancellation of the non-periodic parts from Φ k andΦ k : The top layer wavefunction is obtained easily by the intra-valley inversion symmetry: The Φ k andΦ K of Eqn. (37) carry opposite magnetic fields that cancel with each other, leaving φ b/t k as a Bloch state. Since the crystal momentum (k) dependence, and hence response to an external electric field, is only from the Φ k piece, the wavefunction φ k should have the same topological character as the lowest Landau level wavefunction, according to Laughlin's gauge invariance argument 58 .
To see how this argument applies to our case more explicitly, imagine we apply a time-independent and spatially uniform in-plane external electric field E across the twisted bilayer graphene sample. The Bloch momentum of the electron couples to E through minimal couping, and consequently changes linearly with time: δk ∼ Et. According to Eqn. (27), we know that the zero of Φ k is locked to k, and moves in the direction perpendicular to E. Since zero corresponds to a charge minimum, we conclude that a unit of charge is adiabatically pumped in a direction perpendicular to E during a unit of time. This demonstrates that the zero mode wavefunction Eqn. (1), as a product of a quantum Hall wavefunction and an anti-quantum Hall wavefunction, is indeed a Bloch function which carries Chern number C=1.

C. Zero-structure
We numerically observed that there are multiple zeros occurring at each order of magic angles, as shown in FIG. 1. We now demonstrate they are indeed zeros rather than numerical artifacts.
We noticed that all extra zeros occurring at higher magic angles are located at the reflection symmetric lines. The mirror symmetry M y constrains that ρ(x, y) and ρ * (x, −y) to be the same zero mode solutions. By using the global U (1) phase degree of freedom of the wavefunction, ρ(r) can be chosen to be a purely real function on the reflection symmetric line y=0, the red dot- The zeros of the zero mode wavefunctions are classified into two types by their "movability". One of them is a "movable zero" 59-61 from the quantum Hall wavefunction Φ k , whose location moves linearly with the Bloch wavevector: This zero carries the external Hall response of the Chern band. There are other "frozen zeros" whose locations are fixed and independent of the Bloch momentum k. In particular, among these frozen zeros, one of them is from the anti-quantum Hall state. In FIG. 5, we illustrate the zero-structure, where the black and blue dot represent the movable quantum Hall zero and the frozen anti-quantum Hall zeros. The yellow and red dots are frozen zeros from the function ρ(r). Besides their "movability", zeros are also classified by their "chirality": the wavefunction receives a 2πn phase when the coordinate r encircles the zero once anticlockwise. We numerically noticed that the black and red dots are n=1 zeros, while yellow and blue are n=−1 zeros. Interestingly, as shown in FIG. 1 the center of the unit cells are concentrated with more and more n=−1 zeros at higher magic angles. We discuss the implication for circulating currents in the next section.
We have shown that the zero mode wavefunction shares some similarities with the simple harmonic oscillator system whose eigenstates also have alternating parity and have an increasing number of zeros. In Appendix D, we provide an analytical argument why these features might persist for all higher magic angles by an analogy to the harmonic oscillator 62 . Without loss of generality, we here choose k to be a generic point, K+0.3b1−0.05b2. Since location of the black dot (the zero of the quantum Hall wavefunction Φ k ) is locked with its Bloch momentum according to Eqn. (40), it called a movable zero. The blue dot is the frozen zero from the anti-quantum Hall wavefunction. The yellow and the red dots are the frozen zeros from the function ρ(r). We found the black and red dots are n=1 zeros, while the yellow and blue dots are n=−1 zeros, where n is the 2πn phase that the wavefunction receives when the coordinate r encircles the zero once anticlockwise. The red dotted line is one of the three C3 symmetry related reflection symmetric lines.

VI. EXPERIMENTAL OBSERVATION AND IMPLICATIONS
A. Charge density and scanning tunneling probes One direct consequence of the zeros is a charge density deficiency that can be seen in scanning tunneling spectroscopy experiments 36,63,64 .
We expect a spectroscopy experiment will probe only the top (or bottom) layer, which corresponds to the components φ t (top layer sublattice A wavefunction) and χ t (top layer sublattice B wavefunction). If the spectroscopy measurement has spatial resolution on the level of the atomic spacing, then the sublattice wavefunctions can be probed separately. In this case, fixed zeros in the wavefunction components φ t or χ t correspond to the vanishing of charge density in real space, which will be strongly visible in the spectroscopy experiment. If the ground state is sublattice polarized, which maybe the case on a hexagonal boron nitride substrate, then such spatial resolution is not required to observe the zeros in a spectroscopy experiment. Notice that since the opposite valley wavefunction on the same sublattice is related by T , which acts trivially in real space, we expect the two opposite valley wavefunctions on the same sublattice have the same location of fixed zeros. Therefore, probing the zeros with spectroscopy does not require valley polarization.
The accumulation of zeros at the unit cell center and the unit cell boundary as magic angle order increases, as shown in FIG. 1, should also be visible by spectroscopy with even less atomic resolution. This accumulation will become more prominent at higher magic angles.
Away from the chiral model, the zeros become non-zero minima in the charge density, which we have observed numerically. These will give a less sharp signature in scanning tunneling spectroscopy experiments, but will likely still be observable over some parameter regime.

B. Higher Landau level physics at higher magic angles
As we have seen from Section V, the G(r) piece of the flatband wavefunction in Eqn. (1) has an increasing number of zeros and has an analytical expansion similar to an anti-quantum Hall wavefunction. Consequently, we interpreted the zero mode wavefunction as a product of a higher Landau level anti-quantum Hall state and a lowest Landau level quantum Hall state (Eqn. (37)), where the Landau level index of the former is determined by the order of the magic angle. We also discussed in Section V B that the topological properties of the flatbands are determined by the lowest Landau level quantum Hall piece Φ k since G(r) does not have Bloch momentum dependence.
Nevertheless, we expect the effective interactions projected into the flatbands are modified strongly by both G(r) and Φ k . In particular, we expect the projected interactions would be normalized in a similar way as the normalization effect in higher Landau levels, due to the node-structures in G(r). For example, we expect various kinds of instabilities [65][66][67][68] including charge density waves, bubbles phases, and many other many-body topological phases that are absent at the first magic angle to be possible at flatbands of higher magic angles. Our formulation provides a theoretical and computational pathway towards analyzing interacting physics at different magic angles.

C. Local current and magnetization at higher magic angles
From the charge density of the zero mode wavefunction as plotted in FIG. 1, we observe that for the first magic angle, the charge density maximum occurs at the unit cell center, i.e. the AA stacking point. At higher magic angles, we see an increasing number of zeros appearing at this region. Interestingly, all these zeros are of the same chirality for both layers, while zeros of the opposite chirality are pushed to the boundary of the unit cell. This indicates a stronger phase winding effect and hence circulating currents near the AA stacking region at higher magic angles, which could be experimentally observable nearby the chiral limit.
To see the circulation currents, we first define the following intra-sublattice intra-layer "current operator" J l ss for sublattice s and layer l. The operator J l AA is defined as: operator J l BB (r) is defined in a similar manner but with φ replaced by χ. Here t is the microscopic parameter representing the next-nearest-neighbor hopping strength. We call the above a "current operator" in quot because J l ss is not the current operator of the chiral model, which by definition should be proportional to ∂ k H k , and hence couples distinct sublattices and vanishes within one sublattice. Since the exact flatband wavefunctions are fully sublattice polarized (corresponding physically to hexagonal boron nitride splitting the sublattice degeneracy), the current operator ∂ k H K vanishes within one sublattice polarized state. Nevertheless, we argue that our current operator J l ss has a microscopic origin, and hence should be a physical current operator. The J l ss can be regarded as a continuum version of lattice current i(a † s,i a s,j − a † s,j a s,i ) induced from the next-nearestneighbor hopping process in graphene, where the a s,i are graphene's electron annihilation operators and i, j labels graphene's next-nearest-neighbor sites. In FIG. 6, we plot the real space distribution of J b AA (r), calculated from the bottom layer sublattice-A wavefunction φ b at three different Bloch momenta. According to Ref. (69), the orbital magnetization is dominated by the Gamma point K 0 and Dirac point K in a single valley model, since the bands hybridize most strongly with other bands at these points. Given the strong circulating current present at the second magic angle, it is reasonable to speculate a stronger orbital magnetization 70-76 at higher magic angles than the mag-netization at the first magic angle 69,77,78 for cases close to the chiral limit. Note that the circulation currents are odd under time reversal or a sublattice transformation, hence a strong experimental signal requires valley and sublattice polarization. We leave a detailed exploration of higher magic angle orbital magnetization with more realistic parameters as future work.

VII. CONCLUSION
In this work, we studied the chiral model of twisted bilayer graphene introduced in Ref. (49). We pointed out the intrinsic intra-valley inversion symmetry of the chiral model, protected by the C 2 T crystal symmetry and the linearized Dirac fermion. As a consequence, the energy spectrum is inversion symmetric at all twist angles. Furthermore, zero modes occurring at different magic angles are distinguished by their intra-valley inversion eigenvalue. We numerically found a correspondence of the zero mode inversion parity and the order of magic angles and speculated such an alternating pattern would hold for all magic angles.
We also pointed out the intricate relation between the zero mode wavefunction and the quantum Hall wavefunctions. As guaranteed by intra-valley inversion symmetry, the zero mode wavefunction has an internal spinor structure, and in fact each component can be regarded as a product of a quantum Hall and an anti-quantum Hall wavefunction, which guarantees the zero mode has the periodicity of a Bloch wavefunction. Interestingly, there are an increasing number of zeros occurring in each component at higher magic angles.
In the end, we discussed the implications of our results to realistic systems and observable phenomena. First, these zeros can be detected as charge minima in real space by scanning tunneling spectroscopy. Second, the increasing number of zeros present in the zero mode wavefunction resembles the increasing number of zeros present in higher Landau level wavefunctions. Motivated by this observation, we anticipate higher Landau level physics will occur at the second and higher magic angles. Moreover, we noticed the phase circulation of the flatband wavefunctions at higher magic angles, and anticipate phenomena related to magnetization. We leave more detailed studies on higher magic angles as future work.
Last but not least, it is well known that on a compact manifold, a U (1) magnetic field is subject to a Dirac quantization condition 79 . Our identification of zero mode wavefunctions with two quantum Hall wavefunctions may also shed light on the non-Abelian quantization condition [80][81][82] , where a semi-classical analysis was done recently in Ref. (83).
Note added : during the final stage of the manuscript, we noticed the "unitary particle-hole" symmetry occurring in Ref. (40) which is similar but distinct from our intra-valley inversion symmetry in Eqn. (15); we contrast the difference between the two in the paragraphs under Eqn. (18). We start with setting up the notation of moiré lattice. As mentioned in the main text, we denote the two dimensional lattice vectors as a a=x,y i=1,2 . The area of unit cell is defined to be 2πS: where xy =− yx =1 is the anti-symmetric symbol. The reciprocal basis vectors are: As will be shown in Appendix C, √ S defines an effective magnetic length. We set √ S=1 throughout this work.
Graphene contains A and B sites. As shown in FIG. 2 of the main text, the Dirac points K/K and A/B sites are located at, or, element by element: If we rotate to the chiral basis Ψ c Eqn. (9), the tunneling terms enter in the following way: where the diagonal blocks are: .
(B4) and the off-diagonal blocks are: Using the action of C 2 T in Eqn. (B2), these can be written in terms of only one complex parameter T AB : which satisfy our Lemma 1: where the Pauli matrices τ act on the layer index. In this basis, the chiral matrix σ z that anti-commutes with the Hamiltonian enforces T AA (r)=T BB (r)=0. Notice that Eqn. (B7) also requires linearized Dirac fermion; a quadratic term in the dispersion destroys it. Note that a quadratic term in the dispersion also destroys the exact flatband of the chiral model. We hence demonstrated that for chiral models with linearized Dirac fermion, intra-valley inversion follows from C 2 T symmetry.
Since the lowest Landau level wavefunctions are usually written in terms of holomorphic functions, we start by setting up a notation for complex coordinates. Complex structures ω a=x,y and ω * a=x,y define a one-to-one mapping from two dimensional affine space to the complex plane. We represent the metric and the antisymmetric tensor as They have the properties: ω a =g ab ω b , ω a ω a =0, and ω a ω * a =1. The complex vectors are defined by contracting complex structure with vectors A≡ω a A a , and complex co-vectors as B≡ω a B a .
To distinguish with vectors, complex vectors are unbold. In terms of complex coordinates, the inner product and cross product are respectively A·B≡A a B a =AB * +A * B, A×B≡ ab A a B b =−i(A * B −AB * ). In this work, we took ω x =1 and ω y =i. The quantum Hall system describes two dimensional interacting or noninteracting electrons in a perpendicular magnetic field. In a magnetic field, the electron's coordinate is factorized into the center of its cyclotron motion i.e. guiding center R, and the radius i.e. Landau orbits R: where R commutes withR, but their individual components are noncommutative: (C3) In our case the area of unit cell S plays the same role as magnetic length squared l 2 B = /|eB| where e, B are electron charge and magnetic field strength. When projected into a single Landau level, an electron is fully described by the noncommutative R degrees of freedom. The magnetic translation operator is defined as the following one, which translates the guiding center R by distance d. The magnetic translation algebra is, Due to the single value of wavefunction, any legal wavefunction must transform back to itself after a periodic translation. So we have the boundary condition, where a ∈ A is a lattice vector. From now on we define the whole lattice as A≡{ma 1 +na 2 |m, n ∈ Z}. The phase factor φ a effectively measures the fraction of flux inside the torus. The wavefunctions that satisfy Eqn. (C5) are written in terms of elliptic functions. One choice of elliptic function is the Jacobi theta function 57 . Recently it was also found that the "modified Weierstrass sigma function" is another choice 53,54 . Compared with Jacobi theta function, Weierstrass sigma function has the advantage of being modular invariant.

Modified Weierstrass sigma function
The modified Weierstrass sigma function 53,54 σ(z) is defined as: i.e. a product of the standard Weierstrass sigma functioñ σ(z) and a holomorphic factor e − 1 2Ḡ (A)z 2 , where as will be explained soon the "almost modular form"Ḡ(A) is a modular independent c−number constant that vanishes for square and hexagonal torus. We now introduce thē G(A), and discuss the quasi-periodic property of σ(z).
The standard Weierstrass sigma functionσ(z) has a product series expansion (which is also a fast converging form for numerics), where as defined above, A means the set of lattice points. Clearly, it is modular invariant. It is also quasi-periodic, whereη i is the standard zeta function evaluated at half period, which is related to the k=1 Eisenstein series The Eisenstein series G 2 (a i ) has a highly convergent formula, . (C9) Theη i in addition obey a relation that defines chirality, In our case, the magnetic flux quanta of a unit cell is one, so N φ =1. The (C8) and (C10) suggests a modular independent quantity called "almost modular form", With these formulas in hand, we are ready to get the quasi-periodicity of σ(z): σ(z + a i ) = −e a * i (z+ai/2) σ(z), i = 1, 2.

Quantum Hall wavefunction
The quantum Hall wavefunction is given in Eqn. (31), which we copy below: It has a single zero located at r a k =r a 0 + ab (k − K) b in each unit cell, with r 0 defined in Eqn. (24). Mapping to the complex plane, the zero occurs at z k and its translated counterparts, where z k is: where the first term is zero following from Eqn. (A2) and Eqn. (A3). We used Eqn. (C1) to derive the second term. Since Φ k is not a Bloch function, the "Bloch vector" k here should be understood as labeling the magnetic translation boundary condition Eqn. (C5): t(a 1,2 )Φ k =−e ik·a1,2 Φ k . For a quantum Hall wavefunction, its zero moves linearly with the boundary condition k, reflecting the fact of Chern number C=1 59 . The following diagram FIG. 7 is helpful to quickly figure out r k given the Bloch momentum k.

FIG. 7:
The one-to-one mapping between k and r k , where the first and second letter are r k and k respectively. The figure is constructed by rotating the moiré Brillouin zone by 90 degrees and overlaps with the real space unit cell, precisely because of the mathematical relation r a k =r a 0 + ab (k − K) b . All points in the diagram are illustrated modulo lattice vectors.
Using results derived in the last section, it is easy to find the quasi-periodic boundary condition Φ k satisfies in real and reciprocal space. With i=1,2, they are: The nontrivial phase factors above cannot be removed by a smooth, global, gauge transformation, which reflects the fact that Φ k has a nontrivial Chern number. Technically, these boundary conditions allow one to restrict the discussion to the unit cell and the first Brillouin zone. From now on, we denote k as Bloch momentum inside the first Brillouin zone.
It is straightforward to see how inversion acts on quantum Hall wavefunctions from Eqn. (C13): Φ k (r) = −Φ −k (−r). (C16) We finish this section by showing the boundary condition of G(r), which follows straightforwardly from the periodicity of the zero mode wavefunction in Eqn. (11) and the quantum Hall wavefunction in Eqn. (C14): Appendix D: Analytical Argument for the Nodal Structure The alternating parities and the increasing number of zeros we observed in the chiral model shares many similarities as the simple harmonic oscillator. In this section, we provide an argument for the zero-structure and inversion patterns by making an analogy to simple one-dimensional harmonic oscillators. Specifically, since the additional zeros that occur at higher magic angles occur along a reflection symmetric line, we reduce the zero mode equation to a one variable ordinary differential equation on that line. We can then compare to a harmonic oscillator in one dimension.
The one-dimensional harmonic oscillator is described by the Hamiltonian: whose n th eigenstate φ n (x) satisfies the eigen-equation: where N n =( mω π ) 1 4 (2 n n!) − 1 2 is the normalization factor and H n is the n th Hermite polynomial. Therefore we see that for the harmonic oscillator, the number of zeros of the n-th excited eigenstate is n, and the parity of the n th eigenstate ψ n alternates as (−1) n . Such oscillatory behavior is a generic feature for Sturm-Liouville type differential equations Eqn. (D3) on the interval where p(u) and ω(u) are positive 62 .
We have observed a similar alternating parity and increasing number of zeros of eigenstates at higher magic angles in the chiral model, as discussed in Section III B and Section V C. The problem of the chiral twisted bilayer graphene model is more difficult. One reason is that it is a two-variable differential equation. To make progress, we utilize the symmetry of the problem to reduce the problem to one variable.
We starting by reviewing the zero mode equation, and see how symmetry helps reduce the dimension of the problem. We first recall the zero mode equation from Eqn. (1) and Eqn. (12): −i∂(iG(r)Φ k (r)) = −ηαU φ (r)G(−r)Φ k (r). (D6) By using the lowest Landau level condition that the quantum Hall wavefunction Φ k satisfies, we arrive at the zero mode equation that the function G(r) must satisfy: which is subject to the boundary condition Eqn. (C17). We note that due to the mirror symmetry M y of the problem, both G(x, y) and G * (x, −y) are zero mode solutions of Eqn. (D8). By utilizing the global U (1) phase degree of freedom of wavefunction, one can always choose G(r = 0) to be a purely real number, thereby constraining G(x, 0) to be a real function. We have already used this property for ρ(r) in Section V C, and plotted its real and imaginary part on the reflection symmetric line in FIG. 4. Here we denote the real and imaginary parts of G(x, 0) as R(x) and I(x) respectively. Although the imaginary part vanishes identically at y=0, its y−direction derivative (∂ y I)(x)≡∂ y I(x, y)| y=0 does not. We end up with the following: I(x) = 0, R(x) = 0, (∂ y I)(x) = 0. (D9) The zero mode equation Eqn. (D8) is now rewritten as: subject to the boundary condition Eqn. (C17) which, when reduced to the y=0 line, becomes: where a is the length of the moiré primitive lattice vectors. In the unit S=1 we have been using, its value is a 2 =4π/ √ 3. The derivation so far is exact. The differential equation Eqn. (D10) and its boundary conditions Eqn. (D11)