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}$ (S69)

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. In the case of single FB per
sublattice, the ideal quantum geometry of the FB follows form the fact that
$u_{\mathbf{k}}^{(1)}(\mathbf{r})$ is a holomorphic function of $k$
claassen2015positions,ledwith2020fractionals. However, here the the
orthogonalization of $u_{\mathbf{k}}^{(2)}(\mathbf{r})$ from
$u_{\mathbf{k}}^{(1)}(\mathbf{r})$ gives factors like $\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(2)}\rangle$, which makes
$\tilde{u}_{\mathbf{k}}^{(2)}(\mathbf{r})$ non-holomorphic in $k$. However,
below we show that the two WFs $\tilde{u}_{\mathbf{k}}^{(1)}(\mathbf{r})$ and
$\tilde{u}_{\mathbf{k}}^{(2)}(\mathbf{r})$ together satisfies ideal non-
Abelian quantum geometry.

The non-Abelian Fubiny-Study metric
marzari1997maximallys,resta2011insulatings,marzari2012maximallys,xie2020topologys,ledwith2020fractionals
for the bands polarized on one sublattice is defined as:

$g_{\alpha\beta}^{mn}(\mathbf{k})=\Re\left[\langle\partial_{k_{\alpha}}\tilde{u}_{N,\mathbf{k}}^{(m)}|\left(\mathds{1}-\sum_{n_{1}=1}^{2}|\tilde{u}_{N,\mathbf{k}}^{(n_{1})}\rangle\langle\tilde{u}_{N,\mathbf{k}}^{(n_{1})}|\right)|\partial_{k_{\beta}}\tilde{u}_{N,\mathbf{k}}^{(n)}\rangle\right],$
(S70)

where the sum over $n_{1}$ is restricted to the FBs polarized on one
sublattice, $\tilde{u}_{N,\mathbf{k}}^{(m)}$ are the normalized FB WFs, and
$\langle f|g\rangle\equiv\int_{\text{moir\'{e} unit
cell}}d^{2}\mathbf{r}\,f^{*}(\mathbf{r})g(\mathbf{r})$. Then, the trace of the
non-Abelian Fubini-Study metric can be written in terms of unnormalized WFs as
the following:

$\text{tr}(g_{\alpha\beta}^{mn}(\mathbf{k}))=\sum_{n=1}^{2}\sum_{\alpha\in\\{x,y\\}}g_{\alpha\alpha}^{nn}(\mathbf{k})=\sum_{n=1}^{2}\sum_{\alpha\in\\{x,y\\}}\left(\frac{\langle\partial_{k_{\alpha}}\tilde{u}_{\mathbf{k}}^{(n)}|\partial_{k_{\alpha}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}}-\sum_{n^{\prime}=1}^{2}\frac{\langle\partial_{k_{\alpha}}\tilde{u}_{\mathbf{k}}^{(n)}|\tilde{u}_{\mathbf{k}}^{(n^{\prime})}\rangle\langle\tilde{u}_{\mathbf{k}}^{(n^{\prime})}|\partial_{k_{\alpha}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}||\tilde{u}_{\mathbf{k}}^{(n^{\prime})}||^{2}}\right),$
(S71)

where $||f||^{2}=\langle f|f\rangle$.

Now, we can use the expressions in Eq. (S69), and take advantage of the fact
that $u_{\mathbf{k}}^{(i)}(\mathbf{r})$ is holomorphic function in the
following way. First, writing,
$\partial_{k_{x}}=(\partial_{k}+\overline{\partial_{k}})$ and
$\partial_{k_{y}}=i(\partial_{k}-\overline{\partial_{k}})$ (here
$\partial_{k}=\frac{1}{2}(\partial_{k_{x}}-i\partial_{k_{y}})$ and
$\overline{\partial_{k}}=\frac{1}{2}(\partial_{k_{x}}+i\partial_{k_{y}})$), we
find

$\text{tr}(g_{\alpha\beta}^{mn}(\mathbf{k}))=2\sum_{n=1}^{2}\left(\frac{\langle\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}}+\frac{\langle\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}|\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}}-\sum_{n^{\prime}=1}^{2}\left(\frac{\langle\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}|\tilde{u}_{\mathbf{k}}^{(n^{\prime})}\rangle\langle\tilde{u}_{\mathbf{k}}^{(n^{\prime})}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}||\tilde{u}_{\mathbf{k}}^{(n^{\prime})}||^{2}}+\frac{\langle\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}|\tilde{u}_{\mathbf{k}}^{(n^{\prime})}\rangle\langle\tilde{u}_{\mathbf{k}}^{(n^{\prime})}|\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}||\tilde{u}_{\mathbf{k}}^{(n^{\prime})}||^{2}}\right)\right).$
(S72)

Furthermore, we have that

$\begin{split}&\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(1)}(\mathbf{r})=0\\\
&\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(2)}(\mathbf{r})=\overline{\partial_{k}}\left(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})\right)=-\frac{\langle\partial_{k}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})+\frac{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(2)}\rangle}{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(1)}\rangle^{2}}\langle\partial_{k}u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(1)}\rangle
u_{\mathbf{k}}^{(1)}(\mathbf{r})\\\
&\phantom{\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(2)}(\mathbf{r})}=-\frac{\langle\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}|\tilde{u}_{\mathbf{k}}^{(2)}\rangle}{\langle\tilde{u}_{\mathbf{k}}^{(1)}|\tilde{u}_{\mathbf{k}}^{(1)}\rangle}\tilde{u}_{\mathbf{k}}^{(1)}(\mathbf{r}).\end{split}$
(S73)

The last line also means that
$\langle\tilde{u}_{\mathbf{k}}^{(2)}|\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(2)}\rangle=0$
(sine $\tilde{u}_{\mathbf{k}}^{(2)}$ and $\tilde{u}_{\mathbf{k}}^{(1)}$ are
orthogonal) and
$\langle\tilde{u}_{\mathbf{k}}^{(1)}|\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(2)}\rangle=-\langle\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}|\tilde{u}_{\mathbf{k}}^{(2)}\rangle$.
Moreover, since
$\partial_{k}\tilde{u}_{\mathbf{k}}^{(2)}=\partial_{k}\left(u_{\mathbf{k}}^{(2)}-\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)}\right)=\partial_{k}u_{\mathbf{k}}^{(2)}-\frac{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(2)}\rangle}{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(1)}\rangle}\partial_{k}u_{\mathbf{k}}^{(1)}-\frac{\langle
u_{\mathbf{k}}^{(1)}|\partial_{k}u_{\mathbf{k}}^{(2)}\rangle}{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(1)}\rangle}u_{\mathbf{k}}^{(1)}+\frac{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(2)}\rangle}{\langle
u_{\mathbf{k}}^{(1)}|u_{\mathbf{k}}^{(1)}\rangle^{2}}\langle
u_{\mathbf{k}}^{(1)}|\partial_{k}u_{\mathbf{k}}^{(1)}\rangle
u_{\mathbf{k}}^{(1)}$, we have
$\langle\tilde{u}_{\mathbf{k}}^{(1)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(2)}\rangle=0$.
Plugging all of these into the expression in Eq. (S72), we find the following
simplified expression

$\text{tr}(g_{\alpha\beta}^{mn}(\mathbf{k}))=2\left(\frac{||\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}||^{2}}{||\tilde{u}_{\mathbf{k}}^{(1)}||^{2}}+\frac{||\partial_{k}\tilde{u}_{\mathbf{k}}^{(2)}||^{2}}{||\tilde{u}_{\mathbf{k}}^{(2)}||^{2}}-\frac{|\langle\tilde{u}_{\mathbf{k}}^{(1)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}\rangle|^{2}}{||\tilde{u}_{\mathbf{k}}^{(1)}||^{4}}-\frac{|\langle\tilde{u}_{\mathbf{k}}^{(2)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(2)}\rangle|^{2}}{||\tilde{u}_{\mathbf{k}}^{(2)}||^{4}}-\frac{|\langle\tilde{u}_{\mathbf{k}}^{(2)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}\rangle|^{2}}{||\tilde{u}_{\mathbf{k}}^{(1)}||^{2}||\tilde{u}_{\mathbf{k}}^{(2)}||^{2}}\right).$
(S74)

Clearly, $\text{tr}(g_{\alpha\beta}^{mn}(\mathbf{k}))>0$.

Similarly, the expression for the trace of the non-Abelian Berry curvature is

$\begin{split}\text{tr}(F_{xy}^{mn}(\mathbf{k}))&=\sum_{n=1}^{2}F_{xy}^{nn}(\mathbf{k})=\sum_{n=1}^{2}i\left(\langle\partial_{k_{x}}\tilde{u}_{\mathbf{k}}^{(n)}|\partial_{k_{y}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle-\langle\partial_{k_{y}}\tilde{u}_{{\mathbf{k}}}^{(n)}|\partial_{k_{x}}\tilde{u}_{{\mathbf{k}}}^{(n)}\rangle\right)\\\
&=\sum_{n=1}^{2}i\left[\left(\frac{\langle\partial_{k_{x}}\tilde{u}_{\mathbf{k}}^{(n)}|\partial_{k_{y}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}}-(x\leftrightarrow
y)\right)-\left(\frac{\langle\partial_{k_{x}}\tilde{u}_{\mathbf{k}}^{(n)}|\tilde{u}_{\mathbf{k}}^{(n)}\rangle\langle\tilde{u}_{\mathbf{k}}^{(n)}|\partial_{k_{y}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{4}}-(x\leftrightarrow
y)\right)\right]\\\
&=\sum_{n=1}^{2}2i\left[\left(\frac{\langle\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}}-\frac{\langle\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}|\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{2}}\right)-\left(\frac{\langle\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}|\tilde{u}_{\mathbf{k}}^{(n)}\rangle\langle\tilde{u}_{\mathbf{k}}^{(n)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{4}}-\frac{\langle\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}|\tilde{u}_{\mathbf{k}}^{(n)}\rangle\langle\tilde{u}_{\mathbf{k}}^{(n)}|\overline{\partial_{k}}\tilde{u}_{\mathbf{k}}^{(n)}\rangle}{||\tilde{u}_{\mathbf{k}}^{(n)}||^{4}}\right)\right].\end{split}$
(S75)

Now, using the identities in Eq. (S73) and below it, we simplify the
expression to get the following

$\text{tr}(F_{xy}^{mn}(\mathbf{k}))=2i\left(\frac{||\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}||^{2}}{||\tilde{u}_{\mathbf{k}}^{(1)}||^{2}}+\frac{||\partial_{k}\tilde{u}_{\mathbf{k}}^{(2)}||^{2}}{||\tilde{u}_{\mathbf{k}}^{(2)}||^{2}}-\frac{|\langle\tilde{u}_{\mathbf{k}}^{(1)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}\rangle|^{2}}{||\tilde{u}_{\mathbf{k}}^{(1)}||^{4}}-\frac{|\langle\tilde{u}_{\mathbf{k}}^{(2)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(2)}\rangle|^{2}}{||\tilde{u}_{\mathbf{k}}^{(2)}||^{4}}-\frac{|\langle\tilde{u}_{\mathbf{k}}^{(2)}|\partial_{k}\tilde{u}_{\mathbf{k}}^{(1)}\rangle|^{2}}{||\tilde{u}_{\mathbf{k}}^{(1)}||^{2}||\tilde{u}_{\mathbf{k}}^{(2)}||^{2}}\right).$
(S76)

