# Chiral limit of twisted trilayer graphene

Simon Becker  ETH Zurich, Institute for Mathematical Research, Rämistrasse
101, 8092 Zurich, Switzerland<EMAIL_ADDRESS>,  Tristan Humbert
ENS Paris, Département de Mathématiques et Applications, Rue d’Ulm, Paris,
France<EMAIL_ADDRESS>,  Jens Wittsten  Department of
Engineering, University of Borås, SE-501 90 Borås, Sweden<EMAIL_ADDRESS>and  Mengxuan Yang  Department of Mathematics, University of California,
Berkeley, CA 94720, USA<EMAIL_ADDRESS>
###### Abstract.

We initiate the mathematical study of the Bistritzer-MacDonald Hamiltonian for
twisted trilayer graphene in the chiral limit (and beyond). We develop a
spectral theoretic approach to investigate the presence of flat bands under
specific magic parameters. This allows us to derive trace formulae that show
that the tunnelling parameters that lead to flat bands are nowhere continuous
as functions of the twisting angles.

##  1\. Introduction

Figure 1. (Left): Magic parameters at which the chiral limit of the
Hamiltonian (1.1) exhibits a flat band. $x$-axis is ratio of twisting angles
$\frac{\zeta_{2}}{\zeta_{1}}.$ We assume $\alpha_{23}=\alpha_{21}$. Then
$y$-axis is real part of magic parameter $\operatorname{Re}(\alpha_{12})$ and
color coding is imaginary part $\operatorname{Im}(\alpha_{12}).$ (Right): Only
real magic parameters ($\operatorname{Im}(\alpha_{12})=0$).

Twisted trilayer graphene (TTG) is a material formed by stacking three sheets
of graphene, with slight relative twisting angles $\zeta_{1},\zeta_{2}$
between layers one and two, and two and three, respectively. This stacking
arrangement gives rise to a captivating interference pattern known as a moiré
pattern (see Figure 2). The electron tunnelling described by this pattern
leads to significant modifications in the electronic properties of the
material, such as flat bands at certain magic parameters, see Figure 1.

Each layer of graphene is composed of carbon atoms arranged in a two-
dimensional honeycomb lattice, featuring two distinct atom types, A and B, per
fundamental domain. At sites where neighboring layers (top and middle or
middle and bottom) align, inter-layer interactions between atoms of the same
type, referred to as AA sites, occur. Additionally, tunnelling interactions
between atoms of different types take place at AB and BA sites, where atoms of
type A are stacked over atoms of type B, and vice versa.

The exploration of different graphene layer configurations, including twisted
trilayer graphene, expands on the research conducted on twisted bilayer
graphene [TKV19, BEWZ21, BEWZ22, BHZ1-22, BHZ2-22, BHZ23]. These multilayer
systems offer increased tunability due to a larger set of parameters [KKAJ19,
LVK22]. We also want to mention recent progress on a helical version of
twisted trilayer graphene [Dev23, GMM23a, GMM23b]. The study of the continuum
or Bistritzer-MacDonald (BM) model for twisted trilayer graphene exhibits,
unlike the BM model for bilayer graphene [BM11], commensurable and
incommensurable twisting angles (see Assumption 1). Thus, TTG serves as an
example for a whole range of materials where commensurability matters. The
initial theoretical analysis of twisted trilayer graphene revealed a similar
phenomenon of electronic bands flattening at various _magic_ angles. This
breakthrough subsequently led to experimental observations of correlated
phenomena [PCWT21].

In this paper we study the BM model of twisted trilayer graphene described by
the Hamiltonian

$H(\alpha,\tilde{\alpha})=\begin{pmatrix}W(\tilde{\alpha})&D(\alpha)^{*}\\\
D(\alpha)&W(\tilde{\alpha})\end{pmatrix}$ (1.1)

with

$\begin{split}D(\alpha)&=\begin{pmatrix}2D_{\bar{z}}&\alpha_{12}U(pz)&0\\\
\alpha_{12}U(-pz)&2D_{\bar{z}}&\alpha_{23}U(p\tfrac{\zeta_{2}}{\zeta_{1}}z)\\\
0&\alpha_{23}U(-p\tfrac{\zeta_{2}}{\zeta_{1}}z)&2D_{\bar{z}}\end{pmatrix},\\\
W(\tilde{\alpha})&=\begin{pmatrix}0&\tilde{\alpha}_{12}V(pz)&0\\\
(\tilde{\alpha}_{12}V(pz))^{*}&0&\tilde{\alpha}_{23}V(p\frac{\zeta_{2}}{\zeta_{1}}z)\\\
0&(\tilde{\alpha}_{23}V(p\frac{\zeta_{2}}{\zeta_{1}}z))^{*}&0\end{pmatrix}.\end{split}$
(1.2)

The parameters $\alpha=(\alpha_{12},\alpha_{23})$ and
$\tilde{\alpha}=(\tilde{\alpha}_{12},\tilde{\alpha}_{23})$ describe (rescaled)
hopping amplitudes between layers $i$ and $j$ at AB/BA and AA sites,
respectively. We shall mainly focus on the _chiral limit_ , which discards the
tunnelling at AA sites by setting $\tilde{\alpha}=0$ in (1.1) and gives the
Hamiltonian

$H(\alpha):=H(\alpha,0)=\begin{pmatrix}0&D(\alpha)^{*}\\\
D(\alpha)&0\end{pmatrix}.$

The _anti-chiral limit_ is obtained by setting $\alpha=0$, instead. In Section
2 we provide a brief discussion and derivation of the Hamiltonian (1.1) from
the continuum model for twisted trilayer graphene. In (1.1), the parameters
$\zeta_{1}$ and $\zeta_{2}$ describe the relative (small) twisting angle
between layers. In Figure 2 we show examples of the moiré pattern formed in
TTG for different $\zeta_{1}$ and $\zeta_{2}$, and we will throughout the
paper assume that they are commensurable to analyze the electronic band
structure:

######  Assumption 1 (Commensurable angles).

We assume that $\frac{\zeta_{2}}{\zeta_{1}}\in\mathbb{Q}\setminus\\{0\\}$,
then we can write $\frac{\zeta_{2}}{\zeta_{1}}=3^{j}\frac{r_{1}}{r_{2}}$ with
$r_{1},j\in{\mathbb{Z}}$, $r_{2}\in\mathbb{N}$, and both
$r_{1},r_{2}\not\equiv 0\operatorname{mod}3$. We then set

$p:=\begin{cases}r_{2}&\text{ for }j>0,\\\ 3^{-j}r_{2}&\text{ for }j\leq 0,\\\
\end{cases}\quad\text{ and }\quad q:=\begin{cases}0&\text{ for }j>0,\\\
r_{1}&\text{ for }j\leq 0.\end{cases}$ (1.3)

Under Assumption 1, we have for $j\neq 0$ that the factors $p$ and
$p\frac{\zeta_{2}}{\zeta_{1}}$ in (1.2) satisfy
$p\frac{\zeta_{2}}{\zeta_{1}}\in 3{\mathbb{Z}}+q$ with precisely one of $p$
and $q$ in $3{\mathbb{Z}}$, and therefore also precisely one of $p$ and
$p\frac{\zeta_{2}}{\zeta_{1}}$ in $3{\mathbb{Z}}$, while for $j=0$ we have
neither $p$ nor $q=p\frac{\zeta_{2}}{\zeta_{1}}$ in $3{\mathbb{Z}}$. In
particular, at least one of $p$ and $q$ does not belong to $3{\mathbb{Z}}$,
and if $j\leq 0$ then $p\frac{\zeta_{2}}{\zeta_{1}}=q$.

Figure 2. Moiré patterns for different twisting configurations. The left panel
shows three layers stacked with equal relative twist angles of $4^{\circ}$
between each layer, while in the right panel, the relative twist angles
between layers are $4^{\circ}$ and $7^{\circ}$ degrees.

We introduce $\Gamma:=4\pi i(\omega{\mathbb{Z}}\oplus\omega^{2}{\mathbb{Z}})$,
where $\omega=e^{2\pi i/3}$ is a third root of unity. Then
$\Gamma_{3}:=\Gamma/3$ is the _moiré lattice_ , and $U,V\in
C^{\infty}({\mathbb{C}})$ in (1.2) satisfy

$\displaystyle U(z+na)$ $\displaystyle=\bar{\omega}^{n(a_{1}+a_{2})}U(z),\quad
U(\omega z)=\omega U(z),\quad\overline{U(\bar{z})}=U(z),$ (1.4) $\displaystyle
V(z+na)$ $\displaystyle=\bar{\omega}^{n(a_{1}+a_{2})}V(z),\quad V(\omega
z)=V(z),\quad\overline{V(z)}=V(-z)=V(\bar{z})$ (1.5)

for $n\in\mathbb{N}$ and $a=\frac{4\pi
i}{3}(a_{1}\omega+a_{2}\omega^{2})\in\Gamma_{3}$. In particular, $U$ and $V$
are periodic with respect to $\Gamma$. We provide a complete characterisation
of such potentials in Proposition 2.1. A standard example of $U$ in (1.4) is

$U_{0}(z)=\sum_{j=0}^{2}\omega^{j}e^{\frac{1}{2}(z\bar{\omega}^{j}-\bar{z}\omega^{j})}$
(1.6)

and it will serve as our reference potential for simulations, if nothing else
is mentioned.

Figure 3. Magic angles showing ratio of twisting angles $\zeta_{2}/\zeta_{1}$
($x$-axis) and $\operatorname{Re}(\alpha_{12}/\zeta_{2})$ ($y$-axis) with
color coding given by $\operatorname{Im}(\alpha_{12}/\zeta_{2}).$

We consider now the chiral Hamiltonian $H(\alpha)$, and the eponymous property

$H(\alpha)=-\sigma_{3}H(\alpha)\sigma_{3}.$ (1.7)

We define $\langle
k,z\rangle:=\operatorname{Re}(k\bar{z})=\frac{1}{2}(k\bar{z}+\bar{k}z)$. From
Assumption 1 it follows that the potentials $U(\pm pz)$ and $U(\pm
p\frac{\zeta_{2}}{\zeta_{1}}z)$ in the Hamiltonian
$H(\alpha,\widetilde{\alpha})$ in (1.1) are periodic up to a power of $\omega$
with respect to the moiré lattice $\Gamma_{3}$. Since
$H(\alpha,\widetilde{\alpha})$ commutes with the translation operator

$\mathscr{L}_{a}u(z):=\operatorname{diag}(\omega^{p(a_{1}+a_{2})},1,\bar{\omega}^{q(a_{1}+a_{2})},\omega^{p(a_{1}+a_{2})},1,\bar{\omega}^{q(a_{1}+a_{2})})u(z+a),\quad
a=\tfrac{4}{3}\pi i(\omega a_{1}+\omega^{2}a_{2})\in\Gamma_{3},$

we may enforce Floquet boundary conditions $\mathscr{L}_{a}u=e^{i\langle
a,k\rangle}u$ with $a\in\Gamma_{3}$ for arbitrary $k\in{\mathbb{C}}$, or
equivalently $k\in{\mathbb{C}}/\Gamma_{3}^{*}$. We can then conjugate the
Hamiltonian to obtain

$H_{k}(\alpha,\widetilde{\alpha}):=e^{-i\langle
k,z\rangle}H(\alpha,\widetilde{\alpha})e^{i\langle
k,z\rangle}:H^{1}_{\text{loc}}({\mathbb{C}})\to
L^{2}_{\text{loc}}({\mathbb{C}})$ (1.8)

so that in the chiral limit

$H_{k}(\alpha)=\begin{pmatrix}0&D(\alpha)^{*}+\bar{k}\\\
D(\alpha)+k&0\end{pmatrix}.$ (1.9)

The boundary condition for the operator $H_{k}$ reduces to
$\mathscr{L}_{a}\psi=\psi$ and we denote the associated $L^{2}$ space by
$L^{2}_{0}$ which is a subspace of $L^{2}_{\text{loc}}$ with the $L^{2}$ norm
on ${\mathbb{C}}/\Gamma_{3}$. The Floquet-Bloch decomposition of
$H(\alpha,\widetilde{\alpha})$ then implies that the spectrum can be described
as

$\operatorname{Spec}_{L^{2}({\mathbb{C}})}H(\alpha,\widetilde{\alpha})=\bigcup_{k\in{\mathbb{C}}/\Gamma_{3}^{*}}\operatorname{Spec}_{L^{2}_{0}}H_{k}(\alpha,\widetilde{\alpha})$

where $\Gamma_{3}^{*}$ is the dual lattice of $\Gamma_{3}$ (see (2.5) for a
definition).

For each $k$, the operator $H_{k}(\alpha)$ is an elliptic differential
operator with discrete spectrum. This yields a family of discrete spectra of
$\\{H_{k}(\alpha)\\}_{k\in{\mathbb{C}}/\Gamma_{3}^{*}}$ (cf. [BEWZ22, Section
2C]):

$\begin{gathered}\operatorname{Spec}_{L^{2}_{0}}H_{k}(\alpha)=\\{E_{\pm
j}(k,\alpha)\\}_{j=1}^{\infty},\quad E_{-j}(\alpha,k)=-E_{j}(\alpha,k),\\\
E_{j+1}(k,\alpha)\geq E_{j}(k,\alpha)\geq 0,\quad j\geq 1.\end{gathered}$
(1.10)

We refer to the image of $k\mapsto E_{j}(k,\alpha)$ as a band, and say that
$H(\alpha)$ has a _flat band_ at zero if $E_{1}(k,\alpha)=0$ for all
$k\in{\mathbb{C}}$. We then introduce the set

$\mathcal{A}:=\\{\alpha\in{\mathbb{C}}^{2}:E_{1}(k,\alpha)\equiv 0\text{ for
all }k\in{\mathbb{C}}\\},$ (1.11)

which we refer to as the set of _magic parameters_ or just _magic_ $\alpha$’s.
Note that by the above we have $\alpha\in\mathcal{A}$ if and only if
$0\in\operatorname{Spec}_{L^{2}_{0}}H_{k}(\alpha)$ for all $k\in{\mathbb{C}}$,
i.e.,

$\alpha\in\mathcal{A}\quad\Longleftrightarrow\quad\operatorname{Spec}_{L^{2}_{0}}(D(\alpha))={\mathbb{C}}.$
(1.12)

Just like in [BEWZ22], the key point to obtain a characterisation of magic
parameters is to establish the existence of _protected states_. Indeed, using
the symmetries of the Hamiltonian, we prove in Proposition 3.2 the existence
of three protected states

$\varphi(\alpha,z)\in\ker_{L^{2}_{p,0}}D(\alpha),\ \
\psi(\alpha,z)\in\ker_{L^{2}_{0,0}}D(\alpha),\ \
\rho(\alpha,z)\in\ker_{L^{2}_{-q,0}}D(\alpha),$ (1.13)

analytic in $\alpha$ with $\varphi(0,z)=e_{1}$, $\psi(0,z)=e_{2}$,
$\rho(0,z)=e_{3}$, where $\\{e_{1},e_{2},e_{3}\\}$ is the standard basis of
${\mathbb{C}}^{3}$.

Figure 4. Band structure for different $p,q$: Case I (left panel; $p\equiv
2,q\equiv 0$ mod 3), which shares many qualitative properties with TBG, where
at one Dirac point there is a two fold degenerated Dirac cone (protected
state) and at another Dirac point there is a simple Dirac cone (protected
state); Case II (right panel; $p\equiv q\equiv 1$ mod 3) where there exist
three Dirac cones (protected states) at different Dirac points.

Here, the subspaces $L^{2}_{m,j}$ are associated to irreducible
representations of the symmetry group, see Subsection 2.2.

From Proposition 3.2 we also obtain the following dichotomy for the three
protected states that we illustrate in Figure 4:

  1. Case I:

$-q,p,0$ are not mutually different numbers in $\mathbb{Z}_{3}$. Then two of
the protected states fall into the same subspace. This implies the existence
of overlapping Dirac cones in the Brillouin zone
${\mathbb{C}}/\Gamma_{3}^{*}$, see Figure 4. We also observe that, for
$p\not\equiv 0$ mod 3, we have $-p\in\mathbb{Z}_{3}\setminus\\{0,p,-q\\}$,
which we shall usually assume without loss of generality (see Remark 1).

  2. Case II:

$-q,p,0$ are mutually different numbers in $\mathbb{Z}_{3}$. Then

$\ker_{L^{2}_{p,0}}D(\alpha)\neq\\{0\\},\ \
\ker_{L^{2}_{0,0}}D(\alpha)\neq\\{0\\},\ \
\ker_{L^{2}_{-q,0}}D(\alpha)\neq\\{0\\}.$

This corresponds to three separate Dirac cones in Figure 4.

While Case I behaves in many ways similar to twisted bilayer graphene and most
of the result of [BHZ1-22, BHZ2-22, BHZ23] have an analogue to this case, Case
II exhibits many new phenomena. We note that equal twisted trilayer graphene
introduced in [PT23-1] falls into Case II.

Having these three protected states, we can define their Wronskian

$W(\alpha):=\det[\varphi,\psi,\rho],$ (1.14)

which vanishes if and only if $\alpha$ is magical. A similar treatment has
also been proposed by Popov–Tarnopolsky [PT23-2].

######  Theorem 1 (Flat bands & Wronskian).

The parameter $\alpha\in{\mathbb{C}}^{2}$ is magic if and only if the
Wronskian vanishes:

$\alpha\in\mathcal{A}\iff W(\alpha)=0.$

Following [BEWZ22] and [BHZ1-22], we will use an alternative characterisation,
for $\alpha_{12}\neq 0$, in terms of the compact and non-normal operator which
we call the _Birman-Schwinger operator_ defined by

$B_{k}(\tfrac{\alpha_{23}}{\alpha_{12}}):=R(k)U(-pz)R(k)U(pz)+\Big{(}\frac{\alpha_{23}}{\alpha_{12}}\Big{)}^{2}(R(k)U(p\tfrac{\zeta_{2}}{\zeta_{1}}z)R(k)U(-p\tfrac{\zeta_{2}}{\zeta_{1}}z))$
(1.15)

with $R(k)=(2D_{\bar{z}}+k)^{-1}$. The set of magic parameters is then
characterised by

######  Theorem 2 (Characterisation via the Birman-Schwinger operator).

The parameter $\alpha\in{\mathbb{C}}^{2}\setminus(\\{0\\}\times{\mathbb{C}})$
is magic if and only if
$\frac{1}{\alpha_{12}^{2}}\in\operatorname{Spec}(B_{k}(\tfrac{\alpha_{23}}{\alpha_{12}}))$
for $k\notin\Gamma^{*}.$

The above characterisation is explicit up to fixing the hopping ratio
$\tfrac{\alpha_{23}}{\alpha_{12}},$ i.e., the set of all magic parameters is a
union over all possible ratios $\tfrac{\alpha_{23}}{\alpha_{12}}.$ In
practice, one may just assume that $\alpha_{12}=\alpha_{23}$ for standard TTG.
The Birman-Schwinger operator $B_{k}$ is a Hilbert-Schmidt operator, which
implies that for $\ell\geq 2$ all traces $\operatorname{tr}(B_{k}^{2\ell})$
are well-defined and moreover

$\operatorname{tr}(B_{k}^{\ell})=\sum_{\alpha_{12}\in\mathcal{A}(\frac{\alpha_{23}}{\alpha_{12}})}\frac{1}{\alpha_{12}^{2\ell}},\quad\ell\geq
2$

where we denoted by $\mathcal{A}(\frac{\alpha_{23}}{\alpha_{12}})$ the set of
magic parameters for the fixed hopping ratio $\alpha_{23}/\alpha_{12}$ (see
§3.2 for a precise definition). We can therefore try to understand the
structure of the set of magic parameters
$\mathcal{A}(\frac{\alpha_{23}}{\alpha_{12}})$ by understanding the traces of
the operator $B_{k}$. Similarly to what was done in [BHZ1-22], we provide, in
Theorem 6, semi-explicit formulas for these traces. These expressions are
enough to prove the following result.

######  Theorem 3.

For the potential given in (2.4) with finite Fourier expansion with
coefficients in $c_{n,m}\in\mathbb{Q}(\omega)$ and with the additional
symmetry

$\overline{U(\bar{z})}=U(z)$

, we have for any $\ell\geq 2$ that
$\operatorname{tr}(B_{k}^{\ell})=q_{\ell}\pi/\sqrt{3}$, where
$q_{\ell}\in\mathbb{Q}$. As a consequence, if the set of magic parameters
$\mathcal{A}(\frac{\alpha_{23}}{\alpha_{12}})$ is non-empty, it is infinite.

We will actually prove statements for more general potentials in Theorems 7
and Theorem 8. As non-zero non-normal operators may have only zero in their
spectrum, the existence of a magic parameter from Theorem 2 is not trivial. In
[BEWZ22], the existence was established by computing explicitly the first
trace and proving it was non-zero. In this article, we provide an explicit
calculation of the first trace in Theorem 11 and, as a corollary, we get that
the set of magic parameters is infinite for any twisting angles satisfying
Assumption 1. If we are additionally in Case I, then by adapting the argument
in [BHZ23], we obtain the existence of infinitely many non-simple magic
parameters (see Proposition 3.4 for a precise definition of multiplicity).

######  Corollary 1.1.

Consider the potential $U_{0}$ defined in (1.6) and assume that
$\alpha_{23}/\alpha_{12}\in\mathbb{Q}$ is rational. Then the set of magic
parameters $\mathcal{A}$ is infinite. Assuming Case I, the set of non-simple
magic parameters is infinite.

We refer to Theorem 9 for the most general statement. Numerical experiments,
see Figure 6, suggest that the second part of Corollary 1.1 does not hold for
angles in Case II. Note that Corollary 1.1 implies that when the potential is
the standard one given by (1.6) and $\alpha_{23}/\alpha_{12}\in\mathbb{Q}$,
then there are infinitely many $\alpha$ for which the corresponding chiral
Hamiltonian $H(\alpha)$ has a flat band. This is in contrast to the anti-
chiral model which never exhibits flat bands, see Appendix A.

We then investigate the continuity of magic parameters. Our main result shows
that the set of magic parameters is maximally discontinuous when changing
continuously the twisting angles $\zeta_{1},\zeta_{2}$. By this we mean that
the lowest order trace function computed in the preceding section is
_discontinuous at any point_ in the following sense:

######  Theorem 4.

Let $\zeta_{1}$ and $\zeta_{2}$ with $\zeta_{1}/\zeta_{2}\in\mathbb{Q}$ fixed.
For $r=\alpha_{23}/\alpha_{12}$ fixed, set the trace of unrescaled magic
parameters to be

