$C_{2z}P$ and $C_{2z}PS$ actions are different, hence the conclusion of Ref.
Bernevig  _et al._ (2021a) does not apply to $J^{\prime\prime}$.

###  S3.4 Double hybridization terms (the third and fourth terms)

Here we determine the double hybridization terms, i.e., the third and fourth
terms discussed in Section S3.1. We first consider the
$:\rho_{fc}:V:\rho_{fc}:$ term

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

As explained at the beginning of Section S3.3, since the Wannier functions are
sharp Gaussian functions and $V(\mathbf{r}_{1}-\mathbf{r}_{2})$ decays quickly
as $|\mathbf{r}_{1}-\mathbf{r}_{2}|$ grows, the largest contribution to the
above equation should come from the $\mathbf{R}_{1}=\mathbf{R}_{2}$ terms.
Therefore, in the following we only focus on the
$\mathbf{R}_{1}=\mathbf{R}_{2}$ term. Substituting the Fourier transformation
of the interaction, we can rewrite the interaction as

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

According to Eq. S45, there is

$\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})\
.$ (S217)

Therefore, we have

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

Since
$\eta_{1}\Delta\mathbf{K}_{l_{1}}-\mathbf{Q}_{1}+\eta_{2}\Delta\mathbf{K}_{l_{2}}-\mathbf{Q}_{2}$
is always a moiré reciprocal lattice, there must be
$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$.
We can further simplify the Hamiltonian as

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

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{v}^{(\eta_{2})}_{\alpha_{2}}(\mathbf{k}_{2}+\mathbf{q})|\widetilde{u}^{(\eta_{2})}_{a_{2}}(\mathbf{k}_{2})\rangle$
(S220)

such that the interaction can be written as

$\displaystyle{\hat{H}}_{\widetilde{J}+}=$
$\displaystyle\frac{1}{2N}\sum_{\eta_{1}s_{1}\eta_{2}s_{2}}\sum_{\begin{subarray}{c}\mathbf{R}\alpha_{1}\alpha_{2}\\\
a_{1}a_{2}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|<\Lambda_{c}\\\
|\mathbf{k}_{2}|<\Lambda_{c}\end{subarray}}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}}f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\
.$ (S221)

For low energy states around $\mathbf{k}=0$, we can approximate
$X_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}(\mathbf{k}_{1},\mathbf{k}_{2})$
by $X_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}(0,0)$ (dubbed as
$\widetilde{J}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}$)

$\widetilde{J}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}=\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{v}^{(\eta_{2})}_{\alpha_{2}}(\mathbf{q})|\widetilde{u}^{(\eta_{2})}_{a_{2}}(0)\rangle$
(S222)

We find that, according to the time-reversal symmetry, the matrix elements of
$\widetilde{J}$ are same as those of ${\mathcal{J}}$ (Eq. S189). The $T$
symmetry (Eq. S16) and its representations (Eqs. S36 and S2.2) imply
$\widetilde{v}^{(\eta_{2})}_{\mathbf{Q}\beta,\alpha_{2}}(\mathbf{k})=\widetilde{v}^{(-\eta_{2})*}_{-\mathbf{Q}\beta,\alpha_{2}}(-\mathbf{k})$,
$\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}\beta,a_{2}}(0)=\widetilde{u}^{(-\eta_{2})*}_{-\mathbf{Q}\beta,a_{2}}(0)$
and hence

$\langle\widetilde{v}^{(\eta_{2})}_{\alpha_{2}}(\mathbf{q})|\widetilde{u}^{(\eta_{2})}_{a_{2}}(0)\rangle=\langle\widetilde{u}^{(-\eta_{2})}_{{a}_{2}}(0)|\widetilde{v}^{(-\eta_{2})}_{{\alpha}_{2}}(-\mathbf{q})\rangle\
.$ (S223)

Thus we have

$\displaystyle\widetilde{J}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\alpha_{2}a_{2}}=\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={\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},-\eta_{2}{\alpha}_{2}{a}_{2}}\
.$ (S224)

Thus ${\hat{H}}_{\widetilde{J}+}$ can be written as

$\displaystyle{\hat{H}}_{\widetilde{J}+}=$
$\displaystyle\frac{1}{2N}\sum_{\eta_{1}s_{1}\eta_{2}s_{2}}\sum_{\begin{subarray}{c}\mathbf{R}\alpha_{1}\alpha_{2}\\\
a_{1}a_{2}\end{subarray}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|<\Lambda_{c}\\\
|\mathbf{k}_{2}|<\Lambda_{c}\end{subarray}}{\mathcal{J}}_{\eta_{1}\alpha_{1}a_{1},\eta_{2}\overline{\alpha}_{2}\overline{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}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}\
.$ (S225)

According to Eq. S200, the non-negligible matrix elements of ${\mathcal{J}}$
and $\widetilde{{\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\ ,$
(S226) $\widetilde{{\mathcal{J}}}_{\eta 13,-\eta
13}=\widetilde{{\mathcal{J}}}_{\eta 24,-\eta
24}=-\widetilde{{\mathcal{J}}}_{\eta 13,\eta
24}=-\widetilde{{\mathcal{J}}}_{\eta 24,\eta 13}=J\ .$ (S227)

The interaction Hamiltonian ${\hat{H}}_{\widetilde{J}+}$ can be explicitly
written as

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

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

${\hat{H}}_{\widetilde{J}+}=-\frac{J}{4N}\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}}^{\dagger}c_{\mathbf{k}_{2},\alpha^{\prime}+2,\eta^{\prime}s_{2}}c_{\mathbf{k}_{1},\alpha+2,\eta
s_{1}}$ (S229)

The fourth term in Section S3.1, i.e., $:\rho_{cf}:V:\rho_{cf}:$, is just the
hermitian conjugation of ${\hat{H}}_{\widetilde{J}+}$. Thus the total double
hybridization interaction is given by

${\hat{H}}_{\widetilde{J}}={\hat{H}}_{\widetilde{J}+}+{\hat{H}}_{\widetilde{J}+}^{\dagger}\
.$ (S230)

###  S3.5 Density-hybridization terms (the sixth, seventh, eighth, and ninth
terms)

####  S3.5.1 Vanishing sixth and seventh terms

We now study the sixth and seventh terms discussed in Section S3.1. The sixth
term ($:\hat{\rho}_{ff}:V:\hat{\rho}_{fc}:+(\leftrightarrow)$) has the form

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

We have omitted the overlap of local orbitals at different sites. Changing the
variables $\mathbf{r}_{1}\to\mathbf{r}_{1}+\mathbf{R}$,
$\mathbf{r}_{2}\to\mathbf{r}_{2}+\mathbf{R}$, we can rewrite the interaction
as

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

We have made use of the fact that $\eta\Delta\mathbf{K}_{l}-\mathbf{Q}$ is a
moiré reciprocal lattice and
$e^{i(\eta\Delta\mathbf{K}_{l}-\mathbf{Q})\cdot\mathbf{R}}=1$. We introduce
the matrix

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

such that the interaction can be written as

$\displaystyle{\hat{H}}_{67}=$
$\displaystyle\frac{1}{2\sqrt{N}}\sum_{\mathbf{R}\alpha\eta_{1}s_{1}}\sum_{\eta_{2}s_{2}\alpha^{\prime}a}\sum_{|\mathbf{k}|<\Lambda_{c}}X^{(\eta_{2})}_{\alpha^{\prime},a}(\mathbf{k})e^{i\mathbf{k}\cdot\mathbf{R}}\left(:f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta_{1}s_{1}}:f_{\mathbf{R}\alpha^{\prime}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}a\eta_{2}s_{2}}+f_{\mathbf{R}\alpha^{\prime}\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}a\eta_{2}s_{2}}:f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta_{1}s_{1}}:\right)\
.$ (S234)

Since only the conduction states around $\mathbf{k}=0$ are relevant in the low
energy physics, we will approximate
$X_{\alpha^{\prime},a}^{(\eta_{2})}(\mathbf{k})$ by
$X_{\alpha^{\prime},a}^{(\eta_{2})}(0)$

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

Now we show that $X_{\alpha^{\prime},a}^{(\eta_{2})}(0)$ is guaranteed to be
vanishing by symmetries. According to the $C_{3z}$ symmetry (Eqs. S16 and
S23), there is

$e^{i\frac{2\pi}{3}\eta_{2}[\sigma_{z}]_{\beta\beta}}w^{(\eta_{2})}_{l\beta,\alpha^{\prime}}(\mathbf{r}_{2})=w^{(\eta_{2})}_{l\beta,\alpha^{\prime}}(C_{3z}\mathbf{r}_{2})\zeta_{\eta_{2}\alpha}^{f},\qquad
e^{i\frac{2\pi}{3}\eta_{2}[\sigma_{z}]_{\beta\beta}}\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}\beta,a}(0)=\widetilde{u}^{(\eta_{2})}_{C_{3z}\mathbf{Q}\beta,a}(0)\zeta^{c}_{\eta_{2}a}$
(S236)

with $\zeta_{\eta\alpha}^{f}$ and $\zeta_{\eta\alpha}^{c}$ being $C_{3z}$
eigenvalues given in Eq. S36 and Section S2.2, respectively, i.e.,

$\zeta_{\eta_{2}\alpha}^{f}=e^{i\frac{2\pi}{3}\eta_{2}[\sigma_{z}]_{\alpha\alpha}},\qquad\zeta_{\eta_{2}\alpha}^{c}=\begin{cases}e^{i\frac{2\pi}{3}\eta_{2}[\sigma_{z}]_{aa}},\qquad&a=1,2\\\
1,\qquad&a=3,4\end{cases}\ .$ (S237)

Substituting this condition into Eq. S235, we have

$\displaystyle X_{\alpha^{\prime},a}^{(\eta_{2})}(0)=$
$\displaystyle\zeta_{\eta_{2},\alpha^{\prime}}^{f*}\zeta_{\eta_{2},a}^{c}\frac{1}{\sqrt{\Omega_{0}}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}V(\mathbf{r}_{1}-\mathbf{r}_{2})n_{f}(\mathbf{r}_{1})\sum_{l\beta}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta_{2}}}w^{(\eta_{2})*}_{l\beta,\alpha^{\prime}}(C_{3z}\mathbf{r}_{2})\widetilde{u}^{(\eta_{2})}_{C_{3z}\mathbf{Q}\beta,a}(0)e^{-i\mathbf{Q}\cdot\mathbf{r}_{2}}$
(S238)

Changing the variables as $\mathbf{r}_{1}\to C_{3z}^{-1}\mathbf{r}_{1}$,
$\mathbf{r}_{2}\to C_{3z}^{-1}\mathbf{r}_{2}$, $\mathbf{Q}\to
C_{3z}^{-1}\mathbf{Q}$, we have

$\displaystyle X_{\alpha^{\prime},a}^{(\eta_{2})}(0)=$
$\displaystyle\zeta_{\eta_{2},\alpha^{\prime}}^{f*}\zeta_{\eta_{2},a}^{c}\frac{1}{\sqrt{\Omega_{0}}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}V(C_{3z}^{-1}\mathbf{r}_{1}-C_{3z}^{-1}\mathbf{r}_{2})n_{f}(C_{3z}^{-1}\mathbf{r}_{1})\sum_{l\beta}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta_{2}}}w^{(\eta_{2})*}_{l\beta,\alpha^{\prime}}(\mathbf{r}_{2})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}\beta,a}(0)e^{-i\mathbf{Q}\cdot\mathbf{r}_{2}}$
(S239)

Making use of $V(C_{3z}\mathbf{r})=V(\mathbf{r})$,
$n_{f}(C_{3z}\mathbf{r})=n_{f}(\mathbf{r})$, we have

$X_{\alpha^{\prime},a}^{(\eta_{2})}(0)=\zeta_{\eta_{2},\alpha^{\prime}}^{f*}\zeta_{\eta_{2},a}^{c}X_{\alpha^{\prime},a}^{(\eta_{2})}(0)\
.$ (S240)

Therefore, $X_{\alpha^{\prime},a}^{(\eta_{2})}(0)$ is zero if
$\zeta_{\eta_{2},\alpha^{\prime}}^{f}\neq\zeta_{\eta_{2},a}^{c}$. Due to the
$C_{3z}$ eigenvalues the only possible nonzero matrix elements are
$X_{1,1}^{(\eta_{2})}(0)$ and $X_{2,2}^{(\eta_{2})}(0)$. The above symmetry
analysis also applies for other symmetries, for example the particle-hole
symmetry $P$. As explained in the end of Section S2.2, the $a=1,2$ states have
opposite $P$ eigenvalues with the $\alpha=1,2$ states. Thus
$X_{1,1}^{(\eta_{2})}(0)$ and $X_{2,2}^{(\eta_{2})}(0)$ are also guaranteed to
be zero. Thus the matrix elements in ${\hat{H}}_{67}$ vanish to zeroth order
of $\mathbf{k}$.

Since the seventh term discussed in Section S3.1 is the hermitian conjugation
of the sixth term, its matrix elements also vanish to zeroth order of
$\mathbf{k}$. The sixth and seventh interaction terms will be non-vanishing if
we include higher order terms of $\mathbf{k}$ in the matrix elements, i.e.,
the $\mathbf{k}$-dependence in
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\alpha,a}(\mathbf{k})$. But these matrix
elements will be small for low energy conduction band states around
$\mathbf{k}=0$. We claim that these $\mathbf{k}$-depend matrix elements in the
density-hybridization interaction are irrelevant in the low energy physics.

####  S3.5.2 The eighth and ninth terms

We now study the eighth and ninth terms discussed in Section S3.1. They are
weak because the hybridization, i.e., $\langle f^{\dagger}c\rangle$ and
$\langle c^{\dagger}f\rangle$, is small, as will be discussed in Section S4.
The eighth term ($:\hat{\rho}_{cc}:V:\hat{\rho}_{fc}:+(\leftrightarrow)$) has
the form

$\displaystyle{\hat{H}}_{K+}=$ $\displaystyle\frac{1}{2\Omega_{\rm
tot}^{\frac{3}{2}}}\int
d^{2}\mathbf{r}_{1}d^{2}\mathbf{r}_{2}V(\mathbf{r}_{1}-\mathbf{r}_{2})\sum_{\beta_{1}l_{1}\eta_{1}s_{1}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|,|\mathbf{k}_{1}^{\prime}|<\Lambda_{c}\\\
a_{1}a_{1}^{\prime}\end{subarray}}\sum_{\mathbf{Q}_{1},\mathbf{Q}_{1}^{\prime}\in\mathcal{Q}_{l_{1}\eta_{1}}}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}_{1}^{\prime}\beta_{1},a_{1}^{\prime}}(\mathbf{k}_{1}^{\prime})$
$\displaystyle\times\sum_{\beta_{2}l_{2}\eta_{2}s_{2}}\sum_{\mathbf{R}\alpha}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|<\Lambda_{c}\\\
a_{2}\end{subarray}}\sum_{\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}}w^{(\eta_{2})*}_{l_{2}\beta_{2},\alpha}(\mathbf{r}_{2}-\mathbf{R})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a}(\mathbf{k}_{2})e^{i(\mathbf{k}_{2}-\mathbf{Q}_{2})\cdot\mathbf{r}_{2}}e^{i\eta_{2}\Delta\mathbf{K}_{l_{2}}\cdot\mathbf{R}}$
$\displaystyle\times\left(:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}:f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}+f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}:\right)\
.$ (S241)

Applying the Fourier transformation of $V(\mathbf{r})$, we obtain

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

According to Eq. S45, there are

$\displaystyle\int
d^{2}\mathbf{r}_{2}w^{(\eta_{2})*}_{l_{2}\beta_{2},\alpha}(\mathbf{r}_{2}-\mathbf{R})e^{i(\mathbf{k}_{2}-\mathbf{Q}_{2}-\mathbf{k}_{1}+\mathbf{Q}_{1}+\mathbf{k}_{1}^{\prime}-\mathbf{Q}_{1}^{\prime})\cdot\mathbf{r}_{2}}=$
$\displaystyle=$ $\displaystyle\sqrt{\Omega_{0}}\
e^{i(\mathbf{k}_{2}-\mathbf{Q}_{2}-\mathbf{k}_{1}+\mathbf{Q}_{1}+\mathbf{k}_{1}^{\prime}-\mathbf{Q}_{1}^{\prime})\cdot\mathbf{R}}\widetilde{v}^{(\eta_{2})*}_{\mathbf{Q}_{2}-\mathbf{Q}_{1}+\mathbf{Q}_{1}^{\prime}\beta_{2},\alpha}(\mathbf{k}_{2}-\mathbf{k}_{1}+\mathbf{k}_{1}^{\prime})\
.$ (S243)

Since
$\eta_{2}\Delta\mathbf{K}_{l_{2}}-\mathbf{Q}_{2}+\mathbf{Q}_{1}-\mathbf{Q}_{1}^{\prime}$
is a moiré reciprocal lattice, it must be that
$e^{i(\eta_{2}\Delta\mathbf{K}_{l_{2}}-\mathbf{Q}_{2}+\mathbf{Q}_{1}-\mathbf{Q}_{1}^{\prime})\cdot\mathbf{R}}=1$.
We hence obtain

