# Topologically protected flatness in chiral moiré heterostructures

Valentin Crépel  Center for Computational Quantum Physics, Flatiron Institute,
New York, NY 10010, USA  Peize Ding  Nishchhal Verma  Department of Physics,
Columbia University, New York, NY 10027, USA  Nicolas Regnault  Laboratoire de
Physique de l’École normale supérieure, ENS, Université PSL, CNRS, Sorbonne
Université, Université Paris-Diderot, Sorbonne Paris Cité, 75005 Paris, France
Department of Physics, Princeton University, Princeton, NJ 08544, USA  Raquel
Queiroz  Center for Computational Quantum Physics, Flatiron Institute, New
York, NY 10010, USA  Department of Physics, Columbia University, New York, NY
10027, USA

###### Abstract

The observation of delicate correlated phases in twisted heterostructures of
graphene and transition metal dichalcogenides suggests an inherent resilience
of moiré flat bands against certain types of disorder. We investigate the
robustness of moiré flat bands in the chiral limit of the Bistrizer-MacDonald
model, applicable to both platforms in certain limits. We show a drastic
difference between the protection of the first magic angle and higher magic
angles to chiral symmetric disorder such as random higher moiré potential
harmonics arising, for instance, from lattice relaxation. We find that the
first magic angle is topologically protected by a topological index theorem,
similar to the protection of the zeroth Landau level of Dirac fermions, whose
flatness withstands any chiral symmetric perturbation such as non-uniform
magnetic fields. Focusing on the first magic angle of twisted bilayer
graphene, our analysis reveals a hidden non-local constant of motion that
permits the decomposition of the non-abelian gauge field induced by inter-
layer tunnelings into two decoupled abelian ones, underscoring a topological
mechanism for band flatness. Through numerical simulations, we further show
the strikingly different robustness of flat bands across protected (first) and
unprotected (higher) magic angles in the presence of various types of disorder
and identify the scattering processes that are enhanced or suppressed in the
chiral limit. Interestingly, we find that the suppression of disorder
broadening persists away from the chiral limit and is further accentuated by
isolating a single sublattice polarized flat band in energy. Our analysis
suggests the Berry curvature hotspot at the top of the K and K’ valence band
in the transition metal dichalcogenide monolayers is essential for the
stability of its moiré flat bands and their correlated states.

##  I Introduction

The recent observation of the fractional quantum anomalous Hall effect in
twisted monolayers of $\rm MoTe_{2}$ marks the first unequivocal evidence of
topologically ordered states in quantum material without external magnetic
fields [1, 2, 3, 4]. It paves the way for realizing non-abelian states of
matter through superconducting coupling [5, 6, 7, 8, 9, 10, 11, 12, 13, 14],
which are essential for the development of topologically protected quantum
computing [15]. However, topologically ordered states are typically fragile to
perturbations. The original fractional quantum Hall effect in $\rm GaAs$ was
only observed for extremely high mobility samples [16]. In comparison,
transition metal dichalcogenides (TMDs), such as $\rm MoTe_{2}$, are not
necessarily clean monolayers. Their mobility usually lies two orders of
magnitude below that of typical samples of $\rm GaAs$ that host fractional
quantum Hall states [17, 18, 19, 20, 21]. Additionally, strain, lattice
relaxation, and twist angle variations are expected to further disorder the
heterostructures built from these monolayers [22, 23, 24, 25, 26].
Nonetheless, fractional Chern insulators [27, 28, 29] were observed in both
twisted $\rm MoTe_{2}$ homobilayers and twisted bilayer graphene (TBG) at
finite field [30]. This unexpected observation hints at the possibility that
these heterostructures carry some intrinsic protection to a disorder that
enables the realization of such delicate, strongly correlated phases.

In this work, we ask the question: Can the flatness of moiré flat bands be
itself topologically protected? A hint towards an answer follows from the
pioneering work by Aharonov and Casher [31], who showed that a relativistic
Dirac particle in a magnetic field has exactly $\Phi/\Phi_{0}$ zero modes of
definite chirality, where $\Phi$ is the total magnetic flux through the
lattice, and $\Phi_{0}$ the flux quantum. As a consequence, the zeroth Landau
level of a Dirac particle remains perfectly flat even when the magnetic field
is non-uniform. This perfect immunity against certain forms of disorder is
tied to a chiral anomaly and a topological index that protects the net
chirality of the system [32]. The anticommuting chiral symmetry relates states
of opposite energies, allowing the assignment of a proper chirality to all
zero energy states, such as the zeroth Dirac Landau level. These states with
definite chirality must remain strictly at zero energy, even in the presence
of disorder, as long as the chiral symmetry is preserved. The same anomaly
guarantees that the previous theory cannot be realized on a lattice since the
total chirality of fermions on a lattice must vanish [33]. Therefore, a chiral
Landau level must only appear as a surface theory [34] or within a valley
isolated in momentum space. The latter scenario is relevant for the $K$ valley
in both graphene [35] and TMDs [36]. The chiral anomaly has remarkable
consequences when it comes to the effect of disorder in the Landau levels of
these materials: Provided the disorder does not couple the two Dirac valleys
of opposite chirality, their zeroth Landau level is not broadened by disorder,
unlike every other Landau level at a finite energy [37]. A similar result is
expected for Landau levels of quadratic bands, provided their $g-$factor is
exactly two resulting from the Berry curvature at the band edge [38]. No
anomaly or chiral Landau levels are expected for a gas of free electrons under
a magnetic field, such as in $\rm GaAs$.

With an external magnetic field, which is typically homogeneous, the
additional protection due to the chiral anomaly has limited consequences. Only
a slight broadening difference is seen in experiment between the zeroth and
higher Landau levels in graphene [39]. In fact, in Hall systems, the
protection by a chiral anomaly has little physical importance because the
Landau level spread to cyclotron gap ratio can be efficiently reduced by going
to larger magnetic fields. Contrast this with moiré materials, where there is
no external magnetic field but rather a fictitious gauge field that arises
from highly non-uniform interlayer couplings and where the spectral gap around
the band of interest cannot be arbitrarily increased using an external
experimental knob, and the protection provided by the index theorem can play a
crucial role.

In this paper, we explore the protection of flat bands in moiré bilayers
endowed from the chiral anomaly in the chiral limit of the Bistrizer-MacDonald
model of TBG [40, 41], which has also been shown to be an insightful limit for
twisted TMDs upon adding a mass to the Dirac theory of each monolayer [42],
and may also apply to other types of heterostructures [43, 44, 45, 46, 47].
This limit permits a formal map between the flat bands at the magic angle and
Landau levels under a pseudo, non-abelian gauge field $\vec{\mathcal{A}}$.
Naively, we should expect that these states are _not_ protected by the index
theorem, which is only applicable when the gauge field is abelian (Sec. II).
Remarkably, we discover that for large angles, including the first magic
angle, there is a hidden non-local symmetry that allows for the decomposition
of the $SU(2)$ gauge field into two decoupled abelian components $U(1)\times
U(1)$ (Sec. III). Applied specifically to the chiral limits of TBG and twisted
$\rm MoTe_{2}$, this large-angle abelian regime corresponds to twist angles
$\theta\gtrsim 0.7^{\circ}$ and $\theta\gtrsim 2.1^{\circ}$, respectively. In
the decoupled basis, electrons feel a non-uniform and periodic magnetic field
$B^{\prime}$. The matching of the magnetic length introduced by $B^{\prime}$
and the moiré length approximates well the magic condition, as anticipated in
Ref. [48] where it was understood as a consequence of the topological Atiyah-
Singer index theorem [49, 50]. This leads to striking qualitative differences
in the robustness of the different magic angles against chiral preserving
disorder that we numerically investigate in Sec. IV by introducing random
higher harmonics of the moiré interlayer tunneling potential. We note that our
analysis captures the effect of lattice corrugations and strain introduced at
the moiré unit cell edges, which is particularly important at small twist
angles. We find that the first magic angle is mostly undisturbed by higher
harmonics, which makes it possible to consistently define the same magic angle
condition throughout the sample despite small deformations at sub-moiré length
scales. Importantly, this is not true for higher magic angles, suggesting that
it may be impossible to observe higher magic angles in realistic experimental
settings, which agrees with asymptotic calculations [51].

##  II Index theorem

A highly counter-intuitive property of the zero-energy Landau level of Dirac
particles in a magnetic field was simultaneously discovered by Aharonov and
Casher [31] and Dubrovin and Novikov (who considered the square of the Dirac
Hamiltonian) [38]: it remains perfectly flat and extensively degenerate for
arbitrary non-uniform magnetic fields provided the total magnetic flux
threading the system is a constant. This strikingly contrasts the behavior of
all higher Landau levels of the same system and those stemming from quadratic
bands, which in similar settings acquire a dispersion proportional to the non-
uniformity of the magnetic field.

The perfect immunity of the zeroth Landau level to magnetic field variations
has since then been interpreted as a consequence of a topological index
theorem [32] formulated by Atiyah and Singer (AS) for Dirac operators coupled
to arbitrary (multi-dimensional) gauge fields in even dimensions [49, 50]. The
two-dimensional version with matrix valued gauge field
$\vec{\mathcal{A}}=(\mathcal{A}_{x},\mathcal{A}_{y})$ relevant for the present
study states that the operator

$\mathcal{H}_{\rm D}=\begin{bmatrix}0&\mathcal{D}^{\dagger}\\\
\mathcal{D}&0\end{bmatrix},\quad\mathcal{D}=-2i\mu_{0}\partial_{\bar{z}}+\mathcal{A}_{z},$
(1)

with $2\partial_{\bar{z}}=\partial_{x}+i\partial_{y}$,
$\mathcal{A}_{z}=\mathcal{A}_{x}+i\mathcal{A}_{y}$ and $\mu_{0}$ the identity
with the same size as $\mathcal{A}_{z}$, is subject to the topological anomaly

${\rm dim}~{}({\rm Ker}~{}\mathcal{D})-{\rm dim}~{}({\rm
Ker}~{}\mathcal{D}^{\dagger})=\int_{S}\frac{{\rm Tr}\mathcal{B}}{2\pi},$ (2)

