suppSupplemental References ††thanks: These two authors contributed equally
††thanks: These two authors contributed equally

# Symmetry-based classification of exact flat bands in single and bilayer
moiré systems

Siddhartha Sarkar  Department of Physics, University of Michigan, Ann Arbor,
MI 48109, USA  Xiaohan Wan  Department of Physics, University of Michigan, Ann
Arbor, MI 48109, USA  Shi-Zeng Lin<EMAIL_ADDRESS>Theoretical Division, T-4 and
CNLS, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
Center for Integrated Nanotechnologies (CINT), Los Alamos National Laboratory,
Los Alamos, New Mexico 87545, USA  Kai Sun<EMAIL_ADDRESS>Department of
Physics, University of Michigan, Ann Arbor, MI 48109, USA

###### Abstract

We study the influence of spatial symmetries on the appearance and the number
of exact flat bands (FBs) in single and bilayer systems with Dirac or
quadratic band crossing points, and systematically classify all possible
number of exact flat bands in systems with different point group symmetries.
We find that a maximum of 6 FBs can be protected by symmetries, and show an
example of 6 FBs in a system with QBCP under periodic strain field of
$\mathcal{C}_{6v}$ point group symmetry. All known examples of exact FBs in
single and bilayer systems fall under this classification, including chiral
twisted bilayer graphene, and new examples of exact FBs are found. We show the
construction of wavefunctions for the highly degenerate FBs, and prove that
any such set of FBs are $\mathds{Z}_{2}$ nontrivial, where all WFs polarized
on one sublattice together have Chern number $C=1$ and WFs polarized on the
other sublattice together have $C=-1$. These bands also satisfy ideal non-
Abelian quantum geometry condition. We further show that, just like in TBG,
topological heavy fermion description of the FBs with higher degeneracy is
possible as long as the Berry curvature distribution is peaked around a point
in the BZ.

Introduction.–The interplay between topology and strong interaction makes
topological flat bands(FBs) ideal for nontrivial states such as fractional
Chern insulators, high $T_{c}$ superconductors etc. One well known example is
twisted bilayer graphene (TBG), where the interference between Dirac points
(DPs) in the two layers give rise to FBs, at some “magic” angles, whose
nontrivial topology is a result of nontrivial Berry phase winding of the
constituent DPs of the two layers. Interestingly, Tarnopolosky et al.
Tarnopolsky  _et al._ (2019) showed that, in the chiral limit of TBG, where
the AA tunneling is set to zero, the middle two bands at the Fermi level
become exactly flat at certain magic angles. They also showed that the exact
wave-functions(WFs) of these two bands are complex analytic functions,
analogous to WFs of the lowest Landau level (LLL) on a torus Haldane and
Rezayi (1985). Just like LLL WFs, the analytic WFs of chiral TBG also carry
Chern number $C=\pm 1$. Since then, exact FBs were shown to appear in
continuum models of a number of different systems at chiral limit such as,
twisted bilayer checkerboard lattice (TBCL) Li  _et al._ (2022), TBG with
spatially alternating magnetic field Becker  _et al._ (2022); Le  _et al._
(2022), single layer system with quadratic band crossing point (QBCP) under
periodic strain field Wan  _et al._ (2023), twisted bilayer Fe-based
superlattices Eugenio and Vafek (2023), and TBG with strong second harmonic
tunneling Becker  _et al._ (2023). Curiously, even though the origin of the
exact FBs in all these systems are same, there are some tangible differences
between the exact FBs reported in these articles. Firstly, in some cases the
number of FBs are doubled Becker  _et al._ (2022); Le  _et al._ (2022); Wan
_et al._ (2023); Becker  _et al._ (2023). Second, in case of TBCL, the Chern
number is $C=\pm 2$ rather than $\pm 1$. High number of FBs as well as high
Chern number can potentially give rise to exotic fractional and correlated
states Neupert  _et al._ (2011); Qi (2011); Wang and Ran (2011); Wang  _et
al._ (2012); Yang  _et al._ (2012); Wu  _et al._ (2013); Andrews and Soluyanov
(2020); Liu  _et al._ (2021); Andrews  _et al._ (2021); Zhang  _et al._
(2019).

However, the different number of FBs and different Chern numbers lead to
puzzles. To get a better understanding, let us describe the origin of the
exact FBs in these systems. The general form of the
$\mathbf{k}\cdot\mathbf{p}$ Hamiltonian with chiral symmetry is

$\mathcal{H}(\mathbf{r})=\begin{pmatrix}\mathbf{0}&\mathcal{D}^{\dagger}(\mathbf{r})\\\
\mathcal{D}(\mathbf{r})&\mathbf{0}\end{pmatrix}$ (1)

near a high symmetry momentum (HSM) $\mathbf{k}_{0}$ in the Brillouin zone
(BZ) of the constituent material having a DP or QBCP, at $\mathbf{k}_{0}$,
protected by $n$-fold rotation symmetry $\mathcal{C}_{nz}$ ($n=3$ for DP in
TBG, $n=3$ or $4$ for QBCP) and an antiunitary symmetry $\mathcal{A}$
($=\mathcal{T}$ for QBCP and $=\mathcal{C}_{2z}\mathcal{T}$ for DP, where
$\mathcal{T}$ is time reversal symmetry; $\mathcal{A}\mathbf{r}=\pm\mathbf{r}$
($\mathcal{A}\mathbf{k}=\mp\mathbf{k}$) when $\mathcal{A}=\mathcal{T}$ and
$\mathcal{C}_{2z}\mathcal{T}$, respectively; also
$\mathcal{A}\mathbf{k}_{0}\equiv\mathbf{k}_{0}$). Note that
$\mathcal{D}(\mathbf{r})$ is $l\times l$ matrix, where $l$ is the number of
layers. The operator $\mathcal{D}(\mathbf{r})$ has kinetic part
$\mathcal{D}_{k}(\mathbf{r})$ consisting of spatial derivatives and potential
term $D_{U}(\mathbf{r};\bm{\alpha})$ which arises from moiré periodic
interlayer tunneling Tarnopolsky  _et al._ (2019); Li  _et al._ (2022); Becker
_et al._ (2022); Le  _et al._ (2022); Eugenio and Vafek (2023); Becker  _et
al._ (2023) or externally applied field Le  _et al._ (2022); Wan  _et al._
(2023), and $\bm{\alpha}$ is the vector of tuning parameters of the potential.
Common to all of these systems is the fact that the kinetic part is
antiholomorphic:
$\mathcal{D}_{k}(\mathbf{r})=(-2i\overline{\partial_{z}})^{m}\mathds{1}$,
where $z=x+iy$, $\mathds{1}$ is the identity matrix, ∗ stands for complex
conjugation, and $m=1$ and $2$ for DP and QBCP, respectively. Furthermore,
under $\mathcal{C}_{nz}$ rotation,
$\mathcal{D}(\mathbf{r})=\mathcal{D}_{k}(\mathbf{r})+\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})$
transforms as
$\mathcal{D}(\mathcal{C}_{nz}\mathbf{r})=(\omega^{*})^{2}\mathcal{D}(\mathbf{r})$,
where $\omega=e^{2\pi i(m+3)/n}$. Consequently, the Hamiltonian satisfies
$\text{Diag}\\{\omega\mathds{1},\omega^{*}\mathds{1}\\}\mathcal{H}(\mathbf{r})\text{Diag}\\{\omega^{*}\mathds{1},\omega\mathds{1}\\}=\mathcal{H}(\mathcal{C}_{nz}\mathbf{r})$.
The antiunitary operator $\mathcal{A}$ acts as
$\sigma_{x}\otimes\mathds{1}\mathcal{H}^{*}(\mathbf{r})\sigma_{x}\otimes\mathds{1}=\mathcal{H}(\mathcal{A}\mathbf{r})$
(see Supplemental Material (SM) SM (2) Sec. S-I for more details). An exact FB
of $\mathcal{H}(\mathbf{r})$ at energy $E=0$ with WF
$\Psi_{\mathbf{k}}(\mathbf{r})$ satisfies
$\mathcal{H}(\mathbf{r})\Psi_{\mathbf{k}}(\mathbf{r})=\mathbf{0}$ for all
$\mathbf{k}$. The construction of such a $\Psi_{\mathbf{k}}(\mathbf{r})$ is as
follows. Note that due to $\mathcal{C}_{nz}$ and chiral symmetry
$\mathcal{S}=\sigma_{z}\otimes\mathds{1}$, the two fold degeneracy of DP or
QBCP at HSM $\mathbf{k}_{0}^{m}$ ($\mathbf{k}_{0}$ gets mapped to
$\mathbf{k}_{0}^{m}$ in the moiré Brillouin zone (mBZ)) remains at $E=0$ for
any $\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})$ that keeps $\mathcal{C}_{nz}$
symmetry. This means that there are always two sublattice polarized WFs
$\Psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r})=\\{\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}),\mathbf{0}\\}$
and
$\Psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r})=\\{\mathbf{0},\psi_{\mathbf{k}_{0}^{m}}^{*}(\mathcal{A}\mathbf{r})\\}$
satisfying
$\mathcal{H}(\mathbf{r})\Psi_{\mathbf{k}_{0}^{m},i}(\mathbf{r})=\mathbf{0}$,
or equivalently
$\mathcal{D}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=\mathbf{0}$. If
there are exact FBs, the WFs can be written as
$\\{\psi_{\mathbf{k}}(\mathbf{r}),\mathbf{0}\\}$ and
$\\{\mathbf{0},\psi_{\mathcal{A}\mathbf{k}}^{*}(\mathcal{A}\mathbf{r})\\}$.
Since the kinetic part $\mathcal{D}_{k}(\mathbf{r})$ is antiholomorphic, the
trial WF is naturally
$\psi_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})=f_{\mathbf{k}}(z)\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$,
where $f_{\mathbf{k}}(z)$ is holomorphic function satisfying
$\overline{\partial_{z}}f_{\mathbf{k}}(z)=0$. The function $f_{\mathbf{k}}(z)$
needs to satisfy Bloch periodicity (translation by moiré lattice vector
$\mathbf{a}^{m}$ gives phase shift $e^{i\mathbf{k}\cdot\mathbf{a}^{m}}$).
However, from Louiville’s theorem, such a holomorphic function must have
poles, making $\psi_{\mathbf{k}}(\mathbf{r})$ divergent, unless
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ has a zero that cancels the pole.
Conversely, if $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ has zero at
$\mathbf{r}_{0}$, a Bloch periodic holomorphic function

$\begin{split}&f_{\mathbf{k}}(z;\mathbf{r}_{0})=e^{i(\mathbf{k}\cdot\mathbf{a}^{m}_{1})z/a_{1}}\frac{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}},\tau\right)}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\\\
&\hskip
36.98866pt=e^{i\mathbf{k}\cdot\mathbf{r}}\tilde{f}_{\mathbf{k}}(\mathbf{r};\mathbf{r}_{0}),\\\
&\tilde{f}_{\mathbf{k}}(\mathbf{r};\mathbf{r}_{0})=e^{-i(\mathbf{b}_{2}^{m}\cdot\mathbf{r})k/b_{2}^{m}}\frac{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}},\tau\right)}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\end{split}$
(2)

with a pole at $\mathbf{r}_{0}$ can be constructed. Here $\vartheta(z,\tau)$
is the Jacobi theta function of the first type Ledwith  _et al._ (2020),
$\mathbf{a}^{m}_{i}$ are lattice vectors, $\mathbf{b}^{m}_{i}$ are the
corresponding reciprocal lattice vectors
($\mathbf{a}^{m}_{i}\cdot\mathbf{b}^{m}_{j}=2\pi\delta_{ij}$),
$a_{i}^{m}=(\mathbf{a}^{m}_{i})_{x}+i(\mathbf{a}^{m}_{i})_{y}$,
$b_{i}^{m}=(\mathbf{b}^{m}_{i})_{x}+i(\mathbf{b}^{m}_{i})_{y}$,
$z_{0}=(\mathbf{r}_{0})_{x}+i(\mathbf{r}_{0})_{y}$, $k=k_{x}+ik_{y}$, and
$\tau=a_{2}^{m}/a_{1}^{m}$. Remarkably, at “magic” values of $\bm{\alpha}$,
the WF $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ has a zero, which allows for
such $f_{\mathbf{k}}(z;\mathbf{r}_{0})$, and in turn gives rise to two exact
FBs. Furthermore, the periodic part
$\tilde{f}_{\mathbf{k}}(\mathbf{r};\mathbf{r}_{0})$ is a holomorphic function
of $k$ Kharchev and Zabrodin (2015); Ledwith  _et al._ (2020); this property
along with the presence of the zero in the WF can be used to prove that wave
functions of this form carry Chern number $C=\pm 1$ (see Wang  _et al._
(2021a), also SM SM (2) S-IV for a proof).

It was found in Wan  _et al._ (2023), that in some cases
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ can have two zeros at
$\mathbf{r}_{0}^{(1)}$ and $\mathbf{r}_{0}^{(2)}$ in the unit cell allowing
construction of two holomorphic functions
$f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(1)})$ and
$f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(2)})$ having poles at
$\mathbf{r}_{0}^{(1)}$ and $\mathbf{r}_{0}^{(2)}$, respectively, giving rise
to 4 FBs. Naturally, the following questions arise: (i) how many zeros
(equivalently, how many FBs) can $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ have?
(ii) do all zeros of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ of allow for
construction of exact FB WFs? Related to this was the observation in Le  _et
al._ (2022); four FBs were found even though
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ had a single quadratic zero, and it
was not clear how to construct the FB WFs fully analytically.

Figure 1: Different number of exact flat bands in single layer system with
QBCP under different periodic strain fields $\tilde{A}(\mathbf{r})$ all having
$p6mm$ symmetry. For (a-b),
$\tilde{A}(\mathbf{r})=-\frac{\alpha}{2}\sum_{n=1}^{3}e^{i(4-n)\phi}\cos\left(\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right)$,
whereas for (c),
$\tilde{A}(\mathbf{r})=-\frac{\alpha}{2}\sum_{n=1}^{3}(-e^{i(1-n)\phi}\cos\left(2\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right)+2e^{i(2-n)\phi}\cos\left((\mathbf{b}_{n}^{m}-\mathbf{b}_{n+1}^{m})\cdot\mathbf{r}\right)+2e^{i(4-n)\phi}\cos\left(\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right))$,
where $\phi=2\pi/3$, $\mathbf{b}^{m}_{1}=\frac{4\pi}{\sqrt{3}a^{m}}(0,1)$ and
$\mathbf{b}^{m}_{2,3}=\frac{4\pi}{\sqrt{3}a^{m}}(\mp\sqrt{3}/2,-1/2)$ are the
reciprocal lattice vectors, and $a^{m}$ is the lattice constant of the
superlattice. The vertical axes in (a-c) are normalized eigen-energies
$\tilde{E}=\frac{E}{|\mathbf{b}^{m}|^{2}}$. The “magic” values of $\alpha$ at
which exact FBs appear are
$\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=0.62,-2.88$ and $2.93$ for
(a), (b) and (c), respectively. The number of FBs in (a), (b) and (c) are 2, 4
and 6, respectively. Density plots of $|\psi_{\Gamma^{m}}(\mathbf{r})|$
(normalized by its maximum) are shown at the lower right corner of (a-c). The
dark spots in the density plots denote zeros, the numbers of which in the unit
cell are 1, 2 and 3 for (a), (b) and (c), respectively. The dashed lines in
the density plots mark the boundary of the unit cell. (d-f) Wilson loop
spectrum $\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of
the flat bands of (a-c), respectively.

In this article, we show that not all zeros of
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ allow for construction of exact FB WF,
only the ones at high symmetry points (HSPs) having
$\mathcal{C}_{n^{\prime}z}$ ($n^{\prime}\geq 2$) symmetry allow such
construction. Using this result, we systematically classify systems with
different PG symmetries, and enumerate the number of FBs they can possess
depending on the multiplicity of the HSPs, where the zeros appear, in the unit
cell. We show that all known examples of exact FBs fall under this
classification in Fig. 3(b). Furthermore, we find that a maximum of 6-FBs is
possible in a system with $\mathcal{C}_{6}$ or $\mathcal{C}_{6v}$ PG symmetry
(unless the system is fine-tuned), we show an example of $\mathcal{C}_{6v}$ PG
symmetry in a single layer system with QBCP under periodic strain field.
Additionally, we clarify the topology of these FBs, and explain the
observation of $C=\pm 2$ bands in Li  _et al._ (2022).

Only the zeros of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ at HSPs give
FBs.–For a single layer system, $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ is a
one-component function. To the lowest order,
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}=\mathbf{r}_{0}+\delta(x,y))=\delta(c_{1}x+c_{2}y)+\mathcal{O}(\delta^{2})$
near a generic zero at $\mathbf{r}_{0}$ of $\psi_{\mathbf{k}_{0}^{m}}$,
whereas
$f_{\mathbf{k}}(\mathbf{r}=\mathbf{r}_{0}+\delta(x,y);\mathbf{r}_{0})\sim\frac{c}{\delta(x+iy)}$
since $f_{\mathbf{k}}$ is holomorphic. Hence,
$\psi_{\mathbf{k}+\mathbf{k}_{0}}(\mathbf{r}=\mathbf{r}_{0}+\delta(x,y))\sim\frac{c_{1}x+c_{2}y}{x+iy}$;
such a function is singular at $\mathbf{r}=\mathbf{r}_{0}$ for generic values
of $c_{1}$ and $c_{2}$, hence $\psi_{\mathbf{k}+\mathbf{k}_{0}}(\mathbf{r})$
is not a good WF. At high symmetry points, however, the relative value of
$c_{1}$ and $c_{2}$ is constrained. For example, around
$\mathcal{C}_{n^{\prime}z}$ ($n^{\prime}\geq 2$) symmetric HSP, if
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+\delta\mathcal{C}_{n^{\prime}z}(x,y))=\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+\delta(x,y))$
(see SM SM (2) S-II for the transformation properties of the WF), then
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}=\mathbf{r}_{0}+\delta(x,y))=\mathcal{O}(\delta^{2})$
implying $c_{1}=c_{2}=0$ which makes the WF $\psi_{\mathbf{k}}(\mathbf{r})$
regular (non-singular). In the bilayer case, where the WF
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ is a two-component function, more care
is required to prove similar result; this is done in the SM SM (2) S-II.