$\displaystyle{\hat{H}}_{K+}=$
$\displaystyle\frac{1}{2N^{\frac{3}{2}}\Omega_{0}}\sum_{\beta_{1}l_{1}\eta_{1}s_{1}}\sum_{\begin{subarray}{c}|\mathbf{k}_{1}|,|\mathbf{k}_{1}^{\prime}|<\Lambda_{c}\\\
a_{1}a_{1}^{\prime}\end{subarray}}\sum_{\mathbf{Q}_{1},\mathbf{Q}_{1}^{\prime}\in\mathcal{Q}_{l_{1}\eta_{1}}}V(\mathbf{k}_{1}-\mathbf{Q}_{1}-\mathbf{k}_{1}^{\prime}+\mathbf{Q}_{1}^{\prime})\widetilde{u}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}^{\prime}\beta_{1},a_{1}^{\prime}}(\mathbf{k}_{1}^{\prime})$
$\displaystyle\times\sum_{\beta_{2}l_{2}\eta_{2}s_{2}}\sum_{\mathbf{R}\alpha}\sum_{\begin{subarray}{c}|\mathbf{k}_{2}|<\Lambda_{c}\\\
a_{2}\end{subarray}}\sum_{\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}}\widetilde{v}^{(\eta_{2})*}_{\mathbf{Q}_{2}-\mathbf{Q}_{1}+\mathbf{Q}_{1}^{\prime},\beta_{2};\alpha}(\mathbf{k}_{2}-\mathbf{k}_{1}+\mathbf{k}_{1}^{\prime})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a}(\mathbf{k}_{2})e^{i(\mathbf{k}_{2}-\mathbf{k}_{1}+\mathbf{k}_{1}^{\prime})\cdot\mathbf{R}}$
$\displaystyle\times\left(:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}:f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}+f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}:\right)\
.$ (S244)

We define the matrix

$\displaystyle X_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha
a_{2}}(\mathbf{k}_{1},\mathbf{k}_{1}^{\prime},\mathbf{k}_{2})=$
$\displaystyle\sum_{l_{1}\beta_{1}l_{2}\beta_{2}}\sum_{\mathbf{Q}_{1},\mathbf{Q}_{1}^{\prime}\in\mathcal{Q}_{l_{1}\eta_{1}}}\sum_{\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}}V(\mathbf{k}_{1}-\mathbf{Q}_{1}-\mathbf{k}_{1}^{\prime}+\mathbf{Q}_{1}^{\prime})$
$\displaystyle\times\widetilde{u}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(\mathbf{k}_{1})\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}^{\prime}\beta_{1},a_{1}^{\prime}}(\mathbf{k}_{1}^{\prime})\widetilde{v}^{(\eta_{2})*}_{\mathbf{Q}_{2}-\mathbf{Q}_{1}+\mathbf{Q}_{1}^{\prime},\beta_{2};\alpha}(\mathbf{k}_{2}-\mathbf{k}_{1}+\mathbf{k}_{1}^{\prime})\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a}(\mathbf{k}_{2})$
(S245)

such that the interaction can be written as

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

Since the conduction bands only have low energy states around $\mathbf{k}=0$,
thus in the following we approximate
$X_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha
a_{2}}(\mathbf{k}_{1},\mathbf{k}_{1}^{\prime},\mathbf{k}_{2})$ by
$X_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha a_{2}}(0,0,0)$, dubbed as
${\mathcal{K}}_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha a_{2}}$

$\displaystyle{\mathcal{K}}_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha
a_{2}}=$
$\displaystyle\sum_{l_{1}\beta_{1}l_{2}\beta_{2}}\sum_{\mathbf{Q}_{1},\mathbf{Q}_{1}^{\prime}\in\mathcal{Q}_{l_{1}\eta_{1}}}\sum_{\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}}V(-\mathbf{Q}_{1}+\mathbf{Q}_{1}^{\prime})\widetilde{u}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(0)\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}^{\prime}\beta_{1},a_{1}^{\prime}}(0)\widetilde{v}^{(\eta_{2})*}_{\mathbf{Q}_{2}-\mathbf{Q}_{1}+\mathbf{Q}_{1}^{\prime},\beta_{2};\alpha}(0)\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a}(0)$
$\displaystyle=$
$\displaystyle\sum_{l_{1}\beta_{1}l_{2}\beta_{2}}\sum_{\mathbf{Q}_{1}\in\mathcal{Q}_{l_{1}\eta_{1}}}\sum_{\mathbf{Q}_{2}\in\mathcal{Q}_{l_{2}\eta_{2}}}\sum_{\mathbf{G}}V(\mathbf{G})\widetilde{u}^{(\eta_{1})*}_{\mathbf{Q}_{1}\beta_{1},a_{1}}(0)\widetilde{u}^{(\eta_{1})}_{\mathbf{Q}_{1}-\mathbf{G},\beta_{1};a_{1}^{\prime}}(0)\widetilde{v}^{(\eta_{2})*}_{\mathbf{Q}_{2}-\mathbf{G},\beta_{2};\alpha}(0)\widetilde{u}^{(\eta_{2})}_{\mathbf{Q}_{2}\beta_{2},a}(0)$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{G}}V(\mathbf{G})\langle\widetilde{u}^{(\eta_{1})}_{a_{1}}(0)|\widetilde{u}^{(\eta_{1})}_{a_{1}^{\prime}}(\mathbf{G})\rangle\langle\widetilde{v}^{(\eta_{2})}_{\alpha}(\mathbf{G})|\widetilde{u}^{(\eta_{2})}_{a_{2}}(0)\rangle\
.$ (S247)

The Hamiltonian ${\hat{H}}_{K+}$ can be written in a simple form

$\displaystyle{\hat{H}}_{K+}=$
$\displaystyle\frac{1}{2N^{\frac{3}{2}}\Omega_{0}}\sum_{\eta_{1}s_{1}\eta_{2}s_{2}}\sum_{a_{1}a_{1}^{\prime}a_{2}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{1}^{\prime}|,|\mathbf{k}_{2}|<\Lambda_{c}}\sum_{\mathbf{R}\alpha}e^{i(\mathbf{k}_{2}-\mathbf{k}_{1}+\mathbf{k}_{1}^{\prime})\cdot\mathbf{R}}{\mathcal{K}}_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha
a_{2}}$
$\displaystyle\times\left(:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}:f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}+f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2}a_{2}\eta_{2}s_{2}}:c_{\mathbf{k}_{1}a_{1}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime}a_{1}^{\prime}\eta_{1}s_{1}}:\right)\
.$ (S248)

We now study the ${\mathcal{K}}$ matrix. According to the $C_{2z}P$ symmetry
(Eqs. S36 and S65), there are

$C_{2z}P|\widetilde{v}^{(\eta)}_{\alpha}(\mathbf{k})\rangle=-e^{i\pi\alpha}e^{i\frac{\pi}{2}\eta}|\widetilde{v}^{(-\eta)}_{\overline{\alpha}}(\mathbf{k})\rangle,\qquad
C_{2z}P|\widetilde{u}^{(\eta)}_{a}(\mathbf{k})\rangle=e^{i\pi
a}e^{i\frac{\pi}{2}\eta}|\widetilde{v}^{(-\eta)}_{\overline{a}}(\mathbf{k})\rangle,$
(S249)

where $\overline{\alpha}=2,1$ for $\alpha=1,2$ and $\overline{a}=2,1,4,3$ for
$a=1,2,3,4$. Then if follows that the ${\mathcal{K}}$ matrix satisfies

${\mathcal{K}}_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha
a_{2}}=e^{i\pi(a_{1}^{\prime}-a_{1})}{\mathcal{K}}_{-\eta_{1}\overline{a}_{1}\overline{a}_{1}^{\prime},\eta_{2}\alpha
a_{2}}=-e^{i\pi(a_{2}-\alpha)}{\mathcal{K}}_{\eta_{1}a_{1}a_{1}^{\prime},-\eta_{2}\overline{\alpha}\overline{a}_{2}}$
(S250)

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

${\mathcal{K}}_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha
a_{2}}={\mathcal{K}}_{\eta_{1}\overline{a}_{1}\overline{a}_{1}^{\prime},\eta_{2}\overline{\alpha}\overline{a}_{2}}\
.$ (S251)

Due to the $C_{2z}TP$ symmetry (Eqs. S36 and S65), the ${\mathcal{K}}$ matrix
satisfies

${\mathcal{K}}_{\eta_{1}a_{1}a_{1}^{\prime},\eta_{2}\alpha
a_{2}}=(-1)^{a_{1}+a_{1}^{\prime}}{\mathcal{K}}_{\eta_{1}\overline{a}_{1}^{\prime}\overline{a}_{1},\eta_{2}\alpha
a_{2}}\ .$ (S252)

Numerically we find that there is only one independent non-negligible matrix
element

${\mathcal{K}}_{\eta_{1}14,\eta_{2}24}={\mathcal{K}}_{\eta_{1}23,\eta_{2}13}=-{\mathcal{K}}_{\eta_{1}32,\eta_{2}24}=-{\mathcal{K}}_{\eta_{1}41,\eta_{2}13}=\eta_{1}\eta_{2}K\
,$ (S253)

where $K=4.887$meV for $w_{0}/w_{1}=0.8$. The second largest matrix element is
$0.82$meV. One can verify that the ${\mathcal{K}}$ matrix satisfies the
constraints derived above. Therefore, we can write the Hamiltonian
${\hat{H}}_{K+}$ as

$\displaystyle{\hat{H}}_{K+}=$
$\displaystyle\frac{K}{N^{\frac{3}{2}}\Omega_{0}}\sum_{\eta_{1}s_{1}\eta_{2}s_{2}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{1}^{\prime}|,|\mathbf{k}_{2}|<\Lambda_{c}}\sum_{\mathbf{R}\alpha}e^{i(\mathbf{k}_{2}+\mathbf{k}_{1}^{\prime}-\mathbf{k}_{1})\cdot\mathbf{R}}\eta_{1}\eta_{2}\bigg{(}c_{\mathbf{k}_{1}\overline{\alpha}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime},\alpha+2,\eta_{1}s_{1}}f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2},\alpha+2,\eta_{2}s_{2}}$
$\displaystyle\qquad-
f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2},\alpha+2,\eta_{2}s_{2}}c_{\mathbf{k}_{1},\overline{\alpha}+2,\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime},\alpha,\eta_{1}s_{1}}\bigg{)}\
.$ (S254)

The ninth term discussed in Section S3.1, i.e.,
$:\hat{\rho}_{cc}:V:\hat{\rho}_{cf}:+(\leftrightarrow)$, is just the hermitian
conjugation of ${\hat{H}}_{K+}$. Thus the total density-hybridization
Hamiltonian is given by

${\hat{H}}_{K}={\hat{H}}_{K+}+{\hat{H}}_{K+}^{\dagger}\ .$ (S255)

###  S3.6 Summary of the interaction Hamiltonian and its symmetries

####  S3.6.1 Summary of the interaction Hamiltonian

In summary, the interaction Hamiltonian is given by

${\hat{H}}_{I}={\hat{H}}_{U}+{\hat{H}}_{V}+{\hat{H}}_{W}+{\hat{H}}_{J}+{\hat{H}}_{\widetilde{J}}+{\hat{H}}_{K}\
.$ (S256)

${\hat{H}}_{U}$ (Eq. S157) is the density-density interaction of the local
orbitals

${\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}}:\
,$ (S257)

where $\langle\mathbf{R}\mathbf{R}^{\prime}\rangle$ indexes nearest neighbor
pairs on the triangular lattice, $U_{1}$ is the on-site repulsion interaction,
and $U_{2}$ is the repulsion between nearest neighbors. ${\hat{H}}_{V}$ (Eq.
S179) is the density-density interaction of the conduction bands

${\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}}:\
,$ (S258)

where $V(\mathbf{q})$ (Eq. S145) is the Fourier transformation of the double-
gate-screened Coulomb interaction. ${\hat{H}}_{W}$ (Eq. S172) is the density-
density interaction between the local orbitals and the conduction bands

${\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}}:\
,$ (S259)

where $W_{1}=W_{2}$ and $W_{3}=W_{4}$ due to symmetries. ${\hat{H}}_{J}$ (Eq.
S202) is the exchange interaction between local orbitals and the conduction
bands

$\displaystyle
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}}:\ .$ (S260)

${\hat{H}}_{J}$ can be equivalently written as ferromagnetic coupling between
the flat-U(4) momenta (Eq. S209)

${\hat{H}}_{J}=-J\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})\
,$ (S261)

with the momenta given by Eqs. S203 and S205.
${\hat{H}}_{\widetilde{J}}={\hat{H}}_{\widetilde{J}+}+{\hat{H}}_{\widetilde{J}+}^{\dagger}$
is the double-hybridization interaction, with ${\hat{H}}_{\widetilde{J}+}$
(Eq. S229) given by

$\displaystyle{\hat{H}}_{\widetilde{J}+}=-\frac{J}{4N}\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}}^{\dagger}c_{\mathbf{k}_{2},\alpha^{\prime}+2,\eta^{\prime}s_{2}}c_{\mathbf{k}_{1},\alpha+2,\eta
s_{1}}\ .$ (S262)

Notice that the coefficients $J$ in ${\hat{H}}_{J}$ and
${\hat{H}}_{\widetilde{J}}$ are the same coefficient because the corresponding
matrix elements are related by $T$ symmetry, as explained in Section S3.4.
${\hat{H}}_{K}={\hat{H}}_{K+}+{\hat{H}}_{K+}^{\dagger}$ is the density-
hybridization interaction, with ${\hat{H}}_{K+}$ (Section S3.5.2) given by

$\displaystyle{\hat{H}}_{K+}=$
$\displaystyle\frac{K}{N^{\frac{3}{2}}\Omega_{0}}\sum_{\eta_{1}s_{1}\eta_{2}s_{2}}\sum_{|\mathbf{k}_{1}|,|\mathbf{k}_{1}^{\prime}|,|\mathbf{k}_{2}|<\Lambda_{c}}\sum_{\mathbf{R}\alpha}e^{i(\mathbf{k}_{2}+\mathbf{k}_{1}^{\prime}-\mathbf{k}_{1})\cdot\mathbf{R}}\eta_{1}\eta_{2}\bigg{(}c_{\mathbf{k}_{1}\overline{\alpha}\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime},\alpha+2,\eta_{1}s_{1}}f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2},\alpha+2,\eta_{2}s_{2}}$
$\displaystyle\qquad-
f_{\mathbf{R}\alpha\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}_{2},\alpha+2,\eta_{2}s_{2}}c_{\mathbf{k}_{1},\overline{\alpha}+2,\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}_{1}^{\prime},\alpha,\eta_{1}s_{1}}\bigg{)}\
.$ (S263)

The parameters $U_{1,2}$, $W_{1,3}$, $J$, $K$ for different $w_{0}/w_{1}$ are
given in Table S6.

####  S3.6.2 Symmetries of the interaction Hamiltonian

The continuous BM model (with interaction) has a U(2)$\times$U(2) symmetry
when valley is a good quantum number. Thus the effective topological heavy
fermion model must also have a U(2)$\times$U(2) symmetry. As discussed in
Section S2.7, when $M=0$, the presence of the $C_{2z}PS$ symmetry, where $S$
is the third chiral symmetry, will enhance the U(2)$\times$U(2) symmetry to a
U(4) symmetry, referred to as the flat U(4) symmetry. Recall that the Coulomb
interaction is given by (Eq. S137)

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

with $:\hat{\rho}(\mathbf{r}_{1}):$ being the normal ordered projected density
operator (Section S3.1). Since $C_{2z}P$ (Eqs. S36 and S65) is a symmetry of
the original model, $:\hat{\rho}(\mathbf{r}):$ and hence ${\hat{H}}_{I}$
commute with $C_{2z}P$. (In Ref. Bernevig  _et al._ (2021a), it is $C_{2z}P$
that gives the U(4) symmetry, provided that remote bands are projected out and
the kinetic energy of the flat bands is neglected.) Thus, in order for $H_{I}$
to have the flat-U(4) symmetry, we only need $H_{I}$ to commute with $S$
defined in Eq. S100. One can easily verify that ${\hat{H}}_{U}$,
${\hat{H}}_{V}$, ${\hat{H}}_{W}$, ${\hat{H}}_{J}$, ${\hat{H}}_{\widetilde{J}}$
commute with $S$ whereas ${\hat{H}}_{K}$ anti-commutes with $S$. Therefore,
${\hat{H}}_{U}$, ${\hat{H}}_{V}$, ${\hat{H}}_{W}$, ${\hat{H}}_{J}$,
${\hat{H}}_{\widetilde{J}}$ respect the flat-U(4) symmetry while
${\hat{H}}_{K}$ breaks the flat-U(4) symmetry. (The neglected terms
$J^{\prime}$ and $J^{\prime\prime}$ discussed in Section S3.3.3 commute and
anti-commute with $S$, respectively. Thus $J^{\prime}$ respects the flat-U(4)
symmetry whereas $J^{\prime\prime}$ breaks the flat-U(4) symmetries.) Since
$K$ (and $J^{\prime\prime}$) is small compared to other interaction parameters
(Table S6), flat-U(4) is still an approximate symmetry of the interaction
Hamiltonian. (Including the omitted small matrix elements in the interaction
Hamiltonian and the $\mathbf{k}$-dependence of these matrix elements may yield
more terms that break the flat-U(4). However, we claim these effects are weak
in the low energy physics.)

For the same reason as the flat-U(4) symmetry, $H_{I}$ will have the
chiral-U(4) symmetry as long as it commutes with $C$ defined in Eq. S86. (The
free part ${\hat{H}}_{0}$ anti-commute with $C$ when $v_{\star}^{\prime}=0$.)
One can easily verify that ${\hat{H}}_{U}$, ${\hat{H}}_{V}$, ${\hat{H}}_{W}$,
${\hat{H}}_{J}$, ${\hat{H}}_{\widetilde{J}}$, and ${\hat{H}}_{K}$ all commute
with $C$. Therefore, the chiral-U(4) symmetry is an exact symmetry of the
interaction Hamiltonian summarized in Section S3.6.1. Since for nonzero
$w_{0}/w_{1}$, the continuous model does not really respect the first chiral
symmetry, we expect that including the omitted small terms in ${\hat{H}}_{I}$
will break the chiral symmetry and hence the chiral-U(4). Nevertheless,
chiral-U(4) can be thought as an approximate symmetry of ${\hat{H}}_{I}$.

