$w_{0}/w_{1}$ | $t_{0}$ | $\gamma$ | $M$ | $v_{\star}$ | $v_{\star}^{\prime}$ | $v_{\star}^{\prime\prime}$
---|---|---|---|---|---|---
1.0 | 0.01512 | 4.086 | 4.727 | -3.829 | 1.492 | -0.01237
0.9 | 0.1099 | -10.21 | 4.194 | -4.089 | 1.579 | -0.01680
0.8 | 0.1497 | -24.75 | 3.697 | -4.303 | 1.623 | -0.03320
0.7 | 0.2276 | -39.11 | 3.248 | -4.483 | 1.624 | -0.04012
0.6 | 0.2789 | -52.91 | 2.854 | -4.637 | 1.580 | -0.04009
0.5 | 0.3027 | -65.78 | 2.518 | -4.774 | 1.483 | -0.03509
Units | meV | $\rm eV\cdot\mathring{A}$
Table S4: Parameters of the single-particle Hamiltonian of the topological
heavy fermion model for different $w_{0}/w_{1}$. Since
$v_{\star}^{\prime\prime}$ is extremely small compared to other parameters, in
all the single-particle and many-body calculations we have set
$v_{\star}^{\prime\prime}=0$. Parameters of the BM model used to obtain this
table are $v_{F}=5.944\mathrm{eV\cdot\mathring{A}}$,
$|\mathbf{K}|=1.703\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$,
$\theta=1.05^{\circ}$.

We now derive the form of $H^{(c,\eta)}$ in the valley $\eta=+$ to linear
order of $\mathbf{k}$. The effective Hamiltonian in the other valley can be
obtained through the time-reversal symmetry, i.e.,
$H^{(c,-)}(\mathbf{k})=H^{(c,+)*}(-\mathbf{k})$. For simplicity, in the
following we introduce an additional set of Pauli (identity) matrices
$\zeta_{0,x,y,z}$ to distinguish the $\Gamma_{3}$ space ($a=1,2$) and the
$\Gamma_{1}\oplus\Gamma_{2}$ space ($a=3,4$). The $a=1,2$ and $a=3,4$ span the
subspaces of $\zeta_{z}=1$ and $\zeta_{z}=-1$, respectively. Due to Sections
S2.2 and S65, the symmetry operators in the valley $\eta=+$ are given by

$C_{3z}=e^{i\frac{2\pi}{3}\sigma_{z}}\oplus\sigma_{0},\qquad
C_{2x}=\zeta_{0}\sigma_{x},\qquad C_{2z}T=\zeta_{0}\sigma_{x}K,\qquad
P=-i\zeta_{0}\sigma_{z}\ ,$ (S72)

where $\zeta_{0}\sigma_{x}\equiv\sigma_{x}\oplus\sigma_{x}$,
$\zeta_{0}\sigma_{z}\equiv\sigma_{z}\oplus\sigma_{z}$. We first look at the
subspace $\Gamma_{3}$. According to $C_{2z}T$, only $\sigma_{0,x,y}$ terms are
allowed in the Hamiltonian. According to $P$, $\sigma_{0}$ is odd in $k$ while
$\sigma_{x,y}$ are even (constant). $\sigma_{0}$ is invariant under $C_{3z}$,
hence the coefficient before $\sigma_{0}$ must be invariant under $C_{3z}$.
The lowest odd order polynomial satisfying this condition is in third order of
$k$. Thus the $\sigma_{0}$ term vanishes to linear order of $k$. Since
$\sigma_{x,y}$ are not invariant under $C_{3z}$, their (constant) coefficients
must vanish to linear order of $\mathbf{k}$. (One may show that the lowest
order coefficients of $\sigma_{x,y}$ terms are quadratic in $\mathbf{k}$.)

We then look at the subspace $\Gamma_{1}\oplus\Gamma_{2}$. According to
$C_{2z}T$, only $\sigma_{0,x,y}$ terms are allowed. According to $P$,
$\sigma_{0}$ is odd in $k$ while $\sigma_{x,y}$ are even (constant) in $k$.
$\sigma_{0}$ is invariant under $C_{3z}$, hence the coefficient before
$\sigma_{0}$ must be invariant under $C_{3z}$. The lowest odd order polynomial
satisfying this condition is in third order of $k$. Thus the $\sigma_{0}$ term
vanishes to linear order of $k$. In the subspace $\Gamma_{1}\oplus\Gamma_{2}$,
both $\sigma_{x,y}$ are invariant under $C_{3z}$, hence they are allowed by
$C_{3z}$. However $\sigma_{y}$ is forbidden by the $C_{2x}$ symmetry. Thus the
Hamiltonian in the $\Gamma_{1}\oplus\Gamma_{2}$ subspace must have the form
$M\sigma_{x}$.

Thirdly we look at the coupling between $\Gamma_{3}$ and
$\Gamma_{1}\oplus\Gamma_{2}$. We only need to consider the off-diagonal terms
$\zeta_{x,y}\sigma_{0,x,y,z}$. According to $C_{2z}T$, only
$\zeta_{x}\sigma_{0,x,y}$ and $\zeta_{y}\sigma_{z}$ terms are allowed.
According to $P$, the coefficients of $\zeta_{x}\sigma_{x,y}$ are even
(constant) in $k$ while $\zeta_{x}\sigma_{0}$, $\zeta_{y}\sigma_{z}$ are odd
in $k$. Since $\zeta_{x}\sigma_{x,y}$ are not invariant under $C_{3z}$, their
constant coefficients must vanish. We then study the odd order terms
$\zeta_{x}\sigma_{0}$, $\zeta_{y}\sigma_{z}$. Under $C_{3z}$ they transform as
$C_{3z}\zeta_{x}\sigma_{0}C_{3z}^{-1}=\zeta_{x}\sigma_{0}\cos\frac{2\pi}{3}-\zeta_{y}\sigma_{z}\sin\frac{2\pi}{3}$,
$C_{3z}\zeta_{y}\sigma_{z}\sigma_{0}C_{3z}^{-1}=\zeta_{y}\sigma_{z}\cos\frac{2\pi}{3}+\zeta_{x}\sigma_{0}\sin\frac{2\pi}{3}$.
Thus the $\zeta_{x}\sigma_{0},\zeta_{y}\sigma_{z}$ terms can be chosen as
$v_{\star}k_{x}\zeta_{x}\sigma_{0}-v_{\star}k_{y}\zeta_{y}\sigma_{z}$, which
is also symmetric under $C_{2x}$.

In summary, the effective Hamiltonian $H^{(c,\eta)}$ has the form in Eq. S70.

In Ref. Song  _et al._ (2021) the authors of the present work have shown that
the single-particle Hamiltonian in each valley has a symmetry anomaly of the
$C_{2z}T$ and the particle-hole $P$ symmetries. The anomaly is reflected as
unavoidable $4n+2$ ($n\in\mathbb{N}$) Dirac cones at the Fermi-level. Since
the local orbitals are topologically trivial, the anomaly must be reflected in
the band structure of the conduction bands, meaning that, if we turn off the
coupling between local orbitals and conduction bands, the conduction bands
(given by Eq. S70) must carry the symmetry anomaly by itself and hence be
gapless. The spectrum of Eq. S70 is given by

$\pm\frac{M}{2}\pm\sqrt{\frac{M^{2}}{4}+v_{\star}^{2}\mathbf{k}^{2}}$ (S73)

which has a quadratic touching at $\mathbf{k}=0$. To relate the quadratic
touching to the symmetry anomaly, we consider adding a $C_{3z}$-breaking but
$C_{2z}T$\- and $P$-preserving term $D$ into the Hamiltonian

$\begin{pmatrix}D\sigma_{x}&v_{\star}(k_{x}\sigma_{0}+ik_{y}\sigma_{z})\\\
v_{\star}(k_{x}\sigma_{0}-ik_{y}\sigma_{z})&M\sigma_{x}\end{pmatrix}\ .$ (S74)

The quadratic touching is now split into two linear Dirac points at
$k_{x}=\pm\frac{1}{v_{\star}}\sqrt{DM}$, $k_{y}=0$ ($k_{x}=0$,
$k_{y}=\pm\frac{1}{v_{\star}}\sqrt{DM}$) if $DM>0$ ($DM<0$). The two Dirac
points have the same chirality because they are related by $P$ Song  _et al._
(2021). Thus the quadratic touching is topologically equivalent to two Dirac
points with the same chirality. According to the theorem proven in Ref. Song
_et al._ (2021), the presence of $4n+2$ (here $n=0$) Dirac points is an
equivalent condition to the symmetry anomaly.

###  S2.4 Coupling between topological conduction bands and the local orbitals

We now consider the coupling between the conduction band states and the local
orbitals

${\hat{H}}_{cf}=\sum_{\eta
s}\sum_{a\alpha}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\mathbf{R}}\langle\mathbf{k}a\eta
s|H|W_{\mathbf{R}\alpha\eta s}\rangle c_{\mathbf{k}a\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}+h.c.=\sum_{\eta
s}\sum_{a\alpha}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\mathbf{R}}\langle\mathbf{k}a\eta
s|H|W_{0\alpha\eta s}\rangle e^{-i\mathbf{k}\cdot\mathbf{R}}c_{\mathbf{k}a\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}+h.c.,$ (S75)

where $|W_{\mathbf{R}\alpha\eta s}\rangle$ ($\alpha=1,2$) and
$|\mathbf{k}a\eta s\rangle=c_{\mathbf{k},a,\eta,s}^{\dagger}|0\rangle$
($a=1,2,3,4$) are the single particle states of the local orbitals and the
conduction bands, respectively, $H$ is the BM Hamiltonian (Eq. S5), and
$\Lambda_{c}$ is the cutoff for the conduction bands. Since the Wannier
functions of the local orbitals are localized, the integral
$\langle\mathbf{k}a\eta s|H|W_{0\alpha\eta s}\rangle$ must decay exponentially
with $\mathbf{k}$. However, the particular form of the decay is complicated
due to two reasons. First, as discussed in Section S2.1, the Wannier functions
have two independent components (Eqs. S53 and S54). Thus the overlap between
the conduction band states and the two components will decay differently.
Second, each conduction band state is not a simple plane-wave but a linear
combination of plane-waves ($\mathbf{Q}$-vectors) with different wave-vectors
(Eq. S59), hence the overlap for different plane-wave components will decay
differently. However, the role played by the decay is simply a truncation of
the couplings between the local orbitals and high energy states. We claim that
the particular form of the truncation is not relevant in the low energy
physics. Thus, for simplicity here we assume that the overlap decays as
$e^{-|\mathbf{k}|\lambda^{2}/2}$ with $\lambda$ being the square root of the
spread of the Wannier function (Eq. S50). We hence parameterize the integral
as

$\langle\mathbf{k}a\eta s|H|W_{0\alpha\eta
s}\rangle=\frac{1}{\sqrt{N}}e^{-|\mathbf{k}|^{2}\lambda^{2}/2}H^{(cf,\eta)}_{a\alpha}(\mathbf{k})\
.$ (S76)

Applying the Fourier transformation (Eq. S57) of the local orbitals, we can
rewrite the coupling Hamiltonian as

${\hat{H}}_{cf}=\sum_{\eta
s}\sum_{a\alpha}\sum_{|\mathbf{k}|<\Lambda_{c}}e^{-|\mathbf{k}|^{2}\lambda^{2}/2}H^{(cf,\eta)}_{a\alpha}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}f_{\mathbf{k}\alpha\eta s}+h.c.$ (S77)

For small $\mathbf{k}$, the factor $e^{-|\mathbf{k}|^{2}\lambda^{2}/2}$ is
negligible and we can obtain $H^{(cf,\eta)}_{a\alpha}(\mathbf{k})$ as the
k$\cdot$p Hamiltonian

$H^{(cf,\eta)}_{a\alpha}(\mathbf{k})=\langle\widetilde{u}_{a}^{(\eta)}(0)|h^{(\eta)}(\mathbf{k})|\widetilde{v}_{\alpha}^{(\eta)}(0)\rangle,$
(S78)

where $h^{(\eta)}(\mathbf{k})$ is the BM model (Eq. S5),
$|\widetilde{v}_{\alpha}^{(\eta)}(0)\rangle$ are the Bloch state of the
Wannier functions (Eq. S42), and $|\widetilde{u}_{a}^{(\eta)}(0)\rangle$ are
the Bloch states of the conduction bands (Eq. S59). As explained after Eq.
S69, we always choose $|\widetilde{v}_{\alpha}^{(\eta)}(0)\rangle$ and
$|\widetilde{u}_{a}^{(\eta)}(0)\rangle$ rather than
$|\widetilde{v}_{\alpha}^{(\eta)}(\mathbf{k})\rangle$ and
$|\widetilde{u}_{a}^{(\eta)}(\mathbf{k})\rangle$ as the basis for the
k$\cdot$p expansion. Since $h^{(\eta)}(\mathbf{k})$ only has zeroth and linear
order terms in $\mathbf{k}$, $H^{(cf,\eta)}_{a\alpha}(\mathbf{k})$ also only
has constant and linear terms. As proved in the next paragraph, we find the
coupling Hamiltonian has the form

$H^{(cf,\eta)}=\begin{pmatrix}\gamma\sigma_{0}+v_{\star}^{\prime}(\eta
k_{x}\sigma_{x}+k_{y}\sigma_{y})\\\ v_{\star}^{\prime\prime}(\eta
k_{x}\sigma_{x}-k_{y}\sigma_{y})\end{pmatrix}\ .$ (S79)

The first two-by-two block is the coupling between the local orbitals (forming
$\Gamma_{3}$) and the $\Gamma_{3}$ states of the conduction bands and hence
allows a $\gamma\sigma_{0}$ term, the second two-by-two block the coupling
between the local orbitals and the $\Gamma_{1}\oplus\Gamma_{2}$ states of the
conduction bands and vanishes at $\mathbf{k}=0$. The $\Gamma_{3}$
representation carries angular momenta $L=\pm 1$ and the
$\Gamma_{1}\oplus\Gamma_{2}$ representation carries the angular momentum
$L=0$. The above equation respects the angular momentum conservation up to
$\Delta L=0$ mod 3 because the system only has $C_{3z}$ symmetry. The
parameters obtained from Eq. S78 are tabulated in Table S4. We find that
$v_{\star}^{\prime\prime}$ is very small compared to other velocities. We
hence will omit this term.

Here we give one more justification for the choice of the damping factor
$e^{-|\mathbf{k}|^{2}\lambda^{2}/2}$. First, around $\mathbf{k}=0$, the
damping factor $e^{-|\mathbf{k}|^{2}\lambda^{2}/2}\sim 1$ and works well
because the $H^{(cf,\eta)}(\mathbf{k})$ matrix is by construction (Eq. S78)
accurate at $\mathbf{k}=0$. Second, at large $\mathbf{k}$, the actual damping
factor may be complicated due to the reasons given before Eq. S76, but the
$f-c$ coupling is weak compared to the energy separation of $f$\- and
$c$-bands. Thus, the actual form of damping factor is not important. Our
approximation should have at least captured the main feature of the damping
behavior - it assumes the Wannier function as a single exponential function
and the conduction bands as monochromatic plane waves, which are reasonable
due to the localized and delocalized natures of $f$\- and $c$-electrons. As
will be discussed in Section S2.5 and shown in Fig. S2, the choice

Now we use the symmetry principle to derive the form of
$H^{(cf,+)}(\mathbf{k})$. $H^{(cf,-)}(\mathbf{k})$ is given by
$H^{(cf,+)*}(\mathbf{k})$ due to the time-reversal symmetry. We only keep
zeroth and linear order terms in $\mathbf{k}$. The symmetry operators on the
single-particle states of the conduction bands ($c$) are given in Eq. S72.
According to Eq. S36, the symmetry operators on the single-particle states of
the local orbitals ($f$) are

$C_{3z}=e^{i\frac{2\pi}{3}\sigma_{z}},\qquad C_{2x}=\sigma_{x},\qquad
C_{2z}T=\sigma_{x}K,\qquad P=i\sigma_{z}\ .$ (S80)

We first look at the coupling between the $\Gamma_{3}$ states ($a=1,2$) and
the local orbitals. According to $C_{2z}T$, only $\sigma_{0,x,y}$ and
$i\sigma_{z}$ are allowed. According to $P$, the Hamiltonian satisfy
$(-i\sigma_{z})H^{(cf,\eta)}(-\mathbf{k})(i\sigma_{z})^{\dagger}=-H^{(cf,\eta)}(\mathbf{k})$,
where $-i\sigma_{z}$ and $i\sigma_{z}$ are the representation matrix of $P$
formed by the conduction bands and local orbitals, respectively. Thus the
$\sigma_{x,y}$ terms must be odd in $k$ and $\sigma_{0}$, $i\sigma_{z}$ terms
must be even in $k$. The even order terms are constants to first order of $k$,
i.e., $\gamma\sigma_{0}+i\gamma^{\prime}\sigma_{z}$. $\gamma$ respects both
$C_{3z}$ and $C_{2x}$ symmetries. $\gamma^{\prime}$ respects $C_{3z}$ but
breaks $C_{2x}$, hence the $\gamma^{\prime}$ term is forbidden. The two odd
order terms $\sigma_{x,y}$ transform under $C_{3z}$ as
$C_{3z}\sigma_{x}C_{3z}^{-1}=\sigma_{x}\cos\frac{2\pi}{3}+\sigma_{y}\sin\frac{2\pi}{3}$,
$C_{3z}\sigma_{y}C_{3z}^{-1}=\sigma_{y}\cos\frac{2\pi}{3}-\sigma_{x}\sin\frac{2\pi}{3}$,
thus the odd terms can be chosen as
$v_{\star}^{\prime}(k_{x}\sigma_{x}+k_{y}\sigma_{y})$, which is also symmetric
under $C_{2x}$. Thus the coupling between the $\Gamma_{3}$ states in the QDP
and the local orbitals takes the form
$\gamma\sigma_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}+k_{y}\sigma_{y})$. We
then consider couplings between the $\Gamma_{1}\oplus\Gamma_{2}$ states
($a=3,4$) in the conduction bands and the local orbitals. According to
$C_{2z}T$, only $\sigma_{0,x,y}$ and $i\sigma_{z}$ are allowed. According to
$P$, $\sigma_{x,y}$ terms must be odd in $k$ and $\sigma_{0}$, $i\sigma_{z}$
terms must be even in $k$. The even order terms are constants to first order
of $k$, however they break the $C_{3z}$ symmetry as they couple states with
angular momentum 0 to states with angular momenta $\pm 1$. The two odd order
terms $\sigma_{x,y}$ transform under $C_{3z}$ as
$C_{3z}\sigma_{x}C_{3z}^{-1}=\sigma_{x}\cos\frac{2\pi}{3}-\sigma_{y}\sin\frac{2\pi}{3}$,
$C_{3z}\sigma_{y}C_{3z}^{-1}=\sigma_{y}\cos\frac{2\pi}{3}+\sigma_{x}\sin\frac{2\pi}{3}$,
thus the odd terms can be chosen as
$v_{\star}^{\prime\prime}(k_{x}\sigma_{x}-k_{y}\sigma_{y})$, which is also
symmetric under $C_{2x}$.

###  S2.5 Summary of the single-particle Hamiltonian

According to the discussions in Sections S2.3 and S2.4, we can summarize the
total single-particle Hamiltonian as