Comparing Eq. (S74) and Eq. (S76), we have

$\text{tr}(g_{\alpha\beta}^{mn}(\mathbf{k}))=|\text{tr}(F_{xy}^{mn}(\mathbf{k}))|.$
(S77)

##  S-6 Examples

Below, we show 5 new examples listed in the Fig. 3(b) of the main text: 2 FBs
in systems with QBCP under periodic strain having $p3$ and $p4$ space group
symmetries, 4 FBs in systems with QBCP under periodic strain having $p4$,
$p4mm$ and $p4gm$ space group symmetry.

###  S-6.1 2 flat bands in single layer system with QBCP under moiré potential
with space group symmetry $p3$

Figure S3:  Flat bands in system with QBCP under moiré potential
$\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}=\alpha)=\frac{\alpha}{2}\sum_{n=1}^{3}e^{i(1-n)\phi}\;\;\exp\left(-i(\mathbf{b}_{n}^{m}\cdot\mathbf{r}+\phi_{1})\right)$.
Here, $\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 system for $\phi_{1}\neq 2\pi m/3$ ($m\in\mathds{Z}$) has
$p3$ space group symmetry. For $\phi_{1}=2\pi m/3$ ($m\in\mathds{Z}$), it has
$p31m$ symmetry. (a) Band structure showing 2 exact flat bands at
$\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=4.58$ and
$\phi_{1}=0.7962$ along the high symmetry path in the Moiré Brillouin zone.
The eigen-energy in the vertical axis is normalized as
$\tilde{E}=\frac{E}{|\mathbf{b}^{m}|^{2}}$. (b) Density plots of
$|\psi_{{\Gamma}^{m}}(\mathbf{r})|$. The white dashed line marks the boundary
of the moiré unit cell. The dark points indicate position of zeros of
$\psi_{\Gamma^{m}}(\mathbf{r})$. Clearly, there is only one zero of
$\psi_{\Gamma^{m}}(\mathbf{r})$ at an HSP, namely the corner, in the unit
cell. (c) Wilson loop spectrum
$\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of the flat
bands in (a). (d) Bandwidth of the middle two bands $\ln\tilde{E}_{w}$ as a
function of $\phi_{1}$ and $\tilde{\alpha}$ in polar coordinate of
$\tilde{\alpha}^{2}$ (radius) and $\phi_{1}$ (polar angle). The dark points in
the plot imply flat bands. Since the system has $p3$ symmetry (except the
special lines $\phi_{1}=2\pi m/3$), the co-dimension of the tuning parameter
to obtain flat bands is 2; hence we see flat-bands occurring at isolated
points in the $\tilde{\alpha}^{2}-\phi_{1}$ plane.

###  S-6.2 2 flat bands in single layer system with QBCP under moiré potential
with space group symmetry $p4$

Figure S4:  Flat bands in system with QBCP under moiré potential
$\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}=\alpha)=\alpha\exp(-i\phi)\sum_{n=1}^{2}(-1)^{1-n}\;\;\cos\left(\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right)$
having space group symmetry $p4$. Here,
$\mathbf{b}^{m}_{1}=\frac{2\pi}{a^{m}}(1,0)$ and
$\mathbf{b}^{m}_{2}=\frac{2\pi}{a^{m}}(0,1)$ are the reciprocal lattice
vectors and $a^{m}$ is the lattice constant of the superlattice. (a) Band
structure showing 2 exact flat bands at
$\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=4.83$ and $\phi=1.203067$
along the high symmetry path in the Moiré Brillouin zone. The eigen-energy in
the vertical axis is normalized as $\tilde{E}=\frac{E}{|\mathbf{b}^{m}|^{2}}$.
(b) Density plots of $|\psi_{{\Gamma}^{m}}(\mathbf{r})|$ (normalized by its
maximum). The white dashed line marks the boundary of the moiré unit cell. The
dark points indicate position of zeros of $\psi_{\Gamma^{m}}(\mathbf{r})$.
Clearly, there is only one zero of $\psi_{\Gamma^{m}}(\mathbf{r})$ at an HSP,
namely the corner, in the unit cell. (c) Wilson loop spectrum
$\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of the flat
bands in (a).

###  S-6.3 4 flat bands in single layer system with QBCP under moiré potential
with space group symmetry $p4$

Figure S5:  Flat bands in system with QBCP under moiré potential
$\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}=\alpha)=\alpha\exp(-i\phi)\sum_{n=1}^{2}(-1)^{1-n}\;\;\cos\left(\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right)$
having space group symmetry $p4$. Here,
$\mathbf{b}^{m}_{1}=\frac{2\pi}{a^{m}}(1,0)$ and
$\mathbf{b}^{m}_{2}=\frac{2\pi}{a^{m}}(0,1)$ are the reciprocal lattice
vectors and $a^{m}$ is the lattice constant of the superlattice. (a) Band
structure showing 4 exact flat bands at
$\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=8.41$ and
$\phi=0.98130229$ along the high symmetry path in the Moiré Brillouin zone.
The eigen-energy in the vertical axis is normalized as
$\tilde{E}=\frac{E}{|\mathbf{b}^{m}|^{2}}$. (b) Density plots of
$|\psi_{{\Gamma}^{m}}(\mathbf{r})|$ (normalized by its maximum). The white
dashed line marks the boundary of the moiré unit cell. The dark points
indicate position of zeros of $\psi_{\Gamma^{m}}(\mathbf{r})$. Clearly, there
are two zeros of $\psi_{\Gamma^{m}}(\mathbf{r})$ at HSPs, namely center of the
edges, in the unit cell. (c) Wilson loop spectrum
$\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of the flat
bands in (a).

###  S-6.4 4 flat bands in single layer system with QBCP under moiré potential
with space group symmetry $p4gm$

Figure S6:  Flat bands in system with QBCP under moiré potential
$\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}=\alpha)=i\alpha\sum_{n=1}^{2}(-1)^{1-n}\;\;\cos\left(\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right)$
having space group symmetry $p4gm$. Here,
$\mathbf{b}^{m}_{1}=\frac{2\pi}{a^{m}}(1,0)$ and
$\mathbf{b}^{m}_{2}=\frac{2\pi}{a^{m}}(0,1)$ are the reciprocal lattice
vectors and $a^{m}$ is the lattice constant of the superlattice. Notice that
in addition to $\mathcal{C}_{4z}$, this system has glide symmetry
$\mathcal{G}_{10}=\\{\mathcal{M}_{10}|\frac{1}{2},\frac{1}{2}\\}$:
$\mathcal{D}_{U}(\mathcal{G}_{10}\mathbf{r};\bm{\alpha}=\alpha)=-\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}=\alpha)=\mathcal{D}_{U}^{*}(\mathbf{r};\bm{\alpha}=\alpha)$
for $\alpha\in\mathds{R}$ (by $\mathcal{M}_{10}$ we mean the mirror whose
normal is in the direction of the lattice vector $\mathbf{a}_{1}^{m}$, the
translation part of the glide is $(\mathbf{a}_{1}^{m}+\mathbf{a}_{2}^{m})/2$).
(a) Band structure showing 4 exact flat bands at
$\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=2.24$ along the high
symmetry path in the Moiré Brillouin zone. The eigen-energy in the vertical
axis is normalized as $\tilde{E}=\frac{E}{|\mathbf{b}^{m}|^{2}}$. (b) Density
plots of $|\psi_{{\Gamma}^{m}}(\mathbf{r})|$ (normalized by its maximum). The
white dashed line marks the boundary of the moiré unit cell. The dark points
indicate position of zeros of $\psi_{\Gamma^{m}}(\mathbf{r})$. Clearly, there
are two zeros of $\psi_{\Gamma^{m}}(\mathbf{r})$ at HSPs, namely center of the
edges, in the unit cell. (c) Wilson loop spectrum
$\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of the flat
bands in (a).

###  S-6.5 4 flat bands in single layer system with QBCP under moiré potential
with space group symmetry $p4mm$

Figure S7:  Flat bands in system with QBCP under moiré potential
$\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}=\alpha)=\alpha\sum_{n=1}^{2}(-1)^{1-n}(\cos\left(\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right)-\cos\left(2\mathbf{b}_{n}^{m}\cdot\mathbf{r}\right))$
having space group symmetry $p4mm$. Here,
$\mathbf{b}^{m}_{1}=\frac{2\pi}{a^{m}}(1,0)$ and
$\mathbf{b}^{m}_{2}=\frac{2\pi}{a^{m}}(0,1)$ are the reciprocal lattice
vectors and $a^{m}$ is the lattice constant of the superlattice. Notice that
in addition to $\mathcal{C}_{4z}$, this system has mirror symmetry
$\mathcal{M}_{10}$:
$\mathcal{D}_{U}(\mathcal{M}_{10}\mathbf{r};\bm{\alpha}=\alpha)=\mathcal{D}_{U}(\mathbf{r};\bm{\alpha}=\alpha)=\mathcal{D}_{U}^{*}(\mathbf{r};\bm{\alpha}=\alpha)$
for $\alpha\in\mathds{R}$. (a) Band structure showing 4 exact flat bands at
$\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=-2.62$ along the high
symmetry path in the Moiré Brillouin zone. The eigen-energy in the vertical
axis is normalized as $\tilde{E}=\frac{E}{|\mathbf{b}^{m}|^{2}}$. (b) Density
plots of $|\psi_{{\Gamma}^{m}}(\mathbf{r})|$. The white dashed line marks the
boundary of the moiré unit cell. The dark points indicate position of zeros of
$\psi_{\Gamma^{m}}(\mathbf{r})$. Clearly, there are two zeros of
$\psi_{\Gamma^{m}}(\mathbf{r})$ at HSPs, namely center of the edges, in the
unit cell. (c) Wilson loop spectrum
$\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of the flat
bands in (a).

##  S-7 Twisted bilayer checkerboard lattice (TBCL) is two uncoupled copies of
single layer QBCP system under periodic strain field

Figure S8: TBCL. (a) Moiré unit cell of TBCL system is plotted in blue dashed
line. Moiré unit cell of single layer QBCP system is plotted in red dashed
line. $\mathbf{a}^{m}_{i}$ denotes the Moiré lattice vector of TBCL system. AA
region is marked by magenta disk and AB region is marked by green disk.
$\tilde{\mathbf{a}}^{m}_{i}$ denotes the Moiré lattice vector of the single
layer QBCP Hamiltonians that the TBCL Hamiltonian can be decomposed into (see
Eq. (S84)). (b) Moiré BZ of TBCL system is plotted in blue solid line. Moiré
BZ of single layer QBCP system is plotted in red solid line.
$\mathbf{b}^{m}_{i}$ denotes the Moiré reciprocal lattice vector of TBCL
system. $\tilde{\mathbf{b}}^{m}_{i}$ denotes the Moiré reciprocal lattice
vector of the single layer QBCP system. (c) Band structure of TBCL system with
2 exact flat bands at
$\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=0.13$. (d) Density plot of
$|\psi_{\Gamma^{m}}(\mathbf{r})|$ (normalized by its maximum).
$|\psi_{\Gamma^{m}}(\mathbf{r})|$ has no zero in the unit cell. (e), (f)
Density plot of $|\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})|$ and
$|\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})|$ (normalized by their respective
maximum). See the text below Eq. (S86) for the definition of these two
functions. The zero of $|\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})|$ in (f) is
shifted by
$(\tilde{\mathbf{a}_{1}^{m}}+\tilde{\mathbf{a}_{2}^{m}})/2=\mathbf{a}_{2}^{m}$
from the zero of $|\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})|$ in (e). (g) Berry
Curvature distribution $\Omega(\mathbf{k}_{n})$ (normalized by its average) of
the FBs of TBCL Hamiltonian plotted within the moiré BZ of TBCL Hamiltonian.
Clearly, the Berry Curvature distribution actually has a smaller periodicity
than the moiré BZ of TBCL system.