As discussed in Section S2.9, the presence of both flat-U(4) and chiral-U(4)
implies an SU(4)$\times$SU(4) symmetry. Given that 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 conserves the index $\xi$
($\xi=\eta(-1)^{\alpha-1}$ for $f$-electrons and $\xi=\eta(-1)^{a-1}$ for
$c$-electrons). Since ${\hat{H}}_{U}$, ${\hat{H}}_{V}$, ${\hat{H}}_{W}$,
${\hat{H}}_{J}$, ${\hat{H}}_{\widetilde{J}}$ all have both flat-U(4) and
chiral-U(4) symmetries and $\xi$ is repsected in all of them, they all have
the U(4)$\times$U(4) symmetry. However, as discussed in the last two
paragraphs ${\hat{H}}_{K}$ only has the chiral-U(4) symmetry. In Table S5 we
summarize the local symmetries of all the single-particle and interaction
Hamiltonian terms.

In Eq. S84 we defined the charge-conjugation symmetry

${\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}\
,$ (S265)

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})\ ,$ (S266)

and showed that the single-particle Hamiltonian (Section S2.5) respects the
charge conjugation symmetry. One can show
${\mathcal{P}}_{c}:\hat{\rho}(\mathbf{r}):{\mathcal{P}}_{c}^{-1}=-:\hat{\rho}(\mathbf{r}):$,
which automatically implies that the Coulomb interaction of the form

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

is invariant under the charge conjugation.

###  S3.7 Failure of the strong coupling picture in the first chiral limit of
the BM model ($w_{0}=0$)

The first chiral limit discussed in Section S2.6 is achieved by artificially
enforcing $v_{\star}^{\prime}=0$ in the single-particle Hamiltonian
${\hat{H}}_{0}$, which is obtained from the BM model with nonzero $w_{0}$. In
the actual chiral limit of the BM model, where $w_{0}=0$, the parameter
$v_{\star}^{\prime}$ automatically vanishes. In this subsection we discuss the
heavy fermion model in this actual chiral limit ($w_{0}=0$). We construct the
maximally localized Wannier functions by the same procedure as in Section
S2.1. Since for $w_{0}=0$ the $\Gamma_{3}$ states have energies about
$\pm$100meV, the previously used Wannierization energy window [$-80$meV,
$80$meV] now do not include the $\Gamma_{3}$ states and hence do not support
local Wannier functions. We hence change the energy window to [$-120$meV,
$120$meV] and use all the same other parameters (given around Eq. S42) for the
Wannierization. The resulted Wannier functions still decay exponentially but
are much less localized. The on-site repulsion $U_{1}$ of $f$-electrons
becomes extremely small ($<1$meV) while the hybridization $\gamma$ between
$f$\- and $c$-electrons becomes quite strong ($\sim$100meV). The property
$U_{1}\ll\gamma$ does not change if we use a larger Wannierization energy
window, e.g., [$-160$meV, $160$meV]. Therefore, in the actual first chiral
limit, the strong coupling picture ($U_{1}>\gamma$) is no longer valid, and
the quantum-dot-like behaviors observed in STMs Xie  _et al._ (2019); Wong
_et al._ (2020) cannot be easily explained.

Nevertheless, from a theoretical perspective, as explained in Section S2.10,
our model have the key geometric features of the BM model in the chiral limit:
The flat bands are analytically solvable, and they satisfy the ideal droplet
condition and have relatively flat Berry’s curvature.

##  S4 Correlated insulator states

###  S4.1 Mean field Hamiltonian

${\hat{H}}_{\widetilde{J}}$ will create two particles (or holes) in the local
orbitals, which will cost an energy of the order $U_{1}$. Since $U_{1}$ is the
largest energy scale of the problem, ${\hat{H}}_{\widetilde{J}}$ is an high-
energy process and will be omitted in the following. From a mean field aspect,
the energy contributed by the ${\hat{H}}_{\widetilde{J}}$ term is $\sim
J|\langle f^{\dagger}c\rangle|^{2}$. Since $U_{1}$ is much larger than the
coupling between $f$\- and $c$-electrons, $\langle f^{\dagger}c\rangle$ can be
treated as a small quantity and its second order terms, e.g., $J|\langle
f^{\dagger}c\rangle|^{2}$ can be omitted. We will also omit ${\hat{H}}_{K}$
for a similar reason. From the mean field aspect, the energy contributed by
${\hat{H}}_{K}$ is $\sim K\langle f^{\dagger}c\rangle\langle
c^{\dagger}c\rangle$. According to Section S3.5.2, the involved $\langle
c^{\dagger}c\rangle$ is off-diagonal in the $a$-index, e.g., $\langle
c^{\dagger}_{\mathbf{k}\overline{\alpha}\eta s}c_{\mathbf{k},\alpha+2,\eta
s}\rangle$, which must be small because only diagonal terms are non-vanishing
in the normal state. Thus $\sim K\langle f^{\dagger}c\rangle\langle
c^{\dagger}c\rangle$ can also be treated as a second order small quantity.
Therefore, in the following we approximate the interaction Hamiltonian as
${\hat{H}}_{I}\approx{\hat{H}}_{U}+{\hat{H}}_{W}+{\hat{H}}_{V}+{\hat{H}}_{J}$.

In the following we apply Hartree-Fock (HF) calculation to the remaining
terms. Let us denote the variational slater determinant state as
$|\Psi\rangle$. $|\Psi\rangle$ is supposed to minimize the total energy
$\langle\Psi|{\hat{H}}_{0}+{\hat{H}}_{I}|\Psi\rangle$. To simplify the
calculation, we rewrite the operator ${\hat{H}}_{I}$ as
${\hat{H}}_{I}=:{\hat{H}}_{I}:+{\hat{H}}_{I,MF}$, where $:{\hat{H}}_{I}:$ is
the normal ordered operator of ${\hat{H}}$ with respect to $|\Psi\rangle$ (not
$|G_{0}\rangle$) and ${\hat{H}}_{I,MF}$ is defined as
${\hat{H}}_{I}-:{\hat{H}}_{I}:$. $:{\hat{H}}_{I}:$ by definition will either
annihilate $|\Psi\rangle$ or excite two particles plus two holes upon
$|\Psi\rangle$. Then it follows that
$\langle\Psi|{\hat{H}}_{0}+{\hat{H}}_{I}|\Psi\rangle=\langle\Psi|{\hat{H}}_{0}+{\hat{H}}_{I,MF}|\Psi\rangle$.
Therefore, one only need to minimize
$\langle\Psi|{\hat{H}}_{0}+{\hat{H}}_{I,MF}|\Psi\rangle$.

In this subsection we explicitly derive ${\hat{H}}_{I,MF}$. For simplicity, we
assume that the ground state ($|\Psi\rangle$) preserves the translation
symmetry. For later convenience, we introduce the density matrices

$O^{f}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}=\langle\Psi|f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}|\Psi\rangle=\frac{1}{N}\sum_{\mathbf{k}\in{\rm
MBZ}}\langle\Psi|f_{\mathbf{k}\alpha\eta
s}^{\dagger}f_{\mathbf{k}\alpha^{\prime}\eta^{\prime}s^{\prime}}|\Psi\rangle\
,$ (S268) $O^{c}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}=\frac{1}{N}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(\langle\Psi|c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s^{\prime}}|\Psi\rangle-\frac{1}{2}\delta_{aa^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}\right)$
(S269)
$O^{cf}_{a\eta_{2}s_{2},\alpha\eta_{1}s_{1}}=\frac{1}{\sqrt{N}}\sum_{|\mathbf{k}|<\Lambda_{c}}e^{-i\mathbf{k}\cdot\mathbf{R}}\langle\Psi|c^{\dagger}_{\mathbf{k}a\eta_{2}s_{2}}f_{\mathbf{R}\alpha\eta_{1}s_{1}}|\Psi\rangle\
.$ (S270)

We have defined $O^{c}_{a\eta s,a^{\prime}\eta^{\prime}s^{\prime}}$ as the
expectation of the conduction band density operator with respect to the charge
neutrality point such that it will remain finite even we take the cutoff
$\Lambda_{c}$ to infinity. The filling of local orbitals and conduction bands
are given by

$\nu_{f}=\mathrm{Tr}[O^{f}]-4,\qquad\nu_{c}=\mathrm{Tr}[O^{c}]\ .$ (S271)

We also assume that there is no pairing in $|\Psi\rangle$. The total filling
with respect to the charge neutrality point is given by $\nu=\nu_{f}+\nu_{c}$.
For the normal state $|G_{0}\rangle$ (defined in Section S3.1) at charge
neutrality point, there is $\nu=\nu_{f}=\nu_{c}=0$.

We first decouple the density-density interactions of local orbitals, i.e.,
${\hat{H}}_{U}$ in Eq. S256. To derive the mean field Hamiltonian, we first
divide the density-density interaction into constant, bilinear, and quartic
terms as

$\displaystyle\hat{H}_{U}=$ $\displaystyle
8NU_{1}+48NU_{2}-(4U_{1}+24U_{2})\sum_{\mathbf{R}}\sum_{\alpha\eta
s}f_{\mathbf{R}\alpha\eta s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}+\frac{U_{1}}{2}\sum_{\mathbf{R}}\sum_{\begin{subarray}{c}\alpha\eta s\\\
\beta\eta^{\prime}s^{\prime}\end{subarray}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}$
$\displaystyle+\frac{U_{2}}{2}\sum_{\langle\mathbf{R}\mathbf{R}^{\prime}\rangle}\sum_{\begin{subarray}{c}\alpha\eta
s\\\ \beta\eta^{\prime}s^{\prime}\end{subarray}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}$
$\displaystyle=$ $\displaystyle
8NU_{1}+48NU_{2}-(3.5U_{1}+24U_{2})\sum_{\mathbf{R}}\sum_{\alpha\eta
s}f_{\mathbf{R}\alpha\eta s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}+\frac{U_{1}}{2}\sum_{\mathbf{R}}\sum_{\begin{subarray}{c}\alpha\eta s\\\
\beta\eta^{\prime}s^{\prime}\end{subarray}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}$
$\displaystyle+\frac{U_{2}}{2}\sum_{\langle\mathbf{R}\mathbf{R}^{\prime}\rangle}\sum_{\begin{subarray}{c}\alpha\eta
s\\\ \beta\eta^{\prime}s^{\prime}\end{subarray}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}$ (S272)

The mean field Hamiltonian is obtained by projecting out the terms that create
two particles and two holes upon the ground state $|\Psi\rangle$. According to
the Wick’s theorem we have

$\displaystyle f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}=:f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}:+\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}:f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}:\c{1}f_{\mathbf{R}\alpha\eta
s}}+\wick$ $\displaystyle+$ $\displaystyle\wick{:\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}:}+\wick{:f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}:}+\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{2}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{2}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}}+\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{2}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}\c{2}f_{\mathbf{R}\alpha\eta
s}}$ (S273)

Here $:A:$ is the normal ordered form of $A$ with respect to $|\Psi\rangle$
(not $|G_{0}\rangle$ defined in Section S3.1). The first term will either
annihilate $|\Psi\rangle$ or excite two particle-hole pairs on top of
$|\Psi\rangle$. We will omit the first term and approximate the four fermion
operator by bilinear terms and constant terms

$\displaystyle f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}\approx\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}}+\wick{f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}}+\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}}$ $\displaystyle+$ $\displaystyle\wick{f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}}-\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{2}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{2}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}}-\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{2}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}\c{2}f_{\mathbf{R}\alpha\eta
s}}\ .$ (S274)

Similarly, there are

$\displaystyle f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}\approx\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}}+\wick{f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{1}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}}+\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}}$ $\displaystyle+$ $\displaystyle\wick{f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{1}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}}-\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{2}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{2}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}\c{1}f_{\mathbf{R}\alpha\eta
s}}-\wick{\c{1}f_{\mathbf{R}\alpha\eta
s}^{\dagger}\c{2}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}\c{1}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}\c{2}f_{\mathbf{R}\alpha\eta
s}}\ .$ (S275)

Since the on-site repulsion $U_{1}$ is the largest energy scale and dominates
the correlation physics, we assume that only on-site order parameters, i.e.,
$\langle f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}\rangle$, are
formed. Then we can write the approximate on-site and off-site quartic terms
as

$\displaystyle\sum_{\begin{subarray}{c}\alpha\eta s\\\
\beta\eta^{\prime}s^{\prime}\end{subarray}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}\approx$ $\displaystyle 2\sum_{\alpha\eta
s}\mathrm{Tr}[O^{f}]f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}-2\sum_{\begin{subarray}{c}\alpha\eta
s\\\
\beta\eta^{\prime}s^{\prime}\end{subarray}}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\alpha\eta
s}O^{f}_{\alpha\eta
s,\beta\eta^{\prime}s^{\prime}}-\mathrm{Tr}[O^{f}]\mathrm{Tr}[O^{f}]+\mathrm{Tr}[O^{f}O^{f}]$
(S276)

and

$\displaystyle\sum_{\begin{subarray}{c}\alpha\eta s\\\
\beta\eta^{\prime}s^{\prime}\end{subarray}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}\approx$
$\displaystyle\sum_{\beta\eta^{\prime}s^{\prime}}\mathrm{Tr}[O^{f}]f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}^{\prime}\beta\eta^{\prime}s^{\prime}}+\sum_{\alpha\eta
s}\mathrm{Tr}[O^{f}]f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}-\mathrm{Tr}[O^{f}]\mathrm{Tr}[O^{f}]\
,$ (S277)

respectively. Substituting these equations into the interaction Hamiltonian
${\hat{H}}_{U}$, we obtain

${\hat{H}}_{U}\approx{\hat{H}}_{U,MF}=\overline{H}_{U}-E_{U}$ (S278)

with

$\displaystyle\overline{H}_{U}=\sum_{\mathbf{R}}\sum_{\alpha\eta
s}\left(U_{1}(\nu_{f}+0.5)+6U_{2}\nu_{f}\right)f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}-U_{1}\sum_{\mathbf{R}}O^{f}_{\alpha\eta
s,\beta\eta^{\prime}s^{\prime}}f_{\mathbf{R}\beta\eta^{\prime}s^{\prime}}^{\dagger}f_{\mathbf{R}\alpha\eta
s}\ ,$ (S279)

and

$E_{U}=\frac{NU_{1}}{2}\left(\nu_{f}^{2}+8\nu_{f}-\mathrm{Tr}[O^{f}O^{f}]\right)+3NU_{2}(\nu_{f}^{2}+8\nu_{f}).$
(S280)

Following the same logic as decoupling ${\hat{H}}_{U}$, one can decouple
${\hat{H}}_{W}$ as

$\displaystyle{\hat{H}}_{W}\approx$
$\displaystyle\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{a\eta
s}W_{a}\nu_{f}:c_{\mathbf{k}a\eta s}^{\dagger}c_{\mathbf{k}a\eta
s}:+\sum_{\mathbf{R}a\alpha\eta s}W_{a}:f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta
s}:\nu_{c,a}-N\sum_{a}W_{a}\nu_{f}\nu_{c,a}$
$\displaystyle-\frac{1}{\sqrt{N}}\sum_{\mathbf{R}\alpha\eta_{1}s_{1}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\eta_{2}s_{2}a}W_{a}\left(O^{cf}_{a\eta_{2}s_{2},\alpha\eta_{1}s_{1}}e^{i\mathbf{k}\cdot\mathbf{R}}f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}a\eta_{2}s_{2}}+h.c.\right)-N\sum_{a\alpha}\sum_{\eta_{2}\eta_{1}s_{2}s_{1}}|O^{cf}_{a\eta_{2}s_{2},\alpha\eta_{1}s_{1}}|^{2}W_{a}\
,$ (S281)

where

$\nu_{c,a}=\sum_{\eta s}O^{c}_{a\eta s,a\eta s}\ .$ (S282)

In the first term we have omitted the Umklapp scatterings, i.e.,
$W_{a}\nu_{f}:c_{\mathbf{k}+\mathbf{G}a\eta s}^{\dagger}c_{\mathbf{k}a\eta
s}:$ for $\mathbf{G}\neq 0$, because the state $c_{\mathbf{k}+\mathbf{G}a\eta
s}^{\dagger}$ has a huge kinetic energy and the Umklapp scattering will be
weak. We can further organize the terms as

${\hat{H}}_{W}\approx{\hat{H}}_{W,MF}=\overline{H}_{W}-E_{W}\ ,$ (S283)

with

$\displaystyle\overline{H}_{W}=$
$\displaystyle\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{a\eta
s}\nu_{f}W_{a}c_{\mathbf{k}a\eta s}^{\dagger}c_{\mathbf{k}a\eta
s}+\sum_{\mathbf{R}}\sum_{\alpha\eta s}\nu_{c,a}W_{a}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}$
$\displaystyle-\frac{1}{\sqrt{N}}\sum_{\mathbf{R}\alpha\eta_{1}s_{1}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\eta_{2}s_{2}a}W_{a}\left(O^{cf}_{a\eta_{2}s_{2},\alpha\eta_{1}s_{1}}e^{i\mathbf{k}\cdot\mathbf{R}}f_{\mathbf{R}\alpha\eta_{1}s_{1}}^{\dagger}c_{\mathbf{k}a\eta_{2}s_{2}}+h.c.\right)$
(S284)
$E_{W}=N\sum_{a}\left(W_{a}\nu_{f}\nu_{c,a}+4W_{a}\nu_{c,a}\right)+\sum_{a}\sum_{|\mathbf{k}|<\Lambda_{c}}2W_{a}\nu_{f}-N\sum_{a\alpha}\sum_{\eta_{2}\eta_{1}s_{2}s_{1}}|O^{cf}_{a\eta_{2}s_{2},\alpha\eta_{1}s_{1}}|^{2}W_{a}\
.$ (S285)