$\displaystyle{\hat{H}}_{0}=$ $\displaystyle\sum_{\eta
s}\sum_{aa^{\prime}}\sum_{|\mathbf{k}|<\Lambda_{c}}(H^{(c,\eta)}_{a,a^{\prime}}(\mathbf{k})-\mu\delta_{aa^{\prime}})c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta s}-\mu\sum_{\eta
s}\sum_{\mathbf{R}}f_{\mathbf{R}\alpha\eta s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}$ $\displaystyle+\frac{1}{\sqrt{N}}\sum_{\eta
s}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\mathbf{R}}\sum_{\alpha
a}\left(e^{-|\mathbf{k}|^{2}\lambda^{2}/2-i\mathbf{k}\cdot\mathbf{R}}H^{(cf,\eta)}_{a\alpha}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}+h.c.\right)\ .$ (S81)

$\mu$ is the chemical potential, $\Lambda_{c}$ is the cutoff of the momentum
of the conduction bands, $H^{(c,\eta)}_{a,a^{\prime}}(\mathbf{k})$ and
$H^{(cf)}_{a\alpha}(\mathbf{k})$ are given by

$H^{(c,\eta)}(\mathbf{k})=\begin{pmatrix}0_{2\times 2}&v_{\star}(\eta
k_{x}\sigma_{0}+ik_{y}\sigma_{z})\\\ v_{\star}(\eta
k_{x}\sigma_{0}-ik_{y}\sigma_{z})&M\sigma_{x}\end{pmatrix},$ (S82)
$H^{(cf,\eta)}=\begin{pmatrix}\gamma\sigma_{0}+v_{\star}^{\prime}(\eta
k_{x}\sigma_{x}+k_{y}\sigma_{y})\\\ 0_{2\times 2}\end{pmatrix}\ ,$ (S83)

respectively. The parameters in the single-particle Hamiltonian are given in
Table S4.

In principle, the cutoff $\Lambda_{c}$ should be a small quantity compared to
the Moiré BZ. However, because the states at large $\mathbf{k}$’s have very
high energies, we can formally extend $\Lambda_{c}$ to infinity without
changing the low energy physics. A theoretical advantage of
$\Lambda_{c}\to\infty$ is that the yielded band structure will become periodic
over the moiré BZ. We label the momentum for the conduction bands as
$\mathbf{k}+\mathbf{G}$ with $\mathbf{k}$ being in the moiré BZ and
$\mathbf{G}$ moiré reciprocal lattice. Using a larger $\Lambda_{c}$ implies a
larger lattice of $\mathbf{G}$. In Fig. S2 we have plotted the band structures
of Section S2.5 using the cutoff $\Lambda_{c}=5\sqrt{3}k_{\theta}$.

As shown in Ref. Bernevig  _et al._ (2021a), the anti-unitary single-particle
symmetry $PC_{2z}T$ (anti-commuting with single-particle Hamiltonian) implies
a unitary many-body charge-conjugation symmetry ${\mathcal{P}}_{c}$ that
commutes with the Hamiltonian, i.e.,

${\mathcal{P}}_{c}f_{\mathbf{R}\alpha\eta
s}^{\dagger}{\mathcal{P}}_{c}^{-1}=\sum_{\alpha^{\prime}\eta^{\prime}}D^{f}_{\alpha^{\prime}\eta^{\prime},\alpha\eta}(PC_{2z}T)f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s},\qquad{\mathcal{P}}_{c}c_{\mathbf{k}a\eta
s}^{\dagger}{\mathcal{P}}_{c}^{-1}=\sum_{a^{\prime}\eta^{\prime}}D^{c}_{a^{\prime}\eta^{\prime},a\eta}(PC_{2z}T)c_{-\mathbf{k}a^{\prime}\eta^{\prime}s}\
,$ (S84)

where

$D^{f}(C_{2z}TP)=-\sigma_{y}\tau_{z},\qquad
D^{f}(C_{2z}TP)=(\sigma_{y}\tau_{z})\oplus(\sigma_{y}\tau_{z})\ .$ (S85)

One can verify that ${\mathcal{P}}_{c}$ commutes with ${\hat{H}}_{0}$ in
Section S2.5.

###  S2.6 The (first) chiral limit and the chiral U(4) symmetry

We find that when $v_{\star}^{\prime}=0$ the Hamiltonian (Section S2.5) anti-
commute with the unitary chiral operator

$Cf_{\mathbf{k}\alpha\eta
s}^{\dagger}C^{-1}=\sum_{\alpha^{\prime}\eta^{\prime}}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s}^{\dagger}D^{f}_{\alpha^{\prime}\eta^{\prime},a\eta}(C),\qquad
Cc_{\mathbf{k}a\eta
s}^{\dagger}C^{-1}=\sum_{a^{\prime}\eta^{\prime}}c_{\mathbf{k}a^{\prime}\eta^{\prime}s}^{\dagger}D^{c}_{a^{\prime}\eta^{\prime},a\eta}(C)$
(S86)

with

$D^{f}(C)=\sigma_{z}\tau_{0},\qquad
D^{c}(C)=(-\sigma_{z}\tau_{0})\oplus(\sigma_{z}\tau_{0})\ .$ (S87)

According to Eqs. S36, S2.2 and S65, we find that $C$ satisfies the algebra at
the single-particle level

$C^{2}=1,\qquad[C,T]=0,\qquad[C,C_{3z}]=0,\qquad\\{C,C_{2x}\\}=0,\qquad\\{C,C_{2z}T\\}=0,\qquad[C,P]=0\
,$ (S88)

which are same as the first chiral symmetry discussed in Section S1.2 (Eq.
S21) and in Refs. Song  _et al._ (2021); Bernevig  _et al._ (2021a);
Tarnopolsky  _et al._ (2019); Wang  _et al._ (2021). Therefore, we identify
$C$ as the first chiral symmetry. The band structure after imposing this
chiral symmetry is shown in Fig. S4(b).

Refs. Bernevig  _et al._ (2021a); Bultinck  _et al._ (2020a); Kang and Vafek
(2019) have shown that the presence of chiral symmetry implies a so-called
chiral non-flat U(4) symmetry of the projected interaction Hamiltonian. (A
similar U(4) symmetry was also found in Ref. Kang and Vafek (2019).) Here we
find that $C$ indeed implies a U(4) symmetry of the single-particle
Hamiltonian by using the symmetry $C_{2z}PC$. In next paragraph we will show
the chiral implied U(4) is consistent with the chiral non-flat U(4) discussed
in Refs. Bernevig  _et al._ (2021a); Bultinck  _et al._ (2020a); Kang and
Vafek (2019). In Section S3.6.2 we will show that the chiral U(4) symmetry is
also an approximate symmetry of the interaction Hamiltonian. Since the
Hamiltonian anti-commutes with $P$ and $C$, it must commute with the product
$PC$ as well as $C_{2z}PC$. As discussed in Section S1.2, $C$ is local in real
space, $P$ transforms $\mathbf{r}$ to $-\mathbf{r}$, hence $C_{2z}PC$ is also
local in real space. Its representation matrices are

$D^{f}(C_{2z}PC)=\sigma_{x}\tau_{y},\qquad
D^{c}(C_{2z}PC)=(\sigma_{x}\tau_{y})\oplus(-\sigma_{x}\tau_{y}).$ (S89)

We introduce the continuous symmetry generator

${\hat{\Theta}}_{y0}=\frac{1}{2}\sum_{\mathbf{k}\in\mathrm{MBZ}}\sum_{\alpha\alpha^{\prime}\eta\eta^{\prime}s}[\sigma_{x}\tau_{y}]_{\alpha\eta,\alpha^{\prime}\eta^{\prime}}f_{\mathbf{k}\alpha\eta
s}^{\dagger}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s}+\frac{1}{2}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{aa^{\prime}\eta\eta^{\prime}s}[(\sigma_{x}\tau_{y})\oplus(-\sigma_{x}\tau_{y})]_{a\eta,a^{\prime}\eta^{\prime}}c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s}\ ,$ (S90)

and it commutes with $H_{0}$. Then we find that the Hamiltonian (Section S2.5)
has the continuous symmetry

$[e^{-i\theta{\hat{\Theta}}_{y0}},\hat{H}_{0}]=0\ .$ (S91)

${\hat{\Theta}}_{y0}$ together with the U(2)$\times$U(2) symmetry form the
U(4) symmetry group. The generators of the U(4) group can be explicitly
written as

${\hat{\Theta}}_{\mu\nu}=\frac{1}{2}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\left(\sum_{\mathbf{k}\in\mathrm{MBZ}}\sum_{\alpha\alpha^{\prime}}\Theta^{(\mu\nu,f)}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}f_{\mathbf{k}\alpha\eta
s}^{\dagger}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s^{\prime}}+\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{aa^{\prime}}\Theta^{(\mu\nu,c)}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s^{\prime}}\right)\ ,$ (S92)

where $\mu,\nu=0,x,y,z$ and

$\Theta^{(0\nu,f)}=\sigma_{0}\tau_{0}{\varsigma}_{\nu},\qquad\Theta^{(0\nu,c)}=(\sigma_{0}\tau_{0}{\varsigma}_{\nu})\oplus(\sigma_{0}\tau_{0}{\varsigma}_{\nu})\
,$ (S93)
$\Theta^{(x\nu,f)}=\sigma_{x}\tau_{x}{\varsigma}_{\nu},\qquad\Theta^{(x\nu,c)}=(\sigma_{x}\tau_{x}{\varsigma}_{\nu})\oplus(-\sigma_{x}\tau_{x}{\varsigma}_{\nu})\
,$ (S94)
$\Theta^{(y\nu,f)}=\sigma_{x}\tau_{y}{\varsigma}_{\nu},\qquad\Theta^{(y\nu,c)}=(\sigma_{x}\tau_{y}{\varsigma}_{\nu})\oplus(-\sigma_{x}\tau_{y}{\varsigma}_{\nu})\
,$ (S95)
$\Theta^{(z\nu,f)}=\sigma_{0}\tau_{z}{\varsigma}_{\nu},\qquad\Theta^{(z\nu,c)}=(\sigma_{0}\tau_{z}{\varsigma}_{\nu})\oplus(\sigma_{0}\tau_{z}{\varsigma}_{\nu})\
.$ (S96)

Here $\sigma_{\mu}$, $\tau_{\mu}$, ${\varsigma}_{\mu}$ are Pauli matrices for
the orbital, valley, and spin degree of freedom, respectively. The first and
second eight-by-eight blocks in $\Theta^{(\mu\nu,c)}$ are for the $\Gamma_{3}$
orbitals ($a=1,2$) and the $\Gamma_{1}\oplus\Gamma_{2}$ ($a=3,4$),
respectively.

Figure S4: Band structures in the first, third chiral or flat limits. (a) Band
structure of the original topological heavy fermion model (Section S2.5) at
$w_{0}/w_{1}=0.8$. (b) Band structure in the first chiral symmetry, where
$v_{\star}^{\prime}$ is imposed to zero. (c) Band structure in the third
chiral limit (flat limit), where $M$ is imposed to zero. (d) Band structure in
the flat limit in the absence of the third chiral symmetry, where $M$ is
imposed to zero and a chiral breaking term $A$ is taken into account. (See
Section S2.8). The other parameters are given in Table S4.

In the end we explicitly show that the U(4) symmetry in this chiral limit is
consistent with the (first) chiral non-flat U(4) symmetry in Refs. Bernevig
_et al._ (2021a); Kang and Vafek (2019); Bultinck  _et al._ (2020a). In the
Appendix B2 of Ref. Bernevig  _et al._ (2021a), the sewing matrices of $T$,
$C_{2z}$, and $P$ have been fixed as $\mathbf{k}$-independent matrices

$B(T)=\sigma_{0}\tau_{x},\qquad B(C_{2z})=\sigma_{0}\tau_{x},\qquad
B(P)=-i\sigma_{y}\tau_{z}\ .$ (S97)

The generators of the (first) chiral non-flat U(4) symmetry in Ref. Bernevig
_et al._ (2021a) in this gauge are

$\sigma_{0}\tau_{0}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{x}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{y}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{z}{\varsigma}_{\nu},\qquad(\nu=0,x,y,z)\
,$ (S98)

where $\sigma_{\mu}$, $\tau_{\mu}$, abd ${\varsigma}_{\mu}$ are Pauli matrices
for the band, valley, and spin degrees of freedom, respectively. As shown in
Section S2.1 and Fig. S3, in most area of the moiré BZ, the middle two bands
are mainly contributed by the local orbitals. Thus we expect that, for
$\mathbf{k}\neq 0$ where the overlap between local orbitals and the middle two
bands is finite, the U(4) generators on the local orbitals (after certain
gauge transformation) should be same as those in Eq. S98. We apply the gauge
transformation
$U=e^{i\frac{\pi}{4}\sigma_{0}\tau_{z}}e^{-i\frac{\pi}{4}\sigma_{x}\tau_{z}}$
to the local orbitals such that the transformed representation matrices
$UD^{f}(g)U^{\dagger}$ become same as Eq. S97. (See Eq. S36 for definition of
$D^{f}(g)$’s.) The U(4) generators for the local orbitals in Eqs. S93, S94,
S95 and S96 transform to

$\sigma_{0}\tau_{0}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{y}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{x}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{z}{\varsigma}_{\nu},\qquad(\nu=0,x,y,z)\
.$ (S99)

which are indeed same as Eq. S98. For $\mathbf{k}=0$, our local orbital basis
(forming the representation $\Gamma_{3}$) is different from the basis in Ref.
Bernevig  _et al._ (2021a) (forming the representation
$\Gamma_{1}\oplus\Gamma_{2}$), thus the above comparison is only valid for
$\mathbf{k}\neq 0$.

###  S2.7 The flat limit and the flat U(4) symmetry

Refs. Bernevig  _et al._ (2021a); Bultinck  _et al._ (2020a) discussed a flat
limit where the kinetic energy of the flat middle two bands are discarded. A
so-called non-chiral flat U(4) symmetry emerges in this limit Bernevig  _et
al._ (2021a); Bultinck  _et al._ (2020a). We find that in the topological
heavy fermion model we can achieve an exact flat band if we set $M=0$, i.e.,
the $\Gamma_{1}-\Gamma_{2}$ splitting is set to zero. In the limit $M=0$, the
Hamiltonian (Section S2.5) anti-commutes with the unitary chiral operator

$Sf_{\mathbf{k}\alpha\eta
s}^{\dagger}S^{-1}=\sum_{\alpha^{\prime}\eta^{\prime}}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s}^{\dagger}D^{f}_{\alpha^{\prime}\eta^{\prime},a\eta}(S),\qquad
Sc_{\mathbf{k}a\eta
s}^{\dagger}S^{-1}=\sum_{a^{\prime}\eta^{\prime}}c_{\mathbf{k}a^{\prime}\eta^{\prime}s}^{\dagger}D^{c}_{a^{\prime}\eta^{\prime},a\eta}(S)$
(S100)

with

$D^{f}(S)=\sigma_{0}\tau_{0},\qquad
D^{c}(S)=(-\sigma_{0}\tau_{0})\oplus(\sigma_{0}\tau_{0})\ .$ (S101)

Since $S$ is different from the first and second chiral symmetries discussed
in Ref. Song  _et al._ (2021); Bernevig  _et al._ (2021a), we refer to $S$ as
the third chiral symmetry. The local orbitals and the $a=3,4$ conduction band
basis have the chiral eigenvalue 1, whereas the $a=1,2$ conduction band basis
has the chiral eigenvalue $-1$. Since only hoppings between opposite chiral
eigenvalues are allowed (such that the Hamiltonian anti-commutes with $S$),
the rank of the single-particle Hamiltonian at each $\mathbf{k}$ (in a valley-
spin sector) is at most 4. Therefore, in this limit there must be two exact
flat bands at the zero energy. Correspondingly, we also refer to this chiral
limit as the flat limit. The band structure in this limit is given in Fig.
S4(c).

Similar to the first chiral limit, the third chiral limit also implies a U(4)
symmetry of the single-particle Hamiltonian. In Section S3.6.2 we will show
that the flat U(4) symmetry is also an approximate symmetry of the interaction
Hamiltonian. Since the Hamiltonian anti-commutes with $P$ and $S$, it must
commute with the product $PS$ as well as $C_{2z}PS$. $S$ is local in real
space, while $P$ transforms $\mathbf{r}$ to $-\mathbf{r}$, hence $C_{2z}PS$ is
also local in real space. Its representation matrices are

$D^{f}(C_{2z}PS)=-i\sigma_{y}\tau_{y},\qquad
D^{c}(C_{2z}PS)=(-i\sigma_{y}\tau_{y})\oplus(i\sigma_{y}\tau_{y})\ .$ (S102)

We introduce the continuous symmetry generator

${\hat{\Sigma}}_{y0}=\frac{1}{2}\sum_{\mathbf{k}\in\mathrm{MBZ}}\sum_{\alpha\alpha^{\prime}\eta\eta^{\prime}s}[\sigma_{y}\tau_{y}]_{\alpha\eta,\alpha^{\prime}\eta^{\prime}}f_{\mathbf{k}\alpha\eta
s}^{\dagger}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s}+\frac{1}{2}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{aa^{\prime}\eta\eta^{\prime}s}[(\sigma_{y}\tau_{y})\oplus(-\sigma_{y}\tau_{y})]_{a\eta,a^{\prime}\eta^{\prime}}c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s}\ ,$ (S103)

and it commutes with ${\hat{H}}_{0}$. Then we find that the Hamiltonian
(Section S2.5) has the continuous symmetry

$[e^{-i\theta{\hat{\Sigma}}_{y0}},\hat{H}_{0}]=0\ .$ (S104)

Similar to the case in the first chiral symmetry, the above rotation together
with the U(2)$\times$U(2) symmetry lead to a U(4) symmetry group. The
generators of the U(4) group can be explicitly written as

${\hat{\Sigma}}_{\mu\nu}=\frac{1}{2}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\left(\sum_{\mathbf{k}\in\mathrm{MBZ}}\sum_{\alpha\alpha^{\prime}}\Sigma^{(\mu\nu,f)}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}f_{\mathbf{k}\alpha\eta
s}^{\dagger}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s^{\prime}}+\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{aa^{\prime}}\Sigma^{(\mu\nu,c)}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s^{\prime}}\right)\ ,$ (S105)

where $\mu,\nu=0,x,y,z$ and

$\Sigma^{(0\nu,f)}=\sigma_{0}\tau_{0}{\varsigma}_{\nu},\qquad\Sigma^{(0\nu,c)}=(\sigma_{0}\tau_{0}{\varsigma}_{\nu})\oplus(\sigma_{0}\tau_{0}{\varsigma}_{\nu})\
,$ (S106)
$\Sigma^{(x\nu,f)}=\sigma_{y}\tau_{x}{\varsigma}_{\nu},\qquad\Sigma^{(x\nu,c)}=(\sigma_{y}\tau_{x}{\varsigma}_{\nu})\oplus(-\sigma_{y}\tau_{x}{\varsigma}_{\nu})\
,$ (S107)
$\Sigma^{(y\nu,f)}=\sigma_{y}\tau_{y}{\varsigma}_{\nu},\qquad\Sigma^{(y\nu,c)}=(\sigma_{y}\tau_{y}{\varsigma}_{\nu})\oplus(-\sigma_{y}\tau_{y}{\varsigma}_{\nu})\
,$ (S108)
$\Sigma^{(z\nu,f)}=\sigma_{0}\tau_{z}{\varsigma}_{\nu},\qquad\Sigma^{(z\nu,c)}=(\sigma_{0}\tau_{z}{\varsigma}_{\nu})\oplus(\sigma_{0}\tau_{z}{\varsigma}_{\nu})\
.$ (S109)

Here $\sigma_{\mu}$, $\tau_{\mu}$, ${\varsigma}_{\mu}$ are Pauli matrices for
the orbital, valley, and spin degree of freedom, respectively. The first and
second eight-by-eight blocks in $\Sigma^{(\mu\nu,c)}$ are for the orbitals
$a=1,2$ and $a=3,4$, respectively.

