Quantum valley Hall effect, orbital magnetism, and anomalous Hall effect in twisted multilayer graphene systems

We study the electronic structures and topological properties of $(M+N)$-layer twisted graphene systems. We consider the generic situation that $N$-layer graphene is placed on top of the other $M$-layer graphene, and is twisted with respect to each other by an angle $\theta$. In such twisted multilayer graphene (TMG) systems, we propose that there exists two low-energy bands for each valley emerging from the interface between the $M$ layers and the $N$ layers. These two low-energy bands in the TMG system possess valley Chern numbers that are dependent on both the number of layers and the stacking chiralities, and can be further tuned by an vertical electric field. The valley Chern numbers of the low-energy bands are associated with giant and valley-contrasting orbital magnetizations, which implies that the valley degeneracy can be lifted by weak external magnetic fields or by Coulomb interactions, and a valley-polarized (quantum) anomalous Hall state may be realized.

Twisted bilayer graphene (TBG) has drawn significant attention recently due to the observations of the correlated insulating phases 1-5 and unconventional superconductivity 5,6 . At small twist angles, the lowenergy states of TBG are characterized by four lowenergy bands contributed by the two nearly decoupled monolayer valleys 7,8 . Around the "magic angles", the bandwidths of the four low-energy bands become vanishingly small, and these nearly flat bands are believed to be responsible for the correlated insulating and superconducting properties observed in TBG. Numerous theories have been proposed to understand the electronic structures 9-20 , the correlated insulating phase 9,21-31 , and the mechanism of superconductivity 9,22,24,25,27,[32][33][34][35][36][37] .
On the other hand, the four low-energy bands already exhibits interesting properties at the single-particle level. It has been shown that the four low-energy bands are topologically nontrivial in the sense that they are characterized by odd windings of Wilson loops 13,18,38 , which is an example of the fragile topology 14 . The four flat bands have been further proposed to be equivalent to the zeroth pseudo Landau levels (LLs) with opposite Chern numbers and sublattice polarizations 18 , which is the origin of the nontrivial band topology in the TBG system.
Moreover, recently unconventional ferromagnetic superconductivity and correlated insulating phase have been observed in twisted double bilayer graphene 39,40 . This implies that the low-energy flat bands, which are believed to be responsible for the correlated physics in TBG, may also exist in the twisted double bilayer graphene system. Indeed a recent theoretical study has suggested the presence of flat bands in twisted double bilayer graphene 41 . Motivated by these works, in this paper we study the electronic structures and topological properties of twisted multilayer graphene (TMG). In particular, we consider the general situation that the Nlayer chirally stacked graphene is placed on top of the other M -layer chirally stacked graphene, and they are twisted with respect to each other by a non-vanishing angle θ, as schematically shown in Fig. 1(a) (for the case of M = 2, N = 2). In such a (M + N )-layer TMG sys-tem, we propose that there always exists two low-energy bands (for each valley), and that the bandwidths of the two low-energy bands become vanishingly small at the magic angles of twisted bilayer graphene (TBG) for arbitrary numbers of layers M and N . The flat bands in the TMG system can be interpreted from the pseudo LL representation of TBG 18 , and is protected by an approximate chiral symmetry in chiral graphene multilayers.
Moreover, we also find that there is a Chern-number hierarchy in the (M + N )-layer TMG system. In particular, when the stacking chiralities of the M layers and N layers are the same, the total Chern number of the two low-energy bands for each monolayer valley equals to ±(M − N ) for each spin species 49 . On the other hand, if the stacking chiralities of the M layers and the N layers are opposite, then the total Chern number of the two low-energy bands for each valley is ±(M + N − 2). The valley Chern numbers can be further tuned by an external electric field, leading to gate-tunable quantum valley Hall effect.
The non-vanishing valley Chern numbers of the lowenergy bands are associated with giant and valleycontrasting orbital magnetizations, which implies that the valley degeneracy can be lifted by weak external magnetic fields due to Zeeman couplings, or by Coulomb interactions through spontaneous valley symmetry breaking. In either scenario, a valley-polarized (quantum) anomalous Hall state may be realized. The flat bands at the universal magic angles, together with the Chernnumber hierarchy, make the TMG systems a unique platform to study strongly correlated physics with nontrivial band topology, and may have significant implications on the observed ferromagnetic superconductivity and correlated insulating phase in twisted double bilayer graphene 39,40 .