We emphasize that the term $\sum_{|\mathbf{k}|<\Lambda_{c}}2W_{a}\nu_{f}$
depends on the cutoff ($\Lambda_{c}$) of the conduction bands and diverges as
$\Lambda_{c}\to\infty$. However, the total energy contributed by
${\hat{H}}_{W}$ will converge as $\Lambda_{c}\to\infty$ because the divergent
part of $\sum_{|\mathbf{k}|<\Lambda_{c}}2W_{a}\nu_{f}$ will be canceled by the
divergent part of the first term of $\overline{H}_{W}$.

Then we decouple the exchange interaction ${\hat{H}}_{J}$. Following the same
method as for ${\hat{H}}_{U}$, ${\hat{H}}_{J}$ can be approximated as

$\displaystyle{\hat{H}}_{J}\approx$ $\displaystyle-J\sum_{\mathbf{R}\alpha\eta
s_{1}s_{2}}:f_{\mathbf{R}\alpha\eta s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta
s_{2}}:O^{c}_{\alpha+2\eta s_{2},\alpha+2\eta s_{1}}-J\sum_{\alpha\eta
s_{1}s_{2}}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(O^{f}_{\alpha\eta
s_{1},\alpha\eta
s_{2}}-\frac{1}{2}\delta_{s_{1}s_{2}}\right):c_{\mathbf{k}\alpha+2\eta
s_{2}}^{\dagger}c_{\mathbf{k}\alpha+2\eta s_{1}}:$
$\displaystyle+NJ\sum_{\alpha\eta s_{1}s_{2}}\left(O^{f}_{\alpha\eta
s_{1},\alpha\eta
s_{2}}-\frac{1}{2}\delta_{s_{1}s_{2}}\right)O^{c}_{\alpha+2\eta
s_{2},\alpha+2\eta s_{1}}+J\sum_{\alpha\eta
s_{1}s_{2}}\sum_{|\mathbf{k}|<\Lambda_{c}}O^{f}_{\alpha\eta
s_{1},\overline{\alpha}-\eta
s_{2}}c_{\mathbf{k},\overline{\alpha}+2,-\eta,s_{2}}^{\dagger}c_{\mathbf{k},\alpha+2,\eta,s_{1}}$
$\displaystyle+J\sum_{\alpha\eta
s_{1}s_{2}}\sum_{\mathbf{R}}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}f_{\mathbf{R}\overline{\alpha}-\eta
s_{2}}O_{\overline{\alpha}+2,-\eta,s_{2}\ ;\
\alpha+2,\eta,s_{1}}^{c}-NJ\sum_{\alpha\eta s_{1}s_{2}}O^{f}_{\alpha\eta
s_{1},\overline{\alpha}-\eta s_{2}}O_{\overline{\alpha}+2,-\eta,s_{2}\ ;\
\alpha+2,\eta,s_{1}}^{c}$
$\displaystyle+\frac{J}{\sqrt{N}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha\eta
s_{1}s_{2}}\left(e^{i\mathbf{k}\cdot\mathbf{R}}O^{cf}_{\alpha+2\eta
s_{2},\alpha\eta s_{2}}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}c_{\mathbf{k},\alpha+2,\eta,s_{1}}+h.c.\right)-NJ\sum_{\alpha\eta
s_{1}s_{2}}O^{cf}_{\alpha+2\eta s_{2},\alpha\eta s_{2}}O^{cf*}_{\alpha+2\eta
s_{1},\alpha\eta s_{1}}$
$\displaystyle-\frac{J}{\sqrt{N}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha\eta
s_{1}s_{2}}\left(e^{i\mathbf{k}\cdot\mathbf{R}}O^{cf}_{\overline{\alpha}+2,-\eta,s_{2}\
;\ \overline{\alpha},-\eta,s_{2}}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}c_{\mathbf{k},\alpha+2,\eta,s_{1}}+h.c.\right)+NJ\sum_{\alpha\eta
s_{1}s_{2}}O^{cf}_{\overline{\alpha}+2,-\eta,s_{2}\ ;\
\overline{\alpha},-\eta,s_{2}}O^{cf*}_{\alpha+2\eta s_{1},\alpha\eta s_{1}}\
.$ (S286)

We can organize the terms as

${\hat{H}}_{J}\approx{\hat{H}}_{J,MF}=\overline{H}_{J}-E_{J}\ ,$ (S287)

with

$\displaystyle\overline{H}_{J}=$ $\displaystyle-J\sum_{\mathbf{R}\alpha\eta
s_{1}s_{2}}f_{\mathbf{R}\alpha\eta s_{1}}^{\dagger}f_{\mathbf{R}\alpha\eta
s_{2}}O^{c}_{\alpha+2\eta s_{2},\alpha+2\eta s_{1}}-J\sum_{\alpha\eta
s_{1}s_{2}}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(O^{f}_{\alpha\eta
s_{1},\alpha\eta
s_{2}}-\frac{1}{2}\delta_{s_{1}s_{2}}\right)c_{\mathbf{k}\alpha+2\eta
s_{2}}^{\dagger}c_{\mathbf{k}\alpha+2\eta s_{1}}$
$\displaystyle+J\sum_{\alpha\eta
s_{1}s_{2}}\sum_{|\mathbf{k}|<\Lambda_{c}}O^{f}_{\alpha\eta
s_{1},\overline{\alpha}-\eta
s_{2}}c_{\mathbf{k},\overline{\alpha}+2,-\eta,s_{2}}^{\dagger}c_{\mathbf{k},\alpha+2,\eta,s_{1}}+J\sum_{\alpha\eta
s_{1}s_{2}}\sum_{\mathbf{R}}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}f_{\mathbf{R}\overline{\alpha}-\eta
s_{2}}O_{\overline{\alpha}+2,-\eta,s_{2}\ ;\ \alpha+2,\eta,s_{1}}^{c}$
$\displaystyle+\frac{J}{\sqrt{N}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha\eta
s_{1}s_{2}}\left(e^{i\mathbf{k}\cdot\mathbf{R}}O^{cf}_{\alpha+2\eta
s_{2},\alpha\eta s_{2}}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}c_{\mathbf{k},\alpha+2,\eta,s_{1}}+h.c.\right)$
$\displaystyle-\frac{J}{\sqrt{N}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha\eta
s_{1}s_{2}}\left(e^{i\mathbf{k}\cdot\mathbf{R}}O^{cf}_{\overline{\alpha}+2,-\eta,s_{2}\
;\ \overline{\alpha},-\eta,s_{2}}f_{\mathbf{R}\alpha\eta
s_{1}}^{\dagger}c_{\mathbf{k},\alpha+2,\eta,s_{1}}+h.c.\right)\ ,$ (S288)

and

$\displaystyle E_{J}=$ $\displaystyle-NJ\sum_{\alpha\eta
s_{1}s_{2}}O^{f}_{\alpha\eta s_{1},\alpha\eta s_{2}}O^{c}_{\alpha+2\eta
s_{2},\alpha+2\eta s_{1}}-J\sum_{\alpha\eta
s}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(O^{f}_{\alpha\eta s,\alpha\eta
s}-\frac{1}{2}\right)\frac{1}{2}+NJ\sum_{\alpha\eta
s_{1}s_{2}}O^{f}_{\alpha\eta s_{1},\overline{\alpha}-\eta
s_{2}}O_{\overline{\alpha}+2,-\eta,s_{2}\ ;\ \alpha+2,\eta,s_{1}}^{c}$
$\displaystyle+NJ\sum_{\alpha\eta s_{1}s_{2}}O^{cf}_{\alpha+2\eta
s_{2},\alpha\eta s_{2}}O^{cf*}_{\alpha+2\eta s_{1},\alpha\eta
s_{1}}-NJ\sum_{\alpha\eta s_{1}s_{2}}O^{cf}_{\overline{\alpha}+2,-\eta,s_{2}\
;\ \overline{\alpha},-\eta,s_{2}}O^{cf*}_{\alpha+2\eta s_{1},\alpha\eta
s_{1}}\ .$ (S289)

The second term in $E_{J}$ diverges as $\Lambda_{c}\to\infty$. However, the
total energy contributed by ${\hat{H}}_{J}$ will converge as
$\Lambda_{c}\to\infty$ because the divergent part of the second term in
$E_{J}$ will be canceled by the (divergent part of) second term of
$\overline{H}_{J}$. We can equivalently write ${\hat{H}}_{J}$ in a more
compact form as

$\displaystyle{\hat{H}}_{J}=$
$\displaystyle-\frac{J}{2}\sum_{\mathbf{R}}\sum_{\alpha_{1}\alpha_{2}\eta_{1}\eta_{2}s_{1}s_{2}}(\eta_{1}\eta_{2}+(-1)^{\alpha_{1}+\alpha_{2}})f_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}^{\dagger}f_{\mathbf{R}\alpha_{2}\eta_{2}s_{2}}O^{c}_{\alpha_{2}+2\eta_{2}s_{2},\alpha_{1}+2\eta_{1}s_{1}}$
$\displaystyle-\frac{J}{2}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha_{1}\alpha_{2}\eta_{1}\eta_{2}s_{1}s_{2}}(\eta_{1}\eta_{2}+(-1)^{\alpha_{1}+\alpha_{2}})c_{\mathbf{k}\alpha_{2}+2,\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}\alpha_{1}+2,\eta_{1}s_{1}}\left(O^{f}_{\alpha_{1}\eta_{1}s_{1},\alpha_{2}\eta_{2}s_{2}}-\frac{1}{2}\delta_{\alpha_{1}\alpha_{2}}\delta_{\eta_{1}\eta_{2}}\delta_{s_{1}s_{2}}\right)$
$\displaystyle+\frac{J}{2\sqrt{N}}\sum_{\mathbf{R}}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha_{1}\alpha_{2}\eta_{1}\eta_{2}s_{1}s_{2}}e^{i\mathbf{k}\cdot\mathbf{R}}(\eta_{1}\eta_{2}+(-1)^{\alpha_{1}+\alpha_{2}})\left(O^{cf}_{\alpha_{2}+2\eta_{2}s_{2},\alpha_{2}\eta_{2}s_{2}}f^{\dagger}_{\mathbf{R}\alpha_{1}\eta_{1}s_{1}}c_{\mathbf{k}\alpha_{1}+2\eta_{1}s_{1}}+h.c.\right)\
.$ (S290)

In the end we decouple the Coulomb interaction of the conduction bands.
Usually, it is the Fock term that leads to symmetry breaking in the inner
degrees of freedom. We claim that the symmetry breaking mainly come from the
on-site repulsion of the local orbitals and the interaction of conduction
bands only renormalize the dispersion of the conduction bands. Therefore, for
simplicity, we will omit the Fock channel of this interaction. The Hartree-
mean field is

${\hat{H}}_{V}\approx{\hat{H}}_{V,MF}=\overline{H}_{V}-E_{V}\ ,$ (S291)

with

$\overline{H}_{V}=\frac{V(0)}{\Omega_{0}}\nu_{c}\sum_{\eta
sa}\sum_{|\mathbf{k}|<\Lambda_{c}}c^{\dagger}_{\mathbf{k}a\eta
s}c_{\mathbf{k}a\eta s}$ (S292)

and

$E_{V}=\frac{V(0)}{2\Omega_{0}}N\nu_{c}^{2}+\frac{V(0)}{\Omega_{0}}\sum_{|\mathbf{k}|<\Lambda_{c}}8\nu_{c}\
.$ (S293)

$\overline{H}_{V}$ only shifts the energy the conduction bands.

In summary, the mean field Hamiltonian of the interaction is given by

${\hat{H}}_{I}\approx{\hat{H}}_{I,MF}=\overline{H}_{U}+\overline{H}_{W}+\overline{H}_{V}+\overline{H}_{J}-E_{U}-E_{W}-E_{V}-E_{J}\
.$ (S294)

###  S4.2 Correlated insulator phases at the filling $\nu=0$

Fillings | $\nu=0$ | $\nu=-1$
---|---|---
Phases | VP | K-IVC | IVC | VP | K-IVC | IVC | VP + K-IVC
One-shot HF energy (meV) | -698.645 | -698.702 | -694.888 | -666.326 | -666.359 | -664.607 | -666.360
SCF HF energy (meV) | -704.019 | -704.109 | -700.616 | -671.809 | -671.857 | -670.330 | -671.860
Fillings | $\nu=-2$ | $\nu=-3$
Phases | VP | K-IVC | IVC | VP | K-IVC |  |
One-shot HF energy (meV) | -575.492 | -575.507 | -574.528 | -434.015 | -434.015 |  |
SCF HF energy (meV) | -578.975 | -579.036 | -578.203 | -426.129 | -426.113 |  |
Table S7: The one-shot and self-consistent HF energies of the correlated
insulator phases at integer fillings. Since the energy difference between
flat-U(4) rotation related states are small ($\lesssim 0.1$meV) there may be
strong order-parameter fluctuations at finite temperatures.  Figure S6: HF
band structures of correlated insulator phases at $\nu=0$. (a), (b), (c) are
the one-shot HF band structures of the VP, K-IVC, and IVC phases,
respectively. (d), (e), (f) are the self-consistent HF band structures of the
VP, K-IVC, and IVC phases, respectively. The color represents the composition
of the energy bands, where yellow corresponds to the local orbitals and blue
corresponds to the conduction bands. We have chosen $w_{0}/w_{1}=0.8$ in the
calculation. Other parameters of the single-particle and interaction
hamiltonians are given in Tables S4 and S6.

####  S4.2.1 Numerical results

We consider the parent wave function of valley polarized (VP) state as

$|\mathrm{VP}_{0}^{\nu=0}\rangle=\prod_{\mathbf{R}}f_{\mathbf{R}1+\uparrow}^{\dagger}f_{\mathbf{R}1+\downarrow}^{\dagger}f_{\mathbf{R}2+\uparrow}^{\dagger}f_{\mathbf{R}2+\downarrow}^{\dagger}|\mathrm{FS}\rangle\
.$ (S295)

Here $|\rm FS\rangle$ is the Fermi sea state occupying the lower two bands
(per spin per valley) of $H^{(c,\eta)}(\mathbf{k})$ (Eq. S82), i.e.,

$|\mathrm{FS}\rangle=\prod_{\eta
s}\prod_{n=1,2}\prod_{|\mathbf{k}|<\Lambda_{c}}\left(\sum_{a}\mathcal{U}_{a,n}^{(\eta)}c_{\mathbf{k}a\eta
s}^{\dagger}\right)|0\rangle\ ,$ (S296)

where $\mathcal{U}_{a,n}^{(\eta)}$ ($n=1,2,3,4$) is the $n$-th eigenvector of
$H^{(c,\eta)}(\mathbf{k})$. At $\mathbf{k}=0$, the second and third
eigenstates are degenerate, we should choose $\mathcal{U}_{a,2}^{(\eta)}(0)$
as the vector that smoothly connects to
$\mathcal{U}_{a,2}^{(\eta)}(\mathbf{k}\to 0)$. By definition, the yielded
$O^{c}$ is diagonal in valley and spin indices

$\displaystyle O^{c}_{a\eta s,a^{\prime}\eta^{\prime}s^{\prime}}=$
$\displaystyle\frac{1}{N}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(\langle\mathrm{FS}|c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}^{\prime}a^{\prime}\eta^{\prime}s^{\prime}}|\mathrm{FS}\rangle-\frac{1}{2}\delta_{aa^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}\right)=\frac{1}{N}\delta_{ss^{\prime}}\delta_{\eta\eta^{\prime}}\sum_{|\mathbf{k}|<\Lambda_{c}}\left(\left(\sum_{n=1,2}\mathcal{U}_{a,n}^{(\eta)}(\mathbf{k})\mathcal{U}^{(\eta)*}_{a^{\prime},n}(\mathbf{k})\right)-\frac{1}{2}\delta_{aa^{\prime}}\right)\
.$ (S297)

Given $O^{c}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}=\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}O^{c}_{a\eta
s,a^{\prime}\eta s}$, only the diagonal elements $O^{c}_{a\eta s,a\eta s}$
will enter the mean field Hamiltonian in Sections S4.1, S4.1 and S292. Now we
show that all the diagonal elements are zero. Due to the particle-hole
symmetry $P$ (Eq. S65), we can always pick the gauge
$(-1)^{a-1}\mathcal{U}_{a,n}^{(\eta)}(\mathbf{k})=\mathcal{U}_{a,5-n}^{(\eta)}(-\mathbf{k})$.
Using this condition and the completeness
$\sum_{n=1}^{4}\mathcal{U}_{n}^{(\eta)}(\mathbf{k})\mathcal{U}^{(\eta)\dagger}_{n}(\mathbf{k})=\mathbb{I}_{4\times
4}$, one can show that

$O^{c}_{a\eta s,a^{\prime}\eta s}+(-1)^{a+a^{\prime}}O^{c}_{a\eta
s,a^{\prime}\eta s}=0\ ,$ (S298)

which implies $O^{c}_{a\eta s,a\eta s}=0$. Therefore, $O^{c}$ does not enter
the mean field Hamiltonian because it is valley-diagonal and all its diagonal
elements vanish. The density matrix $O^{cf}$ (Eq. S270) given by the parent VP
state also vanishes because the state respects the particle numbers of local
orbitals and conduction bands separately. The density matrix $O^{f}$ (Eq.
S268) given by the parent VP state is nonzero. We will discuss it in detail in
next subsection. Substituting these density matrices into the mean field
equations, we can obtain the mean field approximation of the total Hamiltonian
as