Now we explicitly show that the U(4) symmetry in the third chiral limit (flat
limit) of our heavy fermion model is consistent with the non-chiral flat U(4)
symmetry in Ref. Bernevig  _et al._ (2021a); Bultinck  _et al._ (2020a), which
commutes with the interaction Hamiltonian but anti-commutes with the single-
particle BM Hamiltonian. It is worth mentioning that our U(4) rotations,
unlike the ones in Ref. Bernevig  _et al._ (2021a), commutes with both the
interaction Hamiltonian (approximately) and the single-particle Hamiltonian
(if $M=0$). In the Appendix B2 of Ref. Bernevig  _et al._ (2021a), the sewing
matrices of $T$, $C_{2z}$, and $P$ have been fixed as $\mathbf{k}$-independent
matrices shown in Eq. S97. The generators of the non-chiral flat U(4) symmetry
in this gauge are

$\sigma_{0}\tau_{0}{\varsigma}_{\nu},\qquad\sigma_{y}\tau_{x}{\varsigma}_{\nu},\qquad\sigma_{y}\tau_{y}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{z}{\varsigma}_{\nu},\qquad(\nu=0,x,y,z)\
,$ (S110)

where $\sigma_{\mu}$, $\tau_{\mu}$, abd ${\varsigma}_{\mu}$ are Pauli matrices
for the band, valley, and spin degrees of freedom, respectively. As shown in
Section S2.1 and Fig. S3, in most area of the moiré BZ, the middle two bands
are mainly contributed by the local orbitals. Thus we expect that, for
$\mathbf{k}\neq 0$ where the overlap between local orbitals and the middle two
bands is finite, the U(4) generators on the local orbitals (after certain
gauge transformation) should be same as those in Eq. S110. We apply the gauge
transformation
$U=e^{i\frac{\pi}{4}\sigma_{0}\tau_{z}}e^{-i\frac{\pi}{4}\sigma_{x}\tau_{z}}$
to the local orbitals such that the transformed representation matrices
$UD^{f}(g)U^{\dagger}$ become same as Eq. S97. (See Eq. S36 for definition of
$D^{f}(g)$’s.) The U(4) generators for the local orbitals in Eqs. S106, S107,
S108 and S109 transform to

$\sigma_{0}\tau_{0}{\varsigma}_{\nu},\qquad-\sigma_{y}\tau_{y}{\varsigma}_{\nu},\qquad\sigma_{y}\tau_{x}{\varsigma}_{\nu},\qquad\sigma_{0}\tau_{z}{\varsigma}_{\nu},\qquad(\nu=0,x,y,z)\
.$ (S111)

which generate the same U(4) rotations as Eq. S110. For $\mathbf{k}=0$, our
local orbital basis (forming the representation $\Gamma_{3}$) is different
from the basis in Ref. Bernevig  _et al._ (2021a) (forming the representation
$\Gamma_{1}\oplus\Gamma_{2}$), thus the above comparison is only valid for
$\mathbf{k}\neq 0$.

As will be revealed by the numerical calculations in Section S4.2, our
flat-U(4) symmetry is a better approximation than our (first) chiral-U(4)
symmetry. A heuristic argument is the following. The flat-U(4) breaking
parameter $M$ only affects low energy single-particle states of $H^{(c,\eta)}$
(Eq. S82) because the $M$ term is just a weak perturbation for high energy
states, while the chiral-U(4) breaking parameter $v_{\star}^{\prime}$ mainly
affects high energy single-particle states because it enters the Hamiltonian
in the form of $v_{\star}^{\prime}k_{x}\sigma_{0}\pm
iv_{\star}^{\prime}k_{y}\sigma_{z}$. Roughly speaking, in the presence of
$v_{\star}^{\prime}$ and $M$, low energy components of the single-particle
Hamiltonian do not commute with the flat-U(4) rotation, while high energy
components do not commute with the chiral-U(4) rotation. Therefore, the energy
change after a chiral-U(4) rotation will be larger than that after a flat-U(4)
rotation, implying that the flat-U(4) symmetry is better approximation of the
Hamiltonian.

In the end of this section we split the flat-U(4) generators to the
fundamental representations of the flat-U(4) group. In Ref. Bernevig  _et al._
(2021a), the fundamental basis is obtained by diagonalizing the Pauli matrix
$\sigma_{y}$ in the band basis that is fixed by the sewing matrices in Eq.
S97. As discussed below Eq. S110, for $\mathbf{k}\neq 0$, the $f$-orbital
basis is related to the basis of Ref. Bernevig  _et al._ (2021a) by the
rotation
$U=e^{i\frac{\pi}{4}\sigma_{0}\tau_{z}}e^{-i\frac{\pi}{4}\sigma_{x}\tau_{z}}$.
Thus now the fundamental U(4) representations should be obtained by
diagonalizing $U^{\dagger}\sigma_{y}U=-\sigma_{z}\tau_{z}$. In the following
paragraph we prove that the eigenvalue $\xi$ of $\sigma_{z}\tau_{z}$ indeed
labels fundamental U(4) representations. It is worth mentioning that the
choice of fundamental U(4) representations is not unique. For example, since
all the U(4) rotations in Eqs. S106, S107, S108 and S109 commute with
$\sigma_{y}$, we can also use eigenvalues of $\sigma_{y}$ in our heavy fermion
basis to label fundamental representations. We choose $\xi$ because (i) it has
the physical meaning of Chern number (Section S4.5) Bernevig  _et al._ (2021a)
and (ii) it simplifies the exchange interaction (Section S3.3).

First, the flat-U(4) operators are already block-diagonal: They do not mix the
eight $f$-orbitals, eight $\Gamma_{3}$ basis of $c$-bands, and eight
$\Gamma_{1}\oplus\Gamma_{2}$ basis of $c$-bands with each other. We need to
further split each block into two fundamental representations of the flat-U(4)
group. We notice that the flat-U(4) generators either do not change the
$\alpha,\eta$ ($a,\eta$) indices of the $f$-electrons ($c$-electrons) or flip
them at the same time. Therefore, the index $\xi=(-1)^{\alpha-1}\eta$
($\xi=(-1)^{a-1}\eta$) of $f$-electrons ($c$-electrons) is unchanged under the
flat-U(4) rotations. (In general, $\xi$ is not conserved by ${\hat{H}}_{0}$.
For example, the $M$ and $v_{\star}^{\prime}$ terms in Eqs. S82 and S83 flip
the $\xi$ index.) Consequently, the four components with the same $\xi$ in
each block form a fundamental representation of the flat-U(4) group. For later
convenience, we introduce the U(4) generator operators for fundamental
representations separately

$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}=$
$\displaystyle\frac{1}{2}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\sum_{\mathbf{k}\in\mathrm{MBZ}}\sum_{\alpha\alpha^{\prime}}\delta_{\xi,(-1)^{\alpha-1}\eta}\delta_{\xi,(-1)^{\alpha^{\prime}-1}\eta^{\prime}}\Sigma^{(\mu\nu,f)}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}f_{\mathbf{k}\alpha\eta
s}^{\dagger}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s^{\prime}}\ ,$ (S112)
$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}=$
$\displaystyle\frac{1}{2}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{aa^{\prime}=1,2}\delta_{\xi,(-1)^{a-1}\eta}\delta_{\xi,(-1)^{a^{\prime}-1}\eta^{\prime}}\Sigma^{(\mu\nu,c)}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s^{\prime}}\ ,$ (S113)
$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}=$
$\displaystyle\frac{1}{2}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{aa^{\prime}=3,4}\delta_{\xi,(-1)^{a-1}\eta}\delta_{\xi,(-1)^{a^{\prime}-1}\eta^{\prime}}\Sigma^{(\mu\nu,c)}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s^{\prime}}\ ,$ (S114)

respectively. The total flat-U(4) generators are given by
${\hat{\Sigma}}_{\mu\nu}=\sum_{\xi}\left({\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}+{\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}+{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}\right)$.

###  S2.8 Flat bands beyond the third chiral symmetry

Consider a system consisting of sublattices $A$ and $B$, if only hoppings
between $A$\- and $B$-sublattices are allowed, then the system has a chiral
symmetry, eigenvalues of which are $1$ and $-1$ for the two sublattices,
respectively. Then the Hamiltonian matrix $H(\mathbf{k})$ satisfies

$H_{\alpha,\alpha^{\prime}}(\mathbf{k})=H_{\beta,\beta^{\prime}}(\mathbf{k})=0\
,$ (S115)

where $\alpha\in A$, $\beta\in B$. We use $L_{A}$ and $L_{B}$ to represent the
number of orbitals (per unit cell) in the $A$\- and $B$-sublattices,
respectively. Without loss of generality, in the following we always label the
sublattice with more (or equal number) orbitals as $A$ such that $L_{A}\geq
L_{B}$. Because the dimension of the Hamiltonian matrix $H(\mathbf{k})$ is
$L_{A}+L_{B}$ and the rank of the Hamiltonian is at most $2L_{B}$, which are
contributed by $L_{B}$ nonzero columns $H_{\alpha,\beta}(\mathbf{k})$ and
$L_{B}$ nonzero rows $H_{\beta,\alpha}(\mathbf{k})$ ($\alpha\in A$, $\beta\in
B$), $H(\mathbf{k})$ will have at least $L_{A}-L_{B}$ zero-energy flat bands
if $L_{A}>L_{B}$. A recent study on flat bands based on the so-called
$S$-matrix Călugăru  _et al._ (2021) pointed out that, in many cases, breaking
the chiral symmetry does not necessarily destroy the flatness of flat bands. A
general condition for the presence of flat bands and their topological
classifications are derived in Ref. Călugăru  _et al._ (2021). Here we only
focus on a simple situation, where intra-sublattice hoppings in the smaller
sublattice ($B$) are taken into account. Adding such hoppings do break the
chiral symmetry but do not change the rank of $H(\mathbf{k})$ — with
$H_{\beta,\beta^{\prime}}(\mathbf{k})$ being nonzero, $H(\mathbf{k})$ still
has at most $L_{B}$ nonzero columns and $L_{B}$ nonzero rows. Therefore, the
flatness is robust against the hoppings within the smaller sublattice.

The above discussion explains the flat bands when we take into account certain
chiral breaking terms in the heavy fermion Hamiltonian. Following the symmetry
argument in Section S2.3, we can add the quadratic term
$A(k_{x}^{2}-k_{y}^{2})\sigma_{x}-2\eta Ak_{x}k_{y}\sigma_{y}$ to the $a=1,2$
block of Eq. S82. This term commutes with $S$ (Eq. S100) and hence breaks the
third chiral symmetry. However, it does not affect the flatness of the flat
bands as long as $M=0$, as shown in Fig. S4(d), where
$A=-30\mathrm{meV}/k_{\theta}^{2}$, $M=0$. Now we explain the flatness in the
presence of $A$ using the $S$-matrix theory Călugăru  _et al._ (2021). In
Section S2.7 we have shown that the eigenvalues of $S$ are $-1$, $+1$, and
$+1$ for the $\Gamma_{3}$ conduction band basis ($a=1,2$),
$\Gamma_{1}\oplus\Gamma_{2}$ conduction band basis ($a=3,4$), and the local
orbitals, respectively. Before we add the $A$ term, only hoppings between
bases with opposite $S$ eigenvalues are present. We thus can label the
$\Gamma_{3}$ conduction band basis ($a=1,2$) and the local orbitals as the
$A$-sublattice, and the $\Gamma_{1}\oplus\Gamma_{2}$ conduction band basis
($a=3,4$) as the $B$-sublattice. There is $L_{A}-L_{B}=2$ and hence there will
be two flat bands. Since the $A$ term is nothing but the hopping within the
smaller ($B$) sublattice, according to the discussion in the last paragraph,
the flatness of flat bands are stable against $A$.

###  S2.9 The U(4)$\times$U(4) symmetry in chiral flat limit

As discussed in Refs. Bernevig  _et al._ (2021a); Bultinck  _et al._ (2020a),
the two-band projected interaction Hamiltonian will have a U(4)$\times$U(4)
group when both the chiral and flat limits are achieved. It should be
emphasized that neither chiral-U(4) nor flat-U(4) is a factor U(4) group of
the total U(4)$\times$U(4) group. In the heavy fermion basis, when both limits
are achieved ($M=v_{\star}^{\prime}=0$), the $\xi$ index defined in Section
S2.7, which equals $(-1)^{\alpha-1}\eta$ and $(-1)^{a-1}\eta$ for $f$\- and
$c$-electrons, respectively, and labels the fundamental flat-U(4)
representations, is conserved in ${\hat{H}}_{0}$. Moreover, every bilinear
term in ${\hat{H}}_{0}$ only involves a single $\xi$ flavor, implying that we
can write ${\hat{H}}_{0}$ as ${\hat{H}}_{0}=\sum_{\xi}{\hat{H}}_{0}^{(\xi)}$
where ${\hat{H}}_{0}^{(\xi)}$ only contains fermion operators of the $\xi$
flavor. Then, not only the total flat-U(4) generators
${\hat{\Sigma}}_{\mu\nu}$, but also the generators acting in each $\xi$
flavor, i.e.,
${\hat{\Sigma}}_{\mu\nu}^{(\xi)}={\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}+{\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}+{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}$
(Eqs. S112, S113 and S114), commute with ${\hat{H}}_{0}$ because there is
$[{\hat{\Sigma}}_{\mu\nu}^{(\xi)},{\hat{H}}_{0}^{(\xi)}]=0$ for the two
$\xi$’s separately. A generic rotation that is commuting with ${\hat{H}}_{0}$
can be parameterized as

$\exp\left(-i\sum_{\mu\nu\xi}\theta_{\mu\nu}^{(\xi)}{\hat{\Sigma}}_{\mu\nu}^{(\xi)}\right)\
,$ (S116)

where $\theta_{\mu\nu}^{(+)}$ can be different with $\theta_{\mu\nu}^{(-)}$.
Therefore, when $M=v_{\star}^{\prime}=0$, ${\hat{H}}_{0}$ has a
U(4)$\times$U(4) symmetry group with the two factor U(4)’s generated by
${\hat{\Sigma}}_{\mu\nu}^{(+)}$ and ${\hat{\Sigma}}_{\mu\nu}^{(-)}$,
respectively.

| ${\hat{H}}_{0}$ | ${\hat{H}}_{0}(M=0)$ | ${\hat{H}}_{0}(v_{\star}^{\prime}=0)$ | ${\hat{H}}_{0}(M=v_{\star}^{\prime}=0)$ | ${\hat{H}}_{U}+{\hat{H}}_{V}+{\hat{H}}_{W}+{\hat{H}}_{J}+{\hat{H}}_{\widetilde{J}}$ | ${\hat{H}}_{K}$
---|---|---|---|---|---|---
Local Symmetry | U(2)$\times$U(2) | flat-U(4) | chiral-U(4) | U(4)$\times$U(4) | U(4)$\times$U(4) | chiral-U(4)
Table S5: Summary of the local symmetries of the single-particle and
interaction Hamiltonians. The single-particle Hamiltonian ${\hat{H}}_{0}$ is
summarized in Section S2.5, and interaction Hamiltonians ${\hat{H}}_{X}$
($X=U,V,W,J,\widetilde{J},K$) are summarized in Section S3.6.1.

Now we prove that, in general, when the single-particle or interaction
Hamiltonian has both the chiral-U(4) and flat-U(4) symmetry, it must at least
have the SU(4)$\times$SU(4) symmetry. Since
$\delta_{\xi,(-1)^{a-1}\eta}\delta_{\xi,(-1)^{a^{\prime}-1}\eta^{\prime}}$,
$\sigma_{y}\tau_{x}$, $\sigma_{x}\tau_{y}$ enter the coefficient matrices of
${\hat{\Sigma}}^{(\xi)}_{x\nu}$, ${\hat{\Sigma}}_{x\nu}$, and
${\hat{\Theta}}_{y\nu}$, respectively, and

$\delta_{\xi,(-1)^{a-1}\eta}\delta_{\xi,(-1)^{a^{\prime}-1}\eta^{\prime}}[\sigma_{y}\tau_{x}]_{\alpha\eta,\alpha^{\prime}\eta^{\prime}}=\frac{1}{2}[\sigma_{y}\tau_{x}+\xi\sigma_{x}\tau_{y}]$
(S117)

we have

${\hat{\Sigma}}^{(\xi)}_{x\nu}=\frac{1}{2}{\hat{\Sigma}}_{x\nu}+\frac{\xi}{2}{\hat{\Theta}}_{y\nu}\
.$ (S118)

Similarly, we have

${\hat{\Sigma}}^{(\xi)}_{y\nu}=\frac{1}{2}{\hat{\Sigma}}_{y\nu}+\frac{\xi}{2}{\hat{\Theta}}_{x\nu}\
.$ (S119)

Therefore, when the Hamiltonian commutes with ${\hat{\Sigma}}_{\mu\nu}$ and
${\hat{\Theta}}_{\mu\nu}$ separately, it must also commute with
${\hat{\Sigma}}^{(\xi)}_{x\nu}$ and ${\hat{\Sigma}}^{(\xi)}_{y\nu}$ for
$\xi=\pm$ separately, which generate a SU(4)$\times$SU(4) group. To see if the
Hamiltonian has a higher U(4)$\times$U(4) symmetry, one has to explicitly
check if the Hamiltonian also commute ${\hat{\Sigma}}^{(\xi)}_{00}$ for
$\xi=\pm$ separately. Since the Hamiltonian must respect charge-U(1), i.e.,
${\hat{\Sigma}}^{(+)}_{00}+{\hat{\Sigma}}^{(-)}_{00}$, one only need to check
if ${\hat{\Sigma}}^{(+)}_{00}-{\hat{\Sigma}}^{(-)}_{00}$, or equivalently the
index $\xi$, is conserved in the Hamiltonian. In summary, the Hamiltonian will
at least have a SU(4)$\times$SU(4) symmetry if it respects both flat-U(4) and
chiral-U(4) symmetries. Given the charge-U(1) is respected, the
SU(4)$\times$SU(4) group will be promoted to U(4)$\times$U(4) if the
Hamiltonian further commutes with
${\hat{\Sigma}}^{(+)}_{00}-{\hat{\Sigma}}^{(-)}_{00}$.

Using the criterion derived in the above paragraph, we find that most terms in
the interaction Hamiltonian also has the U(4)$\times$U(4) symmetry (Section
S3.6.1). We summarize the local symmetries of all the single-particle and
interaction Hamiltonians in Table S5.

###  S2.10 Fubini-Study metric and Berry’s curvature of the flat bands in the
chiral flat limit

It has been shown that the flat (Chern) bands of the BM model in the chiral
limit are analytically solvable Tarnopolsky  _et al._ (2019). Because the
Chern band solutions are holomorphic functions of $k_{x}+ik_{y}$ (or
$k_{x}-ik_{y}$), there is a simple relation between the Berry’s curvature
$\Omega(\mathbf{k})$ and the Fubini-Study metric $g_{ij}(\mathbf{k})$, i.e., ,
$g_{ij}(\mathbf{k})=\frac{1}{2}|\Omega(\mathbf{k})|\delta_{ij}$ Ledwith  _et
al._ (2020), which further leads to $4\det
g(\mathbf{k})=|\Omega(\mathbf{k})|^{2}$, also known as the ideal droplet
condition. The ideal droplet condition together with the flatness of Berry’s
curvature, i.e., $\Omega(\mathbf{k})=const.$, lead to the GMP algebra of
density operators Girvin  _et al._ (1986); Parameswaran  _et al._ (2013); Roy
(2014) and hence can give rise to fractional Chern insulator phases in spin
and valley polarized MATBG Sheffer and Stern (2021); Repellin and Senthil
(2020); Abouelkomsan  _et al._ (2020); Ledwith  _et al._ (2020).