Multiplicity of flat-bands in systems with different PG symmetries.– We
established that among all zeros of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$
only the ones at the HSPs allow construction of exact flat-band WF. However,
notice that there can $m>1$ symmetry related HSPs in one unit cell. Let the
positions of these HSPs be $\mathbf{r}_{0}^{(1)},\dots,\mathbf{r}_{0}^{(m)}$.
Then $m$ zeros would appear in the unit cell together, and hence $m$ complex
analytic Bloch-periodic functions
$f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(1)}),\dots,f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(m)})$
having poles at $\mathbf{r}_{0}^{(1)},\dots,\mathbf{r}_{0}^{(m)}$,
respectively, can be constructed using Eq. 2, and WFs
$\psi_{\mathbf{k}+\mathbf{k}_{0}^{m}}^{(i)}(\mathbf{r})=f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$
($i=1,\dots,m$) form $m$ flat-bands. Moreover, using the antiunitary symmetry
$\mathcal{A}$, we can construct another $m$ flat-band WFs. This implies that
if the HSP, where the zero of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ appears,
has multiplicity $m$, there are $2m$ flat-bands. We show examples of this for
space group $p6mm$ in Fig. 1, where we consider a single layer system with
QBCP at $\mathbf{k}_{0}^{m}=\Gamma$ point under periodic strain field Wan  _et
al._ (2023) with
$\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})=\tilde{A}(\mathbf{r})=A_{x}(\mathbf{r})+iA_{y}(\mathbf{r})$
where $A_{x}=u_{xx}-u_{yy}$ and $A_{y}=u_{xy}$ are moiré periodic shear strain
fields. In Fig. 1(a), (b) and (c), we show that for specific forms of
$\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})$ (see the caption of Fig. 1 for
details), at a critical value of $\alpha$ (the strength of the strain field),
2, 4 and 6 FBs appear, respectively. Correspondingly, the plots of
$\psi_{\Gamma}(\mathbf{r})$ show zeros at HSPs with multiplicity 1, 2 (related
by $\mathcal{C}_{2z}$) and 3 (related by $\mathcal{C}_{3z}$), respectively, in
the unit cell, in agreement with the argument for the number of FBs presented
above.

Figure 2: Four flat bands in TBG Hamiltonian under alternating magnetic field
Le  _et al._ (2022). (a), (b) Moiré unit cell and Brillouin zone,
respectively. (c) Band structure for $\alpha=1.379e^{i0.254\pi}$ showing 4
exact FBs. (d) Wilson loop spectrum
$\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of the 4 flat
bands. (e) Density plot of $|\psi_{K^{m}}(\mathbf{r})|$ (normalized by its
maximum), which has a single quadratic zero at the AB position, as shown in
(f). (g) Density plot of $|\psi_{\Gamma^{m}}(\mathbf{r})|$ (normalized by its
maximum), which has two linear zeros at the AB and BA positions, as shown in
(h). The dashed lines in the density plots mark the boundary of the unit cell.
Figure 3:  (a) Flowchart of the construction of exact flat band wave-functions
starting from the wave-function at $\mathbf{k}_{0}^{m}$. (b) List of all
possible point groups of $\mathbf{k}_{\text{max}}$ for space groups where DP
or QBCP is possible, possible HSPs $\mathbf{r}_{0}$ (having $C_{2}$ or higher
symmetry), their multiplicity in the unit cell, number of flat bands if zeros
appear at these HSPs, the number of parameters need to be varied (co-
dimension) to achieve the flat bands. The last column shows examples of
different cases.

It is important to note that the momentum $\mathbf{k}_{0}^{m}$, where the DP
or QBCP is, may not be the momentum ($\mathbf{k}_{\text{max}}$) with highest
symmetry in the BZ. An example of this is chiral TBG: the $K$ point (where the
DP is) has magnetic little co-group $6^{\prime}$ (generated by
$\mathcal{C}_{3z}$, $\mathcal{C}_{2x}$ and $\mathcal{C}_{2z}\mathcal{T}$).
But, due to a unitary particle-hole symmetry $P=\mathds{1}\otimes i\sigma_{y}$
($P\mathbf{r}=-\mathbf{r}$, $P\mathbf{k}=-\mathbf{k}$ and
$P\mathcal{H}(\mathbf{r})=-\mathcal{H}(-\mathbf{r})P$) Song  _et al._ (2019);
Hejazi  _et al._ (2019) and chiral symmetry
$\mathcal{S}=\sigma_{z}\otimes\mathds{1}$ ($\mathcal{S}\mathbf{r}=\mathbf{r}$,
$\mathcal{S}\mathbf{k}=\mathbf{k}$ and
$\mathcal{S}\mathcal{H}(\mathbf{r})=-\mathcal{H}(\mathbf{r})\mathcal{S}$), an
effective $\mathcal{C}_{2z}$ symmetry can be obtained:
$\mathcal{C}_{2z}=P\mathcal{S}=\sigma_{z}\otimes i\sigma_{y}$ (this symmetry
was found previously in Wang  _et al._ (2021b), and was called intra-valley
inversion symmetry). As a consequence, in chiral TBG,
$\mathbf{k}_{\text{max}}=\Gamma$ with magnetic little co-group $6221^{\prime}$
(generated by $\mathcal{C}_{3z}$, $\mathcal{C}_{2z}$, $\mathcal{C}_{2x}$ and
$\mathcal{T}$; see SM Sec. S-I for a review of symmetries of chiral TBG). In
these cases, where $\mathbf{k}_{0}^{m}\neq\mathbf{k}_{\text{max}}$, if there
is a zero at some HSP $\mathbf{r}_{0}$ of the WF
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$, one should first construct
$\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})=f_{\mathbf{k}_{\text{max}}-\mathbf{k}_{0}^{m}}(z,\mathbf{r}_{0})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$
using the zero, and then check the number of zeros of
$\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$ in the moiré unit cell (since
$\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$ has more symmetry than
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$) to construct the flat-band WFs.
Interestingly, this resolves the ambiguity in the construction of 4 flat-band
WFs in TBG under alternating magnetic field in Le  _et al._ (2022) as is
described below. The $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian around
$\mathbf{k}_{0}=K$ considered in this paper is the same as in Eq. 1 with

$\begin{split}\mathcal{D}(\mathbf{r})&=\begin{pmatrix}-2i\overline{\partial_{z}}&\alpha
U(\mathbf{r})\\\ \alpha
U(-\mathbf{r})&-2i\overline{\partial_{z}}\end{pmatrix},\\\
U(\mathbf{r})&=e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{i\phi}e^{-i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{-i\phi}e^{-i\mathbf{q}_{3}\cdot\mathbf{r}},\end{split}$
(3)

where $\phi=\frac{2\pi}{3}$, $\mathbf{q}_{1}=\frac{4\pi}{3a^{m}}(0,-1)$,
$\mathbf{q}_{2,3}=\frac{2\pi}{3a^{m}}(\pm\sqrt{3},1)$, $a^{m}$ is the moiré
lattice constant (see Fig. 2(a-b) for details of the moiré lattice vectors and
reciprocal lattice vectors). This is exactly the same as chiral TBG model
considered in Tarnopolsky  _et al._ (2019) except that here, unlike
Tarnopolsky  _et al._ (2019), $\alpha$ is allowed to take complex values. It
was found in Le  _et al._ (2022) that, at certain critical complex values of
$\alpha$, there are 4 exact FBs (see Fig. 2(c-d)). At these critical $\alpha$,
the WF $\psi_{\mathbf{k}_{0}^{m}=K}(\mathbf{r})$ has a single quadratic zero
at the AB position $\mathbf{r}_{0}^{(1)}=a^{m}/\sqrt{3}(-1,0)$ as shown in
Fig. 2(e-f). Therefore, it was pointed out in Le  _et al._ (2022) how to
construct the FB WFs using complex analytic function Eq. 2 fully. However, the
WF
$\psi_{\mathbf{k}_{\text{max}}=\Gamma}(\mathbf{r})=f_{\mathbf{k}_{\text{max}}-\mathbf{k}_{0}^{m}}(z,\mathbf{r}_{0}^{(1)})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=f_{\mathbf{q}_{3}}(z,\mathbf{r}_{0}^{(1)})\psi_{K}(\mathbf{r})$
has two linear zeros at AB ($\mathbf{r}_{0}^{(1)}=a^{m}/\sqrt{3}(-1,0)$) and
BA ($\mathbf{r}_{0}^{(2)}=a^{m}/\sqrt{3}(1,0)$) positions, as shown in Fig.
2(g-h). This is due to the fact that this Hamiltonian has the effective
$\mathcal{C}_{2z}=P\mathcal{S}$ as discussed earlier (note that when $\alpha$
is complex, $\mathcal{C}_{2x}$ is broken, hence this system has
$\mathcal{C}_{6}$ point group and $p6$ space group symmetry), hence
$\psi_{\Gamma}(\mathbf{r})$ is $\mathcal{C}_{2z}$ symmetric, and AB and BA
positions are related by $\mathcal{C}_{2z}$. Using these two zeros, 4 FB WFs
can be constructed as follows:

$\begin{split}\Psi_{\mathbf{k}}^{(1)}(\mathbf{r})&=f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(1)})f_{\mathbf{q}_{3}}(z,\mathbf{r}_{0}^{(1)})\begin{bmatrix}\psi_{K}(\mathbf{r})\\\
\mathbf{0}\end{bmatrix},\\\
\Psi_{\mathbf{k}}^{(2)}(\mathbf{r})&=f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(2)})f_{\mathbf{q}_{3}}(z,\mathbf{r}_{0}^{(1)})\begin{bmatrix}\psi_{K}(\mathbf{r})\\\
\mathbf{0}\end{bmatrix},\\\
\Psi_{\mathbf{k}}^{(3)}(\mathbf{r})&=\sigma_{x}\otimes\mathds{1}{\Psi_{\mathcal{A}\mathbf{k}}^{(1)}}^{*}(\mathcal{A}\mathbf{r}),\\\
\Psi_{\mathbf{k}}^{(4)}(\mathbf{r})&=\sigma_{x}\otimes\mathds{1}{\Psi_{\mathcal{A}\mathbf{k}}^{(2)}}^{*}(\mathcal{A}\mathbf{r}).\end{split}$
(4)

We summarize the procedure to find the number of exact FBs from the number of
zeros of $\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$ as well as the
construction of the exact FB WFs in the flow-chart in Fig. 3(a). This insight
also allows us to predict the number of possible FBs in systems with all
different point group symmetries that protect DP or QBCP. This is done via the
following steps: (i) we list all possible point group symmetries (little co-
group Bradley and Cracknell (2010)) of $\mathbf{k}_{\text{max}}$ for space
groups that allow DP or QBCP, (ii) then we list the HSPs (maximal wyckoff
positions) of these space groups, their site-symmetry groups and multiplicity
in the unit cell, (iii) the number of possible symmetry-dictated FBs is then
twice the multiplicity of the HSPs. These results are tabulated in Fig. 3(b).
The table also contains all examples of exact flat bands (in single and
bilayer systems) found in the literature Tarnopolsky  _et al._ (2019); Le  _et
al._ (2022); Becker  _et al._ (2022, 2023); Wan  _et al._ (2023); Eugenio and
Vafek (2023) as well as some new ones: 2 FBs in systems with QBCP under
periodic strain having $p3$ or $p4$ space group symmetries, 4 FBs in systems
with QBCP under periodic strain having $p4$, $p4mm$ and $p4gm$ space group
symmetry (details of these are given in SM SM (2) S-VI), and 6 flat bands in
systems with QBCP under periodic strain having $p6mm$ space group symmetry
(Fig. 1(c)). Furthermore, since in 2D space groups, the maximum multiplicity
of a HSP is 3, the maximum number of FBs can be 6 (not counting fine tuned
cases where zeros can accidentally appear at HSPs not related by symmetries).

Construction of WFs at $\mathbf{k}_{\text{max}}$.–Notice that, in the
construction of $2m$ ($m>1$) FB WFs from zeros of
$\psi_{\mathbf{k}_{\text{max}}}$ at $\mathbf{r}_{0}^{(i)}$, the functions
$f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})$ are independent of each other for all
$\mathbf{k}$ except $\mathbf{k}=\mathbf{0}$ since
$f_{\mathbf{k}=\mathbf{0}}(z;\mathbf{r}_{0}^{(i)})=1$ for all $i=1,\dots,m$.
Therefore, the construction above gives $2m$ independent WFs for
$\mathbf{k}\neq\mathbf{k}_{\text{max}}$, but only 2 independent WFs at
$\mathbf{k}_{\text{max}}$; hence the other $2(m-1)$ WFs at
$\mathbf{k}_{\text{max}}$ need to be obtained separately. To this end, note
that the WFs
$\psi_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(i)}(\mathbf{r})=f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$
are independent, but not mutually orthogonal for
$\mathbf{k}\neq\mathbf{k}_{\text{max}}$. We can modify the later $(m-1)$ WFs
to be orthogonal to $\psi_{\mathbf{k}}^{(1)}(\mathbf{r})$, and define new WFs:

$\tilde{\psi}_{\mathbf{k}}^{(i)}(\mathbf{r})=\begin{cases}\psi_{\mathbf{k}}^{(1)}(\mathbf{r})&\text{
if }i=1,\\\
\psi_{\mathbf{k}}^{(i)}(\mathbf{r})-\frac{\langle\psi_{\mathbf{k}}^{(1)}|\psi_{\mathbf{k}}^{(i)}\rangle}{\langle\psi_{\mathbf{k}}^{(1)}|\psi_{\mathbf{k}}^{(1)}\rangle}\psi_{\mathbf{k}}^{(1)}(\mathbf{r})&\text{
if }i\neq 1,\\\ \end{cases}$ (5)

where $\langle f|g\rangle=\int_{\text{unit
cell}}d^{2}\mathbf{r}f^{*}(\mathbf{r})g(\mathbf{r})$. Remembering
$f_{\mathbf{k}}(z;\mathbf{r}_{0})$ is product of a holomorphic function
$\tilde{f}_{\mathbf{k}}(\mathbf{r};\mathbf{r}_{0})$ of $k=k_{x}+ik_{y}$ and a
phase factor $e^{i\mathbf{k}\cdot\mathbf{r}}$, and
$\tilde{f}_{\mathbf{0}}(\mathbf{r};\mathbf{r}_{0})=1$, we can expand these
functions as
$f_{\mathbf{k}}(z,\mathbf{r}_{0})=1+\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0})k+i(kz^{*}+k^{*}z)/2+\mathcal{O}(k^{2})$,
where
$\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0})\equiv\frac{1}{2}[(\partial_{k_{x}}-i\partial_{k_{y}})\tilde{f}_{\mathbf{k}}(\mathbf{r};\mathbf{r}_{0})]|_{\mathbf{k}=\mathbf{0}}$.
Plugging this back in Eq. 5, we find close to $\mathbf{k}_{\text{max}}$

$\tilde{\psi}_{\mathbf{k}_{\text{max}}+\mathbf{k}}^{(i)}(\mathbf{r})\approx\begin{cases}(1+k\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)})+(kz^{*}+k^{*}z)/2)\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})&\text{
if }i=1,\\\
k(\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(i)})-\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)})-\langle\psi_{\mathbf{k}_{\text{max}}}|(\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(i)})-\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)}))\psi_{\mathbf{k}_{\text{max}}}\rangle)\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})&\text{
if }i\neq 1.\\\ \end{cases}$ (6)

Continuing these WFs to $\mathbf{k}_{\text{max}}$, we get

$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(i)}(\mathbf{r})=\begin{cases}\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})&\text{
if }i=1,\\\
(\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(i)})-\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)})-\langle\psi_{\mathbf{k}_{\text{max}}}|(\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(i)})-\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)}))\psi_{\mathbf{k}_{\text{max}}}\rangle)\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})&\text{
if }i\neq 1.\\\ \end{cases}$ (7)

Note that
$\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(i)})-\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)})$
is a holomorphic function of $z=x+iy$ even though
$\tilde{f}_{\mathbf{k}}(\mathbf{r};\mathbf{r}_{0})$ is not holomorphic, hence
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(i)}(\mathbf{r})$ satisfy
$\mathcal{D}(\mathbf{r})\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(i)}(\mathbf{r})=\mathbf{0}$.
Furthermore, the last $(m-1)$ WFs at $\mathbf{k}_{\text{max}}$ are
automatically orthogonal to $\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$. We
can further mutually orogonalize them using Gram-Schmidt orthogonalization
procedure. The $m$ WFs at $\mathbf{k}_{\text{max}}$ polarized on the other
sublattice are, as usual,
$(\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(i)}(\mathcal{A}\mathbf{r}))^{*}$.
We compare WFs obtained in this way with numerically obtained WFs at
$\mathbf{k}_{\text{max}}=\Gamma$ for the 4 and 6 FBs systems (see SM SM (2)
S-III for details).

Co-dimension of the tuning parameters.–Since
$\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$ is complex $l$ component vector
($l$ is the number of layers), the vector equation
$\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r}_{0})=\mathbf{0}$ is actually a set
of $2l$ equations (real and imaginary part of each component need to be zero).
Therefore, in principle, $2l$ parameters has to be varied to obtain FBs; hence
the co-dimension is generically $2l$. In case of TBG, since one of the
components of $\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$ is automatically
zero at the corners of the moiré unit cell (AB or BA stacking positions) due
to $\mathcal{C}_{3z}$ symmetry (see Tarnopolsky  _et al._ (2019), also SM SM
(2) S-II for details), in all cases, that we are considering, the co-dimension
is 2, unless there are extra symmetries. If the HSP $\mathbf{r}_{0}$ has
mirror-type symmetry $\mathcal{M}$, such that the representation of the
composition $\mathcal{M}\mathcal{A}$ is diagonal with entries $\pm 1$ (by
mirror-type we mean: (i) the symmetry does not commute with $\mathcal{C}_{nz}$
($n\geq 3$), (ii) it commutes with $\mathcal{A}$ whose representation is
$\sigma_{x}\otimes\mathds{1}\mathcal{K}$ and (iii)
$\mathcal{M}^{2}=\mathds{1}$. Note that $C_{2x}$ is also mirror-type according
to this.), then the components of $\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$
are purely real or imaginary at $\mathbf{r}=\mathbf{r}_{0}$, which reduces the
co-dimension to 1. In case of single layer QBCP, mirror symmetry necessarily
has representation $\rho(\mathcal{M})=\sigma_{x}$. In case of chiral TBG, one
can construct an $\mathcal{MA}$ like symmetry
$P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}$ (under which
$P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}(x,y)=(x,-y)$ and
$P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}(k_{x},k_{y})=(-k_{x},k_{y})$,
hence it is a symmetry of the AB and BA stacking point in real space, as well
as a symmetry of the $K^{m}$ point), that has representation
$\rho(P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T})=\sigma_{z}\otimes\sigma_{z}\mathcal{K}$.
Therefore, in QBCP systems with mirror-type symmetry as well as in case of
chiral TBG, the co-dimension is 1. Furthermore, in case of chiral TBG with
alternating magnetic field Le  _et al._ (2022), $\mathcal{C}_{2x}$ is broken,
hence the co-dimension is 2. The co-dimension of tuning parameters for systems
with different symmetries are tabulated in table in Fig. 3(b). Our results for
co-dimensions are in agreement with the ones given in Sheffer  _et al._
(2023).