${\hat{H}}_{0}+{\hat{H}}_{I}\approx{\hat{H}}_{0}+\overline{H}_{U}+\overline{H}_{W}+\overline{H}_{V}+\overline{H}_{J}-E_{U}-E_{W}-E_{V}-E_{J}\
,$ (S299)

where $\overline{H}_{X}$ and $E_{X}$ ($X=U,W,V,J$) are functions the density
matrices (Section S4.1). We refer to the ground state energy of the above mean
field Hamiltonian as the “one-shot” HF energy (Table S7). The spectrum given
by ${\hat{H}}_{0}+\sum_{X}\overline{H}_{X}$ is referred to as the “one-shot”
HF spectrum (Fig. S6(a)). The next step wave function, denoted as
$|\mathrm{VP}_{1}^{\nu=0}\rangle$, occupies the one-shot bands to the filling
of the charge neutrality point. One can then calculate the density matrices of
$|\mathrm{VP}_{1}^{\nu=0}\rangle$ and generate a new mean field Hamiltonian.
The hybridization $O^{cf}$ given by $|\mathrm{VP}_{1}^{\nu=0}\rangle$ is in
general nonzero. Repeating this process until the ground state converges, one
achieve the self-consistent HF ground state, whose energy is called the self-
consistent HF energy. We denote the self-consistent wave function as
$|\mathrm{VP}_{\infty}^{\nu=0}\rangle$ The (one-shot and self-consistent) HF
energies and HF spectra of the VP state are summarized in Table S7 and Fig.
S6, respectively.

In the process of self-consistent iterations, the one-shot spectrum in Fig.
S6(a) continuously changes to the self-consistent spectrum in Fig. S6(d).
Thus, the state $|\mathrm{VP}_{1}^{\nu=0}\rangle$, which occupies the one-shot
bands to the filling $\nu=0$, is adiabatically connected to the self-
consistent state $|\mathrm{VP}_{\infty}^{\nu=0}\rangle$. However, the parent
VP state $|\mathrm{VP}_{0}^{\nu=0}\rangle$ may not adiabatically connected to
them because it occupies the local orbitals and conduction separately but not
the quasi-particle bands. A key difference between
$|\mathrm{VP}_{1}^{\nu=0}\rangle$ and $|\mathrm{VP}_{0}^{\nu=0}\rangle$ is
that $O^{cf}$ given by $|\mathrm{VP}_{1}^{\nu=0}\rangle$ is in general nonzero
while $O^{cf}$ given by $|\mathrm{VP}_{0}^{\nu=0}\rangle$ is zero. In this
work, we refer to $|\mathrm{VP}_{0}^{\nu=0}\rangle$,
$|\mathrm{VP}_{1}^{\nu=0}\rangle$, and $|\mathrm{VP}_{\infty}^{\nu=0}\rangle$
as the parent VP, one-shot VP, and self-consistent VP states, respectively.

The parent wave functions for the inter-valley-coherent (IVC) and Kramers
inter-valley-coherent (K-IVC) states are rotated from
$|\mathrm{VP}_{0}^{\nu=0}\rangle$ by the chiral-U(4) (Eq. S92) and flat-U(4)
(Eq. S105) symmetries as

$|\mathrm{IVC}_{0}^{\nu=0}\rangle=e^{-i\frac{\pi}{2}{\hat{\Theta}}_{x0}}|\mathrm{VP}_{0}^{\nu=0}\rangle=\prod_{\mathbf{R}}\frac{1}{4}(f_{\mathbf{R}1+\uparrow}^{\dagger}-if_{\mathbf{R}2-\uparrow}^{\dagger})(f_{\mathbf{R}1+\downarrow}^{\dagger}-if_{\mathbf{R}2-\downarrow}^{\dagger})(-if_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})(-if_{\mathbf{R}1-\downarrow}^{\dagger}+f_{\mathbf{R}2+\downarrow}^{\dagger})|\rm
FS\rangle\ ,$ (S300)

and

$|\text{K-IVC}_{0}^{\nu=0}\rangle=e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}|\mathrm{VP}_{0}^{\nu=0}\rangle=\prod_{\mathbf{R}}\frac{1}{4}(f_{\mathbf{R}1+\uparrow}^{\dagger}+f_{\mathbf{R}2-\uparrow}^{\dagger})(f_{\mathbf{R}1+\downarrow}^{\dagger}+f_{\mathbf{R}2-\downarrow}^{\dagger})(-f_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})(-f_{\mathbf{R}1-\downarrow}^{\dagger}+f_{\mathbf{R}2+\downarrow}^{\dagger})|\rm
FS\rangle\ ,$ (S301)

respectively. One can verify that the parent IVC state respects the spinless
time-reversal symmetry $T$ (Eq. S36) whereas the parent K-IVC state breaks
$T$. Substituting the parent wave functions into the HF loop, we obtain the
one-shot and self-consistent HF energies (Table S7). The one-shot IVC and
K-IVC states, i.e., $|\mathrm{IVC}^{\nu=0}_{1}\rangle$ and
$|\text{K-IVC}^{\nu=0}_{1}\rangle$, are obtained as Fock states occupying the
corresponding one-shot HF bands to the filling $\nu=0$. The one-shot K-IVC
state has a lower energy.

As shown in Table S7, the energy difference of K-IVC and VP states is only of
the order 0.1meV whereas the energy difference between the IVC state and VP is
of the order 3meV. Since K-IVC and IVC are rotated from VP by the flat-U(4)
and chiral-U(4), respectively, the result suggests that the flat-U(4) symmetry
is a better approximation of the considered Hamiltonian
${\hat{H}}_{0}+\sum_{X}{\hat{H}}_{X}$ ($X=U,W,J,V$). This observation is
consistent with the argument at the end of Section S2.7.

####  S4.2.2 Analytical analysis of the one-shot HF

The one-shot and self-consistent HF band structures of the VP, K-IVC, and IVC
states at $\nu=0$ are summarized in Fig. S6. We can see that the self-
consistent band structures are very close to the one-shot HF band structures.
Here we study the one-shot spectra of the VP and K-IVC states around the
charge neutrality point. We do not study the IVC state here because it has a
higher energy in the mean field. The density matrices given by
$|\text{VP}^{\nu=0}_{0}\rangle$ are

$O^{f}=\sigma_{0}(\frac{\tau_{z}+\tau_{0}}{2}){\varsigma}_{0},\qquad
O^{c}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}=\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}O^{c}_{a\eta
s,a^{\prime}\eta s},\qquad O^{c}_{a\eta s,a\eta s}=0,\qquad O^{cf}=0\ .$
(S302)

Substituting these density matrices into the mean field equations around Eqs.
S279, S4.1, S4.1 and S292, we obtain the one-shot HF mean field and energies
as

$\overline{H}_{U}=-\frac{U_{1}}{2}\sum_{\mathbf{R}}\sum_{\alpha\eta s}\eta
f_{\mathbf{R}\alpha\eta s}^{\dagger}f_{\mathbf{R}\alpha\eta
s},\qquad\overline{H}_{W}=\overline{H}_{V}=0,\qquad\overline{H}_{J}=-\frac{J}{2}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha\eta
s}\eta c_{\mathbf{k}\alpha+2\eta s}^{\dagger}c_{\mathbf{k}\alpha+2\eta s}$
(S303)

and

$E_{U}=-2NU_{1},\qquad E_{W}=E_{V}=E_{J}=0\ ,$ (S304)

respectively. The low energy band structure can be described by the k$\cdot$p
expansion of ${\hat{H}}_{0}+\overline{H}_{U}+\overline{H}_{J}$

$H_{\rm
VP}^{\rm(MF)}(\mathbf{k})\approx\begin{pmatrix}0&v_{\star}(k_{x}\sigma_{0}\tau_{z}+ik_{y}\sigma_{z}\tau_{0})&\gamma\sigma_{0}\tau_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0})\\\
v_{\star}(k_{x}\sigma_{0}\tau_{z}-ik_{y}\sigma_{z}\tau_{0})&M\sigma_{x}\tau_{0}-\frac{J}{2}\sigma_{0}\tau_{z}&0\\\
\gamma\sigma_{0}\tau_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0})&0&-\frac{U_{1}}{2}\sigma_{0}\tau_{z}\end{pmatrix}\otimes{\varsigma}_{0}\
.$ (S305)

For large $\mathbf{k}$ ($v_{\star}|\mathbf{k}|\gg\frac{1}{2}U_{1}$), the low
energy bands are contributed by the local orbitals, which have the energies
$\pm U_{1}/2$. We can see this in Fig. S6 at large $\mathbf{k}$. For
$\mathbf{k}=0$, the energy eigenvalues are given by

$\pm
M\pm\frac{J}{2},\qquad\pm\frac{U_{1}}{4}\pm\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\
.$ (S306)

In Fig. S6(a), (d), the 24 energy levels closest to Fermi level from low to
high are

$-\frac{U_{1}}{4}-\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold}),\qquad\frac{U_{1}}{4}-\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold}),\qquad-M-\frac{J}{2}\;(\text{2-fold}),\qquad
M-\frac{J}{2}\;(\text{2-fold}),$ (S307)
$-M+\frac{J}{2}\;(\text{2-fold}),\qquad
M+\frac{J}{2}\;(\text{2-fold}),\qquad-\frac{U_{1}}{4}+\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold}),\qquad\frac{U_{1}}{4}+\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold})\
.$ (S308)

Notice that in Fig. S6(a), (d) every single band is spin degenerate.

The density matrices given by $|\text{K-IVC}^{\nu=0}_{0}\rangle$ are

$O^{f}=\frac{1}{2}\sigma_{0}\tau_{0}{\varsigma}_{0}-\frac{1}{2}\sigma_{y}\tau_{y}{\varsigma}_{0},\qquad
O^{c}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}=\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}O^{c}_{a\eta
s,a^{\prime}\eta s},\qquad O^{c}_{a\eta s,a\eta s}=0,\qquad O^{cf}=0\ .$
(S309)

According to the mean field equations in Section S4.1, the one-shot HF mean
field and energies are

$\overline{H}_{U}=\frac{U_{1}}{2}\sum_{\mathbf{R}}\sum_{\alpha\eta\alpha^{\prime}\eta^{\prime}s}f_{\mathbf{R}\alpha\eta
s}^{\dagger}[\sigma_{y}]_{\alpha\alpha^{\prime}}[\tau_{y}]_{\eta\eta^{\prime}}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s},\qquad\overline{H}_{W}=\overline{H}_{V}=0\
,$ (S310) $\displaystyle\overline{H}_{J}=$
$\displaystyle-\frac{J}{2}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha_{1}\alpha_{2}\eta_{1}\eta_{2}s_{1}s_{2}}(\eta_{1}\eta_{2}+(-1)^{\alpha_{1}+\alpha_{2}})c_{\mathbf{k}\alpha_{2}+2,\eta_{2}s_{2}}^{\dagger}c_{\mathbf{k}\alpha_{1}+2,\eta_{1}s_{1}}\left(O^{f}_{\alpha_{1}\eta_{1}s_{1},\alpha_{2}\eta_{2}s_{2}}-\frac{1}{2}\delta_{\alpha_{1}\alpha_{2}}\delta_{\eta_{1}\eta_{2}}\delta_{s_{1}s_{2}}\right)$
$\displaystyle=\frac{J}{4}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha_{1}\alpha_{2}\eta_{1}\eta_{2}s}(\eta_{1}\eta_{2}+(-1)^{\alpha_{1}+\alpha_{2}})c_{\mathbf{k}\alpha_{2}+2,\eta_{2}s}^{\dagger}c_{\mathbf{k}\alpha_{1}+2,\eta_{1}s}[\sigma_{y}]_{\alpha_{1}\alpha_{2}}[\tau_{y}]_{\eta_{1}\eta_{2}}$
$\displaystyle=-\frac{J}{2}\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{\alpha_{1}\eta_{1}\alpha_{2}\eta_{2}s}[\sigma_{y}]_{\alpha_{1}\alpha_{2}}[\tau_{y}]_{\eta_{1}\eta_{2}}c_{\mathbf{k}\alpha_{2}+2\eta^{\prime}s}^{\dagger}c_{\mathbf{k}\alpha_{1}+2\eta
s}\ ,$ (S311)

and

$E_{U}=-2NU_{1},\qquad E_{W}=E_{V}=E_{J}=0\ .$ (S312)

Notice that the $\frac{1}{2}\sigma_{0}\tau_{0}{\varsigma}_{0}$ component in
$O^{f}$ does not contribute to the mean field Hamiltonian because they are
canceled by the constant bilinear terms in Eqs. S279 and S4.1, respectively.
The low energy band structure can be described by the k$\cdot$p expansion of
${\hat{H}}_{0}+\overline{H}_{U}+\overline{H}_{J}$.

$H_{\text{K-IVC}}^{\rm(MF)}(\mathbf{k})\approx\begin{pmatrix}0&v_{\star}(k_{x}\sigma_{0}\tau_{z}+ik_{y}\sigma_{z}\tau_{0})&\gamma\sigma_{0}\tau_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0})\\\
v_{\star}(k_{x}\sigma_{0}\tau_{z}-ik_{y}\sigma_{z}\tau_{0})&M\sigma_{x}\tau_{0}-\frac{J}{2}\sigma_{y}\tau_{y}&0\\\
\gamma\sigma_{0}\tau_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0})&0&\frac{U_{1}}{2}\sigma_{y}\tau_{y}\end{pmatrix}\otimes{\varsigma}_{0}\
.$ (S313)

For large $\mathbf{k}$ ($v_{\star}|\mathbf{k}|\gg\frac{1}{2}U_{1}$), the low
energy bands are contributed by the local orbitals, which have the energies
$\pm U_{1}/2$. We can see this in Fig. S6 at large $\mathbf{k}$. For
$\mathbf{k}=0$, the energy eigenvalues are given by

$\pm\sqrt{M^{2}+\frac{J^{2}}{4}},\qquad\pm\frac{U_{1}}{4}\pm\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\
,$ (S314)

where the levels
$\frac{U_{1}}{4}\pm\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}$ are
same as those in the VP state. In Fig. S6(b), (e), the 24 energy levels
closest to Fermi level from low to high are

$-\frac{U_{1}}{4}-\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold}),\qquad\frac{U_{1}}{4}-\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold}),\qquad-\sqrt{M^{2}+\frac{J^{2}}{4}};(\text{4-fold}),$
(S315)
$\sqrt{M^{2}+\frac{J^{2}}{4}}\;(\text{4-fold}),\qquad-\frac{U_{1}}{4}+\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold}),\qquad\frac{U_{1}}{4}+\sqrt{\left(\frac{U_{1}}{4}\right)^{2}+\gamma^{2}}\;(\text{4-fold})\
.$ (S316)

Notice that in Fig. S6(b), (e) every single band is spin degenerate.

###  S4.3 General rules for the ground states

When $M=0$, one can see that the energy levels at $\mathbf{k}=0$ of the VP
(Eq. S306) and K-IVC (Eq. S314) states at $\nu=0$ become the same. A symmetry
reason is that, when $M=0$, the flat-U(4) symmetry discussed in Section S2.7
is an exact symmetry of the considered Hamiltonian. Therefore, when $M=0$, the
VP and K-IVC states must have the same energy and band structure because they
are related by a flat-U(4) rotation. Nonzero $M$ not only makes a difference
in the band structure but also lifts the energy degeneracy between the two
states. Now we argue that nonzero $M$ stabilizes the K-IVC state. As shown in
Eq. S294, the total energy can be calculated as the interaction energy,
$E_{U}+E_{W}+E_{J}+E_{V}$, subtracted from the quasi-particle energy, i.e.,

$-E_{U}-E_{W}-E_{J}-E_{V}+\langle\Psi|{\hat{H}}_{0}+\overline{H}_{U}+\overline{H}_{W}+\overline{H}_{J}+\overline{H}_{V}|\Psi\rangle=-E_{U}-E_{W}-E_{J}-E_{V}+\sum_{\mathbf{k}}\sum_{n\in\mathrm{occ}}\mathcal{E}_{\mathbf{k}n}\
,$ (S317)

where $\mathcal{E}_{\mathbf{k}n}$ is the HF spectrum and $n$ indexes occupied
bands. Since the considered interaction Hamiltonian respects the flat-U(4)
symmetry, the VP and K-IVC states have the same interaction energy (Eqs. S304
and S312). Thus, the state having lower quasi-particle energy will have a
lower total energy. Because the K-IVC state opens a larger gap
($2\sqrt{M^{2}+J^{2}/4}$) than the VP state ($2|M-J/2|$), the K-IVC state has
the lower total energy.

We find that similar arguments apply to other fillings. Here we summarize two
rules for the (translationally invariant) ground states at the one-shot level.
First, as argued at the end of Section S2.7 and numerically confirmed at the
end of Section S4.2.1, the flat-U(4) symmetry is a better approximation (of
the Hamiltonian) than the chiral-U(4) symmetry. Heuristically, the wave
function tends to be symmetric under permutation of flavor indices, e.g., the
four components $\\{\alpha,\eta,s\\}$ fixed $\xi=\eta(-1)^{\alpha-1}$ which
labels the fundamental U(4) representations (Section S2.7), such that the real
space wave function has as many as nodes and hence the Coulomb interaction can
be saved (Hund’s rule). Thus, we summarize the first rule

  1. 1.

When possible, two $f$-electrons tend to occupy one flavor of the flat-U(4)
symmetry, e.g., $f_{\mathbf{R}1\eta s}^{\dagger}f_{\mathbf{R}2\eta
s}^{\dagger}$, such that the two $f$-electrons and their flat-U(4) rotations
form a maximal weight representation of the flat-U(4) group.