In this subsection we will show that in the chiral flat limit
($v_{\star}^{\prime}=M=0$) of our heavy fermion model the Chern band states
are also analytically solvable and are holomorphic functions of $k_{x}+ik_{y}$
or $k_{x}-ik_{y}$ (up to a normalization constant). Thus the ideal droplet
condition of the flat band will be automatically satisfied. (It is worth
noting that this chiral limit is theoretically achieved by imposing
$v_{\star}^{\prime}=M=0$. It shares the same symmetries as the BM model in the
chiral limit, where $w_{0}/w_{1}=0$. But the other parameters, e.g., $M$,
$v_{\star}$, are not required to be obtained at $w_{0}/w_{1}=0$.) According to
Sections S2.7 and S2.6, when both (first) chiral and flat limits are achieved
the single-particle Hamiltonian must commute with $C\cdot S$, representation
of which is given by

$D^{f}(C\cdot S)=\sigma_{z}\tau_{0},\qquad D^{c}(C\cdot
S)=\sigma_{z}\tau_{0}\oplus\sigma_{z}\tau_{0}\ .$ (S120)

Therefore, we can block-diagonalize the Hamiltonian using the eigenvalues
($\zeta$) of $C\cdot S$. The total Hamiltonian in momentum space can be
rewritten as

$\displaystyle{\hat{H}}_{0}=$ $\displaystyle\sum_{\eta
s}\sum_{aa^{\prime}=1,3}\sum_{|\mathbf{k}|<\Lambda_{c}}H^{(c,\eta,+)}_{a,a^{\prime}}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta s}+\sum_{\eta
s}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{a=1,3}\left(e^{-|\mathbf{k}|^{2}\lambda^{2}/2}H^{(cf,\eta,+)}_{a1}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}f_{\mathbf{k}1\eta s}+h.c.\right)$ $\displaystyle+$
$\displaystyle\sum_{\eta
s}\sum_{aa^{\prime}=2,4}\sum_{|\mathbf{k}|<\Lambda_{c}}H^{(c,\eta,-)}_{a,a^{\prime}}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta s}+\sum_{\eta
s}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{a=2,4}\left(e^{-|\mathbf{k}|^{2}\lambda^{2}/2}H^{(cf,\eta,-)}_{a2}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}f_{\mathbf{k}2\eta s}+h.c.\right)$ (S121)

where the two-by-two matrix $H^{(c,\eta,\zeta)}$ and two $H^{(cf,\eta,\zeta)}$
(two-by-one) are given by

$H^{(c,\eta,\zeta)}=v_{\star}\begin{pmatrix}0&\eta k_{x}+i\zeta k_{y}\\\ \eta
k_{x}-i\zeta k_{y}&0\end{pmatrix},\qquad
H^{(cf,\eta,\zeta)}=\begin{pmatrix}\gamma\\\ 0\end{pmatrix}\ .$ (S122)

The Hamiltonian in each $\zeta$ sector still anti-commutes with $S$ and hence
must have flat bands (Section S2.7). As explained in Section S2.5, for the
band structure to be periodic in MBZ, we can formally extend the cutoff of
$c$-bands, $\Lambda_{c}$, to infinity. Here we just assume $\Lambda_{c}$ is
much larger than the size of MBZ. One must be aware that, for a moire
reciprocal lattice $\mathbf{G}$, there is $c_{\mathbf{k}+\mathbf{G},a\eta
s}\neq c_{\mathbf{k},a\eta s}$ and $f_{\mathbf{k}+\mathbf{G},\alpha\eta
s}=f_{\mathbf{k},\alpha\eta s}$ because the former is a continuous field while
the second is a lattice. Then the flat band solution is in the sector
$(\zeta,\eta,s)$ can be written as

$d_{\mathbf{k},\zeta=+,\eta,s}^{\dagger}=\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}f_{\mathbf{k},1,\eta,s}^{\dagger}+\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\sum_{\begin{subarray}{c}\mathbf{G}\\\
|\mathbf{k}+\mathbf{G}|<\Lambda_{c}\end{subarray}}\frac{\gamma/v_{\star}}{\eta(k_{x}+G_{x})+i(k_{y}+G_{y})}c_{\mathbf{k}+\mathbf{G},3,\eta,s}^{\dagger}\
,$ (S123)
$d_{\mathbf{k},\zeta=-,\eta,s}^{\dagger}=\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}f_{\mathbf{k},2,\eta,s}^{\dagger}+\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\sum_{\begin{subarray}{c}\mathbf{G}\\\
|\mathbf{k}+\mathbf{G}|<\Lambda_{c}\end{subarray}}\frac{\gamma/v_{\star}}{\eta(k_{x}+G_{x})-i(k_{y}+G_{y})}c_{\mathbf{k}+\mathbf{G},4,\eta,s}^{\dagger}\
,$ (S124)

where

$\mathcal{N}_{\mathbf{k}}=1+\sum_{\begin{subarray}{c}\mathbf{G}\\\
|\mathbf{k}+\mathbf{G}|<\Lambda_{c}\end{subarray}}\frac{\gamma^{2}/v_{\star}^{2}}{|\mathbf{k}+\mathbf{G}|^{2}}\
.$ (S125)

The states $d_{\mathbf{k},\zeta,\eta,s}^{\dagger}$ are also eigenstate of the
(first) chiral operator $C$ with eigenvalues equal to $\zeta$. Their Chern
numbers are computed to be $\zeta\cdot\eta$.

In order to calculate the quantum metric tensor, we rewrite the flat band
wavefunction in the first-quantized formalism. For simplicity, here we only
focus on the $\zeta=+$, $\eta=+$, ${\varsigma}=\uparrow$ sector and will omit
the indices $\zeta,\eta,{\varsigma}$. There is

$|u_{\mathbf{k}}\rangle=\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}|\widetilde{u}_{\mathbf{k}}\rangle,\qquad|\widetilde{u}_{\mathbf{k}}\rangle=|f_{1},\mathbf{k}\rangle+\sum_{\begin{subarray}{c}\mathbf{G}\\\
|\mathbf{k}+\mathbf{G}|<\Lambda_{c}\end{subarray}}\frac{\gamma/v_{\star}}{(k_{x}+G_{x})+i(k_{y}+G_{y})}|c_{3},\mathbf{k}+\mathbf{G}\rangle\
,$ (S126)

where $|f_{1},\mathbf{k}\rangle$ and $|c_{3},\mathbf{k}+\mathbf{G}\rangle$ are
$f$\- ($\alpha=1$) and $c$-band basis ($a=3$), respectively, and
$|\widetilde{u}_{\mathbf{k}}\rangle$ is a holomorphic function of
$k_{x}+ik_{y}$. The quantum metric can be written as a sum of four terms

$\displaystyle\mathfrak{g}_{ij}(\mathbf{k})=$
$\displaystyle\langle\partial_{k_{i}}u_{\mathbf{k}}|(1-|u_{\mathbf{k}}\rangle\langle
u_{\mathbf{k}}|)|\partial_{k_{j}}u_{\mathbf{k}}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{{\mathcal{N}_{\mathbf{k}}}}\langle\partial_{k_{i}}\widetilde{u}_{\mathbf{k}}|(1-|u_{\mathbf{k}}\rangle\langle
u_{\mathbf{k}}|)|\partial_{k_{j}}\widetilde{u}_{\mathbf{k}}\rangle+\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\left(\partial_{k_{i}}\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\right)\langle\widetilde{u}_{\mathbf{k}}|(1-|u_{\mathbf{k}}\rangle\langle
u_{\mathbf{k}}|)|\partial_{k_{j}}\widetilde{u}_{\mathbf{k}}\rangle$
$\displaystyle+$
$\displaystyle\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\left(\partial_{k_{j}}\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\right)\langle\partial_{k_{i}}\widetilde{u}_{\mathbf{k}}|(1-|u_{\mathbf{k}}\rangle\langle
u_{\mathbf{k}}|)|\widetilde{u}_{\mathbf{k}}\rangle+\left(\partial_{k_{i}}\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\right)\left(\partial_{k_{j}}\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\right)\langle\widetilde{u}_{\mathbf{k}}|(1-|u_{\mathbf{k}}\rangle\langle
u_{\mathbf{k}}|)|\widetilde{u}_{\mathbf{k}}\rangle\ .$ (S127)

Because $(1-|u_{\mathbf{k}}\rangle\langle
u_{\mathbf{k}}|)|\widetilde{u}_{\mathbf{k}}\rangle=0$, the last three terms in
the above equation vanish. Hence we have

$\mathfrak{g}_{ij}(\mathbf{k})=\frac{1}{{\mathcal{N}_{\mathbf{k}}}}\langle\partial_{k_{i}}\widetilde{u}(\mathbf{k})|(1-|u(\mathbf{k})\rangle\langle
u(\mathbf{k})|)|\partial_{k_{j}}\widetilde{u}(\mathbf{k})\rangle\ .$ (S128)

Since $\widetilde{u}_{\mathbf{k}}$ is holomorphic in $k=k_{x}+ik_{y}$, there
are
$\partial_{k_{x}}\widetilde{u}_{\mathbf{k}}=\partial_{k}\widetilde{u}_{\mathbf{k}}$,
$\partial_{k_{y}}\widetilde{u}_{\mathbf{k}}=i\partial_{k}\widetilde{u}_{\mathbf{k}}$,
and the quantum metric must have the form

$\mathfrak{g}_{ij}(\mathbf{k})=\begin{pmatrix}1&i\\\
-i&1\end{pmatrix}f(\mathbf{k})$ (S129)

where
$f(\mathbf{k})=\langle\partial_{k}u(\mathbf{k})|(1-|u(\mathbf{k})\rangle\langle
u(\mathbf{k})|)|\partial_{k}u(\mathbf{k})\rangle$ is a real function of
$\mathbf{k}$. The Fubini-Study metric and Berry’s curvature are given by
$g_{ij}(\mathbf{k})=\mathrm{Re}[\mathfrak{g}_{ij}(\mathbf{k})]=f(\mathbf{k})\delta_{ij}$
and
$\Omega(\mathbf{k})=2\mathrm{Im}[\mathfrak{g}_{xy}(\mathbf{k})]=2f(\mathbf{k})$,
respectively, thus the ideal droplet condition
$g_{ij}(\mathbf{k})=\frac{1}{2}|\Omega(\mathbf{k})|\delta_{ij}$ is satisfied.

We then investigate the flatness of their Berry’s curvatures. We first
consider the parameters $v_{\star}=1.623\mathrm{eV\cdot\mathring{A}}$,
$\gamma=-24.75$meV obtained at $w_{0}/w_{1}=0.8$ (Table S4). The corresponding
Berry’s curvature is shown in Fig. S5(b). We can see that it is not flat at
all. Its deviation from perfect flatness can be estimated by the number

$\delta
C=\sqrt{\int\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\left(\frac{\Omega(\mathbf{k})}{2\pi}-C\right)^{2}}\
,$ (S130)

which is computed to be 2.723 for $\gamma=-24.75$meV. The actual chiral limit
of the BM model, i.e., $w_{0}/w_{1}=0$, leads to a larger gap ($100.8$meV)
between the flat bands the lowest passive bands. That means in the actual
chiral limit $\gamma$ should be $\gamma=-100.8$meV. For comparison, in Fig.
S5(c) and (d) we present the distributions of Berry’s curvatures given by
$\gamma=-50$meV and $\gamma=-100.8$meV, respectively. It is clear that a
larger $|\gamma|$ gives a flatter Berry’s curvature. The derivations from
perfect flatness are 1.1585 and 0.2570 for $\gamma=-50$meV and
$\gamma=-100.8$meV, respectively.

In summary, because of the holomorphic or anti-holomorphic property of the
flat band wavefunction the ideal droplet condition is automatically achieved
in the chiral flat limit of our heavy fermion model, where
$v_{\star}^{\prime}$ and $M$ are artificially imposed to zero. The relative
flatness of the Berry’s curvature can also be reproduced by using parameters
fitting the actual BM model bands in the chiral limit. These lead to the GMP
algebra of density operators. Given that the Hamiltonian (Section S2.10) only
involves $k_{x}+ik_{y}$ or $k_{x}-ik_{y}$ and is periodic in real space, the
Jastrow type wavefunctions for fractional Chern insulators of spin and valley
polarized MATBG Ledwith  _et al._ (2020) should also be applicable to our
model.

Figure S5: Berry’s curvatures of the flat Chern bands in the (first) chiral
flat limit. (a) is an illustration of the moiré Brillouin zone and high
symmetry momenta. (b)-(d) are the Berry’s curvatures in the (first) chiral
flat limit, where $\gamma$ equals to -24.75meV, -50meV, -100.8meV,
respectively, and $v_{\star}^{\prime}=M=0$,
$v_{\star}=1.623\mathrm{eV\cdot\mathring{A}}$. The values of the Berry’s
curvature is normalized in the convention where the area of moiré Brillouin
zone is $4\pi^{2}$.

##  S3 The interaction Hamiltonian and its symmetries

###  S3.1 The projected density operator and Coulomb interaction

We can project the original basis onto the basis of the topological heavy
fermion model as

$c_{\mathbf{k},\mathbf{Q},\beta,\eta,s}^{\dagger}\approx\sum_{\alpha}\widetilde{v}^{(\eta)*}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})f_{\mathbf{k}\alpha\eta
s}^{\dagger}+\sum_{a}\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}$ (S131)

Correspondingly, the projected real space basis (Eq. S11) can be written as

$c_{l\beta\eta
s}^{\dagger}(\mathbf{r})\approx\sum_{\mathbf{R}\alpha}e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})f_{\mathbf{R}\alpha\eta
s}^{\dagger}+\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{a}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})c_{\mathbf{k},a\eta
s}^{\dagger}\ ,$ (S132)

where $\Lambda_{c}$ is the cutoff for the conduction band basis. We emphasize
that even though we formally introduce
$\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})$, in practical
calculations we will only use
$\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(0)$, as we did in the k$\cdot$p
expansion in Section S2.3. The phase factor
$e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}$ comes from the phase shift
gained by the real space basis under translation, as explained around Eq. S43.
We now can write the density operator as

$\displaystyle\hat{\rho}(\mathbf{r})=\sum_{\beta l\eta s}c_{l\beta\eta
s}^{\dagger}(\mathbf{r})c_{l\beta\eta s}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\sum_{\beta l\eta
s}\bigg{[}\sum_{\begin{subarray}{c}\mathbf{R}\alpha\\\
\mathbf{R}^{\prime}\alpha^{\prime}\end{subarray}}e^{i\eta\Delta\mathbf{K}_{l}\cdot(\mathbf{R}-\mathbf{R}^{\prime})}w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})w^{(\eta)}_{l\beta,\alpha^{\prime}}(\mathbf{r}-\mathbf{R}^{\prime})f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta s}+\frac{1}{\Omega_{\rm
tot}}\sum_{\begin{subarray}{c}|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}\\\
aa^{\prime}\\\
\mathbf{Q}\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta}\end{subarray}}e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}}\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})\widetilde{u}^{(\eta)}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k}^{\prime})c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta s}$ $\displaystyle+$
$\displaystyle\frac{1}{\sqrt{\Omega_{\rm tot}}}\sum_{\mathbf{R}\alpha
a}\sum_{\begin{subarray}{c}|\mathbf{k}|<\Lambda_{c}\\\
\mathbf{Q}\in\mathcal{Q}_{l\eta}\end{subarray}}\left(w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}+i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}c_{\mathbf{k}a\eta
s}+w^{(\eta)}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}c^{\dagger}_{\mathbf{k}a\eta
s}f_{\mathbf{R}\alpha\eta s}\right)\bigg{]}$ (S133)

Since $w_{l\beta,\alpha}^{(\eta)}(\mathbf{r}-\mathbf{R})$ is a sharp Gaussian
function centered at $\mathbf{R}$, we can omit the overlap between Wannier
functions at different sites, i.e.,
$w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})w^{(\eta)}_{l\beta,\alpha^{\prime}}(\mathbf{r}-\mathbf{R}^{\prime})$
for $\mathbf{R}\neq\mathbf{R}^{\prime}$. From Eqs. S53 and S54 we find
$\sum_{l\beta}w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})w^{(\eta)}_{l\beta,\alpha^{\prime}}(\mathbf{r}-\mathbf{R})$
is zero if $\alpha\neq\alpha^{\prime}$. This property is guaranteed by the
$PC_{2z}T$ symmetry, which implies the constraint Eq. S49 and hence

$\displaystyle\sum_{l\beta}w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})w^{(\eta)}_{l\beta,\overline{\alpha}}(\mathbf{r}-\mathbf{R})=\sum_{l\beta}il\eta(-1)^{\alpha}w^{(\eta)*}_{-l\overline{\beta},\overline{\alpha}}(\mathbf{r}-\mathbf{R})(-i)l\eta(-1)^{\overline{\alpha}}w^{(\eta)}_{-l\overline{\beta},{\alpha}}(\mathbf{r}-\mathbf{R})$
$\displaystyle=$
$\displaystyle-\sum_{l\beta}w^{(\eta)*}_{-l\overline{\beta},\overline{\alpha}}(\mathbf{r}-\mathbf{R})w^{(\eta)}_{-l\overline{\beta},{\alpha}}(\mathbf{r}-\mathbf{R})=-\sum_{l\beta}w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})w^{(\eta)}_{l\beta,\overline{\alpha}}(\mathbf{r}-\mathbf{R})\
=0\ .$ (S134)

We can understand this property as a consequence of the (single-particle)
anti-unitary particle-hole symmetry ${\mathcal{P}}=PC_{2z}T$. According to the
representations in Eq. S36, $D^{f}(P)=i\sigma_{z}\tau_{z}$,
$D^{f}(C_{2z}T)=\sigma_{x}\tau_{0}$, there is
$D^{f}({\mathcal{P}})=-\sigma_{y}\tau_{z}$. Since ${\mathcal{P}}^{2}=-1$ and
$|W_{\mathbf{R}2\eta s}\rangle=-i\eta{\mathcal{P}}|W_{\mathbf{R}1\eta
s}\rangle$ there must be $\langle W_{\mathbf{R}2\eta s}|O|W_{\mathbf{R}1\eta
s}\rangle=0$ for any hermitian operator that commutes with ${\mathcal{P}}$
(Kramers theorem). Section S3.1 can be obtained if we choose $O$ as the
density operator at $\mathbf{r}$. According to the constraints Eqs. S46 and
S49, the density of the Wannier functions

$n_{f}(\mathbf{r})=\sum_{l\beta}|w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})|^{2}$
(S135)

does not depend on $\alpha$ and $\eta$. Therefore, we can simplify the density
operator as

$\displaystyle\hat{\rho}(\mathbf{r})=\sum_{\eta
s}\bigg{[}\sum_{\mathbf{R}\alpha}n_{f}(\mathbf{r}-\mathbf{R})f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}+\frac{1}{\Omega_{\rm tot}}\sum_{l\beta
aa^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}\end{subarray}}\sum_{\mathbf{Q}\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}}\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})\widetilde{u}^{(\eta)}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k}^{\prime})c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta s}$
$\displaystyle+\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\begin{subarray}{c}l\beta\\\ \mathbf{R}\alpha
a\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}|<\Lambda_{c}\\\
\mathbf{Q}\in\mathcal{Q}_{l\eta}\end{subarray}}\left(w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}+i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}c_{\mathbf{k}a\eta
s}+w^{(\eta)}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}c^{\dagger}_{\mathbf{k}a\eta
s}f_{\mathbf{R}\alpha\eta s}\right)\bigg{]}\ .$ (S136)

We now consider the interaction Hamiltonian of the form

$\hat{H}_{I}=\frac{1}{2}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}V(\mathbf{r}_{1}-\mathbf{r}_{2}):\hat{\rho}(\mathbf{r}_{1})::\hat{\rho}(\mathbf{r}_{2}):$
(S137)