Topology and quantum geometry of the FBs.–As was mentioned earlier, it can be
proven that the sublattice polarized exact FB WFs have Chern number $C=\pm 1$
using the fact that the WFs are (anti)holomorphic in $k$ (the WFs on one
sublattice are holomorphic, the WFs on the other sublattice are
antiholomorphic). This is why the Wilson loop spectra (see SM SM (2) S-IV for
details), for the case where there are only two FBs, have $\pm 1$ winding as
shown in Fig. 1(d). Surprisingly, for the cases of 4 or 6 FBs, the Wilson loop
spectrum still has winding $\pm 1$ as can be seen from Figs. 1(e-f) and 2(d).
This implies that the $m$ WFs polarized to the same sublattice together
actually have Chern number $|C|=1$, instead of $|C|=m$ as we would have
naïvely guessed from the fact that these WFs are holomorphic functions. We
explain this in the following. As mentioned earlier, the sublattice polarized
WFs
$\psi_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(i)}(\mathbf{r})=f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$
are not mutually orthogonal. When we mutually orthogonalize them (similar to
Eq. 5), $m-1$ of the WFs have a singularity; they are proportional to $k$ near
$\mathbf{k}=0$ (Eq. 6). Therefore, they have a Berry phase winding around
$\mathbf{k}=0$, this Berry phase winding cancels the Berry phase winding along
the edge of the BZ to give Chern number $C=0$ (see SM SM (2) S-IV for
details). As a consequence, only one of the sublattice polarized WFs carry
$|C|=1$, the others are topologically trivial. This result is in agreement
with Song  _et al._ (2021), where it was shown that for systems with an
antiunitary particle-hole symmetry $\mathcal{P}$ that satisfies
$\mathcal{P}^{2}=-\mathds{1}$, if a set of bands symmetric around zero energy
has $4p+2$ ($p\in\mathds{N}$) number of Dirac points between them, then they
are $\mathds{Z}_{2}$ nontrivial having Chern number $C=\pm 1$. Indeed, the
systems considered here have $\mathcal{P}$ symmetry; for QBCP systems it can
be constructed as $\mathcal{P}=\mathcal{S}\mathcal{T}$, whereas in TBG it is
$\mathcal{P}=P\mathcal{C}_{2z}\mathcal{T}$. Interestingly, TBCL Li  _et al._
(2022) has 2 FBs with Chern numbers $\pm 2$ which seemingly is in disagreement
with our result. However, it can be shown that the TBCL Hamiltonian is just
two copies of a single layer system with QBCP under a periodic moiré potential
with $p4mm$ space group symmetry unfolded to a BZ twice the size of single
layer moiré BZ, and this BZ doubling gives results in the doubling of the
Chern number (see SM SM (2) S-VII for details). Also, if one replaces the
$p4mm$ symmetric tunneling potential with a $p4gm$ potential in TBCL (this can
be achieved by simply making the amplitude purely imaginary), one finds 4 FBs
with $C=\pm 2$ instead of 2 FBs at some magic angle (see SM SM (2) S-VII);
this doubling of number of FBs in $p4gm$ case is due to “band sticking”
inherent in nonsymmorphic space groups Dresselhaus  _et al._ (2008). We
further show in the  SM SM (2) S-V, that the FBs with higher degeneracy
satisfy ideal non-Abelian trace condition:
$\text{tr}(g^{mn}_{\alpha\beta}(\mathbf{k}))=\sum_{n}^{\prime}\sum_{\alpha\in\\{x,y\\}}g^{nn}_{\alpha\alpha}(\mathbf{k})=|\text{tr}(F_{xy}^{mn}(\mathbf{k}))|=|\sum_{n}^{\prime}F_{xy}^{nn}(\mathbf{k})|$,
where $g_{ij}^{mn}(\mathbf{k})$ is the non-Abelian Fubini-Study metric and
$F_{xy}^{mn}(\mathbf{k})$ is the non-Abelian Berry curvature, and the prime on
the sum over $n$ means that the sum is restricted to FBs polarized on one
sublattice.

Figure 4: Topological heavy fermion model for the system with 6 flat bands in
Fig. 1(c). (a) The band structure at $\tilde{\alpha}=2.53$, which is a little
deviated from the “magic” value $\tilde{\alpha}=2.93$, at which we have 6
exact FBs. The overlaps between the Bloch states and the trial Wannier
functions are represented by the red circles. (b) Band structure given by the
topological heavy fermion model (blue crosses) compared to that of QBCP 6
nearly FBs full Hamiltonian (black lines). (c) The band structure of QBCP 6
nearly FBs model plotted within $\tilde{E}\in(-0.4,0.4)$. The irreps at HSMs
are shown. Total Berry curvature distribution $\Omega(\mathbf{k}_{n})$
(normalized by average Berry curvature) of the 3 FBs polarized on one
sublattice is shown at the bottom left corner of (c).

Topological heavy fermion (THF) model.–Due to the antiunitary particle-hole
symmetry $\mathcal{P}$ discussed above, any set of bands symmetric about the
charge neutrality point is topological Song  _et al._ (2021). As a
consequence, a tight binding description of these bands is never possible.
However, in the case of TBG, due to the fact that the Berry curvature
distribution is peaked at $\Gamma$ point in the mBZ, it was shown in Song and
Bernevig (2022) that hybridization of 2 atomic limit HF bands with 4
topological conduction bands (having nontrivial winding) at $\Gamma$ point can
describe the 2 topological FBs of TBG. This THF model keeps all the relevant
symmetries of TBG and captures the correct topology of the bands. We find that
similar THF description of the high number of FBs discussed in this article is
possible as long as the Berry curvature distribution has a pronounced peak
around some point in the BZ. We show an example of this in Fig. 4 for the
system (space group $G=p6mm$) with 6 FBs in Fig. 1(c). The irreps of the 6 FBs
at the HSMs are $\Gamma_{5}\oplus 2\Gamma_{6}-2M_{1}\oplus 2M_{2}\oplus
M_{3}\oplus M_{4}-K_{1}\oplus K_{2}\oplus 2K_{3}$, which is not a linear
combination of elementary band representations (EBR) Bradlyn  _et al._ (2017).
On the other hand, the two lowest higher energy bands have representations
$\Gamma_{1}$ and $\Gamma_{2}$ at $\Gamma$. Furthermore, replacement of one
$\Gamma_{6}$ with $\Gamma_{1}\oplus\Gamma_{2}$ allows for band representation
$BR=(A_{1}\uparrow G)_{1a}\oplus(A_{2}\uparrow G)_{1a}\oplus(E_{1}\uparrow
G)_{1a}\oplus(E_{2}\uparrow G)_{1a}$ Aroyo  _et al._ (2011, 2006a, 2006b).
This and the fact that the Berry curvature distribution of the 6 FBs being
peaked at $\Gamma$ (Fig 4(c)) suggests a THF model composed of local orbitals
having band representation $BR$ and topological conduction bands (we will
refer to them as $c$-bands) with representation $\Gamma_{6}$. Thus we
construct 6 Gaussian-type trial Wannier function that transform as $A_{1}$,
$A_{2}$, $E_{1}$ and $E_{2}$ representations of PG $C_{6v}$, and calculate the
overlaps of them with the 6 FB WFs. As shown in Fig. 4(a), the Wannier
functions are fully supported by the 6 FBs everywhere in the BZ except near
the $\Gamma$ point, where the $A_{1}$ and $A_{2}$ type Wannier functions are
supported by the lowest higher energy bands. We obtain the maximally localized
Wannier functions (we will refer to them as $f$-orbitals) using the program
Wannier90 Marzari and Vanderbilt (1997); Souza  _et al._ (2001); Pizzi  _et
al._ (2020) (see SM SM (2) S-VIII for details), these MLWFs are extremely
localized as can be seen in Fig. 4(d). Projecting the Hamiltonian into the
space spanned by $f$-orbitals and the remaining two $c$-bands, we obtain the
following single-particle Hamiltonian:

$\displaystyle\hat{\mathcal{H}}=\sum_{|\mathbf{k}|<\Lambda_{c}}H^{(c)}_{ab}(\mathbf{k})c^{\dagger}_{a}(\mathbf{k})c_{b}(\mathbf{k})+\sum_{\mathbf{R}}H^{(f)}_{\alpha\beta}f^{\dagger}_{\alpha}(\mathbf{R})f_{\beta}(\mathbf{R})+$
$\displaystyle\phantom{\hat{\mathcal{H}}}\sum_{|\mathbf{k}|<\Lambda_{c},\mathbf{R}}\left(H^{(cf)}_{a\alpha}(\mathbf{k})e^{-i\mathbf{k}\cdot\mathbf{R}-|\mathbf{k}|^{2}\lambda^{2}/2}c^{\dagger}_{a}(\mathbf{k})f_{\alpha}(\mathbf{R})+\text{h.c.}\right),$
$\displaystyle H^{(c)}(\mathbf{k})\approx
c_{1}(k_{x}^{2}-k_{y}^{2})\sigma_{x}-2c_{1}k_{x}k_{y}\sigma_{y},$
$\displaystyle H^{(f)}\approx\begin{pmatrix}m\sigma_{z}&\mathbf{0}_{2\times
4}\\\ \mathbf{0}_{4\times 2}&\mathbf{0}_{4\times 4}\end{pmatrix},$
$\displaystyle
H^{(cf)}(\mathbf{k})\approx\begin{pmatrix}-ic_{2}(k_{x}-ik_{y})&c_{2}(k_{x}-ik_{y})&0&\gamma&0&0\\\
-ic_{2}(k_{x}+ik_{y})&-c_{2}(k_{x}+ik_{y})&\gamma&0&0&0\\\ \end{pmatrix}$ (8)

where $\alpha\in\\{A_{1},A_{2},E_{1}^{1},E_{1}^{2},E_{2}^{1},E_{2}^{2}\\}$
goes over the $f$-orbitals, $m$ and $\gamma$ set the gap between
$\Gamma_{1}$-$\Gamma_{2}$ and $\Gamma_{6}$-$\Gamma_{6}$ irreps, respectively
(see Fig. 4(b)), $c_{1}$ and $c_{2}$ are real constants. The spectrum of
$\hat{\mathcal{H}}$ matches very well with that of the full Hamiltonian (Fig.
4(b)).

Conclusion.–In this letter, we systematically studied single and bilayer
continuum Hamiltonians under moiré periodic modulations, and classified all
possible number of exact flat bands that can appear in systems with different
point group symmetries, and examined their topology. We not only showed that
all single and bilayer systems in literature (to our knowledge) fall into this
classification, but also found several new examples of model systems that can
host high number of flat bands as well as high Chern number. In the future, it
would be interesting to study the interplay between interactions and high
number of flat bands that can enrich the correlated states in these systems.

Acknowledgements.–This work was supported in part by Air Force Office of
Scientific Research MURI FA9550-23-1-0334 and the Office of Naval Research
MURI N00014-20-1-2479 (XW, SS and KS) and Award N00014-21-1-2770 (XW and KS),
and by the Gordon and Betty Moore Foundation Award N031710 (KS). The work at
LANL (SZL) was carried out under the auspices of the U.S. DOE NNSA under
contract No. 89233218CNA000001 through the LDRD Program, and was supported by
the Center for Nonlinear Studies at LANL, and was performed, in part, at the
Center for Integrated Nanotechnologies, an Office of Science User Facility
operated for the U.S. DOE Office of Science, under user proposals
$\\#2018BU0010$ and $\\#2018BU0083$.

## References

  * Tarnopolsky  _et al._ (2019) G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Physical review letters 122, 106405 (2019).
  * Haldane and Rezayi (1985) F. D. M. Haldane and E. H. Rezayi, Physical Review B 31, 2529 (1985).
  * Li  _et al._ (2022) M.-R. Li, A.-L. He, and H. Yao, Physical Review Research 4, 043151 (2022).
  * Becker  _et al._ (2022) S. Becker, T. Humbert, and M. Zworski, arXiv preprint arXiv:2208.01628 (2022).
  * Le  _et al._ (2022) C. Le, Q. Zhang, C. Fan, X. Wu, and C.-K. Chiu, arXiv preprint arXiv:2210.13976 (2022).
  * Wan  _et al._ (2023) X. Wan, S. Sarkar, S.-Z. Lin, and K. Sun, Physical Review Letters 130, 216401 (2023).
  * Eugenio and Vafek (2023) P. M. Eugenio and O. Vafek, SciPost Physics 15, 081 (2023).
  * Becker  _et al._ (2023) S. Becker, T. Humbert, and M. Zworski, arXiv preprint arXiv:2306.02909 (2023).
  * Neupert  _et al._ (2011) T. Neupert, L. Santos, C. Chamon, and C. Mudry, Physical review letters 106, 236804 (2011).
  * Qi (2011) X.-L. Qi, Physical review letters 107, 126803 (2011).
  * Wang and Ran (2011) F. Wang and Y. Ran, Physical Review B 84, 241103 (2011).
  * Wang  _et al._ (2012) Y.-F. Wang, H. Yao, C.-D. Gong, and D. Sheng, Physical Review B 86, 201101 (2012).
  * Yang  _et al._ (2012) S. Yang, Z.-C. Gu, K. Sun, and S. D. Sarma, Physical Review B 86, 241112 (2012).
  * Wu  _et al._ (2013) Y.-L. Wu, N. Regnault, and B. A. Bernevig, Physical review letters 110, 106802 (2013).
  * Andrews and Soluyanov (2020) B. Andrews and A. Soluyanov, Physical Review B 101, 235312 (2020).
  * Liu  _et al._ (2021) Z. Liu, A. Abouelkomsan, and E. J. Bergholtz, Physical Review Letters 126, 026801 (2021).
  * Andrews  _et al._ (2021) B. Andrews, T. Neupert, and G. Möller, Physical Review B 104, 125107 (2021).
  * Zhang  _et al._ (2019) Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Physical Review B 99, 075127 (2019).
  * SM (2) “See Supplemental Material … for details,” .
  * Ledwith  _et al._ (2020) P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath, Physical Review Research 2, 023237 (2020).
  * Kharchev and Zabrodin (2015) S. Kharchev and A. Zabrodin, Journal of Geometry and Physics 94, 19 (2015).
  * Wang  _et al._ (2021a) J. Wang, J. Cano, A. J. Millis, Z. Liu, and B. Yang, Physical Review Letters 127, 246403 (2021a).
  * Song  _et al._ (2019) Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Physical review letters 123, 036401 (2019).
  * Hejazi  _et al._ (2019) K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Physical Review B 99, 035111 (2019).
  * Wang  _et al._ (2021b) J. Wang, Y. Zheng, A. J. Millis, and J. Cano, Phys. Rev. Res. 3, 023155 (2021b).
  * Bradley and Cracknell (2010) C. Bradley and A. Cracknell,  _The mathematical theory of symmetry in solids: representation theory for point groups and space groups_ (Oxford University Press, 2010).
  * Sheffer  _et al._ (2023) Y. Sheffer, R. Queiroz, and A. Stern, Physical Review X 13, 021012 (2023).
  * Song  _et al._ (2021) Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Phys. Rev. B 103, 205412 (2021).
  * Dresselhaus  _et al._ (2008) M. S. Dresselhaus, G. Dresselhaus, and A. Jorio, “Applications of group theory to the physics of solids,” (2008).
  * Song and Bernevig (2022) Z.-D. Song and B. A. Bernevig, Physical review letters 129, 047601 (2022).
  * Bradlyn  _et al._ (2017) B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Aroyo  _et al._ (2011) M. I. Aroyo, J. M. Perez-Mato, D. Orobengoa, E. Tasci, G. de la Flor, and A. Kirov, Bulg. Chem. Commun 43, 183 (2011).
  * Aroyo  _et al._ (2006a) M. I. Aroyo, J. M. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga, A. Kirov, and H. Wondratschek, Zeitschrift für Kristallographie-Crystalline Materials 221, 15 (2006a).
  * Aroyo  _et al._ (2006b) M. I. Aroyo, A. Kirov, C. Capillas, J. Perez-Mato, and H. Wondratschek, Acta Crystallographica Section A: Foundations of Crystallography 62, 115 (2006b).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Physical review B 56, 12847 (1997).
  * Souza  _et al._ (2001) I. Souza, N. Marzari, and D. Vanderbilt, Physical Review B 65, 035109 (2001).
  * Pizzi  _et al._ (2020) G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, _et al._ , Journal of Physics: Condensed Matter 32, 165902 (2020).

Supplemental Material

###### Contents

  1. S-1 Continuum chiral symmetric moiré Hamiltonian and symmetries
    1. S-1.1 $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian near a Dirac or quadratic band crossing point
    2. S-1.2 Adding moiré potential and Brillouin zone folding
    3. S-1.3 Single layer system with QBCP under periodic potential
    4. S-1.4 Twisted bilayer systems with DP or QBCP
    5. S-1.5 Symmetries of chiral TBG
  2. S-2 Properties of zero energy eigen wave-function $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$
    1. S-2.1 Zero energy eigen wave-function $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ and its transformation under symmetries
    2. S-2.2 Zeros of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$, and constraints from symmetries
      1. S-2.2.1 Single layer system with QBCP
      2. S-2.2.2 Twisted bilayer graphene
  3. S-3 Comparing analytically obtained wave-functions $\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(i)}$ in Eq. 7 of main text with numerical ones for single layer QBCP under periodic strain field with $p6mm$ symmetry
  4. S-4 Topology of the flat-bands
    1. S-4.1 Topology of two flatbands – single flat band on each sub-lattice
    2. S-4.2 Topology of more than two flatbands – multiple flat bands on each sub-lattice
    3. S-4.3 Details of Wilson loop spectrum calculation
  5. S-5 Ideal quantum geometry of the flat-bands
  6. S-6 Examples
    1. S-6.1 2 flat bands in single layer system with QBCP under moiré potential with space group symmetry $p3$
    2. S-6.2 2 flat bands in single layer system with QBCP under moiré potential with space group symmetry $p4$
    3. S-6.3 4 flat bands in single layer system with QBCP under moiré potential with space group symmetry $p4$
    4. S-6.4 4 flat bands in single layer system with QBCP under moiré potential with space group symmetry $p4gm$
    5. S-6.5 4 flat bands in single layer system with QBCP under moiré potential with space group symmetry $p4mm$
  7. S-7 Twisted bilayer checkerboard lattice (TBCL) is two uncoupled copies of single layer QBCP system under periodic strain field
    1. S-7.1 TBCL type Hamiltonian with 4 flat bands
  8. S-8 Details of the topological heavy fermion (THF) model shown in Fig. 4 of the main text
    1. S-8.1 Maximally localized Wannier functions for the $f$-electrons
    2. S-8.2 The $c$-electrons
    3. S-8.3 The single particle Hamiltonian

##  S-1 Continuum chiral symmetric moiré Hamiltonian and symmetries

