# Nature of even and odd magic angles in helical twisted trilayer graphene

Daniele Guerci  Center for Computational Quantum Physics, Flatiron Institute,
New York, New York 10010, USA  Yuncheng Mao  Université Paris Cité, CNRS,
Laboratoire Matériaux et Phénomènes Quantiques, 75013 Paris, France
Christophe Mora  Université Paris Cité, CNRS, Laboratoire Matériaux et
Phénomènes Quantiques, 75013 Paris, France

###### Abstract

Helical twisted trilayer graphene exhibits zero-energy flat bands with large
degeneracy in the chiral limit. The flat bands emerge at a discrete set of
magic twist angles and feature properties intrinsically distinct from those
realized in twisted bilayer graphene. Their degeneracy and the associated band
Chern numbers depend on the parity of the magic angles. Two degenerate flat
bands with Chern numbers $C_{A}=2$ and $C_{B}=-1$ arise at odd magic angles,
whereas even magic angles display four flat bands, with Chern number
$C_{A/B}=\pm 1$, together with a Dirac cone crossing at zero energy. All bands
are sublattice polarized. We demonstrate the structure behind these flat bands
and obtain analytical expressions for the wavefunctions in all cases. Each
magic angle is identified with the vanishing of a zero-mode wavefunction at
high-symmetry position and momentum. The whole analytical structure results
from whether the vanishing is linear or quadratic for the, respectively, odd
and even magic angle. The $C_{3z}$ and $C_{2y}T$ symmetries are shown to play
a key role in establishing the flat bands. In contrast, the particle-hole
symmetry is not essential, except from gapping out the crossing Dirac cone at
even magic angles.

##  I Introduction

Twisted trilayer graphene (TTG) has recently received significant attention
due to its structural similarities to twisted bilayer graphene while
exhibiting distinct features. Different configurations of TTG have been
considered based on the relative twist directions of the top and bottom
layers: (1) mirror symmetrical TTG [1, 2, 3, 4, 5, 6, 7, 8], where the top and
bottom layers are rotated in the same direction by the same angle relative to
the middle layer, (2) twisted monolayer-bilayer graphene [9, 10, 11, 12, 13,
14], where Bernal stacked bilayer graphene is rotated with a small angle
relative to the third layer, and (3) helical TTG [15, 16, 17, 18, 19, 20, 21,
22], where the top and bottom layers are rotated in opposite directions. By
adding an additional graphene layer on top of the bilayer, TTG offers
additional “knobs” to manipulate the system’s physical properties. Besides the
twist angle, the layer shifting and displacement field have also been
identified as key factors for altering the physical properties of TTG. These
features open up promising avenues for studying and controlling the unique
electronic properties of TTG systems.

Recent theoretical studies [18, 23, 24, 19, 20, 21, 22, 25] and experimental
findings [15] have focused on the helical configuration of TTG [16, 17]. This
specific structure is of interest due to the presence of two non-commensurate
moiré interference patterns, resulting in a moiré quasiperiodic crystal [15].
In the case of a small twist angle, the separation between length scales leads
to the introduction of an effective moiré lattice, defined as the closest
commensurate ratio [18], along with a small deviation that gives rise to the
supermoiré lengthscale [18, 19, 25]. The slow supermoiré periodicity can be
seen as a relative shift between the two different moiré scales and
parametrized with a displacement vector [18, 19].

Theoretical studies on helical trilayer graphene [20, 21, 22] have revealed
the significant impact of lattice relaxation effects on the structure at the
supermoiré scale. These studies have shown that energetically favorable ABA
and BAB stacking regions expand at the expense of AAA regions. This leads to a
real-space pattern characterized by triangular domains, where a sizeable gap
separates the two central low-energy bands from the remote ones. The ABA
regions are characterized by a total Chern number $C_{\rm tot}=1$ for the
central bands, while the BAB regions have $C_{\rm tot}=-1$, resulting in a
real-space Chern mosaic [19]. This mosaic is separated by a network of chiral
gapless regions. In the chiral limit [26], the electronic bandstructure
exhibits perfectly flat bands at a discrete series of magic angles. These flat
bands possess interesting properties distinct from those observed in twisted
bilayer graphene [26, 27, 28, 29, 30, 31, 32].

The study of flat bands, or zero modes, in the chiral limit consists in
finding the kernel of an operator with a purely holomorphic derivative and an
abelian or non-abelian periodic potential [33, 34, 35, 29, 31, 36, 37, 38].
Investigating the mathematical structure of this operator for helical trilayer
graphene, Popov and Tarnopolsky [25] have recently identified two possible
scenarios of flat bands for symmetric stackings (AAA, ABA, or BAB). The first
scenario, originally discussed in Ref. [19] (see also Ref. [20]), consists of
two flat bands: a color-entangled Chern 2 band and a Chern -1 band. The color-
entangled band [39, 40, 41, 12, 42, 43, 44, 45] is particularly interesting as
it cannot be simply reduced to a single Landau level. The second scenario,
originally proposed in Ref. [24], features a fourfold degenerate flatband
manifold with an additional Dirac cone crossing the flat bands at the $\Gamma$
point.

Figure 1:  (a) Long wavelength periodicity at the supermoiré lattice. Each
blue and orange point corresponds to a position where a precise AA stacking
occurs between, respectively, the top-middle and the middle-bottom layers. The
blue and orange sets of points show the two incommensurate moiré patterns. (b)
ABA configuration characterized by a relative shift of
$\bm{r}_{0}=(\bm{a}_{1}-\bm{a}_{2})/3$ between the two moiré patterns. (c) The
resulting moiré pattern involves the atomic scale configurations ABA, BAA and
AAB where the atom of one layer lie at the center of two overlying hexagons.
(d) Momentum space Brillouin zone with Dirac zero modes at $K$, $K^{\prime}$
and $\Gamma$, $\bm{q}_{j}$ shows the vectors defining the periodicity at the
moiré scale.

In this work, we demonstrate that helical trilayer graphene with ABA (or BAB)
stackings and equal twist angles displays a systematic series of magic angles,
where the two scenarios are alternatively realized. Odd magic angles feature
twofold degenerate flat bands following the first scenario, while even magic
angles have four degenerate flat bands and a Dirac cone as per the second
scenario. In each case, we prove the emergence of the flat band structure and
derive analytical expression for the zero-energy wavefunctions and the
resulting Chern numbers.

Our theory reveals interesting relations between the dimension of the vector
space spanned by the zero-energy modes and the total Chern number of the band.
It provides an orthogonality relation between chiral and anti-chiral zero
modes and demonstrates how the $C_{3z}$ and $C_{2y}T$ symmetries constrain the
wavefunctions asymptotics and ultimately protect the emergence of the flat
bands, in contrast with particle-hole symmetry which proves inessential.
Finally, we briefly investigate the fate of the flat bands when breaking both
the particle-hole and chiral symmetries.

The plan of the paper is as follows. In Sec. II we introduce the model and the
zero-mode equation in the chiral limit. Then, we discuss the symmetries of the
model protecting three zero-energy Dirac cones at $K$, $K^{\prime}$ and
$\Gamma$ of the mini Brillouin zone in Sec. II.2. Finally, we give in Sec.
II.3 the geometrical relations connecting the chiral and anti-chiral flatband
sectors. In Sec. III we move to the investigation of the zero-modes discussing
the different properties realized at even and odd magic angles. We provide the
analytical solution of the wavefunctions of the different flatbands and we
explain the different properties of odd Sec. III.1 and even Sec. III.2 magic
angles employing symmetry arguments as well as analytical results. In Sec. IV
we discuss the effect of the breaking of the particle-hole symmetry $P$.
Finally, we give a summary of the main results of our work in Sec. V.

##  II The chiral model for helical trilayer graphene

Helical trilayer graphene is formed by stacking three graphene layers in a
staircase configuration with equal twist angles, as depicted in Fig. 1. It
shows two incommensurate moiré interference patterns, resulting in a
superimposed modulation of the relative spatial shift between the two
patterns. This long distance periodic modulation defines a triangular
supermoiré, or moiré-of-moiré lattice. It can be seen as a slow variation of
the atomic registry, first at the moiré scale and then at the supermoiré scale
as shown in Fig. 1. The details of the model [19] we use for the single-
particle Hamiltonian are given in Appendix A. It is a continuum model [46, 47,
48, 49] parametrized by the local relative shift between the two moirés. The
parametrization evolves continuously between AAA, ABA and BAB local stacking
configurations over the supermoiré unit cell.

The atomic relaxation in helical trilayer graphene has been investigated in
recent works [20, 21, 22] where rearrangements of the atomic registry was
demonstrated. Theoretical calculations [20] (see also [21, 22]) performed at
small twist $\approx 1.5^{\circ}$ have revealed that relaxation favours the
formation of a triangular lattice with linear size $\approx 400$nm of large
domains characterized by the energetically favourable ABA/BAB stacking
expanding at the expense of AAA regions, separated by a network of domain
walls hosting chiral edge modes [21].

In this work, we shall focus on the spatially hegemonic ABA stacking shown in
Fig. 1b composed by the periodic modulation of the atomic configurations in
Fig. 1c. The BAB stacking is simply the $C_{2z}T$ symmetric of ABA.
Furthermore, we aim for simplicity and a deeper analytical understanding by
considering the chiral limit [26, 50, 29, 51, 52, 53] of our continuum model,
see Appendix A. Away from the chiral limit, the flat bands become dispersive
but their bandwidth is still much smaller than the energy gap to the remote
bands [19, 20], and they keep their topological features.

###  II.1 Zero modes

As discussed in Appendix A, the Hamiltonian for ABA stacking in the chiral
limit takes the form

$\mathcal{H}_{\rm ABA}(\bm{r})=\begin{pmatrix}0&\mathcal{D}(\bm{r})\\\
\mathcal{D}^{\dagger}(\bm{r})&0\end{pmatrix},$ (1)

for a single graphene valley ($K$), in the “sublattice-Chern” basis
$(\psi_{1},\psi_{2},\psi_{3},\chi_{1},\chi_{2},\chi_{3})$ [29, 51, 52, 53],
where $\psi_{\ell}$ and $\chi_{\ell}$ refer respectively to the $A$ and $B$
sublattice, and $\ell=1,2,3$ labels the three distinct layers from top to
bottom. We have introduced the differential operator

$\mathcal{D}(\bm{r})=-i\sqrt{2}I_{3\times 3}\partial+\mathcal{A}(\bm{r}),$ (2)

with the derivative $\partial=(\partial_{x}-i\partial_{y})/\sqrt{2}$
associated to the complex coordinate $z=(x+iy)/\sqrt{2}$, and the non-abelian
traceless gauge potential

$\mathcal{A}(\bm{r})=\begin{pmatrix}0&\alpha U_{\omega}(\bm{r})&0\\\ \alpha
U_{0}(-\bm{r})&0&\alpha U_{0}(\bm{r})\\\ 0&\alpha
U_{\omega}(-\bm{r})&0\end{pmatrix},$ (3)

resulting from electron tunneling between the layers with the dimensionless
strength $\alpha$. Given the wavevector modulation $\bm{q}_{j+1}=ie^{2i\pi
j/3}$ shown in Fig. 1c, for which we use a complex notation [16], the moiré
potentials are given by

$U_{0}(\bm{r})=\sum_{j=1}^{3}e^{-i\bm{q}_{j}\cdot\bm{r}},$ (4)

$U_{\omega}(\bm{r})=U_{0}(\bm{r}+\bm{r}_{0})$ and
$U_{\omega^{*}}(\bm{r})=U_{0}(\bm{r}-\bm{r}_{0})$ with
$\bm{r}_{0}=(\bm{a}_{1}-\bm{a}_{2})/3$. The moiré lattice vectors are
$\bm{a}_{j}=4\pi e^{i\pi/6}e^{i2\pi(j-1)/3}/3$ with $j=1,2$. In the
expressions above, all energy scales are expressed in unit of
$v_{F}k_{\theta}$, and all momentum (length) scales in units of $k_{\theta}$
($1/k_{\theta}$) with the moiré momentum $k_{\theta}=\theta K_{D}$,
$K_{D}=4\pi/3a_{\rm G}$ and graphene lattice constant $a_{\rm G}\approx
2.46$Å. $\theta$ is the twist angle between consecutive layers and
$\alpha=w_{\rm AB}/v_{F}k_{\theta}$.