where $V$ is the double-gate-screened Coulomb interaction (will be defined in
next paragraph) and $:\hat{\rho}:=\hat{\rho}-\langle
G_{0}|\hat{\rho}|G_{0}\rangle$ is the normal ordered density operator with
respect to a state $|G_{0}\rangle$ at charge neutrality point. We assume
$|G_{0}\rangle$ is such a state that

$\langle
G_{0}|c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}c_{\mathbf{k}^{\prime},\mathbf{Q}^{\prime},\alpha^{\prime},\eta^{\prime},s^{\prime}}|G_{0}\rangle=\frac{1}{2}\delta_{\mathbf{k}\mathbf{k}^{\prime}}\delta_{\mathbf{Q}\mathbf{Q}^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}\
.$ (S138)

One can verify

$\langle G_{0}|f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}|G_{0}\rangle=\frac{1}{2}\delta_{\mathbf{R}\mathbf{R}^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}},\qquad\langle
G_{0}|c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta^{\prime}s^{\prime}}|G_{0}\rangle=\frac{1}{2}\delta_{\mathbf{k}\mathbf{k}^{\prime}}\delta_{aa^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}\
,$ (S139) $\langle G_{0}|f_{\mathbf{R}\alpha\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta^{\prime}s^{\prime}}|G_{0}\rangle=\langle
G_{0}|c_{\mathbf{k}a\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}|G_{0}\rangle=0\
.$ (S140)

Therefore, we can write the normal ordered density operator as

$\displaystyle:\hat{\rho}(\mathbf{r}):=\sum_{\eta
s}\bigg{[}\sum_{\mathbf{R}\alpha}n_{f}(\mathbf{r}-\mathbf{R}):f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}:+\frac{1}{\Omega_{\rm tot}}\sum_{l\beta
aa^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}\\\
\mathbf{Q}\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta}\end{subarray}}e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}}\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})\widetilde{u}^{(\eta)}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k}^{\prime}):c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta s}:$
$\displaystyle+\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\begin{subarray}{c}l\beta\\\ \mathbf{R}\alpha
a\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}|<\Lambda_{c}\\\
\mathbf{Q}\in\mathcal{Q}_{l\eta}\end{subarray}}\left(w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}+i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}c_{\mathbf{k}a\eta
s}+w^{(\eta)}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})\widetilde{u}^{(\eta)*}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k})e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}c^{\dagger}_{\mathbf{k}a\eta
s}f_{\mathbf{R}\alpha\eta s}\right)\bigg{]}\ .$ (S141)

where

$:f_{\mathbf{R}\alpha\eta s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}:=f_{\mathbf{R}\alpha\eta s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}-\frac{1}{2},\qquad:c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta s}:=c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta
s}-\frac{1}{2}\delta_{aa^{\prime}}\ ,$ (S142) $:c^{\dagger}_{\mathbf{k}\eta
s}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}:=c^{\dagger}_{\mathbf{k}\eta
s}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}},\qquad:f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}c^{\dagger}_{\mathbf{k}\eta
s}:=f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}c^{\dagger}_{\mathbf{k}\eta
s}\ .$ (S143)

The double-gate-screened Coulomb interaction is given by

$V(\mathbf{r})=U_{\xi}\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{\sqrt{(\mathbf{r}/\xi)^{2}+n^{2}}}$
(S144)

with $\xi$ being distance between the two gates,
$U_{\xi}=e^{2}/(4\pi\epsilon\xi)$, and $\epsilon\approx 6$ is the dielectric
constant. For $\xi=10$nm there is $U_{\xi}=24$meV. For convenience we
introduce the Fourier transformation of $V(\mathbf{r})$ as

$V(\mathbf{r})=\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}V(\mathbf{q})e^{-i\mathbf{q}\cdot\mathbf{r}},\qquad
V(\mathbf{q})={\pi\xi^{2}U_{\xi}}\cdot\frac{\tanh(\xi|\mathbf{q}|/2)}{\xi|\mathbf{q}|/2}\
.$ (S145)

For simplicity, we denote the four terms in Section S3.1 as $:\rho_{ff}:$,
$:\rho_{cc}:$, $:\rho_{fc}:$, $:\rho_{cf}:$, respectively. Their explicit
forms are given by

$:\rho_{ff}:=\sum_{\eta
s}\sum_{\mathbf{R}\alpha}n_{f}(\mathbf{r}-\mathbf{R}):f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}:\ ,$ (S146)
$:\rho_{cc}:=\frac{1}{\Omega_{\rm tot}}\sum_{\eta s}\sum_{l\beta
aa^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}\\\
\mathbf{Q}\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta}\end{subarray}}e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}}\widetilde{u}^{(\eta)*}_{\mathbf{Q}\beta,a}(\mathbf{k})\widetilde{u}^{(\eta)}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k}^{\prime}):c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta s}:\ ,$ (S147)
$:\rho_{fc}:=\frac{1}{\sqrt{\Omega_{\rm tot}}}\sum_{\eta
s}\sum_{\begin{subarray}{c}l\beta\\\ \mathbf{R}\alpha
a\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}|<\Lambda_{c}\\\
\mathbf{Q}\in\mathcal{Q}_{l\eta}\end{subarray}}w^{(\eta)*}_{l\beta,\alpha}(\mathbf{r}-\mathbf{R})\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}+i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}c_{\mathbf{k}a\eta s},\qquad:\rho_{cf}:=(:\rho_{fc}:)^{\dagger}\
.$ (S148)

There are four diagonal terms in the interaction:

  1. 1.

The term
$:\rho_{ff}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{ff}(\mathbf{r}_{2}):$
is a density-density interaction of the local $f$-orbitals. As will be shown
in Section S3.2, the on-site repulsion ($U_{1}$) is the largest energy scale
of the problem.

  2. 2.

The term
$:\rho_{cc}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{cc}(\mathbf{r}_{2}):$
is the density-density interaction of the conduction bands.

  3. 3.

The term
$:\rho_{fc}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{fc}(\mathbf{r}_{2}):$
creates two particles in the local orbitals and two holes in the conduction
bands. Exciting particles or holes in the local orbitals must overcome a large
repulsion energy of the $f$-electrons. Thus this term is a high energy
process. We will omit this term in the Hartree-Fock calculation.

  4. 4.

The term
$:\rho_{cf}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{cf}(\mathbf{r}_{2}):$
creates two particles in the conduction bands and two holes in the local
orbitals.

There are also $\begin{pmatrix}4\\\ 2\end{pmatrix}=6$ off-diagonal terms in
the interaction:

  1. 5.

The term
$:\rho_{ff}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{cc}(\mathbf{r}_{2}):+(\leftrightarrow)$
is the density-density interaction between the local orbitals and the
conduction bands. This term will effectively shift the energy of conduction
bands if $:\rho_{ff}:\neq 0$.

  2. 6.

The term
$:\rho_{ff}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{fc}(\mathbf{r}_{2}):+(\leftrightarrow)$
creates two particles plus one hole in the local orbitals and a hole in the
conduction bands. As shown in Section S3.5, this term is guaranteed to be
vanishing by the symmetries.

  3. 7.

The term
$:\rho_{ff}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{cf}(\mathbf{r}_{2}):+(\leftrightarrow)$
creates two holes plus a particle in the local orbitals and a particle in the
conduction bands. As the 6th term, this term is guaranteed to be vanishing by
the symmetries.

  4. 8.

The term
$:\rho_{cc}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{fc}(\mathbf{r}_{2}):+(\leftrightarrow)$
creates a particle in the local orbitals and two holes plus one particle in
the conduction bands. Exciting particles or holes in the local orbitals must
overcome a large repulsion energy of the $f$-electrons. Thus this term is a
high energy process. This term is also weak because the hybridization, i.e.,
$\langle f^{\dagger}c\rangle$, is small, as discussed in Section S4. We will
omit this term in the Hartree-Fock calculation.

  5. 9.

The term
$:\rho_{cc}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{cf}(\mathbf{r}_{2}):+(\leftrightarrow)$
creates a hole in the local orbitals and two particles plus one hole in the
conduction bands. Exciting particles or holes in the local orbitals must
overcome a large repulsion energy of the $f$-electrons. Thus this term is a
high energy process. This term is also weak because the hybridization, i.e.,
$\langle c^{\dagger}f\rangle$, is small, as discussed in Section S4. We will
omit this term in the Hartree-Fock calculation.

  6. 10.

The term
$:\rho_{fc}(\mathbf{r}_{1}):V(\mathbf{r}_{1}-\mathbf{r}_{2}):\rho_{cf}(\mathbf{r}_{2}):+(\leftrightarrow)$
is an exchange interaction between local orbitals and conduction bands.

We will derive their explicit forms and calculate their interaction strengths
in Sections S3.2, S3.5, S3.3 and S3.4.

###  S3.2 Density-density terms (the first and fifth terms)

####  S3.2.1 The density-density terms of local orbitals (the first term)

We can write the density-density interaction of the local orbitals, i.e., the
first term discussed in Section S3.1, as

$\displaystyle{\hat{H}}_{U}=$ $\displaystyle\frac{1}{2}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}V(\mathbf{r}_{1}-\mathbf{r}_{2}):\hat{\rho}_{ff}(\mathbf{r}_{1})::\hat{\rho}_{ff}(\mathbf{r}_{2}):=\frac{1}{2}\sum_{\mathbf{R}\mathbf{R}^{\prime}}\sum_{\alpha\eta
s}\sum_{\alpha^{\prime}\eta^{\prime}s^{\prime}}U(\mathbf{R}-\mathbf{R}^{\prime}):f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}::f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}:\
,$ (S149)

where

$U(\mathbf{R})=\int d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}\
V(\mathbf{r}_{1}-\mathbf{r}_{2}-\mathbf{R})n_{f}(\mathbf{r}_{1})n_{f}(\mathbf{r}_{2})\
.$ (S150)

We will not apply the above equation for practical calculation of
$U(\mathbf{R})$ because $V(\mathbf{r})$ is divergent when $\mathbf{r}$
approaches zero. Instead we calculate $U(\mathbf{R})$ in momentum space.
According to the transformation (Eq. S44)

$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})=\frac{1}{\sqrt{N\Omega_{\rm
tot}}}\sum_{\mathbf{k}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{i(\eta\Delta\mathbf{K}_{l}+\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})\
,$ (S151)

where $\Omega_{\rm tot}=N\Omega_{0}$, we can express the density function (Eq.
S135) as

$n_{f}(\mathbf{r})=\frac{1}{N^{2}\Omega_{0}}\sum_{\mathbf{k}\mathbf{k}^{\prime}}\sum_{l\beta}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta}}e^{i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}}\widetilde{v}^{(\eta)*}_{\mathbf{Q}^{\prime}\beta,\alpha}(\mathbf{k}^{\prime})\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})\
.$ (S152)

We define the Fourier transformation of the density function as

$n_{f}(\mathbf{r})=\frac{1}{N\Omega_{0}}\sum_{\mathbf{q},\mathbf{G}}n_{f}(\mathbf{q}+\mathbf{G})e^{-i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}},\qquad
n_{f}(\mathbf{q}+\mathbf{G})=\frac{1}{N}\sum_{\mathbf{k}}\sum_{\mathbf{Q}\beta}\widetilde{v}^{(\eta)*}_{\mathbf{Q}-\mathbf{G}\beta,\alpha}(\mathbf{k}+\mathbf{q})\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})\
.$ (S153)

The function $U(\mathbf{R})$ can be calculated as

$\displaystyle U(\mathbf{R})=$ $\displaystyle\frac{1}{N^{2}\Omega_{0}^{2}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}\
V(\mathbf{r}_{1}-\mathbf{r}_{2}-\mathbf{R})\sum_{\mathbf{q}\mathbf{G}}\sum_{\mathbf{q}^{\prime}\mathbf{G}^{\prime}}e^{-i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}_{1}-i(\mathbf{q}^{\prime}+\mathbf{G}^{\prime})\cdot\mathbf{r}_{2}}n_{f}(\mathbf{q}+\mathbf{G})n_{f}(\mathbf{q}^{\prime}+\mathbf{G}^{\prime})$
$\displaystyle=$ $\displaystyle\frac{1}{N^{2}\Omega_{0}^{2}}\int
d^{2}\mathbf{r}_{2}\
\sum_{\mathbf{q}\mathbf{G}}\sum_{\mathbf{q}^{\prime}\mathbf{G}^{\prime}}V(\mathbf{q}+\mathbf{G})e^{-i(\mathbf{q}+\mathbf{G}+\mathbf{q}^{\prime}+\mathbf{G}^{\prime})\cdot\mathbf{r}_{2}-i(\mathbf{q}+\mathbf{G})\cdot\mathbf{R}}n_{f}(\mathbf{q}+\mathbf{G})n_{f}(\mathbf{q}^{\prime}+\mathbf{G}^{\prime})$
$\displaystyle=$
$\displaystyle\frac{1}{N\Omega_{0}}\sum_{\mathbf{q}\mathbf{G}}V(\mathbf{q}+\mathbf{G})e^{-i\mathbf{q}\cdot\mathbf{R}}n_{f}(\mathbf{q}+\mathbf{G})n_{f}(-\mathbf{q}-\mathbf{G})\
.$ (S154)

For $w_{0}/w_{1}=0.8$, we find that the onsite, nearest neighbor, next-nearest
neighbor, and third nearest neighbor interactions are

$U(0)=57.95\mathrm{meV},\qquad U(\mathbf{a}_{M1})=2.157\mathrm{meV},\qquad
U(\mathbf{a}_{M1}-\mathbf{a}_{M2})=0.1144\mathrm{meV},\qquad
U(2\mathbf{a}_{M1})=0.04145\mathrm{meV}\ ,$ (S155)

respectively. We can see that the interaction decays quickly as $\mathbf{R}$
and the onsite term ($U_{1}=U(0)$) is much larger than the other terms. In
fact, we find that $U_{1}$ is the largest energy scale in this problem. Thus
the onsite interaction will dominate the correlation physics. Nevertheless,
the $U(\mathbf{R})$ ($\mathbf{R}\neq 0$) terms will contribute to a chemical
potential term (for the local orbitals) in the $U_{1}\to\infty$ limit. Suppose
$U_{1}$ has frozen the charge fluctuation of local orbitals, i.e.,
$\forall\mathbf{R}$, $\sum_{\alpha\eta s}\langle f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}\rangle=N_{f}$ is constant, then all
terms in Eq. S149 with $\mathbf{R}\neq\mathbf{R}^{\prime}$ can be decoupled in
the Hartree channel as

$\displaystyle\frac{1}{2}\sum_{\mathbf{R}\neq\mathbf{R}^{\prime}}\sum_{\alpha\eta
s}\sum_{\alpha^{\prime}\eta^{\prime}s^{\prime}}U(\mathbf{R}-\mathbf{R}^{\prime}):f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}::f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}:$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\mathbf{R}\neq\mathbf{R}^{\prime}}\sum_{\alpha\eta
s}\sum_{\alpha^{\prime}\eta^{\prime}s^{\prime}}U(\mathbf{R}-\mathbf{R}^{\prime})\left(2\langle:f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}:\rangle:f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}:-\langle:f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}:\rangle\langle:f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}:\rangle\right)$
$\displaystyle=$ $\displaystyle\left((N_{f}-4)\sum_{\mathbf{R}^{\prime}\neq
0}U(\mathbf{R}^{\prime})\right)\sum_{\mathbf{R}}\sum_{\alpha\eta
s}f_{\mathbf{R}\alpha\eta s}^{\dagger}f_{\mathbf{R}\alpha\eta s}+const.$
(S156)

In order to capture the physics of chemical potential shift in an as simple as
possible model, we consider the approximate Hamiltonian of the form

${\hat{H}}_{U}=\frac{U_{1}}{2}\sum_{\mathbf{R}}\sum_{\alpha\eta
s}\sum_{\alpha^{\prime}\eta^{\prime}s^{\prime}}:f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}::f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}:+\frac{U_{2}}{2}\sum_{\langle\mathbf{R}\mathbf{R}^{\prime}\rangle}\sum_{\alpha\eta
s}\sum_{\alpha^{\prime}\eta^{\prime}s^{\prime}}:f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}::f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\alpha^{\prime}\eta^{\prime}s^{\prime}}:\
,$ (S157)

where $\langle\mathbf{R}\mathbf{R}^{\prime}\rangle$ only sums over nearest
neighbor pairs on the triangular lattice. Rather than naively choosing
$U_{2}=U(\mathbf{a}_{M1})$, we choose $U_{2}$ as such a value that it gives
the same chemical potential shift as Section S3.2.1. In other words, we absorb
all the chemical potential shift effects coming from the $U(\mathbf{R})$
($\mathbf{R}\neq 0$) terms into the nearest neighbor repulsion. $U_{2}$ can be
calculated as

$U_{2}=\frac{1}{6}\sum_{\mathbf{R}\neq
0}U(\mathbf{R})=\frac{1}{6}\left(-U_{1}+\frac{1}{\Omega_{0}}\sum_{\mathbf{G}}V(\mathbf{G})n_{f}(\mathbf{G})n_{f}(-\mathbf{G})\right)\
.$ (S158)

For $w_{0}/w_{1}=0.8$ we find $U_{2}=2.329$meV, which is only slightly larger
than $U(\mathbf{a}_{M1})$.

$w_{0}/w_{1}$ | $U_{1}$ | $U_{2}$ | $W_{1}=W_{2}$ | $W_{3}=W_{4}$ | $J$ | $K$
---|---|---|---|---|---|---
1.0 | 69.35 | 1.914 | 43.46 | 52.90 | 14.48 | 5.122
0.9 | 64.56 | 2.006 | 43.84 | 51.37 | 15.40 | 4.958
0.8 | 57.95 | 2.329 | 44.03 | 50.20 | 16.38 | 4.887
0.7 | 51.72 | 2.656 | 44.05 | 49.33 | 18.27 | 4.897
0.6 | 45.02 | 3.103 | 44.03 | 48.73 | 20.75 | 4.975
0.5 | 37.61 | 3.709 | 44.06 | 48.35 | 23.75 | 5.111
Table S6: Parameters of the interaction Hamiltonian. Parameters of the BM
model used to obtain this table are $v_{F}=5.944\mathrm{eV\cdot\mathring{A}}$,
$|\Delta\mathbf{K}|=1.703\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$,
$\theta=1.05^{\circ}$. Parameters of the double-gate-screened interaction used
to obtain this table are $U_{\xi}=24{\rm meV}$, $\xi=10{\rm nm}$.

####  S3.2.2 The density-density terms between local orbitals and conduction
bands (the fifth term)

Second we consider the interaction between local orbital density and
conduction electron density

$\displaystyle{\hat{H}}_{W}=$ $\displaystyle\frac{1}{\Omega_{\rm tot}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}\sum_{\mathbf{R}\alpha\eta_{1}s_{1}}:f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta_{1}s_{1}}:n_{f}(\mathbf{r}_{1}-\mathbf{R})V(\mathbf{r}_{1}-\mathbf{r}_{2})\sum_{|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}}\sum_{\eta_{2}s_{2}aa^{\prime}}\sum_{l\beta}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta_{2}}}$
$\displaystyle\times\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\beta,a}(\mathbf{k})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k}^{\prime})e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}_{2}}:c_{\mathbf{k}a\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta_{2}s_{2}}:$
(S159)

Applying the transformation (Eq. S145)

$\int
d^{2}\mathbf{r}_{2}V(\mathbf{r}_{1}-\mathbf{r}_{2})e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}_{2}}=V(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}_{1}}$
(S160)

and the variable substitution $\mathbf{r}_{1}\to\mathbf{r}_{1}+\mathbf{R}$, we
can rewrite ${\hat{H}}_{W}$ as