###  S-1.1 $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian near a Dirac or quadratic
band crossing point

A Dirac point (DP) or quadratic band crossing point (QBCP) (at a high symmetry
momentum (HSM) $\mathbf{k}_{0}$) is a two fold degeneracy in the band
structure. Let the corresponding eigenstates be $|\mathbf{k}_{0},A\rangle$ and
$|\mathbf{k}_{0},B\rangle$. In case of graphene, $A$ and $B$ are the
sublattices. However, in general case, they denote the two eigenvectors. These
degeneracies are generally protected by $\mathcal{C}_{nz}$ ($n\in\\{3,4\\}$)
and an antiunitary symmetry $\mathcal{A}$, which in case of QBCP is time
reversal symmetry $\mathcal{T}$ and in case of DP is
$\mathcal{C}_{2z}\mathcal{T}$ ( note that here and elsewhere we will use the
letter $n$ for $n$-fold rotation symmetry $\mathcal{C}_{nz}$ that protects the
DP or QBCP. If there are other rotation symmetries in a particular system, we
will use $n^{\prime}$ or write it explicitly). The eigenstates
$|\mathbf{k}_{0},A\rangle$ and $|\mathbf{k}_{0},B\rangle$ transform as
2-dimensional (co-)irreducible representation (irrep) under these two
symmetries. Furthermore, we will restrict ourselves to chiral ($\mathcal{S}$)
symmetric limit. Here and in the rest of the article, we use the following
gauge of the two eigenstates:

$\begin{split}\mathcal{C}_{nz}\\{|\mathbf{k}_{0},A\rangle,|\mathbf{k}_{0},B\rangle\\}&=\\{|\mathbf{k}_{0},A\rangle,|\mathbf{k}_{0},B\rangle\\}\rho(\mathcal{C}_{nz}),\\\
\mathcal{A}\\{|\mathbf{k}_{0},A\rangle,|\mathbf{k}_{0},B\rangle\\}&=\\{|\mathbf{k}_{0},A\rangle,|\mathbf{k}_{0},B\rangle\\}\rho(\mathcal{A}),\\\
\mathcal{S}\\{|\mathbf{k}_{0},A\rangle,|\mathbf{k}_{0},B\rangle\\}&=\\{|\mathbf{k}_{0},A\rangle,|\mathbf{k}_{0},B\rangle\\}\rho(\mathcal{S}),\\\
\end{split}$ (S1)

where

$\rho(\mathcal{C}_{nz})=\begin{pmatrix}\omega&0\\\
0&\omega^{*}\end{pmatrix},\rho(\mathcal{A})=\sigma_{x},\rho(\mathcal{S})=\sigma_{z},$
(S2)

where $\omega=\exp(i2\pi(m+3)/n)$ and $m$ is an integer which takes values
$m=1$ for DP and $m=2$ for QBCP and $n\in\\{3,4\\}$. Near the HSM, we choose
similar gauge:

$\begin{split}\mathcal{C}_{nz}\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}&=\\{|\mathbf{k}_{0}+\mathcal{C}_{nz}\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathcal{C}_{nz}\mathbf{k},B\rangle\\}\rho(\mathcal{C}_{nz}),\\\
\mathcal{A}\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}&=\\{|\mathbf{k}_{0}+\mathcal{A}\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathcal{A}\mathbf{k},B\rangle\\}\rho(\mathcal{A}),\\\
\mathcal{S}\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}&=\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}\rho(\mathcal{S}),\\\
\end{split}$ (S3)

where $\mathcal{A}\mathbf{k}=\pm\mathbf{k}$ for
$\mathcal{A}=\mathcal{C}_{2z}\mathcal{T}$ and $\mathcal{A}=\mathcal{T}$
respectively. Note that the HSM $\mathbf{k}_{0}$ may have some extra
symmetries such as mirror reflection $\mathcal{M}$ and/or two-fold rotation
symmetry $\mathcal{C}_{2z}$. For mirror symmetry, since
$[\mathcal{M},\mathcal{C}_{nz}]\neq 0$ ($n\in\\{3,4\\}$),
$[\mathcal{M},\mathcal{A}]=0$ and $\mathcal{M}^{2}=\mathds{1}$, the only
representation it can have is $\rho(\mathcal{M})=\sigma_{x}$. On the other
hand, $\mathcal{C}_{2z}$ satisfies $[\mathcal{C}_{2z},\mathcal{C}_{nz}]=0$
($n\in\\{3,4\\}$), $[\mathcal{C}_{2z},\mathcal{A}]=0$ and
$\mathcal{C}_{2z}^{2}=\mathds{1}$. Hence, representation of $\mathcal{C}_{2z}$
can only be $\pm\mathds{1}$. In four-fold rotation symmetric case it is
necessarily $\rho(\mathcal{C}_{2z})=\rho^{2}(\mathcal{C}_{4z})=-\mathds{1}$.
In $6$-fold rotation symmetric case there are two choices:
$\mathcal{C}_{2z}=\mathds{1}$ (example: $\Gamma_{5}$ representation in $p6mm$
space group; here we are following the notation given in Bilbao
Crystallography Server
aroyo2011crystallographys,aroyo2006bilbaoIs,aroyo2006bilbaoIIs) and
$\mathcal{C}_{2z}=-\mathds{1}$ (example: $\Gamma_{6}$ representation in $p6mm$
space group).

Then, the $\mathbf{k}\cdot\mathbf{p}$ Hamiltonian near HSM $\mathbf{k}_{0}$ is
given by:

$\hat{H}=\Omega\sum_{\mathbf{k}}\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}\mathcal{H}(\mathbf{k})\\{\langle\mathbf{k}_{0}+\mathbf{k},A|,\langle\mathbf{k}_{0}+\mathbf{k},B|\\}^{T}.$
(S4)

where $\Omega$ is the volume of the system. Since the Hamiltonian invariant
under $\mathcal{C}_{nz}$ and $\mathcal{A}$, we have

$\begin{split}&\mathcal{C}_{nz}\hat{H}\mathcal{C}_{nz}^{-1}=\Omega\sum_{\mathbf{k}}\mathcal{C}_{nz}\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}\mathcal{H}(\mathbf{k})\\{\langle\mathbf{k}_{0}+\mathbf{k},A|,\langle\mathbf{k}_{0}+\mathbf{k},B|\\}^{T}\mathcal{C}_{nz}^{-1}\\\
&\phantom{\mathcal{C}_{nz}\hat{H}\mathcal{C}_{nz}^{-1}}=\Omega\sum_{\mathbf{k}}\\{|\mathbf{k}_{0}+\mathcal{C}_{nz}\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathcal{C}_{nz}\mathbf{k},B\rangle\\}\rho(\mathcal{C}_{nz})\mathcal{H}(\mathbf{k})\rho(\mathcal{C}_{nz})^{\dagger}\\{\langle\mathbf{k}_{0}+\mathcal{C}_{nz}\mathbf{k},A|,\langle\mathbf{k}_{0}+\mathcal{C}_{nz}\mathbf{k},B|\\}^{T}\\\
&\phantom{\mathcal{C}_{nz}\hat{H}\mathcal{C}_{nz}^{-1}}=\Omega\sum_{\mathbf{k}}\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}\rho(\mathcal{C}_{nz})\mathcal{H}(\mathcal{C}_{nz}^{-1}\mathbf{k})\rho(\mathcal{C}_{nz})^{\dagger}\\{\langle\mathbf{k}_{0}+\mathbf{k},A|,\langle\mathbf{k}_{0}+\mathbf{k},B|\\}^{T}\\\
&\phantom{\mathcal{C}_{nz}\hat{H}\mathcal{C}_{nz}^{-1}}=\hat{H},\\\
\Rightarrow&\rho(\mathcal{C}_{nz})\mathcal{H}(\mathbf{k})\rho(\mathcal{C}_{nz})^{\dagger}=\mathcal{H}(\mathcal{C}_{nz}\mathbf{k}).\end{split}$
(S5)

Similarly, we also have

$\begin{split}\rho(\mathcal{A})\mathcal{H}^{*}(\mathbf{k})\rho(\mathcal{A})^{\dagger}=\mathcal{H}(\mathcal{A}\mathbf{k}),\\\
\end{split}$ (S6)

where the complex conjugation ∗ is due to antiunitarity of $\mathcal{A}$.
Furthermore, since $\hat{H}$ anticommutes with chiral symmetry $\mathcal{S}$,
we have

$\begin{split}\rho(\mathcal{S})\mathcal{H}(\mathbf{k})\rho(\mathcal{S})^{\dagger}=-\mathcal{H}(\mathbf{k}).\\\
\end{split}$ (S7)

These symmetries together imply that to the lowest order the Hamiltonian
matrix has the following form:

$\begin{split}\mathcal{H}_{\text{DP}}(\mathbf{k})&=\begin{pmatrix}0&c^{*}k^{*}\\\
ck&0\end{pmatrix},\\\
\mathcal{H}_{\text{QBCP}}(\mathbf{k})&=\begin{pmatrix}0&c_{1}^{*}(k_{x}^{2}-k_{y}^{2})-2ic_{2}^{*}k_{x}k_{y}\\\
c_{1}(k_{x}^{2}-k_{y}^{2})+2ic_{2}k_{x}k_{y}&0\end{pmatrix},\end{split}$ (S8)

where $k=k_{x}+ik_{y}$, and $c$, $c_{1}$ and $c_{2}$ are complex constants.
Note that for QBCPs protected by $\mathcal{C}_{3z}$ or $\mathcal{C}_{6z}$,
$c_{1}=c_{2}$ making the Hamiltonian isotropic. However, for QBCP protected by
$\mathcal{C}_{4z}$, generically $c_{1}\neq c_{2}$. But, in this article, we
will consider assume isotropic limit for all cases. Then, the complex phase in
$c=|c|e^{i\varphi}$ and $c_{1}=|c_{1}|e^{i\varphi}$ can be absorbed by
rotating the coordinate system $k\rightarrow ke^{-i\varphi}$ and $k\rightarrow
ke^{-i\varphi/2}$ in case of DP and QBCP, respectively. Then dropping the
multiplicative constant $|c|$ and $|c_{1}|$, we obtain the following form of
the Hamiltonian:

$\mathcal{H}_{m}(\mathbf{k})=\begin{pmatrix}0&(k^{m})^{*}\\\
k^{m}&0\end{pmatrix},$ (S9)

with $m=1$ and $m=2$ for DP and QBCP, respectively.

Next, we perform Fourier transform on the basis states
$|\mathbf{k}_{0}+\mathbf{k},A\rangle$ and
$|\mathbf{k}_{0}+\mathbf{k},B\rangle$:

$\begin{split}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle&=\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{r}}|\mathbf{k}_{0}+\mathbf{k},\alpha\rangle,\\\
|\mathbf{k}_{0}+\mathbf{k},\alpha\rangle&=\frac{1}{\Omega}\int
d^{2}\mathbf{r}e^{i\mathbf{k}\cdot\mathbf{r}}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle\text{
}\alpha\in\\{A,B\\}.\end{split}$ (S10)

Note that under lattice translations $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ the
state $|\mathbf{k}_{0}+\mathbf{k},\alpha\rangle$ transforms as
$T_{\mathbf{a}_{i}}|\mathbf{k}_{0}+\mathbf{k},\alpha\rangle=e^{-i(\mathbf{k}_{0}+\mathbf{k})\cdot\mathbf{a}_{i}}|\mathbf{k}_{0}+\mathbf{k},\alpha\rangle$.
This implies

$\begin{split}T_{\mathbf{a}_{i}}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle&=\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{r}}T_{\mathbf{a}_{i}}|\mathbf{k}_{0}+\mathbf{k},\alpha\rangle\\\
&=\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{r}}e^{-i(\mathbf{k}_{0}+\mathbf{k})\cdot\mathbf{a}_{i}}|\mathbf{k}_{0}+\mathbf{k},\alpha\rangle\\\
&=e^{-i\mathbf{k}_{0}\cdot\mathbf{a}}|\mathbf{k}_{0},\mathbf{r}+\mathbf{a},\alpha\rangle\end{split}$
(S11)

In this basis, the Hamiltonian becomes

$\begin{split}\hat{H}&=\Omega\sum_{\mathbf{k}}\\{|\mathbf{k}_{0}+\mathbf{k},A\rangle,|\mathbf{k}_{0}+\mathbf{k},B\rangle\\}\mathcal{H}_{m}(\mathbf{k})\\{\langle\mathbf{k}_{0}+\mathbf{k},A|,\langle\mathbf{k}_{0}+\mathbf{k},B|\\}^{T}\\\
&=\frac{1}{\Omega}\sum_{\mathbf{k}}\int d^{2}\mathbf{r}\int
d^{2}\mathbf{r}^{\prime}\\{|\mathbf{k}_{0},\mathbf{r},A\rangle,|\mathbf{k}_{0},\mathbf{r},B\rangle\\}e^{i\mathbf{k}\cdot\mathbf{r}}\mathcal{H}_{m}(\mathbf{k})e^{-i\mathbf{k}\cdot\mathbf{r}^{\prime}}\\{\langle\mathbf{k}_{0},\mathbf{r}^{\prime},A|,\langle\mathbf{k}_{0},\mathbf{r}^{\prime},B|\\}^{T}\\\
&=\frac{1}{\Omega}\sum_{\mathbf{k}}\int d^{2}\mathbf{r}\int
d^{2}\mathbf{r}^{\prime}\\{|\mathbf{k}_{0},\mathbf{r},A\rangle,|\mathbf{k}_{0},\mathbf{r},B\rangle\\}e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}^{\prime})}\mathcal{H}_{m}(\mathbf{r}^{\prime})\\{\langle\mathbf{k}_{0},\mathbf{r}^{\prime},A|,\langle\mathbf{k}_{0},\mathbf{r}^{\prime},B|\\}^{T}\\\
&=\frac{1}{\Omega}\int d^{2}\mathbf{r}\int
d^{2}\mathbf{r}^{\prime}\\{|\mathbf{k}_{0},\mathbf{r},A\rangle,|\mathbf{k}_{0},\mathbf{r},B\rangle\\}\Omega\delta^{(2)}(\mathbf{r}-\mathbf{r}^{\prime})\mathcal{H}_{m}(\mathbf{r}^{\prime})\\{\langle\mathbf{k}_{0},\mathbf{r}^{\prime},A|,\langle\mathbf{k}_{0},\mathbf{r}^{\prime},B|\\}^{T}\\\
&=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathbf{r},A\rangle,|\mathbf{k}_{0},\mathbf{r},B\rangle\\}\mathcal{H}_{m}(\mathbf{r})\\{\langle\mathbf{k}_{0},\mathbf{r},A|,\langle\mathbf{k}_{0},\mathbf{r},B|\\}^{T},\end{split}$
(S12)

where $\delta^{(2)}(\mathbf{r})$ is the 2-dimensional Dirac delta function and

$\mathcal{H}_{m}(\mathbf{r})=\mathcal{H}_{m}(\mathbf{k}\rightarrow-i\bm{\nabla})=\begin{pmatrix}0&(-2i\partial_{z})^{m}\\\
(-2i\overline{\partial_{z}})^{m}&0\end{pmatrix},$ (S13)

where $\partial_{z}=(\partial_{x}-i\partial_{y})/2$. Note that in this real
space representation the Hamiltonian $\mathcal{H}_{m}(\mathbf{r})$ satisfies
the following symmetry relations:

$\begin{split}\rho(\mathcal{C}_{nz})\mathcal{H}_{m}(\mathbf{r})\rho(\mathcal{C}_{nz})^{\dagger}&=\mathcal{H}_{m}(\mathcal{C}_{nz}\mathbf{r}),\\\
\rho(\mathcal{A})\mathcal{H}_{m}^{*}(\mathbf{r})\rho(\mathcal{A})^{\dagger}&=\mathcal{H}_{m}(\mathcal{A}\mathbf{r}),\\\
\rho(\mathcal{S})\mathcal{H}_{m}(\mathbf{r})\rho(\mathcal{S})^{\dagger}&=-\mathcal{H}_{m}(\mathbf{r}),\\\
\end{split}$ (S14)

where $\mathcal{A}\mathbf{r}=\pm\mathbf{r}$ for $\mathcal{A}=\mathcal{T}$ and
$\mathcal{A}=\mathcal{C}_{2z}\mathcal{T}$, respectively. If there are mirror
or two fold rotation symmetries:

$\begin{split}\rho(\mathcal{C}_{2z})\mathcal{H}_{m}(\mathbf{r})\rho(\mathcal{C}_{2z})^{\dagger}&=\mathcal{H}_{m}(\mathcal{C}_{2z}\mathbf{r}),\\\
\rho(\mathcal{M})\mathcal{H}_{m}(\mathbf{r})\rho(\mathcal{M})^{\dagger}&=\mathcal{H}_{m}(\mathcal{M}\mathbf{r}).\end{split}$
(S15)

###  S-1.2 Adding moiré potential and Brillouin zone folding

Adding a moiré potential whose periodicity (let the moiré lattice vectors be
$\mathbf{a}_{1}^{m}$ and $\mathbf{a}_{2}^{m}$) is much larger than the
periodicity of the lattice (let the lattice vectors be $\mathbf{a}_{1}$ and
$\mathbf{a}_{2}$) makes the moiré Brillouin zone (mBZ) (with moiré reciprocal
lattice vectors $\mathbf{b}_{1}^{m}$ and $\mathbf{b}_{2}^{m}$) much smaller
than the BZ (with reciprocal lattice vectors $\mathbf{b}_{1}$ and
$\mathbf{b}_{2}$) of the lattice. For simplicity, here we are going to
consider the moiré periodicity to be commensurate with the lattice
periodicity. The point $\mathbf{k}_{0}$ in the lattice BZ, gets mapped to some
point
$\mathbf{k}_{0}^{m}=\mathbf{k}_{0}-m_{1}\mathbf{b}_{1}^{m}-m_{2}\mathbf{b}_{2}^{m}$.
We know that under translation by lattice vectors
$n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}$, the basis states transform as
$T_{n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle=e^{-i\mathbf{k}_{0}\cdot(n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2})}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle$.
Under translations by moiré lattice vectors
$n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m}$, they transform as

$\begin{split}T_{n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m}}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle&=e^{-i\mathbf{k}_{0}\cdot(n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m})}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle\\\
&=e^{-i(\mathbf{k}_{0}-m_{1}\mathbf{b}_{1}^{m}-m_{2}\mathbf{b}_{2}^{m})\cdot(n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m})}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle
e^{-i(m_{1}\mathbf{b}_{1}^{m}+m_{2}\mathbf{b}_{2}^{m})\cdot(n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m})}\\\
&=e^{-i\mathbf{k}_{0}^{m}\cdot(n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m})}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle,\end{split}$
(S16)