An inspiring mathematical structure emerges [26, 19] in the chiral limit.
$\mathcal{H}_{\rm ABA}$ hence anticommutes $\\{\mathcal{H}_{\rm
ABA},\Lambda^{z}\\}=0$ with the chiral operator

$\Lambda^{z}=\begin{pmatrix}I_{3\times 3}&0\\\ 0&-I_{3\times 3}\end{pmatrix},$
(5)

and the search for zero-energy modes decomposes into

$\mathcal{D}(\bm{r})\bm{\chi}_{\bm{k}}(\bm{r})=0,\quad\mathcal{D}^{\dagger}(\bm{r})\bm{\psi}_{\bm{k}}(\bm{r})=0,$
(6)

where $\bm{\psi}$ and $\bm{\chi}$ are eigenstates of the chiral operator
$\Lambda^{z}$ with positive $\Lambda^{z}=1$ (chiral sector) and negative
$\Lambda^{z}=-1$ (anti-chiral sector) eigenvalue, respectively. The chiral and
anti-chiral sectors also refer to the $A$-sublattice and $B$-sublattice
polarized states, respectively. The zero modes must also satisfy the Bloch
periodic boundary conditions

$\begin{split}&\bm{\psi}_{\bm{k}}(\bm{r}+\bm{a}_{1/2})=e^{i\bm{k}\cdot\bm{a}_{1/2}}U_{\varphi}\bm{\psi}_{\bm{k}}(\bm{r}),\\\
&\bm{\chi}_{\bm{k}}(\bm{r}+\bm{a}_{1/2})=e^{i\bm{k}\cdot\bm{a}_{1/2}}U_{\varphi}\bm{\chi}_{\bm{k}}(\bm{r}),\end{split}$
(7)

inherited from Eq. 49, with $U_{\varphi}={\text{diag}}[\omega^{*},1,\omega]$.
The remainder of this paper will be devoted to analyzing the properties of the
zero modes, which are solutions of Eq.(6) under the boundary conditions
specified in Eq.(7).

###  II.2 Symmetries and Dirac cones

The Hamiltonian $\mathcal{H}_{\rm ABA}$ (1) is invariant under the spatial
symmetries $C_{3z}$ and $C_{2y}T$ forming the space group $P32^{\prime}1$
($\\#150.27$ in the BNS notation [54, 55]). In the sublattice basis we have:

$C_{3z}\mathcal{H}_{\rm ABA}(\bm{r})C^{-1}_{3z}=\mathcal{H}_{\rm
ABA}(C_{3z}\bm{r}),$ (8)

where

$C_{3z}=\begin{pmatrix}\omega^{*}&0&0\\\ 0&1&0\\\
0&0&\omega^{*}\end{pmatrix}\otimes\begin{pmatrix}\omega&0\\\
0&\omega^{*}\end{pmatrix}.$ (9)

$C_{2y}T$ is the composition of a spinless time-reversal symmetry and a two-
fold rotation around the $y$ axis:

$C_{2y}T\mathcal{H}_{\rm
ABA}(\bm{r})\left(C_{2y}T\right)^{-1}=\mathcal{H}_{\rm ABA}(C_{2y}\bm{r}),$
(10)

where the transformation exchanges top and bottom layer and includes the
complex conjugation operator $\mathcal{K}$:

$C_{2y}T=\begin{pmatrix}0&0&1\\\ 0&1&0\\\ 1&0&0\end{pmatrix}\otimes I_{2\times
2}\mathcal{K}.$ (11)

In addition to these spatial symmetries, the model also exhibits an emerging
particle-hole symmetry [16, 18, 19]

$P=\begin{pmatrix}0&0&1\\\ 0&-1&0\\\ 1&0&0\end{pmatrix}\otimes\sigma^{0},$
(12)

anticommuting with $\mathcal{H}_{\rm ABA}$

$P\mathcal{H}_{\rm ABA}(\bm{r})P^{-1}=-\mathcal{H}_{\rm ABA}(-\bm{r}).$ (13)

In general, i.e. for all twist angles $\theta$, the Hamiltonian ${\cal H}_{\rm
ABA}$ possesses three pairs (chiral and anti-chiral) of zero modes from which
three Dirac cones emerge. They are located at $\Gamma$, $K$ and $K^{\prime}$
of the mini Brillouin zone Fid. 1c corresponding to
$\bm{k}=0,\bm{q}_{1},-\bm{q}_{1}$, and originate from the Dirac cones of the
individual three graphene layers. The zero modes at $\Gamma$ are protected by
particle-hole symmetry $P$ whereas those at $K$ and $K^{\prime}$ are
stabilized by the anti-unitary particle-hole operator $PC_{2y}T$, as further
discussed in Appendix B. We remark that our symmetry analysis is also valid
away from the chiral limit [16, 18].

Magic angles are specific values of the twist angle $\theta$ (or $\alpha$ as
they are related to each other) at which the two central bands of ${\cal
H}_{\rm ABA}$ become perfectly flat. At these angles, the zero modes at
$\Gamma$, $K$, and $K^{\prime}$ are no longer unique, and an extensive
degenerate set of zero modes emerges, forming flat bands.

###  II.3 Geometrical relations

Irrespective of the twist angle and value of $\alpha$, we can derive a set of
identities which shows an interesting structure for the zero modes. Chiral and
anti-chiral zero-energy modes satisfy the relation:

$v(\bm{r})=\bar{\bm{\chi}}_{\bm{k}_{1}}(\bm{r})\cdot\bm{\psi}_{\bm{k}_{2}}(\bm{r})=0,$
(14)

for $\bm{k}_{1}\neq\bm{k}_{2}$ and arbitrary $\bm{r}$ with
$\bar{\bm{\chi}}_{\bm{k}_{2}}(\bm{r})\equiv{\bm{\chi}}^{*}_{\bm{k}_{2}}(\bm{r})$.
The proof of Eq. (14) is straightforward. A direct computation shows that
$\bar{\partial}v(\bm{r})=0\implies v(\bm{r})=v(z)$ by simply using that
${\bm{\chi}}_{\bm{k}_{1}}(\bm{r})$ and $\bm{\psi}_{\bm{k}_{2}}(\bm{r})$ are
zero modes of $\mathcal{D}(\bm{r})$ and $\mathcal{D}^{\dagger}(\bm{r})$,
respectively. $v(z)=0$ then follows from Liouville’s theorem and periodicity
over the moiré unit cell. Eq. (14) tells us that the chiral $\bm{\psi}$ and
anti-chiral $\bar{\bm{\chi}}$ solutions of different momenta at a given
$\bm{r}$ are orthogonal to each other. In addition, we also find that a chiral
zero mode can be generated from a pair of anti-chiral solutions. Namely,

$\bm{\psi}_{-\bm{k}_{1}-\bm{k}_{2}}(\bm{r})=\bar{\bm{\chi}}_{\bm{k}_{1}}(\bm{r})\times\bar{\bm{\chi}}_{\bm{k}_{2}}(\bm{r})$
(15)

solves
$\mathcal{D}^{\dagger}(\bm{r})\bm{\psi}_{-\bm{k}_{1}-\bm{k}_{2}}(\bm{r})=0$.
Similarly,

$\bm{\chi}_{-\bm{k}_{1}-\bm{k}_{2}}(\bm{r})=\bar{\bm{\psi}}_{\bm{k}_{1}}(\bm{r})\times\bar{\bm{\psi}}_{\bm{k}_{2}}(\bm{r}),$
(16)

solves $\mathcal{D}(\bm{r}){\bm{\chi}}_{-\bm{k}_{1}-\bm{k}_{2}}(\bm{r})=0$.
Expressions similar to Eq. (14) and Eq. (15) have been derived in Ref. [25]
but with a different choice of gauge.

Figure 2: Top panel: Wronskian $|W_{A}|$ (orange) and $|W_{B}|$ (blue) for the
chiral and anti-chiral sector as a function of $\alpha$. $|W_{A/B}|$ vanishes
at the magic angles denoted as red and green vertical lines. Bottom panel:
Figs. (a)-(c) show a schematic representation of the zero mode spinors
$\\{K,K^{\prime},\Gamma\\}$ in the chiral and anti-chiral sectors. Panel (a)
represents a generic configuration away from the magic angle where the spinors
are linearly independent and chiral and anti-chiral sectors satisfy Eq. (14),
while (b) and (c) sketch the two scenario realized at the magic angle in the
helical trilayer graphene. In panel (c) $\bm{\phi}_{\Gamma A}$ and
$\bar{\bm{\phi}}_{\Gamma B}$ are the two additional zero modes degenerate with
the flat bands at $\Gamma$.

We introduce the Wronskian [19, 25] of the Dirac spinors in the chiral sector:

$W_{A}(\bm{r})=\bm{\psi}_{\Gamma}(\bm{r})\cdot\left[\bm{\psi}_{K}(\bm{r})\times\bm{\psi}_{K^{\prime}}(\bm{r})\right],$
(17)

which obeys the relation $\bar{\partial}W_{A}=0$. Thus, according to
Liouville’s theorem, we have $W_{A}(\mathbf{r})=W_{A}$. $W_{A}\neq 0$ since,
unlike $v(\mathbf{r})$, Eq. (17) is invariant under translations of moiré
lattice vectors. $W_{A}\neq 0$ implies that the three vectors
$\bm{\psi}_{\Gamma}$, $\bm{\psi}_{K}$, and $\bm{\psi}_{K^{\prime}}$ are
linearly independent. Similarly, we can define the Wronskian in the anti-
chiral sector

$W_{B}(\bm{r})=\bm{\chi}_{\Gamma}(\bm{r})\cdot\left[\bm{\chi}_{K}(\bm{r})\times\bm{\chi}_{K^{\prime}}(\bm{r})\right]$
(18)

satisfying the condition $\partial W_{B}=0$, which implies
$W_{B}(\bm{r})=W_{B}$. Combining the definitions of the Wronskians with
Eqs.(15) and(16), we find that $W_{A}$ and $W_{B}$ are both proportional, up
to a phase, to the scalar products

$\bm{\bar{}}{\chi}_{\Gamma}(\bm{r})\cdot\bm{\psi}_{\Gamma}(\bm{r}),\quad\bm{\bar{}}{\chi}_{K}(\bm{r})\cdot\bm{\psi}_{K}(\bm{r}),\quad\bm{\bar{}}{\chi}_{K^{\prime}}(\bm{r})\cdot\bm{\psi}_{K^{\prime}}(\bm{r}),$
(19)

between the chiral and anti-chiral zero modes at $\Gamma$, $K$, and
$K^{\prime}$. Consequently, $W_{A}$ and $W_{B}$ vanish simultaneously with
these scalar products. As illustrated in the top panel of Fig. 2, the
simultaneous vanishing of $|W_{A}|$ and $|W_{B}|$ defines the positions of the
magic angles [19, 25].

From these different relations, an intuitive picture emerges. At non-magic
angles, the set of chiral vectors ${\cal
C}_{A}\equiv\\{\bm{\psi}_{\Gamma},\bm{\psi}_{K},\bm{\psi}_{K^{\prime}}\\}$ is
linearly independent and span the full three-dimensional space of layers at
each $\bm{r}$. The same holds for the anti-chiral set ${\cal
C}_{B}\equiv\\{\bm{\chi}_{\Gamma},\bm{\chi}_{K},\bm{\chi}_{K^{\prime}}\\}$
whereas chiral and anti-chiral vectors are mutually orthogonal when their
momenta are different as a result of Eq. (14). This is depicted in Fig. 2a.
The situation is quite different at magic angles because the scalar products
of Eq. (19) all vanish. The chiral ${\cal C}_{A}$ and anti-chiral ${\cal
C}_{B}$ sets then form two distinct subspaces which are orthogonal to each
other with the three-dimensional layer space. There are three possible cases,
illustrated in Fig. 2b and c - corresponding to the two scenarios of Ref. [25]
\- (i) the chiral set forms a subspace of dimension 2 and the anti-chiral set
is forced to be one-dimensional, i.e. all three anti-chiral vectors are
collinear. (ii) same as (i) but the role of chiral and anti-chiral are
exchanged. (iii) both chiral and anti-chiral sets have dimension 1.

As we will demonstrate in Sec. III, the chiral and anti-chiral flat bands that
emerge at the magic angle are generated by the two sets ${\cal C}_{A}$ and
${\cal C}_{B}$, leaving the subspace decomposition invariant: $2+1$ for (i),
$1+2$ for (ii), and $1+1$ for (iii). Additionally, for case (iii), we will
show that the third dimension is filled by a pair of additional zero modes at
the $\Gamma$ point, corresponding to a crossing of the flat bands by a Dirac
cone.

By decreasing the twist angle or increasing the value of $\alpha$, we will
uncover in Sec. III a series of magic angles alternating between cases (i) and
(iii) with an even/odd effect. Interestingly, we will also observe that the
dimension or rank of the flat band coincides with the absolute value of the
Chern number, highlighting the relationship between the rank of the flat band
and the number of lowest Landau levels comprising it [42, 43, 44, 45].

Figure 3:  a) Renormalized Fermi velocity at $K$ as a function of $\alpha$.
Vertical red and green lines denote odd $\alpha^{*}_{2n-1}$ and even
$\alpha^{*}_{2n}$ magic angles, respectively. Panel b) and c) show the
bandstructure at the first magic angle $\alpha^{*}_{1}\approx 0.377$
($\theta^{*}_{1}\approx 1.687^{\circ}$) and the second magic angle
$\alpha^{*}_{2}\approx 1.197$ ($\theta^{*}_{2}\approx 0.532^{\circ}$),
respectively. Red lines denote the flat bands, while $\times 2$ and $\times 4$
give the number of zero modes per $\bm{k}$-point for odd and even magic
angles, respectively.