where $S$ is the compact surface on which $\mathcal{H}_{\rm D}$ is defined,
and
$\mathcal{B}=\vec{\nabla}\times\vec{\mathcal{A}}+[\mathcal{A}_{x},\mathcal{A}_{y}]$
denotes the curvature corresponding to the gauge field, i.e. the magnetic
field strength. The block structure in Eq. 1 is provided by the anticommuting
chiral symmetry that Dirac Hamiltonian enjoys in even dimensions.

The topological nature of the AS theorem is clear: it only depends on the
integrated magnetic flux but not on the local details and variation of
$\mathcal{B}$. The right-hand side of Eq. 2 is nothing but the Chern number of
the electromagnetic connection, or gauge field, with which the reader may be
more familiar in the context of Berry connection/curvature [52] (to avoid
confusion, we stress that these two Chern numbers are a priori unrelated and
that the “topology” of the AS theorem is distinct from to band topology). The
most direct physical consequence of the AS theorem is that the number of zero
modes of $\mathcal{H}_{\rm D}$ is independent of local variations of the gauge
potential $\vec{\mathcal{A}}$ and of the metric of the Dirac operator. This is
equivalently stated: the density of states at zero energy is topologically
protected against all perturbations of the Hamiltonian anti-commuting with the
chiral operator. In the case of graphene in a magnetic field, these
perturbations physically correspond to strain and lattice corrugations [37,
53].

The AS index also protects the zeroth Landau level of systems with quadratic
dispersion provided they can be written down as the square of
$\mathcal{H}_{\rm D}$ up to an overall sublattice dependent uniform potential
[38, 54]. For particles in a magnetic field, this condition is equivalent to
having a $g-$factor exactly equal to 2. In other contexts, it can lead to
spin-orbit-like corrections to the quadratic Hamiltonian [42].

##  III Application to chiral TBG

In this section, we introduce the chiral model for TBG and prove that, for
sufficiently large angles, it can be analytically decoupled into two
independent Dirac Hamiltonians of the form of Eq. 1 coupled to abelian gauge
fields. The flat bands appearing in this abelian regime correspond to the
first (experimentally observed) magic angle, whose existence can be understood
as a consequence of the AS theorem. We contrast this behavior with higher
(theoretically predicted) magic angles, which we show are not topologically
protected by the index theorem in Eq. 2 due to their nodal structure. This
difference has important consequences regarding the robustness of the flat
bands, investigated in Sec. IV.

###  III.1 Interlayer tunneling as a fictitious gauge field

We consider the chiral limit of the Bistritzer-MacDonald model for TBG [40,
41], which in the basis $\\{(A,t),(A,b),(B,t),(B,b)\\}$, with $A/B$ and $t/b$
the sublattice and layer degrees of freedom respectively described by
$\vec{\sigma}$ and $\vec{\mu}$ Pauli matrices, and in valley-$K$ reads

$\displaystyle\frac{\mathcal{H}}{|v\kappa|}$
$\displaystyle=\begin{bmatrix}0&\mathcal{D}^{\dagger}\\\
\mathcal{D}&0\end{bmatrix},\quad\mathcal{D}=-2i\mu_{0}\partial_{\bar{z}}+\alpha\mathcal{A}_{z},$
(3) $\displaystyle\mathcal{A}_{z}$
$\displaystyle=A_{z}^{(1)}\mu_{1}+A_{z}^{(2)}\mu_{2},\quad\alpha=\frac{w}{|v\kappa|},$
$\displaystyle A_{z}^{(1)}$
$\displaystyle=\\!\\!\sum_{n=0,1,2}\\!\\!\omega^{n}\cos(\kappa_{n}\cdot
r),\quad A_{z}^{(2)}=\\!\\!\sum_{n=0,1,2}\\!\\!\omega^{n}\sin(\kappa_{n}\cdot
r),$

with $\kappa_{n}$ the $2\pi n/3$-rotation of the $\kappa$ corner of the moiré
Brillouin zone aligned with the negative $y$-axis,
$|\kappa|=2|K|\sin(\theta/2)$, $v$ the Fermi velocity of graphene, and where
$\alpha$ measures the strength of $AB$ tunneling $w$ in unit of the typical
moiré kinetic energy $|v\kappa|$. The Hamiltonian in the other valley is
obtained by time-reversal conjugation.

The suggestive form of Eq. 3 makes it clear that $AB$ interlayer tunnelings in
chiral TBG act as a fictitious gauge potential for the two coupled Dirac cones
with chiral symmetry $\sigma_{3}$ [55]. Furthermore, the chiral model Eq. 3 is
known to host perfectly flat bands at zero energy for
$\alpha=0.586,2.221,3.751,\cdots$ [56]. In that context, it is natural to
wonder whether these exact flat bands are topologically protected by the index
theorem of Eq. 2 or not [57, 48]. At first sight, the answer to this question
seems to be negative. Indeed, $\mathcal{A}_{z}$ physically comes from inter-
layer tunneling and is entirely carried by off-diagonal layer Pauli matrices.
As a result of the trace in Eq. 2, the index corresponding to this gauge field
must vanish. In fact, this is a generic feature of SU($N$) gauge fields with
$N>1$, whose generators are all traceless.

Figure 1: The three components of the SU(2) fictitious magnetic field strength
$\mathcal{B}=B_{1}\mu_{1}+B_{2}\mu_{2}+B_{3}\mu_{3}$ induced by inter-layer
tunneling for a small $\alpha=0.2$ (a), the first magic angle $\alpha=0.586$
(b), and the second magic angle $\alpha=3.751$ (c). The different
$\alpha$-scaling of the abelian
$B_{j=1,2}=\alpha\vec{\nabla}\times\vec{A}^{(j)}$ and non-abelian
$B_{3}=2\alpha^{2}\vec{A}^{(1)}\times\vec{A}^{(2)}$ components causes a
transition from a weakly non-abelian regime ($\alpha\lesssim 1$) to a deeply
non-abelian one ($\alpha\gtrsim 1$). The grey hexagon shows the
$\sqrt{3}\times\sqrt{3}$ enlarged moiré unit cell over which the tunneling
terms are periodic.

We notice that the fictitious gauge field is almost abelian for small values
of $\alpha$ (see Fig. 1). This becomes clear if we split the fictitious field
strength $\mathcal{B}=\mathcal{B}_{\rm ab}+\mathcal{B}_{\rm nab}$ into abelian
$\mathcal{B}_{\rm ab}$ and non-abelian $\mathcal{B}_{\rm nab}$ components,
which respectively read

$\begin{split}\mathcal{B}_{\rm ab}&=B_{1}\mu_{1}+B_{2}\mu_{2},\quad
B_{j=1,2}=\alpha\vec{\nabla}\times\vec{A}^{(j)},\\\ \mathcal{B}_{\rm
nab}&=B_{3}\mu_{3},\quad
B_{3}=2\alpha^{2}\vec{A}^{(1)}\times\vec{A}^{(2)}.\end{split}$ (4)

One of the crucial points of this paper is that the abelian and non-abelian
parts of $\mathcal{B}$ do not scale with the same power of $\alpha$. As a
result, the system is almost abelian for small $\alpha$, i.e. for large twist
angles, while it is almost entirely non-abelian for large $\alpha$, i.e. small
angles. We illustrate this difference in Fig. 1, where we plot the magnetic
strength $B_{1,2,3}$ for a small $\alpha=0.2$, and at the two first magic
angles $\alpha=0.586$ and $3.751$.

The weakly non-abelian case, for $\alpha\ll 1$, bears similarity with that of
electrons evolving in an underlying non-colinear spin texture producing an
effective $SU(2)$ gauge field [58, 59]. In the adiabatic limit where the spin
texture slowly varies in space, the system split into two decoupled
components, respectively aligned and anti-aligned with the underlying texture,
that each experience an abelian gauge potential [60, 61, 62]. Similar
adiabatic approximations have been used in different moiré systems to
approximate the inter-layer tunneling as an effective $U(1)$ abelian gauge
field [63, 64]. For a quantitative description of the physics, the domain of
validity of this adiabatic approximation should be identified, and all the
correcting terms should be explicitly included.

The near-abelian gauge field observed for small $\alpha$ in Fig. 1 hints at a
similar decoupling, thereby separating two distinct regimes in the physics of
TBG: an abelian one for $\alpha\lesssim 1$ in which the AS index theorem can
topologically protect a flat band, and a non-abelian one for $\alpha\gtrsim 1$
where no topological protection occurs and flat bands are accidental.

###  III.2 Hidden constant of motion

We now make the distinction between the abelian and non-abelian regime
analytically precise by proving that an invertible transformation can decouple
the Hamiltonian of chiral TBG (Eq. 3) into two independent abelian $U(1)$
Dirac operators, for which the AS index theorem applies, so long as
$\alpha<0.89$.

Our goal is to find the $\alpha$ for which there exists an invertible
transformation $U$, making the Dirac operator abelian or diagonal; in other
words, we want

$\mathcal{D}^{\prime}=U^{-1}\mathcal{D}U=-2i\mu_{0}\partial_{\bar{z}}+A_{z}^{\prime}\mu_{3},$
(5)

for some $A_{z}^{\prime}$. When such a $U$ exists, we notice that
$[\mu_{3},\mathcal{D}^{\prime}]=0$, such that the original Dirac operator
possesses a hidden constant of motion $[U\mu_{3}U^{-1},\mathcal{D}]=0$.
Conversely, if $\mathcal{D}$ has a constant of motion of the form
$\vec{X}\cdot\vec{\mu}$, then the matrix
$U=N^{-1/2}(\vec{X}\cdot\vec{\mu}+x\mu_{3})$, with $N=2x(x+X_{3})$ and
$x^{2}=\vec{X}\cdot\vec{X}$, provides the desired invertible operator for Eq.
5 to hold. Referring to App. A for a detailed derivation, we find that the
abelian gauge field and field strength are then related to $\vec{X}$ by

$A_{z}^{\prime}=A_{x}^{\prime}+iA_{y}^{\prime}=\frac{\alpha(\vec{X}\cdot\vec{\mathcal{A}})}{x+X_{3}},\quad
B^{\prime}=\vec{\nabla}\times\vec{A}^{\prime}.$ (6)