where in the first equality we used the fact that moiré lattice vectors are
integer multiples of lattice vectors, and in last equality we used
$\mathbf{b}_{i}^{m}.\mathbf{a}_{j}^{m}=2\pi\delta_{ij}$. Therefore, under
moiré translations they transform as basis states at momentum
$\mathbf{k}_{0}^{m}$ in the mBZ.

###  S-1.3 Single layer system with QBCP under periodic potential

If we add a moiré potential to single layer QBCP Hamiltonian that does not
break $\mathcal{C}_{nz}$, $\mathcal{A}=\mathcal{T}$ and $\mathcal{S}$, the
general form of the Hamiltonian becomes:

$\mathcal{H}_{S}(\mathbf{r})=\begin{pmatrix}0&\mathcal{D}_{k}^{\dagger}(\mathbf{r})+\mathcal{D}_{U}^{*}(\mathbf{r};\bm{\alpha})\\\
\mathcal{D}_{k}(\mathbf{r})+\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})&0\end{pmatrix}=\begin{pmatrix}0&(-2i\partial_{z})^{m}+\mathcal{D}_{U}^{*}(\mathbf{r};\bm{\alpha})\\\
(-2i\overline{\partial_{z}})^{m}+\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})&0\end{pmatrix}.$
(S17)

Under moiré periodicities $\mathbf{a}_{1}^{m}$ and $\mathbf{a}_{2}^{m}$, the
basis states $|\mathbf{k}_{0},\mathbf{r},\alpha\rangle$ transform as
$T_{\mathbf{a}_{i}^{m}}|\mathbf{k}_{0},\mathbf{r},\alpha\rangle=e^{-i\mathbf{k}_{0}^{m}\cdot\mathbf{a}_{i}^{m}}|\mathbf{k}_{0},\mathbf{r}+\mathbf{a}_{i}^{m},\alpha\rangle$.
Hence, the representation of the moiré translations are
$\rho(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=e^{-i\mathbf{k}_{0}^{m}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\mathds{1}$,
and
$\mathcal{D}_{U}(\mathbf{r}+\mathbf{a}_{i}^{m};\bm{\alpha})=\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})$.
Due to Eq. (S14), $\mathcal{D}_{U}(\mathbf{r};\bm{\alpha})$ satisfies

$\begin{split}\mathcal{D}_{U}(\mathcal{C}_{nz}\mathbf{r};\bm{\alpha})=(\omega^{*})^{2}\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}),\\\
\end{split}$ (S18)

and if additionally there are mirror and rotation symmetries in Eq. (S15),
then

$\begin{split}\mathcal{D}_{U}(\mathcal{C}_{2z}\mathbf{r};\bm{\alpha})=\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}),\\\
\mathcal{D}_{U}(\mathcal{M}\mathbf{r};\bm{\alpha})=\mathcal{D}_{U}^{*}(\mathbf{r};\bm{\alpha}).\\\
\end{split}$ (S19)

A composition of chiral and time reversal symmetry in these system gives a
anti-unitary particle-hole symmetry $\mathcal{P}=\mathcal{S}\mathcal{T}$, its
representation is
$\rho(\mathcal{P})=\rho(\mathcal{S})\rho(\mathcal{T})=\sigma_{y}$, hence
$\mathcal{P}^{2}=-\mathds{1}$ since
$\rho(\mathcal{P})\rho(\mathcal{P})^{*}=\sigma_{y}\sigma_{y}^{*}=-\mathds{1}$.

###  S-1.4 Twisted bilayer systems with DP or QBCP

In twisted bilayer systems, two layers are twisted from each other by small
angle $\theta$. That creates an offset between the HSMs of the two layers
$\mathbf{k}_{0}^{(1)}=\mathbf{R}(-\theta/2)\mathbf{k}_{0}$ and
$\mathbf{k}_{0}^{(2)}=\mathbf{R}(\theta/2)\mathbf{k}_{0}$:
$\mathbf{k}_{0}^{(1)}-\mathbf{k}_{0}^{(2)}=\mathbf{q}_{1}\approx\mathbf{k}_{0}\times\theta\hat{z}$
(here the numbers in the superscripts denote the layer index). In case of TBG,
$\mathbf{q}_{1}=-(\mathbf{b}_{1}^{m}+\mathbf{b}_{2}^{m})/3$ in the notation
specified in Fig. 2 of the main text. In case of twisted bilayer checkerboard
lattice, $\mathbf{q}_{1}=(\mathbf{b}_{1}^{m}-\mathbf{b}_{2}^{m})/2$ in the
notation specified in Fig. S8. The general form of the twisted bilayer system
in the chiral limit is

$\hat{H}_{TB}=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}\mathcal{H}_{TB}(\mathbf{r})\\{\langle\mathbf{k}_{0}^{(1)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(1)},\mathbf{r},B|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},B|\\}^{T},$
(S20)

with

$\mathcal{H}_{TB}(\mathbf{r})=\left(\begin{array}[]{c|c}\mathcal{H}_{m}(\mathbf{R}(-\theta/2)\mathbf{r})&\begin{array}[]{cc}0&U_{2}(\mathbf{r})\\\
U_{1}(\mathbf{r})&0\end{array}\\\
\hline\cr\begin{array}[]{cc}0&U_{1}^{*}(\mathbf{r})\\\
U_{2}^{*}(\mathbf{r})&0\end{array}&\mathcal{H}_{m}(\mathbf{R}(\theta/2)\mathbf{r})\end{array}\right),$
(S21)

where the zeros in the off-diagonal are due to chiral symmetry, which has
representation $\rho(\mathcal{S})=\mathds{1}\otimes\sigma_{z}$ in this basis.
Note that
$\mathcal{H}_{m}(\mathbf{R}(-\theta/2)\mathbf{r})=e^{im\theta\sigma_{z}/4}\mathcal{H}_{m}(\mathbf{r})e^{-im\theta\sigma_{z}/4}$.
Then transforming
$H_{TB}(\mathbf{r})\rightarrow\text{Diag}\\{e^{-im\theta\sigma_{z}/4}\mathchar
44\relax\penalty
0e^{im\theta\sigma_{z}/4}\\}H_{TB}(\mathbf{r})\text{Diag}\\{e^{im\theta\sigma_{z}/4}\mathchar
44\relax\penalty 0e^{-im\theta\sigma_{z}/4}\\}$, we obtain

$\mathcal{H}_{TB}(\mathbf{r})=\left(\begin{array}[]{c|c}\mathcal{H}_{m}(\mathbf{r})&\begin{array}[]{cc}0&U_{2}(\mathbf{r})\\\
U_{1}(\mathbf{r})&0\end{array}\\\
\hline\cr\begin{array}[]{cc}0&U_{1}^{*}(\mathbf{r})\\\
U_{2}^{*}(\mathbf{r})&0\end{array}&\mathcal{H}_{m}(\mathbf{r})\end{array}\right).$
(S22)

Note that the last transformation does not change the off-diagonal part
because $\\{\sigma_{z},\sigma_{x/y}\\}=0$. Furthermore, the last
transformation changes the basis; however, it can be seen straightforwardly
that the representations of the relevant symmetries remain the same. It is
customary in the literature to change the order of the basis states to write
the Hamiltonian in the following form:

$\hat{H}_{TB}=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}\mathcal{H}_{TB}(\mathbf{r})\\{\langle\mathbf{k}_{0}^{(1)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(1)},\mathbf{r},B|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},B|\\}^{T},$
(S23)

with

$\mathcal{H}_{TB}(\mathbf{r})=\begin{pmatrix}0&0&(-2i\partial_{z})^{m}&U_{2}(\mathbf{r})\\\
0&0&U_{1}^{*}(\mathbf{r})&(-2i\partial_{z})^{m}\\\
(-2i\overline{\partial_{z}})^{m}&U_{1}(\mathbf{r})&0&0\\\
U_{2}^{*}(\mathbf{r})&(-2i\overline{\partial_{z}})^{m}\ &0&0\ \end{pmatrix}$
(S24)

In this basis the representation of the symmetries are
$\rho(\mathcal{C}_{nz})=\text{Diag}\\{\omega,\omega,\omega^{*},\omega^{*}\\}$,
$\rho(\mathcal{A})=\sigma_{x}\otimes\mathds{1}$,
$\rho(\mathcal{S})=\sigma_{z}\otimes\mathds{1}$. The symmetries
$\mathcal{C}_{nz}$ and $\mathcal{A}$ put constraints on $U_{1}(\mathbf{r})$
and $U_{2}(\mathbf{r})$:

$\begin{split}U_{1}(\mathcal{A}\mathbf{r})&=U_{2}^{*}(\mathbf{r})\\\
U_{1}(\mathcal{C}_{nz}\mathbf{r})&=(\omega^{*})^{2}U_{1}(\mathbf{r})\end{split}$
(S25)

There are further constraints on $U_{i}(\mathbf{r})$ from moiré translations.
The basis states for the $\mathbf{k}\cdot\mathbf{p}$ theory of each layer
satisfy:
$T_{\mathbf{a}_{i}}|\mathbf{k}_{0}^{(l)},\mathbf{r},\alpha\rangle=e^{-i\mathbf{k}_{0}^{(l)}\cdot\mathbf{a}_{i}}|\mathbf{k}_{0}^{(l)},\mathbf{r}+\mathbf{a}_{i},\alpha\rangle$.
Then, under a moiré transaltion
$\mathbf{t}^{m}=\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}=n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m}$,
the basis states transform as

$\begin{split}&T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}}\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}\\\
=&\\{|\mathbf{k}_{0}^{(1)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},B\rangle\\}\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})\end{split}$
(S26)

where

$\begin{split}\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})&=\text{Diag}\\{e^{-i\mathbf{k}_{0}^{(1)}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})},e^{-i\mathbf{k}_{0}^{(2)}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})},e^{-i\mathbf{k}_{0}^{(1)}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})},e^{-i\mathbf{k}_{0}^{(2)}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\\}\\\
&=e^{-i\mathbf{k}_{0}^{(1)}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\text{Diag}\\{1,e^{i(\mathbf{k}_{0}^{(1)}-\mathbf{k}_{0}^{(2)})\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})},1,e^{i(\mathbf{k}_{0}^{(1)}-\mathbf{k}_{0}^{(2)})\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\\}\\\
&=e^{-i\mathbf{k}_{0}^{(1)}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\text{Diag}\\{1,e^{i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})},1,e^{i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\\}\\\
&=e^{-i\mathbf{k}_{0}^{m}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\text{Diag}\\{1,e^{i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})},1,e^{i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\\}\end{split}$
(S27)

where we define $\mathbf{k}_{0}^{m}$ as the mapping of $\mathbf{k}_{0}^{(1)}$
in the mBZ after the folding of the BZ. Then, the Hamiltonian of the twisted
bilayer transforms as

$\begin{split}&T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}}\hat{H}_{TB}T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}}^{-1}\\\
=&\int
d^{2}\mathbf{r}\,T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}}\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}\mathcal{H}_{TB}(\mathbf{r})\times\\\
&\hskip
85.35826pt\\{\langle\mathbf{k}_{0}^{(1)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(1)},\mathbf{r},B|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},B|\\}^{T}T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}}^{-1}\\\
=&\int
d^{2}\mathbf{r}\,\\{|\mathbf{k}_{0}^{(1)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},B\rangle\\}\times\\\
&\hskip
85.35826pt\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})\mathcal{H}_{TB}(\mathbf{r})\rho^{\dagger}(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})\times\\\
&\hskip
85.35826pt\\{\langle\mathbf{k}_{0}^{(1)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},A|,\langle\mathbf{k}_{0}^{(1)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},B|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},A|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m},B|\\}^{T}\\\
=&\int
d^{2}\mathbf{r}\,\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}\times\\\
&\hskip
85.35826pt\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})\mathcal{H}_{TB}(\mathbf{r}-\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})\rho^{\dagger}(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})\\{\langle\mathbf{k}_{0}^{(1)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(1)},\mathbf{r},B|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},A|,\langle\mathbf{k}_{0}^{(2)},\mathbf{r},B|\\}^{T}\\\
&=\hat{H}_{TB}\\\
\Rightarrow&\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})\mathcal{H}_{TB}(\mathbf{r})\rho^{\dagger}(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})=\mathcal{H}_{TB}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}).\end{split}$
(S28)

This puts the following constraints on the moiré potentials

$\begin{split}U_{1}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=e^{-i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}U_{1}(\mathbf{r}),\\\
U_{2}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=e^{-i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}U_{2}(\mathbf{r}),\end{split}$
(S29)

###  S-1.5 Symmetries of chiral TBG

In this section we specialize to TBG and summarize what we obtained in the
previous subsection as well as point out some extra symmetries of TBG. The
fact $U_{2}^{*}(\mathbf{r})=U_{1}(-\mathbf{r})$ from Eq. (S25) (since for TBG
$\mathcal{A}=\mathcal{C}_{2z}\mathcal{T}$, we have
$\mathcal{A}\mathbf{r}=-\mathbf{r}$) gives

$\mathcal{H}_{TBG}(\mathbf{r})=\begin{pmatrix}0&0&(-2i\partial_{z})^{m}&U_{1}^{*}(-\mathbf{r})\\\
0&0&U_{1}^{*}(\mathbf{r})&(-2i\partial_{z})^{m}\\\
(-2i\overline{\partial_{z}})^{m}&U_{1}(\mathbf{r})&0&0\\\
U_{1}(-\mathbf{r})&(-2i\overline{\partial_{z}})^{m}\ &0&0\ \end{pmatrix}.$
(S30)

From Eq. (S29) we obtain

$U_{1}(\mathbf{r})=\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}},$
(S31)

where the sum goes over all moiré reciprocal lattice vectors, and
$a_{\mathbf{b}^{m}}$ are Fourier amplitudes. Also from Eq. (S25), noting that
in case of TBG $\omega=e^{2\pi i/3}$, we get

$\begin{split}U_{1}(\mathcal{C}_{3z}\mathbf{r})=e^{2\pi
i/3}U_{1}(\mathbf{r})&\Rightarrow\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathcal{C}_{3z}\mathbf{r}}=e^{2\pi
i/3}\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}}\\\
&\Rightarrow\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathcal{C}_{3z}^{-1}\mathbf{q}_{1}+\mathcal{C}_{3z}^{-1}\mathbf{b}^{m})\cdot\mathbf{r}}=e^{2\pi
i/3}\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}}\\\
&\Rightarrow\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}_{2}^{m}+\mathcal{C}_{3z}^{-1}\mathbf{b}^{m})\cdot\mathbf{r}}=e^{2\pi
i/3}\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}}\\\
&\Rightarrow\sum_{\mathbf{b}^{m}}a_{\mathcal{C}_{3z}(\mathbf{b}^{m}-\mathbf{b}^{m}_{2})}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}}=e^{2\pi
i/3}\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}},\\\
&\Rightarrow a_{\mathcal{C}_{3z}(\mathbf{b}^{m}-\mathbf{b}^{m}_{2})}=e^{2\pi
i/3}a_{\mathbf{b}^{m}}\text{ for all reciprocal lattice vector
}\mathbf{b}^{m}.\end{split}$ (S32)

Eq. (S30) along with Eqs. (S31) and (S32) is the most general form of the TBG
Hamiltonian. Using Eq. (S32), starting from $a_{\mathbf{0}}\equiv\alpha$, we
get

$\begin{split}&a_{-\mathcal{C}_{3z}\mathbf{b}^{m}_{2}}=a_{\mathbf{b}^{m}_{1}}=e^{2\pi
i/3}a_{\mathbf{0}}=e^{2\pi i/3}\alpha,\\\
&a_{\mathcal{C}_{3z}(\mathbf{b}^{m}_{1}-\mathbf{b}^{m}_{2})}=a_{\mathbf{b}^{m}_{2}}=e^{2\pi
i/3}a_{\mathbf{b}^{m}_{1}}=e^{4\pi/3}a_{\mathbf{0}}=e^{4\pi
i/3}\alpha.\end{split}$ (S33)

Hence, if we keep only lowest harmonics the expression for $U_{1}(\mathbf{r})$
becomes

$\begin{split}U_{1}(\mathbf{r})&=a_{\mathbf{0}}(e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{2\pi
i/3}e^{-i(\mathbf{q}_{1}+\mathbf{b}_{1}^{m})\cdot\mathbf{r}}+e^{4\pi
i/3}e^{-i(\mathbf{q}_{1}-\mathbf{b}_{2}^{m})\cdot\mathbf{r}})\\\
&=\alpha(e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{2\pi
i/3}e^{-i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{4\pi
i/3}e^{-i\mathbf{q}_{3}\cdot\mathbf{r}}).\end{split}$ (S34)

This form of $U_{1}(\mathbf{r})$ was used in tarnopolsky2019origins (and in
Eq. 3 in the main text). Similarly, if one starts from
$a_{\mathbf{\mathbf{b}_{1}^{m}+\mathbf{b}_{2}^{m}}}$, one gets relations to
$a_{\mathbf{\mathbf{b}_{2}^{m}-\mathbf{b}_{1}^{m}}}$ and
$a_{\mathbf{\mathbf{b}_{1}^{m}-\mathbf{b}_{2}^{m}}}$. Furthermore, if one sets
$a_{\mathbf{\mathbf{b}_{1}^{m}+\mathbf{b}_{2}^{m}}}=-a_{\mathbf{0}}=-\alpha$,
in that fine-tuned case, $U_{1}(\mathbf{r})$ becomes

$U_{1}(\mathbf{r})=\alpha\left[(e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{2\pi
i/3}e^{-i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{4\pi
i/3}e^{-i\mathbf{q}_{3}\cdot\mathbf{r}})-(e^{2i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{2\pi
i/3}e^{2i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{4\pi
i/3}e^{2i\mathbf{q}_{3}\cdot\mathbf{r}})\right],$ (S35)

this form of $U_{1}(\mathbf{r})$ was used in becker2023degenerates to obtain 4
exact flat bands for some “magic” values of $\alpha$.

Chiral TBG also has the following symmetries:

  * •

$\mathcal{C}_{2x}$ symmetry that swaps
$|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle\leftrightarrow|\mathbf{k}_{0}^{(2)},\mathcal{M}_{y}\mathbf{r},B\rangle$
and
$|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle\leftrightarrow|\mathbf{k}_{0}^{(2)},\mathcal{M}_{y}\mathbf{r},A\rangle$,
where $\mathcal{M}_{y}\mathbf{r}=\mathcal{M}_{y}\\{x,y\\}=\\{x,-y\\}$. This
means the following

$\mathcal{C}_{2x}\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}=\\{|\mathbf{k}_{0}^{(1)},\mathcal{M}_{y}\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathcal{M}_{y}\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathcal{M}_{y}\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathcal{M}_{y}\mathbf{r},B\rangle\\}\rho(\mathcal{C}_{2x})\\\
$ (S36)

where $\rho(\mathcal{C}_{2x})=\sigma_{x}\otimes\sigma_{x}$. This also implies

$\rho(\mathcal{C}_{2x})\mathcal{H}_{TBG}(\mathbf{r})\rho^{\dagger}(\mathcal{C}_{2x})=\mathcal{H}_{TBG}(\mathcal{M}_{y}\mathbf{r})\Rightarrow
U_{1}^{*}(\mathbf{r})=U_{1}(\mathcal{M}_{y}\mathbf{r}).$ (S37)

This and the facts that
$\mathbf{q}\cdot\mathcal{M}_{y}\mathbf{r}=(\mathcal{M}_{y}\mathbf{q})\cdot\mathbf{r}$,
$\mathcal{M}_{y}\mathbf{q}_{1}=-\mathbf{q}_{1}$,
$\mathcal{M}_{y}\mathbf{q}_{2}=-\mathbf{q}_{3}$ and
$\mathbf{M}_{y}\mathbf{q}_{3}=-\mathbf{q}_{2}$ imply $\alpha^{*}=\alpha$ in
Eqs. (S34) and (S35). Applying a magnetic field normal to TBG plane breaks
$\mathcal{C}_{2x}$, this is why in le2022doubles the parameter $\alpha$ was a
complex number under alternating magnetic field.

  * •

Unitary particle-hole symmetry $P$ whose representation is
$\rho(P)=\mathds{1}\otimes i\sigma_{y}$ with $P\mathbf{r}=-\mathbf{r}$ and

$P\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}=\\{|\mathbf{k}_{0}^{(1)},-\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},-\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},-\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},-\mathbf{r},B\rangle\\}\rho(P)\\\
$ (S38)