##  III Even and odd series of magic angles

The sequence of consecutive magic angles is showcased by computing the
renormalized Fermi velocity at $K$ as a function of $\alpha$, see Fig. 3a, in
agreement with the Wronskians shown in Fig. 2. The odd and even magic angles
are highlighted by red and green vertical lines and correspond, respectively,
to the cases (i) and (iii) discussed in Sec. II.3. For odd magic angles
$\alpha^{*}_{2n-1}$, the first being $\alpha^{*}_{1}\approx 0.377$
($\theta^{*}_{1}\approx 1.687^{\circ}$), the spectrum features two degenerate
flat bands as displayed in Fig. 3b. For even magic angles $\alpha^{*}_{2n}$
however, the first one is $\alpha^{*}_{2}\approx 1.197$
($\theta^{*}_{2}\approx 0.532^{\circ}$), four degenerate flat bands arise
coexisting with a Dirac cone crossing them at $\Gamma$, as shown in Fig. 3c.
The structure and degeneracy of the zero-energy bands thus repeats
periodically and depends on the parity of the magic angle label.

Furthermore, inspecting more closely the sequence of magic angles, we find
that the difference between consecutive values rapidly approaches a constant
value
$\alpha^{*}_{2n+1}-\alpha^{*}_{2n-1}\approx\alpha^{*}_{2n+2}-\alpha^{*}_{2n}\simeq
1.214$, as shown in Fig. 4, similarly to the twisted bilayer case [26, 56, 57,
58].

Figure 4: a) Distance between neighboring magic angles in the even and odd
sequence. b) Sequence of magic angles $\alpha_{n}$. Increasing the order $n$
the distance between nearest neighbour magic angles in the even and odd
sectors approaches a constant value $\approx 1.214$ represented by the
horizontal gray line in a).

In the following, we demonstrate that the distinct nature of the even and odd
magic angles arises from symmetry considerations. These constraints dictate
the behavior of the zero-mode wavefunction at the high-symmetry points
together with the identities presented in Sec. II.3.

###  III.1 Zero modes for odd magic angles

To gain analytical insight into the origin of the odd magic angle, we study
the anti-chiral zero mode $\bm{\chi}_{\Gamma}(\bm{r})$ in the vicinity of the
high-symmetry point AA ($\bm{r}=0$). The wavefunctions around AA regions are
not described by pseudo Landau levels [59], reflecting a charge distribution
of the flatbands different from the one realized in twisted bilayer graphene
[20]. As detailed in Appendix C.1, the wavefunction
$\bm{\chi}_{\Gamma}(\bm{r})$ can be formally expanded in powers of $z$ and
$\bar{z}\equiv z^{*}$ close to $r=0$. Enforcing the symmetries $C_{3z}$,
$C_{2y}T$ and $P$ constrains and simplifies the resulting expansion. To
leading order, we find

$\bm{\chi}_{\Gamma}(\bm{r})=\begin{pmatrix}0\\\ \chi_{2}\\\
0\end{pmatrix}+O(\bar{z}),$ (20)

where $C_{2y}T$ (see Eq. (11)) imposes that $\chi_{2}\equiv\chi_{\Gamma 2}(0)$
is a real coefficient. Plotting the real $\chi_{2}$ as a function of $\alpha$
in Fig. 5a, we find that it vanishes for all odd magic angles. The condition
$\chi_{2}=0$, implying $\bm{\chi}_{\Gamma}(0)=0$, thus defines the series of
odd magic angles and yields a vanishing Wronskian from Eq. (18). The vertical
red lines of Fig. 5 exactly matches the one obtained from the Wronskian in
Fig. 2.

Beyond Eq. (20), the next order obeying symmetries is
$\bm{\chi}_{\Gamma}(\bm{r})\simeq\bar{z}(\chi^{\prime}_{1},0,{\chi^{\prime}}^{*}_{1})$
at the odd magic angles. Here, $\chi^{\prime}_{1}$ is an arbitrary complex
coefficient dependent on the magic angle. The presence of this simple zero is
sufficient to predict and explicitly construct [19] the whole anti-chiral flat
band following the seminal reasoning of Tarnopolsky et al. [26] in twisted
bilayer graphene, see also Refs. [27, 29, 30]. The zero-mode wavefunctions
exhibit the analytical expression:

$\bm{\chi}_{\bm{k}}(\bm{r})=\bar{\eta}_{\bm{k}}(\bar{z})\bm{\chi}_{\Gamma}(\bm{r}),$
(21)

where the antiholomorphic $\bar{\eta}_{\bm{k}}(\bar{z})=\eta^{*}_{\bm{k}}(-z)$
is related to the meromorphic function describing the lowest Landau level
(LLL) on a torus [60, 61]:

$\eta_{\bm{k}}(z)=e^{ik_{1}z/a_{1}}\frac{\vartheta_{1}[z/a_{1}-k/b_{2},\omega]}{\vartheta_{1}[z/a_{1},\omega]}$
(22)

with the notation $k_{1}=\bm{k}\cdot\bm{a}_{1}$ and the Jacobi theta-function

$\vartheta_{1}[z,\omega]=\sum_{n\in\mathbb{Z}}e^{i\pi\omega(n+1/2)^{2}}e^{2i\pi(z-1/2)(n+1/2)}$
(23)

which vanishes at $z=0$ and results in a Bloch periodicity (7) for Eq. (21).
Momentum space boundary conditions [30] on the self-periodic part of the Bloch
state
$\bm{u}_{\bar{k}}(\bm{r})=\bm{\chi}_{\bm{k}}(\bm{r})e^{-i\bm{k}\cdot\bm{r}}$
give

${\bm{u}}_{\bar{k}+\bar{b}_{j}}(\bm{r})=e^{-i{\bm{b}}_{j}\cdot\bm{r}}e^{i\phi_{{k,b_{j}}}}{\bm{u}}_{\bar{k}}(\bm{r}),$
(24)

with
$\phi_{{k,b_{1}}}=-2\pi\bar{k}/\bar{b}_{2}+\pi-\pi\bar{b}_{1}/\bar{b}_{2}$ and
$\phi_{{k,b_{2}}}=\pi$ corresponding to a flat band with total Chern number
$C_{B}=-1$ [30].

Thanks to the spatial symmetries, the vanishing of $\chi_{2}$ at the odd magic
angles is thus sufficient to predict a flat band (21) in the anti-chiral
sector with Chern number $-1$. Moreover, employing Eq. (16) also determines
the chiral sector. Eq. (16), evaluated at the origin $\bm{r}=0$ with
$\bm{k}_{1}=K$ and $\bm{k}_{2}=K^{\prime}$, gives

$\bm{\psi}_{K}(0)\times\bm{\psi}_{K^{\prime}}(0)=\bar{\bm{\chi}}_{\Gamma}(0)=0,$
(25)

resulting in the fact that ${\bm{\psi}}_{K}(0)$ and
${\bm{\psi}}_{K^{\prime}}(0)$ are collinear vectors, see Fig. 5a. Expanding
${\bm{\psi}}_{K}$ and ${\bm{\psi}}_{K^{\prime}}$ around $\bm{r}=0$ and
enforcing the symmetries as in Appendix C.1, we find to leading order
$\bm{\psi}_{K}(0)=(\psi_{1},0,\psi_{3})$ and
$\bm{\psi}_{K}(0)=(\psi_{3}^{*},0,\psi_{1}^{*})$ where $\psi_{1/3}$ are
complex coefficients. Using Eq. (25) yields

$\chi_{2}=|\psi_{3}|^{2}-|\psi_{1}|^{2}=0\implies\psi_{1}=\psi_{3}e^{i\varphi}.$
(26)

We choose a gauge with $\varphi=0$ and find

$\bm{\psi}_{K}(0)=\bm{\psi}_{K^{\prime}}(0).$ (27)

As shown in Ref. [19], this identity is sufficient to construct the flat band
wavefunctions for the chiral sector,

${\bm{\psi}}_{\bm{k}}(\bm{r})=a_{k}\eta_{\bm{k}+K^{\prime}}(z){\bm{\psi}}_{K}(\bm{r})+a_{-k}\eta_{\bm{k}+K}(z){\bm{\psi}}_{K^{\prime}}(\bm{r}),$
(28)

satisfying the Bloch periodicity Eq. (7), with the holomorphic function
defined in Eq. (22) and $a_{k}=\vartheta_{1}[(k+K)/b_{2},\omega]$. Momentum
space boundary conditions give $C_{A}=2$ as shown in Ref. [19]. The magic
relation (27) which is not associated to the vanishing of the chiral spinor
$\bm{\psi}$ justifies the charge density distribution homogeneity in the Chern
2 band [20].

Remarkably, all odd magic angles feature a chiral flat band of Chern $2$,
generated by the two vectors ${\bm{\psi}}_{K}(\bm{r})$ and
${\bm{\psi}}_{K^{\prime}}(\bm{r})$, alongside an anti-chiral flat band of
Chern $-1$ where all states align collinearly with
$\bm{\chi}_{\Gamma}(\bm{r})$. Moreover, Eq. (14) reveals that the chiral and
anti-chiral flat band spaces are orthogonal to each other. Specifically,
$\bm{\chi}_{\Gamma}(\bm{r})$ and all the anti-chiral wavefunctions are
oriented in a direction perpendicular to the chiral plane formed by
${\bm{\psi}}_{K}(\bm{r})$ and ${\bm{\psi}}_{K^{\prime}}(\bm{r})$. This
provides a clear understanding of why a chiral Chern $\pm 2$ band is
accompanied by an anti-chiral Chern $\mp 1$ band within a three-layer system.

Finally, we observe that the cross product
$\bm{\psi}_{K}(0)\times\bm{\psi}_{K^{\prime}}(0)$ vanishes also for even magic
angles, see black line in Fig. 5a but for a different reason that we shall
explain in the next Section.

Figure 5: a) Wavefunction absolute value $|\bm{\chi}_{\Gamma}(0)|$ (blue) and
cross product $|\bm{\psi}_{K}(0)\times\bm{\psi}_{K^{\prime}}(0)|$ (black) as a
function of $\alpha$. b) Wavefunction absolute values
$|\bm{\chi}_{K}(\bm{r}_{0})|$ and $|\bm{\psi}_{K^{\prime}}(\bm{r}_{0})|$ as a
function of $\alpha$. The grey solid line shows the renormalized Fermi
velocity $v^{*}/v_{F}$. Vertical green and red lines show the location of odd
and even magic angle, respectively.

###  III.2 Zero-modes for even magic angles