Checkerboard lattice has $\mathcal{C}_{(n=4)z}$ and $\mathcal{T}$ protected
QBCP at the corner, $\mathbf{k}_{0}=M$ point of the Brillouin zone. Upon
twisting the the two layers on top of each other, the $M$ point of one layer
gets mapped to $\mathbf{k}_{0}^{m}=\Gamma^{m}$, the $M$ point of the other
layer gets mapped of $M^{m}$ of the mBZ. Following the derivation of Sec.
SM.1D, we have $\rho(\mathcal{C}_{4z})=\text{Diag}\\{i,i,-i,-i\\}$,
$\rho(\mathcal{T})=\sigma_{x}\otimes\mathds{1}$. Furthermore, the moiré
potential $U_{1}(\mathbf{r})$ in Eq. (S24) satisfies (Eqs. (S25) and (S29))

$U_{1}(\mathcal{C}_{4z}\mathbf{r})=-U_{1}(\mathbf{r}),\;U_{1}(\mathbf{r})=U_{2}^{*}(\mathbf{r}),\;U_{1}(\mathbf{r})=\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}},$
(S78)

where $\mathbf{q}_{1}=(\mathbf{b}_{1}^{m}-\mathbf{b}_{2}^{m})/2$ as shown in
Fig. S8. Together, we have

$\begin{split}U_{1}(\mathcal{C}_{4z}\mathbf{r})=-U_{1}(\mathbf{r})&\Rightarrow\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathcal{C}_{4z}\mathbf{r}}=-\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}_{4z}^{-1}\mathbf{q}_{1}+\mathcal{C}_{4z}^{-1}\mathbf{b}^{m})\cdot\mathbf{r}}=-\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}_{1}^{m}+\mathcal{C}_{4z}^{-1}\mathbf{b}^{m})\cdot\mathbf{r}}=-\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}_{4z}(\mathbf{b}^{m}+\mathbf{b}^{m}_{1})}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}}=-\sum_{\mathbf{b}^{m}}a_{\mathbf{b}^{m}}e^{-i(\mathbf{q}_{1}+\mathbf{b}^{m})\cdot\mathbf{r}},\\\
&\Rightarrow
a_{\mathcal{C}_{4z}(\mathbf{b}^{m}+\mathbf{b}^{m}_{1})}=-a_{\mathbf{b}^{m}}\text{
for all reciprocal lattice vector }\mathbf{b}^{m}.\end{split}$ (S79)

Using above equation, starting from $a_{\mathbf{0}}\equiv\alpha$, we get

$\begin{split}&a_{\mathcal{C}_{4z}\mathbf{b}^{m}_{1}}=a_{\mathbf{b}^{m}_{2}}=-a_{\mathbf{0}}=-\alpha,\\\
&a_{\mathcal{C}_{4z}(\mathbf{b}^{m}_{1}+\mathbf{b}^{m}_{2})}=a_{\mathbf{b}^{m}_{2}-\mathbf{b}^{m}_{1}}=-a_{\mathbf{b}^{m}_{2}}=a_{\mathbf{0}}=\alpha,\\\
&a_{\mathcal{C}_{4z}(\mathbf{b}^{m}_{2})}=a_{-\mathbf{b}^{m}_{1}}=-a_{\mathbf{b}^{m}_{2}-\mathbf{b}^{m}_{1}}=-a_{\mathbf{0}}=-\alpha.\end{split}$
(S80)

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^{-i(\mathbf{q}_{1}+\mathbf{b}_{2}^{m})\cdot\mathbf{r}}+e^{-i(\mathbf{q}_{1}+\mathbf{b}_{2}^{m}-\mathbf{b}^{m}_{1})\cdot\mathbf{r}}-e^{-i(\mathbf{q}_{1}-\mathbf{b}^{m}_{1})\cdot\mathbf{r}})\\\
&=\alpha(e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}-e^{i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{i\mathbf{q}_{1}\cdot\mathbf{r}}-e^{-i\mathbf{q}_{2}\cdot\mathbf{r}})\\\
&=2\alpha(\cos(\mathbf{q}_{1}\cdot\mathbf{r})-\cos(\mathbf{q}_{2}\cdot\mathbf{r})).\end{split}$
(S81)

An additional mirror symmetry $\mathcal{M}_{x}$ with representation
$\rho(\mathcal{M}_{x})=\sigma_{x}\otimes\mathds{1}$, would result in

$\begin{split}&U_{1}(\mathcal{M}_{x}\mathbf{r})=U_{2}(\mathbf{r})=U_{1}^{*}(\mathbf{r})\\\
\Rightarrow&2\alpha(\cos(\mathbf{q}_{1}\cdot\mathcal{M}_{x}\mathbf{r})-\cos(\mathbf{q}_{2}\cdot\mathcal{M}_{x}\mathbf{r}))=2\alpha^{*}(\cos(\mathbf{q}_{1}\cdot\mathbf{r})-\cos(\mathbf{q}_{2}\cdot\mathbf{r}))\\\
\Rightarrow&2\alpha(\cos((\mathcal{M}_{x}\mathbf{q}_{1})\cdot\mathbf{r})-\cos((\mathcal{M}_{x}\mathbf{q}_{2})\cdot\mathbf{r}))=2\alpha^{*}(\cos(\mathbf{q}_{1}\cdot\mathbf{r})-\cos(\mathbf{q}_{2}\cdot\mathbf{r}))\\\
\Rightarrow&2\alpha(\cos(\mathbf{q}_{2}\cdot\mathbf{r})-\cos(\mathbf{q}_{1}\cdot\mathbf{r}))=2\alpha^{*}(\cos(\mathbf{q}_{1}\cdot\mathbf{r})-\cos(\mathbf{q}_{2}\cdot\mathbf{r}))\\\
\Rightarrow&\alpha=-\alpha^{*}\end{split}$ (S82)

Then, replace $\alpha\rightarrow i\alpha$ such that $\alpha\in\mathbb{R}$.
Lastly, performing the following transformation $\text{Diag}\\{e^{\pi
i/4}\mathchar 44\relax\penalty 0e^{\pi i/4}\mathchar 44\relax\penalty 0e^{-\pi
i/4}\mathchar 44\relax\penalty 0e^{-\pi
i/4}\\}\mathcal{H}_{TB}(\mathbf{r})\text{Diag}\\{e^{-\pi i/4}\mathchar
44\relax\penalty 0e^{-\pi i/4}\mathchar 44\relax\penalty 0e^{\pi i/4}\mathchar
44\relax\penalty 0e^{\pi i/4}\\}$ on the Hamiltonian in Eq. S24, we obtain

$\mathcal{H}_{TBCL}(\mathbf{r})=\begin{pmatrix}0&0&i(-2i\partial_{z})^{2}&iU_{1}^{*}(\mathbf{r})\\\
0&0&iU_{1}^{*}(\mathbf{r})&i(-2i\partial_{z})^{2}\\\
-i(-2i\overline{\partial_{z}})^{2}&-iU_{1}(\mathbf{r})&0&0\\\
-iU_{1}(\mathbf{r})&-i(-2i\overline{\partial_{z}})^{2}\ &0&0\
\end{pmatrix}\text{, }U_{1}(\mathbf{r})=2\alpha
i(\cos(\mathbf{q}_{1}\cdot\mathbf{r})-\cos(\mathbf{q}_{2}\cdot\mathbf{r})),\;\alpha\in\mathbb{R},$
(S83)

this is the Hamiltonian considered in li2022magics. Notice that there It was
shown in li2022magics (see also Fig. S8(c)), that for some magic values of
$\alpha$, 2 exact flat bands appear at the charge neutrality point. However,

  * •

these bands have Chern number $C=\pm 2$, which suggest that the FB WFs are not
simply
$f_{\mathbf{k}}(z;\mathbf{r}_{0})\psi_{\mathbf{k}_{0}^{m}}(\mathbf{r})$,
because if it were, then the Chern number would be $C=\pm 1$

  * •

It is also clear from Fig. S8(d), that the WF $\psi_{\Gamma^{m}}$ does not
have a zero at the magic value of $\alpha$ in the moiré unit cell.

These two points makes the FBs in TBCL intriguing. A major clue to solving
this problem comes from the Berry curvature distribution of the sublattice
polarized FB WF in Fig. S8(g). The periodicity of the Berry curvature
distribution in the reciprocal space is smaller than the reciprocal lattice
vectors. This alludes to the possibility that this model is two copies of
Chern number $C=\pm 1$ bands unfolded to a larger Brillouin zone. Below we
show that this is indeed the case. Consider the following transformation

$\mathcal{H}^{(s)}_{TBCL}(\mathbf{r})=U\mathcal{H}_{TBCL}(\mathbf{r})U^{\dagger}=\begin{pmatrix}0&i(-2i\partial_{z})^{2}+iU_{1}^{*}(\mathbf{r})&0&0\\\
-i(-2i\overline{\partial_{z}})^{2}-iU_{1}(\mathbf{r})&0&0&0\\\
0&0&0&i(-2i\partial_{z})^{2}-iU_{1}^{*}(\mathbf{r})\\\
0&0&-i(-2i\overline{\partial_{z}})^{2}+iU_{1}(\mathbf{r})&0\ \end{pmatrix},$
(S84)

where

$U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&1&0&0\\\ 0&0&1&1\\\ -1&1&0&0\\\
0&0&-1&1\end{pmatrix}.$ (S85)

Each of the diagonal blocks clearly corresponds to a single layer QBCP system
with $p4mm$ space group symmetry and moiré lattice vectors
$\tilde{\mathbf{a}}_{i}^{m}$ (and corresponding reciprocal lattice vector
$\tilde{\mathbf{b}}_{i}^{m}$) as shown in Fig. S8(a-b). However, the lattice
vectors $\mathbf{a}_{i}^{m}$ of $\mathcal{H}_{TBCL}$ are smaller (see Fig.
S8(a-b)) because $U_{1}(\mathbf{r}+\mathbf{a}_{i}^{m})=-U_{1}(\mathbf{r})$,
hence

$\mathcal{H}^{(s)}_{TBCL}(\mathbf{r}+\mathbf{a}_{i}^{m})=\sigma_{x}\otimes\mathds{1}\mathcal{H}^{(s)}_{TBCL}(\mathbf{r})\sigma_{x}\otimes\mathds{1}.$
(S86)