such that
$\rho(P)\mathcal{H}_{TBG}(\mathbf{r})=-\mathcal{H}_{TBG}(-\mathbf{r})\rho(P)$.
Note that $P\mathbf{k}=-\mathbf{k}$, but it does not change the valley.

  * •

Composition of $P$ and $\mathcal{S}$ acts like a $\mathcal{C}_{2z}$ symmetry:

    1. 1.

since both $P$ and $\mathcal{S}$ anticommute with $\hat{H}_{TBG}$,
$P\mathcal{S}$ commutes with $\hat{H}_{TBG}$,

    2. 2.

since $P\mathbf{r}=-\mathbf{r}$, $P\mathbf{k}=-\mathbf{k}$, and $\mathcal{S}$
is local in both position and momentum space,
$P\mathcal{S}\mathbf{r}=-\mathbf{r}$ and $P\mathcal{S}\mathbf{k}=-\mathbf{k}$,

    3. 3.

$\rho(P\mathcal{S})$ commutes with $\rho(\mathcal{C}_{3z})$,
$\rho(\mathcal{C}_{2x})$ and $\rho(\mathcal{C}_{2z}\mathcal{T})$.

This has representation $\rho(P\mathcal{S})=\sigma_{z}\otimes i\sigma_{y}$.

  * •

A composition of $P\mathcal{S}$ with $\mathcal{C}_{2x}$ gives a
$\mathcal{C}_{2y}$ type symmetry with representation
$\rho(P\mathcal{SC}_{2x})=i\sigma_{y}\otimes\sigma_{z}$.

  * •

A composition of $P$ with $\mathcal{A}=\mathcal{C}_{2z}\mathcal{T}$ gives an
antiunitary particle-hole symmetry $\mathcal{P}=P\mathcal{C}_{2z}\mathcal{T}$,
which squares to
$\rho(P\mathcal{C}_{2z}\mathcal{T})\rho^{*}(P\mathcal{C}_{2z}\mathcal{T})=(\sigma_{x}\otimes
i\sigma_{y})^{2}=-\mathds{1}$.

  * •

Since $P\mathcal{S}$ is a $\mathcal{C}_{2z}$ type symmetry, a composition of
$P\mathcal{S}$ and $\mathcal{C}_{2z}\mathcal{T}$ gives a time-reversal like
symmetry with representation
$\rho(P\mathcal{S}\mathcal{C}_{2z}\mathcal{T})=i\sigma_{y}\otimes\sigma_{z}$.

  * •

The composition $P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}$
commutes with the Hamiltonian, its representation is diagonal
$\rho(P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T})=\sigma_{z}\otimes\sigma_{z}$,
and it acts on position and momentum as
$P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}\mathbf{r}=P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}\\{x,y\\}=\\{x,-y\\}=\mathcal{M}_{y}\mathbf{r}$
whereas
$P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}\mathbf{k}=P\mathcal{S}\mathcal{C}_{2x}\mathcal{C}_{2z}\mathcal{T}\\{x,y\\}=\\{-k_{x},k_{y}\\}=\mathcal{M}_{x}\mathbf{k}$.
Therefore, it has an effect on position and momentum like
$\mathcal{C}_{2x}\mathcal{T}$.

Because of the $\mathcal{C}_{2z}$ type symmetry $P\mathcal{S}$, the space
group is enhanced in single valley chiral TBG to $p622$ from
$p6^{\prime}2^{\prime}2$ for single-valley non-chiral TBG.

##  S-2 Properties of zero energy eigen wave-function
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$

In the previous section we discussed the moiré Hamiltonians for single and
bilayer DP or QBCP and their symmetry properties. In this section we discuss
the zero energy sublattice polarized eigen wave-functions
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ of the moiré Hamiltonian, how they
transform under the symmetries and how symmetries dictate appearance of zeros
of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$.

###  S-2.1 Zero energy eigen wave-function
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ and its transformation under
symmetries

First, to treat single and bilayer system together, we introduce shorthand
$\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times 2l}$ for the basis states,
where $\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times
2}=\\{|\mathbf{k}_{0},\mathbf{r},A\rangle,|\mathbf{k}_{0},\mathbf{r},B\rangle\\}$
and $\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times
4}=\\{|\mathbf{k}_{0}^{(1)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},A\rangle,|\mathbf{k}_{0}^{(1)},\mathbf{r},B\rangle,|\mathbf{k}_{0}^{(2)},\mathbf{r},B\rangle\\}$
for single and bilayer systems respectively. In this notation, the Hamiltonian
is written as

$\hat{H}_{S/TB}=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times
2l}\mathcal{H}_{S/TB}(\mathbf{r})\\{\langle\mathbf{k}_{0},\mathbf{r}|\\}^{T}_{2l\times
1}.$ (S39)

Before adding any moiré potential, there are two zero energy eigenmodes at
$\mathbf{k}_{0}$ protected by $\mathcal{C}_{nz}$ and $\mathcal{A}$ for both DP
and QBCP. They have $\mathcal{C}_{nz}$ eigenvalues $\omega$ and $\omega^{*}$.
Now, if the moiré potential does not break $\mathcal{C}_{nz}$ and
$\mathcal{A}$, the degeneracy of the DP and QBCP remain. If the moiré
potential does not break $\mathcal{S}$, then the DP or QBCP has to remain at
zero energy since $[\mathcal{C}_{nz},\mathcal{S}]=0$. This means there are two
states $|\Psi_{\mathbf{k}_{0}^{m},1}\rangle$ and
$|\Psi_{\mathbf{k}_{0}^{m},2}\rangle$ such that
$\hat{H}_{S/TB}|\Psi_{\mathbf{k}_{0}^{m},i}\rangle=0$ for any moiré potential
that keeps $\mathcal{C}_{nz}$, $\mathcal{A}$ and $\mathcal{S}$ intact (recall
that $\mathbf{k}_{0}^{m}$ is the momentum in the mBZ where $\mathbf{k}_{0}$
gets mapped to after BZ folding). Furthermore, the irrep of the degeneracy at
$\mathbf{k}_{0}^{m}$ at zero energy before and after adding the moiré
potential is the same. These two states can be chosen to be sublattice
polarized with following properties:

$\begin{split}|\Psi_{\mathbf{k}_{0}^{m},1}\rangle&=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times
2l}\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times 1}\\\
\mathbf{0}_{l\times 1}\end{Bmatrix},\\\
|\Psi_{\mathbf{k}_{0}^{m},2}\rangle&=\mathcal{A}|\Psi_{\mathbf{k}_{0},1}\rangle=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times
2l}\begin{Bmatrix}\mathbf{0}_{l\times 1}\\\
\psi_{\mathbf{k}_{0}^{m}}^{*}(\mathcal{A}\mathbf{r})_{l\times
1}\end{Bmatrix}.\end{split}$ (S40)

Below we discuss the effect of symmetries on the wave-function
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$.

  1. 1.

Due to DP or QBCP protected by $\mathcal{C}_{nz}$, we have

$\begin{split}\mathcal{C}_{nz}|\Psi_{\mathbf{k}_{0}^{m},1}\rangle&=\omega|\Psi_{\mathbf{k}_{0}^{m},1}\rangle\\\
\mathcal{C}_{nz}|\Psi_{\mathbf{k}_{0}^{m},2}\rangle&=\omega^{*}|\Psi_{\mathbf{k}_{0}^{m},2}\rangle.\end{split}$
(S41)

The last two relations mean

$\begin{split}\omega|\Psi_{\mathbf{k}_{0}^{m},1}\rangle=\mathcal{C}_{nz}|\Psi_{\mathbf{k}_{0}^{m},1}\rangle&=\int
d^{2}\mathbf{r}\mathcal{C}_{nz}\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times
2l}\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times 1}\\\
\mathbf{0}_{l\times 1}\end{Bmatrix}\\\ &=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathcal{C}_{nz}\mathbf{r}\rangle\\}_{1\times
2l}\rho(\mathcal{C}_{nz})\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times
1}\\\ \mathbf{0}_{l\times 1}\end{Bmatrix}\\\ &=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathcal{C}_{nz}\mathbf{r}\rangle\\}_{1\times
2l}{\color[rgb]{0,0,0}\text{Diag}\\{\omega,\omega^{*}\\}\otimes\mathds{1}_{l\times
l}}\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times 1}\\\
\mathbf{0}_{l\times 1}\end{Bmatrix}\\\ &=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathcal{C}_{nz}\mathbf{r}\rangle\\}_{1\times
2l}\omega\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times 1}\\\
\mathbf{0}_{l\times 1}\end{Bmatrix}\\\ &=\int
d^{2}\mathbf{r}\\{|\mathbf{k}_{0},\mathbf{r}\rangle\\}_{1\times
2l}\omega\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathcal{C}_{nz}^{-1}\mathbf{r})_{l\times
1}\\\ \mathbf{0}_{l\times 1}\end{Bmatrix}\\\
\Rightarrow\omega\psi_{\mathbf{k}_{0}^{m}}(\mathcal{C}_{nz}^{-1}\mathbf{r})_{l\times
1}&=\omega\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times
1}\Rightarrow\psi_{\mathbf{k}_{0}^{m}}(\mathcal{C}_{nz}\mathbf{r})=\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}).\\\
\end{split}$ (S42)

  2. 2.

If there is a mirror $\mathcal{M}$ (or $\mathcal{C}_{2y}$ like symmetry
$P\mathcal{S}\mathcal{C}_{2x}$ in case of chiral TBG) at $\mathbf{k}_{0}$, we
can choose its sewing matrix
$[B_{\mathbf{k}_{0}^{m}}(\mathcal{M})]_{ij}=\langle\Psi_{\mathbf{k}_{0}^{m},i}|\rho(\mathcal{M})|\Psi_{\mathbf{k}_{0}^{m},j}\rangle=\sigma_{x}$
(see, for example, Bilbao Crystallography server
aroyo2011crystallographys,aroyo2006bilbaoIs,aroyo2006bilbaoIIs). This implies

$\mathcal{M}|\Psi_{\mathbf{k}_{0}^{m},1}\rangle=|\Psi_{\mathbf{k}_{0}^{m},2}\rangle\Rightarrow\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times
1}=\psi_{\mathbf{k}_{0}^{m}}^{*}(\mathcal{MA}\mathbf{r})_{l\times 1}.$ (S43)

This implies that along $\mathcal{MA}$ invariant line (given by
$\mathcal{MA}\mathbf{r}=\mathbf{r}$), the wave-function
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times 1}$ is a real function.

  3. 3.

If there is $\mathcal{C}_{2z}$ symmetry such that
$\mathcal{C}_{2z}\mathbf{k}_{0}=\mathbf{k}_{0}$ (this is true for QBCP in $4$
and $6$-fold rotation symmetry), its representation $\rho(\mathcal{C}_{2z})$
is $\zeta\mathds{1}$ with $\zeta=\pm 1$ (see Sec. S-1.1), and the sewing
matrix
$[B_{\mathbf{k}_{0}^{m}}(\mathcal{C}_{2z})]_{ij}=\langle\Psi_{\mathbf{k}_{0}^{m},i}|\rho(\mathcal{C}_{2z})|\Psi_{\mathbf{k}_{0}^{m},j}\rangle=\zeta\delta_{ij}$
since the irrep of the two degenerate zero energy modes remains the same
before and after adding moiré potential as discussed earlier. From this we
obtain:

$\begin{split}&\mathcal{C}_{2z}|\Psi_{\mathbf{k}_{0}^{m},1}\rangle=\zeta|\Psi_{\mathbf{k}_{0}^{m},1}\rangle\\\
\Rightarrow&\psi_{\mathbf{k}_{0}^{m}}(\mathcal{C}_{2z}\mathbf{r})_{l\times
1}=\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times 1}.\end{split}$ (S44)

  4. 4.

Under moiré lattice translation
$n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m}$, the states
$|\Psi_{\mathbf{k}_{0}^{m},i}\rangle$ should transform as a Bloch states at
$\mathbf{k}_{0}^{m}$, hence

$\begin{split}&T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}}|\Psi_{\mathbf{k}_{0}^{m},i}\rangle=e^{-i\mathbf{k}_{0}^{m}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}|\Psi_{\mathbf{k}_{0}^{m},i}\rangle\\\
\Rightarrow&\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})_{l\times
1}\\\ \mathbf{0}_{l\times
1}\end{Bmatrix}=e^{-i\mathbf{k}_{0}^{m}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\begin{Bmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})_{l\times
1}\\\ \mathbf{0}_{l\times 1}\end{Bmatrix}.\end{split}$ (S45)

The consequences of this for single and twisted bilayer cases are as follows:

     * •

For single layer QBCP systems, as discussed earlier,
$\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})=e^{-i\mathbf{k}_{0}^{m}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\mathds{1}$.
This implies
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$,
or in other words $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ is periodic function
with periodicity of moiré lattice vector.

     * •

For twisted bilayer systems, as discussed earlier,
$\rho(T_{\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m}})=e^{-i\mathbf{k}_{0}^{m}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\text{Diag}\\{1\mathchar
44\relax\penalty
0e^{i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\mathchar
44\relax\penalty 01\mathchar 44\relax\penalty
0e^{i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\\}$. This
implies that the two components of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$
transform differently: The first component satisfies
$\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r})$,
hence it is a periodic function with moiré periodicity. The second component
satisfies
$\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=e^{i\mathbf{q}_{1}\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r})$.
In case of TBG, in the notation of Fig. 2, the phase factor is
$e^{-i(\mathbf{b}_{1}+\mathbf{b}_{2})/3\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}=e^{-2\pi
i(n_{1}+n_{2})/3}$. In case of TBCL, in the notation of Fig. S8, the phase
factor is
${\color[rgb]{0,0,0}e^{i(\mathbf{b}_{1}-\mathbf{b}_{2})/2\cdot(\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})}=e^{\pi
i(n_{1}-n_{2})}}$.

###  S-2.2 Zeros of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$, and constraints
from symmetries

As discussed in the main text, the general construction of exact flat-band
(FB) wave-function (WF) requires a zero in
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ somewhere in the unit cell since that
can cancel the pole of the holomorphic function
$f_{\mathbf{k}}(z;\mathbf{r}_{0})$ such the FB WF is finite everywhere in real
spalce. However, not any zero of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ is
suitable for construction of exact flat band wave-function, only the ones at
high symmetry points (HSPs) are suitable. Below, we present a more detailed
discussion on this for single and bilayer systems.

####  S-2.2.1 Single layer system with QBCP

For single layer systems, $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ is a one
component function. Near a zero at position $\mathbf{r}_{0}$ of
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$
($\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0})=0$), generically the function
behaves as
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+\delta(x,y))=\delta(ax+by)+\mathcal{O}(\delta^{2})$
for small $\delta$, where $a$ and $b$ are complex constants. Since the
holomorphic function $f_{\mathbf{k}}(z;\mathbf{r}_{0})$ has a pole at
$\mathbf{r}_{0}$, near $\mathbf{r}_{0}$ it behave as
$f_{\mathbf{k}}(\mathbf{r}_{0}+\delta(x,y);\mathbf{r}_{0})=\frac{c_{0}}{\delta(x+iy)}+\mathcal{O}(\delta^{0})$.
Hence the function
$\psi_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})=f_{\mathbf{k}}(z;\mathbf{r}_{0})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$
behaves as
$\psi_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+\delta(x,y))=\frac{c_{0}(ax+by)}{x+iy}+\mathcal{O}(\delta)$.
However, this function is singular at $\mathbf{r}_{0}$ for generic values of
$a$ and $b$ since the function oscillates around the point $\mathbf{r}_{0}$.
However, if $\mathbf{r}_{0}$ is a HSP meaning that there exists a rotation
symmetry $\mathcal{C}_{n^{\prime}z}$ ($n^{\prime}\geq 2$) and some moiré
lattice vector $\mathbf{a}^{m}$ such that
$\mathcal{C}_{n^{\prime}z}\mathbf{r}_{0}+\mathbf{a}^{m}=\mathbf{r}_{0}$, then
we know from Sec. S-2.1 that
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+\delta(x,y))=\psi_{\mathbf{k}_{0}^{m}}(\mathcal{C}_{n^{\prime}z}(\mathbf{r}_{0}+\delta(x,y)))=\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}-\mathbf{a}^{m}+\mathcal{C}_{n^{\prime}z}\delta(x,y))=\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+\mathcal{C}_{n^{\prime}z}\delta(x,y)))$.
This implies $a=b=0$. This would make the function
$\psi_{\mathbf{k}+\mathbf{k}_{0}}(\mathbf{r})$ regular.

####  S-2.2.2 Twisted bilayer graphene

In twisted bilayer graphene, $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ is a two
component function. At a zero of both real and imaginary parts of each
component of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ have to be zero. This
means that four equations need to be satisfied at once for a zero. However, it
turns out that due to the twisting between the two layers, one of the
components is automatically zero at some HSPs for any value of the tuning
parameter $\bm{\alpha}$. Then, only the other component of
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ need to be set to zero by tuning the
parameters $\bm{\alpha}$. We show this below.