$\displaystyle{\hat{H}}_{W}=$ $\displaystyle\frac{1}{{\Omega_{\rm tot}}}\int
d^{2}\mathbf{r}_{1}\sum_{\mathbf{R}\alpha\eta_{1}s_{1}}:f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta_{1}s_{1}}:n_{f}(\mathbf{r}_{1})\sum_{|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}}\sum_{\eta_{2}s_{2}aa^{\prime}l}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta_{2}}}e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot(\mathbf{r}_{1}+\mathbf{R})}$
$\displaystyle\times
V(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\beta,a}(\mathbf{k})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k}^{\prime}):c_{\mathbf{k}a\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta_{2}s_{2}}:$
(S161)

We introduce the integral

$\displaystyle
X_{aa^{\prime}}^{(\eta_{2})}(\mathbf{k},\mathbf{k}^{\prime})=\frac{1}{\Omega_{0}}\int
d^{2}\mathbf{r}\
n_{f}(\mathbf{r})\sum_{l\beta}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta_{2}}}V(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})e^{-i(\mathbf{k}-\mathbf{Q}-\mathbf{k}^{\prime}+\mathbf{Q}^{\prime})\cdot\mathbf{r}}\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\beta,a}(\mathbf{k})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(\mathbf{k}^{\prime})$
(S162)

such that the interaction can be written as
($e^{i(\mathbf{Q}-\mathbf{Q}^{\prime})\cdot\mathbf{R}}=1$)

${\hat{H}}_{W}=\frac{1}{N}\sum_{\mathbf{R}\alpha\eta_{1}s_{1}}\sum_{|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}}\sum_{\eta_{2}s_{2}aa^{\prime}}X_{aa^{\prime}}^{(\eta_{2})}(\mathbf{k},\mathbf{k}^{\prime})e^{-i(\mathbf{k}-\mathbf{k}^{\prime})\cdot\mathbf{R}}:f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta_{1}s_{1}}::c_{\mathbf{k}a\eta_{s}s_{2}}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta_{s}s_{2}}:$
(S163)

The conduction bands have low energy states around $\mathbf{k}\approx 0$. For
the low energy states we can approximate the integral as
$X_{aa^{\prime}}^{(\eta_{2})}(\mathbf{k},\mathbf{k}^{\prime})\approx
X_{aa^{\prime}}^{(\eta_{2})}(0,0)$

$\displaystyle X_{aa^{\prime}}^{(\eta_{2})}(0,0)=$
$\displaystyle\frac{1}{\Omega_{0}}\int
d^{2}\mathbf{r}_{1}n_{f}(\mathbf{r})\sum_{l\beta}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta_{2}}}V(\mathbf{Q}-\mathbf{Q}^{\prime})e^{i(\mathbf{Q}-\mathbf{Q}^{\prime})\cdot\mathbf{r}}\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\beta,a}(0)\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}^{\prime}\beta,a^{\prime}}(0)$
$\displaystyle=$
$\displaystyle\frac{1}{\Omega_{0}}\sum_{l\beta}\sum_{\mathbf{Q}\mathbf{G}}n_{f}(\mathbf{G})V(\mathbf{G})\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\beta,a}(0)\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}-\mathbf{G}\beta,a^{\prime}}(0)\
.$ (S164)

Here $\mathbf{Q}$ indexes all vectors in
$\mathcal{Q}_{+}\oplus\mathcal{Q}_{-}$ and $\mathbf{G}$ indexes reciprocal
lattices. Under the crystalline symmetries the function
$n_{f}(\mathbf{G})V(\mathbf{G})$ forms an identity representation, thus the
summation over $\mathbf{Q},\mathbf{G},\beta$ is non-vanishing only if
$\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\beta,a}(0)\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}-\mathbf{G}\beta,a^{\prime}}(0)$
also form an identity representation of the valley-preserving symmetries. To
be concrete, because
$n_{f}(g\mathbf{G})V(g\mathbf{G})=n_{f}(\mathbf{G})V(\mathbf{G})$ for any
crystalline symmetry $g$ that preserves the valley, we can rewrite the $X$
matrix as

$\displaystyle
X_{aa^{\prime}}^{(\eta_{2})}(0,0)=\frac{1}{\Omega_{0}}\sum_{l\beta}\sum_{\mathbf{Q}\mathbf{G}}n_{f}(g\mathbf{G})V(g\mathbf{G})\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\beta,a}(0)\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}-\mathbf{G}\beta,a^{\prime}}(0)\
,$ (S165)

Relabeling $\mathbf{Q}\to g^{-1}\mathbf{Q}$, $\mathbf{G}\to g^{-1}\mathbf{G}$,
we obtain

$\displaystyle X_{aa^{\prime}}^{(\eta_{2})}(0,0)=$
$\displaystyle\frac{1}{\Omega_{0}}\sum_{l\beta}\sum_{\mathbf{Q}\mathbf{G}}n_{f}(\mathbf{G})V(\mathbf{G})\widetilde{u}^{(\eta_{2})*}_{g^{-1}\mathbf{Q}\beta,a}(0)\widetilde{u}^{(\eta_{2})}_{g^{-1}\mathbf{Q}-g^{-1}\mathbf{G}\beta,a^{\prime}}(0)\
.$ (S166)

For simplicity, we denote the symmetry operator matrices in Eqs. S16 and S17
as
$D_{\mathbf{Q}\alpha\eta,\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime}}(g)=\delta_{\mathbf{Q},g\mathbf{Q}^{\prime}}\delta_{\eta\eta^{\prime}}\overline{D}_{\alpha,\alpha^{\prime}}^{(\eta)}(g)$
for $g$ that preserve the valley. Since
$\overline{D}_{\alpha,\alpha^{\prime}}^{(\eta)}(g)$’s are unitary, i.e.,
$\sum_{\alpha}\overline{D}_{\alpha,\beta}^{(\eta)*}(g)\overline{D}_{\alpha,\beta^{\prime}}^{(\eta)}(g)=\delta_{\beta\beta^{\prime}}$,
there are

$\displaystyle\sum_{\beta\mathbf{Q}}\widetilde{u}^{(\eta_{2})*}_{g^{-1}\mathbf{Q}\beta,a}(0)\widetilde{u}^{(\eta_{2})}_{g^{-1}\mathbf{Q}-g^{-1}\mathbf{G}\beta,a^{\prime}}(0)=\sum_{\beta\beta^{\prime}\mathbf{Q}}\sum_{\alpha}\overline{D}_{\alpha,\beta}^{(\eta_{2})*}(g)\widetilde{u}^{(\eta_{2})*}_{g^{-1}\mathbf{Q}\beta,a}(0)\overline{D}_{\alpha,\beta^{\prime}}^{(\eta_{2})}\widetilde{u}^{(\eta_{2})}_{g^{-1}\mathbf{Q}-g^{-1}\mathbf{G}\beta^{\prime},a^{\prime}}(0)$
(S167)

Since the vectors $\widetilde{u}^{(\eta_{2})}$ form irreducible
representations given by Section S2.2, i.e.,

$\sum_{\beta}\overline{D}_{\alpha,\beta^{\prime}}^{(\eta_{2})}(g)\widetilde{u}^{(\eta_{2})}_{g^{-1}\mathbf{Q}\beta^{\prime},a}(0)=\sum_{\beta\mathbf{Q}^{\prime}}{D}_{\mathbf{Q}\alpha\eta_{2},\mathbf{Q}^{\prime}\beta^{\prime}\eta_{2}}\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}^{\prime}\beta^{\prime},a}(0)=\sum_{b}\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}\alpha,a^{\prime}}(0)D_{b\eta_{2},a\eta_{2}}^{(c)}(g)$
(S168)

there is

$\displaystyle\sum_{\beta\mathbf{Q}}\widetilde{u}^{(\eta_{2})*}_{g^{-1}\mathbf{Q}\beta,a}(0)\widetilde{u}^{(\eta_{2})}_{g^{-1}\mathbf{Q}-g^{-1}\mathbf{G}\beta,a^{\prime}}(0)=\sum_{bb^{\prime}}\sum_{\alpha\mathbf{Q}}\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}\alpha,b}(0)D_{b\eta_{2},a\eta_{2}}^{(c)*}(g)\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}-\mathbf{G}\alpha,b^{\prime}}(0)D_{b^{\prime}\eta_{2},a^{\prime}\eta_{2}}^{(c)}(g)$
(S169)

Substituting this relation into Eq. S166, we obtain

$\displaystyle X_{aa^{\prime}}^{(\eta_{2})}(0,0)=$
$\displaystyle\sum_{bb^{\prime}}D_{b\eta_{2},a\eta_{2}}^{(c)*}(g)X_{bb^{\prime}}^{(\eta_{2})}(0,0)D_{b^{\prime}\eta_{2},a^{\prime}\eta_{2}}^{(c)}(g)\
.$ (S170)

In other words, $X$ commutes with all $D^{(c)}(g)$ for crystalline symmetries
$g$ that preserve the valley. These symmetries are generated by $g=C_{3z}$,
$C_{2x}$, $C_{2z}T$. Since $D^{(c)}(g)$ consists of three different
irreducible representations of the crystalline symmetry, i.e.,
$\Gamma_{3}\oplus\Gamma_{1}\oplus\Gamma_{2}$, $X$ must take the following form
in the $\Gamma_{3}\oplus\Gamma_{1}\oplus\Gamma_{2}$ basis according to the
Schur’s lemma

$X^{(\eta)}=\begin{pmatrix}W_{1}^{(\eta)}&0&0&0\\\ 0&W_{1}^{(\eta)}&0&0\\\
0&0&W_{3}^{(\eta)}&W_{3}^{\prime(\eta)}\\\
0&0&W_{3}^{\prime(\eta)}&W_{3}^{(\eta)}\end{pmatrix}$ (S171)

The off diagonal term $W_{3}^{\prime(\eta)}$ is allowed by the crystalline
symmetries because the third and fourth components are mixed states of the
$\Gamma_{1}\oplus\Gamma_{2}$ basis. To be concrete, in the
$\Gamma_{1}\oplus\Gamma_{2}$ subspace in each valley $\eta$, $C_{3z}$,
$C_{2z}$, and $C_{2z}T$ are represented by $\sigma_{0}$, $\sigma_{x}$,
$\sigma_{x}K$ (with $K$ being complex conjugation), respectively, and hence
$W_{3}^{(\eta)}\sigma_{0}+W_{3}^{\prime(\eta)}\sigma_{x}$ commutes with all of
them. However, $W_{3}^{\prime(\eta)}$ is forbidden by the particle-hole
symmetry $P$. Due to $[X^{(\eta)},D^{(c,\eta)}(P)]=0$ and
$D^{(c,\eta)}_{aa^{\prime}}(P)=D^{(c)}_{a\eta,a^{\prime}\eta}(g)=\eta[(-i\sigma_{z})\oplus(-i\sigma_{z})]_{aa^{\prime}}$,
$W_{3}^{\prime(\eta)}$ must vanish and hence $X^{(\eta)}$ is diagonal. Since
$n_{f}(\mathbf{r})$ is a real function, there must be
$n_{f}(\mathbf{Q}-\mathbf{Q}^{\prime})=n_{f}^{*}(\mathbf{Q}^{\prime}-\mathbf{Q})$
and hence $X^{(\eta)}_{aa^{\prime}}(0,0)$ is a hermitian matrix, implying that
$X^{(\eta)}_{aa}(0,0)$ are all real numbers. Due to the time-reversal
symmetry, $X^{(\eta)}_{aa^{\prime}}(0,0)$ does not depend on the $\eta$ index.
We hence write the $X$ matrix as
$X^{(\eta)}_{aa^{\prime}}(0,0)=W_{a}\delta_{aa^{\prime}}$, where
$W_{1}=W_{2}$, $W_{3}=W_{4}$. For $w_{0}/w_{1}=0.8$ we find $W_{1}=44.03$meV
and $W_{3}=50.20$meV. Values of $W_{a}$ at other $w_{0}/w_{1}$ are tabulated
in Table S6.

In summary, the density-density interaction between the local orbitals and the
conduction bands can be written as

${\hat{H}}_{W}=\frac{1}{N}\sum_{\mathbf{R}\alpha\eta_{1}s_{1}}\sum_{|\mathbf{k}|,|\mathbf{k}^{\prime}|<\Lambda_{c}}\sum_{\eta_{2}s_{2}a}W_{a}e^{-i(\mathbf{k}-\mathbf{k}^{\prime})\cdot\mathbf{R}}:f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta_{1}s_{1}}::c_{\mathbf{k}a\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}^{\prime}a\eta_{2}s_{2}}:$
(S172)

####  S3.2.3 Density-density interaction of the conduction bands (the second
term)

We third consider the density-density interaction of the conduction bands

$\displaystyle{\hat{H}}_{V}=$ $\displaystyle\frac{1}{2\Omega_{\rm
tot}^{2}}\sum_{\beta_{1}l_{1}\eta_{1}s_{1}}\sum_{\beta_{2}l_{2}\eta_{2}s_{2}}\sum_{a_{1}a_{1}^{\prime}a_{2}a_{2}^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|,|\mathbf{k}_{1}^{\prime}|<\Lambda_{c}\\\
\mathbf{Q}_{1}\mathbf{Q}_{1}^{\prime}\in\mathcal{Q}_{l_{1}\eta_{1}}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|,|\mathbf{k}_{2}^{\prime}|<\Lambda_{c}\\\
\mathbf{Q}_{2}\mathbf{Q}_{2}^{\prime}\in\mathcal{Q}_{l_{2}\eta_{2}}\end{subarray}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}V(\mathbf{r}_{1}-\mathbf{r}_{2})$
$\displaystyle\times
e^{-i(\mathbf{k}_{1}-\mathbf{Q}_{1}-\mathbf{k}_{1}^{\prime}+\mathbf{Q}_{1}^{\prime})\cdot\mathbf{r}_{1}}\widetilde{u}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}^{\prime}_{1}\beta_{1},a_{1}^{\prime}}(\mathbf{k}_{1}^{\prime})e^{-i(\mathbf{k}_{2}^{\prime}-\mathbf{Q}_{2}^{\prime}-\mathbf{k}_{2}+\mathbf{Q}_{2})\cdot\mathbf{r}_{2}}\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}_{2}^{\prime}\beta_{2},a_{2}^{\prime}}(\mathbf{k}_{2}^{\prime})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a_{2}}(\mathbf{k}_{2})$
$\displaystyle\times:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}::c_{\mathbf{k}_{2}^{\prime}a_{2}^{\prime}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:\
.$ (S173)

After integration over $\mathbf{r}_{1}$, there is

$\displaystyle{\hat{H}}_{V}=$ $\displaystyle\frac{1}{2\Omega_{\rm
tot}^{2}}\sum_{\beta_{1}l_{1}\eta_{1}s_{1}}\sum_{\beta_{2}l_{2}\eta_{2}s_{2}}\sum_{a_{1}a_{1}^{\prime}a_{2}a_{2}^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|,|\mathbf{k}_{1}^{\prime}|<\Lambda_{c}\\\
\mathbf{Q}_{1}\mathbf{Q}_{1}^{\prime}\in\mathcal{Q}_{l_{1}\eta_{1}}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|,|\mathbf{k}_{2}^{\prime}|<\Lambda_{c}\\\
\mathbf{Q}_{2}\mathbf{Q}_{2}^{\prime}\in\mathcal{Q}_{l_{2}\eta_{2}}\end{subarray}}\int
d^{2}\mathbf{r}_{2}V(\mathbf{k}_{1}^{\prime}-\mathbf{Q}_{1}^{\prime}-\mathbf{k}_{1}+\mathbf{Q}_{1})$
$\displaystyle\times
e^{-i(\mathbf{k}_{2}^{\prime}-\mathbf{Q}_{2}^{\prime}-\mathbf{k}_{2}+\mathbf{Q}_{2}-\mathbf{k}_{1}^{\prime}+\mathbf{Q}_{1}^{\prime}+\mathbf{k}_{1}-\mathbf{Q}_{1})\cdot\mathbf{r}_{2}}\widetilde{u}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}^{\prime}_{1}\beta_{1},a_{1}^{\prime}}(\mathbf{k}_{1}^{\prime})\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}_{2}^{\prime}\beta_{2},a_{2}^{\prime}}(\mathbf{k}_{2}^{\prime})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a_{2}}(\mathbf{k}_{2})$
$\displaystyle\times:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}::c_{\mathbf{k}_{2}^{\prime}a_{2}^{\prime}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:\
.$ (S174)

Integration over $\mathbf{r}_{2}$ leads to the momentum conservation
$\mathbf{k}_{1}^{\prime}-\mathbf{Q}_{1}^{\prime}-\mathbf{k}_{1}+\mathbf{Q}_{1}=\mathbf{k}_{2}^{\prime}-\mathbf{Q}_{2}^{\prime}-\mathbf{k}_{2}+\mathbf{Q}_{2}$.
Since $-\mathbf{Q}_{1}^{\prime}+\mathbf{Q}_{1}$ and
$-\mathbf{Q}_{2}^{\prime}+\mathbf{Q}_{2}$ are reciprocal lattices and
$\mathbf{k}_{1}^{\prime}$, $\mathbf{k}_{1}$, $\mathbf{k}_{2}$,
$\mathbf{k}_{2}^{\prime}$ are small momenta around $\Gamma_{M}$, there must be
$\mathbf{k}_{1}^{\prime}-\mathbf{k}_{1}=\mathbf{k}_{2}^{\prime}-\mathbf{k}_{2}$
and
$-\mathbf{Q}_{1}^{\prime}+\mathbf{Q}_{1}=-\mathbf{Q}_{2}^{\prime}+\mathbf{Q}_{2}$
separately. We introduce the momentum
$\mathbf{q}=\mathbf{k}_{1}^{\prime}-\mathbf{k}_{1}=\mathbf{k}_{2}^{\prime}-\mathbf{k}_{2}$,
the reciprocal lattice
$\mathbf{G}=-\mathbf{Q}_{1}^{\prime}+\mathbf{Q}_{1}=-\mathbf{Q}_{2}^{\prime}+\mathbf{Q}_{2}$
and rewrite the interaction as

$\displaystyle{\hat{H}}_{V}=$
$\displaystyle\frac{1}{2\Omega_{0}N}\sum_{\beta_{1}l_{1}\eta_{1}s_{1}}\sum_{\beta_{2}l_{2}\eta_{2}s_{2}}\sum_{a_{1}a_{1}^{\prime}a_{2}a_{2}^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|<\Lambda_{c}\\\
\mathbf{Q}_{1}\in\mathcal{Q}_{l_{1}\eta_{1}}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|<\Lambda_{c}\\\
\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}\end{subarray}}\sum_{\mathbf{G}}\sum_{\begin{subarray}{c}\mathbf{q}\\\
|\mathbf{k}_{1}+\mathbf{q}|,|\mathbf{k}_{2}+\mathbf{q}|<\Lambda_{c}\end{subarray}}V(\mathbf{q}+\mathbf{G})$
$\displaystyle\times\widetilde{u}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}-\mathbf{G}\beta_{1},a_{1}^{\prime}}(\mathbf{k}_{1}+\mathbf{q})\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}_{2}-\mathbf{G}\beta_{2},a_{2}^{\prime}}(\mathbf{k}_{2}+\mathbf{q})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a_{2}}(\mathbf{k}_{2}):c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}+\mathbf{q}a_{1}^{\prime}\eta_{1}s_{1}}::c_{\mathbf{k}_{2}+\mathbf{q}a_{2}^{\prime}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:\
.$ (S175)

We introduce the matrix

$X_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}a_{2}a_{2}^{\prime}}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})=\frac{1}{\Omega_{0}}\sum_{\mathbf{G}}V(\mathbf{q}+\mathbf{G})\langle\widetilde{u}^{(\eta_{1})}_{a_{1}}(\mathbf{k}_{1})|\widetilde{u}^{(\eta_{1})}_{a_{1}^{\prime}}(\mathbf{k}_{1}+\mathbf{q}+\mathbf{G})\rangle\langle\widetilde{u}^{(\eta_{2})}_{a_{2}^{\prime}}(\mathbf{k}_{2}+\mathbf{q}+\mathbf{G})|\widetilde{u}^{(\eta_{2})}_{a_{2}}(\mathbf{k}_{2})\rangle$
(S176)