We now turn to the characterization of the zero-energy modes and flat bands
for even magic angles. In contrast to the previous case,
$\bm{\chi}_{\Gamma}(0)$ remains finite at $\alpha^{*}_{2n}$ as shown in
Fig.5a. Consequently, the zero-mode construction discussed in SectionIII.1
does not apply for even magic angles. To make progress, we examine the
behavior of the zero-mode solution $\bm{\psi}_{K^{\prime}}(\bm{r})$ in the
vicinity of the AB stacking point $\bm{r}_{0}=(\bm{a}_{1}-\bm{a}_{2})/3$. By
expanding the zero-mode equation (6) for the middle layer component
$\psi_{K^{\prime}2}$ to linear order in the deviation $\bm{r}$, we obtain:

$\bar{\partial}^{2}\psi_{K^{\prime}2}-\frac{9\alpha^{2}}{2}\left(\frac{z^{2}}{2}+\frac{\bar{z}}{\sqrt{2}}\right)\psi_{K^{\prime}2}=0,$
(29)

while top and bottom layer amplitudes are given by

$\begin{split}&\bar{\partial}\psi_{K^{\prime}1}=3\alpha
z\psi_{K^{\prime}2}/\sqrt{2},\\\
&\psi_{K^{\prime}3}=\left(i\sqrt{2}\bar{\partial}\psi_{K^{\prime}2}-3\alpha
z\psi_{K^{\prime}1}/\sqrt{2}\right)/(3\alpha).\end{split}$ (30)

We obtain the general solution of Eq. (29):

$\psi_{K^{\prime}2}(\bm{r}+\bm{r}_{0})=\gamma_{A}Ai(\zeta)+\gamma_{B}Bi(\zeta),$
(31)

where $\zeta=(2\alpha)^{2/3}\left(z^{2}/2+\bar{z}/\sqrt{2}\right)$ and we
retain both Airy functions $Ai(z)$ and $Bi(z)$ [62]. A similar behavior of the
zero-mode wavefunction has been obtained for twisted bilayer graphene around
the AB/BA points in Ref. [56]. $C_{3z}$ rotation centered at the high-symmetry
point $\bm{r}_{0}$ implies:

$\bm{\psi}_{K^{\prime}}(C_{3z}\bm{r}+\bm{r}_{0})=\begin{pmatrix}\omega^{*}&0&0\\\
0&\omega^{*}&0\\\
0&0&1\end{pmatrix}\bm{\psi}_{K^{\prime}}(\bm{r}+\bm{r}_{0}),$ (32)

employing the McLaurin expansion of the Airy functions [62] we find
$\gamma_{A}=-3^{2/3}N$, $\gamma_{B}=3^{1/6}N$:

$\bm{\psi}_{K^{\prime}}(\bm{r}+\bm{r}_{0})=\begin{pmatrix}0\\\ 0\\\
\psi_{3}^{(0)}\end{pmatrix}-\frac{3i\alpha\psi_{3}^{(0)}}{\sqrt{2}}\begin{pmatrix}0\\\
\bar{z}\\\ 0\end{pmatrix}+O(z^{2}),$ (33)

where $\psi_{3}^{(0)}\equiv\psi_{K^{\prime}3}(\bm{r}_{0})$. In addition, the
$PC_{2y}T$ symmetry gives:

$\bm{\psi}_{K^{\prime}}(\bm{r}+\bm{r}_{0})=\begin{pmatrix}-1&0&0\\\ 0&1&0\\\
0&0&-1\end{pmatrix}\bm{\psi}^{*}_{K^{\prime}}(C_{2x}\bm{r}+\bm{r}_{0}),$ (34)

implying $\Re\psi^{(0)}_{3}=0$. Consequently,
$\bm{\psi}_{K^{\prime}}(\bm{r}_{0})$ is determined by a single real number
$\Im\psi^{(0)}_{3}$. The condition for the even magic angle is, therefore, the
vanishing of this coefficient, such that
$\bm{\psi}_{K^{\prime}}(\bm{r}_{0})=0$, and the Wronskian $W_{A}$ also becomes
zero as per Eq.(17). The vertical green lines displayed in Fig. 5b align with
those in Fig. 2. At even magic angles, the wavefunction from Eq. (33) further
expands as $\bm{\psi}_{K^{\prime}}(\bm{r}_{0}+\bm{r})\simeq z^{2}(a,b,0)$ at
small $\bm{r}$ (with complex coefficients $a$ and $b$), corresponding to a
double zero. This structure enables the following analytical form for the
chiral flat bands

$\bm{\psi}_{\bm{k}}(\bm{r})=\eta^{(0)}_{\bm{k}^{\prime}}(z)\eta^{(0)}_{\bm{k}-\bm{k}^{\prime}+K}(z)\bm{\psi}_{K^{\prime}}(\bm{r}),$
(35)

with the meromorphic and periodic function

$\eta^{(0)}_{\bm{k}}(z)=e^{ik_{1}z/a_{1}}\frac{\vartheta_{1}[(z-z_{0})/a_{1}-k/b_{2},\omega]}{\vartheta_{1}[(z-z_{0})/a_{1},\omega]},$
(36)

$z_{0}$ denotes the complex representation of $\bm{r}_{0}$. Despite the fact
that $\bm{k}^{\prime}$ is an arbitrary wavevector in Eq. (35), at most two
distinct values of $\bm{k}^{\prime}$ yield independent wavefunctions. Setting
for instance, $\bm{k}^{\prime}=K^{\prime}$ and $\bm{k}^{\prime}=0$, we obtain
two degenerate flat bands described by the wavefunctions
$\bm{\psi}^{(1)}_{\bm{k}}(\bm{r})=\eta^{(0)}_{K^{\prime}}(z)\eta^{(0)}_{K^{\prime}+\bm{k}}(z)\bm{\psi}_{K^{\prime}}(\bm{r})$
and
$\bm{\psi}^{(2)}_{\bm{k}}(\bm{r})=\eta^{(0)}_{\bm{k}+K}(z)\bm{\psi}_{K^{\prime}}(\bm{r})$.

A very similar analytical form was introduced by Popov and Tarnopolsky [24,
25] to describe the fourfold degenerate flat bands but at AAA stacking magic
angles. Mathematically, the construction of Eq. (35) is possible because the
poles brought by the two functions $\eta^{(0)}_{\bm{k}}$ at $z=z_{0}$ are
precisely cancelled by the double zero of $\bm{\psi}_{K^{\prime}}(\bm{r})$ at
$\bm{r}_{0}$ ($z_{0}$). The resulting wavefunctions are finite everywhere,
obey the Bloch periodic boundary conditions Eq. (7), and solve the zero-mode
Eq. (6). All chiral flat band wavefunctions are collinear to
$\bm{\psi}_{K^{\prime}}(\bm{r})$ and thus span a one-dimensional space. It
readily explains the vanishing of the cross product
$\bm{\psi}_{K}(0)\times\bm{\psi}_{K^{\prime}}(0)$ also at even magic angle, as
shown in Fig. 5a. The above derivation relies on the double zero of
$\bm{\psi}_{K^{\prime}}$ at $\bm{r}_{0}$. Alternatively, we can reconstruct
the twofold degenerate flat band from the quadratic vanishing of
$\bm{\psi}^{(1)}_{\Gamma}(\bm{r})$ at $\bm{r}=0$. As discussed in Appendix
C.3, the advantage of this latter approach is that it does not require the
particle-hole symmetry $P$, only $C_{3z}$ and $C_{2y}T$.

Moving to the anti-chiral sector, we can repeat the same analysis for the
wavefunction $\bm{\chi}_{K}(\bm{r})$ in the vicinity of $\bm{r}_{0}$.
$\bm{\chi}_{K}(\bm{r}_{0})=0$ at even magic angles, as shown in Fig. 5b, and
the expansion around $r_{0}$ is quadratic indicating a double zero. The same
construction thus extends to the anti-chiral sector, resulting in two
degenerate flat bands. Additionally, by applying momentum space boundary
conditions, we can determine that each chiral flat band possesses a Chern
number $C_{A}=+1$, whereas the anti-chiral flat bands have $C_{B}=-1$. In
total, we prove that even magic angles feature a two-fold degenerate set of
bands in each sublattice sector, or four flat bands with a vanishing total
Chern number.

The chiral flat bands align with $\bm{\psi}_{K^{\prime}}(\bm{r})$ while the
anti-chiral with $\bm{\chi}_{K}(\bm{r})$, leaving room for a third direction
in the layer space. By following the arguments presented in Ref.[25], one can
analytically demonstrate the existence of a pair of additional zero modes at
$\Gamma$. These zero modes correspond to a Dirac cone crossing the flat bands
in Fig.3c. The proof is outlined as follows. From Eq. (14), we know that
$\bm{\psi}_{K}(\bm{r})$ and $\bm{\chi}_{K^{\prime}}(\bm{r})$ are orthogonal to
each other. Both wavefunctions vanish at $\bm{r}_{0}$ and remain finite
everywhere else. Consequently, we can define a function
$\bm{\phi}_{\Gamma}(\bm{r})$ such that:

$\bm{\psi}_{K^{\prime}}(\bm{r})=\bar{\bm{\phi}}_{\Gamma}(\bm{r})\times\bar{\bm{\chi}}_{K}(\bm{r}),$
(37)

Here, we emphasize that this expression does not uniquely define
$\bm{\phi}_{\Gamma}(\bm{r})$. The subscript $\Gamma$ indicates that
$\bm{\phi}_{\Gamma}(\bm{r})$ has the Bloch periodicity of the $\Gamma$ point,
as inferred from the periodicity of the two other functions. Applying the
operator $\mathcal{D}^{\dagger}(\bm{r})$ to both sides of Eq. (37), we arrive
at

$0=-\left[\mathcal{D}^{*}(\bm{r})\bar{\bm{\phi}}_{\Gamma}(\bm{r})\right]\times\bar{\bm{\chi}}_{K}(\bm{r}),$
(38)

which shows that $\mathcal{D}^{*}(\bm{r})\bar{\bm{\phi}}_{\Gamma}(\bm{r})$
must be proportional to $\bar{\bm{\chi}}_{K}(\bm{r})$, or

$\mathcal{D}^{*}(\bm{r})\bar{\bm{\phi}}_{\Gamma}(\bm{r})=f(\bm{r})\bar{\bm{\chi}}_{K}(\bm{r})$
(39)

with some periodic function $f(\bm{r})$. We introduce the function
$g(\bm{r})$, solution of $\bar{\partial}g(\bm{r})=f(\bm{r})$, and shift
$\bar{\bm{\phi}}_{\Gamma}$ as

$\bar{\bm{\phi}}_{\Gamma}(\bm{r})\to\bar{\bm{\phi}}_{\Gamma}(\bm{r})-ig(\bm{r})\bar{\bm{\chi}}_{K}(\bm{r})/\sqrt{2}$
(40)

to finally obtain

$\mathcal{D}(\bm{r}){\bm{\phi}}_{\Gamma}(\bm{r})=0.$ (41)

This last equation demonstrates that we have constructed an additional anti-
chiral zero-energy solution. Due to its definition in Eq.(37),
$\bm{\phi}_{\Gamma}(\bm{r})$ cannot be proportional to $\bm{\chi}_{K}(\bm{r})$
and therefore lies outside the one-dimensional anti-chiral subspace. Eq.(37)
also implies that
$\bar{\bm{\phi}}_{\Gamma}(\bm{r})\cdot\bm{\psi}_{K^{\prime}}(\bm{r})=0$:
$\bm{\phi}_{\Gamma}(\bm{r})$ and $\bm{\psi}_{K^{\prime}}(\bm{r})$ are
orthogonal. As a result, $\bm{\phi}_{\Gamma}(\bm{r})$ occupies the third
vacant direction. As already discussed earlier, $\bm{\phi}_{\Gamma}(\bm{r})$
satisfies the Bloch periodic condition with momentum at $\Gamma$. A similar
construction yields a second zero-energy state at $\Gamma$ in the chiral
sector, also spanning the third direction orthogonal to both the chiral and
anti-chiral flat bands. This result concludes the characterization of the low-
energy spectrum at even magic angles.

##  IV Breaking the particle-hole symmetry $P$