As was shown in S-2.1, under $\mathcal{C}_{3z}$, each of the two components
satisfies
$\psi_{\mathbf{k}_{0}^{m},i}(\mathcal{C}_{3z}\mathbf{r})=\psi_{\mathbf{k}_{0}^{m},i}(\mathbf{r})$,
whereas under moiré lattice translations they satisfy
$\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r})$
and
$\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}+\sum_{i=1}^{2}n_{i}\mathbf{a}_{i}^{m})=e^{-2\pi
i(n_{1}+n_{2})/3}\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r})$. Near the corners
(BA and AB positions
($\mathbf{r}_{0}=\mp(\mathbf{a}_{2}^{m}-\mathbf{a}_{1}^{m})/3$ respectively,
see Fig. 2 of main text) of the moiré unit cell, we have

$\begin{split}\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+\mathbf{r})&=\psi_{\mathbf{k}_{0}^{m},1}(\mathcal{C}_{3z}\mathbf{r}_{0}+\mathcal{C}_{3z}\mathbf{r})=\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}\pm\mathbf{a}_{2}^{m}+\mathcal{C}_{3z}\mathbf{r})=\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+\mathcal{C}_{3z}\mathbf{r}),\\\
\Rightarrow\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+\mathcal{C}_{3z}\mathbf{r})&=\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+\mathbf{r}),\\\
\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+\mathbf{r})&=\psi_{\mathbf{k}_{0}^{m},2}(\mathcal{C}_{3z}\mathbf{r}_{0}+\mathcal{C}_{3z}\mathbf{r})=\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}\pm\mathbf{a}_{2}^{m}+\mathcal{C}_{3z}\mathbf{r})=e^{\mp
2\pi
i/3}\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+\mathcal{C}_{3z}\mathbf{r}),\\\
\Rightarrow\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+\mathcal{C}_{3z}\mathbf{r})&=e^{\pm
2\pi i/3}\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+\mathbf{r}),\end{split}$
(S46)

for BA and AB positions, respectively. The latter of the two equations implies
that $\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0})=0$ at BA and AB positions
for any $\bm{\alpha}$. The above equation also means the following:

  * •

BA: If for a “magic” value of $\bm{\alpha}$,
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0})=\mathbf{0}$, near the zero we have

$\begin{split}\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+\mathcal{C}_{3z}(\delta
x,\delta y))&=\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+(\delta x,\delta
y)),\\\ \Rightarrow\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+(\delta x,\delta
y))&=c_{1}\delta^{2}(x^{2}+y^{2})+\mathcal{O}(\delta^{3}),\\\
\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+\mathcal{C}_{3z}(\delta x,\delta
y))&=e^{2\pi i/3}\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+(\delta x,\delta
y)),\\\ \Rightarrow\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+(\delta x,\delta
y))&=c_{2}\delta(x+iy)+\mathcal{O}(\delta^{2}).\end{split}$ (S47)

These two together imply that
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+(\delta x,\delta
y))\sim\delta(x+iy)\times\text{(regular function)}$. Then
$f_{\mathbf{k}}(z;\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ is a
regular function near $\mathbf{r}_{0}$, and a FB WF construction is possible.

  * •

AB: The situation is more complicated here. If for a “magic” value of
$\bm{\alpha}$, $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0})=\mathbf{0}$, near
the zero we have generically

$\begin{split}\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+\mathcal{C}_{3z}(\delta
x,\delta y))&=\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+(\delta x,\delta
y)),\\\ \Rightarrow\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+(\delta x,\delta
y))&=c_{1}\delta^{2}(x^{2}+y^{2})+\mathcal{O}(\delta^{3}),\\\
\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+\mathcal{C}_{3z}(\delta x,\delta
y))&=e^{-2\pi i/3}\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+(\delta x,\delta
y)),\\\ \Rightarrow\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+(\delta x,\delta
y))&=c_{2}\delta(x-iy)+\mathcal{O}(\delta^{2}).\end{split}$ (S48)

If this were the case, the function
$f_{\mathbf{k}}(z;\mathbf{r}_{0})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})~{}\sim(x-iy)/(x+iy)$
would be singular, and a FB WF construction would not be possible. However,
the following result (taken from becker2023degenerates) comes to rescue.
Thoerem: If
$\mathcal{D}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=\mathbf{0}$
($\mathcal{D}(\mathbf{r})$ is defined in Eq. 3 of main text), and
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0})=\mathbf{0}$, then
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0}+(x,y))=(x+iy)g(\mathbf{r})$ where
$g(\mathbf{r})$ is a regular two component function. In other words, if
$\mathcal{D}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=\mathbf{0}$, and
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0})=\mathbf{0}$, then
$\overline{\partial_{z}}^{p}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})|_{\mathbf{r}=\mathbf{r}_{0}}=\mathbf{0}$
for any $p\in\mathds{N}$.
Proof: Notice that

$\mathcal{D}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=\mathbf{0}\Rightarrow(2i\overline{\partial_{z}})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=-\begin{pmatrix}0&U_{1}(\mathbf{r})\\\
U_{1}(-\mathbf{r})&0\end{pmatrix}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})\Rightarrow\overline{\partial_{z}}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=\frac{i}{2}\tilde{U}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}).$
(S49)

Then, clearly we have
$\overline{\partial_{z}}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})|_{\mathbf{r}=\mathbf{r}_{0}}=\frac{i}{2}\tilde{U}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}))|_{\mathbf{r}=\mathbf{r}_{0}}=\mathbf{0}$
since $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0})=\mathbf{0}$. Similarly,

$\begin{split}\left.\overline{\partial_{z}}^{2}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})\right|_{\mathbf{r}=\mathbf{r}_{0}}&=\left.\frac{i}{2}\overline{\partial_{z}}(\tilde{U}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}))\right|_{\mathbf{r}=\mathbf{r}_{0}}\\\
&=\left.\frac{i}{2}((\overline{\partial_{z}}\tilde{U}(\mathbf{r}))\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})+\tilde{U}(\mathbf{r})\overline{\partial_{z}}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}))\right|_{\mathbf{r}=\mathbf{r}_{0}}\\\
&=\left.\frac{i}{2}((\overline{\partial_{z}}\tilde{U}(\mathbf{r}))\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})+\frac{i}{2}\tilde{U}^{2}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}))\right|_{\mathbf{r}=\mathbf{r}_{0}}\\\
&=\mathbf{0}.\end{split}$ (S50)

Iterating this procedure of replacing
$\overline{\partial_{z}}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ with
$\frac{i}{2}\tilde{U}(\mathbf{r})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ and
using $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}_{0})=\mathbf{0}$, we see that
$\overline{\partial_{z}}^{p}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})|_{\mathbf{r}=\mathbf{r}_{0}}=\mathbf{0}$
for any $p\in\mathds{N}$.

In fact this theorem implies that, near a zero of
$\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r})$ is locally a holomorphic function.

This theorem together with transformation under $\mathcal{C}_{3z}$ implies
that if there is a zero at AB, the function
$\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ near AB behaves as

$\begin{split}\psi_{\mathbf{k}_{0}^{m},1}(\mathbf{r}_{0}+(\delta x,\delta
y))&=c_{1}\delta^{3}(x+iy)^{3}+\mathcal{O}(\delta^{4}),\\\
\psi_{\mathbf{k}_{0}^{m},2}(\mathbf{r}_{0}+(\delta x,\delta
y))&=c_{2}\delta^{2}(x+iy)^{2}+\mathcal{O}(\delta^{3}).\end{split}$ (S51)

This makes
$f_{\mathbf{k}}(z;\mathbf{r}_{0})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$ a
regular function near AB and, consequently, construction of FB WF possible.

##  S-3 Comparing analytically obtained wave-functions
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(i)}$ in Eq. 7 of main text with
numerical ones for single layer QBCP under periodic strain field with $p6mm$
symmetry

Figure S1: Comparison between analytically obtained wave-functions
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(i)}$ in Eq, 7 of main text (bottom
row) with numerical ones (top row) for single layer QBCP under periodic strain
field with $p6mm$ symmetry. (a-b) correspond to the 4 flat band case in Fig.
1(b) of the main text, (c-f) correspond to 6 flat bands case in Fig. 1(c) of
the main text. Note that in the 6 flat band case,
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(2)}$ and
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(3)}$ from Eq, 7 of main text are not
mutually orthogonal. Therefore, we orthogonalize them by redefining
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(3)}\rightarrow\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(3)}-\frac{\langle\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(2)}|\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(3)}\rangle}{\langle\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(2)}|\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(2)}\rangle}\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(2)}$.
To ensure same gauge between numerics and analytical expression, we further
change the basis such that
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(2)}(\mathcal{C}_{3z}\mathbf{r})=\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})$
and
$\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(3)}(\mathcal{C}_{3z}\mathbf{r})=e^{2\pi
i/3}\tilde{\psi}_{\mathbf{k}_{\text{max}}}^{(3)}(\mathbf{r})$.

##  S-4 Topology of the flat-bands

In the main text, we claim that among all the FBs polarized to the same sub-
lattice, only one carries Chern number $|C|=1$, the others are trivial. Here
we prove this statement. We start by re-deriving the well known result for
when there is a single sub-lattice polarized flat band, and then tackle the
cases when there are more than one FBs polarized to one sub-lattice.

###  S-4.1 Topology of two flatbands – single flat band on each sub-lattice

In the main text, we found that the FB WF can be written as

$\psi_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})=f_{\mathbf{k}}(z;\mathbf{r}_{0})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=e^{i(\mathbf{k}\cdot\mathbf{a}_{1}^{m})z/a_{1}^{m}}\frac{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}},\tau\right)}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}),$
(S52)

and the corresponding Bloch periodic part of the WF

$u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})=e^{-i\mathbf{k}\cdot\mathbf{r}}\psi_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})=e^{-i\mathbf{k}\cdot\mathbf{r}}e^{i(\mathbf{k}\cdot\mathbf{a}_{1}^{m})z/a_{1}^{m}}\frac{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}},\tau\right)}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})=e^{-i(\mathbf{b}_{2}^{m}\cdot\mathbf{r})k/b_{2}^{m}}\frac{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}},\tau\right)}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r}),$
(S53)

where $\vartheta(z,\tau)=-i\sum_{n=-\infty}^{\infty}(-1)^{n}e^{\pi
i\tau(n+1/2)^{2}+\pi i(2n+1)z}$ is the Jacobi theta function of first kind,
$\mathbf{a}_{i}^{m}$ are Moiré lattice vectors, $\mathbf{b}_{i}^{m}$ are the
corresponding reciprocal lattice vectors
($\mathbf{a}_{i}^{m}\cdot\mathbf{b}_{j}^{m}=2\pi\delta_{ij}$),
$a_{i}^{m}=(\mathbf{a}_{i}^{m})_{x}+i(\mathbf{a}_{i}^{m})_{y}$,
$b_{i}^{m}=(\mathbf{b}_{i})_{x}+i(\mathbf{b}_{i})_{y}$,
$z_{0}=(\mathbf{r}_{0})_{x}+i(\mathbf{r}_{0})_{y}$ is complexification of the
position of the zero of $\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$,
$k=k_{x}+ik_{y}$, and $\tau=a_{2}^{m}/a_{1}^{m}$. Clearly, the function
$u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})$ is holomorphic function of $k$
($\partial_{k^{*}}u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})=0$).
Furthermore, the function $u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})$ is
smooth throughout the mBZ but not periodic. The Chern number can be calculated
by taking line integral of Berry connection around the Brillouin zone. To this
end, first notice that the norm of the function $N_{\mathbf{k}}^{-2}=\langle
u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})|u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})\rangle\equiv\int_{\text{moir\'{e}
unit
cell}}d^{2}\mathbf{r}\,u_{\mathbf{k}+\mathbf{k}_{0}^{m}}^{*}(\mathbf{r})u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})\neq
1$. Define normalized function
$u_{N,\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})=N_{\mathbf{k}}u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})$.
Then, the Berry connection is defined as $\mathbf{A}_{\mathbf{k}}\equiv
i\langle
u_{N,\mathbf{k}+\mathbf{k}_{0}}|\bm{\nabla}_{\mathbf{k}}u_{N,\mathbf{k}+\mathbf{k}_{0}}\rangle=i\int_{\text{moir\'{e}
unit
cell}}d^{2}\mathbf{r}\,u_{N,\mathbf{k}+\mathbf{k}_{0}^{m}}^{*}(\mathbf{r})\bm{\nabla}_{\mathbf{k}}u_{N,\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})$.
Hence the Chern number is

$\begin{split}C&=\frac{1}{2\pi}(\int_{3}-\int_{1}+\int_{4}-\int_{2})d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}\\\
&=\frac{1}{2\pi}\int_{\mathbf{0}\rightarrow\mathbf{b}_{1}^{m}}d\mathbf{k}\cdot(\mathbf{A}_{\mathbf{k}}-\mathbf{A}_{\mathbf{k}+\mathbf{b}_{2}^{m}})+\frac{1}{2\pi}\int_{\mathbf{0}\rightarrow\mathbf{b}_{2}^{m}}d\mathbf{k}\cdot(\mathbf{A}_{\mathbf{k}+\mathbf{b}_{1}^{m}}-\mathbf{A}_{\mathbf{k}})\end{split}$
(S54)

where paths 1, 2, 3 and 4 are shown in Fig. S2. However, paths $i$ and $i+2$
are related by a moiré recirprocal lattice vector $\mathbf{b}_{1}^{m}$, hence
we can try to find the relations between
$\mathbf{A}_{\mathbf{k}+\mathbf{b}_{i}^{m}}$ and $\mathbf{A}_{\mathbf{k}}$.
Using the expression for $u_{\mathbf{k}+\mathbf{k}_{0}}(\mathbf{r})$, we find

$\begin{split}u_{\mathbf{k}+\mathbf{k}_{0}^{m}+\mathbf{b}_{2}^{m}}(\mathbf{r})&=e^{-i(\mathbf{b}_{2}^{m}\cdot\mathbf{r})(k+b_{2}^{m})/b_{2}^{m}}\frac{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k+b_{2}^{m}}{b_{2}^{m}},\tau\right)}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})\\\
&=-e^{-i\mathbf{b}_{2}^{m}\cdot\mathbf{r}}u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})\equiv
e^{i\theta_{\mathbf{k},\mathbf{b}_{2}^{\color[rgb]{0,0,0}{m}}}}e^{-i\mathbf{b}_{2}^{m}\cdot\mathbf{r}}u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r}),\end{split}$
(S55)

where we used the fact that $\vartheta(z+1,\tau)=-\vartheta(z,\tau)$. Also,

$\begin{split}&u_{\mathbf{k}+\mathbf{k}_{0}^{m}+\mathbf{b}_{1}^{m}}(\mathbf{r})\\\
&=e^{-i(\mathbf{b}_{2}^{m}\cdot\mathbf{r})(k+b_{1}^{m})/b_{2}^{m}}\frac{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k+b_{1}^{m}}{b_{2}^{m}},\tau\right)}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})\\\
&=e^{-i\mathbf{b}_{1}^{m}\cdot\mathbf{r}}e^{i(\mathbf{b}_{1}^{m}-\frac{b_{1}^{m}}{b_{2}^{m}}\mathbf{b}_{2}^{m})\cdot\mathbf{r}}(-i)\sum_{n=-\infty}^{\infty}(-1)^{n}e^{\pi
i\frac{a_{2}^{m}}{a_{1}^{m}}(n+1/2)^{2}+\pi
i(2n+1)(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k+b_{1}^{m}}{b_{2}^{m}})}\frac{1}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})\\\
&=e^{-i\mathbf{b}_{1}^{m}\cdot\mathbf{r}}e^{i(\mathbf{b}_{1}^{m}-\frac{b_{1}^{m}}{b_{2}^{m}}\mathbf{b}_{2}^{m})\cdot\mathbf{r}}(-i)\sum_{n=-\infty}^{\infty}(-1)^{n}e^{\pi
i\frac{a_{2}^{m}}{a_{1}^{m}}[(n+1/2)^{2}+(2n+1)]+\pi
i(2n+1)(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}})}\frac{1}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})\\\
&=e^{-i\mathbf{b}_{1}^{m}\cdot\mathbf{r}}e^{i(\mathbf{b}_{1}^{m}-\frac{b_{1}^{m}}{b_{2}^{m}}\mathbf{b}_{2}^{m})\cdot\mathbf{r}}(-i)\sum_{n=-\infty}^{\infty}(-1)^{n}e^{\pi
i\frac{a_{2}^{m}}{a_{1}^{m}}[(n+3/2)^{2}-1]+\pi
i(2n+1)(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}})}\frac{1}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})\\\
&=e^{-i\mathbf{b}_{1}^{m}\cdot\mathbf{r}}e^{i(\mathbf{b}_{1}^{m}-\frac{b_{1}^{m}}{b_{2}^{m}}\mathbf{b}_{2}^{m})\cdot\mathbf{r}}e^{-\pi
i\frac{a_{2}^{m}}{a_{1}^{m}}}(-1)e^{-2\pi
i(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}})}(-i)\sum_{n=-\infty}^{\infty}(-1)^{n+1}e^{\pi
i\frac{a_{2}^{m}}{a_{1}^{m}}[(n+3/2)^{2}]+\pi
i(2n+3)(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}})}\frac{\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})}{\vartheta\left(\frac{z-z_{0}}{a_{1}^{m}},\tau\right)}\\\
&=e^{i(\mathbf{b}_{1}^{m}-\frac{b_{1}^{m}}{b_{2}^{m}}\mathbf{b}_{2}^{m})\cdot\mathbf{r}}e^{-\pi
i\frac{a_{2}^{m}}{a_{1}^{m}}}(-1)e^{-2\pi
i(\frac{z-z_{0}}{a_{1}^{m}}-\frac{k}{b_{2}^{m}})}e^{-i\mathbf{b}_{1}^{m}\cdot\mathbf{r}}u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r}),\end{split}$
(S56)

where we used the relations that if $a_{1}^{m}=a^{m}e^{i\phi_{1}}$ and
$a_{2}^{m}=a^{m}e^{i\phi_{2}}$ ($a^{m}$ is length of the vector), we have
$b_{1}^{m}=\frac{2\pi
a^{m}}{|\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m}|}e^{i\phi_{2}-i\pi/2}$ and
$b_{2}^{m}=\frac{2\pi
a^{m}}{|\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m}|}e^{i\phi_{1}+i\pi/2}$, and
hence $b_{2}^{m}/b_{1}^{m}=-a_{1}^{m}/a_{2}^{m}$. To further simplify the
expression, we notice that

$\begin{split}\left(\mathbf{b}_{1}^{m}-\frac{b_{1}^{m}}{b_{2}^{m}}\mathbf{b}_{2}^{m}\right)\cdot\mathbf{r}&=\frac{1}{2}\left(b_{1}^{m}z^{*}+(b_{1}^{m})^{*}z-\frac{b_{1}^{m}}{b_{2}^{m}}(b_{2}^{m}z^{*}+(b_{2}^{m})^{*}z)\right)\\\
&=\frac{1}{2}\left((b_{1}^{m})^{*}-\frac{b_{1}^{m}}{b_{2}^{m}}(b_{2}^{m})^{*}\right)z\\\
&=i\frac{(\mathbf{b}_{1}^{m}\times\mathbf{b}_{2}^{m})\cdot\hat{z}}{b_{2}^{m}}z\\\
&=i(-\frac{2\pi i}{a_{1}^{m}})z=\frac{2\pi z}{a_{1}^{m}},\end{split}$ (S57)