It is, therefore, sufficient to solve for a periodic hidden constant of motion
$\vec{X}$ satisfying $[\vec{X}\cdot\vec{\mu},\mathcal{D}]=0$. We rewrite this
condition as

$\partial_{\bar{z}}\vec{X}=\alpha M\vec{X},\quad
M_{ij}=\varepsilon_{ikj}A_{z}^{(k)},$ (7)

with $\varepsilon_{ijk}$ the Levi-Civita tensor. Note that the matrix $M$ and
the vector $\vec{X}$ have complex entries. The fact that this hidden constant
of motion is complex is imposed by symmetry. If $\vec{X}$ were real, we could
add $(\vec{X}\cdot\mu)\sigma_{3}$ to the Hamiltonian and energetically
separate the two flat bands at the magic angle of chiral TBG while preserving
all the symmetries of the Hamiltonian, in direct contradiction with the
topology of the bands at zero energy [65]. Note also that the antisymmetry of
the Levi-Civita tensor in Eq. 7 implies that $x^{2}$ is uniform, while the
symmetry of the chiral TBG Hamiltonian ensures that it is real; a
demonstration of these facts is given in App. A. Being interested in the small
$\alpha$ limit, we look for solutions as a Taylor series
$\vec{X}=\sum_{n=0}^{\infty}\alpha^{n}\vec{X}^{(n)}$, whose components must
satisfy the recurrence relation $\partial_{\bar{z}}X^{(0)}=0$ and
$\partial_{\bar{z}}X^{(n+1)}=MX^{(n)}$ for $n\geq 0$. Solving this recurrence
relation yields

$X^{(n)}=I^{(n)}(M)C,\quad I(f)=P\left[\int{\rm d}\bar{z}\,(Mf)\right],$ (8)

where $P[f]=f-{\rm mean}(f)$ removes the uniform part of the periodic function
$f$, $I^{(n)}(M)=I\circ I\circ\cdots\circ I\circ I(M)$ denotes $I$ composed
$n$-times (with $I^{(0)}$ the identity), while $C$ is an unconstrained
integration constant. Numerically, we evaluate the matrix multiplication by
$M$ in real space but perform the anti-holomorphic integration in reciprocal
(momentum) space, using fast Fourier transforms to switch from one
representation to another efficiently.

Figure 2: Maximum of the $n$-th coefficient of the $\alpha$-Taylor series of
the constant of motion $\vec{X}$, in logarithmic scale, for different
integration constants $C$. The slopes of the linear fit for $n>10$ (dotted
lines) provide the constant $p_{C}$ that allows to lower bound the radius of
convergence of the Taylor series. These lower bounds are indicated in boxes.

The explicit solution of Eq. 8 allows us to determine for which $\alpha$ the
hidden constant of motion exists by calculating the radius of convergence of
the series that defines $\vec{X}$. To that end, we compute the
$\ell_{\infty}$-norm of the $\vec{X}^{(n)}$ term and find their asymptotic
behavior at large $n$. The results are shown in Fig. 2, where we observe that
the norms asymptotically grow with $n$ as power laws: $\max|\vec{X}^{(n)}|\sim
e^{np_{c}}$ with $p_{C}$ a real coefficient corresponding to the fitted linear
slopes in Fig. 2. The radius of convergence $\alpha_{c}$ of the series
defining $\vec{X}$ is then lower bounded by $e^{-p_{c}}$. Because the uniform
vector $C$ can be arbitrarily chosen, we have repeated the calculation for
$C=[1,0,0]$, $C=[0,1,0]$ and $C=[0,0,1]$. In the two first cases, where the
$C$ is in-plane, we find an identical radius of convergence
$\alpha_{c}^{\parallel}>0.89$; while the third case, where $C$ lies out-of-
plane, gives $\alpha_{c}^{\perp}>0.61$. In practice, any uniform vector $C$
will lead to a constant of motion using Eq. 8, which we can therefore
construct for all $\alpha<\max(\alpha_{c}^{\parallel},\alpha_{c}^{\perp})$.

We have arrived at one of our main results: for $\alpha<0.89$, it is possible
to recast the chiral Hamiltonian of TBG as two decoupled Dirac Hamiltonian
feeling opposite _abelian_ magnetic field $\pm B^{\prime}$ (see Eq. 6). The
effective magnetic field $B^{\prime}$ is represented in Fig. 3a-c for
different values of $\alpha$ using our numerical solution for the constant of
motion $\vec{X}$. Using graphene’s Fermi velocity and typical values of the
inter-layer tunneling and $AB$ sites in TBG, the identified abelian regime
corresponds to twist angles $\theta>0.7^{\circ}$, which only includes the
first magic angle.

###  III.3 Topologically protected magic angle

Within the identified abelian regime ($\alpha<0.89$), we can use the AS
theorem to predict the emergence of topologically protected flat bands. The
only technical subtlety is to identify the minimal compact integration surface
$S$ containing all independent fermionic degrees of freedom of the theory to
avoid double counting the zero modes [48]. Using the translation and $C_{2y}$
symmetries of TBG, we identify this minimal surface with the half moiré unit
cell depicted in Fig. 3c with a dashed triangle.

Figure 3: Magnetic field strength
$B^{\prime}=\vec{\nabla}\times\vec{A}^{\prime}$ after abelianization with
$C=$[0,0,1] for $\alpha=0.25$ (a), $\alpha=0.4$ (b) and $\alpha=0.55$ (c); the
real-space boundaries are the same as in Fig. 1. The minimal integration
surface $S$ containing all independent fermionic degrees of freedom of the
theory is shown as a dotted triangle in (c). d) Flux of $B^{\prime}$ through
$S$ as a function of $\alpha$. The Atiyah-Singer index ensures the existence
of a flat band when this flux $\Phi$ reaches $2\pi$. In our calculations, this
occurs at $\alpha^{*}=0.581$. We used $512$ points along each primitive vector
of the moiré Brillouin zone.

The index theorem (Eq. 2) ensures the existence of a perfectly flat and
topologically protected band when the flux $\Phi$ of $B^{\prime}$ through $S$
reaches an integer multiple of $2\pi$ [48]. In Fig. 3d, we show the
numerically evaluated value of this flux as a function of $\alpha$. We observe
that, although $B^{\prime}$ becomes more and more localized as $\alpha$ grows
(Fig. 3a-c), its maximum value also increases and so does its integrated flux.
The flux $\Phi$ reaches $2\pi$ for $\alpha^{*}=0.581$, which agrees to less
than a percent with the exact value of Ref. [56]. We attribute this difference
to the limited numerical accuracy of our calculations as $\alpha$ approaches
$\alpha_{c}^{\perp}$. Indeed, the number of terms necessary for the Taylor
series to converge increases as $\alpha$ grows, e.g. it reaches $\sim 200$ for
$\alpha=0.5$. At this point, the coefficients admixed when numerically solving
for the recurrence $I^{(n)}(M)$ differ by more than fifteen orders of
magnitudes (which can be seen by extrapolation of Fig. 2 to $n=200$), which
corresponds to our numerical accuracy.

Let us finally mention that the counting of flat bands agrees with the one
expected for twisted bilayer graphene [40, 41]. Indeed, under the decoupling
of Eq. 5, we can split
$\mathcal{D}^{\prime}=\mathcal{D}^{\prime}_{+}\otimes\mathcal{D}^{\prime}_{-}$
into two independent abelian Dirac operators feeling opposite gauge field
$\mathcal{D}^{\prime}_{\pm}=-2i\partial_{\bar{z}}\pm A^{\prime}_{z}$. As a
result, when $\Phi=2\pi$ at $\alpha=\alpha^{*}$, the AS theorem ensures the
existence of two topologically protected flat bands. One is fully localized on
the $A$ sublattice and is spanned by zero modes of $\mathcal{D}^{\prime}_{+}$,
and the other flat band lives on the $B$ sublattice and is spanned by zero
modes of $\mathcal{D}^{\prime\dagger}_{-}$. The two sets of zero modes are
related to one another by the $C_{2y}$ symmetry of the chiral TBG Hamiltonian.

In essence, the invertible transformation $U$ of Eq. 5 allows us to define
layer-spinors that everywhere follow the direction of the fictitious gauge
field $\mathcal{A}_{z}$ created by the $AB$ inter-layer tunneling. This
decoupling generically relies on an adiabatic approximation, whose region of
validity has been analytically identified by explicit computation of the
radius of convergence of the small-$\alpha$ Taylor series defining the layer-
spinors.

###  III.4 No anomaly for higher magic angles

We have explicitly constructed an invertible transformation decomposing the
chiral model of TBG into two independent abelian components for $\alpha<0.89$
by a Taylor expansion in $\alpha$, enabling us to interpret the flat band at
the first magic angle as topologically protected by the AS index theorem. This
construction hints that higher magic angles $\alpha=2.221,3.751,\cdots$ are
not topologically protected as they lie deep in the non-abelian regime (see
Fig. 1). However, it is unclear whether these flat bands are genuinely
accidental, i.e. not protected by an index, or if they are topologically
protected but the abelianization of the theory involves a transformation $U$
that cannot be expressed as a Taylor series in $\alpha$. Here, we show that
the former is true.

Figure 4: Total density in the top layer $\rho_{\rm tot}^{\rm top}(r)$ of the
flat band polarized on the $A$ sublattice at the (a) first $\alpha=0.586$, (b)
second $\alpha=2.221$, and (c) third $\alpha=3.751$ magic angles. The main
difference between (a) and (b-c) is the existence of nodes away from the
origin and corners of the Wigner-Seitz moiré unit cell (shown in black) for
the higher magic angles.