$\mathcal{S}_{4}(r,\zeta_{1},\zeta_{2}):=\bigg{(}\frac{p}{\zeta_{1}}\bigg{)}^{4}\sum_{\alpha_{12}\in\mathcal{A}(r)}\frac{1}{\alpha_{12}^{4}}.$

Then there exists a sequence $(\zeta_{1}^{(n)},\zeta_{2}^{(n)})$ such that for
any $n$, one has $\zeta_{1}^{(n)}/\zeta_{2}^{(n)}\in\mathbb{Q}$ with
$(\zeta_{1}^{(n)},\zeta_{2}^{(n)})\to(\zeta_{1},\zeta_{2})$. Moreover, one has

$\mathcal{S}_{4}(r,\zeta_{1}^{(n)},\zeta_{2}^{(n)})\to+\infty.$

In particular, the trace $\mathcal{S}_{4}(r,\zeta_{1}^{(n)},\zeta_{2}^{(n)})$
is sequentially discontinuous at _any_ point.

In Section 7.1 we show that for a generic choice of potentials (7.1) we can
ensure in Case I that all magic angles are simple within each subspace
$L^{2}_{-p,\ell}$ for $\ell\in{\mathbb{Z}}_{3}.$ If we then choose $V$ as in
(7.3) then the sublattice polarized flat bands will generically have Chern
number $-1$ which is shown in Section 9. Here, we find a very general argument
that should be useful in general for the study of twisted $N$-layer systems.
It is interesting to compare with twisted bilayer graphene where additional
symmetries allowed us to obtain a stronger result under a less restrictive
class of perturbations [BHZ23].

While most of our analysis so far mostly focus on the operator $D(\alpha)$, we
will get back to the full Hamiltonian $H(\alpha)$ in Section 8 by studying its
band structure when $\alpha$ is supposed to be simple, i.e.,
$E_{1}(k,\alpha)=0$ for any $k$ but $E_{2}(k_{0},\alpha)\neq 0$ for some
$k_{0}$. This leads naturally to the question of band touching: namely, does
there exist $k_{1}$ such that $E_{2}(k_{1},\alpha)=0$? In the case of twisted
bilayer graphene, it was shown in [BHZ23, Theo. 4] that there was a spectral
gap for a simple or doubly-degenerate magic angle. Interestingly, in the case
of twisted trilayer graphene, the first two bands always touch. Nevertheless,
if $\alpha$ is simple then the point $K_{0}$ where the two bands touch is
unique (see Section 8 for a more precise statement about the value of
$K_{0}$).

We close out the paper by studying exponential squeezing of bands for the
chiral model: In Theorem 15 in Section 10 we show that as the twisting angles
between the two lattices approach zero, a large number of bands accumulate in
an exponentially small neighbourhood around zero, determined by the twisting
angle ratio. In particular, we observe that even though exponential, the rate
of convergence depends significantly on the ratio of the twisting angles (see
Figure 12).

Notation. We denote by $\sigma_{i}$ the $i$-th Pauli matrix. Let
$\Gamma\subset\mathbb{C}$ be the lattice $\Gamma=4\pi
i(\omega\mathbb{Z}\oplus\omega^{2}\mathbb{Z})$ where $\omega=e^{2\pi i/3}$ is
a third root of unity. The dual lattice $\Gamma^{*}$ is defined as the
$k\in\mathbb{C}$ such that
$\langle\gamma,k\rangle:=\frac{1}{2}(\gamma\bar{k}+\bar{\gamma}k)\in
2\pi\mathbb{Z}$ for all $\gamma\in\Gamma$, and is given by
$\Gamma^{*}=\frac{1}{\sqrt{3}}(\omega\mathbb{Z}\oplus\omega^{2}\mathbb{Z})$.
We will represent $k\in\Gamma^{*}$ as
$k=k(k_{1},k_{2})=\frac{1}{\sqrt{3}}(\omega^{2}k_{1}-\omega k_{2})$. The moiré
lattice and its dual lattice are given by

$\displaystyle\Gamma_{3}:=\tfrac{1}{3}\Gamma=\tfrac{4}{3}\pi i(\omega
a_{1}+\omega^{2}a_{2}),\ a_{1},a_{2}\in{\mathbb{Z}};\ \
\Gamma_{3}^{*}=3\Gamma^{*}=\sqrt{3}(\omega^{2}k_{1}-\omega k_{2}),\
k_{1},k_{2}\in{\mathbb{Z}}.$

Outline of article.

  * •

In Section 2, we derive the full Hamiltonian $H(\alpha,\tilde{\alpha})$ in
(1.1) and study its basic symmetries.

  * •

In Section 3, we discuss the protection of Dirac points in the chiral limit
and its consequences on magic parameters. Then, we state the Birman-Schwinger
principle for twisted trilayer graphene and obtain the scalar operator (1.15)
whose eigenvalues characterize _magic parameters_.

  * •

In Section 4, we study structural properties of trace formulæ that we obtain
from the Birman-Schwinger operator.

  * •

In Section 5, we explicitly work out lowest order trace formulæ that show the
existence of infinitely many magic parameters only assuming Assumption 1.

  * •

In Section 6, we prove the discontinuity of magic parameters with respect to
the twisting angle and the Hölder continuity of the spectrum of the
Hamiltonian in Hausdorff distance.

  * •

In Section 7, we study the generic simplicity of flat bands in different
representations and state conditions for degenerate magic parameters to exist.

  * •

In Section 8, we discuss the theta function construction to construct Bloch
functions associated with flat bands.

  * •

In Section 9, we use results from Section 7 and 8 to compute the Chern number
of the flat bands for a variety of cases.

  * •

In Section 10, we study exponential squeezing of bands in the limit of small
twisting angles.

  * •

In Appendix A, we study the anti-chiral limit ($\alpha\equiv 0$) of twisted
trilayer graphene and show that it does not exhibit any flat bands.

##  2\. The twisted trilayer graphene Hamiltonian

###  2.1. Derivation

We start with the continuum model for twisted trilayer graphene with coupling
parameters $\beta,\tilde{\beta}\in{\mathbb{C}}^{2}$,

$\begin{split}&\mathscr{H}(\beta,\tilde{\beta})=\begin{pmatrix}H_{D}&T_{12}(\zeta_{1}z)&0\\\
T_{12}(\zeta_{1}z)^{*}&H_{D}&T_{23}(\zeta_{2}z)\\\
0&T_{23}(\zeta_{2}z)^{*}&H_{D}\end{pmatrix}\text{ with
}T_{ij}(z)=\begin{pmatrix}\tilde{\beta}_{ij}V(z)&(\beta_{ij}U(-z))^{*}\\\
\beta_{ij}U(z)&\tilde{\beta}_{ij}V(z)\end{pmatrix}.\end{split}$

Here, $H_{D}$ is a Dirac operator
$H_{D}=D_{x_{1}}\sigma_{1}+D_{x_{2}}\sigma_{2}$, and $\mathscr{H}$ acts on a
wavefunction

$\psi=(\psi_{\uparrow,A},\psi_{\uparrow,B},\psi_{\rightarrow,A},\psi_{\rightarrow,B},\psi_{\downarrow,A},\psi_{\downarrow,B})$

with lattice $\uparrow,\rightarrow,\downarrow$ (upper, middle, lower) and
sublattice index $A,B$. Writing this operator in the basis in which
$\psi=(\psi_{\uparrow,A},\psi_{\rightarrow,A},\psi_{\downarrow,A},\psi_{\uparrow,B},\psi_{\rightarrow,B},\psi_{\downarrow,B})$,
we see it is equivalent to a Hamiltonian

$\mathcal{H}(\beta,\tilde{\beta},\zeta)=\begin{pmatrix}\mathcal{W}(\tilde{\beta},\zeta)&\mathcal{D}(\beta,\zeta)^{*}\\\
\mathcal{D}(\beta,\zeta)&\mathcal{W}(\tilde{\beta},\zeta)\end{pmatrix}$

with matrices

$\begin{split}\mathcal{D}(\beta,\zeta)&=\begin{pmatrix}2D_{\bar{z}}&\beta_{12}U(\zeta_{1}z)&0\\\
\beta_{12}U(-\zeta_{1}z)&2D_{\bar{z}}&\beta_{23}U(\zeta_{2}z)\\\
0&\beta_{23}U(-\zeta_{2}z)&2D_{\bar{z}}\end{pmatrix}\text{ and }\\\
W(\tilde{\beta},\zeta)&=\begin{pmatrix}0&\tilde{\beta}_{12}V(\zeta_{1}z)&0\\\
(\tilde{\beta}_{12}V(\zeta_{1}z))^{*}&0&\tilde{\beta}_{23}V(\zeta_{2}z)\\\
0&(\tilde{\beta}_{23}V(\zeta_{2}z))^{*}&0\end{pmatrix},\end{split}$ (2.1)

where $\zeta_{k}=\theta_{k+1}-\theta_{k}$ in the small-angle approximation,
with $\theta_{1},\theta_{2},\theta_{3}\in(0,2\pi)$ denoting the rotation angle
of the three layers, respectively. We then choose $p$ in accordance with
(1.3), and make the change of variables $z\zeta_{1}/p\mapsto z$. By
introducing new coupling parameters
$\alpha_{ij}:=\frac{p}{\zeta_{1}}\beta_{ij}$ and
$\tilde{\alpha}_{ij}:=\frac{p}{\zeta_{1}}\tilde{\beta}_{ij}$, we obtain a new
equivalent Hamiltonian $H(\alpha,\tilde{\alpha})$ (with rescaled energy scale)
given by (1.1). As mentioned in the introduction, the chiral limit is obtained
by setting $\tilde{\alpha}=0$.

######  Remark 1.

We may without loss of generality assume that $p\not\equiv 0$ mod 3, since
otherwise we can just flip the trilayer upside down to effectively replace $p$
with $\lvert q\rvert\not\equiv 0$ mod 3 as follows: Assume that $p\equiv 0$
mod 3. Then, in the notation of Assumption 1, we have
$\zeta_{2}/\zeta_{1}=3^{j}r_{1}/r_{2}$ with $j<0$ and $q=r_{1}\not\equiv 0$
mod 3. Conjugating $\mathcal{D}(\beta,\zeta)$ by a unitary matrix, we get

$\begin{pmatrix}0&0&1\\\ 0&1&0\\\
1&0&0\end{pmatrix}\mathcal{D}(\beta,\zeta)\begin{pmatrix}0&0&1\\\ 0&1&0\\\
1&0&0\end{pmatrix}=\begin{pmatrix}2D_{\bar{z}}&\beta_{23}U(-\zeta_{2}z)&0\\\
\beta_{23}U(\zeta_{2}z)&2D_{\bar{z}}&\beta_{12}U(-\zeta_{1}z)\\\
0&\beta_{12}U(\zeta_{1}z)&2D_{\bar{z}}\end{pmatrix},$

which swaps $\zeta_{1}$ for $-\zeta_{2}$ (as well as $\beta_{12}$ for
$\beta_{23}$), and we have
$\zeta_{1}/\zeta_{2}=3^{-j}r_{2}\operatorname{sgn}(r_{1})/\lvert r_{1}\rvert$
with $\tilde{j}=-j>0$, $\operatorname{sgn}(r_{1})r_{2}\in{\mathbb{Z}}$ and
$\lvert r_{1}\rvert\in\mathbb{N}$. Following the convention in Assumption 1,
we set $\tilde{p}:=\lvert r_{1}\rvert=\lvert q\rvert\not\equiv 0$ mod 3, and
$\tilde{q}=0$. After making the change of variables
$-\zeta_{2}z\mapsto\tilde{p}z$ and setting
$(\alpha_{12},\alpha_{23})=(\beta_{23},\beta_{12})\tilde{p}/(-\zeta_{2})$, we
end up with the equivalent system

$\begin{pmatrix}2D_{\bar{z}}&\alpha_{12}U(\tilde{p}z)&0\\\
\alpha_{12}U(-\tilde{p}z)&2D_{\bar{z}}&\alpha_{23}U(\tilde{p}\frac{\zeta_{1}}{\zeta_{2}}z)\\\
0&\alpha_{23}U(-\tilde{p}\frac{\zeta_{1}}{\zeta_{2}}z)&2D_{\bar{z}}\end{pmatrix}$

where $\tilde{p}\not\equiv 0$ mod 3 and
$\tilde{p}\zeta_{1}/\zeta_{2}=3^{-j}\operatorname{sgn}(r_{1})r_{2}\in
3{\mathbb{Z}}+\tilde{q}$, which proves the claim.

The potentials $U$ and $V$ in (1.2) and (2.1) are smooth functions that
satisfy (1.4)–(1.5). A characterization of such functions $U$ and $V$ is given
in the following proposition. A straightforward proof (which we omit for
brevity) can be obtained by using the fact that a function $u\in
L^{2}({\mathbb{C}}/\Gamma;{\mathbb{C}})$ can be expanded in a Fourier series
with respect to the orthonormal basis

$e_{k}(z):=e^{\frac{i}{2}(z\bar{k}+\bar{z}k)}/\mathrm{Vol}(\mathbb{C}/\Gamma)^{\frac{1}{2}},\quad
k\in\Gamma^{*},$ (2.2)

namely,

$u(z)=\sum_{k\in\Gamma^{*}}c_{k}e_{k}(z),\quad
c_{k}=\int_{{\mathbb{C}}/\Gamma}u(z)\overline{e_{k}(z)}\,dA(z),$

where $dA$ is the Haar measure on the torus $\mathbb{C}/\Gamma$.

######  Proposition 2.1.

Let $a=\frac{4\pi i}{3}(\omega a_{1}+\omega^{2}a_{2})\in\Gamma_{3}$. Let
$j,\ell\in\\{-1,0,1\\}$. Then $u\in L^{2}(\mathbb{C}/\Gamma)$ satisfies

$u(z+a)=\omega^{j(a_{1}+a_{2})}u(z),\quad u(\omega
z)=\bar{\omega}^{\ell}u(z),$ (2.3)

if and only if

$u(z)=\sum_{n,m\in\mathbb{Z}}c_{nm}e^{\frac{i}{2}(-\sqrt{3}(n+m)\operatorname{Re}z-(2j+3(n-m))\operatorname{Im}z)}$
(2.4)

where

$c_{nm}=(u,\phi_{k})_{L^{2}({\mathbb{C}}/\Gamma)},\quad
k={\tfrac{1}{\sqrt{3}}(\omega^{2}(j+3n)-\omega(j-3m))},\quad
n,m\in{\mathbb{Z}},$

satisfies
$c_{nm}=\bar{\omega}^{\ell}c_{(-m)(n-m+j)}=\bar{\omega}^{2\ell}c_{(m-n-j)(-n)}$
for $n,m\in\mathbb{Z}$. If in addition $\overline{u(\bar{z})}=u(z)$ then
$\overline{c_{nm}}=c_{(-m)(-n)}=\bar{\omega}^{\ell}c_{n(n-m+j)}=\bar{\omega}^{2\ell}c_{(m-n-j)m}$
for $n,m\in\mathbb{Z}$, and if also $u(\bar{z})=u(-z)$ (which is only
compatible with $\ell=0$) then $c_{nm}$ is real for all $n,m$.

The conditions on $U$ in (1.4) correspond to taking $j=\ell=-1$ in Proposition
2.1. We illustrate our results for $U$ obtained by choosing $c_{00}=1$, which
then requires $c_{10}=\omega$ and $c_{0(-1)}=\omega^{2}$, and setting all
other coefficients to 0. Noting that $z-\bar{z}=2i\operatorname{Im}z$, while
$z\bar{\omega}-\bar{z}\omega=-i\sqrt{3}\operatorname{Re}z-i\operatorname{Im}z$
and
$z\bar{\omega}^{2}-\bar{z}\omega^{2}=i\sqrt{3}\operatorname{Re}z-i\operatorname{Im}z$,
this potential is given by (1.6). More generally,

$f_{3n+1}(z)=\sum_{j=0}^{2}\omega^{j}e^{\frac{3n+1}{2}(z\bar{\omega}^{j}-\bar{z}\omega^{j})}$

also satisfies (2.3) with $j=\ell=-1$, and so do the linear combinations
$U(z)=\sum_{n\in\mathbb{Z}}b_{3n+1}f_{3n+1}(z)$ that appeared in [BEWZ22].
Observe however that not all potentials in (2.4) for $j=\ell=-1$ are of this
form.

###  2.2. Symmetries of the full Hamiltonian

We will now discuss the symmetries of the twisted trilayer graphene
Hamiltonian (1.1).

First define the following twisted translations

$\mathscr{L}_{a}u(z)=\operatorname{diag}(\omega^{p(a_{1}+a_{2})},1,\bar{\omega}^{q(a_{1}+a_{2})})u(z+a),\
\ a=\tfrac{4}{3}\pi i(\omega a_{1}+\omega^{2}a_{2})\in\Gamma_{3}.$

We then have $[D(\alpha),\mathscr{L}_{a}]=0$ for $D(\alpha)$ defined in (1.2).
Consequently, for $u\in
L^{2}_{\operatorname{loc}}({\mathbb{C}};{\mathbb{C}}^{6})$, the Hamiltonian
commutes with two symmetries

$\begin{split}\mathscr{C}u(z)&:=\operatorname{diag}(1,1,1,\bar{\omega},\bar{\omega},\bar{\omega})u(\omega
z)\\\
\mathscr{L}_{a}u(z)&:=\operatorname{diag}(\omega^{p(a_{1}+a_{2})},1,\bar{\omega}^{q(a_{1}+a_{2})},\omega^{p(a_{1}+a_{2})},1,\bar{\omega}^{q(a_{1}+a_{2})})u(z+a),\end{split}$

such that $\mathscr{C}\mathscr{L}_{a}=\mathscr{L}_{\bar{\omega}a}\mathscr{C}$.
We slightly abuse notation by using $\mathscr{L}_{a}$ and $\mathscr{C}$ on
$L^{2}({\mathbb{C}};{\mathbb{C}}^{6})$ as well as on
$L^{2}({\mathbb{C}};{\mathbb{C}}^{3})$; for $\mathscr{C}$ this means
identifying $u\in L^{2}({\mathbb{C}};{\mathbb{C}}^{3})$ with
$(u,0_{{\mathbb{C}}^{3}})\in L^{2}({\mathbb{C}};{\mathbb{C}}^{6})$.

We introduce

$L^{2}_{k}:=\\{u\in
L^{2}_{\operatorname{loc}}({\mathbb{C}};{\mathbb{C}}^{d}):\mathscr{L}_{a}u=e^{i\langle
a,k\rangle}u,\ a\in\Gamma_{3}\\},\ \ k\in{\mathbb{C}},$

where $\langle a,k\rangle=\tfrac{1}{2}(a\bar{k}+\bar{a}k)$. Note that the
spaces $L^{2}_{k}$ only depend on the equivalence class of
$k\in{\mathbb{C}}/\Gamma_{3}^{*}$, where the dual lattice $\Gamma_{3}^{*}$ is
given by

$\Gamma_{3}^{*}=\\{{\sqrt{3}}(\omega^{2}k_{1}-\omega k_{2}),\
k_{1},k_{2}\in{\mathbb{Z}}\\}.$ (2.5)

We define a map between the spaces
$\\{L^{2}_{k}\\}_{k\in{\mathbb{C}}/\Gamma_{3}^{*}}$:

$\tau(k):L^{2}_{k^{\prime}}\to
L^{2}_{k^{\prime}+k},\quad\tau(k)u(z)=e^{i\langle z,k\rangle}u(z).$ (2.6)

Among the spaces $\\{L^{2}_{k}:k\in{\mathbb{C}}\\}$, we focus on three
distinguished
$k\in\mathcal{K}:=\\{\tfrac{r}{\sqrt{3}}(\omega^{2}-\omega),r\in{\mathbb{Z}}_{3}\\}=\\{0,\pm
i\\}\subset{\mathbb{C}}/\Gamma_{3}^{*}$, as they are the only possible
locations of protected states of the Hamiltonian (cf. Section 3). For
$k\in\mathcal{K}$, in view of the equation

$e^{i\langle a,k\rangle}=\omega^{r(a_{1}+a_{2})},\ \
k=\tfrac{r}{\sqrt{3}}(\omega^{2}-\omega),\ \ r\in{\mathbb{Z}}_{3},\ \
a\in\Gamma_{3},$

we define the subspaces of $L^{2}_{k}$ under the action of $\mathscr{C}$ which
correspond to irreducible representations of the symmetry group (see [BEWZ22,
Section 2.2]),

$L^{2}_{{r,\ell}}=\\{u\in
L^{2}_{\operatorname{loc}}({\mathbb{C}};{\mathbb{C}}^{3}):\mathscr{L}_{a}u=\omega^{r(a_{1}+a_{2})}u,\
\ \mathscr{C}u=\bar{\omega}^{\ell}u,\ \ a\in\Gamma_{3}\\},\ \
r,\ell\in{\mathbb{Z}}_{3}.$

Then for $k\in\mathcal{K}$, we have (see [BHZ2-22, Section 2.1])

$\displaystyle L^{2}_{k}=\bigoplus_{\ell\in{\mathbb{Z}}_{3}}L^{2}_{{r,\ell}},\
\ k=\tfrac{r}{\sqrt{3}}(\omega^{2}-\omega),\ \ r\in{\mathbb{Z}}_{3}.$ (2.7)

######  Remark 2.

A simple computation gives

$\mathcal{K}=\\{-ir:r\in{\mathbb{Z}}_{3}\\},$ (2.8)

thus the correspondence (2.7) can be reformulated as

$L^{2}_{-ir}=\bigoplus_{\ell\in{\mathbb{Z}}_{3}}L^{2}_{{r,\ell}},\ \
r\in{\mathbb{Z}}_{3}.$ (2.9)

Indeed, the points in $\mathcal{K}$ are fixed under rotation by $\omega$
modulo $\Gamma_{3}^{*}$. These three distinguished points in
${\mathbb{C}}/\Gamma_{3}^{*}$ are the locations of protected states in Figure
4111The figure is shown over a fundamental domain of
${\mathbb{C}}/\Gamma_{3}^{*}$ represented in rectangular coordinates
$y=(y_{1},y_{2})$ defined via $z=2i(\omega y_{1}+\omega^{2}y_{2})$, where the
distinguished points $i$, $0$, and $-i$ correspond to $(-1,-1)$, $(0,0)$, and
$(1,1)$ respectively..