such that the interaction can be simplified to

$\displaystyle{\hat{H}}_{V}=$
$\displaystyle\frac{1}{2N}\sum_{\eta_{1}s_{1}a_{1}a_{1}^{\prime}}\sum_{\eta_{2}s_{2}a_{2}a_{2}^{\prime}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}\sum_{\begin{subarray}{c}\mathbf{q}\\\
|\mathbf{k}_{1}+\mathbf{q}|,|\mathbf{k}_{2}+\mathbf{q}|<\Lambda_{c}\end{subarray}}X_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}a_{2}a_{2}^{\prime}}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})$
$\displaystyle\qquad\times:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}+\mathbf{q}a_{1}^{\prime}\eta_{1}s_{1}}::c_{\mathbf{k}_{2}+\mathbf{q}a_{2}^{\prime}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:\
.$ (S177)

For low energy states around $\mathbf{k}=0$, we can approximate the $X$ matrix
by replacing $\widetilde{u}^{(\eta)}_{a}(\mathbf{k})$ with
$\widetilde{u}^{(\eta)}_{a}(0)$

$X_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}a_{2}a_{2}^{\prime}}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})\approx
X_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}a_{2}a_{2}^{\prime}}(0,0;\mathbf{q})=\frac{1}{\Omega_{0}}\sum_{\mathbf{G}}V(\mathbf{q}+\mathbf{G})\langle\widetilde{u}^{(\eta_{1})}_{a_{1}}(0)|\widetilde{u}^{(\eta_{1})}_{a_{1}^{\prime}}(\mathbf{G})\rangle\langle\widetilde{u}^{(\eta_{2})}_{a_{2}^{\prime}}(\mathbf{G})|\widetilde{u}^{(\eta_{2})}_{a_{2}}(0)\rangle\
.$ (S178)

We find that the summation over $\mathbf{G}$ is dominated by the
$\mathbf{G}=0$ component. For example, the diagonal element
$X_{\eta_{1}11,\eta_{2}11}(0,0;0)$ is found to be 49.22meV, and the
$\mathbf{G}=0$ component is $\frac{1}{\Omega_{0}}V(0)=48.33$meV. Therefore, it
is reasonable to keep only the $\mathbf{G}=0$ component. Then we can write
${\hat{H}}_{V}$ as

${\hat{H}}_{V}=\frac{1}{2\Omega_{0}N}\sum_{\eta_{1}s_{1}a_{1}}\sum_{\eta_{2}s_{2}a_{2}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}\sum_{\begin{subarray}{c}\mathbf{q}\\\
|\mathbf{k}_{1}+\mathbf{q}|,|\mathbf{k}_{2}+\mathbf{q}|<\Lambda_{c}\end{subarray}}V(\mathbf{q}):c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}+\mathbf{q}a_{1}\eta_{1}s_{1}}::c_{\mathbf{k}_{2}+\mathbf{q}a_{2}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:\
.$ (S179)

###  S3.3 Exchange terms (the tenth term)

####  S3.3.1 Constraints of the coupling matrix

Now we consider the $:\hat{\rho}_{fc}:V:\hat{\rho}_{cf}:+(\leftrightarrow)$
interaction discussed in Section S3.1 (the tenth term)

$\displaystyle{\hat{H}}_{J}=$ $\displaystyle\frac{1}{2\Omega_{\rm
tot}}\sum_{\begin{subarray}{c}\beta_{1}l_{1}\eta_{1}s_{1}\\\
\beta_{2}l_{2}\eta_{2}s_{2}\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{R}_{1}\alpha_{1}\\\
\mathbf{R}_{2}\alpha_{2}\end{subarray}}\sum_{a_{1}a_{2}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|<\Lambda_{c}\\\
\mathbf{Q}_{1}\in\mathcal{Q}_{l_{1}\eta_{1}}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|<\Lambda_{c}\\\
\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}\end{subarray}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}\ V(\mathbf{r}_{1}-\mathbf{r}_{2})$
$\displaystyle\times
w^{(\eta_{1})*}_{l_{1}\beta_{1},\alpha_{1}}(\mathbf{r}_{1}-\mathbf{R}_{1})w^{(\eta_{2})}_{l_{2}\beta_{2},\alpha_{2}}(\mathbf{r}_{2}-\mathbf{R}_{2})\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})e^{i\eta_{1}\Delta\mathbf{K}_{l_{1}}\cdot\mathbf{R}_{1}+i(\mathbf{k}_{1}-\mathbf{Q}_{1})\cdot\mathbf{r}_{1}}\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}_{2}\beta_{2},a_{2}}(\mathbf{k}_{2})e^{-i\eta_{2}\Delta\mathbf{K}_{l_{2}}\cdot\mathbf{R}_{2}-i(\mathbf{k}_{2}-\mathbf{Q}_{2})\cdot\mathbf{r}_{2}}$
$\displaystyle\times\left(f_{\mathbf{R}_{1}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}_{2}\alpha_{2}\eta_{2}s_{2}}+c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}_{2}\alpha_{2}\eta_{2}s_{2}}f_{\mathbf{R}_{1}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\right)$
(S180)

Now we argue that, since the Wannier functions are well localized, the
$\mathbf{R}_{1}=\mathbf{R}_{2}$ terms will be much stronger than the
$\mathbf{R}_{1}\neq\mathbf{R}_{2}$ terms. In Section S3.2 we find that the
terms $\int d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}\
V(\mathbf{r}_{1}-\mathbf{r}_{2})n_{f}(\mathbf{r}_{1}-\mathbf{R}_{1})n_{f}(\mathbf{r}_{2}-\mathbf{R}_{2})$
with $\mathbf{R}_{1}=\mathbf{R}_{2}$ are about 30 times the terms with
$\mathbf{R}_{1},\mathbf{R}_{2}$ being nearest neighbors. Now
$n_{f}(\mathbf{r}-\mathbf{R})$ is replaced by $w(\mathbf{r}-\mathbf{R})$.
Since $w(\mathbf{r}-\mathbf{R})$ ($\sim\sqrt{n(\mathbf{r}-\mathbf{R})}$) has a
similar profile as $n_{f}(\mathbf{r}-\mathbf{R})$, the
$\mathbf{R}_{1}=\mathbf{R}_{2}$ terms will still be much stronger than the
$\mathbf{R}_{1}\neq\mathbf{R}_{2}$ terms. (We have also numerically verified
this argument. We find that, for $w_{0}/w_{1}=0.8$, the largest coupling for
$|\mathbf{R}_{1}-\mathbf{R}_{2}|=0$ and
$|\mathbf{R}_{1}-\mathbf{R}_{2}|=a_{M}$ are about 16meV and 1meV,
respectively.) Thus we can simplify this term as

$\displaystyle{\hat{H}}_{J}=$ $\displaystyle\frac{1}{2\Omega_{\rm
tot}}\sum_{\begin{subarray}{c}\beta_{1}l_{1}\eta_{1}s_{1}\\\
\beta_{2}l_{2}\eta_{2}s_{2}\end{subarray}}\sum_{\mathbf{R}\alpha_{1}\alpha_{2}}\sum_{a_{1}a_{2}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|<\Lambda_{c}\\\
\mathbf{Q}_{1}\in\mathcal{Q}_{l_{1}\eta_{1}}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|<\Lambda_{c}\\\
\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}\end{subarray}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}\ V(\mathbf{r}_{1}-\mathbf{r}_{2})$
$\displaystyle\times
w^{(\eta_{1})*}_{l_{1}\beta_{1},\alpha_{1}}(\mathbf{r}_{1}-\mathbf{R})w^{(\eta_{2})}_{l_{2}\beta_{2},\alpha_{2}}(\mathbf{r}_{2}-\mathbf{R})e^{i(\eta_{1}\Delta\mathbf{K}_{l_{1}}-\eta_{2}\Delta\mathbf{K}_{l_{2}})\cdot\mathbf{R}}\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})e^{i(\mathbf{k}_{1}-\mathbf{Q}_{1})\cdot\mathbf{r}_{1}}\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}_{2}\beta_{2},a_{2}}(\mathbf{k}_{2})e^{-i(\mathbf{k}_{2}-\mathbf{Q}_{2})\cdot\mathbf{r}_{2}}$
$\displaystyle\times\left(f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}+c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\right)$
(S181)

Substituting the Fourier transformation of the interaction into the above
equation, we obtain

$\displaystyle{\hat{H}}_{J}=$ $\displaystyle\frac{1}{2\Omega_{\rm
tot}}\sum_{\begin{subarray}{c}\beta_{1}l_{1}\eta_{1}s_{1}\\\
\beta_{2}l_{2}\eta_{2}s_{2}\end{subarray}}\sum_{\mathbf{R}\alpha_{1}\alpha_{2}}\sum_{a_{1}a_{2}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|<\Lambda_{c}\\\
\mathbf{Q}_{1}\in\mathcal{Q}_{l_{1}\eta_{1}}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|<\Lambda_{c}\\\
\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}\end{subarray}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\
V(\mathbf{q})e^{-i\mathbf{q}(\mathbf{r}_{1}-\mathbf{r}_{2})}$
$\displaystyle\times
w^{(\eta_{1})*}_{l_{1}\beta_{1},\alpha_{1}}(\mathbf{r}_{1}-\mathbf{R})w^{(\eta_{2})}_{l_{2}\beta_{2},\alpha_{2}}(\mathbf{r}_{2}-\mathbf{R})e^{i(\eta_{1}\Delta\mathbf{K}_{l_{1}}-\eta_{2}\Delta\mathbf{K}_{l_{2}})\cdot\mathbf{R}}\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})e^{i(\mathbf{k}_{1}-\mathbf{Q}_{1})\cdot\mathbf{r}_{1}}\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}_{2}\beta_{2},a_{2}}(\mathbf{k}_{2})e^{-i(\mathbf{k}_{2}-\mathbf{Q}_{2})\cdot\mathbf{r}_{2}}$
$\displaystyle\times\left(f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}+c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\right)$
(S182)

According to Eq. S45, we have

$\int
d^{2}\mathbf{r}_{2}w^{(\eta_{2})}_{l_{2}\beta_{2},\alpha_{2}}(\mathbf{r}_{2}-\mathbf{R})e^{-i(\mathbf{k}_{2}-\mathbf{Q}_{2}-\mathbf{q})\cdot\mathbf{r}_{2}}=\sqrt{\Omega_{0}}\
e^{-i(\mathbf{k}_{2}-\mathbf{Q}_{2}-\mathbf{q})\cdot\mathbf{R}}\widetilde{v}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},\alpha_{2}}(\mathbf{k}_{2}-\mathbf{q})\
.$ (S183) $\int
d^{2}\mathbf{r}_{1}w^{(\eta_{1})*}_{l_{1}\beta_{1},\alpha_{1}}(\mathbf{r}_{1}-\mathbf{R})e^{i(\mathbf{k}_{1}-\mathbf{Q}_{1}-\mathbf{q})\cdot\mathbf{r}_{1}}=\sqrt{\Omega_{0}}\
e^{i(\mathbf{k}_{1}-\mathbf{Q}_{1}-\mathbf{q})\cdot\mathbf{R}}\widetilde{v}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},\alpha_{1}}(\mathbf{k}_{1}-\mathbf{q})\
.$ (S184)

Therefore, we find

$\displaystyle{\hat{H}}_{J}=$
$\displaystyle\frac{1}{2N}\sum_{\begin{subarray}{c}\beta_{1}l_{1}\eta_{1}s_{1}\\\
\beta_{2}l_{2}\eta_{2}s_{2}\end{subarray}}\sum_{\mathbf{R}\alpha_{1}\alpha_{2}}\sum_{a_{1}a_{2}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|<\Lambda_{c}\\\
\mathbf{Q}_{1}\in\mathcal{Q}_{l_{1}\eta_{1}}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|<\Lambda_{c}\\\
\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}\end{subarray}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\
V(\mathbf{q})\widetilde{v}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},\alpha_{1}}(\mathbf{k}_{1}-\mathbf{q})\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})\widetilde{u}^{(\eta_{2})*}_{\mathbf{Q}_{2}\beta_{2},a_{2}}(\mathbf{k}_{2})\widetilde{v}^{(\eta_{2})}_{\mathbf{Q}_{1}\beta_{1},\alpha_{2}}(\mathbf{k}_{2}-\mathbf{q})$
$\displaystyle\times$ $\displaystyle
e^{i(\eta_{1}\Delta\mathbf{K}_{l_{1}}+\mathbf{k}_{1}-\mathbf{Q}_{1}-\eta_{2}\Delta\mathbf{K}_{l_{2}}-\mathbf{k}_{2}+\mathbf{Q}_{2})\cdot\mathbf{R}}\left(f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}+c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\right)$
(S185)

Since
$\eta_{1}\Delta\mathbf{K}_{l_{1}}-\mathbf{Q}_{1}-\eta_{2}\Delta\mathbf{K}_{l_{2}}+\mathbf{Q}_{2}$
is a reciprocal lattice vector of the Moire lattice, there is
$e^{i(\eta_{1}\Delta\mathbf{K}_{l_{1}}-\mathbf{Q}_{1}-\eta_{2}\Delta\mathbf{K}_{l_{2}}+\mathbf{Q}_{2})\cdot\mathbf{R}}=1$
and hence

$\displaystyle{\hat{H}}_{J}=$
$\displaystyle\frac{1}{2N}\sum_{\eta_{1}s_{1}\eta_{2}s_{2}}\sum_{\mathbf{R}\alpha_{1}\alpha_{2}}\sum_{a_{1}a_{2}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\
V(\mathbf{q})\langle\widetilde{v}^{(\eta_{1})}_{\alpha_{1}}(\mathbf{k}_{1}-\mathbf{q})|\widetilde{u}^{(\eta_{1})}_{a_{1}}(\mathbf{k}_{1})\rangle\langle\widetilde{u}^{(\eta_{2})}_{a_{2}}(\mathbf{k}_{2})|\widetilde{v}^{(\eta_{2})}_{\alpha_{2}}(\mathbf{k}_{2}-\mathbf{q})\rangle$
$\displaystyle\times
e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}\left(f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}+c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\right)$
(S186)

We define the matrix

$X_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}(\mathbf{k}_{1},\mathbf{k}_{2})=\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\
V(\mathbf{q})\langle\widetilde{v}^{(\eta_{1})}_{\alpha_{1}}(\mathbf{k}_{1}-\mathbf{q})|\widetilde{u}^{(\eta_{1})}_{a_{1}}(\mathbf{k}_{1})\rangle\langle\widetilde{u}^{(\eta_{2})}_{a_{2}}(\mathbf{k}_{2})|\widetilde{v}^{(\eta_{2})}_{\alpha_{2}}(\mathbf{k}_{2}-\mathbf{q})\rangle$
(S187)

such that the interaction can be written as

$\displaystyle{\hat{H}}_{J}=$
$\displaystyle\frac{1}{2N}\sum_{\begin{subarray}{c}\eta_{1}\alpha_{1}a_{1}\\\
{\eta_{2}\alpha_{2}a_{2}}\end{subarray}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}X_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}(\mathbf{k}_{1},\mathbf{k}_{2})e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}$
$\displaystyle\times\left(f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}+c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\right)\
.$ (S188)

Since the low energy physics only involves conduction band around
$\mathbf{k}_{1,2}\approx 0$, we can approximate the interaction matrix
elements by those at $\mathbf{k}_{1}=\mathbf{k}_{2}=0$, i.e.,

${\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}=X_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}(0,0)=\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\
V(\mathbf{q})\langle\widetilde{v}^{(\eta_{1})}_{\alpha_{1}}(-\mathbf{q})|\widetilde{u}^{(\eta_{1})}_{a_{1}}(0)\rangle\langle\widetilde{u}^{(\eta_{2})}_{a_{2}}(0)|\widetilde{v}^{(\eta_{2})}_{\alpha_{2}}(-\mathbf{q})\rangle$
(S189)

We hence approximate ${\hat{H}}_{J}$ as

$\displaystyle{\hat{H}}_{J}=$
$\displaystyle\frac{1}{2N}\sum_{\begin{subarray}{c}\eta_{1}\alpha_{1}a_{1}\\\
{\eta_{2}\alpha_{2}a_{2}}\end{subarray}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}{\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}$
$\displaystyle\times\left(f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}+c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\right)\
.$ (S190)

After rearranging the creation and annihilation operators, ${\hat{H}}_{J}$ can
be rewritten as (up to a constant)

$\displaystyle{\hat{H}}_{J}=$
$\displaystyle-\frac{1}{N}\sum_{\begin{subarray}{c}\eta_{1}\alpha_{1}a_{1}\\\
{\eta_{2}\alpha_{2}a_{2}}\end{subarray}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}{\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}:f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}::c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}:\
,$ (S191)

where the normal ordered operators are given in Eq. S142.

We now study the ${\mathcal{J}}$ matrix (Eq. S189). One may explicitly apply
all the crystalline symmetries to ${\mathcal{J}}$ as we did around Eq. S166 to
obtain the constraints satisfied by ${\mathcal{J}}$. Here we adopt a
simplified notation, we use $|\widetilde{u}^{(\eta)}_{a}(\mathbf{k})\rangle$
to represent the vector
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\alpha,a}(\mathbf{k})$ and
$D(g)|\widetilde{u}^{(\eta)}_{a}(\mathbf{k})\rangle$ to represent the vector
obtained by acting $D(g)$ (Eqs. S29, S17 and S16) on the vector
$|\widetilde{u}^{(\eta)}_{a}(\mathbf{k})\rangle$. According to the $C_{2z}P$
representations formed by $f$\- and $c$-electrons (Eqs. S36 and S65):

$D^{(f)}_{\alpha\eta,\alpha^{\prime}\eta^{\prime}}(C_{2z}P)=-i[\sigma_{y}]_{\alpha\alpha^{\prime}}[\tau_{y}]_{\eta\eta^{\prime}}$
(S192)
$D^{(c)}_{a\eta,a^{\prime}\eta^{\prime}}(C_{2z}P)=i[\sigma_{y}\oplus\sigma_{y}]_{aa^{\prime}}[\tau_{y}]_{\eta\eta^{\prime}}$
(S193)

there are

$D(C_{2z}P)|\widetilde{v}^{(\eta)}_{\alpha}(\mathbf{k})\rangle=\sum_{\beta\eta^{\prime}}|\widetilde{v}^{(\eta^{\prime})}_{\beta}(\mathbf{k})\rangle
D_{\beta\eta^{\prime},\alpha\eta}^{(f)}(C_{2z}P)=-e^{i\pi\alpha}e^{i\frac{\pi}{2}\eta}|\widetilde{v}^{(-\eta)}_{\overline{\alpha}}(\mathbf{k})\rangle\
,$ (S194)
$D(C_{2z}P)|\widetilde{u}^{(\eta)}_{a}(\mathbf{k})\rangle=\sum_{b\eta^{\prime}}|\widetilde{u}^{(\eta)}_{a}(\mathbf{k})\rangle
D^{(c)}_{b\eta^{\prime},a\eta}(C_{2z}P)=e^{i\pi
a}e^{i\frac{\pi}{2}\eta}|\widetilde{u}^{(-\eta)}_{\overline{a}}(\mathbf{k})\rangle,$
(S195)