We rely on the known physical difference between the first and the higher
magic angle flat bands of chiral TBG: their nodal structure. The total layer-
polarized density of the higher magic angle flat bands is indeed known to
possess nodes away from the center and corners of the Wigner-Seitz moiré unit
cell ($AA/AB/BA$ points), as shown in Fig. 4 [66]. At these points, the flat-
band eigenvectors are fully layer polarized, and the transformation matrix
$U$, if it existed, would be entirely carried by $\mu_{3}$. At a nodal point
$r_{0}$ different of the center and corners of the moiré Wigner-Seitz cell, we
would then have $\vec{X}(r_{0})\cdot\vec{\mu}=x\mu_{3}$, which we show is
impossible in App. A if we want to preserve a uniform value of $x$ since both
$\vec{A}_{1}$ and $\vec{A}_{2}$ are non-zero at generic point $r_{0}$ of the
moiré unit cell. The presence of nodes in the total layer-polarized density of
the flat band of the smaller magic angles prevents the abelianization of the
theory. It ensures that the higher magic $\alpha$ possesses accidental
(unprotected) flat bands. This lack of robustness deep in the non-abelian
regime of chiral TBG ($\alpha>1$) agrees with recent mathematical
large-$\alpha$ asymptotic results [51] and earlier numerical observations
[23].

##  IV Effect of disorder

In the previous section, we identified two types of flat bands in chiral TBG:
the topologically protected ones of the first magic angle ($\alpha<0.89$) and
the accidental ones appearing at higher magic angles ($\alpha>0.89$). Here, we
study their behaviors in the presence of disorder by building a hierarchy
between the various possible disorder channels and analyzing how the flat band
responds to them. Our results, summarized in Tab. 1, are finally substantiated
by numerical simulations.

channel | example | abelian ($\alpha<0.89$) | non-abelian ($\alpha>0.89$)
---|---|---|---
off-resonant |  $AA/BB$-tunneling | protected by spectral gap | smaller gap (exponentially small in $\alpha$ [67])
chiral | homo-strain | protected by AS index | not protected by an index
resonant, non-chiral $\mu_{1,2}$ | AB-tunneling | protected by localization near $AA$ | weakly protected (less localized)
resonant, non-chiral $\mu_{3}$ | hetero-strain | not protected | weakly protected by nodes near $AA$
Table 1: Effect of the different disorder channels on the flat bands in the
abelian ($\alpha<0.89$) and non-abelian ($\alpha>0.89$) regimes. Off-resonant
terms connect the flat bands to degrees of freedom above the spectral gap,
chiral terms anti-commute with the abelian chiral symmetry of both flat bands,
while non-chiral resonant term scatter particles between the two flat bands
and are assumed to preserve $C_{3}$ on average (Fig. 5). The second and third
columns summarize if the flat bands are protected from the different types of
disorder and, if so, by which physical mechanism: either scattering is
suppressed by a spectral gap, the AS index, or because the wavefunctions are
localized in real space, and the disorder matrix elements are suppressed.

###  IV.1 Hierarchy between disorder channels

As represented in Fig. 5, we split the various types of disorder susceptible
to broaden the flat bands of chiral TBG into three categories: the off-
resonant terms scattering particles from the flat-band to degrees of freedom
above or below the spectral gap, the chiral terms that anticommute with the
abelian chiral symmetries of the two flat bands (they can only be defined for
the first magic angle, see Sec. III), and resonant non-chiral scattering
processes directly coupling the two sublattice polarized flat-bands. In this
paragraph, we identify what physical mechanism, if any, protects the flat
bands against those three types of disorders.

Figure 5: Schematic representation of the different disorder channels
considered, e.g. in Tab. 1, using the band structure of chiral TBG at the
first magic angle along a high-symmetry line of the moiré Brillouin zone. The
two zero-energy sublattice-polarized flat bands have been shifted in energy
for visibility.

#### Off-resonant channels protected by the spectral gap:

Disorder channels of the first kind connect the flat bands to degrees of
freedom away from zero energy. Included in this category is, for instance,
inter-layer tunneling that preserves the sublattice. Suppose the disorder
strength is small compared to the spectral gap separating the flat bands from
the rest of the spectrum. In that case, such disorder only contributes to the
bandwidth of the flat band to second order, thereby granting an effective
protection of the flat band against off-resonant scattering terms. Because the
spectral gap separating the flat bands from the rest of the spectrum
exponentially decays with $\alpha$ for large $\alpha$ [67], we expect that
higher magic angles will be more susceptible to off-resonant disorder than the
abelian flat band of the first magic angle. This is especially significant
because the decay of the spectral gap is found to be extremely steep: the gap
above the flat bands of the three first magic angles
$\alpha=0.586,2.221,3.751$ is respectively equal to $0.29$, $0.059$ and
$0.0027$ in units of the typical moiré kinetic energy $|v\kappa|$ (that itself
decreases linearly with $\theta$). The asymptotic scaling of the spectral gap
for large $\alpha$ was recently used to prove the impossibility of observing
very small magic angles in chiral TBG [51].

#### Topologically protected channels:

Another difference between the first and higher magic angles is the presence
of an index theorem, derived in Sec. III, which has direct consequences on the
robustness of the flat bands against another type of disorder (which can be
both off-resonant and resonant). Indeed, the topological nature of the AS
theorem ensures perfect immunity of the flat bands of the first magic angle
against disorder channels that simultaneously anticommute with the chiral
symmetries of the two abelian Dirac operators $\mathcal{D}^{\prime}_{\pm}$
(Sec. II). Because of the $U$ transformation, these two chiral symmetries are
not exactly equal to the chiral operator $\sigma_{3}$ of the original
Hamiltonian Eq. 3, but to its projection onto the two abelian sectors.
Therefore, the perturbations that leave the two flat bands intact must
simultaneously anticommute with $\sigma_{3}$ and commute with the constant of
motion $\vec{X}\cdot\vec{\mu}$. The only channels of disorder complying with
these constraints are $\sigma_{1/2}\mu_{0}$ and
$\sigma_{1/2}(\vec{X}\cdot\vec{\mu})$. The former represents a homo-strain
that identically couples to both layers. The latter is a fine-tuned
combination of interlayer tunneling and heterostrain. In practice, any
realistic disorder will contain some non-zero weight onto these two
topologically protected channels, further lowering the overall strength of the
bandwidth-generating part of the disorder for the first magic angle.

#### Localization protected channels:

We now turn to resonant scattering channels between the two sublattice-
polarized flat bands and first consider those carried by the off-diagonal
layer Pauli matrices $\mu_{1,2}$, e.g. higher harmonics of the $AB$ interlayer
tunneling. Any inter-layer term of this form preserving $C_{3}$ symmetry must
vanish at the $AA$ stacking point. If the disorder only preserves $C_{3}$ on
average, and if it does not fluctuate too wildly over the moiré scale (i.e. if
the disorder auto-correlation distance is not too small compared to the moiré
lattice constant), it will remain strongly suppressed in the $AA$ region. Now,
the wavefunctions of the first magic angle flat bands in chiral TBG are known
to be strongly peaked around the $AA$ point [68, 69, 70]. Due to additional
nodes near the $AA$ region, higher magic angles have their density localized
closer to the regions where the inter-layer resonant disorder is strong [66],
which can be seen in Fig. 4 (note the logarithmic scale). A direct consequence
of this fact is a lower sensitivity of the first magic angle to inter-layer
resonant disorder, i.e. an effective protection from the localization of the
wavefunction.

#### Unprotected channels:

At last, we consider resonant intra-layer disorder, corresponding, for
instance, to hetero-strain. These disorder channels can have large amplitude
near the $AA$ region if they preserve $C_{3}$ on average. The previous
argument is therefore reversed, and we now expect higher magic angles having
nodes near the $AA$ region to enjoy a lower sensitivity to this disorder. On
the other hand, the first magic angle is expected to be very sensitive to
hetero-strain. This has indeed been observed numerically, with the landscape
of quantum phases stabilized by interactions drastically changing upon
including hetero-strain of very weak amplitude [71, 72, 73, 74].

To summarize (see Tab. 1), the first magic angle possesses an intrinsic
robustness to a broad class of disorder channels thanks to its topological
origin, large gap, and the localization of its wavefunctions. It is, however,
sensitive to hetero-strain, which can lead to substantial changes in the
quantitative physics of TBG [75]. On the other hand, higher magic angles are
expected to strongly broaden in the presence of any disorder except intra-
layer resonant channels.

###  IV.2 Numerical simulations

To solidify these observations, we numerically simulate the effect of disorder
in the chiral model of TBG. The role of the AS index in the robustness against
disorder can be further highlighted by focusing on disorder in the $AB$ inter-
layer tunneling channel, for which all magic angles possess some protection
due to the localization of their wavefunctions and where the first magic angle
enjoys an additional protection due to its abelian nature (third line of Tab.
1).

To be more concrete, we introduce a $C_{3}$ symmetric disorder by adding
higher-order harmonics to the inter-layer tunneling potential. This amounts to
the substitution $\mathcal{A}_{z}\to\mathcal{A}_{z}+\delta\mathcal{A}_{z}$ in
Eq. 3 with

$\displaystyle\delta\mathcal{A}_{z}=\sum_{m=1}^{m_{c}}\mu_{m}(a_{m,z}^{(1)}\mu_{1}+a_{m,z}^{(2)}\mu_{2}),$
(9) $\displaystyle a_{m,z}^{(1)}+ia_{m,z}^{(2)}=\sum_{n}e^{i(\kappa_{m,n}\cdot
r+\phi_{m,n})},$

where $\kappa_{m,0}$ is the $m$-th smallest vector in the first trine of the
plane that connects the two corners $\kappa$ and $\kappa^{\prime}$ of the
moiré Brillouin zone (up to reciprocal lattice vectors), $\kappa_{m,n}$
denotes its rotation by $2n\pi/3$, and $\mu_{m}$ stands for the amplitude of
the $m$-th harmonic. The $\phi_{m,n}$ phases are included such that
$\delta\mathcal{A}_{z}$ does not break $C_{3}$ and the $C_{2y}$ symmetry of
Eq. 3. The precise definition of all these terms is given in App. B. We only
consider the first $m_{c}=10$ harmonics and discretize $\mathcal{D}$ on a
momentum-space lattice containing all plane waves with wavevectors of norm
smaller than a certain radius $R_{k}$.