Define the involutive $\mathcal{P}\mathcal{T}$ symmetry and the involutive
mirror symmetry $\mathscr{M}$ by

$(\mathcal{P}\mathcal{T}u)(z)=\begin{pmatrix}0&\operatorname{id}\\\
\operatorname{id}&0\end{pmatrix}\overline{u(-z)},\quad\mathscr{M}u(z)=\begin{pmatrix}0&\mathscr{E}\\\
\mathscr{E}&0\end{pmatrix}u(-\bar{z}),$

with $\mathscr{E}=\operatorname{diag}(1,-1,1)$. Using the symmetry properties
(1.4)–(1.5) of $U$ and $V$ it is then easy to check (assuming
$\tilde{\alpha}\in{\mathbb{R}}^{2}$) that

$(\mathcal{P}\mathcal{T})H(\alpha,\tilde{\alpha})(\mathcal{P}\mathcal{T})=H(\alpha,\tilde{\alpha}),\quad(\mathcal{P}\mathcal{T})\mathscr{C}=\omega\mathscr{C}(\mathcal{P}\mathcal{T})\quad\text{and}\quad(\mathcal{P}\mathcal{T})\mathscr{L}_{a}=\mathscr{L}_{-a}(\mathcal{P}\mathcal{T})$

and

$\mathscr{M}H(\alpha,\tilde{\alpha})\mathscr{M}=-H(\alpha,\tilde{\alpha}),\quad\mathscr{M}\mathscr{C}=\bar{\omega}\mathscr{C}^{-1}\mathscr{M},\quad\text{and}\quad\mathscr{M}\mathscr{L}_{a}=\mathscr{L}_{-\bar{a}}\mathscr{M}.$

Then the joint application of mirror and $\mathcal{P}\mathcal{T}$ symmetry
leave $L^{2}_{j,\ell}$ spaces invariant:

$L^{2}_{j,\ell}\xrightarrow[]{{\mathscr{M}}}L^{2}_{j,-\ell+1}\xrightarrow[]{{\mathcal{P}\mathcal{T}}}L^{2}_{j,\ell}.$

Indeed, if $u\in L^{2}_{j,\ell}$ then $\mathscr{M}u\in L^{2}_{j,-\ell+1}$
since

$\mathscr{L}_{a}\mathscr{C}(\mathscr{M}u)=\bar{\omega}\mathscr{M}(\mathscr{L}_{-\bar{a}}\mathscr{C}^{-1}u)=\omega^{j(a_{1}+a_{2})}\bar{\omega}^{-\ell+1}\mathscr{M}u,$

and if $v\in L^{2}_{j,-\ell+1}$ then $\mathcal{P}\mathcal{T}v\in
L^{2}_{j,\ell}$ since

$\mathscr{L}_{a}\mathscr{C}(\mathcal{P}\mathcal{T}v)=\bar{\omega}\mathcal{P}\mathcal{T}(\mathscr{L}_{-a}\mathscr{C}v)=\bar{\omega}\mathcal{P}\mathcal{T}(\omega^{-j(a_{1}+a_{2})}\bar{\omega}^{-\ell+1}v)=\omega^{j(a_{1}+a_{2})}\bar{\omega}^{\ell}\mathcal{P}\mathcal{T}v.$

In particular, defining the map

$Q:L_{\operatorname{loc}}^{2}({\mathbb{C}};{\mathbb{C}}^{3})\to
L_{\operatorname{loc}}^{2}({\mathbb{C}};{\mathbb{C}}^{3})\text{ with
}Qv(z)=\overline{v(-z)},$ (2.10)

one can check that $QD(\alpha)Q=D(\alpha)^{*}$ and $Q:L^{2}_{j,\ell}\to
L^{2}_{j,-\ell}$. The symmetries can then give the symmetry of the spectrum of
$H$ in each representation.

######  Proposition 2.2.

The densely defined operator $H(\alpha,\tilde{\alpha}):L^{2}_{j,\ell}\to
L^{2}_{j,\ell}$, $j,\ell\in{\mathbb{Z}}_{3}$ satisfies

$\operatorname{Spec}_{L^{2}_{j,\ell}}H(\alpha,\tilde{\alpha})=-\operatorname{Spec}_{L^{2}_{j,\ell}}H(\alpha,\tilde{\alpha}).$

###### Proof.

This follows from conjugating $H(\alpha,\tilde{\alpha})$ by
$\mathscr{M}(\mathcal{P}\mathcal{T}).$ ∎

###  2.3. Floquet theory

Since $U$ is periodic with respect to $\Gamma,$ we can also perform the Bloch-
Floquet transform on $H$ in (1.1) with respect to $\Gamma$ and standard
$\Gamma$-lattice translations. This way, we obtain the same fibre operators
(1.8) but with periodic boundary conditions on ${\mathbb{C}}/\Gamma.$ Thus, we
have the equality of sets

$\operatorname{Spec}_{L^{2}({\mathbb{C}})}H(\alpha,\widetilde{\alpha})=\bigcup_{k\in{\mathbb{C}}/\Gamma_{3}^{*}}\operatorname{Spec}_{L^{2}_{0}}H_{k}(\alpha,\widetilde{\alpha})=\bigcup_{k\in{\mathbb{C}}/\Gamma^{*}}\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}H_{k}(\alpha,\widetilde{\alpha}).$
(2.11)

However, working on $\Gamma$ rather than $\Gamma_{3}$ has the disadvantage
that, due to the larger lattice, it artificially inflates the number of bands
by a factor $9$, since
$\lvert{\mathbb{C}}/\Gamma\rvert=9\lvert{\mathbb{C}}/\Gamma_{3}\rvert$. This
is thoroughly explained in [BHZ2-22, (2.6)], see also [BHZ2-22, Figure $4$].
We shall thus mostly refrain from working on $\Gamma$ when referring to flat
bands. One thing we do want to mention is that one can argue analogously to
(1.12) and conclude that

$\alpha\in\mathcal{A}\quad\Longleftrightarrow\quad\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}(D(\alpha))={\mathbb{C}}.$
(2.12)

##  3\. Protected states and a Birman-Schwinger principle

In this section, we consider the Hamiltonian $H(\alpha)$ in the chiral limit
$\tilde{\alpha}=0$ satisfying equation (1.7) on $L^{2}({\mathbb{C}}/\Gamma)$.
The main purpose is to prove the existence of protected states at zero energy
(see Figure 4) in Proposition 3.2, and to establish a spectral
characterisation in Theorem 5 of the set of magic parameters $\mathcal{A}$
defined in (1.11), as well as a characterisation of $\mathcal{A}$ in terms of
the Wronskian of protected states, see Theorem 1.

###  3.1. Protected states at zero energy

For $\alpha_{12}=\alpha_{23}=0$, we have
$\operatorname{ker}_{L^{2}({\mathbb{C}}/\Gamma)}H(0)=\operatorname{Span}\\{e_{1},\ldots,e_{6}\\}$
by Liouville’s theorem, where $\\{e_{j}\\}$ is the standard basis on
${\mathbb{C}}^{6}$. Moreover, for $a\in\Gamma_{3}$,

$\begin{split}\mathscr{L}_{a}e_{1}&=\omega^{p(a_{1}+a_{2})}e_{1},\quad\mathscr{L}_{a}e_{2}=e_{2},\quad\mathscr{L}_{a}e_{3}=\bar{\omega}^{q(a_{1}+a_{2})}e_{3},\\\
\mathscr{L}_{a}e_{4}&=\omega^{p(a_{1}+a_{2})}e_{4},\quad\mathscr{L}_{a}e_{5}=e_{5},\quad\mathscr{L}_{a}e_{6}=\bar{\omega}^{q(a_{1}+a_{2})}e_{6},\\\
\mathscr{C}e_{1}&=e_{1},\quad\mathscr{C}e_{2}=e_{2},\quad\mathscr{C}e_{3}=e_{3},\quad\mathscr{C}e_{4}=\bar{\omega}e_{4},\quad\mathscr{C}e_{5}=\bar{\omega}e_{5},\quad\mathscr{C}e_{6}=\bar{\omega}e_{6},\end{split}$

which implies that

$\begin{split}e_{1}\in L^{2}_{p,0},\ e_{2}\in L^{2}_{0,0},\ e_{3}\in
L^{2}_{-q,0},\ e_{4}\in L^{2}_{p,1},\ e_{5}\in L^{2}_{0,1},\ e_{6}\in
L^{2}_{-q,1}.\end{split}$

By Assumption 1, there is at least one element in $\\{p,-q,0\\}$ which is
distinct from the others. Hence, there is an $n\in\\{p,-q,0\\}$ such that

$\dim\operatorname{ker}_{L^{2}_{n,0}}H(0,0)=\dim\operatorname{ker}_{L^{2}_{n,1}}H(0,0)=1.$

It follows from Proposition 2.2 that

$\operatorname{ker}_{L^{2}_{n,0}}H(\alpha,\widetilde{\alpha})\neq\emptyset,\
\operatorname{ker}_{L^{2}_{n,1}}H(\alpha,\widetilde{\alpha})\neq\emptyset,$
(3.1)

for all $\alpha\in{\mathbb{C}}^{2}$, $\widetilde{\alpha}\in{\mathbb{R}}^{2}$
as $H(\alpha)$ depends analytically on
$\alpha_{12},\alpha_{23}\in{\mathbb{C}}$ so that the one-dimensional
eigenspace in $\ker_{L^{2}_{n,0}}H(\alpha,\widetilde{\alpha})$ and
$\ker_{L^{2}_{n,1}}H(\alpha,\widetilde{\alpha})$ cannot split for $\alpha\neq
0$ or $\widetilde{\alpha}\neq 0$. In the chiral limit, $\widetilde{\alpha}=0$,
since
$\operatorname{ker}_{L^{2}_{r,\ell}}H(\alpha)=\operatorname{ker}_{L^{2}_{r,\ell}}D(\alpha)\oplus\\{0_{{\mathbb{C}}^{3}}\\}+\\{0_{{\mathbb{C}}^{3}}\\}\oplus\operatorname{ker}_{L^{2}_{r,\ell}}D(\alpha)^{*}$
we find that for this $n$ we have for all $\alpha\in{\mathbb{C}}^{2}$ that

$\operatorname{ker}_{L^{2}_{n,0}}D(\alpha)\neq\\{0\\}.$ (3.2)

Note also that (3.1) implies the following:

######  Corollary 3.1.

For all $\alpha\in\mathbb{C}^{2},\widetilde{\alpha}\in\mathbb{R}^{2}$,
$0\in\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}H(\alpha,\widetilde{\alpha})$.
In particular,
$0\in\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}H(\alpha).$

We now focus on the chiral limit. The following proposition gives the
existence of three protected states in the corresponding representation spaces
for all $\alpha\in{\mathbb{C}}^{2}$.

######  Proposition 3.2 (Existence of protected states).

For any $\alpha\in{\mathbb{C}}^{2}$ and $p,q$ as in Assumption 1 with
$p\not\equiv 0\bmod 3$, we have three protected states with their location
differing:

  1. Case I:

If $-q,p,0$ are not mutually different numbers in $\mathbb{Z}_{3}$, i.e.,
$\\{-q,p,0\\}=\\{p,0\\}$ as a set, then there exists $m\in\\{p,0\\}$ and
$m\neq n\in\\{p,0\\}$ such that

$\dim\ker_{L^{2}_{m,0}}D(\alpha)\geq 2,\ \
\ker_{L^{2}_{n,0}}D(\alpha)\neq\\{0\\};$

  2. Case II:

If $-q,p,0$ are mutually different numbers in $\mathbb{Z}_{3},$ then

$\ker_{L^{2}_{p,0}}D(\alpha)\neq\\{0\\},\ \
\ker_{L^{2}_{0,0}}D(\alpha)\neq\\{0\\},\ \
\ker_{L^{2}_{-q,0}}D(\alpha)\neq\\{0\\}.$

###### Proof.

Note that $H(\alpha)$ and $D(\alpha)$ depend analytically in
$\alpha_{12},\alpha_{23}\in{\mathbb{C}}$. In Case I, the decomposition into
representations $L^{2}_{r,\ell},\ r=p,0,-q$ and $\ell=0,1$ does not
distinguish the elements in the nullspace, i.e., there are $n\neq m,\text{
where }m,n\in{\mathbb{Z}}_{3}$ such that

$\operatorname{dim}\ker_{L^{2}_{n,0}}H(0)=1\text{ and
}\operatorname{dim}\ker_{L^{2}_{m,0}}H(0)=2.$

In view of (3.2) we only have to show that
$\dim\ker_{L^{2}_{m,0}}D(\alpha)\geq 2$.

For $\alpha\notin\mathcal{A}$, by equation (2.12), there exists
$k\in{\mathbb{C}}$ such that $(D(\alpha)-k)^{-1}:L^{2}({\mathbb{C}}/\Gamma)\to
L^{2}({\mathbb{C}}/\Gamma)$ is compact so that
$\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}D(\alpha)$ is discrete. The
rotational symmetry

$D(\alpha)\mathscr{C}=\bar{\omega}\mathscr{C}D(\alpha)$

yields that

$(D(\alpha)-k)\mathscr{C}u=\bar{\omega}\mathscr{C}(D(\alpha)-\omega k)u.$

This means that $\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}D(\alpha)$ is
invariant under multiplication by $\omega$. As
$\dim\ker_{L^{2}({\mathbb{C}}/\Gamma)}D(0)=3$ and one of the elements in
$\ker_{L^{2}({\mathbb{C}}/\Gamma)}D(\alpha)$ is protected (as it corresponds
to the protected state in $\ker_{L^{2}_{n,0}}D(\alpha)$), the
$\omega$-rotational symmetry of
$\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}D(\alpha)$ implies that
$\dim\ker_{L^{2}_{m,0}({\mathbb{C}}/\Gamma)}D(\alpha)=2$ for
$\alpha\notin\mathcal{A}$. Therefore, since $\mathcal{A}$ is discrete, we have
$\dim\ker_{L^{2}_{m,0}}D(\alpha)\geq 2$ for $\alpha\in{\mathbb{C}}^{2}$.

For Case II, the translation operator $\mathscr{L}_{a}$ splits all functions
$e_{1},e_{2},e_{3}$ in three different representation spaces
$L^{2}_{p,0},L^{2}_{0,0},L^{2}_{-q,0}$, which gives

$\dim\ker_{L^{2}_{p,0}}H(0)=\dim\ker_{L^{2}_{0,0}}H(0)=\dim\ker_{L^{2}_{-q,0}}H(0)=1.$

They correspond to three distinct points in Figure 4. Since they are simple
and the spectrum of the Hamiltonian
$\operatorname{Spec}_{L^{2}_{r,0}}H(\alpha)$ is even for $r=p,0,-q$ due to
Proposition 2.2, the analyticity of $H(\alpha)$ in
$\alpha_{12},\alpha_{23}\in{\mathbb{C}}$ yields that the eigenvalues are
protected for all $\alpha\in{\mathbb{C}}^{2}$. ∎

###  3.2. Spectral characterisation of magic parameters

Following the strategy of [BEWZ22], we now give a spectral characterisation of
the set of magic parameters using a compact operator $B_{k}$, which we will
refer as the Birman-Schwinger operator. For $\alpha=(\alpha_{12},\alpha_{23})$
with $\alpha_{12}\neq 0$, we shall for the moment consider the ratio
$\alpha_{23}/\alpha_{12}$ to be fixed but arbitrary, and view

$D(\alpha)=\begin{pmatrix}2D_{\bar{z}}&\alpha_{12}U(pz)&0\\\
\alpha_{12}U(-pz)&2D_{\bar{z}}&\alpha_{12}\tfrac{\alpha_{23}}{\alpha_{12}}U(\tfrac{p\zeta_{2}}{\zeta_{1}}z)\\\
0&\alpha_{12}\tfrac{\alpha_{23}}{\alpha_{12}}U(-\tfrac{p\zeta_{2}}{\zeta_{1}}z)&2D_{\bar{z}}\end{pmatrix}$

as a family of operators depending on a single complex parameter
$\alpha_{12}\in{\mathbb{C}}$. The non-self-adjoint operator family
${\mathbb{C}}\ni\alpha_{12}\mapsto D(\alpha):H^{1}({\mathbb{C}}/\Gamma)\to
L^{2}({\mathbb{C}}/\Gamma)$ is a holomorphic family of elliptic Fredholm
operators of index 0.

Let $H_{k}(\alpha)$ be as in (1.8), so that

$H_{k}(\alpha)=\begin{pmatrix}0&D(\alpha)^{*}+\bar{k}\\\
D(\alpha)+k&0\end{pmatrix}:H^{1}({\mathbb{C}}/\Gamma;{\mathbb{C}}^{6})\to
L^{2}({\mathbb{C}}/\Gamma;{\mathbb{C}}^{6}).$

Note that, with $e_{k}\in L^{2}({\mathbb{C}}/\Gamma;{\mathbb{C}})$ as in (2.2)
for $k\in\Gamma^{*}$, we have $2D_{\bar{z}}e_{k}=ke_{k}$, which implies that

$\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}(2D_{\bar{z}})=\Gamma^{*}.$

Therefore, we consider the identity

$D(\alpha)+k=(2D_{\bar{z}}+k)(\mathrm{Id}+\alpha_{12}T_{k})\text{ for
}k\notin\Gamma^{*},$ (3.3)

where we have introduced the operator

$T_{k}=\begin{pmatrix}0&R(k)U(pz)&0\\\
R(k)U(-pz)&0&\frac{\alpha_{23}}{\alpha_{12}}R(k)U(p\tfrac{\zeta_{2}}{\zeta_{1}}z)\\\
0&\frac{\alpha_{23}}{\alpha_{12}}R(k)U(-p\tfrac{\zeta_{2}}{\zeta_{1}}z)&0\end{pmatrix}$
(3.4)

with $R(k)=(2D_{\bar{z}}+k)^{-1}.$ We conclude that
$0\in\bigcap_{k\in{\mathbb{C}}}\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}H_{k}(\alpha)$
with $\alpha_{12}\in\mathbb{C}\setminus\\{0\\}$, $\alpha_{23}\in\mathbb{C}$ if
and only if

$-\alpha_{12}^{-1}\in\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}T_{k}$

for some $k\notin\Gamma^{*}$. Introduce the set

$\mathcal{A}_{k}=\mathcal{A}_{k}(\tfrac{\alpha_{23}}{\alpha_{12}}):=1/(\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}T_{k}\backslash\\{0\\})\text{
for }k\notin\Gamma^{*}.$

It can be shown as in [BEWZ22, Proposition 3.1] (see also [BHZ2-22,
Proposition 2.2] for a more detailed argument) that the set $\mathcal{A}_{k}$
is independent of $k\notin\Gamma^{*}$, and that for any fixed
$\tfrac{\alpha_{23}}{\alpha_{12}}$ we have

$\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}(D(\alpha))=\begin{cases}\Gamma^{*},&\alpha\notin\mathcal{A},\\\
{\mathbb{C}},&\alpha\in\mathcal{A},\end{cases}$ (3.5)

where $\mathcal{A}:=\\{(\alpha_{12},\alpha_{23})\in{\mathbb{C}}^{2}:\
\alpha_{12}\in\mathcal{A}(\tfrac{\alpha_{23}}{\alpha_{12}})\\}.$ Note that
this definition agrees with the definition of $\mathcal{A}$ in (1.11) in view
of (2.12).

######  Remark 3.

If we go back to the rescaling done for (2.1), (unrescaled) magic parameters
are defined as

$\mathcal{B}:=\\{\beta\in\mathbb{C}^{2}:\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}({\mathcal{D}}(\beta,\zeta))={\mathbb{C}}\\}.$

We note that the condition
$\operatorname{Spec}_{L^{2}({\mathbb{C}}/\Gamma)}(D(\alpha))={\mathbb{C}}$ is
invariant under rescaling of the $z$ variable. This means that the set
$\mathcal{B}$ is obtained by a rescaling of the set $\mathcal{A}$, namely

$\mathcal{B}=\frac{\zeta_{1}}{p}\mathcal{A}.$

In this article, we will for simplicity study the set of rescaled magic
parameters $\mathcal{A}$, with the exception of section §6, where we will
investigate the continuity of the set of _unrescaled_ magic parameters with
respect to the twisting angles $\zeta_{1},\zeta_{2}$.

The following proposition gives symmetries of the set $\mathcal{A}$ of magic
$\alpha$’s.

######  Proposition 3.3.

We have

$\operatorname{Spec}D(\alpha)=\operatorname{Spec}D(-\alpha)=\operatorname{Spec}D(\bar{\alpha})\quad\text{
and hence }\quad\mathcal{A}=-\mathcal{A}=\overline{\mathcal{A}}.$

###### Proof.

Conjugating $D(\alpha)$ by $\operatorname{diag}(1,-1,1)$ shows that
$\operatorname{Spec}D(\alpha)=\operatorname{Spec}D(-\alpha)$.

Applying the map $Fv(z)=\overline{v(\bar{z})}$ we find that
$(D_{\bar{z}}Fv)(z)=(D_{z}\bar{v})(\bar{z})=-(\overline{D_{\bar{z}}v})(\bar{z})=-(FD_{\bar{z}}v)(z)$
and

$D(\alpha)(Fv)=-F(D(-\bar{\alpha})v)$

which implies that
$\operatorname{Spec}(D(\alpha))=-\overline{\operatorname{Spec}(D(-\bar{\alpha}))}$.
By (3.5) we see that $\operatorname{Spec}D(\alpha)$ is invariant under
additive inversion and complex conjugation, which gives

$\operatorname{Spec}(D(\alpha))=-\overline{\operatorname{Spec}(D(-\bar{\alpha}))}=\operatorname{Spec}(D(-\bar{\alpha}))=\operatorname{Spec}(D(\bar{\alpha}))$

where the last identity follows from the first part of the proof. ∎

For fixed $\frac{\alpha_{23}}{\alpha_{12}}$, we define the _multiplicity_ of
$\alpha\in\mathcal{A}$ as

$m(\alpha):=\min_{k\notin\Gamma^{*}}\\{m(\alpha^{-1}_{12}):\alpha^{-1}_{12}\in\operatorname{Spec}_{L^{2}_{0}}T_{k}\\}.$

In particular, when $\alpha^{-1}_{12}$ is a simple eigenvalue of $T_{k}$, we
call $\alpha\in\mathcal{A}$ a _simple_ magic parameter. From the identity
(3.3), we obtain

######  Proposition 3.4.

In the notation of (1.10), the following statements are equivalent:

  1. (1)