The twist-angle dependency in the Pauli matrices obtained by replacing
$\bm{\sigma}$ with $\bm{\sigma}_{\pm\theta}$, see Eq. 45, breaks the particle-
hole symmetry $P$ (12). The symmetry breaking is negligible in the small twist
angle regime $\theta$ and has been ignored in the previous analysis. In the
following we first discuss the stability of the zero modes to this
perturbation. Then, we move to consider the effect of the perturbation away
from the chiral limit.

###  IV.1 Robustness of the flat bands

In the chiral limit the twist angle difference between top and bottom layer
$\bm{\sigma}_{\pm\theta}$ enters in the Hamiltonian $H_{\rm ABA}$ (1) as:

$\mathcal{D}(\bm{r})=-i\sqrt{2}M_{-\theta}\partial+\mathcal{A}(\bm{r}),$ (42)

where $M_{\theta}={\text{diag}}\left(e^{i\theta},1,e^{-i\theta}\right)$.
Differently from twisted bilayer graphene [26] the layer dependent phases
$M_{\theta}$ in Eq. (42) cannot be gauged away since
$[M_{\theta/2},\mathcal{A}(\bm{r})]\neq 0$. As mentioned before, the layer
dependent phase breaks the particle-hole symmetry $P$ [18, 20] while $C_{3z}$
and $C_{2y}T$ are still symmetries of $\mathcal{H}_{\rm ABA}$. In the chiral
limit $w_{\rm AA}=0$ the Hamiltonian $\mathcal{H}_{\rm ABA}$ is characterized
by three Dirac cones at $K$, $K^{\prime}$ and $\Gamma$ that are protected by
$C_{3z}$ and the chiral symmetry $\Lambda_{z}$. From Eq. (42) we readily
realize that the orthogonality relation (14) transforms into the identity:

$v_{\theta}(\bm{r})=\bar{\bm{\chi}}_{\bm{k}_{1}}(\bm{r})\cdot\left[M_{\theta}\bm{\psi}_{\bm{k}_{2}}(\bm{r})\right]=0$
(43)

for $\bm{k}_{1}\neq\bm{k}_{2}$ and arbitrary $\bm{r}$. We emphasize that this
is no longer a vector product as it is not definite positive. Relations (15)
and (16) are also modified accordingly and become
$\bm{\psi}_{-\bm{k}_{1}-\bm{k}_{2}}=[M_{\theta}\bar{\bm{\chi}}_{\bm{k}_{1}}]\times[M_{\theta}\bar{\bm{\chi}}_{\bm{k}_{2}}]$
and
$\bm{\chi}_{-\bm{k}_{1}-\bm{k}_{2}}=[M_{-\theta}\bar{\bm{\psi}}_{\bm{k}_{1}}]\times[M_{-\theta}\bar{\bm{\psi}}_{\bm{k}_{2}}]$,
while Eqs. (17) and (18) are still satisfied.

At odd magic angles, the linear vanishing of $\bm{\chi}_{\Gamma}(\bm{r})$
around $\bm{r}=0$ still holds, protected by the symmetries $C_{3z}$ and
$C_{2y}T$. The positions of the magic angles are only slightly shifted by the
particle-hole symmetry breaking term, $\alpha_{1}^{*}\approx 0.3772$ for
instance. Moreover, since

$\bar{\bm{\chi}}_{\Gamma}(0)=0=[M_{\theta}{\bm{\psi}}_{K}]\times[M_{\theta}{\bm{\psi}}_{K^{\prime}}],$
(44)

we still obtain that $\bm{\psi}_{K}(0)=\bm{\psi}_{K^{\prime}}(0)$ at the odd
magic angle, yielding the Chern $2$ zero-mode chiral solution of Eq. (28)
together with the anti-chiral Chern $1$ band. The structure of the flat bands
are thus fully stable under particle-hole symmetry breaking as confirmed by
our numerical calculations.

At even magic angles, as discussed in Appendix C.3, the $C_{3z}$ and $C_{2y}T$
symmetries alone are sufficient to preserve a double zero in one of the two
wavefunctions $\bm{\psi}_{\Gamma}(\bm{r})$ at $\bm{r}=0$, which automatically
yields the two flat bands in the anti-chiral sector. The analytical
expressions for these bands are explicitly given in Eq. (64). It is worth
noting that the choice of $\bm{\psi}_{\Gamma}(\bm{r})$ vanishing at $\bm{r}=0$
is not continuously connected to the zero-mode solution $\bm{\psi}_{\Gamma}$
away from magic angles; instead, it requires the admixture with the first
excited band, which occurs at even magic angles. In the anti-chiral sectors,
the two flat bands are protected by the conservation of ${\rm
Tr}\Lambda^{z}=0$ (see Eq. (5)) within the zero-energy manifold. This
conservation condition dictates that the two chiral flat bands must be
accompanied by two anti-chiral flat bands. Consequently, we find that the
fourfold degenerate flat band structure is maintained at even magic angles,
even in the presence of particle-hole symmetry breaking. Our numerical
calculations confirm this protection. Finally, as already noted in Ref. [20],
the additional Dirac cone at $\Gamma$ is gapped by breaking $P$. The gap is
however quite small, on the order of the energy scale $\sim\theta
v_{F}k_{\theta}$.

###  IV.2 Away from the chiral limit

Figure 6: Separated middle bands of the first (left panel) and the second
(right panel) magic angles, labeled with the Chern number associated to each
band. The separated band are obtained by breaking simultaneously the particle-
hole and the chiral symmetry. To amplify the effect of the particle-hole
symmetry breaking the twist angle in $\bm{\sigma}_{\pm\theta}$ is multiplied
by a factor 20.

The primary emphasis of this work has been on the chiral limit, which allowed
us to provide analytical wavefunctions for the flat bands. We now introduce a
slight deviation from this limit by incorporating a finite corrugation factor
$w_{AA}/w_{AB}=0.05$, while still accounting for the particle-hole symmetry
breaking term. In this case, all zero-energy modes are lifted and gaps open
between bands as illustrated in Fig. 6 in the vicinity of the first two magic
angles. At twist angles different from magic values, the Dirac cones at
$\Gamma$, $K$ and $K^{\prime}$ become gapped when both the chiral and
particle-hole symmetries are broken. At magic angles, whether even or odd, the
flat bands acquires a finite bandwidth, on the order of a few meV for
$w_{AA}/w_{AB}=0.05$, and separate from each other.

As discussed in Sec. III.1, the odd magic angles in the chiral limit exhibit a
flat band with Chern number $+2$ ($-1$) which is fully polarized on the A(B)
sublattice. Remarkably, even when departing from the chiral limit, where the
two sublattices become mixed, the resulting two bands still retain their Chern
numbers of $+2$ and $-1$, as shown in Fig.6 for the first magic angle. These
Chern numbers are in fact robust and persist up to significant values of the
corrugation as the gap to the remote bands is large.

In contrast to this, as discussed in Sec. III.2, we identified a set of four
flat bands with Chern numbers $+1,+1,-1,-1$, and complete sublattice
polarization in the chiral limit for even magic angles. The Dirac cone
crossing at $\Gamma$, as shown in the previous section, is gapped out by
particle-hole symmetry breaking. When we move outside the chiral limit, the
mixing of sublattices leads to the rearrangement of Chern numbers to
$+1,0,0,-1$, as depicted in Fig.6, for the four split bands at a small
corrugation of $w_{AA}/w_{AB}=0.05$. Increasing the corrugation further leads
to topological transitions that involve the Dirac cone at $\Gamma$ and change
the topological properties of the flat bands. Their study is beyond the scope
of this work.

##  V Conclusions

This work analyses the mathematical structure of flat bands in equal-twist
helical trilayer graphene in the chiral limit and for ABA stacking. We
determine analytical expressions for all zero-energy wavefunctions at magic
angles together with the band Chern numbers. We derive an orthogonality
relation between chiral and anti-chiral zero modes, which constrains the
number of generators in the zero-energy manifold and reveals a connection
between the dimensionality of the vector space spanned by the zero modes and
the total Chern number of the band.

In contrast to twisted bilayer graphene, helical twisted trilayer graphene
exhibits an even/odd variation in the composition and features of the flat
bands. At odd magic angles, we find a twofold degenerate zero-energy manifold
at each $\bm{k}$ point, comprising an anti-chiral flat band with Chern number
$C_{B}=-1$ and a chiral flat band with Chern number $C_{A}=2$. The latter,
being generated by two linearly independent spinors, cannot be reduced to the
lowest Landau level, leading to interesting implications on the properties of
the correlated ground state [19], which will be the topic of future studies.

Even magic angles, on the other hand, are distinguished by a four-fold
degenerate manifold, where both chiral and anti-chiral flat bands have pairs
of, or double, zeros, resulting in two zero-modes in each sector. The two
sublattice-polarized bands in a given sector are all collinear to a single
space-dependent spinor and carry the same Chern number, $C_{A}=+1$ for the
chiral $A$-polarized bands, $C_{B}=-1$ for the anti-chiral $B$-polarized
bands. We also demonstrate that in addition to the two flat bands, there must
be two additional degenerate zero modes within the zero-energy manifold,
thereby explaining the presence of a Dirac cone crossing the flat bands at
$\Gamma$.

We also investigate the stability of flat bands and zero modes at magic angles
under a weak breaking of particle-hole symmetry. We find that all the features
listed above remain valid, with the exception of the Dirac cone, which becomes
slightly gapped by the perturbation. Only the joint breaking of particle-hole
and chiral symmetries gaps out and splits all bands. Interestingly, at odd
magic angle, the resulting isolated bands retain the Chern numbers $+2,-1$
analytically determined in the chiral limit. In the light of of Ref. [63], it
is natural to ask which features protect the emergence and properties of the
flat bands. Our analysis highlights the importance of the structure of the
differential operator $\mathcal{D}=-i\sqrt{2}\partial+\mathcal{A}$, where
$\mathcal{A}$ is a non-abelian traceless SU(3) gauge potential and there is
clearly a natural generalization to the SU(N) case for multilayer stackings
with $N>3$, where Chern bands with $C>2$ are expected. In addition, the
symmetries $C_{3z}$ and $C_{2y}T$ appear to be crucial for maintaining exactly
flat bands tuned only by the twist angle, whereas particle-hole symmetry is
not needed. We emphasize that, in contrast to twisted bilayer graphene, the
$C_{2z}T$ symmetry is broken here and the model belongs to the Altland-
Zirnbauer [64] class AIII instead of CI, although both classes can have
protected Dirac cone on a 2D surface.

Finally, the emergence of isolated bands with non-zero Chern values in a more
realistic context, where both chiral and particle-hole symmetries are broken,
opens up exciting possibilities for realizing an anomalous quantum Hall effect
without the need for an almost aligned hBN substrate. Moreover, the potential
for fractional Hall states consisting of bands with Chern numbers of $\pm 2$
represents a promising direction for future investigations.

###### Acknowledgements.

We acknowledge discussions with Jie Wang, Jen Cano and Eslam Khalaf. C.M. and
Y.M. acknowledge support by the French National Research Agency (project
TWISTGRAPH, ANR-21-CE47-0018). D.G. acknowledges support from the Flatiron
Institute, a division of the Simons Foundation.

##  Appendix A Local Hamiltonian

We review here some fundamental properties of the local Hamiltonian describing
equal twist helical trilayer graphene. A more general formulation, which
includes non-equal twist configurations, is provided in Ref [18, 25]. In the
basis ${\bm{\Psi}}=(\psi_{1},\chi_{1},\psi_{2},\chi_{2},\psi_{3},\chi_{3})$
where $\psi,\,\chi$ correspond to the wave function amplitude on the A and B
sublattices, respectively, the Hamiltonian near valley K reads [18, 19]

$H_{\rm
eTTG}(\bm{r};{\bm{\phi}})=\begin{pmatrix}v_{F}\hat{\bm{k}}\cdot\bm{\sigma}_{\theta}&T(\bm{r},{\bm{\phi}})&0\\\
h.c.&v_{F}\hat{\bm{k}}\cdot\bm{\sigma}&T(\bm{r},-{\bm{\phi}})\\\
0&h.c.&v_{F}\hat{\bm{k}}\cdot\bm{\sigma}_{-\theta}\end{pmatrix},$ (45)