For the considered $C_{3}$-preserving inter-layer tunneling disorder, the
chiral model of TBG still possesses exact magic angles where the bands at zero
energy become perfectly degenerate [76, 56]. The precise value of these magic
$\alpha$ can be efficiently determined using the method of Refs. [56, 67],
where all magic angles are determined as eigenvalues of a carefully designed
operator $\Lambda$, whose explicit expression is given in App. B for
comprehensiveness. We first use this method to characterize the sensitivity of
the magic angles to the different harmonics in Eq. 9. For this purpose, we
choose one specific harmonic $m$ and compute the change $\delta\alpha$ in the
first, second, and third magic angles when $\mu_{m}=\delta\mu$ is set to a
small but finite value (keeping all other $\mu_{n\neq m}=0$). For small
disorder strength, this corresponds to taking the derivatives
$\partial\alpha/\partial\mu$ of $\Lambda$’s eigenvalues with respect to
$\mu_{m}$. We numerically evaluate this derivative using a small enough
$\delta\mu=0.02$, and plot our results for $|\partial\alpha/\partial\mu|$ in
Fig. 6a. We observe that the $\alpha$ deviations in the first magic angle
exponentially decay with $m$ and are negligible for $m>5$, indicating that the
first magic angle is robust against microscopic details on scales smaller than
the moiré length. In contrast, the deviations for the higher magic angle can
reach up to ten times the amplitude of the added perturbation and do not decay
with $m$. This suggests that the higher magic angles are easily affected by
microscopic disorder at all scales down to graphene’s lattice constant.
Perturbations on those scales, such as lattice relaxations, are inevitable and
necessarily spoil the prediction for higher magic angles [41, 56].
Perturbations on scales larger than the moiré length are also studied in App.
C.

Figure 6: a) Variation $|\partial\alpha/\partial\mu|$ of the first (blue),
second (orange), and third (green) magic angle computed for a small
$\delta\mu=0.02$ perturbation in the $m$-th harmonics of the $AB$ inter-layer
tunneling (see Eq. 9). To obtain converged results for large $m$, we used a
large truncation radius $R_{k}=25\sqrt{3}|\kappa|$. b) Histogram of the magic
angles for $500$ realization of the inter-layer tunneling disorder with all
$\mu_{m}$ drawn uniformly and normalized as
$\sqrt{\sum_{m}\mu_{m}^{2}}=\delta\mu=0.1$. c) $E_{\gamma}$, the energy of the
top active band crossing zero energy in chiral TBG at the center of the moiré
Brillouin zone ($\gamma$), averaged over $500$ disorder realizations for the
same distribution of $\mu_{m}$ as in (b). The dashed line shows the result in
the absence of disorder. $R_{k}=25\sqrt{3}|\kappa|$ for (b) and
$R_{k}=10\sqrt{3}|\kappa|$ for (c) is used. Colors are provided as a guide to
the eye, allowing us to identify the regions closest to 0.586 (blue), 2.221
(orange), and 3.751 (green), which correspond to the magic angles in the
absence of disorder.

To understand the physical implications of the difference in $\delta\alpha$
for first and higher magic angles observed in Fig. 6a, we consider realistic
disorder realization containing all harmonics $m\leq m_{c}$ with amplitudes
$\mu_{m}$ drawn uniformly in $[-1,1]$ with the overall normalization
$\sqrt{\sum_{m}\mu_{m}^{2}}=\delta\mu=0.1$. We plot the distribution of the
first three magic angles for 500 realizations of such disorder in Fig. 6b,
which shows a similar trend as Fig. 6a: the distribution of magic $\alpha$ is
sharply peaked around the original value 0.586 for the first magic angle. At
the same time, it is extremely broad and featureless for higher magic angles.
To interpret this result, picture a sample in which different patches of the
sample (different unit cells) are characterized by slightly different atomic
arrangements and lattice relaxation profiles that locally change the values of
the $\mu_{m}$. We can assign an effective magic angle $\alpha$ to each of
these regions. When TBG is close to the first unperturbed magic angle, the
distribution of these magic $\alpha$ is peaked, and the magic condition is
almost simultaneously satisfied for all patches (see Fig. 6b). This is not
true for higher magic angles, such that the density of states at zero energy
is smeared in many patches, yielding a strong broadening of the original flat
bands.

We can estimate the flat band broadening by computing $E_{\gamma}$ the average
energy of the topmost band at the center of the moiré Brillouin zone ($\gamma$
point) for 500 realizations of the same realistic disorder pattern for a fixed
twist angle, i.e. as a function $\alpha$. The results presented in Fig. 6c
pictorially capture our previous statement: the higher magic angles are
completely washed out due to their high sensitivity to disorder. The average
$\gamma$-point energy only displays an apparent non-monotonic behavior with an
almost vanishing bandwidth at the first magic angle, which we have seen to be
protected against disorder (Tab. 1).

Our numerical calculations clearly show that only the flat bands at the first
magic angle retain a large density of states in the presence of disorder. As
explained in Tab. 1, this is due to the combination of the topological
protection from the AS theorem, large spectral gap, and strong localization of
the wavefunctions in the flat band. Higher magic angles are sensitive to
lattice scale details and should be interpreted as accidental artifacts of the
continuum model rather than experimentally observable features.

###  IV.3 Away from the chiral limit

We finally check whether our results, which find their roots in the chiral
limit of TBG, also hold in the more realistic situation where non-zero $AA$
tunnelings are included in the Bistritzer-MacDonald model [40]. When the $AA$
tunneling amplitude $w_{0}$ is non-zero, we cannot identify the magic angle
with the point where the $\gamma$-point energy vanishes, as we did in Fig. 6c,
since the model never features perfectly flat bands. Instead, we locate the
magic angle as the points for which the Fermi velocity at the corner of the
moiré Brillouin zone vanishes [76, 40]. We have numerically computed this
Fermi velocity $v_{m}(\alpha)$ averaged over 500 realizations of the generic
disorder pattern described above as a function of $\alpha$, both in the chiral
limit $w_{0}=0$ and for the realistic value $w_{0}=0.8w$ [40]. Our results are
presented in Fig. 7, where, for better comparison, we have chosen
$\delta\mu=0.1$ for the chiral limit and $\delta\mu=0.2$ in the realistic
regime, leading to the similar average velocity at the first magic angle. In
both the chiral and realistic (nonchiral) limits, disorder barely affects the
Fermi velocity (blue lines) from its clean values (dotted lines) as long as we
are in the abelian regime ($\alpha<0.86$). In the non-abelian regime, the
average velocity deviations become significant. Consequently, only the first
magic angle still strongly suppresses the Fermi velocity $v_{m}(\alpha)$. In
contrast, all higher magic angles are washed out by disorder – we interpret
this as a direct consequence of the robustness granted by the AS topological
index described in the second line of Tab. 1.

To understand why a stronger disorder strength is necessary for the realistic
limit ($\delta\mu=0.2$) to obtain the same average value of $v_{m}(\alpha)$ at
the first magic angle than in the chiral limit ($\delta\mu=0.1$), we should
take a closer look at the effect of the non-chiral part of the disorder. More
precisely, it is known that the wavefunction of the flat bands is more
localized near the $AA$ region in the realistic limit $w_{0}=0.8w$ than in the
chiral limit [41]. Therefore, the third line of Tab. 1 qualitatively captures
the better protection of the flat band against inter-layer tunneling disorder.

Figure 7: Fermi velocity at the corner of the moiré Brillouin zone
$v_{m}(\alpha)$ averaged over $500$ disorder realizations for the same
distribution of $\mu_{m}$ as in Fig. 6, both in the chiral limit (a) and the
realistic value for $AA$ tunneling strength $w_{0}=0.8w$ (b). To obtain
similar results at the first magic angle, allowing for a fair comparison, we
have used $\delta\mu=0.1$ for the chiral limit (a) and $\delta\mu=0.2$ in the
realistic regime (b). The dashed lines show the same quantity in the absence
of disorder. We have normalized the results by the Fermi velocity of graphene
$v$ and shown them on a log scale to better discern the effect of the weak
disorder.

##  V Implications for TMDs

While TBG has been the main focus up to now, let us expand on the implications
of our results to twisted TMDs, which share the same chiral limit as TBG (Eq.
3) up to the addition of a mass term that describes the gap of the
semiconducting monolayers [42]. Most of our analysis, particularly the
decoupling into abelian components from Sec. III, goes through in the presence
of this mass term. For twisted $\rm MoTe_{2}$ homobilayers, similar estimates
as above show that the physics is effectively abelian for $\theta>2^{\circ}$,
which lead to topologically protected flat bands at the first magic angles
$\theta\sim 3.2^{\circ}$ [42]. Away from the chiral limit, the resulting flat
band will acquire some dispersion due to additional orbital-preserving terms
in the Hamiltonian but remain robust against disorder if the results of Fig. 7
can be generalized to twisted $\rm MoTe_{2}$ homobilayers. The mass term and
the strong spin-orbit coupling of TMDs [36] strongly suppress all resonant
non-chiral disorder – the only ones that spoil the protection of the flat
bands, see Tab 1 – as they fully lift the degeneracy of the flat band to leave
a single topologically protected valence flat band per spin flavor. Twisted
TMDs, therefore, are ideal hosts for delicate correlated topological phases of
matter. They host weakly dispersing Chern bands fully protected against all
forms of disorder thanks to the combinations of the AS topological index
theorem and their spectral gap.

We also note that other approaches aim at interpreting the physics of twisted
TMDs as particles moving in an abelian gauge field [63, 64, 77, 78, 79, 80].
Our work precisely determines in which regimes and for which angles such a
picture is valid.

##  VI Conclusion