The magic parameter $\alpha\in\mathcal{A}$ has multiplicity $m$.

  2. (2)

$\min_{k\in{\mathbb{C}}/\Gamma_{3}^{*}}\\{\dim\ker_{L_{0}^{2}}(D(\alpha)+k)\\}=m$.

  3. (3)

$E_{m}(k,\alpha)=0$ for any $k\in{\mathbb{C}}$ and $E_{m+1}(k_{0},\alpha)\neq
0$ for some $k_{0}\in{\mathbb{C}}$.

We will now establish a spectral characterization of the set $\mathcal{A}$ of
magic $\alpha$’s in terms of a scalar operator. Let $T_{k}$ be the operator in
(3.4) and note that

$\begin{split}T_{k}^{2}&=A_{k}\oplus B_{k}\text{ with
}B_{k}=(T_{k})_{21}(T_{k})_{12}+(T_{k})_{23}(T_{k})_{32}\text{ and }\\\
A_{k}&=\begin{pmatrix}(T_{k})_{12}(T_{k})_{21}&(T_{k})_{12}(T_{k})_{23}\\\
(T_{k})_{32}(T_{k})_{21}&(T_{k})_{32}(T_{k})_{23}\end{pmatrix}=\begin{pmatrix}(T_{k})_{12}\\\
(T_{k})_{32}\end{pmatrix}\begin{pmatrix}(T_{k})_{21}&(T_{k})_{23}\end{pmatrix}=:S^{(1)}_{k}S^{(2)}_{k}.\end{split}$
(3.6)

We observe that $B_{k}=S^{(2)}_{k}S^{(1)}_{k}$. Let $(\lambda,\varphi)$ be an
eigenpair of $B_{k}$. Then
$(\sqrt{\lambda},((T_{k})_{12}\varphi,\sqrt{\lambda}\varphi,(T_{k})_{32}\varphi)^{t})$
is an eigenpair of $T_{k}$. That this is no coincidence, follows from the
Schur complement formula, which shows for $\lambda\neq 0$ that

$\displaystyle\begin{pmatrix}1&S^{(1)}_{k}/\lambda\\\
S^{(2)}_{k}/\lambda&1\end{pmatrix}$ $\displaystyle=\begin{pmatrix}1&0\\\
S^{(2)}_{k}/\lambda&1\end{pmatrix}\begin{pmatrix}1&0\\\
0&1-S^{(2)}_{k}S^{(1)}_{k}/\lambda^{2}\end{pmatrix}\begin{pmatrix}1&S^{(1)}_{k}/\lambda\\\
0&1\end{pmatrix}$ $\displaystyle=\begin{pmatrix}1&S^{(1)}_{k}/\lambda\\\
0&1\end{pmatrix}\begin{pmatrix}1-S^{(1)}_{k}S^{(2)}_{k}/\lambda^{2}&0\\\
0&1\end{pmatrix}\begin{pmatrix}1&0\\\ S^{(2)}_{k}/\lambda&1\end{pmatrix}.$

This shows that
$\operatorname{Spec}(A_{k})\setminus\\{0\\}=\operatorname{Spec}(S^{(1)}_{k}S^{(2)}_{k})\setminus\\{0\\}=\operatorname{Spec}(S^{(2)}_{k}S^{(1)}_{k})\setminus\\{0\\}=\operatorname{Spec}(B_{k})\setminus\\{0\\}$
such that
$\operatorname{Spec}(T_{k})\setminus\\{0\\}=\pm\sqrt{\operatorname{Spec}(B_{k})}\setminus\\{0\\}.$

Therefore, we have a spectral characterization of the set of magic parameters
in terms of the Hilbert-Schmidt operator $B_{k}$.

######  Theorem 5.

Let $\alpha\in{\mathbb{C}}^{2}\setminus(\\{0\\}\times{\mathbb{C}})$, then with
$B_{k}$ defined in (1.15), $\alpha$ is magic if and only if
$\frac{1}{\alpha_{12}^{2}}\in\operatorname{Spec}(B_{k}(\tfrac{\alpha_{23}}{\alpha_{12}}))$
for $k\notin\Gamma^{*}.$

Finally, we observe that

$\operatorname{Spec}_{L^{2}_{r}}(T_{k})=\operatorname{Spec}_{L^{2}_{s}}(T_{k}),\quad
k\notin\Gamma^{*}.$

To see this, we note that by equation (2.6),

$\tau({s-r}):L^{2}_{r}\to
L^{2}_{s},\quad\tau(s-r)^{*}T_{k}\tau(s-r)=T_{k+s-r}.$ (3.7)

Hence, since the spectrum of $T_{k}$ is independent of $k$, it follows that if
there exists a flat band, i.e., $(D(\alpha)+k)u_{k}=0$ with $0\neq u_{k}\in
L^{2}_{r}$ for all $k\in\mathbb{C}$, then there are also $0\neq v_{k}\in
L^{2}_{s}$ such that $(D(\alpha)+k)v_{k}=0$ for all $k$.

###  3.3. Theta function argument

To construct the entire flat band from a single Bloch function with a zero,
one can use the structure of $D(\alpha)$ and employ theta functions for this
construction. To this end, we let $\theta$ be the Jacobi theta function
$\theta(\zeta):=\theta_{1}(\zeta|\omega):=-\theta_{\frac{1}{2},\frac{1}{2}}(\zeta|\omega)$
given by

$\begin{gathered}\theta(\zeta)=-\sum_{n\in\mathbb{Z}}\exp(\pi
i(n+\tfrac{1}{2})^{2}\omega+2\pi i(n+\tfrac{1}{2})(\zeta+\tfrac{1}{2})),\ \ \
\theta(-\zeta)=-\theta(\zeta),\\\ \theta(\zeta+m)=(-1)^{m}\theta(\zeta),\ \
\theta(\zeta+n\omega)=(-1)^{n}e^{-\pi in^{2}\omega-2\pi i\zeta
n}\theta(\zeta),\end{gathered}$

so that $\theta$ has simple zeros at ${\mathbb{Z}}+\omega{\mathbb{Z}}$ (and no
other zeros) – see [Mu83].

We now define

$F_{k}(z):=e^{-\frac{i}{2}(\bar{\omega}z+\bar{z})k}\frac{\theta(\frac{3z}{4\pi
i\omega}+\frac{k}{\sqrt{3}\omega})}{\theta\Big{(}\frac{3z}{4\pi
i\omega}\Big{)}}.$

A computation shows that $F_{k}(z+a)=F_{k}(z)$ for $a\in\Gamma_{3}$. We also
define

$G_{k}(z):=e^{\frac{i}{2}(\bar{k}-k\bar{\omega})z}\frac{\theta(\frac{3z}{4\pi
i\omega}+\frac{k}{\sqrt{3}\omega})}{\theta\Big{(}\frac{3z}{4\pi
i\omega}\Big{)}}$

satisfying $G_{k}(z+a)=e^{i\langle a,k\rangle}G_{k}(z)$ for $a\in\Gamma_{3}$.
Note that in particular we have $\tau(k)F_{k}(z)=G_{k}(z)$, where $\tau$ is
defined in equation (2.6). If $u\in L^{2}_{k^{\prime}}$ is a solution to
$D(\alpha)u=0$ with $u(z_{*})=0$ (cf. Lemma 8.1), then for any
$k\in{\mathbb{C}}$ we can construct

$(D(\alpha)+k)(F_{k}(z-z_{*})u(z))=0\text{ and
}D(\alpha)(G_{k}(z-z_{*})u(z))=0.$ (3.8)

In particular, multiplication by $F_{k}$ leaves $L^{2}_{k^{\prime}}$ invariant
and multiplication by $G_{k}$ maps $L^{2}_{k^{\prime}}$ to
$L^{2}_{k^{\prime}+k}.$

######  Remark 4.

We shall refer to equation (3.8) as the _theta function argument_.

###  3.4. Flat bands via Wronskians

Here we prove Theorem 1, which gives an alternative characterisation of magic
parameters $\mathcal{A}$ as the zero set of the Wronskian $W(\alpha)$ in
(1.14) of the protected states.

We devote the rest of this section to proving the theorem. First note that the
Wronskian $W(\alpha)$ is constant in $z$, as it is periodic and satisfies
equation $D_{\bar{z}}W=0$. Thus we can consider the Wronskian $W(\alpha)$ at
distinguished points which are invariant under rotation by $\omega$ modulo
$\Gamma_{3}$:

$z_{S}:=\tfrac{1}{3}(\gamma_{2}-\gamma_{1})=\tfrac{4\sqrt{3}}{9},\quad\gamma_{j}=\tfrac{4\pi
i}{3}\omega^{j}\in\Gamma_{3},\ \ j=1,2.$

In fact, $\omega(\pm z_{S})=\pm z_{S}\mp\gamma_{2}$ shows that $\\{0,\pm
z_{S}\\}$ is fixed under rotation by $\omega$ as a set of elements in
${\mathbb{C}}/\Gamma_{3}$, and we will refer to them as _stacking points_.
These are precisely the points in ${\mathbb{C}}/\Gamma_{3}$ corresponding to
the distinguished points
$\mathcal{K}=\\{-ir:r\in{\mathbb{Z}}_{3}\\}\subset{\mathbb{C}}/\Gamma_{3}^{*}$
given in (2.8). Indeed,
$\Gamma_{3}^{*}=\sqrt{3}(\omega{\mathbb{Z}}\oplus\omega^{2}{\mathbb{Z}})$, so
${\mathbb{C}}/\Gamma_{3}^{*}$ is obtained from ${\mathbb{C}}/\Gamma_{3}$
through multiplication by $3\sqrt{3}/(4\pi i)$. Hence, $\pm z_{S}$ corresponds
to

$\pm z_{S}\cdot\frac{3\sqrt{3}}{4\pi i}=\mp i\in{\mathbb{C}}/\Gamma_{3}^{*}.$

Now we consider $u\in L^{2}_{r,0}$ for some $r\in{\mathbb{Z}}_{3}$, then

$\begin{split}u(\pm z_{S})&=\mathscr{C}u(\pm z_{S})=u(\pm\omega z_{S})=u(\pm
z_{S}\mp\gamma_{2})\\\ &=\operatorname{diag}(\omega^{\pm p},1,\omega^{\mp
q})\mathscr{L}_{\mp\gamma_{2}}u(\pm z_{S})=\operatorname{diag}(\omega^{\mp
r\pm p},\omega^{\mp r},\omega^{\mp r\mp q})u(\pm z_{S}).\end{split}$ (3.9)

In Case I,

  1. (1)

if $p\equiv-q$, where we have $\varphi,\rho\in\ker_{L^{2}_{p,0}}D(\alpha)$,
$\psi\in\ker_{L^{2}_{0,0}}D(\alpha)$, using (3.9) we obtain

$\varphi_{2}(\pm z_{S})=0,\ \rho_{2}(\pm z_{S})=0,\ \psi_{1}(\pm
z_{S})=\psi_{3}(\pm z_{S})=0,$

which yields

$W(\alpha)=\psi_{2}(\varphi_{1}\rho_{3}-\varphi_{3}\rho_{1})(\pm z_{S});$

  2. (2)

if $q\equiv 0$, where we have $\varphi\in\ker_{L^{2}_{p,0}}D(\alpha)$,
$\psi,\rho\in\ker_{L^{2}_{0,0}}D(\alpha)$, we obtain

$\varphi_{2}(\pm z_{S})=\varphi_{3}(\pm z_{S})=0,\ \psi_{1}(\pm z_{S})=0,\
\rho_{1}(\pm z_{S})=0,$

which yields

$W(\alpha)=\varphi_{1}(\psi_{2}\rho_{3}-\psi_{3}\rho_{2})(\pm z_{S}).$

In Case II, we obtain

$\varphi_{2}(\pm z_{S})=\varphi_{3}(\pm z_{S})=0,\ \ \psi_{1}(\pm
z_{S})=\psi_{3}(\pm z_{S})=0,\ \ \rho_{1}(\pm z_{S})=\rho_{2}(\pm z_{S})=0,$

therefore

$W(\alpha)=\varphi_{1}(\pm z_{S})\psi_{2}(\pm z_{S})\rho_{3}(\pm z_{S}).$

This leads to the following result.

######  Proposition 3.5.

In the notation of equation (1.13) and (1.14), we have

$W(\alpha)=0\quad\Longrightarrow\quad\alpha\in\mathcal{A}.$

###### Proof.

For Case I, we only discuss the subcase $p\equiv-q$ as the other subcase is
similar. Assume $W(\alpha)=0$. Then if $\psi(\pm z_{S})=0$, we can apply the
theta function argument (cf. equation (3.8)) to $\psi$ to conclude that
$\alpha\in\mathcal{A}$. Otherwise if $\varphi_{1}\rho_{3}-\varphi_{3}\rho_{1}$
vanishes at $\pm z_{S}$, we have that
$c_{1}\varphi+c_{2}\rho\in\ker_{L^{2}_{p,0}}D(\alpha)$ vanishes at $\pm z_{S}$
for some $c_{1},c_{2}\in{\mathbb{C}}$ not both equal to zero. Applying the
theta function argument (3.8) to $c_{1}\varphi+c_{2}\rho$ we find that
$\alpha\in\mathcal{A}$.

For Case II, we have

$W(\pm z_{S})=\varphi_{1}(\pm z_{S})\psi_{2}(\pm z_{S})\rho_{3}(\pm z_{S}).$

Therefore, $W=0$ implies that one of $\varphi,\psi,\rho$ vanishes at $z_{S}$
and one of $\varphi,\psi,\rho$ vanishes at $-z_{S}$. In either case, the theta
function argument (3.8) gives $\alpha\in\mathcal{A}$. ∎

The kernel of $D(\alpha)$ provides an expression for the inverse of
$D(\alpha)-k$ with $k\notin\Gamma^{*}$ and $\alpha\notin\mathcal{A}$. In
particular, the following result shows that $\alpha\in\mathcal{A}$ implies
that $W(\alpha)=0$. Combined with Proposition 3.5, this proves Theorem 1.

######  Proposition 3.6.

In view of the notation in (1.13) and (1.14), we define a three-by-three
matrix

$\mathbf{W}(\alpha):=[\varphi,\psi,\rho],\ \
W(\alpha)=\det\mathbf{W}(\alpha).$

Then $W(\alpha)\neq 0$ and $k\notin\Gamma^{*}$ imply that, with the cofactor
matrix denoted by $\operatorname{adj}$,

$(D(\alpha)-k)^{-1}=\frac{1}{W(\alpha)}\mathbf{W}(\alpha)(2D_{\bar{z}}-k)^{-1}({\rm{adj}}\mathbf{W}(\alpha)).$
(3.10)

For a fixed $k\notin\Gamma^{*}$, $\alpha_{12}\mapsto(D(\alpha)-k)^{-1}$ is a
meromorphic family of compact operators with poles of finite rank at
$\alpha_{12}\in\mathcal{A}(\frac{\alpha_{23}}{\alpha_{12}})$.

###### Proof.

If $W(\alpha)\neq 0$, then
$\mathbf{W}(\alpha)^{-1}=(W(\alpha))^{-1}{\rm{adj}}\mathbf{W}(\alpha)$. Write
$D(\alpha)=2D_{\bar{z}}+\mathbf{U}_{\alpha}(z)$, and note that

$2D_{\bar{z}}[\mathbf{W}(\alpha)]=-\mathbf{U}_{\alpha}(z)\mathbf{W}(\alpha)\quad\Longleftrightarrow\quad\mathbf{W}(\alpha)=e^{-\frac{i}{2}\int\mathbf{U}_{\alpha}(z)d\bar{z}}.$

Hence, $\mathbf{W}(\alpha)^{-1}$ is the integrating factor of the equation

$2D_{\bar{z}}u+(\mathbf{U}_{\alpha}(z)-k)u=v$

so that

$(2D_{\bar{z}}-k)(\mathbf{W}(\alpha)^{-1}u)=\mathbf{W}(\alpha)^{-1}v.$

This yields equation (3.10). In view of equation (3.3),

$(D(\alpha)-k)^{-1}=(I+\alpha_{12}T_{k})^{-1}(D(0)-k)^{-1},$

where, using analytic Fredholm theory (see for instance [DyZw19, Theorem
C.8]), we see that $\alpha_{12}\mapsto(I+\alpha_{12}T_{k})^{-1}$ is a
meromorphic family of operators with poles of finite rank. ∎

##  4\. Trace formulæ

In this section, we use the Birman-Schwinger operator $B_{k}$, defined in
equation (1.15), to analyse the set of magic parameters
$\mathcal{A}(\alpha_{23}/\alpha_{12})$ for a fixed hopping ratio. More
precisely, using Theorem 2, we see that the traces of $B_{k}$ are actually
equal to the sums of powers of magic parameters. Following [BHZ1-22, Theorem
4], our goal is to provide a semi-explicit formula for

$\operatorname{tr}(B_{k}^{\ell})=\sum_{\alpha_{12}\in\mathcal{A}(\frac{\alpha_{23}}{\alpha_{12}})}\frac{1}{\alpha_{12}^{2\ell}},\quad\ell\geq
2.$

We will prove the following result, which is a slight variation of [BHZ1-22,
Theorem 4].

######  Theorem 6.

Let $B_{k}:L^{2}_{0}\to L^{2}_{0}$ be a meromorphic family of Hilbert-Schmidt
operators defined for $k\notin\Gamma^{*}$. We suppose that $B_{k}$ satisfies
the three properties stated in Lemma 4.1. Then one has, for any $\ell\geq 2$,

$\begin{split}\operatorname{tr}(B_{k}^{\ell})&=\frac{2i\pi\omega}{3\sqrt{3}}\sum_{n\in\mathbb{Z}}n\left[\sum_{m\in\mathbb{Z}}\mathrm{Res}\big{(}\langle
B_{k}^{\ell}e_{i},e_{i}\rangle_{L^{2}},\sqrt{3}(m\omega^{2}-n\omega)\big{)}\right],\end{split}$
(4.1)

where all sums are finite sums.

###### Proof.

The only difference compared to the study of twisted bilayer graphene
[BHZ1-22, Theo.4] is that the set of forbidden values for the parameter $k$ is
$\Gamma^{*}$ instead of $(3\Gamma^{*}+i)\cup(3\Gamma^{*}-i)$. This causes no
harm to the proof as we can, by the residue theorem, translate the contour
integral, used in the proof, by a small parameter $\varepsilon>0$ without
changing the final value. ∎

######  Remark 5.

Using the conjugation (3.7), we see that the trace, for the operator $B_{k}$
defined in (1.15), is the same on any space $L^{2}_{ir},$ where
$r\in{\mathbb{Z}}_{3}$. The trace on the whole $L^{2}$ space is then $9$ times
the result of (4.1), see equation (2.11) and the discussion thereafter. Most
of the results derived from these trace formulas will only use algebraic
properties of the traces and will thus not be affected by the possible
rescaling of (4.1).

The proof now boils down to proving the following lemma, which is a slight
adaptation of [BHZ1-22, Lemma 2.2].

######  Lemma 4.1.

Let $e_{k}(z)$ be as in (2.2) for $k\in\Gamma^{*}$. Consider a potential $U\in
C^{\infty}(\mathbb{C}/\Gamma;\mathbb{C})$ satisfying the first two symmetries
of (1.4) with a finite number of non-zero Fourier mode in its decomposition
(2.4). For $\ell\geq 2$, one has:

  * •

The trace is constant in $k$,

$\operatorname{tr}(B_{k}^{\ell})=\tau_{\ell}\text{ independent of
}k\in\mathbb{C}\setminus\Gamma^{*}.$

  * •

The function $\mathbb{C}\setminus\Gamma^{*}\ni k\mapsto\langle
B_{k}^{\ell}e_{m},e_{m}\rangle_{L^{2}}$ is a finite sum of rational fractions
on the complex plane $\mathbb{C}$ with degree equal to $-2\ell$ and with $($a
finite number of$)$ poles contained in $\Gamma^{*}$.

  * •

For any $\gamma\in\Gamma^{*}$ and for any $k\notin\Gamma^{*}$, we have

$\langle B_{k}^{\ell}e_{3\gamma+i},e_{3\gamma+i}\rangle_{L^{2}}=\langle
B_{k-3\gamma}^{\ell}e_{i},e_{i}\rangle_{L^{2}}.$

###### Proof.

The first point is a consequence of the fact that the spectrum of $B_{k}$
doesn’t depend on $k$, which follows from Theorem 2. For the last two points,
we prove by induction that $k\mapsto B_{k}^{\ell}e_{3\gamma+i}$, where
$\gamma\in\Gamma^{*}$, is of the form

$B_{k}^{\ell}e_{3\gamma+i}=\sum_{\nu\in F}R_{\nu+3\gamma}(k)e_{\nu+3\gamma},$
(4.2)

where $F\subset\Gamma^{*}$ is a finite set and $R_{\nu}(k)$ is a sum of
rational fractions of degree $-2\ell$ with poles located on $\Gamma^{*}$.
Moreover, we will prove that one has the relation
$R_{\nu+3\gamma}(k)=R_{\nu}(k-3\gamma).$

The result is clear for $\ell=0$. Suppose the result is true for $\ell$, and
let’s prove it holds for $\ell+1$. The main observation is that multiplication
by $U(\pm z)$ acts as a shift on the Fourier basis. The multiplication by
$U(-z)$ sends $e_{\nu}$ to a linear combination of $e_{\ell}$ for
$\ell\in\Gamma^{*}$. Then applying $(D(0)-k)^{-1}$ multiplies the coefficient
of $e_{\ell}$ by $(\ell-k)^{-1}$. Multiplying by $U(z)$ gives back a linear
combination of $e_{\nu}$ with $\nu\in\Gamma^{*}$. Finally, applying
$(D(0)-k)^{-1}$ multiplies the coefficient of $e_{\nu}$ by $(\nu-k)^{-1}$.
This means that, using the induction hypothesis (4.2),

$B_{k}^{\ell+1}e_{3\gamma+i}=\sum_{\nu\in F}R_{\nu+3\gamma}(k)\sum_{\eta\in
L}\sum_{\beta\in
L}\frac{a_{\eta}}{k-(\nu+\beta-\eta+3\gamma)}\frac{a_{\beta}}{k-(\nu+\beta+3\gamma)}e_{\nu+\beta-\eta+3\gamma},$

where $L\subset\Gamma^{*}$ is a finite subset that depends only on $U$ and
$a_{\bullet}$ are constants. Thus, it is clear from this formula that the
induction carries on to $\ell+1$. This concludes the proof of the lemma. ∎