A. The lattice structures
We consider the most generic case of chirally stacked twisted multilayer graphene, i.e., we place N chiral graphene multilayers on top of M chiral graphene multilayers, and twist them with respect to each other by an angle θ. This is schematically shown in Fig. 1(a) for the case of M = 2, N = 2. Similar to the case of TBG, commensurate moiré supercells are formed when the twist angle θ(m) obeys the condition cos θ(m) = (3m 2 + 3m + 1/2)/(3m 2 + 3m + 1) 42 , where m is a positive integer. The lattice vectors of the moiré superlattice are expressed as t 1 = (− √ 3L s /2, L s /2), and t 2 = (0, L s ), where L s = |t 1 | = a/(2 sin (θ/2)) is the size of the moiré supercell, and a = 2.46Å is the lattice constant of graphene. In TBG it is well known that there are atomic corrugations, i.e., the variation of interlayer distances on the moiré length scale. In particular, in the AB(BA) region of TBG, the interlayer distance d AB ≈ 3.35Å while in the AA-stacked region the interlayer distance d AA ≈ 3.6Å 43 . Such atomic corrugations may be modeled as 11 where g 1 , g 2 , and g 3 = g 1 +g 2 are the three reciprocal lattice vectors of the moire supercell. We take d 0 = 3.433Å and d 1 = 0.0278Å in order to reproduce the interlayer distances in AA-and AB-stacked bilayer graphene. In this paper, the atomic corrugations of the two twisted layers at the interface (between the M layers and the N layers) is also be modeled by Eq. (1). On the other hand, the interlayer distances within the untwisted M layers and the untwisted N layers are set to the interlayer distance of Bernal bilayer graphene d AB = 3.35Å. At a small twist angle θ, the Brillouin zone (BZ) of the moiré supercell has been significantly reduced compared with those of the untwisted multilayers as shown in Fig. 1(b).

B. The effective Hamiltonian
The low-energy effective Hamiltonian of the twisted (M + N )-layer TMG of the K valley is expressed as where H K α (M ) and H + α (N ) are the effective Hamiltonians for the M -layer and N -layer graphene with stacking chiralities α, α = +/−. In particular, where h 0 (k) = − v F (k−K M )·σ stands for the low-energy effective Hamiltonian for monolayer graphene near the Dirac point K M , and h α is the interlayer hopping, with and h − = h † + . The off-diagonal term U represents the coupling between the twisted M layers and N layers. Here we assume that there is only the nearest neighbor interlayer coupling, i.e., the topmost layer of the M -layer graphene is only coupled with the bottom-most layer of the N -layer graphene, thus where the 2×2 matrix U describes the tunneling between the Dirac states of the twisted bilayers 8,11 where r AB = ( √ 3L s /3, 0), u 0 and u 0 denote the intersublattice and intrasublattice interlayer tunneling amplitudes, with u 0 ≈ 0.098 eV, and u 0 ≈ 0.078 eV 11 . u 0 is smaller than u 0 due to the effects of atomic corrugations 11,18 is the shift between the Dirac points of the N layers and the M layers. The phase factor g(r) is defined as . It worth to note that Eq. (2) is the effective Hamiltonian for the K valley. The Hamiltonian for the K valley is readily obtained by applying a time-reversal operation to H K α,α (M + N ).