In this work, we have demonstrated that the flat band appearing at the first
magic angle in twisted graphene and transition metal dichalcogenides bilayers,
both in the chiral limit, can be understood as originating from Dirac
particles evolving in an abelian gauge field. Their zero energy Landau levels
become perfectly flat and topologically protected against disorder thanks to a
topological index theorem when the moiré length and the effective magnetic
length of the system match. This description is only possible when the inter-
layer tunneling acting as a fictitious gauge field can be adiabatically
decoupled into abelian components. Such decoupling occurs for large twist
angles where the typical kinetic energy on the moiré scale is not negligible
compared to the strength of tunneling and corresponds to $\theta>0.7^{\circ}$
for twisted bilayer graphene and $\theta>2^{\circ}$ for twisted $\rm
MoTe_{2}$. The topological nature of the zeroth Landau level of chiral Dirac
operators has critical consequences on the robustness of the magic angle
against perturbations. In particular, we find that only the first magic angle
benefits from the enhanced robustness to disorder granted by the Atiyah-Singer
topological index. Therefore, we conclude that flat bands appearing at smaller
twist angles are accidental, in the sense that they do not have an index and
are strongly susceptible to spreading out in energy even for weak disorder on
scales smaller than the moiré pattern. Let us finally highlight that all the
flat bands at all magic angles are characterized by the same band topology
(Chern number, Wannier obstructions, etc…) [65, 81], but only for those
appearing at the first magic angle, the full Hamiltonian of TBG features an
additional Atiyah-Singer topological index protecting the flat band.

##  VII Acknowledgement

Research on topological properties of moiré materials is supported as part of
Programmable Quantum Materials, an Energy Frontier Research Center funded by
the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences
(BES), under award DE-SC0019443. The Flatiron Institute is a division of the
Simons Foundation. N.R. acknowledges support from the QuantERA II Programme
which has received funding from the European Union’s Horizon 2020 research and
innovation program under Grant Agreement No 101017733.