The trace formula is not fully explicit but it is sufficient to prove that the
traces, for the standard potential, are of the form $q_{\ell}\pi/\sqrt{3}$,
for $q_{l}$ a rational number. This extraordinary rational condition seems to
reflect some _hidden integrability_ of the set of magic parameters, which we
do not fully understand yet. However, just like in [BHZ1-22], this condition
implies that the set of magic parameters is infinite. Now that we have the
trace formula of Theorem 6, the proofs of [BHZ1-22, Theorems 5 and 6] carry
over to our setting and we get

######  Theorem 7.

Consider a potential $U\in C^{\infty}(\mathbb{C}/\Gamma;\mathbb{C})$
satisfying the first two symmetries of (1.4) with finitely many non-zero
Fourier modes $c_{n}\in\mathbb{Q}(\omega/\sqrt{3})$ appearing in the
decomposition (2.4). Suppose moreover that
$\alpha_{23}/\alpha_{12}\in\mathbb{Q}$. Then for any $\ell\geq 2$, one has
$\tau_{\ell}\in\pi\mathbb{Q}(\omega/\sqrt{3})$. If $U$ also has the third
symmetry of (1.4) then the traces are real and thus
$\tau_{\ell}\in\pi\mathbb{Q}/\sqrt{3}$. In particular, for all potentials
satisfying all three symmetries in (1.4), one has

$\forall\ell\geq
2,\quad\operatorname{tr}(B_{k}^{\ell})=\sum_{\alpha_{12}\in\mathcal{A}(U,\alpha_{23}/\alpha_{12})}\alpha_{12}^{-2\ell}=\frac{\pi}{\sqrt{3}}q_{\ell},\;\;\;q_{\ell}\in\mathbb{Q},$

where $\mathcal{A}(U,\alpha_{23}/\alpha_{12})$ is the set of magic parameters
counting multiplicity for a potential $U$ and fixed hopping ratio
$\alpha_{23}/\alpha_{12}$.

We deduce also the infiniteness of magic parameters, for our canonical choice
of potential.

######  Theorem 8.

Under the assumptions and with the same notation as in Theorem 7 one has the
implication

$|\mathcal{A}(U,\alpha_{23}/\alpha_{12})|>0\implies|\mathcal{A}(U,\alpha_{23}/\alpha_{12})|=+\infty.$

Let $N\geq 0$, for a tuple $c=(c_{n})_{\\{n:\|n\|_{\infty}\leq N\\}}$, define
$U$ to be the potential defined by (2.4). Then the above implication holds for
a generic $($in the sense of Baire$)$ set of coefficients
$c=(c_{n})_{\\{n:\|n\|_{\infty}\leq N\\}}\in\mathbb{C}^{(2N+1)^{2}}$ that
contains $\mathbb{Q}(\omega/\sqrt{3})^{(2N+1)^{2}}.$

Finally, we extend the result of [BHZ23] on the existence of non-simple magic
parameters, see Figure 5.

Figure 5. Fixing $\alpha_{32}=\alpha_{12}$ and $\zeta=(1,3)$ (left) and
$\zeta=(2,3)$ (right) we plot magic parameters $\alpha_{12}$. Non-simple
eigenvalues are indicated by numbers. The multiplicity of the flat bands of
the Hamiltonian is twice the multiplicity indicated in the figures.

For the rest of this subsection, we will assume we are in Case I, i.e.,
$p,-q,0$ are not mutually different modulo $3$. The argument in [BHZ23] was
based on the possibility of taking $k=0$ on a specific translational invariant
subspace $L^{2}_{-ir}$ for some $r\in{\mathbb{Z}}_{3}$. Then, using a theta
function argument, it was shown (in [BHZ23, Theorem 1]) that the eigenfunction
of a simple magic parameter $\alpha$ on this subspace $L^{2}_{-ir}$ has to
live in the smaller rotational invariant subspace $L^{2}_{r,2}$ (cf. (2.9)).
The existence of non-simple magic parameters was then shown by proving that
the traces of the Birman-Schwinger operator $B_{0}$ at $k=0$ on $L^{2}_{r,1}$
did not all vanish. In the present case, we recall that the Birman-Schwinger
operator, acting on $L^{2}_{0}$ of functions invariant by the translations
$\mathscr{L}_{a}$222We define the action of $\mathscr{L}_{a}$ on scalar valued
functions by projecting on the third component, see (3.6)., is given by

$B_{k}(\tfrac{\alpha_{23}}{\alpha_{12}}):=R(k)_{p}^{-1}U(-pz)R(k)_{0}U(pz)+(\tfrac{\alpha_{23}}{\alpha_{12}})^{2}R(k)_{p}U(p\tfrac{\zeta_{2}}{\zeta_{1}}z)R(k)_{p+p\frac{\zeta_{2}}{\zeta_{1}}}^{-1}U(-p\tfrac{\zeta_{2}}{\zeta_{1}}z).$
(4.3)

Here, we have denoted by $R(k)_{p}$ the restriction of the resolvent
$R(k)=(2D_{\bar{z}}+k)^{-1}$ on the translational invariant subspace
$L^{2}_{p}$. The presence of the factor $R(k)_{0}$ makes it impossible to put
$k=0$. Note that $p\frac{\zeta_{2}}{\zeta_{1}}\equiv q\bmod 3$, which means
that $p,0,p+p\frac{\zeta_{2}}{\zeta_{1}}$ are not mutually different modulo
$3$. Recall that we defined

$L^{2}_{ir}:=\\{u\in
L^{2}(\mathbb{C}/\Gamma;\mathbb{C}):\mathscr{L}_{a}u=e^{i\langle
ir,a\rangle}u=\bar{\omega}^{r}u,\ a\in\Gamma/3\\},\quad r\in\mathbb{Z}_{3}.$

The restrictions of $B_{k}$ on $L^{2}_{ip}$ are conjugated by the operator
$\tau(\cdot)$ defined in (2.6). We see that, for $r\in\mathbb{Z}_{3}$,

$(B_{k})_{L^{2}_{ir}}:=R(k)_{p-r}^{-1}U(-pz)R(k)_{-r}U(pz)+(\tfrac{\alpha_{23}}{\alpha_{12}})^{2}R(k)_{p-r}U(p\tfrac{\zeta_{2}}{\zeta_{1}}z)R(k)_{p+p\frac{\zeta_{2}}{\zeta_{1}}-r}^{-1}U(-p\tfrac{\zeta_{2}}{\zeta_{1}}z).$

######  Remark 6.

In particular, there is a space $\mathcal{H}=L^{2}_{ir}$ such that one can
take $k=0$ for $(B_{k})_{\mathcal{H}}$. Because the spectrum of $T_{k}$ is the
same on each of these spaces, we can suppose, without loss of generality that
$p=-q\equiv-1$ modulo $3$. This means we can choose $r=1$ (see Remark 2), i.e
$\mathcal{H}=L^{2}_{i}$, and we study $(B_{0})_{L^{2}_{i}}$, and we use
Proposition 7.1 to get

$\alpha\in\mathcal{A}\text{ is simple}\implies
u_{i}\in\ker_{L^{2}_{i}}D(\alpha)\implies
u_{i}\in\ker_{L^{2}_{-1,-1}}D(\alpha).$

Just like in the proof of [BHZ23, Theorem 1], the existence of a non-simple
magic parameter is obtained by showing that
$\operatorname{tr}\big{(}(B_{0})^{\ell}_{L^{2}_{-1,1}}\big{)}$ is non-zero for
a well chosen value of $\ell$. The proof is exactly the same as in the case of
twisted bilayer graphene, where the argument was based on the fact that one
can write

$\operatorname{tr}\big{(}(B_{0})^{\ell}_{L^{2}_{-1,1}}\big{)}=\frac{1}{3}\operatorname{tr}\big{(}(B_{0})^{\ell}_{L^{2}_{i}}\big{)}+\mathcal{R}_{-1,1}^{\ell},$

where
$\operatorname{tr}\big{(}(B_{0})^{\ell}_{L^{2}_{i}}\big{)}=q_{\ell}\pi/\sqrt{3}$
for a non-zero rational number $q_{\ell}$, which means that the first term is
transcendental. Now, we were able to show in [BHZ23] that the remainder
$\mathcal{R}_{-1,1}^{\ell}$ was algebraic, thus proving that the trace is non-
zero. We get the following refinement of Theorem 8.

######  Theorem 9.

Suppose that we are in Case I. Then, under the assumptions and with the same
notation as in Theorem 7, one has the implication

$|\mathcal{A}(U,\alpha_{23}/\alpha_{12})|>0\implies|\mathcal{A}_{m}(U,\alpha_{23}/\alpha_{12})|=+\infty.$

Here, we denoted by $\mathcal{A}_{m}(U,\alpha_{23}/\alpha_{12})$ the set of
non-simple magic parameters and by $\mathcal{A}(U,\alpha_{23}/\alpha_{12})$
the set of all magic parameters for the potential $U$. In particular, the set
of magic parameters for our canonical potential $U=U_{0}$ defined in (1.6) is
infinite.

Let $N\geq 0$, and for a tuple $a=(a_{n})_{\\{n:\|n\|_{\infty}\leq N\\}}$,
define $U_{a}$ to be the potential defined by (2.4) with coefficients $a$.
Then the above implication holds for a generic $($in the sense of Baire$)$ set
of coefficients $a=(a_{n})_{\\{n:\|n\|_{\infty}\leq
N\\}}\in\mathbb{C}^{(2N+1)^{2}}$ that contains
$(\mathbb{Q}(\omega/\sqrt{3}))^{(2N+1)^{2}}.$

######  Remark 7.

In Case II, the above argument fails as we cannot consider $k=0$ on any of the
rotational invariant subspaces. We note that numerical evidence suggest that
there are choices of $\zeta_{1},\zeta_{2}$ for which all magic parameters are
simple, see Figure 6.

Figure 6.  For $\zeta=(1,4)$, $\alpha_{12}=\alpha_{23}$, all magic parameters
seem to be simple. For the largest (with positive real part) magic parameter a
contour plot of the element of $\ker(D(\alpha))$ is shown.

##  5\. Traces for powers of order 2

Note that none of the results above proves the existence of a magic parameter.
For the bilayer case, the existence, for the Bistritzer–MacDonald potential,
was proven by explicitly computing the trace for $\ell=2$ and showing that it
is non-zero. In the following section, we provide for $U=U_{0}$ explicit
formulas for $\operatorname{tr}(B_{k}^{2})$ that depend on the choice of
parameters $\zeta_{1},\zeta_{2}$, and $\alpha_{23}/\alpha_{12}$.

For our numerics, it is convenient to use rectangular coordinates $z=2i(\omega
y_{1}+\omega^{2}y_{2})$, see [BEWZ22, §3.3] for details. In these coordinates,
we may introduce

$\begin{gathered}{\mathcal{D}_{k}}:={\omega^{2}(D_{y_{1}}+k_{1})-\omega(D_{y_{2}}+k_{2})},\\\
{\mathcal{V}(y)}:=\sqrt{3}(e^{-i(y_{1}+y_{2})}+\omega
e^{i(2y_{1}-y_{2})}+\omega^{2}e^{i(-y_{1}+2y_{2})}),\end{gathered}$

with periodic periodic boundary conditions (for $y\mapsto y+2\pi n$,
$n\in{\mathbb{Z}}^{2}$). In the following, we shall write
$\widehat{\mathcal{V}_{\pm}}(y):=\widehat{\mathcal{V}(\pm y)}.$ The operator
$\ B_{k}$, defined in (4.3), reads in the new coordinates

$\begin{gathered}\hat{B}_{k}(\alpha_{23}/\alpha_{12}):=\widehat{\mathcal{D}}_{k}^{-1}\widehat{\mathcal{V}}_{-}(p)\widehat{\mathcal{D}}_{k}^{-1}\widehat{\mathcal{V}}_{+}(p)+(\alpha_{23}/\alpha_{12})^{2}\left(\widehat{\mathcal{D}}_{k}^{-1}\widehat{\mathcal{V}}_{+}(\tilde{p})\widehat{\mathcal{D}}_{k}^{-1}\widehat{\mathcal{V}}_{-}(\tilde{p})\right),\end{gathered}$

where we denoted by $\tilde{p}=p\zeta_{2}/\zeta_{1}$ and we introduced

$\begin{gathered}\widehat{\mathcal{D}}_{k}:={\omega^{2}(D+k_{1})-\omega(D+k_{2})},\text{
with }D=\operatorname{diag}(\ell)_{\ell\in{\mathbb{Z}}}\\\
\widehat{\mathcal{V}}_{\pm}(y):=\sqrt{3}\left(J^{\pm}\otimes J^{\pm}+\omega
J^{\mp 2}\otimes J^{\pm}+\omega^{2}J^{\pm}\otimes J^{\mp
2}\right),\end{gathered}$

where $J$ is the right-shift $J((a_{n})_{n})=(a_{n+1})_{n}$ – see [BEWZ22,
(3.17)]. The space $L^{2}_{i}(\mathbb{C}/\Gamma;\mathbb{C})$ corresponds to

$\begin{gathered}\ell_{(1,1)}^{2}:=\\{f\in\ell^{2}(\mathbb{Z}^{2}):\forall\,n\notin\left(3\mathbb{Z}+1\right)\times\left(3\mathbb{Z}+1\right),\>\>f_{n}=0\\}.\end{gathered}$

As in [BEWZ22, §3.3], we introduce auxiliary operators $J^{p,q}:=J^{p}\otimes
J^{q},\>\>p,q\in\mathbb{Z}$. For a diagonal matrix
$\Lambda=(\Lambda_{i,j})_{i,j\in\mathbb{Z}}$ acting on
$\ell^{2}(\mathbb{Z}^{2})$, we define a new diagonal matrix

$\Lambda_{p,q}:=(\Lambda_{i+p,j+q})_{i,j\in\mathbb{Z}}.$

We recall the following properties [BEWZ22, (3.24)]

$J^{p,q}\Lambda
J^{p^{\prime},q^{\prime}}=\Lambda_{p,q}J^{p+p^{\prime},q+q^{\prime}}=J^{p+p^{\prime},q+q^{\prime}}\Lambda_{-p^{\prime},-q^{\prime}}.$

Denoting the inverse of $\widehat{\mathcal{D}}_{k}^{-1}$ by

$\Lambda=\Lambda_{k}:=\widehat{\mathcal{D}}_{k}^{-1},\ \ \
\Lambda_{m,n}=\frac{1}{\omega^{2}(m+k_{1})-\omega(n+k_{2})},\ \ \
(k_{1},k_{2})\notin\mathbb{Z}^{2},$

we see that

$\displaystyle\frac{1}{3}\widehat{B}_{k}$
$\displaystyle=\Lambda\left(J^{-p,-p}+\omega
J^{2p,-p}+\omega^{2}J^{-p,2p}\right)\Lambda\left(J^{p,p}+\omega
J^{-2p,p}+\omega^{2}J^{p,-2p}\right)$
$\displaystyle\quad+(\alpha_{23}/\alpha_{12})^{2}\Lambda\left(J^{\tilde{p},\tilde{p}}+\omega
J^{-2\tilde{p},\tilde{p}}+\omega^{2}J^{\tilde{p},-2\tilde{p}}\right)\Lambda\left(J^{-\tilde{p},-\tilde{p}}+\omega
J^{2\tilde{p},-\tilde{p}}+\omega^{2}J^{-\tilde{p},2\tilde{p}}\right).$

This means that when expanding $\widehat{B}_{k}^{\ell},$ we get a diagonal
term of the following form:

$\begin{gathered}((\widehat{B}_{k}^{\ell})_{ii})_{i\in{\mathbb{Z}}}=\sum_{\pi\in\Theta_{\ell}}a_{\pi}\prod_{i=1}^{\ell}\Lambda_{\tilde{\alpha}_{i},\tilde{\beta}_{i}}\Lambda_{{\tilde{\gamma}_{i}},{\tilde{\delta}_{i}}},\\\
\pi:=\left[\big{(}(\alpha_{1},\beta_{1}),n_{1}\big{)},\big{(}(\gamma_{1},\delta_{1}),n_{1}\big{)},\big{(}(\alpha_{2},\beta_{2}),n_{2}\big{)},...,\big{(}(\gamma_{\ell},\delta_{\ell}),n_{l}\big{)}\right],\end{gathered}$

where

$\begin{gathered}\tilde{\alpha}_{i}=\sum_{j=1}^{i-1}n_{j}(\alpha_{j}+\gamma_{j})\quad\tilde{\beta}_{i}=\sum_{j=1}^{i-1}n_{j}(\beta_{j}+\delta_{j}),\quad{\tilde{\gamma}_{i}}=n_{i}\alpha_{i}+\sum_{j=1}^{i-1}n_{j}(\alpha_{j}+\gamma_{j}),\\\
\tilde{\delta}_{i}=n_{i}\beta_{i}+\sum_{j=1}^{i-1}n_{j}(\beta_{j}+\delta_{j}).\end{gathered}$

Here, the tuple $\pi$ is such that for each $1\leq k\leq\ell,$ we have either

$n_{k}=p,\,\,\text{and
}\,-(\alpha_{k},\beta_{k}),(\gamma_{k},\delta_{k})\in\\{(1,1),(-2,1),(1,-2)\\},$

or we have

$n_{k}=\tilde{p},\,\,\text{and
}\,(\alpha_{k},\beta_{k}),-(\gamma_{k},\delta_{k})\in\\{(1,1),(-2,1),(1,-2)\\}.$

We call such a tuple an admissible tuple. The sum is taken over the
combinatorial set

$\Theta_{\ell}:=\\{\pi:\pi\text{ is admissible and
}\widetilde{\alpha_{\ell+1}}=\widetilde{\beta_{\ell+1}}=0\\}.$ (5.1)

In other words, we consider all tuples of the form

$\pi:=\left[\big{(}(\alpha_{1},\beta_{1}),n_{1}\big{)},\big{(}(\gamma_{1},\delta_{1}),n_{1}\big{)},\big{(}(\alpha_{2},\beta_{2}),n_{2}\big{)},...,\big{(}(\gamma_{\ell},\delta_{\ell}),n_{l}\big{)}\right],$

admissible in the previous sense and such that the total sum of the (weighted)
couples is zero. Finally, the coefficient $a_{\pi}$ is given by

$a_{\pi}=3^{\ell}(\alpha_{23}/\alpha_{12})^{2s_{\pi}}\omega^{m_{\pi}},\quad
m_{\pi}:=\tfrac{2}{3}\sum_{i=1}^{\ell}(\gamma_{i}+\beta_{i}),$

where $s_{\pi}$ denotes the number of $k$’s such that $n_{k}=\tilde{p}$. We
get the equivalent trace formula given in Theorem 6, in these Fourier
coordinates. The translation between the two coordinate systems is exactly the
same as in [BHZ1-22, Theorem 7]. To state the result we introduce the notation
$\gamma_{(a,b)}=\omega^{2}a-\omega b$, together with the orbits of

$\sigma:\mathbb{Z}\times\mathbb{Z}\to\mathbb{Z},\>\>\>\sigma(m,n)=(-(n+m),m),$
(5.2)

which satisfies

$\gamma_{\sigma(m,n)}=\omega\gamma_{(m,n)}.$

This implies that $\sigma^{3}=\operatorname{id}$ and that the orbits under
$\sigma$ are of cardinality $3$. For a tuple
$\pi=\left[((\alpha_{1},\beta_{1}),n_{1}),((\gamma_{1},\delta_{1}),m_{1}),((\alpha_{2},\beta_{2}),n_{2}),((\gamma_{2},\delta_{2}),m_{2})\right]\in\Theta_{2}$,
we define

$\sigma(\pi)=\left[(\sigma(\alpha_{1},\beta_{1}),n_{1}),(\sigma(\gamma_{1},\delta_{1}),m_{1}),(\sigma(\alpha_{2},\beta_{2}),n_{2}),(\sigma(\gamma_{2},\delta_{2}),m_{2})\right].$

######  Theorem 10.

Let ${\ell}\geq 2$ and $\Theta_{\ell}$ be as in (5.1) with coefficients
$\tilde{\alpha}_{i},..,\tilde{\delta}_{i},m_{\pi}$ as defined above. Then the
traces are given by

${\operatorname{tr}}\big{(}\widehat{B}_{k}^{\ell}\big{)}=\frac{-2i\omega\pi}{9}\sum_{\pi\in\Theta_{\ell}}(\alpha_{23}/\alpha_{12})^{2s_{\pi}}\mathcal{C}(\pi),$

where we define the _contribution_ of the tuple $\pi$ as

$\mathcal{C}(\pi)=\sum_{(\eta_{i},\varepsilon_{i})\in
P(f_{\sigma(\pi)})}{\operatorname{Res}}(f_{\sigma(\pi)},-\gamma_{(\eta_{i},\varepsilon_{i})})\varepsilon_{i}.$

Here, the summation is taken on
$P(f_{\pi}):=\\{-{({\tilde{\alpha}_{i}},{\tilde{\beta}_{i}})},-{({\tilde{\gamma}_{i}},{\tilde{\delta}_{i}})},1\leqslant
i\leqslant 2\\},$ the set of indices of poles of $f_{\pi}$, with

$f_{\pi}(k):=3^{\ell}\omega^{m_{\pi}}\prod_{i=1}^{\ell}\frac{1}{(k+\gamma_{({\tilde{\alpha}_{i}},{\tilde{\beta}_{i}})}+\mu)(k+\gamma_{({\tilde{\gamma}_{i}},{\tilde{\delta}_{i}})}+\mu)},\quad\mu:=\omega^{2}-\omega.$

Thus, there are two components in the trace formula. The first one is the set
$\Theta_{\ell}$. Note that this set only depends on $\zeta_{2}/\zeta_{1}$ but
_not_ on $\alpha_{23}/\alpha_{12}$. The second component is given by the
coefficients $a_{\pi}$ which depend on both parameters, with a polynomial
dependence in $\alpha_{23}/\alpha_{12}$. The full computation of
$\Theta_{\ell}$ is the main difficulty to make the trace formula from Theorem
6 fully explicit. This combinatorial set describes how the Jordan blocks in
$\hat{B}_{k}$ interact with each other when computing the powers of
$\hat{B}_{k}$ and thus reflects the _non-normal_ nature of the operator
$B_{k}$. However, for $\ell=2$ and for $U=U_{0}$, the set $\Theta_{\ell}$ is
sufficiently small to be computable by hand explicitly. We get different
formulas for different values of $\zeta_{2}/\zeta_{1}$. We sum up the
computations in the following theorem. When studying the continuity of the set
of magic parameters with respect to parameters $\zeta_{1},\zeta_{2}$ (see
Section 6), we will study the set of _unrescaled magic parameters_. Using the
rescaling of equation (1.1), we see that this set is defined to be