where $v_{F}\approx 10^{6}$m/s is the graphene velocity. The set of phases
${\bm{\phi}}=(\phi_{1},\phi_{2},\phi_{3})$ parametrizes the position on the
supermoiré lattice. With the choice of gauge $\phi_{1}=0$,

$\bm{R}=\frac{\phi_{2}}{\pi}{\mathbf{a}}^{\rm
MM}_{1}+\frac{\phi_{3}}{\pi}{\mathbf{a}}^{\rm MM}_{2},$ (46)

with the supermoiré lattice vectors ${\mathbf{a}}^{\rm
MM}_{1/2}=\frac{4\pi}{3\theta k_{\theta}}e^{\mp i\pi/3}$. ${\bm{\phi}}$ also
controls the relative shift between the two moiré patterns and a change of
gauge shifts all phases $\phi_{j}$ by the same amount. $\bm{R}=0$ defines the
AAA stacking, whereas ABA and BAB correspond to $\bm{R}=({\mathbf{a}}^{\rm
MM}_{2}-{\mathbf{a}}^{\rm MM}_{1})/3$ and $\bm{R}=({\mathbf{a}}^{\rm
MM}_{1}-{\mathbf{a}}^{\rm MM}_{2})/3$ parametrized by
${\bm{\phi}}=\pm(0,-\pi/3,+\pi/3)$.

${\bm{\sigma}}$ is the vector of Pauli matrices in the sublattice space,
$\bm{\sigma}_{\theta}\equiv
e^{i\theta\sigma^{z}/2}\bm{\sigma}e^{-i\theta\sigma^{z}/2}$ and
$\hat{\bm{k}}=-i\nabla_{\bm{r}}$. The tunneling between different layers is
described by the moiré potential:

$T(\bm{r},{\bm{\phi}})=\sum^{3}_{j=1}T_{j}e^{-i\bm{r}\cdot\bm{q}_{j}}e^{-i\phi_{j}},$
(47)

where $T_{j+1}=w_{\rm AA}\sigma^{0}+w_{\rm AB}[\sigma^{x}\cos 2\pi
j/3+\sigma^{y}\sin 2\pi j/3]$, $w_{\rm AB}=110$meV and $w_{\rm AA}=rw_{\rm
AB}$ with $r$ dimensionless parameter quantifying atomic corrugation, using
complex notation $\bm{q}_{j+1}=ik_{\theta}e^{2i\pi j/3}$ [16] with
$k_{\theta}=\theta K_{D}$, $K_{D}=4\pi/3a_{\rm G}$ and $a_{\rm G}\approx
2.46$Å. The moiré lattice is characterized by the reciprocal lattice vectors
$\bm{b}_{1/2}=\bm{q}_{1}-\bm{q}_{2/3}$ and primitive vectors $\bm{a}_{1/2}$.
Ignoring the twist angle dependency in the Pauli matrices
$\bm{\sigma}_{\pm\theta}$ the Hamiltonian (45) becomes invariant under the
particle-hole symmetry:

$PH_{\rm eTTG}(\bm{r};{\bm{\phi}})P^{-1}=-H_{\rm eTTG}(-\bm{r};{\bm{\phi}}),$
(48)

with $P$ given in Eq. (12). Under moiré lattice translations we have:

$H_{\rm eTTG}(\bm{r}+\bm{a}_{1/2};\bm{\phi})=U_{\varphi}H_{\rm
eTTG}(\bm{r};\bm{\phi})U^{\dagger}_{\varphi},$ (49)

where the matrix
$U_{\varphi}={\text{diag}}[\omega^{*},1,\omega]\otimes\sigma^{0}$ corresponds
to a layer dependent phase factor $\omega=\exp(2\pi i/3)$.

$H_{\rm eTTG}$ also exhibits a supermoiré periodicity with $\bm{R}$ (or
${\bm{\phi}}$). One can show [19] the following identity

$H_{\rm
eTTG}\left(\bm{r}+\frac{\bm{a}_{l}}{2};{\bm{\phi}}+\Delta{\bm{\phi}}_{l}\right)=\tilde{U}H_{\rm
eTTG}(\bm{r};{\bm{\phi}})\tilde{U}^{\dagger},$ (50)

where $\tilde{U}={\text{diag}}(1,1,\omega)$. The phase shifts
$\Delta{\bm{\phi}}_{1}=(0,\pi,0)$, $\Delta{\bm{\phi}}_{2}=(0,0,\pi)$
correspond respectively to $\bm{R}\to\bm{R}+{\mathbf{a}}^{\rm MM}_{1}$ and
$\bm{R}\to\bm{R}+{\mathbf{a}}^{\rm MM}_{2}$. As a result, the AAA stacking
points are periodically replicated, forming a triangular lattice generated by
${\mathbf{a}}^{\rm MM}_{1/2}$ and characterized by ABA and BAB domains [19]
(see also [20, 21]).

As noted above, we focus on the ABA stacking configuration by setting
${\bm{\phi}}=(0,2\pi/3,-2\pi/3)$ and the chiral limit $w_{AA}=0$ (suppressed
tunneling between A and A orbitals). The resulting Hamiltonian is given by Eq.
(1) in the sublattice-Chern basis.

##  Appendix B Particle-hole symmetry $P$ protected Dirac cones

We discuss the symmetries protecting the Dirac cones at the high-symmetry
points $K$, $K^{\prime}$ and $\Gamma$ in the ABA stacking configuration.

The irreducible representations of the space group composed by $C_{3z}$ and
$C_{2y}T$ at the high-symmetry momenta $\Gamma$ and $K$ ($K^{\prime}$) verify
the $C_{3}$ point group character table and are all one-dimensional. This can
be understood more directly by considering an eigenstate of $C_{3z}$,
$C_{3z}|\omega\rangle=\omega|\omega\rangle$. From the relation

$(C_{2y}T)C_{3z}(C_{2y}T)^{-1}=C_{3z}^{-1},$ (51)

we obtain $C_{3z}C_{2y}T|\omega\rangle=\omega\,C_{2y}T|\omega\rangle$.
$C_{2y}T$ thus does not circulate between the eigenstates of $C_{3z}$ and
cannot protect a twofold degeneracy. The three zero-energy Dirac cones at $K$,
$K^{\prime}$ and $\Gamma$ arising in the band spectrum of ABA trilayer
graphene [18, 19] are therefore stable in the presence of the particle-hole
symmetry $P$ (12). In the chiral limit $w_{\rm AA}=0$, the Dirac cones are
protected by $\Lambda^{z}$ and persists even in the absence of $P$.

Since $P$ and $C_{3z}$ commute, if the spectrum at $\Gamma$ hosts two states
with $C_{3z}$ eigenvalues $\omega$, $\omega^{*}$ near charge neutrality (and
no other states), particle-hole symmetry $P$ automatically pins these two
states at zero energy. This is proven by contradiction: if we assume that the
two states sit at opposite non-vanishing energies, then $P$ permutes them.
This is however impossible since $P$ cannot change the $C_{3z}$ eigenvalue
which completes the proof. In fact, $P$ restricted to these two states must be
the identity as it commutes with $C_{3z}$. It further shows that breaking
$C_{3z}$ does not lift the Dirac crossings as the trivial (identity)
representation of $P=I_{2\times 2}$ cannot deform continuously to the
traceless $\sigma^{x}$ matrix permuting states with opposite non-zero
energies. $K$ and $K^{\prime}$ are however not stable under $P$ \- which
permutes $K$ and $K^{\prime}$ \- and the stability of their Dirac cones must
come from a different operator. $C_{2y}T$ sends $\bm{k}\to-C_{2y}\bm{k}$ and
therefore also permutes $K$ and $K^{\prime}$. Combining $P$ and $C_{2y}T$ we
find:

$PC_{2y}T\mathcal{H}_{\rm
ABA}(\bm{r})\left(PC_{2y}T\right)^{-1}=-\mathcal{H}_{\rm ABA}(-C_{2y}\bm{r}),$
(52)

with

$P^{\prime}=PC_{2y}T=\begin{pmatrix}1&0&0\\\ 0&-1&0\\\
0&0&1\end{pmatrix}\otimes I_{2\times 2}\mathcal{K},$ (53)

which leaves $K$ and $K^{\prime}$ invariant and acts as a (anti-unitary)
particle-hole operator. $P^{\prime}C_{3z}(P^{\prime})^{-1}=C^{-1}_{3z}$ being
anti-unitary $P^{\prime}$ cannot permute the eigenvalues $(\omega,\omega^{*})$
which implies that an isolated pair of states at $K$ and $K^{\prime}$ is
degenerate and pinned at zero energy, in the two-dimensional subspace
$P^{\prime}=I_{2\times 2}\mathcal{K}$.

##  Appendix C Symmetries of the zero-mode wavefunction

The symmetries of the Hamiltonian $\mathcal{H}_{\rm ABA}$ discussed in Section
(II.2) constraints the properties of the zero-mode wavefunction and gives
insight on the nature of the different magic angles. Considering the symmetry
$g$ the A sublattice polarized solution $\bm{\psi}$ we have:

$\bm{\psi}_{\bm{k}}(\bm{r})=e^{i\Xi_{g}}M_{g}\bm{\psi}_{g\bm{k}}(g\bm{r}),$
(54)

where the action of the symmetry $g$ on layer, momentum and space degrees of
freedom is given in table 1.

Symmetry | $\bm{k}$ | $\bm{r}$ | $M_{g}$
---|---|---|---
$C_{3z}$ | $C_{3z}\bm{k}$ | $C_{3z}\bm{r}$ | ${\text{diag}}(\omega,1,\omega)$
$P$ | $-\bm{k}$ | $-\bm{r}$ | $\begin{pmatrix}0&0&-1\\\ 0&1&0\\\ -1&0&0\end{pmatrix}$
$C_{2y}T$ | $-C_{2y}\bm{k}$ | $C_{2y}\bm{r}$ | $\begin{pmatrix}0&0&1\\\ 0&1&0\\\ 1&0&0\end{pmatrix}\mathcal{K}$
$PC_{2y}T$ | $C_{2y}\bm{k}$ | $-C_{2y}\bm{r}$ | ${\text{diag}}(-1,1,-1)\mathcal{K}$
Table 1: Action of the symmetries of $\mathcal{H}_{\rm ABA}$ on the zero-mode
$\bm{\psi}$. Rows refer to the different symmetries, from top to bottom
$C_{3z}$, $P$, $C_{2y}T$ and $PC_{2y}T$ the last two involving the complex
conjugation $\mathcal{K}$. Columns show the action of the symmetry on the
momentum $\bm{k}$, space $\bm{r}$ coordinates and $M_{g}$ is the
representation of the symmetry acting on the three dimensional layer degree of
freedom.

Similar expressions are also obtained for the B sublattice zero mode
$\bm{\chi}_{\bm{k}}$. The phase $\Xi_{g}$ is fixed by taking the decoupled
limit $\alpha=0$ and depends on $\bm{k}$.

We will now employ the symmetry relations in Eq. (54) of the model to
constrain the zero mode wavefunctions in the ${\cal C}_{A}$ and ${\cal C}_{B}$
sectors around high-symmetry points $\bm{r}$ in the moiré unit cell.

###  C.1 Odd magic angles

Thanks to Eq. (25) odd magic angles are determined by looking at the behavior
of $\bm{\chi}_{\Gamma}(\bm{r})$ around $\bm{r}\approx 0$ where the kernel
$\mathcal{D}(\bm{r})$ takes the form:

$\mathcal{D}(\bm{r})\simeq\begin{pmatrix}-i\sqrt{2}\partial&-3\alpha
z/\sqrt{2}&0\\\ 3\alpha&-i\sqrt{2}\partial&3\alpha\\\ 0&3\alpha
z/\sqrt{2}&-i\sqrt{2}\partial\end{pmatrix}.$ (55)