where we used
$b_{2}^{m}(b_{1}^{m})^{*}-b_{1}^{m}(b_{2}^{m})^{*}=2i(\mathbf{b}_{1}^{m}\times\mathbf{b}_{2}^{m})\cdot\hat{z}$
and
$(\mathbf{b}_{1}^{m}\times\mathbf{b}_{2}^{m})\cdot\hat{z}=\frac{4\pi^{2}}{\hat{z}\cdot(\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m})}$.
These together imply that

$\begin{split}u_{\mathbf{k}+\mathbf{k}_{0}^{m}+\mathbf{b}_{1}^{m}}(\mathbf{r})=-e^{-\pi
i\frac{a_{2}^{m}}{a_{1}^{m}}}e^{2\pi i\frac{z_{0}}{a_{1}^{m}}}e^{2\pi
i\frac{k}{b_{2}^{m}}}e^{-i\mathbf{b}_{1}^{m}\cdot\mathbf{r}}u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r})\equiv
e^{i\theta_{\mathbf{k},\mathbf{b}_{1}^{\color[rgb]{0,0,0}{m}}}}e^{-i\mathbf{b}_{1}^{m}\cdot\mathbf{r}}u_{\mathbf{k}+\mathbf{k}_{0}^{m}}(\mathbf{r}).\end{split}$
(S58)

Eqs. S58 and S55 imply that

$\begin{split}\frac{N_{\mathbf{k}+\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}}{N_{\mathbf{k}}}&=e^{\Im(\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}})}.\\\
u_{N,\mathbf{k}+\mathbf{k}_{0}^{m}+\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}&=e^{i\Re(\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}})}e^{-i\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}\cdot\mathbf{r}}u_{N,\mathbf{k}+\mathbf{k}_{0}^{m}}.\end{split}$
(S59)

Using this, one finds
$\mathbf{A}_{\mathbf{k}+\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}=\mathbf{A}_{\mathbf{k}}-\bm{\nabla}_{\mathbf{k}}\Re(\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}})$.
Notice that in our case,

$\begin{split}\bm{\nabla}_{\mathbf{k}}\Re(\theta_{\mathbf{k},\mathbf{b}_{2}^{\color[rgb]{0,0,0}{m}}})&=0\\\
\bm{\nabla}_{\mathbf{k}}\Re(\theta_{\mathbf{k},\mathbf{b}_{1}^{\color[rgb]{0,0,0}{m}}})&=\frac{\hat{z}\cdot(\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m})}{a^{m}}(-\sin\phi_{1},\cos\phi_{1}),\end{split}$
(S60)

where we used the expressions $k=k_{x}+ik_{y}$ and $b_{2}^{m}=\frac{2\pi
a^{m}}{|\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m}|}e^{i\phi_{1}+i\pi/2}$.
Therefore, using Eq. S54 the Chern number is

$\begin{split}C&=\frac{1}{2\pi}\int_{\mathbf{0}\rightarrow\mathbf{b}_{2}^{m}}d\mathbf{k}\cdot\left(-\frac{\hat{z}\cdot(\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m})}{a^{m}}(-\sin\phi_{1},\cos\phi_{1})\right)\\\
&=\frac{1}{2\pi}\int_{0}^{1}dx\frac{2\pi
a^{m}}{|\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m}|}(-\sin\phi_{1},\cos\phi_{1})\cdot\left(-\frac{\hat{z}\cdot(\mathbf{a}_{1}^{m}\times\mathbf{a}_{2}^{m})}{a^{m}}(-\sin\phi_{1},\cos\phi_{1})\right)\\\
&=-1,\end{split}$ (S61)

where we parametrized the line $\mathbf{0}\rightarrow\mathbf{b}_{2}^{m}$ by
$x$. Therefore, the Chern number is always $|C|=1$ (whereas the sign of the
Chern number depends on the convention).

###  S-4.2 Topology of more than two flatbands – multiple flat bands on each
sub-lattice

Figure S2: (a) If a WF is smooth at all $\mathbf{k}$ in the Brillouin zone
defined by the reciprocal lattice vectors $\mathbf{b}_{i}^{m}$, its Chern
number can be evaluated from the line integrals of the Berry connections along
the paths $1,2,3$ and $4$ as shown in Eq. S54. If, in some gauge, the WF is
singular at $\mathbf{k}_{\text{max}}$, the Brillouin zone needs to be divided
into patches $N$ and $S$, where $S$ contains $\mathbf{k}_{\text{max}}$. A
different gauge for the WF is needed in $S$ that is smooth in $S$. (b) A
noncontractible loop parallel to $\mathbf{b}_{2}^{m}$ with base momentum
$\mathbf{k}_{\text{base}}$. The Wilson loop is evaluated along such loops. For
numerical evaluation of Wilson loop, the noncontractible loop needs to be
discretized. The Wilson loop is then a function of $\mathbf{k}_{\text{base}}$,
which can be varied parallel to $\mathbf{b}_{1}^{m}$.

As was described in the main text, if there are multiple zeros of
$\psi_{\mathbf{k}_{\text{max}}}$ at
$\mathbf{r}_{0}^{(1)},\dots,\mathbf{r}_{0}^{(m)}$, then $m$ FB WFs
$\mathbf{f}_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$
can be constructed on one sublattice. These WFs are smooth functions of
$\mathbf{k}$, and are independent (but not mutually orthogonal) of each other
everywhere except at $\mathbf{k}=0$ since
$\mathbf{f}_{\mathbf{0}}(z;\mathbf{r}_{0}^{(i)})=1$ for all $i$. For
concreteness and simplicity, let us focus on $m=2$, but the procedure can be
straightforwardly extended to any integer value of $m$. Similar to what was
shown in the main text, after orthogonalization,

$\tilde{u}_{\mathbf{k}}^{(i)}(\mathbf{r})=\begin{cases}u_{\mathbf{k}}^{(1)}(\mathbf{r})&\text{
if }i=1,\\\ u_{\mathbf{k}}^{(2)}(\mathbf{r})-\frac{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(2)}\rangle}{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(1)}\rangle}u_{\mathbf{k}}^{(1)}(\mathbf{r})&\text{
if }i=2,\\\ \end{cases}$ (S62)

where
$u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(i)}(\mathbf{r})=e^{-i\mathbf{k}\cdot\mathbf{r}}f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})=\tilde{f}_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$
is the periodic part of the full Bloch function. The function
$\tilde{u}_{\mathbf{k}}^{(1)}(\mathbf{r})$ is exactly the same as the one
discussed for two FB case, and has Chern number -1. The function
$\tilde{u}_{\mathbf{k}}^{(2)}(\mathbf{r})$ is more interesting. It is smooth
at any $\mathbf{k}$ except $\mathbf{k}_{\text{max}}$ due to the reason
mentioned above. Furthermore, around $\mathbf{k}_{\text{max}}$ the function
$\tilde{u}_{\mathbf{k}}^{(2)}(\mathbf{r})$ has the form

$\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})\approx
k(\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(2)})-\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)})-\langle\psi_{\mathbf{k}_{\text{max}}}|(\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(2)})-\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(1)}))\psi_{\mathbf{k}_{\text{max}}}\rangle)\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$
(S63)

where $k=k_{x}+ik_{y}$ and
$\tilde{f}_{\mathbf{0}}^{\prime}(\mathbf{r};\mathbf{r}_{0}^{(i)})=\left.\partial_{k}\tilde{f}_{\mathbf{k}}(\mathbf{r};\mathbf{r}_{0}^{(i)})\right|_{\mathbf{k}=\mathbf{0}}$,
where $\partial_{k}=(\partial_{k_{x}}-i\partial_{k_{y}})/2$. Its phase winds
around $\mathbf{k}_{\text{max}}$. However, the gauge
$e^{-i\varphi_{\mathbf{k}}}\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})$
(where $\varphi_{\mathbf{k}}=\text{arg}(k_{x}+ik_{y})$) is smooth. However,
the Berry curvature $F_{xy}(\mathbf{k})$ is gauge independent. So, while
calculating the Chern number using the usual formula
$C=\frac{1}{2\pi}\int_{\text{mBZ}}d^{2}\mathbf{k}F_{xy}(\mathbf{k})$ (where
$F_{xy}(\mathbf{k})=\partial_{k_{x}}A_{y}(\mathbf{k})-\partial_{k_{y}}A_{x}(\mathbf{k})$),
we can break the mBZ into two regions $S$ and $N$ (see Fig. S2), and use
$\tilde{u}^{S}_{\mathbf{k}+\mathbf{k}_{\text{max}}}(\mathbf{r})=e^{-i\varphi_{\mathbf{k}}}\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})$
and
$\tilde{u}^{N}_{\mathbf{k}+\mathbf{k}_{\text{max}}}(\mathbf{r})=\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})$
in the two regions respectively. Using the normalization factor
$N_{\mathbf{k}}$
($N_{\mathbf{k}}^{-2}=\langle\tilde{u}^{S}_{\mathbf{k}+\mathbf{k}_{\text{max}}}|\tilde{u}^{S}_{\mathbf{k}+\mathbf{k}_{\text{max}}}\rangle=\langle\tilde{u}^{N}_{\mathbf{k}+\mathbf{k}_{\text{max}}}|\tilde{u}^{N}_{\mathbf{k}+\mathbf{k}_{\text{max}}}\rangle$),
we define $\mathbf{A}^{S/N}_{\mathbf{k}}\equiv i\langle
N_{\mathbf{k}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{S/N}|\bm{\nabla}_{\mathbf{k}}N_{\mathbf{k}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{S/N}\rangle$.
Then, in the overlap region between the $S$ and $N$ regions, we have
$\mathbf{A}^{S}_{\mathbf{k}}=\mathbf{A}^{N}_{\mathbf{k}}+\bm{\nabla}_{\mathbf{k}}\varphi_{\mathbf{k}}$.
Using all of these, we get the expression for the Chern number of the second
band to be

$\begin{split}C&=\frac{1}{2\pi}\int_{\text{mBZ}}d^{2}\mathbf{k}F_{xy}(\mathbf{k})\\\
&=\frac{1}{2\pi}\left(\int_{S}d^{2}\mathbf{k}F_{xy}^{S}(\mathbf{k})+\int_{N}d^{2}\mathbf{k}F_{xy}^{N}(\mathbf{k})\right)\\\
&=\frac{1}{2\pi}\left(\int_{\partial
S}d\mathbf{k}\cdot(\mathbf{A}_{\mathbf{k}}^{S}-\mathbf{A}_{\mathbf{k}}^{N})+\int_{\partial\text{mBZ}}d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}\right)\\\
&=\frac{1}{2\pi}\left(\int_{\partial
S}d\mathbf{k}\cdot\bm{\nabla}_{\mathbf{k}}\varphi_{\mathbf{k}}+\int_{\partial\text{mBZ}}d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}\right)\\\
&=\frac{1}{2\pi}\int_{\partial\text{mBZ}}d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}+1.\end{split}$
(S64)

Now the integral over the edge of the mBZ is very similar to what we did in
the previous subsection. We know that

$\begin{split}&\frac{1}{2\pi}\int_{\partial\text{mBZ}}d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}\\\
=&\frac{1}{2\pi}(\int_{3}-\int_{1}+\int_{4}-\int_{2})d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}\\\
=&\frac{1}{2\pi}\int_{-\mathbf{b}_{1}^{m}/2\rightarrow\mathbf{b}_{1}^{m}/2}d\mathbf{k}\cdot(\mathbf{A}_{\mathbf{k}}^{N}-\mathbf{A}_{\mathbf{k}+\mathbf{b}_{2}^{m}}^{N})+\frac{1}{2\pi}\int_{-\mathbf{b}_{2}^{m}/2\rightarrow\mathbf{b}_{2}^{m}/2}d\mathbf{k}\cdot(\mathbf{A}_{\mathbf{k}+\mathbf{b}_{1}^{m}}^{N}-\mathbf{A}_{\mathbf{k}}^{N}).\end{split}$
(S65)

Now, using Eq. S58 and Eq. S55 we get

$\begin{split}&\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{N}(\mathbf{r})\\\
=&\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{(2)}(\mathbf{r})\\\
=&u_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{(2)}(\mathbf{r})-\frac{\langle
u_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{(1)}|u_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{(2)}\rangle}{\langle
u_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{(1)}|u_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{(1)}\rangle}u_{\mathbf{k}+\mathbf{k}_{\text{max}}+\mathbf{b}_{i}^{m}}^{(1)}(\mathbf{r})\\\
=&e^{-i\mathbf{b}_{i}^{m}\cdot\mathbf{r}}\left[e^{i\theta_{\mathbf{k},\mathbf{b}_{1}^{\color[rgb]{0,0,0}{m}}}^{(2)}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})-\frac{\langle
e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(1)}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}|e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(2)}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}\rangle}{\langle
e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(1)}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}|e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(1)}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}\rangle}e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(1)}}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}(\mathbf{r})\right]\\\
=&e^{-i\mathbf{b}_{i}^{m}\cdot\mathbf{r}}e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(2)}}\left[u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})-\frac{\langle
u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}|u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}\rangle}{\langle
u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}|u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}\rangle}u_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(1)}(\mathbf{r})\right]\\\
=&e^{-i\mathbf{b}_{i}^{m}\cdot\mathbf{r}}e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(2)}}\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{(2)}(\mathbf{r})\\\
=&e^{-i\mathbf{b}_{i}^{m}\cdot\mathbf{r}}e^{i\theta_{\mathbf{k},\mathbf{b}_{i}^{\color[rgb]{0,0,0}{m}}}^{(2)}}\tilde{u}_{\mathbf{k}+\mathbf{k}_{\text{max}}}^{N}(\mathbf{r}),\end{split}$
(S66)

which means that the ratio of wave-functions at the two sides of the mBZ
related by reciprocal lattice vector remains the same before and after the
orthogonalization in Eq. S62. But, we know from the previous subsection that
the integral
$\frac{1}{2\pi}\int_{\partial\text{mBZ}}d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}$
only depends on this ratio. Therefore, using our knowledge from the previous
subsection, we immediately obtain
$\frac{1}{2\pi}\int_{\partial\text{mBZ}}d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}=-1$.
Therefore, the Chern number of the second band is

$C=\frac{1}{2\pi}\int_{\partial\text{mBZ}}d\mathbf{k}\cdot\mathbf{A}_{\mathbf{k}}^{N}+1=-1+1=0.$
(S67)

A similar calculation shows that the Chern number of the third band would also
be zero in the case of $m=3$. This proves that the total Chern number of all
FB WFs polarized on the same sublattice is $C=\pm 1$.

###  S-4.3 Details of Wilson loop spectrum calculation

To numerically capture the topology of the flat bands, we evaluate Wilson loop
spectrum as shown in Figs. 1-2 of main text as well as Figs. S3-S7, S9. Here
we give details of the Wilson loop calculation. The Wilson loop at point
$\mathbf{k}_{\text{base}}$ along a noncontractible loop in the
$\mathbf{b}_{2}^{m}$ direction (see Fig. S2(b)) is defined as

$W_{\mathbf{k}_{\text{base}}+\mathbf{b}_{2}^{m}\leftarrow\mathbf{k}_{\text{base}}}^{mn}=\lim_{N\rightarrow\infty}\sum_{i_{1},\dots,i_{N-1}\in\text{FBs}}\langle
u_{\mathbf{k}_{\text{base}}+\mathbf{b}_{2}^{m}}^{m}|u_{\mathbf{k}_{\text{base}}+\frac{N-1}{N}\mathbf{b}_{2}^{m}}^{i_{N-1}}\rangle\langle
u_{\mathbf{k}_{\text{base}}+\frac{N-1}{N}\mathbf{b}_{2}^{m}}^{i_{N-1}}|\dots|u_{\mathbf{k}_{\text{base}}+\frac{1}{N}\mathbf{b}_{2}^{m}}^{i_{1}}\rangle\langle
u_{\mathbf{k}_{\text{base}}+\frac{1}{N}\mathbf{b}_{2}^{m}}^{i_{1}}|u_{\mathbf{k}_{\text{base}}}^{n}\rangle.$
(S68)

In numerical evaluations, $N$ is finite. The quantities
$\theta_{i}(\mathbf{k})=\Im\log\left(\text{eig}_{i}(W_{\mathbf{k}+\mathbf{b}_{2}^{m}\leftarrow\mathbf{k}})\right)$
are plotted as vary $\mathbf{k}$ between
$\mathbf{0}\rightarrow\mathbf{b}_{1}^{m}$ in the figures mentioned above. Due
to the anti-unitary particle-hole symmetry $\mathcal{P}$
($\mathcal{P}=\mathcal{S}\mathcal{T}$ in QBCP systems and
$\mathcal{P}=P\mathcal{C}_{2z}\mathcal{T}$ in TBG (see Sec. S-I)) that
satisfies $\mathcal{P}^{2}=-\mathds{1}$, the Wilson loop spectrum
$\theta_{i}(\mathbf{k})$ is (Kramers) doubly degenerate at
$\mathbf{k}=\mathbf{0},\frac{1}{2}\mathbf{b}^{m}$ and $\mathbf{b}^{m}$
TBGIIBernevigs. Furthermore, for $\mathcal{C}_{2z}\mathcal{T}$ symmetric
systems, the degeneracies at $\theta_{i}(\mathbf{k})=0,\pi$ are protected. In
all examples discussed here (except twisted bilayer checkerboard lattice,
which is discussed later), the winding of the Wilson loop spectrum is $\pm 1$,
which is consistent with the analysis done earlier in this section.

##  S-5 Ideal quantum geometry of the flat-bands

In this section we prove that the set of FBs polarized on each sublattice
satisfy ideal non-Abelian quantum geometry.

As was described in the main text, if there are multiple zeros of
$\psi_{\mathbf{k}_{\text{max}}}$ at
$\mathbf{r}_{0}^{(1)},\dots,\mathbf{r}_{0}^{(m)}$, then $m$ FB WFs
$\mathbf{f}_{\mathbf{k}}(z;\mathbf{r}_{0}^{(i)})\psi_{\mathbf{k}_{\text{max}}}(\mathbf{r})$
can be constructed on one sublattice. These WFs are smooth functions of
$\mathbf{k}$, and are independent (but not mutually orthogonal) of each other
everywhere except at $\mathbf{k}=0$ since
$\mathbf{f}_{\mathbf{0}}(z;\mathbf{r}_{0}^{(i)})=1$ for all $i$. For