The state $f_{\mathbf{R}1\eta s}^{\dagger}f_{\mathbf{R}2\eta s}^{\dagger}$ and
its flat-U(4) rotations form a degenerate U(4) multiplet if $M=0$. The
multiplet form the irreducible representation $[2]_{4}$. (Here
$[\lambda_{1},\lambda_{2}\cdots]_{4}$ is the Young tableau notation for SU(4)
irreducible representations.) Readers may refer to Ref. Lian  _et al._ (2021)
for how the U(4) irreps are used to label the ground states of MATBG. We then
need to consider how nonzero $M$ will split the U(4) multiplet and select the
ground state. Because the considered interaction Hamiltonian respects the
flat-U(4) symmetry, states in the multiplet must have the same interaction
energy. Thus, it is the single-particle energy of the conduction bands that
determines the ground state - the phase opening the largest gap has the lowest
energy. However, since there are many levels, i.e., 8 from local orbitals and
16 from conduction bands, it is still not immediate to determine the gap at
generic fillings. We now assume that the $\Gamma_{3}$ states from the local
orbitals and the $\Gamma_{3}$ states ($a=1,2$) from the conduction bands are
all at high energies. This is a reasonable assumption based on the BM model
because previous studies Lian  _et al._ (2021); Kang and Vafek (2019);
Bultinck  _et al._ (2020a) has obtained the ground states by projecting the
Hamiltonian into the topological flat bands, which only has the
$\Gamma_{1}\oplus\Gamma_{2}$ states ($a=3,4$). We have also confirmed that the
lowest (closest to the gap) energy levels are indeed contributed by the
$\Gamma_{1}\oplus\Gamma_{2}$ ($a=3,4$) states for all the integer fillings, in
agreement with those in Refs. Bernevig  _et al._ (2021b); Bultinck  _et al._
(2020a); Vafek and Kang (2021); Cea and Guinea (2020); Zhang  _et al._ (2020);
Kang  _et al._ (2021). Thus we only need to look at the gap formed by the
$\Gamma_{1}\oplus\Gamma_{2}$ states ($a=3,4$). The effective Hamiltonian for
the $\Gamma_{1}\oplus\Gamma_{2}$ states ($a=3,4$) can be obtained by
restricting the $a$ index of the one-shot mean field Hamiltonian
${\hat{H}}_{0}+\overline{H}_{U}+\overline{H}_{W}+\overline{H}_{J}+\overline{H}_{V}$
(Sections S2.5, S279, S4.1, S4.1 and S292) into the subspace spanned by
$c_{\mathbf{k}a\eta s}$ ($a=3,4$). It reads

$\displaystyle H_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}^{(\Gamma_{1}\oplus\Gamma_{2})}(\mathbf{k}=0)=$
$\displaystyle\nu_{f}W_{3}\delta_{aa^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}+M[\sigma_{x}]_{a-2,a^{\prime}-2}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}-J\delta_{aa^{\prime}}\delta_{\eta\eta^{\prime}}\left(O_{a^{\prime}-2,\eta^{\prime},s^{\prime}\
;\ a-2,\eta,s}^{f}-\frac{1}{2}\delta_{ss^{\prime}}\right)$
$\displaystyle+J\delta_{\overline{a},a^{\prime}}\delta_{-\eta,\eta^{\prime}}O^{f}_{\overline{a}-2,-\eta,s^{\prime}\
;\ a-2,\eta,s}\ ,$ (S318)

where $a,a^{\prime}=3,4$. We can equivalently write it in a more compact form
as

$H^{(\Gamma_{1}\oplus\Gamma_{2})}(\mathbf{k}=0)=\nu_{f}W_{3}\sigma_{0}\tau_{0}{\varsigma}_{0}+M\sigma_{x}\tau_{0}{\varsigma}_{0}-\frac{J}{2}\tau_{z}(O^{fT}-\frac{1}{2}\sigma_{0}\tau_{0}{\varsigma}_{0})\tau_{z}-\frac{J}{2}\sigma_{z}(O^{fT}-\frac{1}{2}\sigma_{0}\tau_{0}{\varsigma}_{0})\sigma_{z}\
.$ (S319)

The second rule then follows

  1. 2.

For a set of flat-U(4) rotation related states, the state minimizing the
energy of Eq. S319 is the ground state.

The onsite repulsion $U_{1}$, and the couplings between local orbitals and
conduction bands, i.e., the $\gamma$ and $v_{\star}$ terms (Eq. S83), do not
enter the second rule because they only affect states at higher energies
contributed by the $f$-orbitals and the $\Gamma_{3}$ $c$-band basis but not
the low energy states contributed by the $\Gamma_{1}\oplus\Gamma_{2}$ $c$-band
basis, as exampled in Section S4.2.2.

In the next subsection, we will apply the two rules to the ground states at
other integer fillings.

###  S4.4 Correlated insulator phases at the fillings $\nu=-1,-2,-3$

Figure S7: HF band structures of correlated insulator phases at $\nu=-1$. (a),
(b), (c), (d) are the one-shot HF band structures of the VP, K-IVC, IVC, and
mixed VP + K-IVC phases, respectively. (e), (f), (g), (h) are the self-
consistent HF band structures of the VP, K-IVC, IVC, and mixed VP + K-IVC
phases, respectively. The color represents the composition of the energy
bands, where yellow corresponds to the local orbitals and blue corresponds to
the conduction bands. We have chosen $w_{0}/w_{1}=0.8$ in the calculation.
Other parameters of the single-particle and interaction hamiltonians are given
in Tables S4 and S6.

####  S4.4.1 Filling $\nu=-1$

We consider the parent wave function of (fully) valley polarized (VP) and
(partially) spin polarized state at $\nu=-1$ as

$|\mathrm{VP}_{0}^{\nu=-1}\rangle=\prod_{\mathbf{R}}f_{\mathbf{R}1+\uparrow}^{\dagger}f_{\mathbf{R}1+\downarrow}^{\dagger}f_{\mathbf{R}2+\uparrow}^{\dagger}|\mathrm{FS}\rangle\
,$ (S320)

where $|\rm FS\rangle$ is the Fermi sea state of the conduction bands at the
charge neutrality point (Eq. S296). The parent wave functions for the IVC and
K-IVC states are rotated from $|\mathrm{VP}_{0}^{\nu=-1}\rangle$ by the
chiral-U(4) (Eq. S92) and flat-U(4) (Eq. S105) symmetries as

$|\mathrm{IVC}_{0}^{\nu=-1}\rangle=e^{-i\frac{\pi}{2}{\hat{\Theta}}_{x0}}|\mathrm{VP}_{0}^{\nu=-1}\rangle=\prod_{\mathbf{R}}\frac{1}{2\sqrt{2}}(f_{\mathbf{R}1+\uparrow}^{\dagger}-if_{\mathbf{R}2-\uparrow}^{\dagger})(f_{\mathbf{R}1+\downarrow}^{\dagger}-if_{\mathbf{R}2-\downarrow}^{\dagger})(-if_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})|\rm
FS\rangle\ ,$ (S321)

and

$|\text{K-IVC}_{0}^{\nu=-1}\rangle=e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}|\mathrm{VP}_{0}^{\nu=-1}\rangle=\prod_{\mathbf{R}}\frac{1}{2\sqrt{2}}(f_{\mathbf{R}1+\uparrow}^{\dagger}+f_{\mathbf{R}2-\uparrow}^{\dagger})(f_{\mathbf{R}1+\downarrow}^{\dagger}+f_{\mathbf{R}2-\downarrow}^{\dagger})(-f_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})|\rm
FS\rangle\ ,$ (S322)

respectively. Ref. Lian  _et al._ (2021) found that a mixed state of VP and
K-IVC has the lowest energy at $\nu=-1$. The parent wave function of the mixed
state can be written as

$|\text{VP+K-IVC}_{0}^{\nu=-1}\rangle=\prod_{\mathbf{R}}f_{\mathbf{R}1+\downarrow}^{\dagger}e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}f_{\mathbf{R}1+\uparrow}^{\dagger}f_{\mathbf{R}2+\uparrow}^{\dagger}e^{i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}|\mathrm{FS}\rangle=\prod_{\mathbf{R}}\frac{1}{2}f_{\mathbf{R}1+\downarrow}^{\dagger}(f_{\mathbf{R}1+\uparrow}^{\dagger}+f_{\mathbf{R}2-\uparrow}^{\dagger})(-f_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})|\rm
FS\rangle\ .$ (S323)

According to Ref. Lian  _et al._ (2021), the two states
($f_{\mathbf{R}1+\uparrow}^{\dagger}f_{\mathbf{R}2+\uparrow}^{\dagger}$) from
the same valley-spin flavor are rotated to a K-IVC state, whereas the state
$f_{\mathbf{R}1+\downarrow}^{\dagger}$ remains valley and spin polarized.
Feeding the parent states into the HF loop, we obtain the one-shot and self-
consistent HF band structures (Fig. S7) and energies (Table S7). We find that
the mixed state indeed has the lowest energy.

All the one-shot states, i.e., $|X^{\nu=-1}_{1}\rangle$ ($X=$VP, IVC, K-IVC,
VP+K-IVC), and self-consistent states, i.e., $|X^{\nu=-1}_{\infty}\rangle$,
are gapped Fock states occupying the HF bands (Fig. S7) up to the filling
$\nu=-1$. As will be explained in Section S4.5, all these states have the
Chern number $C=1$. It is worth mentioning that the parent states
$|X^{\nu=-1}_{0}\rangle$, even though lead to gapped one-shot states with
nonzero Chern numbers, do not have Chern numbers because they occupy gapless
conduction bands plus a select number of local orbitals.

We now apply the two rules given in Section S4.3 to discuss the one-shot
ground states $|X^{\nu=-1}_{1}\rangle$ ($X=$VP, IVC, K-IVC, VP+K-IVC). VP,
K-IVC and VP+K-IVC states satisfy the first rule, i.e., two $f$-electrons tend
to occupy a flavor flat-U(4) symmetry. Here the two $f$-electrons are
$f_{\mathbf{R}1+\uparrow}^{\dagger}f_{\mathbf{R}2+\uparrow}^{\dagger}$ or its
flat-U(4) partners. We need to compare the quasi-particle energies of the
three states to determine the ground states. The density matrices given by
$|\text{VP}^{\nu=-1}_{0}\rangle$ are

$O^{f}=\sigma_{0}(\frac{\tau_{z}+\tau_{0}}{2})(\frac{{\varsigma}_{0}+{\varsigma}_{z}}{2})+(\frac{\sigma_{0}+\sigma_{z}}{2})(\frac{\tau_{0}+\tau_{z}}{2})(\frac{{\varsigma}_{0}-{\varsigma}_{z}}{2}),\quad
O^{c}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}=\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}O^{c}_{a\eta
s,a^{\prime}\eta s},\quad O^{c}_{a\eta s,a\eta s}=0,\quad O^{cf}=0\ .$ (S324)

The density matrix $O^{f}$ is diagonal in spin. The two spin blocks of $O^{f}$
are

$O^{f(s=\uparrow)}=\sigma_{0}(\frac{\tau_{z}+\tau_{0}}{2}),\qquad
O^{f(s=\downarrow)}=(\frac{\sigma_{0}+\sigma_{z}}{2})(\frac{\tau_{0}+\tau_{z}}{2})\
.$ (S325)

Substituting $O^{f(s)}$ into Eq. S319, we obtain the two spin blocks of
$H^{(\Gamma_{1}\oplus\Gamma_{2})}(\mathbf{k}=0)$ as

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\uparrow)}_{\rm
VP}=-W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}-\frac{J}{2}\tau_{z}\frac{\sigma_{0}\tau_{z}}{2}\tau_{z}-\frac{J}{2}\sigma_{z}\frac{\sigma_{0}\tau_{z}}{2}\sigma_{z}=-W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}-\frac{J}{2}\sigma_{0}\tau_{z}\
,$ (S326)

and

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\downarrow)}_{\rm
VP}=-W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}-\frac{J}{2}\tau_{z}X\tau_{z}-\frac{J}{2}\sigma_{z}X\sigma_{z}=-W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}-\frac{J}{2}\mathrm{diag}([1,-1,-1,-1])\
,$ (S327)

respectively, where
$X=(\frac{\sigma_{0}+\sigma_{z}}{2})(\frac{\tau_{0}+\tau_{z}}{2})-\frac{1}{2}\sigma_{0}\tau_{0}=\frac{1}{2}\mathrm{diag}([1,-1,-1,-1])$.
Here the subscript “VP” indicates that the corresponding density matrices are
generated by the parent VP state. The lowest three levels at $\mathbf{k}=0$ of
the VP state are

$-W_{3}-\sqrt{M^{2}+\frac{J^{2}}{4}},\qquad-W_{3}-M-\frac{J}{2},\qquad-
W_{3}+M-\frac{J}{2}\ ,$ (S328)

which come from the $s=\downarrow,\eta=+$, $s=\uparrow,\eta=+$, and
$s=\uparrow,\eta=+$ flavors, respectively. (We have made use of $J/2>M$.)

Applying a flat-U(4) rotation $e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}$ to the
order parameter, we obtain the order parameters of the K-IVC state as

$O^{f(s=\uparrow)}=\frac{1}{2}\sigma_{0}\tau_{0}-\frac{1}{2}\sigma_{y}\tau_{y},\qquad
O^{f(s=\downarrow)}=\begin{pmatrix}0&0&0&1\\\ 0&-1&0&0\\\ 0&0&-1&0\\\
1&0&0&0\end{pmatrix}\ .$ (S329)

Substituting $O^{f(s)}$ into Eq. S319, we obtain the two spin blocks of
$H^{(\Gamma_{1}\oplus\Gamma_{2})}(\mathbf{k}=0)$ as

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\uparrow)}_{\rm
K-IVC}=-W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}-\frac{J}{2}\sigma_{y}\tau_{y},\qquad
H^{(\Gamma_{1}\oplus\Gamma_{2},s=\downarrow)}_{\rm
K-IVC}=-W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}+\frac{J}{2}\begin{pmatrix}0&0&0&1\\\
0&1&0&0\\\ 0&0&1&0\\\ 1&0&0&0\end{pmatrix}\ ,$ (S330)

where the subscript “K-IVC” indicates that the corresponding density matrices
are generated by the parent K-IVC state. The lowest three levels at
$\mathbf{k}=0$ of the K-IVC state are all
$-W_{3}-\sqrt{M^{2}+\frac{J^{2}}{4}}$, where two come from the $s=\uparrow$
flavor and one comes from the $s=\downarrow$ flavor. Clearly, the K-IVC state
has lower quasi-particle energy than the VP state. In the end we consider the
one-shot HF mean field Hamiltonian spanned by the $\Gamma_{1}\oplus\Gamma_{2}$
for the VP+K-IVC state

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\uparrow)}_{\rm
VP+K-IVC}=-W_{3}+M\sigma_{x}\tau_{0}-\frac{J}{2}\sigma_{y}\tau_{y},\qquad
H^{(\Gamma_{1}\oplus\Gamma_{2},s=\downarrow)}_{\rm
VP+K-IVC}=-W_{3}+M\sigma_{x}\tau_{0}-\frac{J}{2}\mathrm{diag}([1,-1,-1,-1])\
,$ (S331)

where the subscript “VP+K-IVC” indicates that the corresponding density
matrices are generated by the parent mixed (VP+K-IVC) state. The lowest three
levels at $\mathbf{k}=0$ of the VP+K-IVC state are also all
$-W_{3}-\sqrt{M^{2}+\frac{J^{2}}{4}}$, which are same as those of the K-IVC
state. Thus the rules in Section S4.3 cannot tell us whether K-IVC or VP+K-IVC
has a lower energy. In fact, the one-shot and self-consistent HF energies of
the two states are indeed very close to each other (Table S7).

Taking into account the bands at $\mathbf{k}=0$ contributed by the
$\Gamma_{3}$ basis still cannot distinguish the energies of the K-IVC and the
VP+K-IVC states. Following a similar calculation as in Section S4.2.2, one
will find that the energy levels from the $\Gamma_{3}$ basis are the same for
both states. The energy difference must be contributed by the bands at finite
$\mathbf{k}$.

As shown in Fig. S7, one-shot and self-consistent HF bands of all the
different states at $\nu=-1$ have a feature in common: There is a set of flat
bands above the zero energy. Such flat bands are also observed in one of our
previous studies Bernevig  _et al._ (2021b) but have not been understood. (See
the particle excitation spectra in Fig. 11 of Ref. Bernevig  _et al._
(2021b).) We now can explain the origin of the flatness through our
topological heavy fermion model. We consider the k$\cdot$p expansion of the
one-shot mean-field Hamiltonian. For simplicity, here we mainly focus on the
VP state. Following the same procedure we have done to obtain the one-shot
mean field Hamiltonian at $\nu=0$ (Eq. S305), we obtain

$H_{\rm
VP}^{\rm(MF)}(\mathbf{k})\approx\begin{pmatrix}-W_{1}\sigma_{0}\tau_{0}{\varsigma}_{0}&v_{\star}(k_{x}\sigma_{0}\tau_{z}+ik_{y}\sigma_{z}\tau_{0}){\varsigma}_{0}&\gamma\sigma_{0}\tau_{0}{\varsigma}_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0}){\varsigma}_{0}\\\
v_{\star}(k_{x}\sigma_{0}\tau_{z}-ik_{y}\sigma_{z}\tau_{0}){\varsigma}_{0}&-W_{3}\sigma_{0}\tau_{0}{\varsigma}_{0}+M\sigma_{x}\tau_{0}{\varsigma}_{0}-\frac{J}{2}\sigma_{0}\tau_{z}{\varsigma}_{0}&0\\\
\gamma\sigma_{0}\tau_{0}{\varsigma}_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0}){\varsigma}_{0}&0&-(U_{1}+6U_{2})\sigma_{0}\tau_{0}{\varsigma}_{0}-U_{1}(O^{f}-\frac{1}{2}\sigma_{0}\tau_{0}{\varsigma}_{0})\end{pmatrix}.$
(S332)