C. The emergence of two flat bands and the universal magic angles
We continue to study the electronic structures of the (M + N )-layer TMG systems using the effective Hamiltonian given by Eq.  Fig. 2(a)-(d) respectively 50 . Clearly there are two lowenergy flat bands marked by the red lines that are separated from the other bands. The two low-energy bands are almost exactly flat at θ = 1.05 • for all these TMG systems with different layers, indicating that the magic angle of TBG is universal for the TMG systems regardless the number of layers. It turns out that the two flat bands in TMG originates from the twisted bilayer at the interface, and they remain flat even after being coupled with the other graphene layers due to an (approximate) chiral symmetry of Eq. (2). More details about the origin of the flat bands in the TMG systems can be found in Appendix A. It worth to note that in realistic situations there are further neighbor interlayer hoppings in graphene multilayers, which would break the chiral symmetry of the effective Hamiltonian in Eq. (2). The remote interlayer hoppings in graphene multilayers may broaden the flat bands, and we will discus in detail in Sec. III.

II. THE CHERN-NUMBER HIERARCHY AND QUANTUM VALLEY HALL EFFECT
A. The Chern-number hierarchy The flat bands at the universal magic angle make the TMG systems a perfect platform to study the strongly correlated physics. In addition to the flat bands and the universal magic angles, the low-energy bands in the TMG systems also exhibit unusual topological properties with non-vanishing valley Chern numbers. To be specific, when the stacking chiralities of the M layers and the N layers are the same , the total Chern number of the two low-energy bands for each monolayer valley equals to ±(M − N ). On the other hand, if the stacking chiralities of the M layers and the N layers are opposite, then the total Chern number of the two flat bands for each valley equals to ±(M + N − 2). Such a Chern-number hierarchy is more concisely summerized in the following equation where C K α,α (C K α,α ) denotes the total Chern number of the two low-energy flat bands for the K (K ) valley, and the subscripts α, α = ± represent the stacking chiralities of the M layers and N layers. We would like to emphasize that the total Chern number of the two flat bands (per valley per spin) is a more robust quantity than the Chern number of each individual flat band. This is because the former is protected by the energy gaps between the two flat bands and the other high-energy bands, while the latter is crucially dependent on how the gap between the two flat bands is opened up. Eq. (7) has been numerically verified using the effective Hamiltonian of TMG shown in Eq. (2). In particular, in Fig. 3(a) we plot the Wilson-loop eigenvalues (denoted as w(k)) of the (2 + 2) TMG (M = 2, N = 2) at the first magic angle with the same stacking chirality. The blue circles and red diamonds represent the Wilson loops of the K and K valleys respectively. As clearly shown in the figure, for each valley the total Chern number of the two flat bands vanishes. In Fig. 3(b) we still plot the Wilson loops of the (2 + 2) TMG at the first magic angle, but with opposite stacking chiralities. It is clearly seen that for the K valley (blue circles) the two Wilson loops carry the same Chern number −1, giving rise to a total Chern number of −2 for the K valley (+2 for the K valley), which is consistent with Eq. (7). In Fig. 3(c) and (d) we plot the Wilson-loop eigenvalues for the the (2+4) TMG (M = 2, N = 4) at the first magic angle. When the stacking chiralities are the same, the total Chern number of the two flat bands for the K (K ) valley equals to −2 (+2); while if the stacking chiralities are opposite, then the total Chern number of the two flat bands equals to −4 (+4) for the K (K ) valley. Again, this is in perfect agreement with Eq. (7).
In order to understand the origin of the Chern-number hierarchy shown in Eq. (7), we manually tune the strength of the nearest neighbor interlayer hopping t ⊥ , and inspect the evolution of the band structures of the TMG system. In Fig. 4 we show the band structures of (2+2)-TMG (for the K valley) at θ = 1.05 • with different interlayer hopping strengths t ⊥ . When t ⊥ = 0 ( Fig. 4(a)), the (2+2)-TMG system consists of three mutually decoupled subsystems: the TBG at the interface, one graphene monolayer below the TBG, and one graphene monolayer above the TBG. The TBG at the interface would give rise to two nearly flat bands for each valley at the magic angle (marked by red lines in Fig. 4(a)), whereas the bottom and top graphene monolayers contribute two Dirac cones with giant Fermi velocities around the K s and K s point respectively, as shown in Fig. 4(a). If a small t ⊥ is turned on, small gaps would be opened up at the two Dirac points of the bottom and top graphene monolayers, as shown in Fig. 4(b) for t ⊥ = 0.01 eV. As a result, the valence and/or conduction band of each Dirac cone would acquire a Berry phase of ±π due to the pseudospin winding around the Dirac point, and has a Chern number ±1/2. If the two Dirac cones were decoupled from the TBG at the interface, then the Chern numbers of the valence band and the conduction band of each Dirac cone must be opposite. However, if the two Dirac cones were coupled with the flat bands from the TBG at the interface, then the conduction and valence bands of each Dirac cone could carry the same Chern number with a total Chern number of ±1. If the total Chern numbers of the two Dirac cones are the same, say, both are +1, then the two flat bands from the interface TBG must have the opposite Chern number -2, because the total Chern number of all the three subsystems must vanish. On the other hand, if the total Chern numbers of the two Dirac cones are opposite, then the total Chern number of the two flat bands has to vanish. As t ⊥ continues to increase, the gaps of two the two Dirac cones become larger and larger as shown in Fig. 4(c) and (d) for t ⊥ = 0.08 eV and t ⊥ = 0.2 eV. Eventually the Dirac cones from the top and bottom graphene monolayers are pushed to high energies due to the strong interlayer hopping t ⊥ ≈ 0.39 eV. On the other hands, by virtue of the approximate chiral symmetry, the two flat bands originating from the TBG at the interface would remain flat at the magic angles regardless of the magnitude of t ⊥ .
It is straightforward to generalize the above argument to (M + N )-layer TMG with arbitrary number of layers M and N . When t ⊥ = 0, the low-energy band structure of the (M + N )−layer TMG system at the magic angle includes the two flat bands contributed by the interface TBG, the (M − 1) Dirac cones at K s from the graphene monolayers below the interface, and (N − 1) Dirac cones at K s from the graphene monolayers above the interface 51 . As t ⊥ increases, the three subsystems are coupled together, and gaps are opened up at the two Dirac points of the M − 1 and N − 1 graphene monolayers. As a result, the Chern numbers of the valence and/or conduction bands of the M − 1 and N − 1 Dirac cones would be ±(M −1)/2 and ±(N −1)/2. Again, since the two sets of graphene monolayers are coupled with the TBG at the interface, the conduction and valence bands are allowed to have the same Chern numbers for each set of Dirac cones, such that the totl Chern number for the M − 1 (N − 1) graphene monolayers equals to ±(M − 1) (± (N −1)). On the other hands, the Chern number of the two flat bands originating from the interface TBG has to cancel those of the two sets of Dirac cones, which lead to two situations: (i) if the Chern numbers of the two sets of Dirac cones have the same sign and add together, then the total Chern number of the two flat bands equals to ±(M + N − 2) (ii) if the Chern numbers of the two sets of Dirac cones have opposite signs, then the the total Chern number of the two flat bands equals to ±(M − N ) . This concludes our explanation to the Chern number hierarchy given by Eq. (7).