## References

  * Cai  _et al._ [2023] J. Cai, E. Anderson, C. Wang, X. Zhang, X. Liu, W. Holtzmann, Y. Zhang, F. Fan, T. Taniguchi, K. Watanabe, _et al._ , Nature 622, 63 (2023).
  * Xu  _et al._ [2023] F. Xu, Z. Sun, T. Jia, C. Liu, C. Xu, C. Li, Y. Gu, K. Watanabe, T. Taniguchi, B. Tong, J. Jia, Z. Shi, S. Jiang, Y. Zhang, X. Liu, and T. Li, Physical Review X 13, 031037 (2023).
  * Zeng  _et al._ [2023] Y. Zeng, Z. Xia, K. Kang, J. Zhu, P. Knüppel, C. Vaswani, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, Nature 622, 69 (2023).
  * Park  _et al._ [2023] H. Park, J. Cai, E. Anderson, Y. Zhang, J. Zhu, X. Liu, C. Wang, W. Holtzmann, C. Hu, Z. Liu, _et al._ , Nature 622, 74 (2023).
  * Mong  _et al._ [2014] R. S. Mong, D. J. Clarke, J. Alicea, N. H. Lindner, P. Fendley, C. Nayak, Y. Oreg, A. Stern, E. Berg, K. Shtengel, _et al._ , Physical Review X 4, 011036 (2014).
  * Vaezi [2014] A. Vaezi, Physical Review X 4, 031009 (2014).
  * Repellin  _et al._ [2018] C. Repellin, A. M. Cook, T. Neupert, and N. Regnault, npj Quantum Materials 3, 14 (2018).
  * Barkeshli [2016] M. Barkeshli, Physical Review Letters 117, 096803 (2016).
  * Barkeshli and Wen [2011] M. Barkeshli and X.-G. Wen, Physical Review B 84, 10.1103/physrevb.84.115121 (2011).
  * Barkeshli and Qi [2012] M. Barkeshli and X.-L. Qi, Physical Review X 2, 031013 (2012).
  * Crépel  _et al._ [2019] V. Crépel, N. Regnault, and B. Estienne, Physical Review B 100, 125128 (2019).
  * Katzir  _et al._ [2020] B. A. Katzir, A. Stern, E. Berg, and N. H. Lindner, arXiv preprint arXiv:2011.13950 (2020).
  * Reddy  _et al._ [2024] A. P. Reddy, N. Paul, A. Abouelkomsan, and L. Fu, arXiv preprint arXiv:2403.00059 (2024).
  * Crépel and Regnault [2024] V. Crépel and N. Regnault, arXiv preprint arXiv:2403.05622 (2024).
  * Kitaev [2003] A. Y. Kitaev, Annals of physics 303, 2 (2003).
  * Stormer  _et al._ [1999] H. L. Stormer, D. C. Tsui, and A. C. Gossard, Reviews of Modern Physics 71, S298 (1999).
  * Movva  _et al._ [2015] H. C. Movva, A. Rai, S. Kang, K. Kim, B. Fallahazad, T. Taniguchi, K. Watanabe, E. Tutuc, and S. K. Banerjee, ACS nano 9, 10402 (2015).
  * Chuang  _et al._ [2016] H.-J. Chuang, B. Chamlagain, M. Koehler, M. M. Perera, J. Yan, D. Mandrus, D. Tomanek, and Z. Zhou, Nano letters 16, 1896 (2016).
  * Xu  _et al._ [2017] S. Xu, J. Shen, G. Long, Z. Wu, Z.-q. Bao, C.-C. Liu, X. Xiao, T. Han, J. Lin, Y. Wu, _et al._ , Physical Review Letters 118, 067702 (2017).
  * Pisoni  _et al._ [2019] R. Pisoni, T. Davatz, K. Watanabe, T. Taniguchi, T. Ihn, and K. Ensslin, Physical Review Letters 123, 117702 (2019).
  * Zhao  _et al._ [2023] S. Zhao, J. Huang, V. Crépel, X. Wu, T. Zhang, H. Wang, X. Han, Z. Li, C. Xi, S. Pan, _et al._ , arXiv preprint arXiv:2308.02821 (2023).
  * Wilson  _et al._ [2020] J. H. Wilson, Y. Fu, S. Das Sarma, and J. H. Pixley, Physical Review Research 2, 023325 (2020).
  * Carr  _et al._ [2019] S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, Physical Review Research 1, 013001 (2019).
  * Nakatsuji and Koshino [2022] N. Nakatsuji and M. Koshino, Physical Review B 105, 245408 (2022).
  * Sainz-Cruz  _et al._ [2021] H. Sainz-Cruz, T. Cea, P. A. Pantaleón, and F. Guinea, Physical Review B 104, 075144 (2021).
  * Thomson and Alicea [2021] A. Thomson and J. Alicea, Physical Review B 103, 125138 (2021).
  * Sheng  _et al._ [2011] D. Sheng, Z.-C. Gu, K. Sun, and L. Sheng, Nature communications 2, 389 (2011).
  * Neupert  _et al._ [2011] T. Neupert, L. Santos, C. Chamon, and C. Mudry, Physical Review Letters 106, 236804 (2011).
  * Regnault and Bernevig [2011] N. Regnault and B. A. Bernevig, Physical Review X 1, 021014 (2011).
  * Xie  _et al._ [2021a] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, _et al._ , Nature 600, 439 (2021a).
  * Aharonov and Casher [1979] Y. Aharonov and A. Casher, Physical Review A 19, 2461 (1979).
  * Katsnelson and Novoselov [2007] M. Katsnelson and K. Novoselov, Solid State Communications 143, 3 (2007).
  * Nielsen and Ninomiya [1981] H. Nielsen and M. Ninomiya, Physics Letters B 105, 219 (1981).
  * Zhang  _et al._ [2012] F. Zhang, C. L. Kane, and E. J. Mele, Physical Review B 86, 081303 (2012).
  * Novoselov  _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. Dubonos, and A. Firsov, nature 438, 197 (2005).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Physical Review Letters 108, 196802 (2012).
  * Kawarabayashi  _et al._ [2009] T. Kawarabayashi, Y. Hatsugai, and H. Aoki, Physical Review Letters 103, 156804 (2009).
  * Dubrovin  _et al._ [1980] B. Dubrovin, S. Novikov, _et al._ , Pisma V Zhurnal Eksperimentalnoi I Teoreticheskoi Fiziki 79, 1006 (1980).
  * Giesbers  _et al._ [2007] A. Giesbers, U. Zeitler, M. Katsnelson, L. Ponomarenko, T. Mohiuddin, and J. Maan, Physical Review Letters 99, 206803 (2007).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Physical Review Letters 122, 106405 (2019).
  * Crépel  _et al._ [2023a] V. Crépel, N. Regnault, and R. Queiroz, arXiv preprint arXiv:2305.10477 (2023a).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Physical Review B 100, 085109 (2019).
  * Guerci  _et al._ [2022] D. Guerci, P. Simon, and C. Mora, Physical Review Research 4, L012013 (2022).
  * Wan  _et al._ [2023] X. Wan, S. Sarkar, S.-Z. Lin, and K. Sun, Physical Review Letters 130, 216401 (2023).
  * Wang and Liu [2022] J. Wang and Z. Liu, Physical Review Letters 128, 176403 (2022).
  * Crépel  _et al._ [2023b] V. Crépel, A. Dunbrack, D. Guerci, J. Bonini, and J. Cano, Physical Review B 108, 075126 (2023b).
  * Parhizkar and Galitski [2023] A. Parhizkar and V. Galitski, arXiv preprint arXiv:2301.00824 (2023).
  * Atiyah and Singer [1968] M. F. Atiyah and I. M. Singer, Annals of mathematics , 484 (1968).
  * Atiyah and Singer [1984] M. Atiyah and I. M. Singer, Proceedings of the National Academy of Sciences 81, 2597 (1984).
  * Becker  _et al._ [2024] S. Becker, I. Oltman, and M. Vogel, arXiv preprint arXiv:2402.12799 (2024).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Physical Review Letters 49, 405 (1982).
  * Katsnelson and Prokhorova [2008] M. Katsnelson and M. Prokhorova, Physical Review B 77, 205424 (2008).
  * Estienne  _et al._ [2023] B. Estienne, N. Regnault, and V. Crépel, Physical Review Research 5, L032048 (2023).
  * Liu  _et al._ [2019] J. Liu, J. Liu, and X. Dai, Physical Review B 99, 155415 (2019).
  * Becker  _et al._ [2022] S. Becker, M. Embree, J. Wittsten, and M. Zworski, Probability and Mathematical Physics 3, 69 (2022).
  * Sheffer and Stern [2021] Y. Sheffer and A. Stern, Physical Review B 104, L121405 (2021).
  * Ye  _et al._ [1999] J. Ye, Y. B. Kim, A. J. Millis, B. I. Shraiman, P. Majumdar, and Z. Tešanović, Physical Review Letters 83, 3737 (1999).
  * Bruno  _et al._ [2004] P. Bruno, V. K. Dugaev, and M. Taillefumier, Physical Review Letters 93, 096806 (2004).
  * Martin and Batista [2008] I. Martin and C. Batista, Physical Review Letters 101, 156402 (2008).
  * Nakatsuji  _et al._ [2015] S. Nakatsuji, N. Kiyohara, and T. Higo, Nature 527, 212 (2015).
  * Verma  _et al._ [2022] N. Verma, Z. Addison, and M. Randeria, Science Advances 8, eabq2765 (2022).
  * Wu  _et al._ [2019] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Physical Review Letters 122, 086402 (2019).
  * Morales-Durán  _et al._ [2023] N. Morales-Durán, N. Wei, and A. H. MacDonald, arXiv preprint arXiv:2308.03143 (2023).
  * Xie  _et al._ [2021b] F. Xie, N. Regnault, D. Călugăru, B. A. Bernevig, and B. Lian, Physical Review B 104, 115167 (2021b).
  * Wang  _et al._ [2021] J. Wang, Y. Zheng, A. J. Millis, and J. Cano, Physical Review Research 3, 023155 (2021).
  * Becker  _et al._ [2023] S. Becker, I. Oltman, and M. Vogel, arXiv preprint arXiv:2309.02701 (2023).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Physical Review X 8, 031087 (2018).
  * Zang  _et al._ [2022] J. Zang, J. Wang, A. Georges, J. Cano, and A. J. Millis, arXiv preprint arXiv:2210.11573 (2022).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Physical Review Letters 129, 047601 (2022).
  * Zhang  _et al._ [2019] Y.-H. Zhang, H. C. Po, and T. Senthil, Physical Review B 100, 125104 (2019).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Physical Review X 11, 041063 (2021).
  * Parker  _et al._ [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, N. Bultinck, _et al._ , Physical Review Letters 127, 027601 (2021).
  * Wagner  _et al._ [2022] G. Wagner, Y. H. Kwan, N. Bultinck, S. H. Simon, and S. Parameswaran, Physical Review Letters 128, 156401 (2022).
  * Liu  _et al._ [2022] X. Liu, G. Farahi, C.-L. Chiu, Z. Papic, K. Watanabe, T. Taniguchi, M. P. Zaletel, and A. Yazdani, Science 375, 321 (2022).
  * Sheffer  _et al._ [2023] Y. Sheffer, R. Queiroz, and A. Stern, Physical Review X 13, 021012 (2023).
  * Crépel and Millis [2024] V. Crépel and A. Millis, arXiv preprint arXiv:2403.15546 (2024).
  * Li  _et al._ [2024] H. Li, Y. Su, Y. B. Kim, H.-Y. Kee, K. Sun, and S.-Z. Lin, arXiv preprint arXiv:2402.02251 (2024).
  * Kang and Vafek [2023] J. Kang and O. Vafek, Physical Review B 107, 075408 (2023).
  * Vafek and Kang [2023] O. Vafek and J. Kang, Physical Review B 107, 075123 (2023).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Physical Review Letters 123, 036401 (2019).

##  Appendix A Invertible transformation

In this appendix, we assume that a constant of motion
$\vec{X}\cdot\vec{\mu}=X_{1}\mu_{1}+X_{2}\mu_{2}+X_{3}\mu_{3}$ solution of

$\left[\vec{X}\cdot\vec{\mu},\mathcal{D}\right]=0\quad\Longleftrightarrow\quad\partial_{\bar{z}}\vec{X}+\alpha\vec{X}\times\vec{\mathcal{A}}=0,$
(10)

exist in order to

  * •

show that $\vec{X}\cdot\vec{X}$ is uniform in space and real.

  * •

derive the form of the effective gauge potential $A_{z}^{\prime}$ appearing in
Eq. 6.

  * •

prove that $\vec{X}$ cannot be carried by $\mu_{3}$ only away from high-
symmetry points of the moiré unit cell.

###  A.1 Uniform and real $x^{2}$

We first demonstrate that $x^{2}=\vec{X}\cdot\vec{X}$ is uniform and real. We
can first easily show that it is holomorphic

$\partial_{\bar{z}}(\vec{X}\cdot\vec{X})=2\vec{X}\cdot\partial_{\bar{z}}\vec{X}=-2\alpha\vec{X}\cdot(\vec{X}\times\mathcal{A})=0.$
(11)

Because $x$ is also additionally periodic (due to the periodicity of
$\mathcal{A}$), it must be constant (uniform in space). Then, we use the
$C_{2y}$ symmetry of the chiral Hamiltonian for TBG. This symmetry exchanges
the layers and complex conjugates the gauge fields
$A_{z}^{(j)}\to\left[A_{z}^{(j)}\right]^{*}$ for $j=1,2$, as can be derived
from their explicit representation (Eq. 3). As a result, we find that
$\vec{X}(-r_{1},r_{2})=\vec{X}^{*}(r_{1},r_{2})$, which ensures that all
components of $\vec{X}$ are real on the $y$-axis where $r_{1}=0$ [we use
$(r_{1},r_{2})$ for the position instead of $(x,y)$ to avoid confusion of the
horizontal position with the $x^{2}=\vec{X}\cdot\vec{X}$], which results in a
real $x^{2}$.

###  A.2 Abelianization

We then perform the calculation of the transformed Dirac operators
$\mathcal{D}^{\prime}=U^{-1}\mathcal{D}U$ with

$U=U^{-1}=\frac{1}{\sqrt{N}}(\vec{X}\cdot\vec{\mu}+x\mu_{3}),\quad
N=2x(x+X_{3}).$ (12)

A direct calculation gives

$\begin{split}\mathcal{D}^{\prime}&=U^{-1}\mathcal{D}U,\\\
&=\frac{1}{N}[(\vec{X}\cdot\vec{\mu})\mathcal{D}(\vec{X}\cdot\vec{\mu})+x^{2}\mu_{3}\mathcal{D}\mu_{3}]+\frac{x}{N}[\mu_{3}\mathcal{D}(\vec{X}\cdot\vec{\mu})+(\vec{X}\cdot\vec{\mu})\mathcal{D}\mu_{3}]+\sqrt{N}(-2i\partial_{\bar{z}}N^{-1/2})\mu_{0}\\\
&=\frac{1}{N}[\underbrace{(\vec{X}\cdot\vec{\mu})^{2}}_{x^{2}\mu_{0}}\mathcal{D}+x^{2}\mu_{3}\mathcal{D}\mu_{3}]+\frac{x}{N}[\mu_{3}(\vec{X}\cdot\vec{\mu})\mathcal{D}+(\vec{X}\cdot\vec{\mu})\mathcal{D}\mu_{3}]+\frac{i}{N}(\partial_{\bar{z}}N)\mu_{0}\\\
&=\frac{x^{2}}{N}[2\cdot(-2i\partial_{\bar{z}})\mu_{0}]+\frac{x}{N}[\underbrace{\\{\mu_{3},\vec{X}\cdot\vec{\mu}\\}}_{2X_{3}}(-2i\partial_{\bar{z}})+\underbrace{[\mu_{3},\vec{X}\cdot\vec{\mu}]}_{2i(\hat{u}_{3}\times\vec{X})\cdot\vec{\mu}}(\alpha\vec{\mathcal{A}}\cdot\vec{\mu})]+\frac{2ix}{N}(\partial_{\bar{z}}X_{3})\mu_{0}\\\
&=\frac{2x(x+X_{3})}{N}(-2i\partial_{\bar{z}})\mu_{0}-\frac{2ix}{N}[(\partial_{\bar{z}}X_{3})\mu_{0}+i\alpha(\vec{X}\cdot\vec{\mathcal{A}})\mu_{3}]+\frac{2ix}{N}(\partial_{\bar{z}}X_{3})\mu_{0}\\\
&=(-2i\partial_{\bar{z}})\mu_{0}+\frac{\alpha(\vec{X}\cdot\vec{\mathcal{A}})}{x+X_{3}}\mu_{3},\end{split}$
(13)

where the first line is a simple re-ordering, the second line uses the fact
that $(\vec{X}\cdot\mu)$ commutes with $\mathcal{D}$, the third uses that
$\vec{\mathcal{A}}$ has no $\mu_{3}$ component such that
$\mu_{3}(\vec{\mathcal{A}}\cdot\vec{\mu})\mu_{3}=-(\vec{\mathcal{A}}\cdot\vec{\mu})$,
the fourth and fifth simplify the products of Pauli matrices. This represents
two decoupled copies subject, each subject to an abelian $U(1)$ gauge field

$\pm A^{\prime}_{z}=\pm\frac{\alpha(\vec{X}\cdot\vec{\mathcal{A}})}{x+X_{3}},$
(14)

that is opposite for the two copies.

###  A.3 Out-of-plane constant of motion

We finally prove by contradiction that $\vec{X}\cdot\vec{\mu}$ cannot be
purely carried by $\mu_{3}$ at a point $r_{0}$ away from the origin and
corners of the Wigner-Seitz moiré unit cell. We will rely on the fact that $x$
is constant and reason by contradiction. Assume that we can write
$\vec{X}(r_{0})\cdot\vec{\mu}=x\mu_{3}$. Then, consider a small displacement
$\eta=(\eta_{x},\eta_{y})$ away from this point. The equation defining the
constant of motion Eq. 10 gives $X_{1}(r_{0}+\eta)=-x\alpha
A_{z}^{(1)}(r_{0})\bar{\eta}$, $X_{2}(r_{0}+\eta)=x\alpha
A_{z}^{(2)}(r_{0})\bar{\eta}$, and $X_{3}(r_{0}+\eta)=x$ with
$\bar{\eta}=\eta_{x}-i\eta_{y}$. For the norm $x^{2}$ to be constant, we
therefore need that
$x^{2}=x^{2}[1+\alpha^{2}\bar{\eta}^{2}[(A_{z}^{(1)})^{2}(r_{0})+(A_{z}^{(2)})^{2}(r_{0}))]$,
which requires $(A_{z}^{(1)})^{2}(r_{0})+(A_{z}^{(2)})^{2}(r_{0})=0$. This is
only true for the field in Eq. 3 at the origin and corners of the Wigner-Seitz
moiré unit cell, providing the desired contradiction.

##  Appendix B Numerical determination of magic angles

In this Appendix, we provide details on our numerical simulation of disorder
in the chiral TBG Hamiltonian. We use the method introduced in Ref. [56],
which we briefly summarize here for completeness.

Figure 8: a) Definition of the angle $\phi_{m,n}$ as the polar angle of
$\kappa_{m,n}$, measured from the negative $y$-axis, represented for the the
$m=3$ harmonics (red arrows). We also show the lowest harmonics corresponding
to the tunneling in the Bistritzer-MacDonald model in green. b) Definition of
the $\kappa_{m,n}$ vectors, which start from the origin and point to the red
dots. c) The flow of magic angles on the complex plane as a function of
increasing disorder strength $|\mu|$ (color bar) while $\hat{\mu}$ is fixed.
We put disorder in the first $6$ channels and use a truncation
$R_{k}=20\sqrt{3}\kappa$.