The density matrix $O^{f}$ is given in Eq. S324. Comparing it to the one-shot
mean field Hamiltonian at $\nu=0$ (Eq. S305), there are three additional
terms, i.e., $-W_{1}\sigma_{0}\tau_{0}{\varsigma}_{0}$,
$-W_{3}\sigma_{0}\tau_{0}{\varsigma}_{0}$, and
$-(U_{1}+6U_{2})\sigma_{0}\tau_{0}{\varsigma}_{0}$ in the three diagonal
blocks. These three terms come from the Hartree channel terms $\nu_{f}W_{1}$,
$\nu_{f}W_{3}$, $\nu_{f}(U_{1}+6U_{2})$ of $\overline{H}_{W}$ (Section S4.1)
and $\overline{H}_{U}$ (Eq. S279) and only shift the energies of the three
blocks. Without the energy shift and the couplings
($\gamma,v_{\star}^{\prime}$) between $c$\- and $f$-electrobs, the $c$-bands
(the first two blocks) would have a quadratic touching at zero energy and the
$f$-levels (the third block) have energies $\pm U_{1}/2$, as illustrated by
the red bands in Fig. 3(a) in the main text. Using the parameters obtained at
$w_{0}/w_{1}=0.8$ (Table S6), the average energy shift of the $c$-electron
bands (the first two blocks) is $-(W_{1}+W_{3})/2\approx-47$meV, and the
energy shift of the $f$-electron bands is $-(U_{1}+6U_{2})\approx-72$meV.
Thus, the relative energy shift of $f$-electrons with respect to $c$-electrons
is given by $\delta E\approx(-72+47)\mathrm{meV}\approx-25$meV and is
approximately $-U_{1}/2\approx 29$meV. That means, if we turn off the
hybridizations, the upper branch of the $f$-electron levels will be shifted to
the quadratic touching point of the $c$-electrons, as shown by the red bands
in Fig. 3(b) in the main text. Turning on the hybridizations will gap out the
quadratic touching point, then the upper branch of $f$-electron bands form an
isolated set of flat bands.

####  S4.4.2 Filling $\nu=-2$

Figure S8: HF band structures of correlated insulator phases at $\nu=-2$. (a),
(b), (c) are the one-shot HF band structures of the VP, K-IVC, and IVC phases,
respectively. (d), (e), (f) are the self-consistent HF band structures of the
VP, K-IVC, and IVC phases, respectively. The color represents the composition
of the energy bands, where yellow corresponds to the local orbitals and blue
corresponds to the conduction bands. We have chosen $w_{0}/w_{1}=0.8$ in the
calculation. Other parameters of the single-particle and interaction
hamiltonians are given in Tables S4 and S6.

We consider the parent wave function of valley polarized (VP) state at
$\nu=-2$ as

$|\mathrm{VP}_{0}^{\nu=-2}\rangle=\prod_{\mathbf{R}}f_{\mathbf{R}1+\uparrow}^{\dagger}f_{\mathbf{R}2+\uparrow}^{\dagger}|\mathrm{FS}\rangle\
,$ (S333)

where $|\rm FS\rangle$ is the Fermi sea state of the conduction bands at the
charge neutrality point (Eq. S296). The parent wave functions for the IVC and
K-IVC states are rotated from $|\mathrm{VP}_{0}^{\nu=-1}\rangle$ by the
chiral-U(4) (Eq. S92) and flat-U(4) (Eq. S105) symmetries as

$|\mathrm{IVC}_{0}^{\nu=-2}\rangle=e^{-i\frac{\pi}{2}{\hat{\Theta}}_{x0}}|\mathrm{VP}_{0}^{\nu=-1}\rangle=\prod_{\mathbf{R}}\frac{1}{2}(f_{\mathbf{R}1+\uparrow}^{\dagger}-if_{\mathbf{R}2-\uparrow}^{\dagger})(-if_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})|\rm
FS\rangle\ ,$ (S334)

and

$|\text{K-IVC}_{0}^{\nu=-2}\rangle=e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}|\mathrm{VP}_{0}^{\nu=-1}\rangle=\prod_{\mathbf{R}}\frac{1}{2}(f_{\mathbf{R}1+\uparrow}^{\dagger}+f_{\mathbf{R}2-\uparrow}^{\dagger})(-f_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})|\rm
FS\rangle\ ,$ (S335)

respectively. Feeding the parent states into the HF loops, we obtain the HF
band structures and energies as shown in Fig. S8 and Table S7, respectively.
We find that the K-IVC state has the lowest energy.

We now apply the two rules given in Section S4.3 to discuss the one-shot
states, i.e., $|X^{\nu=-2}_{1}\rangle$ ($X=$VP, IVC, K-IVC), which are Fock
states occupying the one-shot HF bands (Fig. S8) to the filling $\nu=-2$. The
VP and K-IVC states satisfy the first rule, i.e., two $f$-electrons tend to
occupy a flavor flat-U(4) symmetry. Here the two $f$-electrons are
$f_{\mathbf{R}1+\uparrow}^{\dagger}f_{\mathbf{R}2+\uparrow}^{\dagger}$ or its
flat-U(4) partners. We need to compare the quasi-particle energies of the
three states to determine the ground states. The density matrices given by
$|\text{VP}^{\nu=-2}_{0}\rangle$ are

$O^{f}=\sigma_{0}(\frac{\tau_{z}+\tau_{0}}{2})(\frac{{\varsigma}_{0}+{\varsigma}_{z}}{2}),\qquad
O^{c}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}=\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}O^{c}_{a\eta
s,a^{\prime}\eta s},\qquad O^{c}_{a\eta s,a\eta s}=0,\qquad O^{cf}=0\ .$
(S336)

According to Eq. S319 and the density matrices, the one-shot HF mean field
Hamiltonian spanned by the $\Gamma_{1}\oplus\Gamma_{2}$ states at
$\mathbf{k}=0$, i.e., $H^{(\Gamma_{1}\oplus\Gamma_{2})}$, is diagonal in the
spin index. The spin-up and down blocks are given by

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\uparrow)}_{\rm
VP}=-2W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}-\frac{J}{2}\sigma_{0}\tau_{z},\qquad
H^{(\Gamma_{1}\oplus\Gamma_{2},s=\downarrow)}_{\rm
VP}=-2W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}+\frac{J}{2}\sigma_{0}\tau_{0}\
,$ (S337)

where the subscript “VP” indicates that the corresponding density matrices are
generated by parent VP state. The lowest two levels at $\mathbf{k}=0$ of the
VP state are

$-2W_{3}-M-\frac{J}{2},\qquad-2W_{3}+M-\frac{J}{2}\ ,$ (S338)

both of which come from the $s=\uparrow$ flavor. (We have making use of
$J/2>M$.) Applying a flat-U(4) rotation
$e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}$ to the order parameter, we obtain the
one-shot HF mean field Hamiltonian spanned by the $\Gamma_{1}\oplus\Gamma_{2}$
(Eq. S319) for the K-IVC state

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\uparrow)}_{\rm
K-IVC}=-2W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}-\frac{J}{2}\sigma_{y}\tau_{y},\qquad
H^{(\Gamma_{1}\oplus\Gamma_{2},s=\downarrow)}_{\rm
K-IVC}=-2W_{3}\sigma_{0}\tau_{0}+M\sigma_{x}\tau_{0}+\frac{J}{2}\sigma_{0}\tau_{0}\
,$ (S339)

where the subscript “K-IVC” indicates that the corresponding density matrices
are generated by parent VP state. The lowest two levels at $\mathbf{k}=0$ of
the K-IVC state are

$-2W_{3}-\sqrt{M^{2}+\frac{J^{2}}{4}},\qquad-2W_{3}-\sqrt{M^{2}+\frac{J^{2}}{4}}\
,$ (S340)

both of which come from the $s=\uparrow$ flavor. Clearly, the K-IVC state,
where the $M\sigma_{x}\tau_{0}$ and $\frac{J}{2}\sigma_{y}\tau_{y}$ terms in
the spin-up sector anti-commute with each other, has lower quasi-particle
energy than the VP state, where the $M\sigma_{x}\tau_{0}$ and
$\frac{J}{2}\sigma_{0}\tau_{z}$ terms in the spin-up sector commute with each
other.

####  S4.4.3 Filling $\nu=-3$

Figure S9: HF band structures of correlated insulator phases at $\nu=-3$. (a),
(b) are the one-shot HF band structures of the VP and K-IVC phases,
respectively. (c), (d) are the self-consistent HF band structures of the VP
and K-IVC phases, respectively. The color represents the composition of the
energy bands, where yellow corresponds to the local orbitals and blue
corresponds to the conduction bands. We have chosen $w_{0}/w_{1}=0.8$ in the
calculation. Other parameters of the single-particle and interaction
hamiltonians are given in Tables S4 and S6.

We consider the parent wave function of valley polarized (VP) and spin
polarized state at $\nu=-3$ as

$|\mathrm{VP}_{0}^{\nu=-3}\rangle=\prod_{\mathbf{R}}f_{\mathbf{R}1+\uparrow}^{\dagger}|\mathrm{FS}\rangle\
,$ (S341)

where $|\rm FS\rangle$ is the Fermi sea state of the conduction bands at the
charge neutrality point. The parent wave functions for the IVC and K-IVC
states are rotated from $|\mathrm{VP}_{0}^{\nu=-3}\rangle$ by the chiral-U(4)
(Eq. S92) and flat-U(4) (Eq. S105) symmetries as

$|\mathrm{IVC}_{0}^{\nu=-3}\rangle=e^{-i\frac{\pi}{2}{\hat{\Theta}}_{x0}}|\mathrm{VP}_{0}^{\nu=-3}\rangle=\prod_{\mathbf{R}}\frac{1}{\sqrt{2}}(f_{\mathbf{R}1+\uparrow}^{\dagger}-if_{\mathbf{R}2-\uparrow}^{\dagger})|\rm
FS\rangle\ ,$ (S342)

and

$|\text{K-IVC}_{0}^{\nu=-3}\rangle=e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}|\mathrm{VP}_{0}^{\nu=-3}\rangle=\prod_{\mathbf{R}}\frac{1}{\sqrt{2}}(f_{\mathbf{R}1+\uparrow}^{\dagger}+f_{\mathbf{R}2-\uparrow}^{\dagger})|\rm
FS\rangle\ ,$ (S343)

respectively. We can see that $|\mathrm{IVC}_{0}^{\nu=-3}\rangle$ and
$|\text{K-IVC}_{0}^{\nu=-3}\rangle$ are related by a valley-U(1) rotation,
which adds the factor $i$ to the $f$-electron
($f_{\mathbf{R}2-\uparrow}^{\dagger}$) in the valley $\eta=-$. (One should
notice that, at the filling $\nu=-2$, the K-IVC (Eq. S335) and IVC (Eq. S334)
parent states are not related by a valley-U(1) rotation. If one applies the
same rotation as above, one will obtain
$\prod_{\mathbf{R}}\frac{1}{2}(f_{\mathbf{R}1+\uparrow}^{\dagger}+f_{\mathbf{R}2-\uparrow}^{\dagger})(f_{\mathbf{R}1-\uparrow}^{\dagger}+f_{\mathbf{R}2+\uparrow}^{\dagger})|\text{FS}\rangle$,
which is still different from the K-IVC parent state (Eq. S335).) Thus the two
states must have the same band structure and the same energy. Hence in the
following we will only discuss the VP and the K-IVC states. Feeding the parent
states into the HF loops, we obtain the HF band structures and energies as
shown in Fig. S9 and Table S7, respectively. We find that the VP state has the
lowest energy.

All the one-shot states, i.e., $|X^{\nu=-3}_{1}\rangle$ ($X=$VP, K-IVC), and
self-consistent states, i.e., $|X^{\nu=-3}_{\infty}\rangle$, are Fock states
occupying the HF bands (Fig. S7) up to the filling $\nu=-3$. Even though they
do not have indirect gaps, they do have direct gaps at every momentum. Hence
Chern numbers are still well defined for the direct gaps. As will be explained
in Section S4.5, all these states have the Chern number $C=1$. It is worth
mentioning that the parent states $|X^{\nu=-3}_{1}\rangle$, even though lead
to one-shot states with nonzero Chern numbers, do not have Chern numbers
because they occupy gapless conduction bands plus one local orbital.

We now apply the two rules given in Section S4.3 to discuss the ground states.
We need to compare the quasi-particle energies of the two states to determine
the ground states. The density matrices given by
$|\text{VP}^{\nu=-3}_{0}\rangle$ are

$O^{f}=(\frac{\sigma_{0}+\sigma_{z}}{2})(\frac{\tau_{z}+\tau_{0}}{2})(\frac{{\varsigma}_{0}+{\varsigma}_{z}}{2}),\qquad
O^{c}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}=\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}O^{c}_{a\eta
s,a^{\prime}\eta s},\qquad O^{c}_{a\eta s,a\eta s}=0,\qquad O^{cf}=0\ .$
(S344)

According to Eq. S319 and the density matrices, the one-shot HF mean field
Hamiltonian spanned by the $\Gamma_{1}\oplus\Gamma_{2}$ states at
$\mathbf{k}=0$, i.e., $H^{(\Gamma_{1}\oplus\Gamma_{2})}$, is diagonal in the
spin index. The spin-up and down blocks of $\Gamma_{1}\oplus\Gamma_{2}$ states
are given by

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\uparrow)}_{\rm
VP}=-3W_{3}+M\tau_{0}\sigma_{x}-\frac{J}{2}\mathrm{diag}([1,-1,-1,-1]),\qquad
H^{(\Gamma_{1}\oplus\Gamma_{2},s=\downarrow)}_{\rm
VP}=-3W_{3}+M\sigma_{x}+\frac{J}{2}\sigma_{0}\ ,$ (S345)

where the subscript “VP” indicates that the input density matrices are
generated by the parent VP state. The lowest level at $\mathbf{k}=0$ of the VP
state is

$-3W_{3}-\sqrt{M^{2}+\frac{J^{2}}{4}},$ (S346)

which comes from the $s=\uparrow$ flavor. Applying a flat-U(4) rotation
$e^{-i\frac{\pi}{2}{\hat{\Sigma}}_{x0}}$ to the order parameter, we obtain the
one-shot HF mean field Hamiltonian spanned by the $\Gamma_{1}\oplus\Gamma_{2}$
(Eq. S319) for the K-IVC state

$H^{(\Gamma_{1}\oplus\Gamma_{2},s=\uparrow)}_{\rm
K-IVC}=-2W_{3}+M\tau_{0}\sigma_{x}+\frac{J}{2}\begin{pmatrix}0&0&0&1\\\
0&1&0&0\\\ 0&0&1&0\\\ 1&0&0&0\end{pmatrix},\qquad
H^{(\Gamma_{1}\oplus\Gamma_{2},s=\downarrow)}_{\rm
K-IVC}=-2W_{3}+M\sigma_{x}+\frac{J}{2}\sigma_{0}\ ,$ (S347)

where the subscript “K-IVC” indicates that the input density matrices are
generated by the parent K-IVC state. The lowest level at $\mathbf{k}=0$ of the
K-IVC state is still $-3W_{3}-\sqrt{M^{2}+\frac{J^{2}}{4}}$. Thus the second
rule given in Section S4.3 cannot tell us which has lower energy. In fact, the
one-shot and self-consistent HF energies of the two states are indeed very
close to each other (Table S7). As discussed in the end of Section S4.4.1, the
energy difference between the two states must be contributed by bands at
finite $\mathbf{k}$.

###  S4.5 Chern numbers of the ground states

We first consider the Chern numbers of the insulator ground states at the
filling $\nu=0$. Since the Chern numbers are stable against perturbations, we
can simplify the problem by (i) looking at the one-shot mean-field Hamiltonian
instead of the self-consistent mean-field Hamiltonian and (ii) deforming the
one-shot mean-field Hamiltonian adiabatically. As shown in previous
subsections, flat-U(4) is a good approximation of the Hamiltonian. Tuning $M$
to zero, which makes flat-U(4) an exact symmetry of the considered
Hamiltonian, will not close the gap of the one-shot VP ($2|M-J/2|$), gap of
one-shot K-IVC ($2\sqrt{M^{2}+J^{2}/4}$), and gaps of their U(4) partner
states, provided that $M<J/2$. Therefore, the VP and K-IVC states must have
the same Chern number since they are related by flat-U(4) rotation. We rewrite
the kp expansion of the one-shot mean-field Hamiltonian for the VP state (Eq.
S305) as

$H^{({\rm MF})}_{\rm
VP}(\mathbf{k})=\begin{pmatrix}A(\mathbf{k})&S(\mathbf{k})\\\
S^{\dagger}(\mathbf{k})&B\end{pmatrix}\otimes{\varsigma}_{0}$ (S348)

with

$A(\mathbf{k})=\begin{pmatrix}0&v_{\star}(k_{x}\tau_{z}\sigma_{0}+ik_{y}\tau_{0}\sigma_{z})\\\
v_{\star}(k_{x}\tau_{z}\sigma_{0}-ik_{y}\tau_{0}\sigma_{z})&M\tau_{0}\sigma_{x}-\frac{J}{2}\tau_{z}\sigma_{0}\end{pmatrix},\qquad
S(\mathbf{k})=\begin{pmatrix}\gamma\sigma_{0}\tau_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0})\\\
0\end{pmatrix}\ ,$ (S349)