where $\overline{\alpha}=2,1$ for $\alpha=1,2$ and $\overline{a}=2,1,4,3$ for
$a=1,2,3,4$. Inserting $D^{\dagger}(C_{2z}P)D(C_{2z}P)=1$ into the inner
products in Eq. S189, it follows that

${\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}=-e^{i\pi(\alpha_{2}-a_{2})}{\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},-\eta_{2}\overline{\alpha}_{2}\overline{a}_{2}}=-e^{i\pi(a_{1}-\alpha_{1})}{\mathcal{J}}_{-\eta_{1}\overline{\alpha}_{1}\overline{a}_{1},\eta_{2}\alpha_{2}a_{2}}$
(S196)

Similarly, due to the $C_{2x}$ symmetry (Eqs. S36 and S2.2), the
${\mathcal{J}}$ matrix satisfies

${\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}={\mathcal{J}}_{\eta_{1}\overline{\alpha}_{1}\overline{a}_{1},\eta_{2}\overline{\alpha}_{2}\overline{a}_{2}}\
,$ (S197)

due to the $C_{2y}$ ($=C_{2x}C_{2z}$) symmetry (Eqs. S36 and S2.2), the
${\mathcal{J}}$ matrix satisfies

${\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}={\mathcal{J}}_{-\eta_{1}\alpha_{1}a_{1},-\eta_{2}\alpha_{2}a_{2}}\
,$ (S198)

and due to the $C_{3z}$ symmetry (Eqs. S36 and S2.2), the ${\mathcal{J}}$
matrix satisfies

${\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}=\zeta_{\eta_{1}\alpha_{1}}^{*}\zeta_{\eta_{1}a_{1}}\zeta_{\eta_{2}a_{2}}^{*}\zeta_{\eta_{2}\alpha_{2}}{\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}\
.$ (S199)

Here $\zeta_{\eta\alpha}=D^{(f)}_{\alpha\eta,\alpha\eta}(C_{3z})$ and
$\zeta_{\eta a}=D^{(c)}_{a\eta,a\eta}(C_{3z})$ are the $C_{3z}$ eigenvalues.

####  S3.3.2 Leading order coupling

For $w_{0}/w_{1}=0.8$, we find that the largest matrix elements of
${\mathcal{J}}$ are $16.38$meV and the second largest matrix elements are
$4.28$meV. We will only keep the largest matrix elements for simplicity. Due
to the constraints derived above, we find the largest elements of
${\mathcal{J}}$ are given by

${\mathcal{J}}_{\eta 13,\eta 13}={\mathcal{J}}_{\eta 24,\eta
24}=-{\mathcal{J}}_{\eta 13,-\eta 24}=-{\mathcal{J}}_{\eta 24,-\eta 13}=J\ .$
(S200)

where $J$ is the value of the largest matrix element. The values of $J$ for
other $w_{0}/w_{1}$ are tabulated in Table S6. The interaction Hamiltonian can
be simplified as

$\displaystyle{\hat{H}}_{J}=$
$\displaystyle-\frac{J}{N}\sum_{\mathbf{R}\alpha\eta
s_{1}s_{2}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}\big{(}:f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta
s_{2}}::c_{\mathbf{k}_{2},\alpha+2,\eta
s_{2}}^{\dagger}c_{\mathbf{k}_{1},\alpha+2,\eta
s_{1}}:-f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}f_{\mathbf{R},\overline{\alpha},-\eta
s_{2}}c_{\mathbf{k}_{2},\overline{\alpha}+2,-\eta
s_{2}}^{\dagger}c_{\mathbf{k}_{1},\alpha+2,\eta s_{1}}\big{)}\ .$ (S201)

We can write ${\hat{H}}_{J}$ in a more compact form as

$H_{J}=-\frac{J}{2N}\sum_{\mathbf{R}s_{1}s_{2}}\sum_{\alpha\alpha^{\prime}\eta\eta^{\prime}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}(\eta\eta^{\prime}+(-1)^{\alpha+\alpha^{\prime}}):f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s_{2}}::c_{\mathbf{k}_{2},\alpha^{\prime}+2,\eta^{\prime}s_{2}}^{\dagger}c_{\mathbf{k}_{1},\alpha+2,\eta
s_{1}}:$ (S202)

In Section S2.7 we have block-diagonalized the flat-U(4) generators (Eqs.
S112, S113 and S114) such that each block carries a fundamental U(4)
representation. We find that ${\hat{H}}_{J}$ is nothing but the coupling
between two fundamental flat-U(4) “spins”. We first define the fundamental
flat-U(4) generators at each site (momentum) for $f$-electrons ($c$-electrons)
as

$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R})=$
$\displaystyle\frac{1}{2}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\sum_{\alpha\alpha^{\prime}}\delta_{\xi,(-1)^{\alpha-1}\eta}\delta_{\xi,(-1)^{\alpha^{\prime}-1}\eta^{\prime}}\Sigma^{(\mu\nu,f)}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}\ ,$ (S203)
$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}(\mathbf{q})=$
$\displaystyle\frac{1}{2N}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}|<\Lambda_{c}\\\
|\mathbf{k}+\mathbf{q}|<\Lambda_{c}\end{subarray}}\sum_{aa^{\prime}=1,2}\delta_{\xi,(-1)^{a-1}\eta}\delta_{\xi,(-1)^{a^{\prime}-1}\eta^{\prime}}\Sigma^{(\mu\nu,c)}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}+\mathbf{q}\alpha\eta
s}^{\dagger}c_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s^{\prime}}\ ,$ (S204)
$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q})=$
$\displaystyle\frac{1}{2N}\sum_{\eta\eta^{\prime}}\sum_{ss^{\prime}}\sum_{\begin{subarray}{c}|\mathbf{k}|<\Lambda_{c}\\\
|\mathbf{k}+\mathbf{q}|<\Lambda_{c}\end{subarray}}\sum_{aa^{\prime}=3,4}\delta_{\xi,(-1)^{a-1}\eta}\delta_{\xi,(-1)^{a^{\prime}-1}\eta^{\prime}}\Sigma^{(\mu\nu,c)}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}+\mathbf{q}\alpha\eta
s}^{\dagger}c_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s^{\prime}}\ .$ (S205)

The $\Sigma^{(\mu\nu,f)}$ and $\Sigma^{(\mu\nu,c)}$ matrices are given in Eqs.
S106, S107, S108 and S109. Eqs. S112, S113 and S114 are related to the above
equations as
${\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}=\sum_{\mathbf{R}}{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R})$,
${\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}=N{\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}(\mathbf{q}=0)$,
${\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}=N{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q}=0)$.
Now we consider the product of flat-U(4) generators
$\sum_{\mu\nu\xi}e^{-i\mathbf{q}\cdot\mathbf{R}}{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R}){\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q})$.
Substituting the explicit forms of the $\Sigma^{(\mu\nu,f)}$ and
$\Sigma^{(\mu\nu,c)}$ matrices into the product and expanding the summation
over $\mu$, we obtain

$\displaystyle\sum_{\mu\nu\xi}e^{-i\mathbf{q}\cdot\mathbf{R}}{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R}){\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q})$
$\displaystyle=$
$\displaystyle\frac{1}{4N}\sum_{\xi}\sum_{\nu=0,x,y,z}\sum_{\eta\alpha}^{\prime}\sum_{\eta^{\prime},a=3,4}^{\prime}\sum_{s_{1}s_{2}s_{3}s_{4}}\sum_{\mathbf{k}}^{\prime}e^{-i\mathbf{q}\cdot\mathbf{R}}\bigg{(}f_{\mathbf{R},\alpha,\eta,s_{1}}^{\dagger}[{\varsigma}_{\nu}]_{s_{1}s_{2}}f_{\mathbf{R},\alpha,\eta
s_{2}}c_{\mathbf{k}+\mathbf{q},a,\eta^{\prime},s_{3}}^{\dagger}[{\varsigma}_{\nu}]_{s_{3}s_{4}}c_{\mathbf{k},a,\eta^{\prime},s_{4}}$
$\displaystyle+\left((-1)^{\alpha+a}-(-1)^{\alpha+a}\eta\eta^{\prime}\right)f_{\mathbf{R},\overline{\alpha},-\eta
s_{1}}^{\dagger}[{\varsigma}_{\nu}]_{s_{1}s_{2}}f_{\mathbf{R},\alpha,\eta
s_{2}}c_{\mathbf{k}+\mathbf{q},\overline{a},-\eta^{\prime},s_{3}}^{\dagger}[{\varsigma}_{\nu}]_{s_{3}s_{4}}c_{\mathbf{k},a,\eta^{\prime},s_{4}}$
$\displaystyle+\eta\eta^{\prime}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}[{\varsigma}_{\nu}]_{s_{1}s_{2}}f_{\mathbf{R}\alpha\eta
s_{2}}c_{\mathbf{k}+\mathbf{q}a\eta^{\prime}s_{3}}^{\dagger}[{\varsigma}_{\nu}]_{s_{3}s_{4}}c_{\mathbf{k}a\eta^{\prime}s_{4}}\bigg{)}\
,$ (S206)

where $\sum_{\eta\alpha}^{\prime}$ sums over $\eta,\alpha$ satisfying
$\xi=(-1)^{\alpha-1}\eta$, $\sum_{\eta,a=3,4}^{\prime}$ sums over $\eta,a$
satisfying $\xi=(-1)^{a-1}\eta$, and $\sum_{\mathbf{k}}^{\prime}$ sums over
$\mathbf{k}$ satisfying $|\mathbf{k}|<\Lambda_{c}$,
$|\mathbf{k}+\mathbf{q}|<\Lambda_{c}$. Due to the identity
$\sum_{\nu}[{\varsigma}_{\nu}]_{s_{1}s_{2}}[{\varsigma}_{\nu}]_{s_{3}s_{4}}=2\delta_{s_{1}s_{4}}\delta_{s_{2}s_{3}}$,
there is

$\displaystyle\sum_{\mu\nu\xi}e^{-i\mathbf{q}\cdot\mathbf{R}}{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R}){\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q})$
$\displaystyle=$
$\displaystyle\frac{1}{2N}\sum_{\xi}\sum_{\eta\alpha}^{\prime}\sum_{\eta^{\prime},a=3,4}^{\prime}\sum_{s_{1}s_{2}}\sum_{\mathbf{k}}^{\prime}e^{-i\mathbf{q}\cdot\mathbf{R}}\bigg{(}f_{\mathbf{R},\alpha,\eta,s_{1}}^{\dagger}f_{\mathbf{R},\alpha,\eta
s_{2}}c_{\mathbf{k}+\mathbf{q},a,\eta^{\prime},s_{2}}^{\dagger}c_{\mathbf{k},a,\eta^{\prime},s_{1}}+\eta\eta^{\prime}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta
s_{2}}c_{\mathbf{k}+\mathbf{q}a\eta^{\prime}s_{2}}^{\dagger}c_{\mathbf{k}a\eta^{\prime}s_{1}}$
$\displaystyle+\left((-1)^{\alpha+a}-(-1)^{\alpha+a}\eta\eta^{\prime}\right)f_{\mathbf{R},\overline{\alpha},-\eta
s_{1}}^{\dagger}f_{\mathbf{R},\alpha,\eta
s_{2}}c_{\mathbf{k}+\mathbf{q},\overline{a},-\eta^{\prime},s_{2}}^{\dagger}c_{\mathbf{k},a,\eta^{\prime},s_{1}}\bigg{)}\
.$ (S207)

The first line after the equal sign vanishes unless $\eta=\eta^{\prime}$,
which, according to $\xi=(-1)^{\alpha-1}\eta=(-1)^{a-1}\eta$, also implies
$a=\alpha+2$. The second line after the equal sign vanishes unless
$\eta=-\eta^{\prime}$, which, according to
$\xi=(-1)^{\alpha-1}\eta=(-1)^{a-1}\eta$, also implies
$a=\overline{\alpha}+2$. Therefore, the above equation can be simplified to

$\displaystyle\sum_{\mu\nu\xi}e^{-i\mathbf{q}\cdot\mathbf{R}}{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R}){\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q})$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\xi}\sum_{\eta\alpha}^{\prime}\sum_{s_{1}s_{2}}\sum_{\mathbf{k}}^{\prime}e^{-i\mathbf{q}\cdot\mathbf{R}}\bigg{(}f_{\mathbf{R},\alpha,\eta,s_{1}}^{\dagger}f_{\mathbf{R},\alpha,\eta
s_{2}}c_{\mathbf{k}+\mathbf{q},\alpha+2,\eta,s_{2}}^{\dagger}c_{\mathbf{k},\alpha+2,\eta,s_{1}}$
$\displaystyle-f_{\mathbf{R},\overline{\alpha},-\eta
s_{1}}^{\dagger}f_{\mathbf{R},\alpha,\eta
s_{2}}c_{\mathbf{k}+\mathbf{q},\alpha+2,\eta,s_{2}}^{\dagger}c_{\mathbf{k},\overline{\alpha}+2,-\eta,s_{1}}\bigg{)}\
.$ (S208)

The summation $\sum_{\xi}\sum_{\eta\alpha}^{\prime}$ can be equivalently
replaced by $\sum_{\eta\alpha}$ that has no restriction on $\eta\alpha$.
Comparing it with ${\hat{H}}_{J}$, we find ${\hat{H}}_{J}$ equals to

${\hat{H}}_{J}=-J\sum_{\mathbf{R}\mathbf{q}}\sum_{\mu\nu\xi}e^{-i\mathbf{q}\cdot\mathbf{R}}:{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R})::{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q}):\
.$ (S209)

####  S3.3.3 Other coupling terms

For $w_{0}/w_{1}=0.8$, we find that the second largest and third largest
matrix elements are given by

${\mathcal{J}}_{\eta 14,\eta 14}={\mathcal{J}}_{\eta 23,\eta
23}={\mathcal{J}}_{\eta 14,-\eta 23}={\mathcal{J}}_{\eta 14,-\eta
13}=J^{\prime}=\mathrm{4.28meV}\ ,$ (S210)

and

$\displaystyle{\mathcal{J}}_{\eta 12,\eta 23}={\mathcal{J}}_{\eta 14,\eta
21}={\mathcal{J}}_{\eta 21,\eta 14}={\mathcal{J}}_{\eta 23,\eta 12}$
$\displaystyle=$ $\displaystyle{\mathcal{J}}_{\eta 12,-\eta
14}={\mathcal{J}}_{\eta 14,-\eta 12}={\mathcal{J}}_{\eta 21,-\eta
23}={\mathcal{J}}_{\eta 23,-\eta 21}=-J^{\prime\prime}=-\mathrm{1.17meV}\ ,$
(S211)

respectively. Following the calculation around Section S3.3.2 in last
subsubsection, we find that the $J^{\prime}$ term can be rewritten as a ferro-
magnetic coupling btween $f$-electrons and $c$-electrons

$-J^{\prime}\sum_{\mathbf{R}\mathbf{q}}\sum_{\mu\nu\xi}e^{-i\mathbf{q}\cdot\mathbf{R}}:{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R})::{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,-\xi)}(\mathbf{q}):\
.$ (S212)

Different with $J$, which couples U(4) moments of $f$\- and $c$-electrons with
the same $\xi$ index, $J^{\prime}$ couples U(4) moments $f$\- and
$c$-electrons with opposite $\xi$ indices. The $J^{\prime}$ term respects the
flat-U(4) symmetry because it is a product of flat-U(4) moments, but it breaks
the U(4)$\times$U(4) symmetry.

The $J^{\prime\prime}$ term can be rewritten as

$J^{\prime\prime}\frac{1}{N}\sum_{\mathbf{R}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}\sum_{ss^{\prime}\eta}e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}\left(f_{\mathbf{R}\eta
1s}^{\dagger}c_{\mathbf{k}_{1}\eta 2s}+f_{\mathbf{R}-\eta
2s}^{\dagger}c_{\mathbf{k}_{1}-\eta 1s}\right)\left(c_{\mathbf{k}_{2}\eta
3s^{\prime}}^{\dagger}f_{\mathbf{R}\eta 2s^{\prime}}+c_{\mathbf{k}_{2}-\eta
4s^{\prime}}^{\dagger}f_{\mathbf{R}-\eta 1s^{\prime}}\right)+h.c.$ (S213)

$J^{\prime\prime}$ cannot be written as a product of flat-U(4) moments and it
breaks the flat-U(4) symmetry. To see the symmetry breaking, we consider to
apply the valley rotation $e^{-i\frac{\pi}{2}\hat{\Sigma}_{x0}}$ (Eq. S105).
Using the explicit U(4) representations defined in Eqs. S106, S107, S108 and
S109, we find that the $J^{\prime\prime}$ term transforms to

$\to-J^{\prime\prime}\frac{1}{N}\sum_{\mathbf{R}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{2}|<\Lambda_{c}}\sum_{ss^{\prime}\eta}e^{i(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{R}}\left(f_{\mathbf{R}\eta
1s}^{\dagger}c_{\mathbf{k}_{1}\eta 2s}+f_{\mathbf{R}-\eta
2s}^{\dagger}c_{\mathbf{k}_{1}-\eta 1s}\right)\left(c_{\mathbf{k}_{2}\eta
3s^{\prime}}^{\dagger}f_{\mathbf{R}\eta 2s^{\prime}}+c_{\mathbf{k}_{2}-\eta
4s^{\prime}}^{\dagger}f_{\mathbf{R}-\eta 1s^{\prime}}\right)+h.c.$ (S214)

Thus the $J^{\prime\prime}$ term is odd under
$e^{-i\frac{\pi}{2}\hat{\Sigma}_{x0}}$. Since $J^{\prime\prime}=-1.17$meV is
small compared to other energy scales, e.g., $U_{1}=57.95$meV, $J=16.38$meV,
flat-U(4) is still an approximate symmetry of the full interaction
Hamiltonian.

It was shown in Ref. Bernevig  _et al._ (2021a) that in the flat nonchiral
limit the flat-U(4) is an exact symmetry of the projected Coulomb Hamiltonian.
But here we find that the flat-U(4) symmetry is broken by $J^{\prime\prime}$.
The reason for this seemingly contradiction is that the $\Gamma_{3}$ states,
which are absent in the projected Coulomb Hamiltonian, are added in the heavy
fermion model. Notice that the inter-valley rotation in Ref. Bernevig  _et
al._ (2021a) is generated by the unitary operator $C_{2z}P$, while the inter-
valley rotation here is generated by $C_{2z}PS$, where $S$ is the third-chiral
operator defined in Section S2.7. $S$ is needed in our model because otherwise
the inter-valley rotation would not commute with the kinetic energy of the
dispersive $c$-bands. The representations of the two operators are given by
$D^{f}(C_{2z}P)=-i\sigma_{y}\tau_{y}$,
$D^{c}(C_{2z}P)=(i\sigma_{y}\tau_{y})\oplus(i\sigma_{y}\tau_{y})$ and
$D^{f}(C_{2z}PS)=-i\sigma_{y}\tau_{y}$,
$D^{c}(C_{2z}PS)=(-i\sigma_{y}\tau_{y})\oplus(i\sigma_{y}\tau_{y})$,
respectively, where the first and second blocks of $D^{c}$ act on $\Gamma_{3}$
and $\Gamma_{1}\oplus\Gamma_{2}$ states, respectively. Since $J$ and
$J^{\prime}$ terms only involve $f$-electrons and $\Gamma_{1}\oplus\Gamma_{2}$
$c$-electrons, on which $C_{2z}P$ and $C_{2z}PS$ actions are the same, the
conclusion of Ref. Bernevig  _et al._ (2021a) applies and hence $J$,
$J^{\prime}$ terms respect the flat-U(4) symmetry. However the
$J^{\prime\prime}$ term involves $\Gamma_{3}$ $c$-electrons, on which