In fact, the diagonal blocks are exactly the same as the Hamiltonian that was
reported to host exact FBs for magic values of $\alpha$ in eugenio2022twisteds
(the Hamiltonian in eugenio2022twisteds is written in a $\pi/4$ rotated
coordinate system than the one here, which results in a difference in factor
of imaginary $i$ between the two). The magic value reported in
eugenio2022twisteds are the same as that reported in li2022magics after
nondimensionalization. Indeed, when we transform the TBCL WF
$\Psi_{\Gamma^{m}}(\mathbf{r})=\\{\psi_{\Gamma^{m}}(\mathbf{r}),\mathbf{0}\\}^{T}$
to $\Psi_{\Gamma^{m}}^{(s)}(\mathbf{r})=U\Psi_{\Gamma^{m}}(\mathbf{r})$, the
nonzero components of it,
$\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})=(\psi_{\Gamma^{m},1}(\mathbf{r})+\psi_{\Gamma^{m},2}(\mathbf{r}))/\sqrt{2}$
and
$\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})=(-\psi_{\Gamma^{m},1}(\mathbf{r})+\psi_{\Gamma^{m},2}(\mathbf{r}))/\sqrt{2}$)
(recall
$\psi_{\Gamma^{m}}(\mathbf{r})=\\{\psi_{\Gamma^{m},1}(\mathbf{r}),\psi_{\Gamma^{m},2}(\mathbf{r})\\}$
is a two component function for twisted bilayer systems) have zeros in the
unit cell defined by vectors $\tilde{\mathbf{a}}_{i}^{m}$ shifted from each
other by $\mathbf{a}_{1}^{m}$ for magic values of $\alpha$ as shown in Fig.
S8(e-f).

Clearly, from the Eq. S45 and subsequent discussion, we know that
$\psi_{\Gamma^{m},1}(\mathbf{r})$ has periodicity $\mathbf{a}_{i}^{m}$,
whereas $\psi_{\Gamma^{m},2}(\mathbf{r})$ has periodicity
$\tilde{\mathbf{a}}_{i}^{m}$; hence $\Psi_{\Gamma^{m},1/3}^{(s)}(\mathbf{r})$
have periodicity $\tilde{\mathbf{a}}_{i}^{m}$. Also, since
$\psi_{\Gamma^{m},2}(\mathbf{r}-\mathbf{a}_{1}^{m})=-\psi_{\Gamma^{m},2}(\mathbf{r})$,
we have
$\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}-\mathbf{a}_{1}^{m})=\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})$.
Using the zeros of $\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})$ and
$\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})$ at $\mathbf{r}_{0}=\mathbf{0}$ and
$\mathbf{r}_{0}=-\mathbf{a}_{1}^{m}$ respectively (see Fig. S8(e-f)), we can
construct the FB WF of $\mathcal{H}^{(s)}(\mathbf{r})$ as

$\Psi_{\mathbf{k}}^{(s)}(\mathbf{r})=\begin{Bmatrix}f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\\\
0\\\
e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix}\text{,
}f_{\mathbf{k}}(z;\mathbf{0})=e^{i(\mathbf{k}\cdot\tilde{\mathbf{a}}_{1}^{m})z/\tilde{a}_{1}^{m}}\frac{\vartheta\left(\frac{z}{\tilde{a}_{1}^{m}}-\frac{k}{\tilde{b}_{2}^{m}},\tilde{a}_{2}^{m}/\tilde{a}_{1}^{m}\right)}{\vartheta\left(\frac{z}{\tilde{a}_{1}^{m}},,\tilde{a}_{2}^{m}/\tilde{a}_{1}^{m}\right)}$
(S87)

where $f_{\mathbf{k}}(z;\mathbf{r}_{0})$ satisfies
$f_{\mathbf{k}}(z+\tilde{a}_{i}^{m};\mathbf{r}_{0})=e^{i\mathbf{k}\cdot\tilde{\mathbf{a}}_{i}^{m}}f_{\mathbf{k}}(z;\mathbf{r}_{0})$,
and
$f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})=f_{\mathbf{k}}(z;-\mathbf{a}_{1}^{m})$.
It can be easily checked that $\Psi_{\mathbf{k}}^{(s)}(\mathbf{r})$ satisfies
the Bloch periodicity
$\Psi_{\mathbf{k}}^{(s)}(\mathbf{r}+\mathbf{a}_{i}^{m})=e^{i\mathbf{k}\cdot\mathbf{a}_{i}^{m}}\sigma_{x}\otimes\mathds{1}\Psi_{\mathbf{k}}^{(s)}(\mathbf{r})$
corresponding to Eq. (S86):

$\displaystyle\Psi_{\mathbf{k}}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})$
$\displaystyle=$
$\displaystyle\begin{Bmatrix}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\\\
e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+2a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+2\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix}$ $\displaystyle=$
$\displaystyle\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+2a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+2\mathbf{a}_{1}^{m})\\\
0\\\
f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix}$ $\displaystyle=$
$\displaystyle\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+\tilde{a}_{1}^{m}-\tilde{a}_{2}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\tilde{\mathbf{a}}_{1}^{m}-\tilde{\mathbf{a}}_{2}^{m})\\\
0\\\
f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix}$ $\displaystyle=$
$\displaystyle\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}e^{i\mathbf{k}\cdot(\tilde{\mathbf{a}}_{1}^{m}-\tilde{\mathbf{a}}_{2}^{m})}f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\\\
0\\\
f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix},\begin{array}[]{c}\text{ since
}f_{\mathbf{k}}(z+\tilde{a}_{1}^{m}-\tilde{a}_{2}^{m};\mathbf{0})=e^{i\mathbf{k}\cdot(\tilde{\mathbf{a}}_{1}^{m}-\tilde{\mathbf{a}}_{2}^{m})}f_{\mathbf{k}}(z;\mathbf{0})\\\
\text{ and
}\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\tilde{\mathbf{a}}_{1}^{m}-\tilde{\mathbf{a}}_{2}^{m})=\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\end{array}$
(S88c) $\displaystyle=$ $\displaystyle
e^{i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\\\
0\\\
e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix},\text{ since
}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}e^{i\mathbf{k}\cdot(\tilde{\mathbf{a}}_{1}^{m}-\tilde{\mathbf{a}}_{2}^{m})}=e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}e^{i\mathbf{k}\cdot
2\tilde{\mathbf{a}}_{1}^{m}}=e^{i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}$
$\displaystyle=$ $\displaystyle
e^{i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}\sigma_{x}\otimes\mathds{1}\Psi_{\mathbf{k}}^{(s)}(\mathbf{r}),$
$\displaystyle\Psi_{\mathbf{k}}^{(s)}(\mathbf{r}+\mathbf{a}_{2}^{m})$
$\displaystyle=$
$\displaystyle\begin{Bmatrix}f_{\mathbf{k}}(z+a_{2}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{2}^{m})\\\
0\\\
e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+a_{1}^{m}+a_{2}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m}+\mathbf{a}_{2}^{m})\\\
0\end{Bmatrix}$ $\displaystyle=$
$\displaystyle\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+a_{1}^{m}+a_{2}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m}+\mathbf{a}_{2}^{m})\\\
0\\\
f_{\mathbf{k}}(z+a_{2}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{2}^{m})\\\
0\end{Bmatrix}$ $\displaystyle=$
$\displaystyle\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+\tilde{a}_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\tilde{\mathbf{a}}_{1}^{m})\\\
0\\\
f_{\mathbf{k}}(z+a_{1}^{m}+a_{2}^{m}-a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m}+\mathbf{a}_{2}^{m}-\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix}$ $\displaystyle=$
$\displaystyle\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}e^{i\mathbf{k}\cdot\tilde{\mathbf{a}}_{1}^{m}}f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\\\
0\\\
f_{\mathbf{k}}(z+a_{1}^{m}+\tilde{a}_{2}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m}+\tilde{\mathbf{a}}_{2}^{m})\\\
0\end{Bmatrix},\begin{array}[]{c}\text{ since
}f_{\mathbf{k}}(z+\tilde{a}_{1}^{m};\mathbf{0})=e^{i\mathbf{k}\cdot\tilde{\mathbf{a}}_{1}^{m}}f_{\mathbf{k}}(z;\mathbf{0})\\\
\text{ and
}\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\tilde{\mathbf{a}}_{1}^{m})=\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\end{array}$
(S88f) $\displaystyle=$
$\displaystyle\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}e^{i\mathbf{k}\cdot\mathbf{a}_{2}^{m}}f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\\\
0\\\
e^{i\mathbf{k}\cdot\tilde{\mathbf{a}}_{2}^{m}}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix},\begin{array}[]{c}\text{ since
}f_{\mathbf{k}}(z+a_{1}^{m}+\tilde{a}_{2}^{m};\mathbf{0})=e^{i\mathbf{k}\cdot\tilde{\mathbf{a}}_{2}^{m}}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\\\
\text{,
}\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{2}^{m}+\tilde{\mathbf{a}}_{1}^{m})=\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{2}^{m})\text{,
and
}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}e^{i\mathbf{k}\cdot\tilde{\mathbf{a}}_{1}^{m}}=e^{i\mathbf{k}\cdot\mathbf{a}_{2}^{m}}\end{array}$
(S88i) $\displaystyle=$ $\displaystyle
e^{i\mathbf{k}\cdot\mathbf{a}_{2}^{m}}\sigma_{x}\otimes\mathds{1}\begin{Bmatrix}f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})\\\
0\\\
e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\end{Bmatrix},\text{ since
}e^{i\mathbf{k}\cdot\tilde{\mathbf{a}}_{2}^{m}}=e^{i\mathbf{k}\cdot\mathbf{a}_{2}^{m}}e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}$
$\displaystyle=$ $\displaystyle
e^{i\mathbf{k}\cdot\mathbf{a}_{2}^{m}}\sigma_{x}\otimes\mathds{1}\Psi_{\mathbf{k}}^{(s)}(\mathbf{r}).$
(S88j)

Hence the WF

$\Psi_{\mathbf{k}}(\mathbf{r})=U^{\dagger}\Psi_{\mathbf{k}}^{(s)}(\mathbf{r})=\frac{1}{\sqrt{2}}\begin{Bmatrix}f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})-e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
f_{\mathbf{k}}(z;\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})+e^{-i\mathbf{k}\cdot\mathbf{a}_{1}^{m}}f_{\mathbf{k}}(z+a_{1}^{m};\mathbf{0})\psi_{\Gamma^{m},1}^{(s)}(\mathbf{r}+\mathbf{a}_{1}^{m})\\\
0\\\ 0\end{Bmatrix}$ (S89)

satisfies
$\mathcal{H}_{TBCL}(\mathbf{r})\Psi_{\mathbf{k}}(\mathbf{r})=\mathbf{0}$ and
have the correct Bloch periodicity. The FB WF polarized on the other
sublattice can be obtained as
$\sigma_{x}\otimes\mathds{1}\Psi_{\mathbf{k}}^{*}(\mathbf{r})$.