$B=-\frac{U_{1}}{2}\sigma_{0}\tau_{z}$. Since $U_{1}$ is large, we can
integrate out the $f$-orbital levels to derive a four-by-four Hamiltonian to
describe the low energy bands around $\mathbf{k}=0$. Applying a second order
perturbation theory, or a Schrieffer-Wolf transformation that decouples $f$\-
and $c$-electrons, we obtain

$H^{({\rm eff})}_{\rm VP}(\mathbf{k})\approx
A(\mathbf{k})+S(\mathbf{k})B^{-1}S^{\dagger}(\mathbf{k})=\begin{pmatrix}\frac{2\gamma^{2}}{U_{1}}\tau_{z}\sigma_{0}+\frac{4\gamma
v_{\star}^{\prime}}{U_{1}}(k_{x}\tau_{0}\sigma_{x}+k_{y}\tau_{z}\sigma_{y})&v_{\star}(k_{x}\tau_{z}\sigma_{0}+ik_{y}\tau_{0}\sigma_{z})\\\
v_{\star}(k_{x}\tau_{z}\sigma_{0}-ik_{y}\tau_{0}\sigma_{z})&M\tau_{0}\sigma_{x}-\frac{J}{2}\tau_{z}\sigma_{0}\end{pmatrix}\otimes{\varsigma}_{0}\
.$ (S350)

Since the “integrated out” $f$-orbitals are topologically trivial, there
should be no worry about missing topological bands at higher energies. Without
closing the gap, we can continuously change the parameters as
$2\gamma^{2}/U\to m$, $v_{\star}^{\prime}\to 0$, $M\to 0$, $J/2\to-m$ with
$m>0$ such that the Hamiltonian becomes

$m\zeta_{z}\tau_{z}\sigma_{0}+v_{\star}k_{x}\zeta_{x}\tau_{z}\sigma_{0}-v_{\star}k_{y}\zeta_{y}\tau_{0}\sigma_{z}\
.$ (S351)

Here $\zeta_{x,y,z}$ are Pauli matrices introduced for the $\Gamma_{3}$
($\zeta_{z}=1$) and the $\Gamma_{1}\oplus\Gamma_{2}$ ($\zeta_{z}=-1$)
subspaces. We can view this Hamiltonian as four independent gapped Dirac
points in the eigenspaces of $\tau_{z}$ and $\sigma_{z}$. To be concrete, we
use $\eta=\pm 1$ and $(-1)^{\alpha-1}$ ($\alpha=1,2$) to represent the
eigenvalues of $\tau_{z}$ and $\sigma_{z}$, respectively. Then the Dirac
Hamiltonian can be written as

$m\eta\zeta_{z}+v_{\star}k_{x}\eta\zeta_{x}-v_{\star}k_{y}(-1)^{\alpha-1}\zeta_{y}.$
(S352)

The Chern number contributed by each Dirac point is given by
$-\frac{1}{2}\eta\cdot\eta\cdot(-1)^{\alpha-1}=-\frac{1}{2}(-1)^{\alpha-1}$.
The total Chern number is hence zero.

We then consider the Chern numbers of the insulator states at other fillings.
We can simplify the problem by deforming the one-shot mean-field Hamiltonian
adiabatically in the sense that the direct gap at each $\mathbf{k}$ does not
close. We first let $M\to 0$ to recover the flat-U(4) symmetry. The states do
not close their gaps ($2|J/2-M|$ or $\sqrt{M^{2}+J^{2}/4}$) in this process
because $M$ is small compared to $J/2$. Thus all the phases related by
flat-U(4) rotations must have the same topology. For a set of flat-U(4)
related states, we only need to calculate the Chern number of the simplest one
of them, i.e., the valley, spin, and orbital polarized state. We hence
consider the parent wave function

$|{\rm trial}\rangle=\prod_{\mathbf{R}}\prod_{\alpha\eta
s}(f_{\mathbf{R}\alpha\eta s}^{\dagger})^{n_{\alpha\eta
s}}|\mathrm{FS}\rangle\ ,$ (S353)

where $n_{\alpha\eta s}=1$ or 0 are the occupation numbers of the $\alpha$-th
($\alpha=1,2$) $f$-orbitals in the valley $\eta$ ($=\pm$) and spin $s$
($=\uparrow\downarrow$) sector. The density matrix and filling of local
orbitals are given by

$O_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}^{f}=\delta_{\alpha,\alpha^{\prime}}\delta_{\eta,\eta^{\prime}}\delta_{s,s^{\prime}}n_{\alpha\eta
s},\qquad\nu_{f}=-4+\sum_{\alpha\eta s}n_{\alpha\eta s},\qquad\nu_{c}=0\ .$
(S354)

The k$\cdot$p expansion of the one-shot mean-field Hamiltonian (from Sections
S2.5, S279, S4.1, S4.1 and S292) for the parent state reads

$H^{\rm(MF)}(\mathbf{k})\approx\begin{pmatrix}\nu_{f}W_{1}&v_{\star}(k_{x}\sigma_{0}\tau_{z}+ik_{y}\sigma_{z}\tau_{0}){\varsigma}_{0}&\gamma\sigma_{0}\tau_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0}){\varsigma}_{0}\\\
v_{\star}(k_{x}\sigma_{0}\tau_{z}-ik_{y}\sigma_{z}\tau_{0}){\varsigma}_{0}&\nu_{f}W_{3}+M\sigma_{x}\tau_{0}{\varsigma}_{0}-J(O^{f}-\frac{1}{2}\sigma_{0}\tau_{0}{\varsigma}_{0})&0\\\
\gamma\sigma_{0}\tau_{0}{\varsigma}_{0}+v_{\star}^{\prime}(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}\tau_{0}){\varsigma}_{0}&0&\nu_{f}(U_{1}+6U_{2})-U_{1}(O^{f}-\frac{1}{2}\sigma_{0}\tau_{0}{\varsigma}_{0})\end{pmatrix}\
.$ (S355)

We further deform the Hamiltonian $H^{\rm(MF)}(\mathbf{k})$ adiabatically to
simplify the problem. First, we can continuously turn off the $\nu_{f}W_{1}$,
$\nu_{f}W_{3}$, $\nu_{f}W_{1}(U_{1}+6U_{2})$ terms without changing the direct
gap at $\mathbf{k}=0$. A simple argument to justify this process is the
following. As mentioned in Section S4.3, for all the fillings the first block
($\Gamma_{3}$ representation from conduction bands) and the third block
($\Gamma_{3}$ representation from local orbitals) are at high energies and the
levels closest to the gap are always from the second block
($\Gamma_{1}\oplus\Gamma_{2}$ representation from the conduction bands), i.e.,
$\nu_{f}W_{3}\pm M\pm J/2$. In the process of turning off the three terms, the
first block and the third block continue to be at high energies and do not
enter the low energy physics around $\mathbf{k}=0$. This is confirmed by the
calculated levels in in the following (Eqs. S357 and S358). Second, we tune
$v_{\star}^{\prime}\to 0$, $M\to 0$, which is also adiabatic. The deformed
Hamiltonian is diagonal in the sublattice ($\sigma_{z}$), valley ($\tau_{z}$),
and spin (${\varsigma}_{z}$) indices, i.e.,

$H^{(\mathrm{MF},\alpha,\eta,s)}(\mathbf{k})\approx\begin{pmatrix}0&v_{\star}(\eta
k_{x}+i(-1)^{\alpha-1}k_{y})&\gamma\\\ v_{\star}(\eta
k_{x}-i(-1)^{\alpha-1}k_{y})&-J(n_{\alpha\eta s}-\frac{1}{2})&0\\\
\gamma&0&-U_{1}(n_{\alpha\eta s}-\frac{1}{2})\end{pmatrix}\ .$ (S356)

For $n_{\alpha\eta s}=1$, given
$|\frac{U_{1}}{4}\pm\sqrt{\gamma^{2}+U_{1}^{2}/4}|>J/2$, the energy levels
from low to high are

$-\frac{U_{1}}{4}-\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}},\qquad-\frac{J}{2},\qquad-\frac{U_{1}}{4}+\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}}\
.$ (S357)

For $n_{\alpha\eta s}=0$, the energy levels from low to high are

$\frac{U_{1}}{4}-\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}},\qquad\frac{J}{2},\qquad\frac{U_{1}}{4}+\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}}\
.$ (S358)

We can choose the chemical potential at zero: For $n_{\alpha\eta s}=1$, the
$-\frac{U_{1}}{4}-\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}}$ and $-J/2$ levels
are occupied, while $-\frac{U_{1}}{4}+\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}}$
is empty; For $n_{\alpha\eta s}=0$, the
$\frac{U_{1}}{4}-\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}}$ level is occupied,
while the $J/2$ and $\frac{U_{1}}{4}+\sqrt{\gamma^{2}+\frac{U_{1}^{2}}{16}}$
levels are empty. The total number of occupied levels (subtracting 12
contributed by $c$-electrons) equals to $\nu_{f}$. Now we integrate out the
third block (local orbitals), as we did in Eq. S350, to derive an effective
Dirac Hamiltonian. For $n_{\alpha\eta s}=1$, there is

$H^{({\rm
eff},\alpha,\eta,s)}(\mathbf{k})\approx\begin{pmatrix}\frac{2\gamma^{2}}{U_{1}}&v_{\star}(k_{x}\eta\sigma_{0}+i(-1)^{\alpha-1}k_{y}\tau_{0})\\\
v_{\star}(k_{x}\eta\sigma_{0}-i(-1)^{\alpha-1}k_{y})&-\frac{J}{2}\end{pmatrix}\
.$ (S359)

Tuning $\frac{2\gamma^{2}}{U_{1}}\to m$, $\frac{J}{2}\to-m$ ($m>0$), the above
Hamiltonian is adiabatically deformed to $m\xi_{z}+v_{\star}\eta
k_{x}\xi_{x}-v_{\star}(-1)^{\alpha-1}\xi_{y}$. Hence the contributed Chern
number by this Hamiltonian is $-\frac{1}{2}\eta(-1)^{\alpha-1}$. For
$n_{\alpha\eta s}=0$, there is

$H^{({\rm
eff},\alpha,\eta,s)}(\mathbf{k})\approx\begin{pmatrix}-\frac{2\gamma^{2}}{U_{1}}&v_{\star}(k_{x}\eta\sigma_{0}+i(-1)^{\alpha-1}k_{y}\tau_{0})\\\
v_{\star}(k_{x}\eta\sigma_{0}-i(-1)^{\alpha-1}k_{y})&+\frac{J}{2}\end{pmatrix}\
.$ (S360)

Tuning $\frac{2\gamma^{2}}{U_{1}}\to m$, $\frac{J}{2}\to-m$ ($m>0$), the above
Hamiltonian is adiabatically deformed to $-m\xi_{z}+v_{\star}\eta
k_{x}\xi_{x}-v_{\star}(-1)^{\alpha-1}\xi_{y}$. Hence the contributed Chern
number by this Hamiltonian is $\frac{1}{2}\eta(-1)^{\alpha-1}$. Therefore, for
the parent state of the form Eq. S353, the total Chern number contributed by
the $n_{\alpha\eta s}=1$ and $n_{\alpha\eta s}=0$ flavors is given by

$C=\sum_{\alpha\eta s}\left(-\frac{1}{2}\eta(-1)^{\alpha-1}n_{\alpha\eta
s}+\frac{1}{2}\eta(-1)^{\alpha-1}(1-n_{\alpha\eta s})\right)=-\sum_{\alpha\eta
s}\eta(-1)^{\alpha-1}n_{\alpha\eta s}\ .$ (S361)

According to the $C_{3z}$ representation matrix of the $f$-orbitals,
$D^{f}(C_{3z})=e^{\frac{2\pi}{3}\tau_{z}\sigma_{z}}$ (Eq. S36), the factor
$\eta(-1)^{\alpha-1}$ is nothing but the angular momentum of the corresponding
$f$-orbital.

Substituting Eq. S361 into the VP states at the fillings $\nu=0,-1,-2,-3$, we
find their Chern numbers as

$C_{\nu=0}=0,\qquad C_{\nu=-1}=-1,\qquad C_{\nu=-2}=0,\qquad C_{\nu=-3}=-1\ .$
(S362)

The flat-U(4) related K-IVC and VP+K-IVC states have the same Chern numbers.
These results are same as Ref. Lian  _et al._ (2021).

##  S5 More discussions

##### Insights for the heavy fermion physics from experiments

— Our model fits precisely with the experimental facts of Ref. Rozen  _et al._
(2021) which writes “The correlated state features an unusual combination of
seemingly contradictory properties, some associated with itinerant electrons -
such as the absence of a thermodynamic gap, metallicity and a Dirac-like
compressibility - and others associated with localized moments, such as a
large entropy and its disappearance under a magnetic field.”

##### Charge-2 excitations

— As an example of the improved understanding of the MATBG physics through our
model, as shown in Ref. Bernevig  _et al._ (2021b), there is a high-energy
“flat-band” charge-2 collective mode upon the ground state at $\nu=0$ (see
Fig. 4 of Ref. Bernevig  _et al._ (2021b)). Since its energy is about 120meV
$\sim 2U_{1}$, it is likely that the charge-2 collective mode is just an
excitation of two $f$-electrons at the same site.

##### Position-dependence of STM spectrum and Landau level quantization

— One key feature of our heavy fermion model is that the quasi-particle
excitation has a minimal band gap at the $\Gamma_{M}$ point, and the
corresponding wavefunctions are contributed by $c$-electrons. This feature
should be reflected in position dependence of the spectrum: At an AA-stacking
site, the local spectrum (density of states) is mainly contributed by
$f$-electrons and hence should have strong weight at the large gap ($\sim U$)
states but less weight at the minimal gap ($\sim J$) states. On the contrary,
at an AB-stacking site, where c-electron dominates, the local spectrum should
have more weight at the minimal gap ($\sim J$) states. In addition, since
$c$-electrons are more delocalized, their response to magnetic field will be
much more significant than that of $f$-electrons. As a result, the Landau
level (LL) quantization in the presence of a magnetic field should be easier
to be observed at AB-stacking sites.

The position-dependences of local spectrum (density of states) and LLs have
been seen in STM experiments, e.g., Refs. Choi  _et al._ (2021); Wong  _et
al._ (2020). Comparing Fig. 1c and Extended Fig. 2a of Ref. Choi  _et al._
(2021) one can see that, in the same sample, when the large gap ($\sim U$) is
well developed, there are additional peaks inside the large gap ($\sim U$) at
AB-stacking sites but none at AA-stacking sites. This phenomenon is consistent
with the above analysis. Ref. Choi  _et al._ (2021) also found that LLs at AA-
stacking sites are harder to resolve than those at AB-stacking sites, which
are also consistent the analysis above. Comparing Fig. 1c and Fig. 3a of Ref.
Wong  _et al._ (2020), one can see the same behavior: Spectrum at AB-stacking
sites have additional peaks inside the large gap ($\sim U$). (But Fig. 1c and
Fig. 3a of Ref. Wong  _et al._ (2020) are from two samples. We did not find
spectra at both AA and AB sites of the same sample in Ref. Wong  _et al._
(2020).)

The quasi-particle spectra shown in Fig. 2 of the main text also explain the
Landau fans observed in transport experiments Wu  _et al._ (2021); Saito  _et
al._ (2021a); Das  _et al._ (2021); Park  _et al._ (2021). Around integer
fillings $\nu\geq 0$ ($\nu\leq 0$), the quasi-particles (quasi-holes) are
contributed by c-electrons and will form Landau levels in the presence of
magnetic field, while quasi-holes (quasi-particles) are contributed by
f-electrons and will form Landau levels. Therefore, Landau fans are observed
at $\nu=n+\delta$ ($\nu=-n-\delta$) for $n=0,1,2,3$ and $0<\delta\ll 1$.

##### Ground states at fillings $\nu=\pm 3$

— Transport experiments Yankowitz  _et al._ (2019); Lu  _et al._ (2019)
observed resistivity peaks at the fillings $\nu=\pm 3$, suggesting existence
of strongly correlated phases. Then, STM Nuckolls  _et al._ (2020); Choi  _et
al._ (2021) and transport Wu  _et al._ (2021); Saito  _et al._ (2021a); Das
_et al._ (2021); Park  _et al._ (2021) experiments found that a finite
magnetic field can stabilize a correlated Chern insulator phase with Chern
numbers $\nu=\pm 1$ at $\nu=\pm 3$. These experiments suggest that the Chern
insulator phase is at least a metastable state slightly above the actual
ground state.

A perturbation Lian  _et al._ (2021) and an ED calculation Xie  _et al._
(2021) by some of the authors found that the ground states at $\nu=\pm 3$
would be correlated Chern insulator with Chern number $\pm 1$ if the so-called
flat-metric-condition is imposed to the model. These results are also
consistent analysis in Ref. Bultinck  _et al._ (2020a). However, with
realistic parameters, there is no indirect gap between the quasi-hole and
electron bands - even though there is a direct gap at each momentum Bernevig
_et al._ (2021b). Thus, with realistic parameters, the correlated Chern
insulator phase becomes unstable Bernevig  _et al._ (2021b). The gapless
spectrum obtained in the current work Fig. S9 is fully consistent with the
results in Refs. Bernevig  _et al._ (2021b); Xie  _et al._ (2021).

DMRG studies Kang and Vafek (2020); Soejima  _et al._ (2020) found that, using
realistic parameters, gapless $C_{2z}T$-symmetric nematic or gapped
$C_{2z}T$-symmetric stripe states may have lower energies than the Chern