The solution of the zero-mode equation is obtained performing a Taylor
expansion in $z$ and $\bar{z}$. Employing Eq. (54) we readily realize that
$\chi_{\Gamma 2}(C_{3z}\bm{r})=\chi_{\Gamma 2}(\bm{r})$ and $\chi_{\Gamma
1/3}(C_{3z}\bm{r})=\omega^{*}\chi_{\Gamma 1/3}(\bm{r})$. In addition we also
have $\chi_{\Gamma 2}(C_{2y}\bm{r})=\chi^{*}_{\Gamma 2}(\bm{r})$ and
$\chi_{\Gamma 1/3}(C_{2y}\bm{r})=\chi^{*}_{\Gamma 3/1}(\bm{r})$ implying:

$\bm{\chi}_{\Gamma}(\bm{r})\simeq\chi_{2}\begin{pmatrix}3i\alpha z^{2}/4\\\
1\\\ -3i\alpha z^{2}/4\end{pmatrix}+\begin{pmatrix}\chi^{\prime}_{1}\bar{z}\\\
3\sqrt{2}\alpha\Im\chi^{\prime}_{1}z\bar{z}\\\
{\chi^{\prime}}^{*}_{1}\bar{z}\end{pmatrix},$ (56)

where $\chi_{2}\equiv\chi_{\Gamma 2}(0)\in\mathbb{R}$ while
$\chi^{\prime}_{1}\equiv\bar{\partial}\chi_{\Gamma 1}|_{0}\in\mathbb{C}$.
Thus, $C_{3z}$ and $C_{2y}T$ reduces the perfect flatness of the entire band
to the vanishing of a single real number $\chi_{2}$ [63]. Notice that if we
further impose the particle-hole symmetry $P$ we have
$\Re\chi^{\prime}_{1}=0$. We readily realize that the expression (56) solves
$\mathcal{D}\bm{\chi}_{\Gamma}=0$ up to small terms of the order
$z^{2}\bar{z}$. At the magic angle we have $\chi_{2}=0$ and
$\bm{\chi}_{\Gamma}$ has a simple zero for $\bm{r}\to 0$:

$\bm{\chi}_{\Gamma}(\bm{r}\to
0)\sim\bar{z}\left(\chi^{\prime}_{1},0,{\chi^{\prime}}^{*}_{1}\right)^{T}.$
(57)

We emphasize that the simple zero at the odd magic angles where $\chi_{2}=0$
persists also in the absence of the particle-hole symmetry $P$. The vanishing
of $\bm{\chi}_{\Gamma}(0)$ implies
$\bm{\psi}_{K}(0)\times\bm{\psi}_{K^{\prime}}(0)=0$ resulting in the fact that
$\bm{\psi}_{K}(0)$ and $\bm{\psi}_{K^{\prime}}(0)$ are colinear. Eq. (54) for
$C_{3z}$ gives
$\psi_{K/K^{\prime}2}(C_{3z}\bm{r})=\omega\psi_{K/K^{\prime}2}(\bm{r})$ and
$\psi_{K/K^{\prime}1/3}(C_{3z}\bm{r})=\psi_{K/K^{\prime}1/3}(\bm{r})$.
Furthermore, $K$ and $K^{\prime}$ zero modes are related by
$\psi_{K^{\prime}1/3}(\bm{r})=\psi^{*}_{K3/1}(C_{2y}\bm{r})$ and
$\psi_{K^{\prime}2}(\bm{r})=\psi^{*}_{K2}(C_{2y}\bm{r})$. These symmetries
reduces $\bm{\psi}_{K}(0)\times\bm{\psi}_{K^{\prime}}(0)=0$ to
$\bm{\psi}_{K}(0)=\bm{\psi}_{K^{\prime}}(0)$, where the identity holds up to a
phase, which gives rise to the the Chern 2 zero mode wavefunction in Eq. (28).

###  C.2 Even magic angles

Even magic angles realizes the a 1+1 decomposition corresponding to one-
dimensional flat bands in both chiral and anti-chiral sectors. The four-fold
degeneracy of the flat band manifold, see Fig. 3c, originates from a double
zero in the spinors $\bm{\psi}_{K^{\prime}}(\bm{r})$ and
$\bm{\chi}_{K}(\bm{r})$ at $\bm{r}_{0}$. The origin of this double zero can be
explained by symmetry reasoning. To start with we expand
$\mathcal{D}^{\dagger}(\bm{r})$ around the AB stacking point $\bm{r}_{0}$
finding:

$\mathcal{D}^{\dagger}(\bm{r}+\bm{r}_{0})\simeq\begin{pmatrix}-i\sqrt{2}\bar{\partial}&-3\alpha
z/\sqrt{2}&0\\\ 3\alpha z/\sqrt{2}&-i\sqrt{2}\bar{\partial}&3\alpha\\\
0&-3\alpha\bar{z}/\sqrt{2}&-i\sqrt{2}\bar{\partial}\end{pmatrix}.$ (58)

Focusing on the chiral sector and fixing the center of the $C_{3z}$ rotation
around $\bm{r}_{0}$ we find Eq. (32) which implies:

$\bm{\psi}_{K^{\prime}}(\bm{r}+\bm{r}_{0})\simeq\psi^{(0)}_{3}\begin{pmatrix}0\\\
-i3\alpha\bar{z}/\sqrt{2}\\\
1\end{pmatrix}+\frac{z^{2}}{2}\begin{pmatrix}\psi^{(0)^{\prime\prime}}_{1}\\\
\psi^{(0)^{\prime\prime}}_{2}\\\ 0\end{pmatrix},$ (59)

where
$\psi^{(0)^{\prime\prime}}_{1/2}\equiv\partial^{2}\psi_{K^{\prime}1/2}|_{\bm{r}_{0}}$
and we have included terms up to second order in $\bm{r}$. Imposing, the
$PC_{2y}T$ symmetry (34) implies $\Re\psi^{(0)}_{3}=0$,
$\Re\psi^{(0)^{\prime\prime}}_{1}=0$ and $\Im\psi^{(0)^{\prime\prime}}_{2}=0$.
Thus, at the magic angle where $\psi^{(0)}_{3}=0$ we have

$\bm{\psi}_{K^{\prime}}(\bm{r}\to 0+\bm{r}_{0})\sim
z^{2}\left(\psi^{(0)^{\prime\prime}}_{1},\psi^{(0)^{\prime\prime}}_{2},0\right)^{T}/2,$
(60)

enabling to attach two lowest Landau levels with simple pole at $\bm{r}_{0}$
to the zero mode spinor $\bm{\psi}_{K^{\prime}}(\bm{r})$ (35). Following a
similar line or reasoning one could show that the anti-chiral zero mode
$\bm{\chi}_{K}$ shows a double zero at $\bm{r}_{0}$.

###  C.3 Alternative derivation for even magic angles

We provide here an alternative argument for the protection of a twofold
degenerate flat band at even magic angles. The construction in Sec. III.2
relies on the quadratic vanishing of $\bm{\psi}_{K^{\prime}}(\bm{r})$ at
$\bm{r}=\bm{r}_{0}$ derived in Sec. C.2 which uses the $PC_{2y}T$ symmetry,
and therefore $P$. However, the zero-mode flat bands also host a specific
wavefunction
$\bm{\psi}^{(1)}_{\Gamma}(\bm{r})=\eta^{(0)}_{K^{\prime}}(z)\eta^{(0)}_{K^{\prime}}(z)\bm{\psi}_{K^{\prime}}(\bm{r})$
which exhibits a double zero at $\bm{r}=0$. As we show below, this double zero
is protected solely by the $C_{3z}$ and $C_{2y}T$ symmetries. To linear order
in $z$, the zero mode equation takes the form:

$\begin{split}&-i\sqrt{2}\bar{\partial}\psi_{\Gamma
2}-3\alpha\bar{z}(\psi_{\Gamma 1}-\psi_{\Gamma 3})/\sqrt{2},\\\
&-i\sqrt{2}\bar{\partial}(\psi_{\Gamma 1}-\psi_{\Gamma 3})=0,\\\
&-i\sqrt{2}\bar{\partial}(\psi_{\Gamma 1}+\psi_{\Gamma 3})+6\alpha\psi_{\Gamma
2}=0.\end{split}$ (61)

Solving the system of differential equations and imposing $C_{3z}$, Eq. (54),
we find that the zero mode behaves as:

$\bm{\psi}_{\Gamma}(\bm{r})\simeq\psi_{2}\begin{pmatrix}0\\\ 1\\\
0\end{pmatrix}-\frac{3\alpha i\bar{z}}{\sqrt{2}}\psi_{2}\begin{pmatrix}1\\\
0\\\ 1\end{pmatrix}+z^{2}\begin{pmatrix}\psi^{\prime\prime}_{1}/2\\\ 0\\\
\psi^{\prime\prime}_{3}/2\end{pmatrix}$ (62)

where $\psi_{2}\equiv\psi_{\Gamma 2}(0)$ and
$\psi^{\prime\prime}_{1/3}=\partial^{2}\psi_{\Gamma 1/3}|_{0}$. $C_{2y}T$
gives the additional condition $\psi_{\Gamma
1/3}(C_{2y}\bm{r})=\psi^{*}_{\Gamma 3/1}(\bm{r})$ and $\psi_{\Gamma
2}(C_{2y}\bm{r})=\psi^{*}_{\Gamma 2}(\bm{r})$ implying $\psi_{2}\in\mathbb{R}$
and $\psi^{\prime\prime}_{1}={\psi^{\prime\prime}}^{*}_{3}$. In summary, a
double zero is expected,

$\bm{\psi}_{\Gamma}(\bm{r}\to 0)\sim
z^{2}\left(\psi^{{}^{\prime\prime}}_{1},0,{\psi^{\prime\prime}}^{*}_{1}\right)^{T}/2$
(63)

as soon as the real coefficient $\psi_{2}$ vanishes. It corresponds to the
specific solution $\bm{\psi}^{(1)}_{\Gamma}(\bm{r})$ introduced above.
Breaking $P$ but keeping $C_{3z}$ and $C_{2y}T$ intact, one only moves the
magic angle for which $\psi_{2}=0$ since $\psi_{2}$ remains a real number. The
twofold degenerate flat band in the chiral sector is a direct consequence of
Eq. (63), with the analytical structure

$\bm{\psi}_{\bm{k}}(\bm{r})=\eta^{(0)}_{\bm{k}^{\prime}}(z)\eta^{(0)}_{\bm{k}-\bm{k}^{\prime}}(z)\bm{\psi}_{\Gamma}(\bm{r}),$
(64)

alternative to Eq. (35).