Lastly, we can calculate Chern number of the FB WF in the gauge of
$\Psi^{(s)}_{\mathbf{k}}(\mathbf{r})$. The two nonzero components both have
$f_{\mathbf{k}}(z;\mathbf{0})$, and hence would give Chern number $C=-1$ if we
integrate over the BZ given by $\tilde{\mathbf{b}}_{1}^{m}$ and
$\tilde{\mathbf{b}}_{2}^{m}$. But, due to Eq. (S86), the BZ over which we have
to integrate is $\mathbf{b}_{1}^{m}$ and $\mathbf{b}_{2}^{m}$, which is twice
as big as the one given by $\tilde{\mathbf{b}}_{1}^{m}$ and
$\tilde{\mathbf{b}}_{2}^{m}$. This is why the Chern number of these two FBs
are $C=\pm 2$.

###  S-7.1 TBCL type Hamiltonian with 4 flat bands

Figure S9: TBCL-4FBs. (a)Band structure of TBCL type system with 4 exact FBs
at $\tilde{\alpha}=\frac{\alpha}{|\mathbf{b}^{m}|^{2}}=0.56i$. (b) Zoom in of
the 4 FBs in (a), which shows "band sticking" along $X^{m}-Y^{m}$. (c) Density
plot of $|\psi_{\Gamma^{m}}(\mathbf{r})|$ (normalized by its maximum).
Clearly, there is only one zero of $\psi_{\Gamma^{m}}(\mathbf{r})$ at an HSP,
namely the corner, in the unit cell. Moiré unit cell of TBCL system is plotted
in white dashed line. (d),(e) Density plot of
$|\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})|$ and
$|\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})|$ (normalized by their respective
maximum). See the text below Eq. (S91) for the definition of these two
functions. Moiré unit cell of single layer QBCP system is plotted in white
dashed line. Clearly, there are two zeros of
$|\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})|$/$|\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})|$
at HSPs, namely center of the edges, in the unit cell. (f) Wilson loop
spectrum $\tilde{\theta}({\mathbf{k}})=\frac{\theta({\mathbf{k}})}{2\pi}$ of
the flat bands in (a).

We end the section by showing the example of a closely related Hamiltonian to
that of TBCL. Consider the Hamiltonian that is obtained by replacing
$U_{1}(\mathbf{r})\rightarrow iU_{1}(\mathbf{r})$ in Eq. (S83)

$\mathcal{H}(\mathbf{r})=\begin{pmatrix}0&0&i(-2i\partial_{z})^{2}&U_{1}^{*}(\mathbf{r})\\\
0&0&U_{1}^{*}(\mathbf{r})&i(-2i\partial_{z})^{2}\\\
-i(-2i\overline{\partial_{z}})^{2}&U_{1}(\mathbf{r})&0&0\\\
U_{1}(\mathbf{r})&-i(-2i\overline{\partial_{z}})^{2}\ &0&0\
\end{pmatrix}\text{, }U_{1}(\mathbf{r})=2\alpha
i(\cos(\mathbf{q}_{1}\cdot\mathbf{r})-\cos(\mathbf{q}_{2}\cdot\mathbf{r})),\;\alpha\in\mathbb{R}.$
(S90)

Remarkably, this Hamiltonian at a “magic” value of $\alpha$ has 4 FBs (see
Fig. S9(a-b)) with bands polarized to each sublattice possessing Chern number
$C=\pm 2$ as can be seen the Wilson loop spectrum in Fig. S9(f). Even more
curious is the the fact that the wave-function $\psi_{\Gamma^{m}}(\mathbf{r})$
has a single zero at corner of the unit cell (Fig. S9(c)). The Chern number as
well as the number of zeros are, once again, seemingly in disagreement with
the construction of FB WFs discussed in the main text. This can be understood
using the decomposition that was used for TBCL and symmetries of the
decomposed single layer Hamiltonian. Once we transform the Hamiltonian to

$\mathcal{H}^{(s)}(\mathbf{r})=U\mathcal{H}_{TBCL}(\mathbf{r})U^{\dagger}=\begin{pmatrix}0&i(-2i\partial_{z})^{2}+U_{1}^{*}(\mathbf{r})&0&0\\\
-i(-2i\overline{\partial_{z}})^{2}+U_{1}(\mathbf{r})&0&0&0\\\
0&0&0&i(-2i\partial_{z})^{2}-U_{1}^{*}(\mathbf{r})\\\
0&0&-i(-2i\overline{\partial_{z}})^{2}-U_{1}(\mathbf{r})&0\ \end{pmatrix}.$
(S91)

We find a new symmetry for each diagonal block. Same as in TBCL, each of the
diagonal blocks clearly corresponds to a single layer QBCP system with moiré
lattice vectors $\tilde{\mathbf{a}}_{i}^{m}$, whereas the lattice vectors
$\mathbf{a}_{i}^{m}$ of $\mathcal{H}$ are smaller (see Fig. S8(a-b)) because
$\mathcal{H}^{(s)}(\mathbf{r}+\mathbf{a}_{i}^{m})=\sigma_{x}\otimes\mathds{1}\mathcal{H}^{(s)}(\mathbf{r})\sigma_{x}\otimes\mathds{1}$.
Notice that each block has a glide symmetry
$\mathcal{G}_{10}=\\{\mathcal{M}_{10}|\frac{1}{2},\frac{1}{2}\\}$, where
$\mathcal{G}_{10}\mathbf{r}=\mathcal{G}_{10}(x,y)=\mathcal{M}_{10}(x,y)+\frac{1}{2}\tilde{\mathbf{a}}_{1}^{m}+\frac{1}{2}\tilde{\mathbf{a}}_{2}^{m}=(-y,-x)+\frac{1}{2}\tilde{\mathbf{a}}_{1}^{m}+\frac{1}{2}\tilde{\mathbf{a}}_{2}^{m}$
(we denote the mirror as $\mathcal{M}_{10}$ since its normal is the direction
$\tilde{\mathbf{a}}_{1}^{m}$ direction):
$\mathcal{H}^{(s)}(\mathcal{G}_{10}\mathbf{r})=\mathds{1}\otimes\sigma_{x}\mathcal{H}^{(s)}(\mathbf{r})\mathds{1}\otimes\sigma_{x}$
(which is due to the fact that
$U_{1}(\mathcal{G}_{10}\mathbf{r})=U_{1}^{*}(\mathbf{r})$ for
$\alpha=\alpha^{*}$). Therefore, each block forms a QBCP system under periodic
potential with $p4gm$ symmetry. However, it is known that in $p4gm$ lattice,
high symmetry points have multiplicity of 2 (see, for example, BCS
aroyo2011crystallographys,aroyo2006bilbaoIs,aroyo2006bilbaoIIs for a list of
high symmetry points in the unit cell for any space group) in the unit cell.
Therefore, if the components of $\Psi_{\Gamma^{m}}^{(s)}=U\Psi_{\Gamma^{m}}$,
$\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})=(\psi_{\Gamma^{m},1}(\mathbf{r})+\psi_{\Gamma^{m},2}(\mathbf{r}))/\sqrt{2}$
and
$\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})=(-\psi_{\Gamma^{m},1}(\mathbf{r})+\psi_{\Gamma^{m},2}(\mathbf{r}))/\sqrt{2}$,
have zeros, they come in pairs in the unit cell defined by lattice vectors
$\tilde{\mathbf{a}}_{i}^{m}$. Indeed, there are two zeros in each of
$\Psi_{\Gamma^{m},1}^{(s)}(\mathbf{r})$ and
$\Psi_{\Gamma^{m},3}^{(s)}(\mathbf{r})$ at “magic” value of $\alpha$ as can be
seen from Fig. S9(d-e). Because of this, for each diagonal block of
$\mathcal{H}^{(s)}(\mathbf{r})$, once can construct 2 FB WFs polarized to each
sublattice (indeed each block of $\mathcal{H}^{(s)}(\mathbf{r})$ is equivalent
to the Hamiltonian we consider in Fig. S6 where we find 4 FBs at “magic”
$\alpha$). This degeneracy can also be understood from “band sticking” effect
in nonsymmorphic lattices. The two zeros allow for defining two independent
holomorphic functions $f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(1)})$ and
$f_{\mathbf{k}}(z;\mathbf{r}_{0}^{(2)})$; this, together with the construction
of Eqs. (S87) and (S89) for each holomorphic function, give the 2 FB wFs
polarized to each sublattice (considering two sublattices, total 4 FBs) for
the Hamiltonian $\mathcal{H}(\mathbf{r})$ in Eq. (S90). Furthermore, we know (
from Sec. S-IVB) that the total Chern number of the multiple bands polarized
on the same sublattice is $C=1$ when evaluated over Brillouin zone defined by
$\tilde{\mathbf{b}}_{i}^{m}$. However, due to
$\mathcal{H}^{(s)}(\mathbf{r}+\mathbf{a}_{i}^{m})=\sigma_{x}\otimes\mathds{1}\mathcal{H}^{(s)}(\mathbf{r})\sigma_{x}\otimes\mathds{1}$,
the Brillouin zone is now twice large (defined by reciprocal lattice
$\mathbf{b}_{i}^{m}$). This is why the Chern number of two bands polarized on
each sublattice is $C=\pm 2$.

##  S-8 Details of the topological heavy fermion (THF) model shown in Fig. 4
of the main text

As was discussed in the main text, due to the antiunitary particle-hole
symmetry $\mathcal{P}$, any set of bands symmetric about the charge neutrality
point is topological TBGIIBernevigs. 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^{m}$
point in the mBZ, it was shown in song2022magics that hybridization of 2
atomic limit HF bands with 4 topological conduction bands (having nontrivial
winding) at $\Gamma^{m}$ 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 mBZ. We showed an
example of this in Fig. 4 of the main text for the system (space group
$G=p6mm$) with 6 FBs in Fig. 4(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) bradlyn2017topologicals. 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}$ (we use the same notation as Topological
Quantum Chemistry section of Bilbao Crystallography Server (BCS)
aroyo2011crystallographys,aroyo2006bilbaoIs,aroyo2006bilbaoIIs). This and the
fact that the Berry curvature distribution of the 6 FBs being peaked at
$\Gamma^{m}$ (Fig 4(c)) suggest a THF model composed of local orbitals having
band representation $BR$ (we will refer to them as $f$-electrons) and
topological conduction bands (we will refer to them as $c$-bands) with
representation $\Gamma_{6}$. The details of the construction of localized
Wannier functions as well as the single particle THF Hamiltonian are shown
below. This section follows song2022magics.

###  S-8.1 Maximally localized Wannier functions for the $f$-electrons

We start by recalling the basis states $|\mathbf{k}_{0}+\mathbf{k},A/B\rangle$
of the QBCP described in Sec. S-I. In systems with 3-fold rotation symmetry,
QBCP can only appear at the $\Gamma$ point, hence $\mathbf{k}_{0}=\mathbf{0}$
in our example. The basis in the real space then is

$|\mathbf{r},\alpha\rangle=\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{r}}|\mathbf{k},\alpha\rangle=\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i(\mathbf{k}+\mathbf{b}^{m})\cdot\mathbf{r}}|\mathbf{k},\mathbf{b}^{m},\alpha\rangle,\text{,
}\alpha\in\\{A,B\\},$ (S92)