$\mathcal{B}=\frac{\zeta_{1}}{p}\mathcal{A},$

see Remark 3. Using the same notation $\mathcal{B}(r)$, we deduce the
following explicit formulas for

$\mathcal{S}_{2\ell}(\beta_{23}/\beta_{12})=\sum_{\beta_{12}\in\mathcal{B}(\beta_{23}/\beta_{12})}\frac{1}{\beta_{12}^{2\ell}},\quad\ell\geq
2,$

for $\ell=2$.

######  Theorem 11.

Consider $U=U_{0}$ defined in (1.6). Denote by
$\mathcal{B}(\beta_{23}/\beta_{12},\zeta_{1},\zeta_{2})$ the set of unrescaled
magic parameters for the choice of parameters $\beta_{23}/\beta_{12}$,
$\zeta_{1}$, $\zeta_{2}$, and let
$\mathcal{S}_{2\ell}(\beta_{23}/\beta_{12},\zeta_{1},\zeta_{2})$ denote the
corresponding trace.

  * •

If $\zeta_{2}/\zeta_{1}=-1,$ then we have

$\mathcal{S}_{4}(\beta_{23}/\beta_{12},\zeta_{1},-\zeta_{1})=\left(1+\left(\frac{\beta_{23}}{\beta_{12}}\right)^{2}\right)^{2}\frac{4\pi
p^{2}}{9\sqrt{3}\zeta_{1}^{4}}.$

  * •

If $\zeta_{2}/\zeta_{1}\neq\pm 1,$ then we have

$\mathcal{S}_{4}(\beta_{23}/\beta_{12},\zeta_{1},\zeta_{2})=\frac{4\pi
p^{2}}{9\sqrt{3}\zeta_{1}^{4}}\left(\frac{1}{(\zeta_{2}/\zeta_{1})^{2}}\left(\frac{\beta_{23}}{\beta_{12}}\right)^{4}+\frac{3}{1-\zeta_{2}/\zeta_{1}+(\zeta_{2}/\zeta_{1})^{2}}\left(\frac{\beta_{23}}{\beta_{12}}\right)^{2}+1\right).$

  * •

If $\zeta_{2}/\zeta_{1}=1,$ then we have

$\mathcal{S}_{4}(\beta_{23}/\beta_{12},\zeta_{1},\zeta_{1})=\left(1-\left(\frac{\beta_{23}}{\beta_{12}}\right)^{2}+\left(\frac{\beta_{23}}{\beta_{12}}\right)^{4}\right)\frac{4\pi
p^{2}}{9\sqrt{3}\zeta_{1}^{4}}.$

In particular, if the hopping ratio $\beta_{23}/\beta_{12}\in\mathbb{R}$, then
$\mathcal{S}_{4}>0$ and there exist magic parameters.

###### Proof.

We will prove the equivalent formulas for $\mathcal{A}$ to keep the notation
of the article consistent. First note that the case

$\zeta_{2}/\zeta_{1}=-1,$

reduces to the bilayer case and the formula thus follows from [BEWZ22]. We
list a few tricks that make the calculation easier.

  1. (1)

First, it is easy to see that $\Theta_{2}$ is stable under the map $\sigma$ in
(5.2). Moreover, we have $m_{\sigma(\pi)}\equiv m_{\pi}-2\;\mathrm{mod}\;3.$
Because we have the relation $\gamma_{\sigma(m,n)}=\omega\gamma_{(m,n)}$, we
see that

${\operatorname{Res}}(f_{\sigma(\pi)},\gamma_{\sigma(m,n)})=\omega^{-2}\omega^{2+1}{\operatorname{Res}}(f_{\pi},\gamma_{m,n})=\omega{\operatorname{Res}}(f_{\pi},\gamma_{m,n}).$

In other words, to compute the contributions on all the orbit, it is
sufficient to compute only one residue out of the three, namely:

$\mathcal{C}(\pi)+\mathcal{C}(\sigma(\pi))+\mathcal{C}(\sigma^{2}(\pi))=-2i\pi\omega\\!\\!\sum_{(\eta_{i},\varepsilon_{i})\in
P(f_{\pi})}{\operatorname{Res}}(f_{\pi},-\gamma_{(\eta_{i},\varepsilon_{i})})(\varepsilon_{i}+\eta_{i}\omega-(\eta_{i}+\varepsilon_{i})\omega^{2}).$

We will use the notation $\mathcal{C}([\pi])$ for the contribution of the
whole orbit.

  2. (2)

The invariance by cyclic permutation of the trace implies that

$\mathcal{C}((\alpha_{1},\beta_{1}),(\gamma_{1},\delta_{1}),(\alpha_{2},\beta_{2}),(\gamma_{2},\delta_{2}))=\mathcal{C}((\alpha_{2},\beta_{2}),(\gamma_{2},\delta_{2}),(\alpha_{1},\beta_{1}),(\gamma_{1},\delta_{1})).$

  3. (3)

Moreover, we also get the following relation

$\mathcal{C}((\alpha_{1},\beta_{1}),(\gamma_{1},\delta_{1}),(\alpha_{2},\beta_{2}),(\gamma_{2},\delta_{2}))=\mathcal{C}((-\gamma_{1},-\delta_{1}),(-\alpha_{1},-\beta_{1}),(-\gamma_{1},-\delta_{1})(-\alpha_{2},-\beta_{2})).$

Indeed, using $\gamma_{(-m,-n)}=-\gamma_{(m,n)}$, we see that we get a minus
sign in the residue part of the contribution, but we also get a minus sign in
the “epsilon” part so the resulting contribution is unchanged.

All of this will reduce the number of residues one will have to compute. What
is left is to compute the set $\Theta_{2}$. Consider

$\pi:=\left[\big{(}(\alpha_{1},\beta_{1}),n_{1}\big{)},\big{(}(\gamma_{1},\delta_{1}),n_{1}\big{)},\big{(}(\alpha_{2},\beta_{2}),n_{2}\big{)},\big{(}(\gamma_{2},\delta_{2}),n_{2}\big{)}\right]\in\Theta_{2}.$

Let’s start by introducing some terminology. We’ll call $A$, $B$ and $C$ the
lines through the origin and with direction $(1,1)$, $(-2,1)$ and $(1,-2)$
respectively. We’ll say that a couple $(\alpha,\beta)\in\mathbb{Z}^{2}$ is of
type $A$, $B$ or $C$ if the corresponding point in the $\mathbb{R}^{2}$ plane
is on the corresponding line. The sum of the four (weighted) couples is zero
by definition. By the pigeonhole principle, at least two of the couples have
the same type. Without loss of generality, let’s suppose it is $A$. There are
a few cases to cover.

The case $p\neq\pm\tilde{p}$:

  * •

If $(\alpha_{1},\beta_{1})=-(\gamma_{1},\delta_{1})$ (i.e.,
$(\alpha_{1},\beta_{1})$ and $(\gamma_{1},\delta_{1})$ have same type) then we
must also have $(\alpha_{2},\beta_{2})=-(\gamma_{2},\delta_{2})$ and the
converse also holds. This gives our first family of tuples:

$\displaystyle\pi=\left[\big{(}(\alpha_{1},\beta_{1}),p\big{)},\big{(}(\gamma_{1},\delta_{1}),p\big{)},\big{(}(\alpha_{2},\beta_{2}),\tilde{p}\big{)},\big{(}(\gamma_{2},\delta_{2}),\tilde{p}\big{)}\right],\,\,\,$
$\displaystyle(\alpha_{1},\beta_{1})=-(\gamma_{1},\delta_{1}),\,\,\,\,\,\
(\alpha_{2},\beta_{2})=-(\gamma_{2},\delta_{2}).\quad\quad\quad$

  * •

If $(\alpha_{1},\beta_{1})=-(\alpha_{2},\beta_{2})=(1,1)$, then
$p(\alpha_{1},\beta_{1})+\tilde{p}(\alpha_{2},\beta_{2})=(p-\tilde{p})(1,1).$
This is a non-zero vector on the line $A$. If one of the last two couples is
of type $A$, the so is the last one, so this case is included in the previous
one. We can thus suppose that $(\gamma_{1},\delta_{1})$ is of another type,
say $B$. Because $p\neq\tilde{p}$ and $(1,1)+(-2,1)+(1,-2)=0$, this implies
that $p-\tilde{p}=-p=-\tilde{p}$ which implies $p=0$ and gives a
contradiction. Thus, there are no additional tuple if $p\neq\tilde{p}$.

  * •

The last case also does not produce any other tuples in the case
$p\neq\tilde{p}$.

We know from the previous three points that we only need to consider the
residue for

$\displaystyle\pi_{1}$
$\displaystyle=[((-1,-1),p),((1,1),p),((1,1),\tilde{p}),((-1,-1),\tilde{p})],$
$\displaystyle\pi_{2}$
$\displaystyle=[((-1,-1),p),((1,1),p),((-2,1),\tilde{p}),((2,-1),\tilde{p})],$
$\displaystyle\pi_{3}$
$\displaystyle=[((-1,-1),p),((1,1),p),((1,-2),\tilde{p}),((-1,2),\tilde{p})].$

A straightforward computation of the residues then gives

$\operatorname{tr}(B_{k}^{2})(p,\alpha_{23}/\alpha_{12},\zeta_{2}/\zeta_{1})=\frac{4\pi}{9\sqrt{3}p^{2}}\left(\frac{1}{(\zeta_{2}/\zeta_{1})^{2}}\left(\frac{\alpha_{23}}{\alpha_{12}}\right)^{4}+\frac{3}{1-\zeta_{2}/\zeta_{1}+(\zeta_{2}/\zeta_{1})^{2}}\left(\frac{\alpha_{23}}{\alpha_{12}}\right)^{2}+1\right).$

The case $p=\tilde{p}$: The condition becomes

$(\alpha_{1},\beta_{1})+(\gamma_{1},\delta_{1})+(\alpha_{2},\beta_{2})+(\gamma_{2},\delta_{2})=0,$

where we have

$-(\alpha_{1},\beta_{1}),(\gamma_{1},\delta_{1}),(\alpha_{2},\beta_{2}),-(\gamma_{2},\delta_{2})\in\\{(1,1),(-2,1),(1,-2)\\}.$

The only solutions are given by

$(\alpha_{1},\beta_{1})=-(\alpha_{2},\beta_{2}),\,\,(\gamma_{1},\delta_{1})=-(\gamma_{2},\delta_{2}),$

or by

$(\alpha_{1},\beta_{1})=-(\gamma_{1},\delta_{1}),\,\,(\gamma_{2},\delta_{2})=-(\alpha_{2},\beta_{2}).$

We can compute the residues, using again the three remarks, which gives the
formula

$\operatorname{tr}(B_{k}^{2})(p,\alpha_{23}/\alpha_{12},1)=\left(1-\left(\frac{\alpha_{23}}{\alpha_{12}}\right)^{2}+\left(\frac{\alpha_{23}}{\alpha_{12}}\right)^{4}\right)\frac{4\pi}{9\sqrt{3}p^{2}}.$

This completes the proof. ∎

By combining Theorems 8 and 11 we obtain Corollary 1.1 from the introduction.

##  6\. Continuity of bands & discontinuity of magic parameters

In Proposition 6.1 we show the Hölder 1/2-continuity of the spectrum of the
Hamiltonian in Hausdorff-distance. This implies that the bands do not change
dramatically when changing the twisting angles. This is in contrast to Theorem
4, which says that the coupling coefficients giving rise to flat bands are
rather ill-behaved when changing the twisting angles.

######  Proposition 6.1.

The map
$(\beta,\tilde{\beta},\zeta)\in{\mathbb{C}}^{2}\times{\mathbb{C}}^{2}\times\mathbb{R}^{2}\mapsto\operatorname{Spec}_{L^{2}({\mathbb{C}};{\mathbb{C}}^{4})}(\mathcal{H}(\beta,\zeta))$
is locally Hölder 1/2 continuous in Hausdorff distance.

###### Proof.

Let $L>0$ to be fixed later on. If
$d(\lambda,\operatorname{Spec}_{L^{2}({\mathbb{C}};{\mathbb{C}}^{4})}(\mathcal{H}(\beta,\zeta)))<\varepsilon$
then there is $\varphi_{\varepsilon}$ normalized such that
$\|(\mathcal{H}(\beta,\zeta)-\lambda)\varphi_{\varepsilon}\|\leq\varepsilon\|\varphi_{\varepsilon}\|$.
We then define $\eta_{L}(z):=\chi(z/L)$ with cut-off function $\chi\in
C_{c}^{\infty}({\mathbb{C}};[0,1])$ such that $\chi(z)=1$ for
$|z|_{\infty}\leq 1/2.$ We then set $\eta_{j,L}(z):=\eta_{L}(z-j)$ for
$j\in{\mathbb{Z}}+i{\mathbb{Z}}$ and define
$\alpha^{\uparrow}_{L}:=\sup_{z}\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}|\eta_{j,L}(z)|^{2}=\mathcal{O}(L^{2})$
and
$\alpha_{L}^{\downarrow}:=\inf_{z}\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}|\eta_{j,L}(z)|^{2}=\mathcal{O}(L^{2})$
with $\alpha^{\uparrow}_{L}/\alpha_{L}^{\downarrow}=\mathcal{O}(1).$ Thus, we
have that

$\begin{split}\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|\eta_{j,L}(\mathcal{H}(\beta,\zeta)-\lambda)\varphi_{\varepsilon}\|^{2}&\leq\alpha^{\uparrow}_{L}\|(\mathcal{H}(\beta,\zeta)-\lambda)\varphi_{\varepsilon}\|^{2}\lesssim\alpha_{L}^{\uparrow}\varepsilon^{2}\|\varphi_{\varepsilon}\|^{2}\lesssim\varepsilon^{2}\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|\eta_{j,L}\varphi_{\varepsilon}\|^{2}.\end{split}$

Let $\delta>0$, we then have that using
$\|u+v\|^{2}\leq(1+\delta)\|v\|^{2}+(1+\delta^{-1})\|w\|^{2}$

$\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|(\mathcal{H}(\beta,\zeta)-\lambda)\eta_{j,L}\varphi_{\varepsilon}\|^{2}\leq(1+\delta^{-1})\varepsilon^{2}\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|\eta_{j,L}\varphi_{\varepsilon}\|^{2}+(1+\delta)\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|[\eta_{j,L},\mathcal{H}(\beta,\zeta)]\varphi_{\varepsilon}\|^{2}.$

Now,
$\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|[\eta_{j,L},\mathcal{H}(\beta,\zeta)]\varphi_{\varepsilon}\|^{2}\lesssim\|\varphi_{\varepsilon}\|^{2}\lesssim\frac{1}{L^{2}}\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|\eta_{j,L}\varphi_{\varepsilon}\|^{2}.$
This implies that

$\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|(\mathcal{H}(\beta,\zeta)-\lambda)\eta_{j,L}\varphi_{\varepsilon}\|^{2}\lesssim\Big{(}(1+\delta^{-1})\varepsilon^{2}+\tfrac{1+\delta}{L^{2}}\Big{)}\sum_{j\in{\mathbb{Z}}+i{\mathbb{Z}}}\|\eta_{j,L}\varphi_{\varepsilon}\|^{2}.$

Since this summed inequality holds, there is in particular some $j$ such that

$\|(\mathcal{H}(\beta,\zeta)-\lambda)\eta_{j,L}\varphi_{\varepsilon}\|\lesssim\sqrt{(1+\delta^{-1})\varepsilon^{2}+\tfrac{(1+\delta)}{L^{2}}}\|\eta_{j,L}\varphi_{\varepsilon}\|.$

This implies that for new parameters $\beta^{\prime},\zeta^{\prime}$

$\|(\mathcal{H}(\beta^{\prime},\zeta^{\prime})-\lambda)\eta_{j,L}\varphi_{\varepsilon}\|\lesssim\Big{(}\sqrt{(1+\delta^{-1})\varepsilon^{2}+\tfrac{1+\delta}{L^{2}}}+(|\beta-\beta^{\prime}|+\|U\|_{C^{1}}|\zeta-\zeta^{\prime}|L)\Big{)}\|\eta_{j,L}\varphi_{\varepsilon}\|.$

Choosing
$L=(|\beta-\beta^{\prime}|+\|U\|_{C^{1}}|\alpha-\alpha^{\prime}|)^{-1/2},$ the
result follows. ∎

While the previous implies that the spectrum is well-behaved (in particular
also as we approach incommensurable angles), we now prove Theorem 4, stated in
the introduction saying that the set of magic parameters $\mathcal{B}$ is
_discontinuous_ for varying parameters $\zeta_{1}$ and $\zeta_{2}$.

###### Proof of Theorem 4.

We start by introducing
$h:=\frac{\zeta_{2}}{\zeta_{1}}=3^{k}\frac{r_{1}}{r_{2}},\text{ with
}r_{1},r_{2}\not\equiv 0\operatorname{mod}3.$ We shall then choose a sequence
$(\zeta_{1}^{(n)},\zeta_{2}^{(n)})$ such that
$\frac{\zeta_{1}^{(n)}}{\zeta_{2}^{(n)}}$ yields integers $p_{n}$ in the trace
formula diverging to infinity. More precisely, we choose

$\zeta_{1}^{(n)}=\zeta_{1},\quad\zeta_{2}^{(n)}=3^{k}\frac{r_{1}}{r_{2}}\times\frac{3^{n}}{3^{n}-1}\zeta_{1}^{(n)}=:h_{n}\zeta_{1}^{(n)}\
\Longrightarrow\
(\zeta_{1}^{(n)},\zeta_{2}^{(n)})\to(\zeta_{1},\zeta_{2}),\quad
p_{n}=r_{2}(3^{n}-1).$

Using Theorem 11, then for $n$ large enough, we have
$\zeta_{1}^{(n)}\neq\zeta_{2}^{(n)}$ and thus

$\mathcal{S}_{4}(r,\zeta_{1}^{(n)},\zeta_{2}^{(n)})=\frac{4\pi
p_{n}^{2}}{9\sqrt{3}\zeta_{1}^{4}}\left(\frac{1}{h_{n}^{2}}r^{4}+\frac{3}{1-h_{n}+h_{n}^{2}}r^{2}+1\right).$

This means that we have

$\mathcal{S}_{4}(r,\zeta_{1}^{(n)},\zeta_{2}^{(n)})\sim\underbrace{\frac{4\pi}{9\sqrt{3}\zeta_{1}^{4}}\left(\frac{1}{h^{2}}r^{4}+\frac{3}{1-h+h^{2}}r^{2}+1\right)}_{>0}p_{n}^{2}\to+\infty.\qed$

Since the operator $B_{k}$ depends analytically on
$\frac{\alpha_{23}}{\alpha_{12}}$, we immediately conclude:

######  Proposition 6.2.

Let $\zeta_{1}/\zeta_{2}\in\mathbb{Q}\setminus\\{0\\}$ fixed. The magic
parameters depend continuously on the ratio $\frac{\alpha_{23}}{\alpha_{12}}$.
In Case I, the same is true for magic parameters corresponding to eigenvalues
of $B_{k}$ on subspaces $L^{2}_{-p,\ell}$ for $\ell\in{\mathbb{Z}}_{3}.$

The latter part implies that flat bands of multiplicity $2$ are not destroyed
by varying the ratio $\zeta_{1}/\zeta_{2}.$ This condition is actually sharp,
as numerical experiments suggest that magic parameters of multiplicity 2 in
the case of equal twisting angles $\zeta_{1}=\zeta_{2}$ with
$\alpha_{23}=\alpha_{12}$ are immediately destroyed when varying this ratio.

##  7\. Generic simplicity

In this section, we argue that for $\zeta_{1},\zeta_{2}$ as in Case I, the
magic parameters are either simple or two-fold degenerate for a generic set of
tunnelling potentials satisfying the translational and rotational symmetry of
the tunnelling potentials in the chiral limit.

Since we assume Case I, Remark 6 gives that $T_{0}$ is well-defined on
$L^{2}_{ip}$.

We start by noticing that the existence of an eigenfunction of $T_{0}$ in a
representation $L^{2}_{j,\ell}$ implies the existence of another eigenfunction
in $L^{2}_{j,-\ell+1}$ for $\ell\in\\{0,1\\}.$ In particular, spectrum of
$T_{0}$ in representations $L^{2}_{j,0}$ or $L^{2}_{j,1}$ always gives rise to
degenerate flat bands.

######  Proposition 7.1.

Assume Case I, then $T_{0}:L^{2}_{-p,\ell}\to L^{2}_{-p,\ell}$ is well-defined
for $\ell\in{\mathbb{Z}}_{3}$, and in addition,

$\operatorname{Spec}_{L^{2}_{-p,\ell}}(T_{0})=\operatorname{Spec}_{L^{2}_{-p,-\ell+1}}(T_{0})$

with equality of geometric multiplicity. In particular, if
$1/\alpha\in\operatorname{Spec}_{L^{2}_{-p,\ell}}(T_{0})$ then the
multiplicity of the magic angle is at least $2$. Let
$u_{i}\in\ker_{L^{2}_{-p,i}}(D(\alpha))\setminus\\{0\\}$, then

$u_{0}(z)=(z\pm z_{S})w_{0}(z),\quad u_{1}(z)=z^{2}w_{1}(z),\quad
u_{2}(z)=zw_{2}(z)\quad\text{with }w_{i}\in
C^{\omega}({\mathbb{C}};{\mathbb{C}}).$

If $u\in\ker_{L^{2}_{-p,1}}(D(\alpha))\setminus\\{0\\}$ then
$\wp(z,\frac{4}{3}\pi iw,\frac{4}{3}\pi
iw^{2})u\in\ker_{L^{2}_{-p,0}}(D(\alpha))$, where $\wp(z)$ is the
$\Gamma_{3}$-periodic Weierstrass $p$-function with a double pole at zero.

###### Proof.

We take $v\in L^{2}_{-p,\ell}$ with $T_{0}v=-\lambda v$ and $\lambda\neq 0$.
Multiplying by $2D_{\bar{z}}$ we find that $D(1/\lambda)v=0$. Using the anti-
linear symmetry $Q$, defined in (2.10), we have that $Qv\in L^{2}_{-p,-\ell}$
and $D(1/\lambda)^{*}Qv=0.$ This implies that $(1-1/\lambda
T_{0}^{*})(2D_{\bar{z}}Qv)=0,$ where $(2D_{\bar{z}}Qv)\in L^{2}_{-p,-\ell+1}.$