B. Gate tunable quantum valley Hall effect
The valley Chern numbers given by Eq. (7) can be further tuned by applying an vertical electric potential V ⊥ . Taking the case of (2 + 2)-layer TMG as an example, we study the dependence of the valley Chern numbers on the vertical electric potential V ⊥ at θ = 1.05 • . We consider three different fillings: (i) the chemical potential is in the single-particle gap below the two flat bands, with the Chern number of all the filled bands denoted by C K −1 ; (ii) the chemical potential is at the charge neutrality point with one flat band being filled, and the Chern number of the filled bands is denoted by C K 0 ; (iii) the chemical potential is above the two flat bands but below the other high-energy bands, and the Chern number of the filled bands is denoted by C K +1 . The superscript K refers to the K valley.
In Table I we show the dependence of C K −1 , C K 0 , and C K +1 on the vertical electric potential V ⊥ (in units of meV) for (2 + 2)-layer TMG with opposite stacking chiralities at θ = 1.05 • . We see that at small electric potential |V ⊥ | 36 meV, the Chern numbers of the K valley are unchanged. In particular, the Chern number carried by the two flat bands is always C +1 − C −1 = −2. However, when V ⊥ 36 meV, the Chern numbers at the three different fillings become zeros.
In Table II still we show C K −1 , C K 0 , and C K +1 vs. V ⊥ for (2 + 2)-layer TMG at θ = 1.05 • , but with the same stacking chirality. When the vertical potential V ⊥ = 0, both C K +1 and C K −1 are vanishing, and C K 0 is ill defined since the two flat bands touch each other at K s and K s points. However, once a small V ⊥ is applied, the Chern number at the charge neutrality point becomes C K 0 = +2(−2) for negative (positive) V ⊥ . As |V ⊥ | is further increased, C K 0 remains as constant ±2, but there are additional topological transitions between the flat bands and the other high-energy bands such that the Chern numbers below and above the flat bands C K −1 and C ( +1 K) are changed from ±1 to ±2.