We begin by describing the chiral model without disorder. The Hamiltonian can
be written as

$\mathcal{H}_{\rm D}=\begin{pmatrix}0&\mathcal{D}\\\
\mathcal{D}^{\dagger}&0\end{pmatrix},\quad{\rm
where}\quad\mathcal{D}=\begin{pmatrix}\hat{D}&\alpha\hat{T}_{-}\\\
\alpha\hat{T}_{+}&\hat{D}\end{pmatrix}$ (15)

and $\hat{T}_{\pm}=\hat{T}(\pm z)$ is the inter-layer hopping written in terms
of the complex variable $z=x+iy$. The tunneling potential $\hat{T}$ is related
to the gauge field $\mathcal{A}_{z}$ defined in Eq. 3 of the main text by

$\hat{T}(z)=A_{z}^{(1)}+iA_{z}^{(2)}=\sum_{n=0,1,2}\omega^{n}\exp(i\kappa_{n}\cdot
r)$ (16)

with $\omega=e^{i2\pi/3}$, and where $\kappa_{n}\equiv\kappa_{1,n}$ are the
three smallest vectors connecting the two corners of the moiré Brillouin zone,
represented with green arrows in Fig. 8a. We consider the effects of disorder
in the inter-layer tunneling potential $\hat{T}_{\pm}$ through the addition of
a perturbation that we decompose in harmonics labeled by $m$

$\hat{T}\rightarrow\hat{T}+\sum\limits_{m}\mu_{m}\hat{T}_{m}$ (17)

where $\\{\mu_{m}\\}$ is a random vector of small norm
$|\mu|=[\sum_{m}\mu_{m}^{2}]^{1/2}\ll 1$. The wave-vectors corresponding to
the smallest ten harmonics ($m\leq 10$) are depicted in Fig. 8b. These
channels were used in Fig. 6 of the main text as well. Notice that, for each
$m$, there are several wave-vectors of equal norm related by $C_{3}$ and/or
$C_{2y}$ symmetry that contribute to $\hat{T}_{m}$. We denote them as
$\kappa_{m,n}$, with $n$ an auxiliary index taking either three or more
distinct values depending on the number of nonequivalent wave-vectors in the
$m$-th harmonic (compare for instance $m=2$ and $m=3$ in Fig. 8b). In the main
text, we considered perturbations $\hat{T}_{m}$ that preserve the $C_{3}$ and
$C_{2y}$ symmetry of the chiral Hamiltonian for TBG (Eq. 3), which constrains
the functional form of these inter-layer tunneling terms to

$\hat{T}_{m}=a_{m,z}^{(1)}+ia_{m,z}^{(2)}=\sum_{n}e^{i(\kappa_{m,n}\cdot
r+\phi_{m,n})}.$ (18)

As illustrated in Fig. 8a and b, these terms can be interpreted as long-range
tunnelings in momentum space, with an amplitude equal to $\mu_{m}$, a
direction given by $\kappa_{m,n}$, and a phase $\phi_{m,n}$. Amongst all gauge
equivalent choices for the $\phi_{m,n}$ phases leading to $C_{3}$ and $C_{2y}$
symmetric $\hat{T}_{m}$, we choose the one where $\phi_{m,n}$ is the polar
angle of $\kappa_{m,n}$ measured from the negative $y$-axis, as illustrated in
Fig. 8a. The main advantage of this prescription is that it equally applies to
harmonics containing three and six wave-vectors (Fig. 8b).

The problem of finding the magic angle corresponds to finding vectors
$|\psi\rangle$ and $|\chi\rangle$ such that

$\begin{pmatrix}\hat{D}&\alpha\hat{T}_{-}\\\
\alpha\hat{T}_{+}&\hat{D}\end{pmatrix}\begin{pmatrix}|\psi\rangle\\\
|\chi\rangle\end{pmatrix}=0\quad\Rightarrow\quad\Lambda|\psi\rangle=\dfrac{1}{\alpha^{2}}|\psi\rangle,\quad\Lambda=\hat{D}^{-1}\hat{T}_{+}\hat{D}^{-1}\hat{T}_{-},$
(19)

which itself reduces to finding the eigenvalues $\\{\lambda_{i}\\}$ of the
(non-Hermitian) operator $\Lambda$. The magic angles $\\{\alpha_{i}\\}$ are
then calculated as $\alpha_{i}=1/\sqrt{\lambda_{i}}$. The operators $\hat{D}$
and $\hat{T}_{\pm}$ can be represented as sparse matrices on a momentum space
grid and are, in practice, truncated by removing all momenta with norm higher
than a certain cutoff $R_{k}$.

As an illustration of the method, we introduce disorder in the first $m_{c}=6$
channels using $\vec{\mu}=|\mu|\hat{\mu}$, varying the overall strength of the
disorder $|\mu|$ for a fixed random vector normalized to unity $\hat{\mu}$. In
Fig. 8c, we show the evolution of the magic angles$\\{\alpha_{i}\\}$ on the
complex plane as a function of increasing disorder strength $|\mu|$. We see
that only the first magic angle is robust against the perturbation, while all
the other magic angles are dramatically changed as disorder strength
increases.

##  Appendix C Small Momentum Perturbation

The inter-layer tunneling decomposed into higher harmonic $m\geq 1$ of Eq. 17
can only describe disorder with a correlation length smaller than moiré unit
cell. In this Appendix, we explore the opposite regime where the disorder
correlation length is larger than the moiré unit cell. Physically, twist angle
variations across the sample belong to this class of perturbations.

Introducing harmonics that are periodic on an enlarged moiré unit cell in real
space requires access to wave-vectors that are fractions of the original
$\kappa_{1,n}$ appearing in the Bistritzer-MacDonald model (Eq. 3). To achieve
this numerically, we rescale the momentum-space lattice of Fig. 8b by a factor
$r=1/M$, with $M=3p+1$ and $p\in\mathbb{Z}$. This choice ensures that the
reduced lattice still contains the $\kappa_{1,n}$. By convention, rescaling by
a negative value $r$ corresponds to a rescaling by $|r|>0$ followed by a
$180^{\circ}$ degree rotation, as illustrated in Fig. 9a for the case $r=-1/2$
(corresponding to $M=-2$, $p=-1$). The inter-layer harmonics
$\hat{T}_{m}^{\prime}$ on the rescaled lattice with wave-vectors
$\kappa_{m,n}^{\prime}$ can now be used to introduce disorder with correlation
length larger than the moiré lattice constant. The definitions of
$\hat{T}_{m}^{\prime}$ and $\kappa_{m,n}^{\prime}$ are the same as
$\hat{T}_{m}$ and $\kappa_{m,n}$ (see App. B), albeit on the rescaled
honeycomb lattice. We also highlight that
$\kappa_{1,n}=\kappa_{|M|,n}^{\prime}$, which is also illustrated in Fig. 9a.

Figure 9: a) Illustration of the rescaled honeycomb lattice for $M=-2$ (grey)
lines. The rescaled lattice is generated by the basis vectors
$\kappa_{1,n}^{\prime}$ with $|\kappa^{\prime}_{1,n}|=|\kappa|/2$, and
contains the wave-vectors $|\kappa_{1,n}|$ appearing in the Bistritzer-
Macdonald model. The original moiré Brillouin zone is drawn in black. b)
Variation $|\partial\alpha/\partial\mu|$ of the first (blue), second (orange),
and third (green) magic angle, numerically evaluated using $\delta\mu=0.02$,
as a function of $M=\left|\frac{\kappa}{q_{\rm
pert}}\right|=\left|\frac{\kappa}{\kappa_{1,0}^{\prime}}\right|$. To obtain
converged results for large $m$, we used a large truncation radius
$R_{k}=35\sqrt{3}|\kappa|$.

We reproduce the calculation of Fig. 6a for the fractional wave-vectors
$\kappa_{m,n}^{\prime}$ and consider an inter-layer tunneling of the form

$\hat{T}^{\prime}=\hat{T}^{\prime}_{|M|}+\delta\mu\hat{T}_{1}^{\prime},$ (20)

where $\hat{T}^{\prime}_{|M|}$ is nothing but the terms appearing in the un-
perturbed chiral Hamiltonian for TBG (Eq. 3), while $\hat{T}_{1}^{\prime}$ is
the lowest accessible fractional harmonic on the rescaled lattice. The latter
realizes a perturbation of momentum $q_{\rm pert}=|\kappa/M|$. In Fig. 9b, we
show the derivative of the magic angles $|\partial\alpha/\partial\mu|$,
numerically computed using a small disorder strength $\delta\mu=0.02$ (see
Sec. IV.2), for the first three magic angles as a function of
$M=\left|\frac{\kappa}{q_{\rm
pert}}\right|=\left|\frac{\kappa}{\kappa_{1,0}^{\prime}}\right|$. We observe
that, unlike the perturbation of larger momentum where only the first magic
angle enjoyed protection, all magic angles respond similarly to disorder.
Thus, there is no special protection of the first magic angle against
perturbations of correlation length longer than the moiré lattice vector.