The expansion of the $u_{i}$ at the zeros follows as in Section $6$ of
[BHZ23], cf. in particular [BHZ23, Theorem $5$]. By the usual theta function
argument, the number of flat bands is equal to two (cf. Lemma 8.4). ∎

We now provide an argument showing that the spectrum of $T_{0}$ in each
representation $L^{2}_{j,\ell}$ is simple.

###  7.1. Generalized potentials

We consider the general class of potentials of the form

$\begin{split}V&=\begin{pmatrix}0&U_{+}&0\\\ U_{-}&0&Y_{+}\\\
0&Y_{-}&0\end{pmatrix}\text{ where }U_{\pm}(\omega z)=\omega
U_{\pm}(z),Y_{\pm}(\omega z)=\omega Y_{\pm}(z)\\\
U_{\pm}&(z+na)=\bar{\omega}^{\pm pn(a_{1}+a_{2})}U_{\pm}(z),\quad
Y_{\pm}(z+na)=\bar{\omega}^{\pm\frac{p\zeta_{2}}{\zeta_{1}}n(a_{1}+a_{2})}Y_{\pm}(z),\end{split}$
(7.1)

with $a=\frac{4\pi i}{3}(a_{1}\omega+a_{2}\omega^{2})\in\Gamma_{3}$. We do not
however assume $\overline{U_{\pm}(z)}=U_{\pm}(\bar{z})$ and
$\overline{Y_{\pm}(z)}=Y_{\pm}(\bar{z})$. For the topology defining the notion
of _genericity_ , it is convenient to use the following Hilbert space of real
analytic potentials defined using the following norm: for fixed $\delta>0$,

$\|V\|_{\delta}^{2}:=\sum_{\pm}\sum_{k\in\Gamma^{*}}(|a^{\pm}_{k}|^{2}+|b^{\pm}_{k}|^{2})e^{2|k|\delta},\
\ \ \ U_{\pm}(z)=\sum_{k\in\Gamma^{*}}a^{\pm}_{k}e^{i\langle z,k\rangle},\ \ \
Y_{\pm}(z)=\sum_{k\in\Gamma^{*}}b^{\pm}_{k}e^{i\langle z,k\rangle}.$

Then we define $\mathscr{V}=\mathscr{V}_{\delta}$ by

$V\in\mathscr{V}\ \Longleftrightarrow\ \text{ $V$ satisfies \eqref{eq:newU},
}\ \|V\|_{\delta}<\infty.$ (7.2)

We also define the more general potential

$V=\begin{pmatrix}W_{1}&U_{+}&X_{+}\\\ U_{-}&W_{2}&Y_{+}\\\
X_{-}&Y_{-}&W_{3}\end{pmatrix}$ (7.3)

with $U_{\pm},Y_{\pm}$ as before and

$W_{i}(\omega z)=\omega W_{i}(z)\text{ and }X_{\pm}(\omega z)=\omega
X_{\pm}(z)$

as well as

$W_{i}(z+na)=W_{i}(z)\text{ and }X_{\pm}(z+na)=\bar{\omega}^{\pm
p(1+\zeta_{2}/\zeta_{1})n(a_{1}+a_{2})}X_{\pm}(z).$

The set of such $V$ with entries satisfying (7.2) are denoted by
$\mathscr{V}_{\text{full}}.$ Such $V$ still satisfy the relevant symmetries
$[\mathscr{L}_{a},V]=0\text{ and }[\mathscr{C},V]=0.$ The potentials $W_{i}$
describe magnetic potentials or strain fields and $X_{\pm}$ describes
tunnelling between top and bottom layers. It is worth noticing that our
argument to show generic simplicity relies on the matrix set
$\mathscr{V}_{\text{full}}$ which requires more potentials than there are
tunnelling potentials in chiral limit TTG.

Before proceeding with the results, we shall explain the simplifications made
along the way. We start with the simplifications needed already in the case of
TBG [BHZ23]:

  * •

We give up on the symmetry $\overline{U_{\pm}(z)}=U_{\pm}(\bar{z}),$ since a
potential satisfying this symmetry does not describe a dense set of functions
in $L^{2}_{j,\ell}$ spaces.

  * •

We do not assume that $U_{+}(z)=U_{-}(-z)$ as is assumed in the idealized
chiral BM model. This choice would introduce a coupling between the components
of the potential that is not general enough to split eigenvalues.

Sticking to these simplifications in the case of TTG allows us to show that
for a generic potential all eigenvalues in each representation are simple.
Giving up on the symmetry $U_{+}(z)=U_{-}(-z)$ breaks the $C_{2}\mathcal{T}$
symmetry of the Hamiltonian which is given by the operator $Q$ in (2.10). This
is not an issue in the case of TBG, since one still has the symmetry
$(Av)(z):=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}\overline{v(z)}$ such that
$AD_{\text{TBG}}(\alpha)A=-D_{\text{TBG}}(\alpha)^{*}$ relating
$D_{\text{TBG}}(\alpha)$ to its adjoint. This symmetry is no longer available
in the case of TTG.

One can still show, along the lines of what has been done in TBG, that for a
generic perturbation the eigenvalues of $T_{k}$ are simple in each subspace
$L^{2}_{-p,\ell},$ but a splitting of representations is not possible. Indeed,
one can only infer that if $T_{0}$ has a simple eigenvalue in each
representation that under a generic perturbation they either split (leading to
flat bands of Chern number -1) right away or move together but are related by
$\theta$ functions. This is then enough to conclude that the Chern number has
to be $-1$ under generic perturbations, see Corollary 9.2.

Our main result of this section is the following theorem:

######  Theorem 12.

Assume Case I. There exists a generic subset (an intersection of open dense
sets), $\mathscr{V}_{0}\subset\mathscr{V}$, with $\mathscr{V}$ defined in
(7.2), such that if $V\in\mathscr{V}_{0}$ then at any magic parameter
$\alpha\in{\mathbb{C}}^{2}$, the operator $T_{0}$ has at most a simple
eigenvalue in each subspace $L^{2}_{-p,\ell}.$ For a generic set
$\mathscr{V}_{0}\subset\mathscr{V}_{\text{full}}$, the flat bands are either
simple bands in $L^{2}_{-p}$ such that

$\dim\ker_{L^{2}_{-p,2}}(D(\alpha))=1\text{ and
}\dim\ker_{L^{2}_{-p,0}}(D(\alpha))=\dim\ker_{L^{2}_{-p,1}}(D(\alpha))=0,$

two-fold degenerate flat bands with

$\dim\ker_{L^{2}_{-p,2}}(D(\alpha))=0\text{ and
}\dim\ker_{L^{2}_{-p,0}}(D(\alpha))=\dim\ker_{L^{2}_{-p,1}}(D(\alpha))=1,$

or three-fold degenerate flat bands

$\dim\ker_{L^{2}_{-p,2}}(D(\alpha))=\dim\ker_{L^{2}_{-p,0}}(D(\alpha))=\dim\ker_{L^{2}_{-p,1}}(D(\alpha))=1,$

where for $u\in\ker_{L^{2}_{-p,2}}(D(\alpha))$, and for precisely one of
$\\{\pm z_{S}\\}$,

$u(z)=zw(z)\text{ and }u(z)=(z\pm z_{S})^{2}v(z)\text{ for }w,v\in
C^{\omega}({\mathbb{C}}).$

In addition, we have in the three-fold degenerate case

$\ker_{L^{2}_{-p,1}}D(\alpha)=\mathbb{C}\wp(z\pm z_{S},\tfrac{4}{3}\pi
iw,\tfrac{4}{3}\pi iw^{2})u(z)$

and

$\ker_{L^{2}_{-p,0}}D(\alpha)=\mathbb{C}\frac{u(z)}{\wp(z\mp
z_{S},\tfrac{4}{3}\pi iw,\tfrac{4}{3}\pi iw^{2})}.$

###  7.2. Proof of generic simplicity

Our proof of Theorem 12 is an adaptation of the argument for generic
simplicity of resonances by Klopp–Zworski [KZ95] – see also [DyZw19, §4.5.5].

Note that we are in Case I, i.e.,
$\\{-q,p,0\\}\equiv\\{p,0\\}\bmod{\mathbb{Z}}_{3}$. The key is to study
simplicity of the eigenvalues of $T_{0}$ on a rotational decomposition of
$L^{2}_{j}$. Hence, let

$R_{k}:=(2D_{\bar{z}}-k)^{-1}:L^{2}_{ip}\to L^{2}_{ip}\cap H^{1}.$

In particular, since $-p\notin\\{p,0\\}$, we may set $k=0$ and see that

$R:=R_{\mathbf{0}}:L^{2}_{ip}\to L^{2}_{ip}\cap H^{1},$ (7.4)

is well defined together with $T_{0}=RV$. We need to decompose (cf. equation
(2.9))

$L^{2}_{ip}=\bigoplus_{\ell=0}^{2}L^{2}_{{-p,\ell}},\ \ \
L^{2}_{{-p,\ell}}\simeq L^{2}(F),$

where $F$ is a fixed fundamental domain of $L^{2}_{-p,\ell}$, i.e.,
${\mathbb{C}}$ quotient by the joint group action $\mathscr{L}$ and
$\mathscr{C}$. For $V\in\mathscr{V}$ and $\ell\in{\mathbb{Z}}_{3}$

$RV:L^{2}_{{-p,\ell}}\to L^{2}_{{-p,\ell}}.$

The next lemma shows that the eigenvalues of $RV$ restricted to each of the
three representations $L^{2}_{-p,\ell}$ with $\ell\in{\mathbb{Z}}_{3}$ are
simple. This argument is analogous to the proof of [BHZ23, Lemma $5.2$] and in
our proof we shall only point out the main differences to the argument
provided there.

######  Lemma 7.2.

In the notation of (7.2) and (7.4), there exists a generic subset
$\mathscr{V}_{\ell}$ of $\mathscr{V}$ such that for $V\in\mathscr{V}_{\ell}$,
the eigenvalues of $RV|_{L^{2}_{{-p,\ell}}}$ are simple for each
$\ell\in{\mathbb{Z}}_{3}.$

###### Proof.

The proof follows the same string of arguments as [BHZ23, Lemma $5.2$] in the
case of twisted bilayer graphene: Consider operators $RV$ acting on
$\mathscr{H}:=L^{2}_{{-p,\ell}}$. The condition that one has to verify is, in
the notation of [BHZ23, (5.14),(5.16)], that for
$w\in\ker_{L^{2}_{-p,\ell}}(T_{k}-\lambda)\setminus\\{0\\}$ and
$\widetilde{w}\in\ker_{L^{2}_{-p,\ell+1}}(T_{k}^{*}-\bar{\lambda})\setminus\\{0\\}$

$\forall\,V\in\mathscr{V}\ \ \langle V{w},R^{*}\widetilde{w}\rangle=0\implies
w\text{ or }\widetilde{w}=0.$

In our case, this leads to the component-wise constraint (with $\circ_{j}$
denoting components of $\bullet=w,\widetilde{w}$)

$\begin{split}&\langle
U_{+}w_{2},R^{*}\widetilde{w}_{1}\rangle_{L^{2}(F)}+\langle
U_{-}w_{1},R^{*}\widetilde{w}_{2}\rangle_{L^{2}(F)}+\langle
Y_{+}w_{3},R^{*}\widetilde{w}_{2}\rangle_{L^{2}(F)}+\langle
Y_{-}w_{2},R^{*}\widetilde{w}_{3}\rangle_{L^{2}(F)}=0,\end{split}$ (7.5)

where $F$ is a fundamental domain. Since $V$ is arbitrary in $F$, this implies
that ${w}$ or $R^{*}\widetilde{w}\equiv 0$ since the elements
$w,\widetilde{w}$ are, by ellipticity of the $D(-1/\lambda)$ operator, real-
analytic functions with full support componentwise. Since it is easy to see
that $\widetilde{w}$ cannot be constant, we readily conclude that indeed $w$
or $\widetilde{w}\equiv 0.$ ∎

Next we show that eigenvalues within different representations can be split.
For twisted bilayer graphene this has been addressed in [BHZ23, Lemma $6.3$].

######  Lemma 7.3.

Assume Case I. For $W$ with $[\mathscr{L}_{a},W]=0$ and $[\mathscr{C},W]=0,$
we have

$\operatorname{Spec}_{L^{2}_{{-p,\ell}}}(RW)\cap
D(\lambda_{0},2r)=\\{\lambda_{0}\\},\ \ \ \ell\in{\mathbb{Z}}_{3},\ r>0,$

and $\lambda_{0}$ is a simple eigenvalue of all $RW|_{L^{2}_{{-p,\ell}}}$.
Then, for every $\varepsilon>0$ there exists $V\in\mathscr{V}_{\text{full}}$,
$\|V\|_{\delta}<\varepsilon$, such that

$\operatorname{Spec}_{L^{2}_{{-p,\ell}}}(R(W+V))\cap
D(\lambda_{0},r)=\\{\widetilde{\lambda}_{\ell}\\},\ \
\widetilde{\lambda}_{0}=\widetilde{\lambda}_{1}\neq\widetilde{\lambda}_{2}$
(7.6)

or

$\operatorname{Spec}_{L^{2}_{{-p,1}}}(R(W+V))\cap
D(\lambda_{0},r)=\operatorname{Spec}_{L^{2}_{{-p,2}}}(R(W+V))\cap
D(\lambda_{0},r)=\\{\widetilde{\lambda}\\}$

with $R(W+V)u=\widetilde{\lambda}u$ for $u\in L^{2}_{-p,2}$, where for
precisely one of $\\{\pm z_{S}\\}$

$u(z)=zw(z)\text{ and }u(z)=(z\pm z_{S})^{2}v(z)\text{ for }w,v\in
C^{\omega}({\mathbb{C}}).$

In addition, we have

$\ker_{L^{2}_{-p,1}}(R(W+\lambda)-\widetilde{\lambda})=\mathbb{C}\wp(z\pm
z_{S},\tfrac{4}{3}\pi iw,\tfrac{4}{3}\pi iw^{2})u(z)$

and

$\ker_{L^{2}_{-p,0}}(R(W+\lambda)-\widetilde{\lambda})=\mathbb{C}\frac{u(z)}{\wp(z\mp
z_{S},\tfrac{4}{3}\pi iw,\tfrac{4}{3}\pi iw^{2})}.$

###### Proof.

We start by observing that

$\begin{split}w_{2}(z\pm z_{S})&=\bar{\omega}w_{2}(\omega(z\pm
z_{S}))=\bar{\omega}w_{2}(\omega z\pm z_{S}\mp\gamma_{2})\\\
&=\bar{\omega}\operatorname{diag}(\omega^{\pm p},1,\omega^{\mp
q})\mathscr{L}_{\mp\gamma_{2}}w_{2}(\omega z\pm z_{S})\\\
&=\bar{\omega}\operatorname{diag}(\omega^{\pm p\mp j},\omega^{\mp
j},\omega^{\mp q\mp j})w_{2}(\omega z\pm z_{S}).\end{split}$ (7.7)

If $q=0$ and $p=1$ then

$\bar{\omega}\operatorname{diag}(\omega^{2p},\omega^{p},\omega^{-q+p})=(\omega,1,1)$

showing that the zero at $z_{S}$ is at least of second order while at $-z_{S}$
we get $(1,\omega,\omega).$ To summarize, if $w_{2}$ vanishes at either $\pm
z_{S}$ then it vanishes to second order.

As in [BHZ23, Lemma $6.3$] we have ${w}_{k},\widetilde{w}_{k}\in L^{2}_{j,k}$,
such that $\langle w_{k},\widetilde{w}_{k}\rangle=1$, and

$(2\lambda_{0}D_{\bar{z}}-W){w}_{k}=0,\ \ \
(2\bar{\lambda}_{0}D_{z}-W^{*})R^{*}\widetilde{w}_{k}=0.$

We can split an eigenvalue with eigenvectors $w_{k}$, if we can find $V$ such
that (see (LABEL:eq:UpUm) for the notation)

$\begin{gathered}\langle
Vw_{0},R^{*}\widetilde{w}_{0}\rangle_{L^{2}(F)}\neq\langle
Vw_{2},R^{*}\widetilde{w}_{2}\rangle_{L^{2}(F)}.\end{gathered}$

If for all (analytic) $V$ the terms were equal it would follow that

$(w_{0})_{\nu}\overline{(R^{*}\widetilde{w}_{0})}_{\ell}=(w_{2})_{\nu}\overline{(R^{*}\widetilde{w}_{2})}_{\ell}.$

The left hand side vanishes at both $\pm z_{S}$ and therefore so does the
right-hand side for every $\nu,\ell.$ If one component $(w_{2})_{\nu}$ of
$w_{2}$ does not vanish at $z_{S}$, then $R^{*}\widetilde{w}_{2}(z_{S})=0$ and
vice versa. We conclude that $w_{2}$ vanishes at either $z_{S}$ or $-z_{S}$
and if not, then $R^{*}\widetilde{w}_{2}$ vanishes at both.

If $w_{2}$ or $R^{*}\widetilde{w}_{2}$ vanishes at both $\pm z_{S}$, then that
function has $5$ zeros, and we have split the eigenvalue as indicated in
(7.6). As the theta function argument would then allow us to construct $5$
flat bands, this would contradict the assumption that
$\ker_{L^{2}_{ip}}(D(\alpha))$ is three-dimensional.

Thus, we assume that neither of the two has a zero at both $\pm z_{S}$. We
shall focus here on $w_{2}$ having a zero at $z_{S}$ or $-z_{S}$ (but not
both) as the case for $R^{*}\widetilde{w}_{2}$ is completely analogous. By the
symmetry argument at the beginning of this proof, the entire vector $w_{2}$
then vanishes at $\pm z_{S}$ to second order.

If $w_{2}$ vanishes at either $\pm z_{S}$ to second order, we can define

$v_{1}(z)=w_{2}(z)\wp(z\pm z_{S},\tfrac{4}{3}\pi iw,\tfrac{4}{3}\pi
iw^{2})\in\ker_{L^{2}_{-p,1}}(D(\alpha)).$

This function vanishes at $0$ to second order and at $\mp z_{S}$ to first
order. Finally,

$v_{0}(z)=v_{1}(z)\wp(z,\tfrac{4}{3}\pi iw,\tfrac{4}{3}\pi
iw^{2})\in\ker_{L^{2}_{-p,0}}(D(\alpha))$

vanishes to first order at $\pm z_{S}$ and to second order at $\mp z_{S}.$ ∎

We find from Theorem 12 that the flat bands in $\ker(D(\alpha)+k)$ of
multiplicity $m$ are the vector spaces $V_{m}$ given for
$w_{i}\in\ker_{L^{2}_{0,i}}(D(\alpha))$ by

$\begin{split}V_{1}(k)&=\\{\zeta F_{k}(z)w_{2}(z);\zeta\in{\mathbb{C}}\\}\\\
V_{2}(k)&=\\{\zeta_{1}F_{k_{1}}(z+z_{S})w_{0}(z)+\zeta_{2}F_{k_{1}}(z-z_{S})w_{0}(z);\zeta_{1},\zeta_{2}\in{\mathbb{C}}\\}\text{
and }\\\
V_{3}(k)&=\\{\zeta_{1}F_{k_{1}}(z-z_{S})w_{0}(z)+\zeta_{2}F_{k_{2}}(z)w_{1}(z)+\zeta_{3}F_{k_{3}}(z)w_{2}(z);\zeta_{1},\zeta_{2},\zeta_{3}\in{\mathbb{C}}\\}.\end{split}$

We can now finish

###### Proof of Theorem 12.

We restrict us to explaining the second part of the proof, as the first part
is even simpler since it does not require the splitting of representations in
Lemma 7.2. We have simplicity of the spectrum of $T_{0}$ on $L^{2}_{ip}$,
modulo the necessary multiplicity of coupled spectrum in
$L^{2}_{-p,0},L^{2}_{-p,1}$. Using then Lemma 7.2 (strictly speaking its
proof) and Lemma 7.3, we conclude that for every $r>0$, the set

$\mathscr{V}_{r}:=\\{V:R(W+V)\text{ satisfies the conclusion of Lemma
}\ref{l:split}\text{ in }{\mathbb{C}}\setminus D(0,r)\\}$

is open and dense. We then obtain $\mathscr{V}_{0}$ by taking the intersection
of $\mathscr{V}_{1/n}$. ∎

##  8\. Band touching for simple magic parameters

Throughout this section, we shall assume without loss of generality that
$p\notin 3\mathbb{Z}$ (see Remark 1). We recall that $\alpha$ is a simple
magic parameter if $E_{1}(k,\alpha)=0$ for any $k$ but
$E_{2}(k_{0},\alpha)\neq 0$ for some $k_{0}$. This leads naturally to the
question of band touching: namely, does there exist
$K_{0}\in{\mathbb{C}}/\Gamma_{3}^{*}$ such that $E_{2}(K_{0},\alpha)=0$? It
turns out this is impossible for twisted bilayer graphene: it was shown in
[BHZ23, Theorem 4] that there was a spectral gap for a simple or doubly-
degenerate magic angle. By this, we mean that the first non-zero band doesn’t
touch the flat bands at any point. Interestingly, in the case of twisted
trilayer graphene, the first two bands always touch. Nevertheless, if we
suppose that the magic parameter $\alpha$ is simple, the point where the two
bands touch is unique. This is the statement of the next theorem and is the
main result of this section.

######  Theorem 13.

For $\alpha\in\mathcal{A}$ simple, there exists a unique
$K_{0}\in\\{-ip,0,iq\\}$ such that

$\dim\ker_{L^{2}_{0}}(D(\alpha)+K_{0})\geq
2,\quad\dim\ker_{L^{2}_{0}}(D(\alpha)+k)=1$

for any $k\notin K_{0}+\Gamma_{3}^{*}$. In other words, the first two bands
touch at only one point $K_{0}\in{\mathbb{C}}/\Gamma_{3}^{*}$. Moreover:

  * •

For Case I, we have $K_{0}=-ip$ if $p=-q$ and $K_{0}=0$ if $q=0$, according to
Proposition 8.5.

  * •

For Case II, we have $K_{0}\equiv-i(m-1)$ where the integer $m$ is
characterized in Proposition 8.6.