where we broke down the sum over the $\mathbf{k}$ into sum over $\mathbf{k}$
in the mBZ and sum over moiré reciprocal lattice vectors and
$|\mathbf{k},\mathbf{b}^{m},\alpha\rangle\equiv|\mathbf{k}+\mathbf{b}^{m},\alpha\rangle$.
Recall that these basis states have the following transformation properties

$\begin{split}\mathcal{C}_{3z}\\{|\mathbf{r},A\rangle,|\mathbf{r},B\rangle\\}&=\\{|\mathcal{C}_{3z}\mathbf{r},A\rangle,|\mathcal{C}_{3z}\mathbf{r},B\rangle\\}\rho(\mathcal{C}_{3z}),\,\rho(\mathcal{C}_{3z})=\text{Diag}\\{e^{4\pi
i/3},e^{2\pi i/3}\\},\\\
\mathcal{C}_{2z}\\{|\mathbf{r},A\rangle,|\mathbf{r},B\rangle\\}&=\\{|\mathcal{C}_{2z}\mathbf{r},A\rangle,|\mathcal{C}_{2z}\mathbf{r},B\rangle\\}\rho(\mathcal{C}_{2z}),\,\rho(\mathcal{C}_{2z})=\mathds{1},\\\
\mathcal{M}_{x}\\{|\mathbf{r},A\rangle,|\mathbf{r},B\rangle\\}&=\\{|\mathcal{M}_{x}\mathbf{r},A\rangle,|\mathcal{M}_{x}\mathbf{r},B\rangle\\}\rho(\mathcal{M}_{x}),\,\rho(\mathcal{M}_{x})=\sigma_{x},\\\
\mathcal{T}\\{|\mathbf{r},A\rangle,|\mathbf{r},B\rangle\\}&=\\{|\mathbf{r},A\rangle,|\mathbf{r},B\rangle\\}\rho(\mathcal{M}_{x}),\,\rho(\mathcal{T})=\sigma_{x},\\\
T_{\mathbf{R}}\\{|\mathbf{r},A\rangle,|\mathbf{r},B\rangle\\}&=\\{|\mathbf{r}+\mathbf{R},A\rangle,|\mathbf{r}+\mathbf{R},B\rangle\\}\end{split}$
(S93)

where we chose $\rho(\mathcal{C}_{2z})=\mathds{1}$ to specify that the irrep
label of the QBCP is $\Gamma_{5}$ in the notation of BCS (we could have just
as easily chose $\rho(\mathcal{C}_{2z})=-\mathds{1}$, then irrep would have
been $\Gamma_{6}$). Also, here,
$\mathbf{R}=n_{1}\mathbf{a}_{1}^{m}+n_{2}\mathbf{a}_{2}^{m}$ is a moiré
lattice vector.

We want to construct trial Wannier functions that transform as $A_{1}$
($s$-type orbital), $A_{2}$, $E_{1}$ ($p$-type orbitals) and $E_{2}$ ($d$-type
orbitals) representations of $\mathcal{C}_{6v}$ at $1a$ Wyckoff position or
the center of the Wigner-Seitz unit cell. The simplest one of these is the
construction of the $E_{2}$ reps:

$\begin{split}|W^{\prime}_{\mathbf{R},E_{2},d_{x^{2}-y^{2}}+id_{2xy}}\rangle&=\frac{1}{\Omega\sqrt{2\pi\lambda_{0}^{2}}}\int
d^{2}\mathbf{r}e^{-(\mathbf{r}-\mathbf{R})^{2}/2\lambda_{0}^{2}}|\mathbf{r},A\rangle=\frac{1}{{\color[rgb]{0,0,0}\Omega}\sqrt{2\pi\lambda_{0}^{2}}}\int
d^{2}\mathbf{r}e^{-(\mathbf{r}-\mathbf{R})^{2}/2\lambda_{0}^{2}}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i(\mathbf{k}+\mathbf{b}^{m})\cdot\mathbf{r}}|\mathbf{k},\mathbf{b}^{m},A\rangle\\\
&=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}|\mathbf{k},\mathbf{b}^{m},A\rangle,\\\
|W^{\prime}_{\mathbf{R},E_{2},d_{x^{2}-y^{2}}-id_{2xy}}\rangle&=\frac{1}{\Omega\sqrt{2\pi\lambda_{0}^{2}}}\int
d^{2}\mathbf{r}e^{-(\mathbf{r}-\mathbf{R})^{2}/2\lambda_{0}^{2}}|\mathbf{r},B\rangle=\frac{1}{{\color[rgb]{0,0,0}\Omega}\sqrt{2\pi\lambda_{0}^{2}}}\int
d^{2}\mathbf{r}e^{-(\mathbf{r}-\mathbf{R})^{2}/2\lambda_{0}^{2}}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i(\mathbf{k}+\mathbf{b}^{m})\cdot\mathbf{r}}|\mathbf{k},\mathbf{b}^{m},B\rangle\\\
&=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}|\mathbf{k},\mathbf{b}^{m},B\rangle,\end{split}$
(S94)

because the basis states $|\mathbf{r},\alpha\rangle$ transform as $\Gamma_{5}$
rep, which are also $d$-type. This part is exactly the same as constructing
$p_{x}\pm ip_{y}$ orbitals in TBG song2022magics. However, constructing the
other 4 Wannier functions is new in this system compared to TBG. To create the
$E_{1}$ rep or the $p$-orbitals, all we need is to have an extra negative sign
under rotation, that can be done in the following way

$\begin{split}|W^{\prime}_{\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle&=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}ie^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}|\mathbf{k},\mathbf{b}^{m},A\rangle,\\\
|W^{\prime}_{\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle&=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}ie^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}|\mathbf{k},\mathbf{b}^{m},B\rangle,\end{split}$
(S95)

where
$\theta_{\mathbf{k}+\mathbf{b}^{m}}=\text{arg}((k_{x}+b^{m}_{x})+i(k_{y}+b^{m}_{y}))$.
We can easily verify