III. VERIFICATIONS USING AN EMPIRICAL TIGHT-BINDING MODEL
The flat bands at the first magic angle (Fig. 2) and the Chern-number hierarchy (Eq. (7)) in TMG have been verified using a realistic microscopic tight-binding (TB) model. To be specific, the hopping parameter between two p z orbitals at different carbon sites i and j in the multilayer system is expressed in the Slater-Koster form where V σ = V 0 σ e −(r−dc)/δ0 , and V π = V 0 π e −(r−a0)/δ0 . d = (d x , d y , d z ) is the displacement vector between the two carbon sites. a 0 = a/ √ 3 = 1.42Å, d c = 3.35Å is the interlayer distance in AB-stacked bilayer graphene, and δ 0 = 0.184 a. V 0 σ = 0.48 eV and V 0 π = −2.7 eV. The atomic corrugations at the interface between the M layers and the N layers in (M + N ) TMG are modeled by Eq. (1), and their effects can be taken into account by plugging Eq. (1) into the hopping parameter shown in Eq. (8).
The bandstructures calculated using the Slater-Koster TB model at θ ≈ 1.08 • for the (2 + 2) TMG are shown in Fig. 4. To be specific, the bandstructure of (2+2) TMG with opposite and the same stacking chiralities are shown in Fig. 4(a) and (c) respectively. It is evident that there are four low-energy bands (contributed by the two valleys) which are separated from the other highenergy bands, and the bandwidths are on the order of 10-15 meV, which are greater than those from the continuum model (see Fig. 2(a). This is because in the continuum model only the nearest neighbor interlayer hopping is kept (see Eq. (A3)), which imposes a chiral symmetry to the Hamiltonian of Eq. (2). As argued in Appendix A, the chiral symmetry of Eq. (2) would pin the zeroth pseudo LLs emerging from the twisted bilayer at the interface to zero energy, leading to almost vanishing band- width as shown in Fig. 2. However, such a chiral symmetry is broken in the realistic Slater-Koster TB model, and the bandwidths are expected to be enhanced due to the further neighbor interlayer hoppings. In Fig. 4(b) and (d) we show the Wilson loops of the four low-energy bands of (2+2) TMG at θ ≈ 1.08 • calculated using the Slater-Koster TB model. Fig. 4(b) ((d)) denotes the case with opposite (the same) stacking chiralities with the valley Chern number being ±2 (0). It is evident that the Wilson loops shown in Fig. 4(b) and (d) are topologically equivalent to those calculated by the continuum model as shown in Fig. 3(a)-(b), which provide a strong numerical support to the Chern-number hierarchy given by Eq. (7).

IV. VALLEY-CONTRASTING ORBITAL MAGNETIZATIONS AND VALLEY-POLARIZED STATES
The valley Chern numbers given by Eq. (7) implies opposite orbital magnetizations for the two monolayer valleys K and K . In particular, according to the "modern theory" of orbital magnetizations [44][45][46][47] , the orbital magnetization of the (M + N )-layer TMG for either K or K valley can be expressed as where nk and |u nk are the eigenenergies and (the periodic part of) Bloch eigenstates of a (M + N )-layer TMG Hamiltonian (denoted by H) for either K or K valley, µ is the chemical potential, and H k = e −ik·r He ik·r . Since the K and K valleys are transformed to each other by a time-reversal operation, it is naturally expected that the valley-contrasting Chern numbers shown in Eq. (7) would lead to opposite orbital magnetization for the two valleys. The valley-contrasting orbital magnetizations further suggests that the valley degeneracy may be lifted by external magnetic fields or by Coulomb interactions, leading to a valley-polarized state with possibly quantized anomalous Hall conductivity due to the non-vanishing Chern number. We have exploited this idea using the continuum model given by Eq. (2). To be explicit, we consider the case of (2 + 2)-layer TMG at the first magic angle θ = 1.05 • . We have calculated the orbital magnetizations (M z ) of the two flat bands for the K valley using Eq. (9). The dependence of M z on the chemical potential µ is plotted in Fig. 5(a). The blue circles in Fig. 5(a) represent the situation that the bottom bilayer and the top bilayer have opposite stacking chiralities with the valley Chern number −2 (see Eq. (7)). In this case the magnitude of M z is giant, and has a dip at the charge neutrality point µ = 0 due to the twofold degeneracy at the Dirac points. When µ > 0, M z increases almost linearly with µ due to the non-vanishing Chern number 46 . At the charge neutrality point µ = 0 , our calculations show that M z is as large as −11.375 µ B (for the K valley) per moiré supercell. The orbital magnetizations for the K valley are just opposite to those of the K valley. On the other hand, the red diamonds in Fig. 4(a) denote the case that the bottom bilayer and the top bilayer have the same stacking chirality with vanishing valley Chern number. In this case the orbital magnetization vanishes identically for any chemical potential due to the presence of a PT symmetry (P is inversion).
The giant orbital magnetizations associated with the non-vanishing valley Chern numbers implies that the val-ley degeneracy can be easily lifted by a weak external magnetic field or by Coulomb interactions, leading to a valley-polarized quantum anomalous Hall state. In Fig. 5(b) we plot the Zeeman-splitted bandstructure of (2 + 2) TMG with opposite stacking chiralities (with valley Chern number ±2) at the first magic angle, with an external magnetic field B = 2.5 T pointing along the −ẑ direction (pointing "down"). The chemical potential is set to the charge neutrality point (µ = 0) such that it will not be changed after applying the magnetic field. As clearly shown in Fig. 4(b), both the valley and the spin degrees of freedom are splitted by the weak external magnetic field, but the valley splitting is much greater than the spin splitting due to the giant orbital magnetization. An external field B = 2.5 T is already strong enough to open a gap between the states of the K and K valleys, leading to a quantum anomalous Hall state with Chern number -2 for each spin species. We expect that such a valley-polarized state may also be realized by the spontaneous valley symmetry breaking due to exchange interactions, which we will leave for future study. The same kind of idea can also be applied to other (M + N )-layer TMG with non-vanishing valley Chern numbers.

V. SUMMARY
To summarize,we have studied the electronic structures and topological properties of the (M + N )-layer TMG system. We have proposed that, with the chiral approximation, there always exists two low-energy bands whose bandwidths become minimally small at the magic angle of TBG. We have further shown that the two flat bands in the TMG system are topologically nontrivial, and exhibits a Chern-number hierarchy. In particular, when the stacking chiralities of the M layers and the N layers are opposite, the total Chern number of the two low-energy bands for each valley equals to ±(M + N − 2) (per spin). On the other hand, if the stacking chiralities of the M layers and the N layers are the same, then the total Chern number of the two low-energy bands for each valley is ±(M − N ) (per spin). The non-vanishing valley Chern numbers are associated with giant and valley-contrasting orbital magnetizations along z direction, which implies that the valley degeneracy can be lifted by weak external magnetic fields or by Coulomb interactions, leading to a valleypolarized (quantum) anomalous Hall state. Specifically, for the (2+2)-layer TMG with opposite stacking chiralities, we have shown that when the chemical potential is at the charge neutrality point, an external magnetic field ∼ 2.5 T is already strong enough to lift the valley degeneracy, and may lead to a quantum anomalous Hall state. Our work is a step forward in understanding the electronic properties of twisted multilayer graphene. The universal magic angles and the Chern-number hierarchy proposed in this work make the TMG system a perfect platform to study the interplay between electrons' Coulomb correlations and nontrivial band topology.
Note added : During the preparation of our manuscript, we note a very recently posted work Ref. 48, in which the electronic structures, superconductivity, and correlated insulating phase have been discussed for twisted double bilayer graphene with the same stacking chirality (AB-AB stacking).
Note that the diagonal blocks h 0 (k ± e A) = − v F (k ± eA/ ) are equivalent to the Dirac fermions coupled with opposite pseudo magnetic fields, which would generate pseudo LLs of opposite Chern numbers ±1 18 . In particular, the two zeroth pseudo LLs have opposite sublattice polarizations, thus the intrasublattice coupling term ±3iu 0 in Eq. (A3) vanishes in the subspace of the zeroth pseudo LLs 18 . Therefore, if we re-write Eq. (A3) in a mixed basis consisted of the Dirac states from the first and fourth layers, and the zeroth pseudo LLs from the second and third layers, then H K +,+ (2 + 2) becomes Each element of Eq. (A4) is a 3 × 3 matrix. In particular, and where ηt ⊥ denotes the coupling between the zeroth pseudo LL and the Dirac states of the first (fourth) layer, which can be expressed as some integral over the eigenfunctions of the zeroth pseudo LLs. Note that we have dropped the higher pseudo LLs in Eq. (A4). The Diagonalizations of Eq. (A6) and Eq. (A5) yield two decoupled zero modes. These two zero modes originate from the two zeroth pseudo LLs contributed by the M th and (M + 1)th twisted bilayer (M = 2 for Eq. (A4)), and they stay at zero energy even after being coupled with the other layers due to the chiral symmetry of the effective Hamiltonian of TMG: all the matrix elements in Eq. (2) and Eq. (A1)-Eq. (A4) are intersublattice. As a consequence, if we apply the gauge transformation such that all the basis functions at the B sublattice changes sign, then the total Hamiltonian and eigenenergies would change sign as well. However, the eigenenergies are supposed to be invariant under such a gauge transformation, which thus enforces that both E(k) and −E(k) have to be the eigenenergies of the Hamiltonian. Therefore, a zero mode would stay at zero energy as long as the chiral symmetry is preserved. Similar argument applies to any (M + N )-layer TMG systems with either opposite or the same stacking chiralities. As long as the chiral symmetry is preserved, the zeroth pseudo LLs emerging from the interface between the M layers and N layers would be pinned to zero energy.
On the other hands, it is well known that at the magic angles of TBG the bandwidth of the two low-energy bands for each valley is minimal. From the perspective of the pseudo LLs 18 , it means that at the magic angles, the states within the zeroth pseudo LLs are minimally coupled with each other (and to the higher pseudo LLs), thus they are almost exactly flat. As discussed above, by virtue of the chiral symmetry, these zeroth pseudo LLs that are maximally flat at the magic angles would remain flat even after being coupled with the other layers in the TMG systems. It follows that the magic angles in TBG should be universal in the TMG systems by virtue of the chiral symmetry in Eq. (2).