We will split the proofs between Case I and II, see the end of the section for
the proofs of Propositions 8.5 and 8.6. The main technique is the theta
function argument, explained in Section 3.3. To use it properly, we first need
to study the zeroes of eigenfunctions for a simple magic parameter.

###  8.1. Vanishing of eigenfunctions for simple magic parameters

We start with a simple lemma characterizing the zeros of zero modes of the
Hamiltonian:

######  Lemma 8.1.

Suppose that $w\in C^{\infty}(\mathbb{C};\mathbb{C}^{2})$ and that
$(D(\alpha)+k)w=0$ for some $k$ and that $w(z_{0})=0$. Then
$w(z)=(z-z_{0})w_{0}(z)$, where $w_{0}\in
C^{\infty}(\mathbb{C};\mathbb{C}^{2})$.

###### Proof.

Since the solution to the elliptic equation is (real) analytic, it suffices to
show that $(2D_{\bar{z}})^{\ell}w(z_{0})=0$ for all $\ell$. Then, writing
$(2D_{\bar{z}})^{\ell}w(z)=(2D_{\bar{z}})^{\ell-1}\left[\left(U-k\right)w\right](z),$
the claim follows by induction on $\ell$. ∎

Using the previous lemma, we find the order of vanishing of elements in
$L^{2}_{r,\ell}$ for $\ell=0,1,2$:

######  Lemma 8.2.

Let $u\in\ker_{L^{2}_{r,\ell}}(D(\alpha))$ for some $r\in{\mathbb{Z}}_{3}$.

  1. (1)

If $\ell=2$, then $u(z)=zw(z)$ for some $w\in
C^{\infty}(\mathbb{C};\mathbb{C}^{3})$.

  2. (2)

If $\ell=1$, then $u(z)=z^{2}w(z)$ for some $w\in
C^{\infty}(\mathbb{C};\mathbb{C}^{3})$.

  3. (3)

If $\ell=0$, and $u(0)=0$, then $u(z)=z^{3}w(z)$ some $w\in
C^{\infty}(\mathbb{C};\mathbb{C}^{3})$.

###### Proof.

For the first case, since $u(\omega z)=\bar{\omega}^{2}u(z)=\omega u(z)$ we
find $u(0)=0$. The result then follows from Lemma 8.1.

For the second case, we also have $u(0)=0$ since $u(\omega
z)=\bar{\omega}u(z)$. Since by Lemma 8.1, we have $u(z)=zu_{1}(z)$ and
$D(\alpha)u=0$ implies $D(\alpha)u_{1}=0$. Using the identity

$\bar{\omega}(zu_{1}(z))=\bar{\omega}u(z)=u(\omega z)=(\omega z)u_{1}(\omega
z),$

we obtain $u_{1}(\omega z)=\omega u_{1}(z)$. This implies $u_{1}(0)=0$. As
$D(\alpha)(zu_{1})=zD(\alpha)u_{1}$, applying Lemma 8.1 again yields
$u_{1}(z)=zu_{2}(z)$ for some $u_{2}\in
C^{\infty}(\mathbb{C};\mathbb{C}^{3})$.

The third case follows from an analogous iteration but with one more step than
the proof of the second case. ∎

We record another useful lemma.

######  Lemma 8.3.

Let $u\in\ker_{L^{2}_{k}}D(\alpha)$ for some $k\in\\{0,-ip,iq\\}$, i.e., at
points in ${\mathbb{C}}/\Gamma^{*}_{3}$ where there exists protected states.
Then $u(0)=0$ implies $\dim\ker_{L^{2}_{k}}D(\alpha)\geq 2$.

###### Proof.

Note that $L^{2}_{k}=\bigoplus_{\ell=0}^{2}L^{2}_{{r,\ell}}$ for
$r=ik\in{\mathbb{Z}}_{3}$. Let $u\in\ker_{L^{2}_{k}}D(\alpha)$ with $u(0)=0$.
Note that there exists a protected state
$\varphi\in\ker_{L^{2}_{r,0}}D(\alpha)$ as $k\in\\{0,-ip,iq\\}$. If
$u,\varphi$ are linearly independent, the lemma is obviously true.

Otherwise $u=c\varphi\in L^{2}_{r,0}$ with $u(0)=0$, so by Lemma 8.2 we obtain
$u=z^{3}w(z)$ for some $w(z)\in C^{\infty}(\mathbb{C};\mathbb{C}^{2})$. This
means that

$\tilde{u}(z)=\wp^{\prime}(z,\tfrac{4}{3}\pi iw,\tfrac{4}{3}\pi iw^{2})u(z)\in
L^{2}_{r,0},\ \ \hat{u}(z)=\wp(z,\tfrac{4}{3}\pi iw,\tfrac{4}{3}\pi
iw^{2})u(z)\in L^{2}_{r,2}$

satisfies $\tilde{u},\hat{u}\in\ker_{L^{2}_{k}}D(\alpha)$, where $\wp(z)$ as
before is the $\Gamma_{3}$-periodic Weierstrass $p$-function with a double
pole at zero and $\wp^{\prime}(z)$ denotes its derivative. This implies that
$\dim\ker_{L^{2}_{k}}D(\alpha)\geq 3$ in this case. ∎

We give a simple but useful lemma, especially in constructing flat bands and
dealing with multiplicity. This reduces the multiplicity of flat bands to the
number of zeroes of eigenfunctions.

######  Lemma 8.4.

For $\\{z_{i}\\}_{1\leq i\leq N}\subset{\mathbb{C}}$, consider
$\\{F_{k}(z-z_{i}):1\leq i\leq N\\}$. If there exists
$k_{0}\in{\mathbb{C}}\setminus\Gamma_{3}^{*}$ such that
$\\{F_{k_{0}}(z-z_{i}):1\leq i\leq N\\}$ is linearly dependent, then there
exist $i\neq j$ such that $z_{i}-z_{j}\in\Gamma_{3}$.

###### Proof.

Note that $F_{k}(z)\in L^{2}({\mathbb{C}}/\Gamma_{3})$. By [BHZ2-22],
$F_{k}(z)$ gives Green kernel of $2D_{\bar{z}}+k$ with

$(2D_{\bar{z}}+k)F_{k}(z)=c(k)\delta_{0}(z),\ \ c(k)=2\pi
i\theta(z(k))/\theta^{\prime}(0)\neq 0\text{ when }k\notin\Gamma_{3}^{*}.$

If there exists $k_{0}\in{\mathbb{C}}\setminus\Gamma_{3}^{*}$ such that
$\sum_{i=1}^{N}c_{i}F_{k_{0}}(z-z_{i})=0$ for some non-trivial $\\{c_{i}\\}$,
applying $2D_{\bar{z}}+k$ yields

$\sum_{i=1}^{N}c_{i}c(k_{0})\delta_{0}(z-z_{i})=0,$

which may only be true if there exist $i\neq j$ such that
$z_{i}=z_{j}\in{\mathbb{C}}/\Gamma_{3}$. ∎

###  8.2. Proof of theorem 13

We first show that the second band always touches the first one. It will be
convenient to split the proof between Case I and Case II, beginning with Case
I. We will use the various lemmas proved in the previous subsection, and we
will also get some additional information on the vanishing point of the
eigenfunctions. Our focus on Case I is no coincidence as in Case II the
location of the zero is not uniquely determined, see Figure 9.

######  Proposition 8.5.

Let $\alpha\in\mathcal{A}$ be simple and assume Case I.

  1. (1)

If $-q\equiv p$ $\mathrm{mod}\ 3$, then the unique protected state
$\psi(z)\in\ker_{L^{2}_{0,0}}D(\alpha)$ has a simple zero at
$-pz_{S}+\Gamma_{3}$; there exists $u\in
E=\mathrm{span}\\{\varphi,\rho\\}\setminus\\{0\\}$ such that $u(z)$ has a
simple zero at $pz_{S}+\Gamma_{3}$, while any $v\in E$ linearly independent of
$u$ has no zero.

  2. (2)

If $q\equiv 0$ $\mathrm{mod}\ 3$, then the protected state
$\varphi(z)\in\ker_{L^{2}_{p,0}}D(\alpha)$ only vanishes simply at
$pz_{S}+\Gamma_{3}$; there exists $u\in
E=\mathrm{span}\\{\psi,\rho\\}\setminus\\{0\\}$ such that $u$ only vanishes
simply at $-pz_{S}+\Gamma_{3}$, while any $v\in E$ linearly independent of $u$
has no zero.

In both cases $u$ is unique (up to a rescaling).

###### Proof.

We only show the first case; the second case is essentially the same. By the
symmetry under $\mathscr{L}_{a},a\in\Gamma_{3}$, we only consider points in
${\mathbb{C}}/\Gamma_{3}$. Recall that by the proof of Proposition 3.5, there
exists $u(z),u^{\prime}(z)\in E$ such that at least one of $\psi(z),u(z)$
vanishes at $z_{S}$ and at least one of $\psi(z),u^{\prime}(z)$ vanishes at
$-z_{S}$.

Recall that by the argument provided in [BHZ2-22, Theorem 3] we know that if
$\alpha\in\mathcal{A}$ is simple, then no eigenstate in
$\ker_{L^{2}_{r,\ell}}D(\alpha)$ with $r,\ell\in{\mathbb{Z}}_{3}$ can vanish
at $z_{0}\notin\\{\pm z_{S},0\\}.$

We also recall that none of the above eigenstates can vanish at $0$, leaving
only $\pm z_{S}$ as possible zeros. Indeed, assume $u\in L^{2}_{r,0}$ for some
$r\in{\mathbb{Z}}_{3}$ and $u(0)=0$. By Lemma 8.2, $u(z)=z^{3}w(z)$ for some
$w\in C^{\infty}({\mathbb{C}};{\mathbb{C}}^{3})$. Using the theta function
argument we construct

$\\{F_{k+ir}(z-z_{S})\wp(z)\tau(ir)u(z),F_{k+ir}(z+z_{S})\wp(z)\tau(ir)u(z)\\}\subset\ker_{L^{2}_{0}}(D(\alpha)+k)$

where $\wp$ is the Weierstrass $p$-function with $\Gamma_{3}$-periodicity,
which has a double pole at zero and simple zeroes at $\pm z_{S}$. The function
$\tau(ir)$ above has the mapping property $L^{2}_{-ir}\to L^{2}_{0}$ by
equation (2.6). By Lemma 8.4, $\dim\ker_{L^{2}_{0}}(D(\alpha)+k)\geq 2$, which
gives a contradiction to $\alpha\in\mathcal{A}$ being simple.

Now we show that $\psi$ cannot vanish at $pz_{S}$. Assume $\psi(pz_{S})=0$.
Write $\psi=(\psi_{1},\psi_{2},\psi_{3})^{T}$. Since $\psi\in L^{2}_{0,0}$, we
have

$\psi_{1}(z+a)=\bar{\omega}^{p(a_{1}+a_{2})}\psi_{1}(z),\
\psi_{2}(z+a)=\psi_{2}(z),\
\psi_{3}(z+a)={\omega}^{q(a_{1}+a_{2})}\psi_{3}(z).$

Define $v(z):=\psi(z+pz_{S})$. Then $v(0)=0$, and

$\begin{split}\mathscr{C}v(z)&=v(\omega z)=\psi(\omega
z+pz_{S})=\psi(z+\bar{\omega}pz_{S})\\\
&=\psi(z+pz_{S}+p\gamma_{1})=v(z+p\gamma_{1})=(\bar{\omega}^{p^{2}}v_{1},v_{2},\omega^{pq}v_{3})^{T}.\end{split}$

As in the proof of Lemma 8.2, this implies that $v(z)=z^{2}w(z)$ for some
$w\in C^{\infty}({\mathbb{C}})$, i.e., $\psi(z)=(z-pz_{S})^{2}w(z-pz_{S})$ has
a double zero at $pz_{S}$. Using the argument in the previous step we conclude
that $\alpha\in\mathcal{A}$ has at least multiplicity two, which is a
contradiction.

Therefore, $\psi$ can only vanish at $-pz_{S}$, and only of order one. Indeed,
if it has a zero of multiplicity two, the theta function argument would imply
that $\alpha\in\mathcal{A}$ has at least multiplicity two, which is again a
contradiction. The same argument shows that there exists $u(z)\in
E=\mathrm{span}\\{\varphi(z),\rho(z)\\}$ such that $u(z)$ only vanishes at
$pz_{S}$ of order one.

It remains to show that any $v\in E$ linearly independent from $u$ has no
zero. As $E\subset L^{2}_{p,0}$, we can show as before that $v(z)$ can at most
admit a simple zero only at $pz_{S}$, as otherwise we get a multiplicity two
$\alpha\in\mathcal{A}$ using the theta function argument. This is in fact also
impossible, as otherwise we have

$\\{F_{k+ip}(z-pz_{S})\tau(ip)u(z),F_{k+ip}(z-pz_{S})\tau(ip)v(z)\\}\subset\ker_{L^{2}_{0}}(D(\alpha)+k)$

and thus $\dim\ker_{L^{2}_{0}}(D(\alpha)+k)\geq 2$, which contradicts the fact
that $\alpha\in\mathcal{A}$ is simple. ∎

We now turn to Case II. The fact that the two bands touch then essentially
follows from lemma 8.3. The point $K_{0}=-i(m-1)$ at which the bands touch is
characterized by which eigenfunction vanishes at the stacking point $-z_{S}$.
This integer $m$ is not made explicit and numerical evidence suggest that it
depends on the magic parameter $\alpha$, see Figure 9.

######  Proposition 8.6.

If $\alpha\in\mathcal{A}$ simple, then there exists $m\in{\mathbb{Z}}_{3}$ for
which we can find $u\in L^{2}_{m,0}$ among $\\{\varphi,\psi,\rho\\}$ with
$u(-z_{S})=0$ and $u(-z_{S})\neq 0$. Moreover, we have
$\dim\ker_{L^{2}_{0}}(D(\alpha)-i(m-1))\geq 2$.

###### Proof.

By Proposition 3.6, for $\alpha\in\mathcal{A}$ we have at least one of the
protected states $\varphi,\psi,\rho$ vanishes at $z_{S}$ and at least one of
$\varphi,\psi,\rho$ vanishes at $-z_{S}$.

We first note that no protected state can vanish at $\pm z_{S}$
simultaneously. Indeed, if there exists a single protected state, say $\psi\in
L^{2}_{0,0}$ (construction for other protected states are similar), that
vanishes at $\pm z_{S}$ simultaneously, then, using the theta function
argument, we can for all $k$ construct

$\displaystyle v_{1,k}=F_{k}(z-z_{S})\psi(z),\
v_{2,k}=F_{k}(z+z_{S})\psi(z)\in\ker_{L^{2}_{0}}(D(\alpha)+k)$

such that $\dim\ker_{L^{2}_{0}}(D(\alpha)+k)\geq 2$ by Lemma 8.4, which
contradicts the assumption that $\alpha$ is simple. Therefore, one of
$\\{\varphi,\psi,\rho\\}$ has to vanish at $z_{S}$ and another one has to
vanish at $-z_{S}$. This proves the first part of the result.

Assume $u\in\\{\varphi,\psi,\rho\\}$ with $u\in L^{2}_{m,0}$ and
$u(-z_{S})=0$. Using the theta function argument in the form of $G_{k}(z)$
(cf. equation (3.8)), we define

$u_{k}(z):={G}_{k+im}(z+z_{S})u(z)\in\ker_{L^{2}_{k}}D(\alpha),\ \forall
k\in{\mathbb{C}}.$

In particular, for $k=-i(m-1)$ we have by [BHZ2-22, Lemma 3.2] that

$\displaystyle u_{-i(m-1)}\in\ker_{L^{2}_{m-1,2}}D(\alpha),\ \
u_{-i(m-1)}(0)=0.$

By Proposition 3.2 we have $\ker_{L^{2}_{m-1,0}}D(\alpha)\neq\\{0\\}$, so
there is a $v_{-i(m-1)}\in L^{2}_{-i(m-1)}$ with $v_{-i(m-1)}\neq
u_{-i(m-1)}$. From the definition of $\tau(k)$ in (2.6) follows that

$v\in L^{2}_{k},\ \ D(\alpha)v=0\implies\tau(-k)v\in L^{2}_{0},\ \
(D(\alpha)+k)v=0.$

Applying this for $k=-i(m-1)$ to $u_{-i(m-1)}$ and $v_{-i(m-1)}$ we obtain the
second part of the result. ∎

Figure 7.  For $\zeta=(1,1)$ and $\alpha_{12}=\alpha_{23}$ all magic
parameters are at least two-fold degenerate. For the largest (with positive
real part) magic parameter the two elements of $\ker_{L^{2}_{0,0}}(D(\alpha))$
and $\ker_{L^{2}_{0,2}}(D(\alpha))$ are shown.

Figure 8. Let $\alpha_{12}=\alpha_{23}.$ For $\zeta=(1,3)$ there are
infinitely many magic parameters that are simple or two-fold degenerate,
respectively. A contour plot of the largest two-fold degenerate zero modes
with $\alpha\approx 1.12+0.65i$ in $\ker_{L^{2}_{0,0}}(D(\alpha))$ (left) and
$\ker_{L^{2}_{0,2}}(D(\alpha))$ (right) is shown.

Figure 9. Let $\alpha_{12}=\alpha_{23}$. For $\zeta_{1}=4$ and $\zeta_{2}=7$
we have simple magic parameters at $1.8999$ and $1.9288$ with simple element
of $\ker_{L^{2}_{0,0}}(D(\alpha))$ illustrated. We observe that the unique
zero appears at two different points $\pm z_{S}.$

Figure 10. Let $\alpha_{12}=\alpha_{23}$. For $\zeta=(1,1)$ (left) we have
two-fold degenerate magic parameter at $\alpha\approx 0.82825$ with simple
wavefunction in $\ker_{L^{2}_{0,0}}(D(\alpha))$ illustrated and for
$\zeta=(1,7)$ (right) we have a two-fold degenerate magic parameter at
$\alpha\approx 1.0179+0.9281i$ with simple element of
$\ker_{L^{2}_{0,0}}(D(\alpha))$ illustrated. There is a double zero at one of
$\\{\pm z_{S}\\}.$

Proposition 3.2 shows that the second band can touch the first band at at
least one point and the location of this touching point is given by
Proposition 8.5 and 8.6 for Case I and II respectively. We now show that in
both Case I and II, when $\alpha\in\mathcal{A}$ is simple, the second band can
only touch the flat band at a unique point among
$\\{-ip,0,iq\\}\subset{\mathbb{C}}/\Gamma_{3}^{*}$ given by Proposition 8.5
and 8.6.

######  Proposition 8.7.

For $\alpha\in\mathcal{A}$ simple, there exists a unique
$K_{0}\in\\{-ip,0,iq\\}$ such that

$\dim\ker_{L^{2}_{0}}(D(\alpha)+k)=1$

for any $k\notin K_{0}+\Gamma_{3}^{*}$.

###### Proof.

Note that $\alpha\in\mathcal{A}$ is simple. Without loss of generality, we can
work within ${\mathbb{C}}/\Gamma_{3}^{*}$. By Proposition 8.5 and 8.6, there
exists $K_{0}\in\\{-ip,0,iq\\}$ such that

$\dim\ker_{L^{2}_{0}}(D(\alpha)+K_{0})=2.$

If there exists $k_{0}\neq K_{0}$ such that

$\dim\ker_{L^{2}_{0}}(D(\alpha)+k_{0})=2,$

then we can take a basis $\\{\varphi_{1},\varphi_{2}\\}$ of the space
$\ker_{L^{2}_{0}}(D(\alpha)+k_{0})$ such that $\varphi_{1}(z)$ has a unique
zero $z_{1}\in{\mathbb{C}}/\Gamma_{3}$ and $\varphi_{2}(z)$ does not vanish.
By Proposition 8.5 and 8.6, we can pick
$\varphi_{K_{0}}\in\ker_{L^{2}_{0}}(D(\alpha)+K_{0})$ such that
$\varphi_{K_{0}}(z)$ does not vanish.

We consider the Wronskian

$W=\det[\varphi_{1},\varphi_{2},\varphi_{K_{0}}].$

Note that the Wronskian satisfies

$(2D_{\bar{z}}+2k_{0}+K_{0})W=0,\ \
W(z+\gamma)=e^{-i\langle\gamma,i(q-p)\rangle}W(z),\ \gamma=\tfrac{4\pi
i}{3}(a_{1}\omega+a_{2}\omega^{2}).$

This yields $W\equiv 0$. Therefore there are $f_{1}(z),f_{2}(z)$ periodic with
respect to $\Gamma_{3}$ such that

$\varphi_{1}=f_{1}(z)\varphi_{2}+f_{2}(z)\varphi_{K_{0}}.$

Applying $D(\alpha)+k_{0}$ to the above equation we obtain

$D_{\bar{z}}f_{1}(z)\cdot\varphi_{2}+(D_{\bar{z}}+k_{0}-K_{0})f_{2}(z)\cdot\varphi_{K_{0}}=0.$

If $\\{\varphi_{1},\varphi_{2}\\}$ is linearly independent for all
$z\in{\mathbb{C}}/\Gamma_{3}$, then we must have have $D_{\bar{z}}f_{1}=0$ and
$(D_{\bar{z}}+k_{0}-K_{0})f_{2}=0$. This yields $f_{2}\equiv 0$ as
$D_{\bar{z}}+K_{0}-k_{0}:L^{2}({\mathbb{C}}/\Gamma_{3})\to
L^{2}({\mathbb{C}}/\Gamma_{3})$ is invertible for $k_{0}\neq K_{0}$ and
$f_{1}(z)$ is a meromorphic function over the torus $\Gamma_{3}$. As
$\varphi_{2}$ has no zero and $\varphi_{1}$ is smooth, $f_{1}$ must be
holomorphic therefore it is a constant. This contradicts to the linear
independence assumption. Therefore, there exists some $z_{*}$ such that
$c_{1}\varphi_{2}(z_{*})+c_{2}\varphi_{K_{0}}(z_{*})=0$ for some $c_{1},c_{2}$
non-zero, as neither $\varphi_{2}$ nor $\varphi_{K_{0}}$ vanishes.

To obtain a contradiction, we can then define the flat band eigenfunctions
$u_{k}\in\ker_{L^{2}_{0}}(D(\alpha)+k)$ as follows:

$u_{k}(z)=c_{1}F_{k-k_{0}}(z-z_{*})\theta(k-K_{0})\varphi_{2}(z)+c_{2}F_{k-K_{0}}(z-z_{*})\theta(k-k_{0})\varphi_{K_{0}}(z).$

Note that this functions is smooth in $z$ even though