## References

  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021).
  * Kim  _et al._ [2022] H. Kim, Y. Choi, C. Lewandowski, A. Thomson, Y. Zhang, R. Polski, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Evidence for unconventional superconductivity in twisted trilayer graphene, Nature 606, 494 (2022).
  * Kim  _et al._ [2023] H. Kim, Y. Choi, Étienne Lantagne-Hurtubise, C. Lewandowski, A. Thomson, L. Kong, H. Zhou, E. Baum, Y. Zhang, L. Holleis, K. Watanabe, T. Taniguchi, A. F. Young, J. Alicea, and S. Nadj-Perge, Imaging inter-valley coherent order in magic-angle twisted trilayer graphene (2023), arXiv:2304.10586 [cond-mat.str-el] .
  * Călugăru  _et al._ [2021] D. Călugăru, F. Xie, Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Twisted symmetric trilayer graphene: Single-particle and many-body hamiltonians and hidden nonlocal symmetries of trilayer moiré systems with and without displacement field, Phys. Rev. B 103, 195411 (2021).
  * Xie  _et al._ [2021] F. Xie, N. Regnault, D. Călugăru, B. A. Bernevig, and B. Lian, Twisted symmetric trilayer graphene. ii. projected hartree-fock study, Physical Review B 104, 10.1103/physrevb.104.115167 (2021).
  * Lei  _et al._ [2021] C. Lei, L. Linhart, W. Qin, F. Libisch, and A. H. MacDonald, Mirror symmetry breaking and lateral stacking shifts in twisted trilayer graphene, Physical Review B 104, 10.1103/physrevb.104.035139 (2021).
  * Christos  _et al._ [2022] M. Christos, S. Sachdev, and M. S. Scheurer, Correlated insulators, semimetals, and superconductivity in twisted trilayer graphene, Phys. Rev. X 12, 021018 (2022).
  * Guerci  _et al._ [2022] D. Guerci, P. Simon, and C. Mora, Higher-order van hove singularity in magic-angle twisted trilayer graphene, Phys. Rev. Res. 4, L012013 (2022).
  * Chen  _et al._ [2020] S. Chen, M. He, Y.-H. Zhang, V. Hsieh, Z. Fei, K. Watanabe, T. Taniguchi, D. H. Cobden, X. Xu, C. R. Dean, and M. Yankowitz, Electrically tunable correlated and topological states in twisted monolayer–bilayer graphene, Nature Physics 17, 374 (2020).
  * He  _et al._ [2021] M. He, Y.-H. Zhang, Y. Li, Z. Fei, K. Watanabe, T. Taniguchi, X. Xu, and M. Yankowitz, Competing correlated states and abundant orbital magnetism in twisted monolayer-bilayer graphene, Nature Communications 12, 10.1038/s41467-021-25044-1 (2021).
  * Polshyn  _et al._ [2021] H. Polshyn, Y. Zhang, M. A. Kumar, T. Soejima, P. Ledwith, K. Watanabe, T. Taniguchi, A. Vishwanath, M. P. Zaletel, and A. F. Young, Topological charge density waves at half-integer filling of a moiré superlattice, Nature Physics 18, 42 (2021).
  * Ledwith  _et al._ [2022a] P. J. Ledwith, A. Vishwanath, and E. Khalaf, Family of ideal chern flatbands with arbitrary chern number in chiral twisted graphene multilayers, Phys. Rev. Lett. 128, 176404 (2022a).
  * Yang [2023] M. Yang, Flat bands and high chern numbers in twisted multilayer graphene (2023), arXiv:2303.00103 [math-ph] .
  * Suárez Morell  _et al._ [2013] E. Suárez Morell, M. Pacheco, L. Chico, and L. Brey, Electronic properties of twisted trilayer graphene, Phys. Rev. B 87, 125414 (2013).
  * Uri  _et al._ [2023] A. Uri, S. C. de la Barrera, M. T. Randeria, D. Rodan-Legrain, T. Devakul, P. J. D. Crowley, N. Paul, K. Watanabe, T. Taniguchi, R. Lifshitz, L. Fu, R. C. Ashoori, and P. Jarillo-Herrero, Superconductivity and strong interactions in a tunable moiré quasiperiodic crystal (2023), arXiv:2302.00686 [cond-mat.mes-hall] .
  * Mora  _et al._ [2019] C. Mora, N. Regnault, and B. A. Bernevig, Flatbands and perfect metal in trilayer moiré graphene, Phys. Rev. Lett. 123, 026402 (2019).
  * Zhu  _et al._ [2020] Z. Zhu, S. Carr, D. Massatt, M. Luskin, and E. Kaxiras, Twisted trilayer graphene: A precisely tunable platform for correlated electrons, Phys. Rev. Lett. 125, 116404 (2020).
  * Mao  _et al._ [2023] Y. Mao, D. Guerci, and C. Mora, Supermoiré low-energy effective theory of twisted trilayer graphene, Phys. Rev. B 107, 125423 (2023).
  * Guerci  _et al._ [2023] D. Guerci, Y. Mao, and C. Mora, Chern mosaic and ideal flat bands in equal-twist trilayer graphene (2023), arXiv:2305.03702 [cond-mat.mes-hall] .
  * Devakul  _et al._ [2023] T. Devakul, P. J. Ledwith, L.-Q. Xia, A. Uri, S. de la Barrera, P. Jarillo-Herrero, and L. Fu, Magic-angle helical trilayer graphene (2023), arXiv:2305.03031 [cond-mat.str-el] .
  * Nakatsuji  _et al._ [2023] N. Nakatsuji, T. Kawakami, and M. Koshino, Multi-scale lattice relaxation in general twisted trilayer graphenes (2023), arXiv:2305.13155 [cond-mat.mes-hall] .
  * Foo  _et al._ [2023] D. C. W. Foo, Z. Zhan, M. M. A. Ezzi, L. Peng, S. Adam, and F. Guinea, Extended magic phase in twisted graphene multilayers (2023), arXiv:2305.18080 [cond-mat.str-el] .
  * Dunbrack and Cano [2023] A. Dunbrack and J. Cano, Intrinsically multilayer moiré heterostructures, Phys. Rev. B 107, 235425 (2023).
  * Popov and Tarnopolsky [2023a] F. K. Popov and G. Tarnopolsky, Magic angles in equal-twist trilayer graphene (2023a), arXiv:2303.15505 [cond-mat.str-el] .
  * Popov and Tarnopolsky [2023b] F. K. Popov and G. Tarnopolsky, Magic angle butterfly in twisted trilayer graphene (2023b), arXiv:2305.16385 [cond-mat.str-el] .
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of magic angles in twisted bilayer graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Ledwith  _et al._ [2020] P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath, Fractional chern insulator states in twisted bilayer graphene: An analytical approach, Phys. Rev. Research 2, 023237 (2020).
  * Popov and Milekhin [2021] F. K. Popov and A. Milekhin, Hidden wave function of twisted bilayer graphene: The flat band as a landau level, Phys. Rev. B 103, 155150 (2021).
  * Ledwith  _et al._ [2021] P. J. Ledwith, E. Khalaf, and A. Vishwanath, Strong coupling theory of magic-angle graphene: A pedagogical introduction, Annals of Physics 435, 168646 (2021).
  * Wang  _et al._ [2021] J. Wang, Y. Zheng, A. J. Millis, and J. Cano, Chiral approximation to twisted bilayer graphene: Exact intravalley inversion symmetry, nodal structure, and implications for higher magic angles, Physical Review Research 3, 10.1103/physrevresearch.3.023155 (2021).
  * Sheffer and Stern [2021] Y. Sheffer and A. Stern, Chiral magic-angle twisted bilayer graphene in a magnetic field: Landau level correspondence, exact wave functions, and fractional chern insulators, Phys. Rev. B 104, L121405 (2021).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Magic-angle twisted bilayer graphene as a topological heavy fermion problem, Phys. Rev. Lett. 129, 047601 (2022).
  * Atiyah and Singer [1984] M. Atiyah and I. M. Singer, Dirac operators coupled to vector potentials, Proceedings of the National Academy of Sciences 81, 2597 (1984).
  * San-Jose  _et al._ [2012] P. San-Jose, J. González, and F. Guinea, Non-abelian gauge potentials in graphene bilayers, Phys. Rev. Lett. 108, 216802 (2012).
  * Katsnelson and Prokhorova [2008] M. I. Katsnelson and M. F. Prokhorova, Zero-energy states in corrugated bilayer graphene, Phys. Rev. B 77, 205424 (2008).
  * Ledwith  _et al._ [2022b] P. J. Ledwith, A. Vishwanath, and D. E. Parker, Vortexability: A unifying criterion for ideal fractional chern insulators (2022b), arXiv:2209.15023 [cond-mat.str-el] .
  * Gao  _et al._ [2022] Q. Gao, J. Dong, P. Ledwith, D. Parker, and E. Khalaf, Untwisting moiré physics: Almost ideal bands and fractional chern insulators in periodically strained monolayer graphene (2022), arXiv:2211.00658 [cond-mat.mes-hall] .
  * Parhizkar and Galitski [2023] A. Parhizkar and V. Galitski, A generic topological criterion for flat bands in two dimensions (2023), arXiv:2301.00824 [cond-mat.mes-hall] .
  * Barkeshli and Qi [2012] M. Barkeshli and X.-L. Qi, Topological nematic states and non-abelian lattice dislocations, Phys. Rev. X 2, 031013 (2012).
  * Wu  _et al._ [2013] Y.-L. Wu, N. Regnault, and B. A. Bernevig, Bloch model wave functions and pseudopotentials for all fractional chern insulators, Phys. Rev. Lett. 110, 106802 (2013).
  * Wu  _et al._ [2014] Y.-L. Wu, N. Regnault, and B. A. Bernevig, Haldane statistics for fractional chern insulators with an arbitrary chern number, Phys. Rev. B 89, 155113 (2014).
  * Wang  _et al._ [2022] J. Wang, S. Klevtsov, and Z. Liu, Origin of model fractional chern insulators in all topological ideal flatbands: Explicit color-entangled wavefunction and exact density algebra (2022), arXiv:2210.13487 [cond-mat.mes-hall] .
  * Dong  _et al._ [2022] J. Dong, P. J. Ledwith, E. Khalaf, J. Y. Lee, and A. Vishwanath, Exact many-body ground states from decomposition of ideal higher chern bands: Applications to chirally twisted graphene multilayers, arXiv:2210.13477 (2022).
  * Mera and Ozawa [2023] B. Mera and T. Ozawa, Uniqueness of landau levels and their analogs with higher chern numbers (2023), arXiv:2304.00866 [cond-mat.mes-hall] .
  * Estienne  _et al._ [2023] B. Estienne, N. Regnault, and V. Crépel, Ideal chern bands are landau levels in curved space (2023), arXiv:2304.01251 [cond-mat.mes-hall] .
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011), https://www.pnas.org/content/108/30/12233.full.pdf .
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Phys. Rev. B 86, 155449 (2012).
  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene bilayer with a twist: Electronic structure, Phys. Rev. Lett. 99, 256802 (2007).
  * Mele [2010] E. J. Mele, Commensuration and interlayer coherence in twisted bilayer graphene, Phys. Rev. B 81, 161405 (2010).
  * Vafek and Kang [2020] O. Vafek and J. Kang, Renormalization group study of hidden symmetry in twisted bilayer graphene with coulomb interactions, Phys. Rev. Lett. 125, 257602 (2020).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Phys. Rev. X 10, 031034 (2020).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Charged skyrmions and topological origin of superconductivity in magic-angle graphene, Science Advances 7, 10.1126/sciadv.abf5299 (2021).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. iii. interacting hamiltonian and exact symmetries, Phys. Rev. B 103, 205413 (2021).
  * Aroyo  _et al._ [2006a] M. I. Aroyo, J. M. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga, A. Kirov, and H. Wondratschek, Bilbao crystallographic server: I. databases and crystallographic computing programs, Zeitschrift für Kristallographie - Crystalline Materials 221, 15 (2006a).
  * Aroyo  _et al._ [2006b] M. I. Aroyo, A. Kirov, C. Capillas, J. M. Perez-Mato, and H. Wondratschek, Bilbao Crystallographic Server. II. Representations of crystallographic point groups and space groups, Acta Crystallographica Section A 62, 115 (2006b).
  * Ren  _et al._ [2021] Y. Ren, Q. Gao, A. H. MacDonald, and Q. Niu, Wkb estimate of bilayer graphene’s magic twist angles, Phys. Rev. Lett. 126, 016404 (2021).
  * Becker  _et al._ [2022] S. Becker, M. Embree, J. Wittsten, and M. Zworski, Mathematics of magic angles in a model of twisted bilayer graphene, Probability and Mathematical Physics 3, 69 (2022).
  * Navarro-Labastida and Naumis [2023] L. A. Navarro-Labastida and G. G. Naumis, 3/2 magic angle quantization rule of flat bands in twisted bilayer graphene and its relationship to the quantum hall effect, Phys. Rev. B 107, 155428 (2023).
  * Liu  _et al._ [2019] J. Liu, J. Liu, and X. Dai, Pseudo landau level representation of twisted bilayer graphene: Band topology and implications on the correlated insulating phase, Phys. Rev. B 99, 155415 (2019).
  * Haldane and Rezayi [1985] F. D. M. Haldane and E. H. Rezayi, Periodic laughlin-jastrow wave functions for the fractional quantized hall effect, Phys. Rev. B 31, 2529 (1985).
  * Haldane [2018] F. D. M. Haldane, A modular-invariant modified weierstrass sigma-function as a building block for lowest-landau-level wavefunctions on the torus, Journal of Mathematical Physics 59, 071901 (2018).
  * [62] DLMF, NIST Digital Library of Mathematical Functions, https://dlmf.nist.gov/, Release 1.1.10 of 2023-06-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Sheffer  _et al._ [2023] Y. Sheffer, R. Queiroz, and A. Stern, Symmetries as the guiding principle for flattening bands of dirac fermions, Phys. Rev. X 13, 021012 (2023).
  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Phys. Rev. B 55, 1142 (1997).