$\displaystyle\mathcal{C}_{3z}\\{|W^{\prime}_{\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{C}_{3z}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{C}_{3z}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathcal{C}_{3z}\mathbf{k},\mathcal{C}_{3z}\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathcal{C}_{3z}\mathbf{k},\mathcal{C}_{3z}\mathbf{b}^{m},B\rangle\\}\text{Diag}\\{e^{4\pi
i/3},e^{2\pi i/3}\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathcal{C}_{3z}^{-1}\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathcal{C}_{3z}^{-1}\mathbf{k}+\mathcal{C}_{3z}^{-1}\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathcal{C}_{3z}^{-1}(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathcal{C}_{3z}^{-1}(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}\text{Diag}\\{e^{4\pi
i/3},e^{2\pi i/3}\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathcal{C}_{3z}\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i(\theta_{\mathbf{k}+\mathbf{b}^{m}}-2\pi/3)}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i(\theta_{\mathbf{k}+\mathbf{b}^{m}}-2\pi/3)}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}\text{Diag}\\{e^{4\pi
i/3},e^{2\pi i/3}\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathcal{C}_{3z}\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}\text{Diag}\\{e^{4\pi
i/3}e^{-2\pi i/3},e^{2\pi i/3}e^{2\pi i/3}\\}$
$\displaystyle=\\{|W^{\prime}_{\mathcal{C}_{3z}\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathcal{C}_{3z}\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}\text{Diag}\\{e^{2\pi
i/3},e^{4\pi i/3}\\},$ (S96a)
$\displaystyle\mathcal{C}_{2z}\\{|W^{\prime}_{\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{C}_{2z}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{C}_{2z}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathcal{C}_{2z}\mathbf{k},\mathcal{C}_{2z}\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathcal{C}_{2z}\mathbf{k},\mathcal{C}_{2z}\mathbf{b}^{m},B\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathcal{C}_{2z}^{-1}\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathcal{C}_{2z}^{-1}\mathbf{k}+\mathcal{C}_{2z}^{-1}\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathcal{C}_{2z}^{-1}(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathcal{C}_{2z}^{-1}(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathcal{C}_{2z}\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i(\theta_{\mathbf{k}+\mathbf{b}^{m}}-\pi)}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i(\theta_{\mathbf{k}+\mathbf{b}^{m}}-\pi)}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathcal{C}_{2z}\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}\text{Diag}\\{e^{-\pi
i},e^{\pi i}\\}$
$\displaystyle=\\{|W^{\prime}_{\mathcal{C}_{2z}\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathcal{C}_{2z}\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}\text{Diag}\\{-1,-1\\},$
(S96b)
$\displaystyle\mathcal{M}_{x}\\{|W^{\prime}_{\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{M}_{x}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{M}_{x}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathcal{M}_{x}\mathbf{k},\mathcal{M}_{x}\mathbf{b}^{m},B\rangle,e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathcal{M}_{x}\mathbf{k},\mathcal{M}_{x}\mathbf{b}^{m},A\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathcal{M}_{x}^{-1}\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathcal{M}_{x}\mathbf{k}+\mathcal{M}_{x}\mathbf{b}^{m})^{2}}\\{e^{i\theta_{\mathcal{M}_{x}^{-1}(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},B\rangle,e^{-i\theta_{\mathcal{M}_{x}^{-1}(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},A\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathcal{M}_{x}\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{i(\pi-\theta_{\mathbf{k}+\mathbf{b}^{m}})}|\mathbf{k},\mathbf{b}^{m},B\rangle,e^{-i(\pi-\theta_{\mathbf{k}+\mathbf{b}^{m}})}|\mathbf{k},\mathbf{b}^{m},A\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathcal{M}_{x}\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},B\rangle,e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},A\rangle\\}\text{Diag}\\{-1,-1\\}$
$\displaystyle=\\{|W^{\prime}_{\mathcal{M}_{x}\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle,|W^{\prime}_{\mathcal{M}_{x}\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle\\}\text{Diag}\\{-1,-1\\}$
$\displaystyle=\\{|W^{\prime}_{\mathcal{M}_{x}\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathcal{M}_{x}\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}(-\sigma_{x}),$
(S96c)
$\displaystyle\mathcal{T}\\{|W^{\prime}_{\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}$
$\displaystyle=-i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{T}|\mathbf{k},\mathbf{b}^{m},A\rangle,e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\mathcal{T}|\mathbf{k},\mathbf{b}^{m},B\rangle\\}$
$\displaystyle=-i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|-\mathbf{k},-\mathbf{b}^{m},B\rangle,e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|-\mathbf{k},-\mathbf{b}^{m},A\rangle\\}$
$\displaystyle=-i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{i(-\mathbf{k})\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(-\mathbf{k}-\mathbf{b}^{m})^{2}}\\{e^{-i\theta_{-(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},B\rangle,e^{i\theta_{-(\mathbf{k}+\mathbf{b}^{m})}}|\mathbf{k},\mathbf{b}^{m},A\rangle\\}$
$\displaystyle=-i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{-i(\pi+\theta_{\mathbf{k}+\mathbf{b}^{m}})}|\mathbf{k},\mathbf{b}^{m},B\rangle,e^{i(\pi+\theta_{\mathbf{k}+\mathbf{b}^{m}})}|\mathbf{k},\mathbf{b}^{m},A\rangle\\}$
$\displaystyle=i\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\\{e^{-i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},B\rangle,e^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},A\rangle\\}$
$\displaystyle=\\{|W^{\prime}_{\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle\\}\sigma_{x}.$
(S96d)

Similiarly, one can check that the following trial Wannier functions transform
as $A_{1}$ and $A_{2}$ rep of $C_{6v}$

$\begin{split}|W^{\prime}_{\mathbf{R},A_{1}}\rangle&=\frac{\sqrt{\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}(e^{2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},A\rangle+e^{-2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},B\rangle),\\\
|W^{\prime}_{\mathbf{R},A_{2}}\rangle&=-\frac{\sqrt{\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}i(e^{2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},A\rangle-e^{-2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}|\mathbf{k},\mathbf{b}^{m},B\rangle).\end{split}$
(S97)

In the basis
$\\{|W^{\prime}_{\mathbf{R},A_{1}}\rangle,|W^{\prime}_{\mathbf{R},A_{2}}\rangle,|W^{\prime}_{\mathbf{R},E_{1},p_{x}+ip_{y}}\rangle,|W^{\prime}_{\mathbf{R},E_{1},p_{x}-ip_{y}}\rangle,|W^{\prime}_{\mathbf{R},E_{2},d_{x^{2}-y^{2}}+d_{2xy}}\rangle,|W^{\prime}_{\mathbf{R},E_{2},d_{x^{2}-y^{2}}-d_{2xy}}\rangle\\}$
the representations of the symmetries are

$\begin{split}\rho^{f}(\mathcal{C}_{3z})&=\begin{pmatrix}1&0&0&0&0&0\\\
0&1&0&0&0&0\\\ 0&0&e^{2\pi i/3}&0&0&0\\\ 0&0&0&e^{4\pi i/3}&0&0\\\
0&0&0&0&e^{4\pi i/3}&0\\\ 0&0&0&0&0&e^{2\pi i/3}\\\
\end{pmatrix},\rho^{f}(\mathcal{C}_{2z})=\begin{pmatrix}1&0&0&0&0&0\\\
0&1&0&0&0&0\\\ 0&0&-1&0&0&0\\\ 0&0&0&-1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&1\\\
\end{pmatrix},\rho^{f}(\mathcal{M}_{x})=\begin{pmatrix}1&0&0&0&0&0\\\
0&-1&0&0&0&0\\\ 0&0&0&-1&0&0\\\ 0&0&-1&0&0&0\\\ 0&0&0&0&0&1\\\ 0&0&0&0&1&0\\\
\end{pmatrix}\\\ \rho^{f}(\mathcal{T})&=\begin{pmatrix}1&0&0&0&0&0\\\
0&1&0&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&0&0&1\\\ 0&0&0&0&1&0\\\
\end{pmatrix},\rho^{f}(\mathcal{S})=\begin{pmatrix}0&-i&0&0&0&0\\\
i&0&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&-1&0&0\\\ 0&0&0&0&1&0\\\ 0&0&0&0&0&-1\\\
\end{pmatrix},\rho^{f}(\mathcal{P})=\rho^{f}(\mathcal{ST})=\begin{pmatrix}0&-i&0&0&0&0\\\
i&0&0&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&-1&0&0&0\\\ 0&0&0&0&0&1\\\ 0&0&0&0&-1&0\\\
\end{pmatrix},\end{split}$ (S98)

Next we calculate the overlap between these trial Wannier functions and the
energy eigenstates. Denoting the numerically obtained energy eigenstates as
$|\psi_{\mathbf{k},n}\rangle$, we define the overlap matrix as

$\displaystyle A_{n,1}(\mathbf{k})$
$\displaystyle\equiv\langle\psi_{\mathbf{k},n}|W^{\prime}_{\mathbf{0},A_{1}}\rangle=\frac{\sqrt{\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}(e^{2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},A\rangle+e^{-2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},B\rangle)$
$\displaystyle A_{n,2}(\mathbf{k})$
$\displaystyle\equiv\langle\psi_{\mathbf{k},n}|W^{\prime}_{\mathbf{0},A_{2}}\rangle=\frac{\sqrt{\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}i(e^{2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},A\rangle-e^{-2i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},B\rangle)$
$\displaystyle A_{n,3}(\mathbf{k})$
$\displaystyle\equiv\langle\psi_{\mathbf{k},n}|W^{\prime}_{\mathbf{0},E_{1},p_{x}+ip_{y}}\rangle=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}ie^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},A\rangle$
$\displaystyle A_{n,4}(\mathbf{k})$
$\displaystyle\equiv\langle\psi_{\mathbf{k},n}|W^{\prime}_{\mathbf{0},E_{1},p_{x}-ip_{y}}\rangle=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}ie^{i\theta_{\mathbf{k}+\mathbf{b}^{m}}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},B\rangle$
$\displaystyle A_{n,5}(\mathbf{k})$
$\displaystyle\equiv\langle\psi_{\mathbf{k},n}|W^{\prime}_{\mathbf{0},E_{2},d_{x^{2}-y^{2}}+id_{2xy}}\rangle=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},A\rangle$
$\displaystyle A_{n,6}(\mathbf{k})$
$\displaystyle\equiv\langle\psi_{\mathbf{k},n}|W^{\prime}_{\mathbf{0},E_{2},d_{x^{2}-y^{2}}-id_{2xy}}\rangle=\frac{\sqrt{2\pi\lambda_{0}^{2}}}{\Omega}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}+\mathbf{b}^{m})^{2}}\langle\psi_{\mathbf{k},n}|\mathbf{k},\mathbf{b}^{m},B\rangle$
(S99)

Figure S10: Overlap $|A_{n,\alpha}(\mathbf{k})|$ (Eq. (S-8.1)) between the
trial Wannier functions and the energy eigenstates are plotted in red circles
on the energy bands. The top row shows all the 8 lowest bands, whereas the
bottom row zooms into the the lowest 6 bands for better visualization. (a-b)
show the overlap $|A_{n,5}(\mathbf{k})|=|A_{n,6}(\mathbf{k})|$ (the equality
is due to chiral symmetry $\mathcal{S}$) of the Wannier functions
corresponding to $E_{2}$ ($d$-type) irrep. (c-d) show the overlap
$|A_{n,3}(\mathbf{k})|=|A_{n,4}(\mathbf{k})|$ (the equality is due to chiral
symmetry $\mathcal{S}$) of the Wannier functions corresponding to $E_{1}$
($p$-type) irrep. (e-f) show the overlap $|A_{n,1}(\mathbf{k})|$ of the
Wannier function corresponding to $A_{1}$ ($s$-type) irrep. (g-h) show the
overlap $|A_{n,2}(\mathbf{k})|$ of the Wannier function corresponding to
$A_{2}$ irrep.

We plot these overlap functions for $n=\pm 1,\pm 2,\pm 3,\pm 4$ (the bands are
numbered away from the charge neutrality as $\pm 1,\pm 2,\dots$) in Fig. S10.
Clearly the Wannier functions are completely supported by the middle six bands
everywhere in the mBZ except at $\Gamma^{m}$, where the $A_{1}$ and $A_{2}$
type Wannier functions are supported by the lowest higher bands at
$\Gamma^{m}$. We feed the overlap matrix $A_{n,\alpha}(\mathbf{k})$ ($n=\pm
1,\dots,\pm 4,\alpha=1,\dots,6$) into the machinary of Wannier90
marzari1997maximallys,souza2001maximallys,pizzi2020wannier90s to construct the
Maximally localized Wannier functions (MLWFs). We chose $\lambda_{0}=a^{m}/10$
for the numerical calculation, and used a $20\times 20$ grid to discretize the
mBZ, and chose an energy window such that only the lowest 8 bands fall inside
the window for the disentanglement and Wannierization procedure. Wannier90
returns MLWFs in the plane wave basis
$|\mathbf{k},\mathbf{b}^{m},\beta\rangle$ ($\beta=A,B$) as

$|W_{\mathbf{R},\alpha}\rangle=\frac{1}{\Omega\sqrt{N}}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}|\mathbf{k},\mathbf{b}^{m},\beta\rangle
e^{-i\mathbf{k}\cdot\mathbf{R}}\tilde{v}_{\mathbf{b}^{m}\beta,\alpha}(\mathbf{k}),$
(S100)

where $N$ is the number of moiré unit cell. We can write the MLWFs in the real
space basis

$w_{\beta,\alpha}(\mathbf{r}-\mathbf{R})=\langle\mathbf{r},\beta|W_{\mathbf{R},\alpha}\rangle=\frac{1}{\Omega\sqrt{N}}\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m}}e^{i(\mathbf{k}+\mathbf{b}^{m})\cdot(\mathbf{r}-\mathbf{R})}\tilde{v}_{\mathbf{b}^{m}\beta,\alpha}(\mathbf{k}).$
(S101)

The density plots of $\sum_{\beta}|w_{\beta,\alpha}|^{2}$ are shown in Fig.
4(d) of the main text. Note that since the $|W_{\mathbf{R},A_{1}}\rangle$ and
$|W_{\mathbf{R},A_{2}}\rangle$ are related by chiral symmetry, their density
plot look the same. Also, since $p_{x}\pm ip_{y}$ are related by
$\mathcal{T}$, $\sum_{\beta}|w_{\beta,3}|^{2}=\sum_{\beta}|w_{\beta,4}|^{2}$;
this is why we only plotted $\sum_{\beta}|w_{\beta,3}|^{2}$. Similarly,
$\sum_{\beta}|w_{\beta,5}|^{2}=\sum_{\beta}|w_{\beta,6}|^{2}$ since
$d_{x^{2}-y^{2}}\pm id_{2xy}$ are related by $\mathcal{T}$, and we plotted
$\sum_{\beta}|w_{\beta,5}|^{2}$. Clearly, they are well localized within the
unit cell.

A creation operators of the Wannier states can be introduced as

$\begin{split}f_{\mathbf{R},\alpha}^{\dagger}&=\sum_{\beta}\int
d^{2}\mathbf{r}\langle\mathbf{r},\beta|W_{\mathbf{R},\alpha}\rangle
c_{\beta}^{\dagger}(\mathbf{r})=\sum_{\beta}\int
d^{2}\mathbf{r}\,w_{\beta,\alpha}(\mathbf{r}-\mathbf{R})c_{\beta}^{\dagger}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum_{\mathbf{k},\mathbf{b}^{m},\beta}e^{-i\mathbf{k}\cdot\mathbf{R}}\tilde{v}_{\mathbf{b}^{m}\beta,\alpha}(\mathbf{k})c_{\mathbf{k},\mathbf{b}^{m},\beta}^{\dagger},\\\
f_{\mathbf{k},\alpha}^{\dagger}&=\sum_{\mathbf{b}^{m},\beta}\tilde{v}_{\mathbf{b}^{m}\beta,\alpha}(\mathbf{k})c_{\mathbf{k},\mathbf{b}^{m},\beta}^{\dagger},\end{split}$
(S102)

where $c_{\mathbf{k},\mathbf{b}^{m},\beta}^{\dagger}$ is the creation operator
of the plane-wave state $|\mathbf{k},\mathbf{b}^{m},\beta\rangle$.

###  S-8.2 The $c$-electrons

The Wannier functions span
$\Gamma_{1}\oplus\Gamma_{2}\oplus\Gamma_{5}\oplus\Gamma_{6}$ representations
among $\Gamma_{1}\oplus\Gamma_{2}\oplus\Gamma_{5}\oplus 2\Gamma_{6}$
representation formed by the middle 8 bands. However, the middle six bands
actually have representations $\Gamma_{5}\oplus 2\Gamma_{6}$. Thus to get the
correct band structure at the $\Gamma^{m}$ point, we need two additional
degrees of freedom that form $\Gamma_{6}$ representation. Following
song2022magics, they can be formally written as

$c_{\mathbf{k},a}^{\dagger}=\sum_{\mathbf{b}^{m},\beta}\tilde{u}_{\mathbf{b}^{m}\beta,a}(\mathbf{k})c_{\mathbf{k},\mathbf{b}^{m},\beta}^{\dagger},$
(S103)

where $\tilde{u}_{\mathbf{b}^{m}\beta,a}(\mathbf{k})$ will be determined
below. Note that in the plane wave basis the single layer QBCP Hamiltonian in
Eq. (S17), can be written as

$\begin{split}\hat{H}&=\sum_{\mathbf{k}\in\text{mBZ}}\sum_{\mathbf{b}^{m},{\mathbf{b}^{m}}^{\prime}}\sum_{\alpha,\beta\in\\{A,B\\}}[h_{\mathbf{b}^{m},{\mathbf{b}^{m}}^{\prime}}(\mathbf{k})]_{\alpha,\beta}c_{\mathbf{k},\mathbf{b}^{m},\alpha}^{\dagger}c_{\mathbf{k},{\mathbf{b}^{m}}^{\prime},\beta},\\\
[h_{\mathbf{b}^{m},{\mathbf{b}^{m}}^{\prime}}(\mathbf{k})]&=\begin{bmatrix}0&(k^{*}+{b^{m}}^{*})^{2}\delta_{\mathbf{b}^{m},{\mathbf{b}^{m}}^{\prime}}+\tilde{A}^{*}(\mathbf{b}^{m}-{\mathbf{b}^{m}}^{\prime})\\\
(k+b^{m})^{2}\delta_{\mathbf{b}^{m},{\mathbf{b}^{m}}^{\prime}}+\tilde{A}(\mathbf{b}^{m}-{\mathbf{b}^{m}}^{\prime})&0\end{bmatrix},\end{split}$
(S104)

where $\tilde{A}(\mathbf{b}^{m})$ are the Fourier components of the periodic
field $\mathcal{D}_{U}(\mathbf{r};\mathbf{\alpha})\equiv\tilde{A}(\mathbf{r})$
in Eq. (S17), $k=k_{x}+ik_{y}$ and $b^{m}=b^{m}_{x}+ib^{m}_{y}$. We can
diagonalize the Hamiltonian

$\sum_{{\mathbf{b}^{m}}^{\prime},\beta}[h_{\mathbf{b}^{m},{\mathbf{b}^{m}}^{\prime}}(\mathbf{k})]_{\alpha,\beta}u_{\mathbf{k},{\mathbf{b}^{m}}^{\prime},\beta,n}=\epsilon_{n}(\mathbf{k})u_{\mathbf{k},{\mathbf{b}^{m}},\alpha,n},$
(S105)

where $\\{u_{\mathbf{k},n}\\}$ is the $n$-th eigenvector of matrix
$[h(\mathbf{k})]$ with eigenvalue $\epsilon_{n}(\mathbf{k})$. We denote the
eigenvalues as
$\dots\leq\epsilon_{-2}(\mathbf{k})\leq\epsilon_{-1}(\mathbf{k})\leq\epsilon_{1}(\mathbf{k})\leq\epsilon_{2}(\mathbf{k})\leq\dots$,
where $\epsilon_{-1}(\mathbf{k})$ and $\epsilon_{1}(\mathbf{k})$ are the
eigenvalues with lowest magnitude. With this, we can define the projector to
the lowest 8 bands

$P_{{\mathbf{b}^{m}}\alpha,{\mathbf{b}^{m}}^{\prime}\beta}(\mathbf{k})=\sum_{n=\pm
1,\pm 2,\pm 3,\pm
4}u_{\mathbf{k},\mathbf{b}^{m},\alpha,n}u_{\mathbf{k},{\mathbf{b}^{m}}^{\prime},\beta,n}^{*}.$
(S106)

On the other hand, the projector to the 6 Wannier states are

$Q_{{\mathbf{b}^{m}}\alpha,{\mathbf{b}^{m}}^{\prime}\beta}(\mathbf{k})=\sum_{\gamma=1,\dots,6}\tilde{v}_{\mathbf{b}^{m}\alpha,\gamma}(\mathbf{k})\tilde{v}_{{\mathbf{b}^{m}}^{\prime}\beta,\gamma}^{*}(\mathbf{k}).$
(S107)

Since, by construction, the Wannier states are linear combinations of the 8
lowest bands, we have $P(\mathbf{k})Q(\mathbf{k})P(\mathbf{k})=Q(\mathbf{k})$.
Then, the projector to the remaining states is given by
$P(\mathbf{k})-Q(\mathbf{k})$. The eigenvectors of
$P(\mathbf{k})-Q(\mathbf{k})$ with eigenvalue 1 are
$\tilde{u}_{\mathbf{b}^{m}\beta,a}(\mathbf{k})$. We fix the gauge of these two
vectors by requiring the representations of the symmetries to be the following

$\rho^{c}(\mathcal{C}_{3z})=\begin{pmatrix}e^{4\pi i/3}&0\\\ 0&e^{2\pi
i/3}\end{pmatrix},\rho^{c}(\mathcal{C}_{2z})=\begin{pmatrix}-1&0\\\
0&-1\end{pmatrix},\rho^{c}(\mathcal{M}_{x})=\begin{pmatrix}0&-1\\\
-1&0\end{pmatrix},\rho^{c}(\mathcal{T})=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\rho^{c}(\mathcal{S})=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}.$ (S108)

###  S-8.3 The single particle Hamiltonian

After obtaining the $f$-electron and the $c$-band basis, we are at a position
to obtain the single-particle effective topological heavy fermion Hamiltonian.
To this end, we define two matrices

$\begin{split}[U(\mathbf{k})]&=[\\{\tilde{v}_{1}(\mathbf{k})\\},\\{\tilde{v}_{2}(\mathbf{k})\\},\\{\tilde{v}_{3}(\mathbf{k})\\},\\{\tilde{v}_{4}(\mathbf{k})\\},\\{\tilde{v}_{5}(\mathbf{k})\\},\\{\tilde{v}_{6}(\mathbf{k})\\},\\{\tilde{u}_{1}(\mathbf{k})\\},\\{\tilde{u}_{2}(\mathbf{k})\\}],\\\
[U_{C}(\mathbf{k})]&=[\dots,\\{u_{\mathbf{k},-7}\\},\\{u_{\mathbf{k},-6}\\},\\{u_{\mathbf{k},-5}\\},\\{u_{\mathbf{k},5}\\},\\{u_{\mathbf{k},6}\\},\\{u_{\mathbf{k},7}\\},\dots]\end{split}$
(S109)

Then, we project the hamiltonian matrix $[h(\mathbf{k})]$ into the lowest 8
bands for small $|\mathbf{k}|$ (accurate to the second order in expansion
w.r.t. $|\mathbf{k}|$) in the following way

$\begin{split}[h_{P}(\mathbf{k})]&=[\tilde{h}(\mathbf{k})]-[C(\mathbf{k})]^{\dagger}[\tilde{h}_{C}(\mathbf{k})][C(\mathbf{k})]=\left[\begin{array}[]{c|c}H^{(f)}(\mathbf{k})&H^{(fc)}(\mathbf{k})\\\
\hline\cr H^{(cf)}(\mathbf{k})&H^{(c)}(\mathbf{k})\end{array}\right],\\\
[\tilde{h}(\mathbf{k})]&=[U(\mathbf{0})]^{\dagger}[h(\mathbf{k})][U(\mathbf{0})],\\\
[\tilde{h}_{C}(\mathbf{k})]&=[U_{C}(\mathbf{0})]^{\dagger}[h(\mathbf{k})][U_{C}(\mathbf{0})],\\\
[C(\mathbf{k})]&=[U_{C}(\mathbf{0})]^{\dagger}[h(\mathbf{k})][U(\mathbf{0})].\end{split}$
(S110)

Due to the choice of the gauge for the $c$-bands in the previous subsection,
$H^{(c)}(\mathbf{k})$ has the form

$H^{(c)}(\mathbf{k})=\begin{pmatrix}0&c_{1}(k^{*})^{2}\\\
c_{1}k^{2}&0\end{pmatrix},$ (S111)

where $k=k_{x}+ik_{y}$ and $c_{1}$ is a real constant. For the $f$-electrons,
since they are localized, we get

$H^{(f)}(\mathbf{k})\approx\begin{pmatrix}m\sigma_{z}&\mathbf{0}_{2\times
4}\\\ \mathbf{0}_{4\times 2}&\mathbf{0}_{4\times 4}\end{pmatrix},$ (S112)

where $2|m|$ sets the gap between the $\Gamma_{1}$ and $\Gamma_{2}$ reps (see
Fig. 4(b)). The coupling between the $c$-bands and the $f$-electrons (keeping
only the lowest order terms)

$H^{(cf)}(\mathbf{k})\approx\begin{pmatrix}-ic_{2}(k_{x}-ik_{y})&c_{2}(k_{x}-ik_{y})&0&\gamma&0&ic_{3}(k_{x}+ik_{y})\\\
-ic_{2}(k_{x}+ik_{y})&-c_{2}(k_{x}+ik_{y})&\gamma&0&ic_{3}(k_{x}-ik_{y})&0\\\
\end{pmatrix},$ (S113)

where $\gamma$, $c_{2}$ and $c_{3}$ are real constants. For the Flat bands, we
find $c_{3}$ to be negligible. Furthermore, $2|\gamma|$ sets the gap between
the two $\Gamma_{6}$ reps. Since the $f$-electrons are localized, the integral
$\langle\mathbf{k},a|\hat{H}|W_{\mathbf{0},\alpha}\rangle$ ($\hat{H}$ is the
QBCP Hamiltonian in Eq. (S104)) should decay exponentially with
$|\mathbf{k}|$. For simplicity, we choose the decay factor to be
$e^{-\lambda^{2}|\mathbf{k}|^{2}/2}$ with $\lambda$ being the spread of the
Wannier functions; this is same as what was done in the case of TBG
song2022magics. Lastly, since at large $\mathbf{k}$ has huge kinetic energy,
we put a cutoff $|\mathbf{k}|<\Lambda_{c}$ for the $c$-electron momentum. All
these considerations together give the single particle THF Hamiltonian (Eq.
(8) of main text)

$\begin{split}&\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})+\\\
&\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),\\\
&H^{(c)}(\mathbf{k})\approx
c_{1}(k_{x}^{2}-k_{y}^{2})\sigma_{x}-2c_{1}k_{x}k_{y}\sigma_{y},\\\
&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},\\\
&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}.\end{split}$ (S114)

apsrev4-1 ref.bib