# Topological Mixed Valence Model in Magic-Angle Twisted Bilayer Graphene

Yantao Li  Department of Physics, Kent State University, Kent, Ohio 44242, USA
Benjamin M. Fregoso  Department of Physics, Kent State University, Kent, Ohio
44242, USA  Maxim Dzero  Department of Physics, Kent State University, Kent,
Ohio 44242, USA

###### Abstract

We develop a model to describe the mixed valence regime in magic-angle twisted
bilayer graphene (MATBG) using the recently developed heavy-fermion framework.
By employing the large-$N$ slave-boson approach, we derive the self-consistent
mean field equations and solve them numerically. We find that the SU(8)
symmetry constraint moiré system exhibits novel mixed-valence properties which
are different from conventional heavy-fermions systems. We find the solutions
describing the physics at the filling near the Mott insulator regime in the
limit of strong Coulomb interactions between the flat-band fermions. Our model
can provide additional insight into the possible microscopic origin of
unconventional superconductivity in MATBG.

_Introduction_.—The discovery of correlated electronic phases including
superconductivity in magic-angle twisted bilayer graphene (MATBG) Cao  _et
al._ (2018a, b) has stimulated research efforts to explore various electronic
properties in graphene-based multilayer structures Yankowitz  _et al._ (2019);
Sharpe  _et al._ (2019); Fleischmann  _et al._ (2020); Jones  _et al._ (2020);
Cao  _et al._ ; Xu  _et al._ (2021); Choi  _et al._ (2021); Park  _et al._
(2022); Uri  _et al._ (2023) as well as in van der Waals heterostructures and
other platforms Naik and Jain (2018); Tang  _et al._ (2020); Wang  _et al._
(2020); Zhang  _et al._ (2020); Shabani  _et al._ (2021); Xu  _et al._ (2022);
Xiong  _et al._ (2022); Meng  _et al._ (2023). As a result, a new field
"twistronic physics" has emerged Kennes  _et al._ (2021) which focuses on
theoretical aspects of these systems and covers both static Xu and Balents
(2018); Po  _et al._ (2018); Wu  _et al._ (2018a); Kang and Vafek (2018,
2019); Koshino (2019); Hejazi  _et al._ (2019); Chebrolu  _et al._ (2019); Cea
_et al._ (2019); Liu  _et al._ (2019); Liang  _et al._ (2020); Tran  _et al._
(2020); Tritsaris  _et al._ (2020); Ramires and Lado (2021); Lake and Senthil
(2021); Qin and MacDonald (2021); Eaton  _et al._ (2022); Zhou  _et al._
(2022) and non-equilibrium properties Topp  _et al._ (2019); Li  _et al._
(2020); Katz  _et al._ (2020); Vogl  _et al._ (2020, 2021).

The MATBG system consists of two single graphene sheets which are twisted
relative to each other at certain angles called magic angles Bistritzer and
MacDonald (2011); Lopes dos Santos  _et al._ (2007); Shallcross  _et al._
(2008, 2010); Mele (2010, 2011); Carr  _et al._ (2017). It is believed that it
is the flat bands that appear at such magic angles are the main driver for the
exotic physical phenomena which were experimentally observed in these systems.
The origin of magic angles is related to the case when the electron tunneling
in the AA-region of MATBG is neglected corresponding to the chiral limit
Tarnopolsky  _et al._ (2019); Khalaf  _et al._ (2019). It can be theoretically
shown that one does not need to use the chiral limit for other twisted
graphene stacks to exhibit such flat bands associated with Dirac cones Li  _et
al._ (2022). Furthermore, the discovery of superconductivity in MATBG
demonstrates yet another example of superconductivity emerging from the
’strange metal’ phase Cao  _et al._ (2020) and, as such, is reminiscent of the
physics of the high-Tc copper-based and heavy-fermion superconductors. In
passing we note, that several works have recently attempted to explain the
origin of superconductivity in MATBG from various perspectives Wu  _et al._
(2018b); Lian  _et al._ (2019); González and Stauber (2019); Christos  _et
al._ (2022); Fischer  _et al._ (2022).

Most recently an alternative viewpoint has emerged. Specifically, focusing on
a first magic angle, Song and Bernevig showed that a model for twisted bilayer
graphene can be mapped to the heavy-fermion model Song and Bernevig (2022).
This elegant theory is based on the experimental fact that AA-region in MATBG
exhibits a quantum dot-like behavior Tilak  _et al._ (2021) and so one can
describe the physics of this region using a model with the flat band electrons
($f$-electrons). The electrons in the AB/BA regions play the role of
conduction electrons ($c$-electrons). By mapping the Bistritzer-MacDonald
model Bistritzer and MacDonald (2011) to the periodic Anderson model at the
first magic angle, Song and Bernevig created a way to bridge the MATBG with
the more conventional heavy-fermion systems. However, in striking contrast
with the conventional heavy-fermion system, the MATBG hosts SU(8) symmetry due
to spin and valley degrees of freedom and two central flat bands. Shortly
after several important works appeared, which addressed various aspects of
this unique moiré heavy-fermion-like system Chou and Sarma (2022); Yu  _et
al._ (2023); Hu  _et al._ (2023a, b); Zhou and Song (2023); Huang  _et al._
(2023); C?lug?ru  _et al._ (2023); Singh  _et al._ (2023); Chou and Sarma
(2023).

In conventional heavy-fermion systems, one usually distinguishes between the
so-called local moment (or Kondo) regime and the mixed-valent one. In the
local moment regime, the energy of the flat ($f$-orbital) band lies well below
the Fermi energy of the conduction electrons, while in the mixed-valent
regime, it lies close to the Fermi energy. In order to describe both of these
regimes on a technical level, one considers the limit when the local Coulomb
repulsion is taken to infinity. Then one introduces the projection operators
along with the Lagrange multiplier to enforce the constraint of the single
occupancy on the $f$-levels. In the mean-field approximation, one replaces the
projection operators and constraint fields with the $c$-numbers which are
computed self-consistently and describe the renormalization of the $f$-energy
level, hybridization between $c$ and $f$-electrons as well as renormalization
of the chemical potential.

Early works on the application of the heavy-fermion model to MATBG have
typically focused on the Kondo regime and have not solved the problem self-
consistently. In this paper, we attempt to resolve this issue. In what
follows, we present a new heavy-fermion model for the MATBG and discuss the
physics both in the Kondo regime and the mixed valence regime self-
consistently. Our self-consistent solutions pave the way to go beyond the
mean-field approximation and consider the effects of fluctuations on the
competing phases of MATBG.

An important aspect of our model is that we consider the strong limit of
Coulomb interactions, this could provide another way to explain the
microscopic origin of unconventional ($d$-wave) superconductivity, which can
be shown to arise from the quantum mechanical fluctuations in the number of
the $f$-electrons and is purely electronic in origin. Indeed, as the very
recent experiment suggests Oh  _et al._ (2021), the MATBG does indicate the
$d$-wave superconductivity similar to the Cooper pair symmetry of the
superconductivity in the conventional heavy fermion systems Lavagna  _et al._
(1987).

_Setup and Topological Heavy Fermion Model_.—The MATBG setup is illustrated in
Fig. 1(a). Two single graphene layers are stacked together and they are
twisted relative to each other with angle $\theta$. The AA regions are
designated by white spots (with red arrows depicting flat-band electrons) and
AB/BA regions are dark green (with blue spheres depicting conduction
electrons).

Next, we consider the first magic angle $\theta_{m}=1.05^{\circ}$. At this
magic angle, the MATBG can be mapped from the Bistritzer-MacDonald model to a
special SU(8) periodic Anderson model, i.e. the Song-Bernevig model (SB) Song
and Bernevig (2022). We choose the SB model as the starting point. The
Hamiltonian is

$\hat{H}=\hat{H}_{0,\textrm{c}}+\hat{H}_{0,\textrm{f}}+\hat{H}_{0,\textrm{cf}}+\hat{H}_{\textrm{U}},$
(1)

where

$\hat{H}_{0\textrm{c}}=\sum_{a,a^{\prime},\eta,s}\sum_{\mathbb{p}}h^{(c,\eta
s)}_{aa^{\prime}}(\mathbb{p})c^{\dagger}_{\mathbb{p},a,\eta,s}c_{\mathbb{p},a^{\prime},\eta,s},$
(2)

is the Hamiltonian of the conduction ($c$-) electrons in AB/BA moiré lattice
sites,

$\hat{H}_{0\textrm{f}}=\sum_{\alpha,\alpha^{\prime},\eta,s}\sum_{\mathbb{k}}h^{(f,\eta
s)}_{\alpha\alpha^{\prime}}(\mathbb{k})f^{\dagger}_{\mathbb{k},\alpha,\eta,s}f_{\mathbb{k},\alpha^{\prime},\eta,s},$
(3)

is the Hamiltonian of flat band ($f$-) electrons in AA moiré lattice sites,

$\displaystyle\hat{H}_{0\textrm{cf}}=$
$\displaystyle\sum_{\alpha,a,\eta,s}\sum_{\mathbb{G}}\sum_{\mathbb{k}\in\text{mBZ}}[V^{(\eta
s)}_{\alpha
a}(\mathbb{k}+\mathbb{G})f^{\dagger}_{\mathbb{k},\alpha,\eta,s}c_{\mathbb{k}+\mathbb{G},a,\eta,s}$
(4) $\displaystyle+\text{H.c.}],$

accounts for the hybridization between the $c$ and $f$-electrons and

$\hat{H}_{\textrm{U}}=\frac{U}{2}\sum_{\mathbb{R}}:\hat{n}^{f}_{\mathbb{R}}::\hat{n}^{f}_{\mathbb{R}}:,$
(5)

is the Hamiltonian describing the local Coulomb repulsion between the $f$\-
electrons.

Figure 1: (a) Sketch of the setup. Two single graphene layers twist relative
to each other with angle $\theta=\theta_{m}$. The AA regions behave like
quantum dots (red arrows) and AB/BA regions behave like conduction states
(blue balls). (b) Moiré momentum space with 7 mBZs. The blue dashed lines mean
the path of the spectrum. (c) and (d) are spectrum with parameters $\rho=0$
(before hybridization) and $\rho=0.5$ (full hybridization) separately. Also,
both of them are set with $\mu_{f}=\lambda=0$ and other parameters can be seen
in the main text.

In the expressions above $\alpha=1,2$, $a=1,2,3,4$, $\eta=\pm$, and
$s=\uparrow,\downarrow$ are flat band, conduction band, valley, and spin
indexes correspondingly, $\mathbb{p}=\mathbb{k}+\mathbb{G}$, $\mathbb{G}$ is
the reciprocal lattice vectors in the moiré momentum space,
$\hat{n}^{f}_{\mathbb{R}}$ is the on-site density operator of $f$-electrons,
$U$ is the strength of the Coulomb repulsion. Lastly, the matrices which
appear in the expressions above are defined according to

$\begin{split}\hat{h}^{(c,\eta
s)}(\mathbb{p})&=\left[\begin{matrix}-\mu_{c}\hat{\sigma}_{0}&\nu_{\star}(\eta
p_{x}\hat{\sigma}_{0}+ip_{y}\hat{\sigma}_{z})\\\ \nu_{\star}(\eta
p_{x}\hat{\sigma}_{0}-ip_{y}\hat{\sigma}_{z})&M\hat{\sigma}_{x}-\mu_{c}\hat{\sigma}_{0}\end{matrix}\right],\\\
V^{(\eta
s)}(\mathbb{p})&=e^{\frac{-|\mathbb{p}|^{2}\lambda^{2}}{2}}\begin{bmatrix}\gamma\sigma_{0}+\nu_{\star}^{\prime}(\eta
p_{x}\sigma_{x}+p_{y}\sigma_{y})\\\ 0_{2\times 2}\\\
\end{bmatrix},\end{split}$ (6)

and $h^{(f,\eta s)}=(\epsilon_{f_{0}}-\mu_{f})\hat{\sigma}_{0}$, where
$\mu_{c}$ and $\mu_{f}$ are the chemical potentials for $c$ and $f-$ electrons
accordingly, $\hat{\sigma}_{0}=I_{2\times 2}$ is the unit matrix,
$\hat{\sigma}_{j}$ ($j=x,y,z$) are the Pauli matrices. Note that the
Hamiltonian is expressed in the moiré momentum space using the plane wave
approximation. The size of moiré momentum space (per valley per spin) is
$2+4N_{G}$ which means $2$ flat bands and $4$ conduction bands with $N_{G}$
moiré Brillouin zones (mBZ). The values of the parameters are
$\nu_{\star}=-4.303$ eVÅ, $M=3.697$ meV, $\gamma=-24.75$ meV,
$\nu_{\star}^{\prime}=1.622$ eVÅ, and the dampling factor $\lambda=0.3375$
$a_{M}$, where $a_{M}$ is the moiré lattice constant. Note that all these
parameters correspond to ’magic angle’ $\theta_{m}=1.05^{\circ}$,
$U_{0}=W_{AA}/W_{AB}=0.8$, $W_{AA}$ and $W_{AB}$ are interlayer hopping
amplitudes in AA regions, and AB/BA regions separately, and the velocity of
the electron in single layer graphene $v_{F}=5.94$ eVÅ. We note that when
$M=0$ the system reaches the limit of the flat band and we do not consider
this case in this paper.

_Slave-Boson approach and Mean Field Equations_.—To handle the SB model in the
mixed valence region, we use the slave-boson approach. We extend the number of
orbital, spin, and valley for both $c$\- and $f$-electrons to $N$ flavors and
set the interaction

$U=\infty$ (7)

to exclude the double occupancy. Introducing the slave-boson operators
$b^{\dagger}_{\mathbb{R}}$ and $b_{\mathbb{R}}$ at each AA site in real space
the constraint becomes

$Q=\sum^{N}_{l=1}\sum_{\alpha}(f^{\dagger}_{\mathbb{R},\alpha,l}f_{\mathbb{R},\alpha,l}+b^{\dagger}_{\mathbb{R}}b_{\mathbb{R}}).$
(8)

Note that $N=4$ includes spin and valley for two flat band orbitals. Here, we
suppose different valley has the same band structure and we choose the valley
$\eta=+$ in the above Hamiltonian. We stress that there are two index spaces
in the Hamiltonian. One is the extended $N$ space with $1/N$ expansion
corresponding to the index $l$. Another is the moiré momentum space which has
a size of $2+4N_{G}$ corresponding to the index $\alpha$, $\alpha^{\prime}$,
$a$, and $a^{\prime}$ separately

We introduce Lagrangian multipliers $\lambda_{\mathbb{R}}$ to ensure the
number of $f$ electrons is $Q=1$. We rewrite $Q\rightarrow q_{0}N$,
$b_{\mathbb{R}}\rightarrow b_{\mathbb{R}}\sqrt{N}$, and $V_{\alpha
a}\rightarrow V_{\alpha a}/\sqrt{N}$. The local gauge transformation is
$b_{\mathbb{R}}=\rho_{\mathbb{R}}\exp(i\theta_{\mathbb{R}})$,
$f_{\mathbb{R}}=f^{\prime}_{\mathbb{R}}\exp(i\theta_{\mathbb{R}})$ and
$\lambda_{\mathbb{R}}=\lambda^{\prime}_{\mathbb{R}}-\theta_{\mathbb{R}}$. We
then rewrite $f^{\prime}_{\mathbb{R}}$ and $\lambda^{\prime}_{\mathbb{R}}$ to
$f_{\mathbb{R}}$ and $\lambda_{\mathbb{R}}$. The partition function is

$Z=\int\mathscr{D}(cc^{\dagger}ff^{\dagger}\rho\lambda)\exp(-S),$ (9)

where the action is $S=\int_{0}^{\beta}L(\tau)d\tau$, and

$\displaystyle L=\sum^{N}_{\begin{subarray}{c}a,a^{\prime};l=1\\\
\mathbb{p}\end{subarray}}(\partial_{\tau}+h^{(c)}_{aa^{\prime}}(\mathbb{p}))c^{\dagger}_{\mathbb{p},a,l}c_{\mathbb{p},a^{\prime},l}(\tau)+\sum_{\begin{subarray}{c}\alpha,\alpha^{\prime}\\\
\mathbb{\mathbb{k},\mathbb{k}^{\prime}\in\text{mBZ}}\end{subarray}}(\partial_{\tau}+h^{(f)}_{\alpha\alpha^{\prime}}(\mathbb{k})\delta_{\mathbb{k},\mathbb{k}^{\prime}}+\frac{i\lambda(\mathbb{k}-\mathbb{k}^{\prime};\tau)}{\sqrt{N_{L}}})f^{\dagger}_{\mathbb{k},\alpha,l}f_{\mathbb{k}^{\prime},\alpha^{\prime},l}(\tau)+\frac{1}{\sqrt{N_{L}}}\sum^{N}_{\begin{subarray}{c}\mathbb{G},\alpha,a;l=1\\\
\mathbb{\mathbb{k},\mathbb{k}^{\prime}\in\text{mBZ}}\end{subarray}}$ (10)
$\displaystyle\left.[V_{\alpha
a}(\mathbb{k}+\mathbb{G})c^{\dagger}_{\mathbb{k}+\mathbb{G},a,l}(\tau)f_{\mathbb{k}^{\prime},\alpha,l}(\tau)\rho^{\dagger}(\mathbb{k}-\mathbb{k}^{\prime};\tau)+\text{H.c.}]\right.+\left.\frac{iN}{\sqrt{N_{L}}}\sum_{\mathbb{k},\mathbb{k}^{\prime}\in\text{mBZ}}\rho({\mathbb{k}};\tau)\lambda(\mathbb{k}^{\prime}-\mathbb{k};\tau)\rho(-\mathbb{k}^{\prime};\tau)\right.\left.-iq_{0}N\sqrt{N_{L}}\lambda(0;\tau)\right.,$

where $\beta=\frac{1}{T}$, $T$ is the temperature, and $N_{L}$ is the number
of lattice site in moiré real space. We set
$\lambda(\mathbb{k;\tau})=\frac{1}{T}\bar{\lambda}\delta_{\mathbb{k},0}$,
$\rho(\mathbb{k};\tau)=\frac{1}{T}\bar{\rho}\delta_{\mathbb{k},0}$ to get the
mean field action. We rewrite
$\frac{\bar{\lambda}}{\sqrt{N_{L}}}\rightarrow\bar{\lambda}$ and
$\frac{\bar{\rho}}{\sqrt{N_{L}}}\rightarrow\bar{\rho}$. After applying
$\partial S_{0}/\partial\bar{\rho}=0$, $\partial
S_{0}/\partial\bar{\lambda}=0$, and summing over the Matsubara frequency
$\omega_{n}$, we end up with three mean field equations as follows

$\bar{\rho}\bar{\lambda}=\frac{i}{2N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{\bar{\rho}}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j}),$ (11)
$q_{0}-\bar{\rho}^{2}=-\frac{i}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{\bar{\lambda}}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j}),$ (12)
$n_{t}=q_{0}-\bar{\rho}^{2}+\frac{1}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{c}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j})-2N_{G},$ (13)

where $n_{t}$ is the total filling,
$n_{F}(\epsilon)=(\exp(\epsilon/T)+1)^{-1}$ is the Fermi-Dirac distribution,
$\mathcal{E}_{j}$ is the eigenvalues of matrix
$A_{0}=\\{\\{\hat{h}^{(c)},\hat{V}\\},\\{\hat{V}^{\dagger},\hat{h}^{(f)}+i\bar{\lambda}\sigma_{0}\\}\\}$,
$A^{\bar{\rho}}_{0}=\partial{(-i\omega_{n}+A_{0})}/\partial\bar{\rho}$,
$A^{\bar{\lambda}}_{0}=\partial{(-i\omega_{n}+A_{0})}/\partial\bar{\lambda}$,
$A^{c}_{0}=-\partial{(-i\omega_{n}+A_{0})}/\partial\mu_{c}$, and $A_{0}$
expands in the moiré momentum space. Note that $P$ can be constructed by the
eigenvectors of $A_{0}$ and
$P=(c_{1},c_{2},\cdots,c_{(2+4N_{G})})_{2+4N_{G}\times 2+4N_{G}}$, where
$c_{j}$ are the eigenvectors of $A_{0}$. (See Supplemental Material Note1 for
the detailed derivation). The above three self-consistent mean field equations
are one of our main results.

_Numerics_.—Now, we will numerically solve the equations Eq. 11, Eq. 12, and
Eq. 13. To solve them self-consistently, we set an error bar
$errs=\sum^{3}_{n=1}(l_{n}-r_{n})^{2}$, where $l_{n}$ and $r_{n}$ represent
left and right sides of $n$th mean field equation separately. We set
$\mu_{c}=\mu_{f}$ and go through the parameters regions $\rho\in[0,0.5]$,
$\mu_{f}\in[-100,100]$ meV, and $i\lambda\in[-100,100]$ meV. Since we set the
interaction $U=\infty$, $Q=1$, so we have $q_{0}=1/4$. To reach the mixed
valence region, we also set the total filling $n_{t}=0.83\times q_{0}$. We
find the solutions to make the $errs\approx 0$. There exist two solutions: one
is positive $\mu_{f}$, and another is negative $\mu_{f}$. We note that
$\mathcal{E}_{j}(\bar{\lambda},\bar{\rho},\mu_{f},\mathbb{k})$ is numerically
calculated and depends on the chemical potential $\mu_{f}$ and momentum
$\mathbb{k}$ with $\mathbb{k}\in\text{mBZ}$.

We substitute the solutions to $A_{0}$ and we get the spectrums, see Fig. 2.
We also plot the variation of parameters in the mixed valence region as a
function of the temperature $T$, see Fig. 2.

Figure 2: (a) Spectrum of the self-consistent solutions with parameters
$\mu_{f}=19$ meV, $i\lambda=22$ meV, $\rho=0.4$, and $T=0.01K$. (b) The
spectrum of the self-consistent solutions with parameters $\mu_{f}=-5$ meV,
$i\lambda=-5$ meV, $\rho=0.1$, and $T=0.01K$. (c) Self-consistent solutions
$\mu_{f}$ and $i\lambda$ vs $T$. (d) Self-consistent solutions $\rho$,
$n_{f}$, and $n_{c}$ vs $T$ with $n_{t}=0.8q_{0}$, where $T_{k}\approx 24K$ is
the Kondo temperature.

_Discussion_.—Different from the conventional heavy fermion system, the
topological mixed valence model in MATBG has SU(8) symmetry. The two central
flat bands for each valley and spin are particle-hole symmetric in the chiral
limit. Although the particle-hole symmetry is broken at the experimental
lattice relaxation range $U_{0}=0.8$, the valence bands host a similar band
structure as the conduction bands before the hybridization. This makes the
difficult to reach the Kondo region just by pushing the flat bands lower away
from the conduction bands as it does in the conventional heavy fermion system.
One needs to consider many valence bands together. We also note that it might
be interesting to consider different flavors of slave boson Dorin and
Schlottmann (1993) since we have a total of 8 flat bands; they are not
degenerate at the experimental range. We leave this task for future work.

_Summary_.—We have introduced a new model to describe the mixed valence region
of the magic-angle twisted bilayer graphene with infinite Coulomb interaction.
We start from the SB model and use the slave-boson method in large-$N$
expansion. We derive a new group of mean field equations to describe the mixed
valence regions of twisted bilayer graphene. The solutions can catch the
physics of the filling near the strong correlation, which is at the edge of
the Mott insulator, and then could be approaching the unconventional
superconductivity. Our topological mixed valence model paves the way to study
the possible origin of superconductivity in twisted bilayer graphene. We hope
our model could stimulate further research in the mixed valence region in
various related van der Waals heterostructures materials or platforms.

_Note added_.—A related paper comes out recently Lau and Coleman (2023), which
deals with the mixed valence model in twisted bilayer graphene with finite
Coulomb interaction.

###### Acknowledgements.

_Acknowledgments._ We would like to acknowledge the very useful discussions
with Yang-Zhi Chou. This work was financially supported by the National
Science Foundation Grants No. NSF-DMR-2002795 (Y. L. and M. D.) and NSF-
DMR-2015639 (B.M.F.). Parts of this paper were written during the Aspen Center
of Physics 2023 summer program on “New Directions on Strange Metals in
Correlated Systems" (M. D.), which was supported by the National Science
Foundation Grant No. PHY-2210452.

## References

  * Cao  _et al._ (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018a).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Yankowitz  _et al._ (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Sharpe  _et al._ (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Fleischmann  _et al._ (2020) M. Fleischmann, R. Gupta, F. Wullschläger, S. Theil, D. Weckbecker, V. Meded, S. Sharma, B. Meyer, and S. Shallcross, Nano Letters 20, 971 (2020).
  * Jones  _et al._ (2020) A. J. H. Jones, R. Muzzio, P. Majchrzak, S. Pakdel, D. Curcio, K. Volckaert, D. Biswas, J. Gobbo, S. Singh, J. T. Robinson, K. Watanabe, T. Taniguchi, T. K. Kim, C. Cacho, N. Lanata, J. A. Miwa, P. Hofmann, J. Katoch, and S. Ulstrup, Advanced Materials 32, 2001656 (2020).
  * (7) Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, “Large pauli limit violation and reentrant superconductivity in magic-angle twisted trilayer graphene,” arXiv:2103.12083 .
  * Xu  _et al._ (2021) S. Xu, M. M. Al Ezzi, N. Balakrishnan, A. Garcia-Ruiz, B. Tsim, C. Mullan, J. Barrier, N. Xin, B. A. Piot, T. Taniguchi, K. Watanabe, A. Carvalho, A. Mishchenko, A. K. Geim, V. I. Fal’ko, S. Adam, A. H. C. Neto, K. S. Novoselov, and Y. Shi, Nature Physics 17, 619 (2021).
  * Choi  _et al._ (2021) Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Nature 589, 536 (2021).
  * Park  _et al._ (2022) J. M. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature Materials 21, 877 (2022).
  * Uri  _et al._ (2023) A. Uri, S. C. de la Barrera, M. T. Randeria, D. Rodan-Legrain, T. Devakul, P. J. D. Crowley, N. Paul, K. Watanabe, T. Taniguchi, R. Lifshitz, L. Fu, R. C. Ashoori, and P. Jarillo-Herrero, “Superconductivity and strong interactions in a tunable moiré quasiperiodic crystal,” (2023), arXiv:2302.00686 [cond-mat.mes-hall] .
  * Naik and Jain (2018) M. H. Naik and M. Jain, Physical Review Letters 121, 266401 (2018).
  * Tang  _et al._ (2020) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Nature 579, 353 (2020).
  * Wang  _et al._ (2020) L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Nature Materials 19, 861 (2020).
  * Zhang  _et al._ (2020) L. Zhang, Z. Zhang, F. Wu, D. Wang, R. Gogna, S. Hou, K. Watanabe, T. Taniguchi, K. Kulkarni, T. Kuo, S. R. Forrest, and H. Deng, Nature Communications 11, 5888 (2020).
  * Shabani  _et al._ (2021) S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, Nature Physics 17, 720 (2021).
  * Xu  _et al._ (2022) Y. Xu, K. Kang, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, Nature Nanotechnology 17, 934 (2022).
  * Xiong  _et al._ (2022) Y. Xiong, Y. Wang, R. Zhu, H. Xu, C. Wu, J. Chen, Y. Ma, Y. Liu, Y. Chen, K. Watanabe, T. Taniguchi, M. Shi, X. Chen, Y. Lu, P. Zhan, Y. Hao, and F. Xu, Science Advances 8, eabo0375 (2022), https://www.science.org/doi/pdf/10.1126/sciadv.abo0375 .
  * Meng  _et al._ (2023) Z. Meng, L. Wang, W. Han, F. Liu, K. Wen, C. Gao, P. Wang, C. Chin, and J. Zhang, Nature 615, 231 (2023).
  * Kennes  _et al._ (2021) D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio, Nature Physics 17, 155 (2021).
  * Xu and Balents (2018) C. Xu and L. Balents, Physical Review Letters 121, 087001 (2018).
  * Po  _et al._ (2018) H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Physical Review X 8, 031089 (2018).
  * Wu  _et al._ (2018a) F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. 121, 026402 (2018a).
  * Kang and Vafek (2018) J. Kang and O. Vafek, Physical Review X 8, 031088 (2018).
  * Kang and Vafek (2019) J. Kang and O. Vafek, Physical Review Letters 122, 246401 (2019).
  * Koshino (2019) M. Koshino, Physical Review B 99, 235406 (2019).
  * Hejazi  _et al._ (2019) K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Phys. Rev. B 99, 035111 (2019).
  * Chebrolu  _et al._ (2019) N. R. Chebrolu, B. L. Chittari, and J. Jung, Phys. Rev. B 99, 235417 (2019).
  * Cea  _et al._ (2019) T. Cea, N. R. Walet, and F. Guinea, Phys. Rev. B 100, 205113 (2019).
  * Liu  _et al._ (2019) J. Liu, Z. Ma, J. Gao, and X. Dai, Phys. Rev. X 9, 031021 (2019).
  * Liang  _et al._ (2020) X. Liang, Z. A. H. Goodwin, V. Vitale, F. Corsetti, A. A. Mostofi, and J. Lischner, Physical Review B 102, 155146 (2020).
  * Tran  _et al._ (2020) K. Tran, J. Choi, and A. Singh, 2D Materials 8, 022002 (2020).
  * Tritsaris  _et al._ (2020) G. A. Tritsaris, S. Carr, Z. Zhu, Y. Xie, S. B. Torrisi, J. Tang, M. Mattheakis, D. T. Larson, and E. Kaxiras, 2D Materials 7, 035028 (2020).
  * Ramires and Lado (2021) A. Ramires and J. L. Lado, Phys. Rev. Lett. 127, 026401 (2021).
  * Lake and Senthil (2021) E. Lake and T. Senthil, Phys. Rev. B 104, 174505 (2021).
  * Qin and MacDonald (2021) W. Qin and A. H. MacDonald, Phys. Rev. Lett. 127, 097001 (2021).
  * Eaton  _et al._ (2022) A. Eaton, Y. Li, H. A. Fertig, and B. Seradjeh, Phys. Rev. B 106, 045117 (2022).
  * Zhou  _et al._ (2022) Y. Zhou, D. N. Sheng, and E.-A. Kim, Phys. Rev. Lett. 128, 157602 (2022).
  * Topp  _et al._ (2019) G. E. Topp, G. Jotzu, J. W. McIver, L. Xian, A. Rubio, and M. A. Sentef, Physical Review Research 1, 023031 (2019).
  * Li  _et al._ (2020) Y. Li, H. A. Fertig, and B. Seradjeh, Physical Review Research 2, 043275 (2020).
  * Katz  _et al._ (2020) O. Katz, G. Refael, and N. H. Lindner, Physical Review B 102, 155123 (2020).
  * Vogl  _et al._ (2020) M. Vogl, M. Rodriguez-Vega, and G. A. Fiete, Phys. Rev. B 101, 235411 (2020).
  * Vogl  _et al._ (2021) M. Vogl, M. Rodriguez-Vega, B. Flebus, A. H. MacDonald, and G. A. Fiete, Phys. Rev. B 103, 014310 (2021).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Lopes dos Santos  _et al._ (2007) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Shallcross  _et al._ (2008) S. Shallcross, S. Sharma, and O. A. Pankratov, Phys. Rev. Lett. 101, 056803 (2008).
  * Shallcross  _et al._ (2010) S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Phys. Rev. B 81, 165105 (2010).
  * Mele (2010) E. J. Mele, Phys. Rev. B 81, 161405 (2010).
  * Mele (2011) E. J. Mele, Phys. Rev. B 84, 235439 (2011).
  * Carr  _et al._ (2017) S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Phys. Rev. B 95, 075420 (2017).
  * Tarnopolsky  _et al._ (2019) G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Khalaf  _et al._ (2019) E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019).
  * Li  _et al._ (2022) Y. Li, A. Eaton, H. A. Fertig, and B. Seradjeh, Phys. Rev. Lett. 128, 026404 (2022).
  * Cao  _et al._ (2020) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020).
  * Wu  _et al._ (2018b) F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121, 257001 (2018b).
  * Lian  _et al._ (2019) B. Lian, Z. Wang, and B. A. Bernevig, Phys. Rev. Lett. 122, 257002 (2019).
  * González and Stauber (2019) J. González and T. Stauber, Phys. Rev. Lett. 122, 026801 (2019).
  * Christos  _et al._ (2022) M. Christos, S. Sachdev, and M. S. Scheurer, Phys. Rev. X 12, 021018 (2022).
  * Fischer  _et al._ (2022) A. Fischer, Z. A. H. Goodwin, A. A. Mostofi, J. Lischner, D. M. Kennes, and L. Klebl, npj Quantum Materials 7, 5 (2022).
  * Song and Bernevig (2022) Z.-D. Song and B. A. Bernevig, Phys. Rev. Lett. 129, 047601 (2022).
  * Tilak  _et al._ (2021) N. Tilak, X. Lai, S. Wu, Z. Zhang, M. Xu, R. d. A. Ribeiro, P. C. Canfield, and E. Y. Andrei, Nature Communications 12, 4180 (2021).
  * Chou and Sarma (2022) Y.-Z. Chou and S. D. Sarma, “Kondo lattice model in magic-angle twisted bilayer graphene,” (2022), arXiv:2211.15682 [cond-mat.str-el] .
  * Yu  _et al._ (2023) J. Yu, M. Xie, B. A. Bernevig, and S. D. Sarma, “Magic-angle twisted symmetric trilayer graphene as topological heavy fermion problem,” (2023), arXiv:2301.04171 [cond-mat.mes-hall] .
  * Hu  _et al._ (2023a) H. Hu, G. Rai, L. Crippa, J. Herzog-Arbeitman, D. C?lug?ru, T. Wehling, G. Sangiovanni, R. Valenti, A. M. Tsvelik, and B. A. Bernevig, “Symmetric kondo lattice states in doped strained twisted bilayer graphene,” (2023a), arXiv:2301.04673 [cond-mat.str-el] .
  * Hu  _et al._ (2023b) H. Hu, B. A. Bernevig, and A. M. Tsvelik, “Kondo lattice model of magic-angle twisted-bilayer graphene: Hund’s rule, local-moment fluctuations, and low-energy effective theory,” (2023b), arXiv:2301.04669 [cond-mat.str-el] .
  * Zhou and Song (2023) G.-D. Zhou and Z.-D. Song, “Kondo phase in twisted bilayer graphene – a unified theory for distinct experiments,” (2023), arXiv:2301.04661 [cond-mat.str-el] .
  * Huang  _et al._ (2023) C. Huang, X. Zhang, G. Pan, H. Li, K. Sun, X. Dai, and Z. Meng, “Evolution from quantum anomalous hall insulator to heavy-fermion semimetal in magic-angle twisted bilayer graphene,” (2023), arXiv:2304.14064 [cond-mat.str-el] .
  * C?lug?ru  _et al._ (2023) D. C?lug?ru, M. Borovkov, L. L. H. Lau, P. Coleman, Z.-D. Song, and B. A. Bernevig, “Tbg as topological heavy fermion: Ii. analytical approximations of the model parameters,” (2023), arXiv:2303.03429 [cond-mat.str-el] .
  * Singh  _et al._ (2023) K. Singh, A. Chew, J. Herzog-Arbeitman, B. A. Bernevig, and O. Vafek, “Topological heavy fermions in magnetic field,” (2023), arXiv:2305.08171 [cond-mat.str-el] .
  * Chou and Sarma (2023) Y.-Z. Chou and S. D. Sarma, “Scaling theory of intrinsic kondo and hund’s rule interactions in magic-angle twisted bilayer graphene,” (2023), arXiv:2306.03121 [cond-mat.str-el] .
  * Oh  _et al._ (2021) M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Nature 600, 240 (2021).
  * Lavagna  _et al._ (1987) M. Lavagna, A. J. Millis, and P. A. Lee, Phys. Rev. Lett. 58, 266 (1987).
  * Dorin and Schlottmann (1993) V. Dorin and P. Schlottmann, Phys. Rev. B 47, 5095 (1993).
  * Lau and Coleman (2023) L. L. H. Lau and P. Coleman, “Topological mixed valence model for twisted bilayer graphene,” (2023), arXiv:2303.02670 [cond-mat.str-el] .

Supplemental Material for “Topological Mixed Valence Model in Magic-Angle
Twisted Bilayer Graphene”

##  I Hamiltonian and Parameters

The periodic Anderson model in twisted bilayer graphene Song and Bernevig
(2022) can be expressed as

$\hat{H}=\hat{H}_{0,c}+\hat{H}_{0,f}+\hat{H}_{0,cf}+\hat{H}_{U},$ (S1)

where $\hat{H}_{0,c}$ is the Hamiltonian of the conduction ($c$-) electrons in
AB/BA moiré lattice sites, $\hat{H}_{0,f}$ is the Hamiltonian of flat band
($f$-) electrons in AA moiré lattice sites, $\hat{H}_{0,cf}$ is the coupling
Hamiltonian of them, and $\hat{H}_{U}$ is the Hamiltonian of Coulomb
interaction among the $f$\- electrons. The explicit form of them are

$\hat{H}_{0,c}=\sum_{a,a^{\prime},\eta,s}\sum_{\mathbb{p}}h^{(c,\eta
s)}_{aa^{\prime}}(\mathbb{p})c^{\dagger}_{\mathbb{p},a,\eta,s}c_{\mathbb{p},a^{\prime},\eta,s}=\sum_{a,a^{\prime},\eta,s}\sum_{\mathbb{G}}\sum_{\mathbb{k}\in\text{mBZ}}h^{(c,\eta
s)}_{aa^{\prime}}(\mathbb{k}+\mathbb{G})c^{\dagger}_{\mathbb{k}+\mathbb{G},a,\eta,s}c_{\mathbb{k}+\mathbb{G},a^{\prime},\eta,s},$
(S2)

$\hat{H}_{0,f}=\sum_{\alpha,\alpha^{\prime},\eta,s}\sum_{\mathbb{k}}h^{(f,\eta
s)}_{\alpha\alpha^{\prime}}(\mathbb{k})f^{\dagger}_{\mathbb{k},\alpha,\eta,s}f_{\mathbb{k},\alpha^{\prime},\eta,s},$
(S3)

$\hat{H}_{0,cf}=\sum_{\alpha,a,\eta,s}\sum_{\mathbb{G}}\sum_{\mathbb{k}\in\text{mBZ}}[V^{(\eta
s)}_{\alpha
a}(\mathbb{k}+\mathbb{G})f^{\dagger}_{\mathbb{k},\alpha,\eta,s}c_{\mathbb{k}+\mathbb{G},a,\eta,s}+\text{H.c.}],$
(S4)

and

$\hat{H}_{U}=\frac{U}{2}\sum_{\mathbb{R}}:\hat{n}^{f}_{\mathbb{R}}::\hat{n}^{f}_{\mathbb{R}}:,$
(S5)

where $\mathbb{p}=\mathbb{k}+\mathbb{G}$, $\hat{n}^{f}_{\mathbb{R}}$ is the
on-site density operator of $f$-electrons and $U$ is Coulomb interaction among
them. The single particle Hamiltonian of the model can also be expressed as

$\begin{aligned} \hat{H_{0}}&=\begin{bmatrix}h^{(c,\eta
s)}(\mathbb{p})&V^{(\eta s)}(\mathbb{p})\\\ (V^{(\eta
s)}(\mathbb{p}))^{\dagger}&h^{(f,\eta s)}(\mathbb{k})\\\ \end{bmatrix}\\\
&=\begin{bmatrix}0_{2\times 2}-\mu\sigma_{0}&\nu_{\star}(\eta
p_{x}\sigma_{0}+ip_{y}\sigma_{z})&\gamma\sigma_{0}+\nu_{\star}^{\prime}(\eta
p_{x}\sigma_{x}+p_{y}\sigma_{y})\\\ \nu_{\star}(\eta
p_{x}\sigma_{0}-ip_{y}\sigma_{z})&M\sigma_{x}-\mu\sigma_{0}&0_{2\times 2}\\\
\gamma\sigma_{0}+\nu_{\star}^{\prime}(\eta
p_{x}\sigma_{x}+p_{y}\sigma_{y})&0_{2\times
2}&(\epsilon_{f_{0}}-\mu)\sigma_{0}\end{bmatrix}\\\ \\\ \end{aligned},$ (S6)

where

$\displaystyle h^{(c,\eta s)}(\mathbb{p})$
$\displaystyle=\begin{bmatrix}0_{2\times 2}-\mu\sigma_{0}&\nu_{\star}(\eta
p_{x}\sigma_{0}+ip_{y}\sigma_{z})\\\ \nu_{\star}(\eta
p_{x}\sigma_{0}-ip_{y}\sigma_{z})&M\sigma_{x}-\mu\sigma_{0}\\\ \end{bmatrix},$
(S7)

$\displaystyle V^{(\eta s)}(\mathbb{p})$
$\displaystyle=e^{\frac{-|\mathbb{p}|^{2}\lambda^{2}}{2}}\begin{bmatrix}\gamma\sigma_{0}+\nu_{\star}^{\prime}(\eta
p_{x}\sigma_{x}+p_{y}\sigma_{y})\\\ 0_{2\times 2}\\\ \end{bmatrix},$ (S8)

and $h^{(f,\eta s)}=(\epsilon_{f_{0}}-\mu)\sigma_{0}$.

The Hamiltonian is expressed in the moiré momentum space using the plane wave
approximation. The size of moiré momentum space (per valley per spin) is
$2+4N_{G}$ which means $2$ flat bands and $4$ conduction bands and conduction
bands are expanded to $N_{G}$ moiré Brillouin zone (mBZ), see Fig. S1. Then
the single particle Hamiltonian can be written in matrix form as

$\displaystyle\hat{H}_{0}=\left(\begin{array}[]{cccc}\hat{H}_{0,c}&\hat{H}_{0,cf}\\\
(\hat{H}_{0,cf})^{\dagger}&\hat{H}_{0,f}\\\ \end{array}\right),$ (S9)

where

$\displaystyle\hat{H}_{0,c}=\left(\begin{array}[]{cccc}h^{(c,\eta
s)}(\mathbb{k}+\mathbb{G_{1}})&0&\cdots&0\\\ 0&h^{(c,\eta
s)}(\mathbb{k}+\mathbb{G_{2}})&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&h^{(c,\eta s)}(\mathbb{k}+\mathbb{G_{N_{G}}})\end{array}\right),$
(S10)

$\displaystyle\hat{H}_{0,cf}=\left(\begin{array}[]{cccc}V^{(\eta s)}_{\alpha
a}(\mathbb{k}+\mathbb{G_{1}})&V^{(\eta s)}_{\alpha
a}(\mathbb{k}+\mathbb{G_{2}})&\cdots&V^{(\eta s)}_{\alpha
a}(\mathbb{k}+\mathbb{G_{N_{G}}})\end{array}\right),$ (S11)

and $\hat{H}_{0,f}=h^{(f,\eta s)}.$ Note that $\mathbb{G_{n}=\Gamma_{n}}$ and
$\mathbb{\Gamma_{n}}$ is the position of symmetry point in $n$-th moiré
Brillouin zone, for example, $\mathbb{G_{1}=\Gamma_{1}}=0$.

The values of the parameters are $\nu_{\star}=-4.303$ eVÅ, $M=3.697$ meV,
$\gamma=-24.75$ meV, $\nu_{\star}^{\prime}=1.622$ eVÅ, and $\lambda=0.3375$
$a_{M}$, where $a_{M}$ is the moiré lattice constant. Note that all these
parameters are corresponding to $\theta_{m}=1.05^{\circ}$,
$U_{0}=W_{AA}/W_{AB}=0.8$, and the velocity of electron in single layer
graphene $v_{F}=5.94$ eVÅ.

Figure S1: Sketch of moiré momentum space which includes $N_{G}=37$ moiré
Brillouin zones. $K$, $K^{\prime}$ and $\Gamma$ are the symmetry points in the
moiré Brillouin zone. In numerical calculation, we increase $N_{G}$ from $9$
to $19$ and $37$ to check the convergence.

##  II Slave-boson Mean field equations in Large-N Expansion

We extend the number of orbital, spin, and valley for both $c$\- and
$f$-electrons to $N$ and set

$U=\infty.$ (S12)

By introducing the slave-boson operators $b^{\dagger}_{\mathbb{R}}$ and
$b_{\mathbb{R}}$ at each AA site to exclude the double occupancy. The
Hamiltonian in the large-N approximation is

$\displaystyle\hat{H}$
$\displaystyle=\hat{H}_{0,c}+\hat{H}_{0,f}+\hat{H}_{0,cf}$ (S13)
$\displaystyle=\sum_{l}^{N}\left\\{\sum_{a,a^{\prime}}\sum_{\mathbb{p}}h^{(c)}_{aa^{\prime}}(\mathbb{p})c^{\dagger}_{\mathbb{p},a,l}c_{\mathbb{p},a^{\prime},l}+\sum_{\alpha,\alpha^{\prime}}\sum_{\mathbb{k}}h^{(f)}_{\alpha\alpha^{\prime}}(\mathbb{k})f^{\dagger}_{\mathbb{k},\alpha,l}f_{\mathbb{k},\alpha^{\prime},l}\right\\}$
$\displaystyle+\frac{1}{\sqrt{N_{L}}}\sum_{l}^{N}\left\\{\sum_{\alpha,a}\sum_{\mathbb{G}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\mathbb{q}}[V_{\alpha
a}(\mathbb{k}+\mathbb{G})c^{\dagger}_{\mathbb{k}+\mathbb{G},a,l}f_{\mathbb{k}+\mathbb{q},\alpha,l}b^{\dagger}_{\mathbb{q}}+\text{H.c.}]\right\\}.$

The constraint is

$Q=\sum^{N}_{l=1}\sum_{\alpha}(f^{\dagger}_{\mathbb{R},\alpha,l}f_{\mathbb{R},\alpha,l}+b^{\dagger}_{\mathbb{R}}b_{\mathbb{R}}).$
(S14)

Note that $N$ includes orbital, spin, and valley. Here, we suppose different
valley has the same band structure and we choose the valley $\eta=+$ in the
above Hamiltonian. In our case, $N=4$ and $N=\infty$ is the mean field value.
We stress that there are two index spaces in the Hamiltonian. One is the
extended $N$ space with $1/N$ expansion corresponding to the index $l$.
Another is the moiré momentum space which has size of $N_{G}$ corresponding to
the index $a,a^{\prime},\alpha$, and $\alpha^{\prime}$. We only consider one
species of the slave boson for simplification. Since we have two flat bands,
two species of the slave bosons Dorin and Schlottmann (1993) might be needed
and we will consider this case in future works.

We introduce Lagrangian multipliers $\lambda_{\mathbb{R}}$ to ensure the
number of $f$ electrons is $Q$ or $N_{f}$, where $Q\equiv N_{f}$. We rewrite
$Q\rightarrow q_{0}N$, $b_{\mathbb{R}}\rightarrow b_{\mathbb{R}}\sqrt{N}$, and
$V_{\alpha a}\rightarrow V_{\alpha a}/\sqrt{N}$. The local gauge
transformation is
$b_{\mathbb{R}}=\rho_{\mathbb{R}}\exp(i\theta_{\mathbb{R}})$,
$f_{\mathbb{R}}=f^{\prime}_{\mathbb{R}}\exp(i\theta_{\mathbb{R}})$ and
$\lambda_{\mathbb{R}}=\lambda^{\prime}_{\mathbb{R}}-\theta_{\mathbb{R}}$. We
then rewrite $f^{\prime}_{\mathbb{R}}$ and $\lambda^{\prime}_{\mathbb{R}}$ to
$f_{\mathbb{R}}$ and $\lambda_{\mathbb{R}}$. The partition function is

$Z=\int\mathscr{D}(cc^{\dagger}ff^{\dagger}\rho\lambda)\exp(-S),$ (S15)

where the action is $S=\int_{0}^{\beta}Ld\tau$

$\displaystyle L=\sum^{N}_{\begin{subarray}{c}a,a^{\prime};l=1\\\
\mathbb{p}\end{subarray}}(\partial_{\tau}+h^{(c)}_{aa^{\prime}}(\mathbb{p}))c^{\dagger}_{\mathbb{p},a,l}c_{\mathbb{p},a^{\prime},l}(\tau)+\sum_{\begin{subarray}{c}\alpha,\alpha^{\prime}\\\
\mathbb{\mathbb{k},\mathbb{k}^{\prime}\in\text{mBZ}}\end{subarray}}(\partial_{\tau}+h^{(f)}_{\alpha\alpha^{\prime}}(\mathbb{k})\delta_{\mathbb{k},\mathbb{k}^{\prime}}+\frac{i\lambda(\mathbb{k}-\mathbb{k}^{\prime};\tau)}{\sqrt{N_{L}}})f^{\dagger}_{\mathbb{k},\alpha,l}f_{\mathbb{k}^{\prime},\alpha^{\prime},l}(\tau)+\frac{1}{\sqrt{N_{L}}}\sum^{N}_{\begin{subarray}{c}\mathbb{G},\alpha,a;l=1\\\
\mathbb{\mathbb{k},\mathbb{k}^{\prime}\in\text{mBZ}}\end{subarray}}$ (S16)
$\displaystyle\left.[V_{\alpha
a}(\mathbb{k}+\mathbb{G})c^{\dagger}_{\mathbb{k}+\mathbb{G},a,l}(\tau)f_{\mathbb{k}^{\prime},\alpha,l}(\tau)\rho^{\dagger}(\mathbb{k}-\mathbb{k}^{\prime};\tau)+\text{H.c.}]\right.+\left.\frac{iN}{\sqrt{N_{L}}}\sum_{\mathbb{k},\mathbb{k}^{\prime}\in\text{mBZ}}\rho({\mathbb{k}};\tau)\lambda(\mathbb{k}^{\prime}-\mathbb{k};\tau)\rho(-\mathbb{k}^{\prime};\tau)\right.\left.-iq_{0}N\sqrt{N_{L}}\lambda(0;\tau)\right.$
,

where $\beta=\frac{1}{T}$ and $T$ is the temperature. We set
$\lambda(\mathbb{k;\tau})=\frac{1}{T}\bar{\lambda}\delta_{\mathbb{k},0}$,
$\rho(\mathbb{k};\tau)=\frac{1}{T}\bar{\rho}\delta_{\mathbb{k},0}$, then the
mean field action is

$\displaystyle S_{0}=$
$\displaystyle-N\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\ln^{\det{(-i\omega_{n}+A_{0})}}+\frac{iN}{T}(\bar{\rho}^{2}/\sqrt{N_{L}}-q_{0}\sqrt{N_{L}})\bar{\lambda}$
(S17) $\displaystyle=$
$\displaystyle-N\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum_{j=1}^{2+4N_{G}}\ln^{(-i\omega_{n}+\mathcal{E}_{j})}+\frac{iN}{T}(\bar{\rho}^{2}/\sqrt{N_{L}}-q_{0}\sqrt{N_{L}})\bar{\lambda}$
$\displaystyle=$
$\displaystyle-N\sum_{\mathbb{k}\in\text{mBZ}}\sum_{j=1}^{2+4N_{G}}\ln^{(1+e^{\mathcal{E}_{j}/T})}+\frac{iN}{T}(\bar{\rho}^{2}/\sqrt{N_{L}}-q_{0}\sqrt{N_{L}})\bar{\lambda},$

where $\mathcal{E}_{j}$ is the eigenvalues of $A_{0}$ and $A_{0}$ is a matrix
with size of $2+4N_{G}$. We rewrite
$\frac{\bar{\lambda}}{\sqrt{N_{L}}}\rightarrow\bar{\lambda}$ and
$\frac{\bar{\rho}}{\sqrt{N_{L}}}\rightarrow\bar{\rho}$, then we have

$\displaystyle S_{0}=$
$\displaystyle-N\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\ln^{\det{(-i\omega_{n}+A_{0})}}+\frac{iN}{T}N_{L}(\bar{\rho}^{2}-q_{0})\bar{\lambda}$
(S18) $\displaystyle=$
$\displaystyle-N\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum_{j=1}^{2+4N_{G}}\ln^{(-i\omega_{n}+\mathcal{E}_{j})}+\frac{iN}{T}N_{L}(\bar{\rho}^{2}-q_{0})\bar{\lambda}$
$\displaystyle=$
$\displaystyle-N\sum_{\mathbb{k}\in\text{mBZ}}\sum_{j=1}^{2+4N_{G}}\ln^{(1+e^{\mathcal{E}_{j}/T})}+s\frac{iN}{T}N_{L}(\bar{\rho}^{2}-q_{0})\bar{\lambda},$
$\displaystyle
A_{0}=\begin{bmatrix}\begin{array}[]{ccccc}h^{(c)}(\mathbb{k}+\mathbb{G_{1}})&0&\cdots&0&\bar{\rho}V(\mathbb{k}+\mathbb{G_{1}})\\\
0&h^{(c)}(\mathbb{k}+\mathbb{G_{2}})&\cdots&0&\bar{\rho}V(\mathbb{k}+\mathbb{G_{2}})\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&h^{(c)}(\mathbb{k}+\mathbb{G_{N_{G}}})&\bar{\rho}V(\mathbb{k}+\mathbb{G_{N_{G}}})\\\
(\bar{\rho}V(\mathbb{k}+\mathbb{G_{1}}))^{\dagger}&(\bar{\rho}V(\mathbb{k}+\mathbb{G_{2}}))^{\dagger}&\cdots&(\bar{\rho}V(\mathbb{k}+\mathbb{G_{N_{G}}}))^{\dagger}&h^{(f)}(\mathbb{k})+i\bar{\lambda}\sigma_{0}\\\
\end{array}\end{bmatrix},$ (S19)

we relabel

$\displaystyle\hat{h}^{(c)}=\left(\begin{array}[]{cccc}h^{(c)}(\mathbb{k}+\mathbb{G_{1}})&0&\cdots&0\\\
0&h^{(c)}(\mathbb{k}+\mathbb{G_{2}})&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&h^{(c)}(\mathbb{k}+\mathbb{G_{N_{G}}})\end{array}\right)$ (S20)

and

$\displaystyle\hat{V}=\left(\begin{array}[]{cccc}\bar{\rho}V(\mathbb{k}+\mathbb{G_{1}})&\bar{\rho}V(\mathbb{k}+\mathbb{G_{2}})&\cdots&\bar{\rho}V(\mathbb{k}+\mathbb{G_{N_{G}}})\end{array}\right).$
(S21)

So, we can rewrite

$\displaystyle
A_{0}=\begin{bmatrix}\begin{array}[]{ccccc}\hat{h}^{(c)}&\hat{V}\\\
\hat{V}^{\dagger}&h^{(f)}+i\bar{\lambda}\sigma_{0}\\\
\end{array}\end{bmatrix},$ (S22)

After applying $\partial S_{0}/\partial\bar{\rho}=0$ and $\partial
S_{0}/\partial\bar{\lambda}=0$, we have the saddle point mean field equations
as follows

$\bar{\lambda}=-i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\frac{1}{\det{(-i\omega_{n}+A_{0})}}\frac{\partial(\det(-i\omega_{n}+A_{0}))}{(2\bar{\rho})\partial\bar{\rho}},$
(S23)

$q_{0}-\bar{\rho}^{2}=i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\frac{1}{\det{(-i\omega_{n}+A_{0})}}\frac{\partial(\det(-i\omega_{n}+A_{0}))}{\partial\bar{\lambda}}.$
(S24)

Using the following formula

$\frac{\partial{(\det{(C}))}}{\partial
t}=(\det{(C}))\cdot\mathrm{Tr}((C)^{-1}\frac{\partial{C}}{\partial t}),$ (S25)

where $C$ is a matrix and $t$ is a scalar. We have

$\frac{\partial{(\det{(-i\omega_{n}+A_{0}}))}}{\partial\bar{\rho}}=(\det{(-i\omega_{n}+A_{0}}))\cdot\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}\frac{\partial{(-i\omega_{n}+A_{0})}}{\partial\bar{\rho}})$
(S26)

and

$\frac{\partial{(\det{(-i\omega_{n}+A_{0}}))}}{\partial\bar{\lambda}}=(\det{(-i\omega_{n}+A_{0}}))\cdot\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}\frac{\partial{(-i\omega_{n}+A_{0})}}{\partial\bar{\lambda}})$
(S27)

.

Then the mean-field equations become

$\bar{\lambda}=-i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}\frac{\partial{(-i\omega_{n}+A_{0})}}{(2\bar{\rho})\partial\bar{\rho}}),$
(S28)

$q_{0}-\bar{\rho}^{2}=i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}\frac{\partial{(-i\omega_{n}+A_{0})}}{\partial\bar{\lambda}}).$
(S29)

We define

$\displaystyle
A^{\bar{\rho}}_{0}=\frac{\partial{(-i\omega_{n}+A_{0})}}{\partial\bar{\rho}}=\begin{bmatrix}\begin{array}[]{ccccc}0&0&\cdots&0&V(\mathbb{k}+\mathbb{G_{1}})\\\
0&0&\cdots&0&V(\mathbb{k}+\mathbb{G_{2}})\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&0&V(\mathbb{k}+\mathbb{G_{N_{G}}})\\\
(V(\mathbb{k}+\mathbb{G_{1}}))^{\dagger}&(V(\mathbb{k}+\mathbb{G_{2}}))^{\dagger}&\cdots&(V(\mathbb{k}+\mathbb{G_{N_{G}}}))^{\dagger}&0\\\
\end{array}\end{bmatrix}$ (S30)

and

$\displaystyle
A^{\bar{\lambda}}_{0}=\frac{\partial{(-i\omega_{n}+A_{0})}}{\partial\bar{\lambda}}=\begin{bmatrix}\begin{array}[]{ccccc}0&0&\cdots&0&0\\\
0&0&\cdots&0&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\cdots&0&0\\\
0&0&\cdots&0&i\sigma_{0}\\\ \end{array}\end{bmatrix}.$ (S31)

The mean field equations are written as

$\bar{\lambda}=-i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}\frac{A^{\bar{\rho}}_{0}}{2\bar{\rho}}),$
(S32)

$q_{0}-\bar{\rho}^{2}=i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}A^{\bar{\lambda}}_{0}).$
(S33)

Actually, we can get the above equations from the action directly before we do
$\mathrm{Tr}\ln=\ln\det$. We introduce an invertible matrix $P$ to diagonal
the matrix $A_{0}$. Then $P$ can be constructed by the eigenvectors of $A_{0}$
and $P=(c_{1},c_{2},\cdots,c_{(2+4N_{G})})_{2+4N_{G}\times 2+4N_{G}}$, where
$c_{j}$ are eigenvectors of $A_{0}$. Since

$\displaystyle\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}A^{\bar{\rho}}_{0})=\mathrm{Tr}(P^{-1}(-i\omega_{n}+A_{0})^{-1}PP^{-1}A^{\bar{\rho}}_{0}P)=\sum^{2+4N_{G}}_{j}((-i\omega_{n}+\mathcal{E}_{j})^{-1}(P^{-1}A^{\bar{\rho}}_{0}P)_{jj})$
(S34)

and

$\begin{aligned}
&\mathrm{Tr}((-i\omega_{n}+A_{0})^{-1}A^{\bar{\lambda}}_{0})=\mathrm{Tr}(P^{-1}(-i\omega_{n}+A_{0})^{-1}PP^{-1}A^{\bar{\lambda}}_{0}P)=\sum^{2+4N_{G}}_{j}((-i\omega_{n}+\mathcal{E}_{j})^{-1}(P^{-1}A^{\bar{\lambda}}_{0}P)_{jj})\\\
\end{aligned},$ (S35)

then the mean-field equations become

$\bar{\lambda}\bar{\rho}=\frac{i}{2}\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum^{2+4N_{G}}_{j}\frac{(P^{\dagger}A^{\bar{\rho}}_{0}P)_{jj}}{i\omega_{n}-\mathcal{E}_{j}},$
(S36)

$q_{0}-\bar{\rho}^{2}=-i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum^{2+4N_{G}}_{j}\frac{(P^{\dagger}A^{\bar{\lambda}}_{0}P)_{jj}}{i\omega_{n}-\mathcal{E}_{j}}.$
(S37)

Note that $P$ is an unitary matrix and $P^{-1}=P^{\dagger}$. We also have the
third mean field equation,

$n_{t}=q_{0}-\bar{\rho}^{2}+\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum^{2+4N_{G}}_{j}\frac{(P^{\dagger}A^{c}_{0}P)_{jj}}{i\omega_{n}-\mathcal{E}_{j}}-2N_{G},$
(S38)

where $n_{t}=\text{constant}$ is the total number of flat band and conduction
band electrons and we set $n_{t}=1.66$ below. We define

$\displaystyle
A^{c}_{0}=\begin{bmatrix}\begin{array}[]{ccccc}\mathbb{1}_{4\times
4}&0&\cdots&0&0\\\ 0&\mathbb{1}_{4\times 4}&\cdots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\cdots&\mathbb{1}_{4\times 4}&0\\\
0&0&\cdots&0&0_{2\times 2}\\\ \end{array}\end{bmatrix}.$ (S39)

After summing over the Matsubara frequency $\omega_{n}$, we have

$\bar{\lambda}\bar{\rho}=\frac{i}{2N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{\bar{\rho}}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j}),$ (S40)

$q_{0}-\bar{\rho}^{2}=-\frac{i}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{\bar{\lambda}}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j}),$ (S41)

$n_{t}=q_{0}-\bar{\rho}^{2}+\frac{1}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{c}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j})-2N_{G},$ (S42)

where $n_{F}(\epsilon)=\frac{1}{\exp(\epsilon/T)+1}$ is the Fermi-Dirac
distribution. Note that
$\mathcal{E}_{j}(\bar{\lambda},\bar{\rho},\mu,\mathbb{k})$ is numerically
calculated and depends on the chemical potential $\mu$ and momentum
$\mathbb{k}$ with $\mathbb{k}\in\text{mBZ}$.

##  III Green Functions

Now, let us think about the Green functions.

$G_{0}^{-1}=i\omega_{n}-A_{0}.$ (S43)

So,

$\displaystyle G^{ff}_{0}$
$\displaystyle=-i\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum^{2+4N_{G}}_{j}\frac{(P^{\dagger}A^{\bar{\lambda}}_{0}P)_{jj}}{i\omega_{n}-\mathcal{E}_{j}},$
(S44) $\displaystyle G^{cc}_{0}$
$\displaystyle=\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum^{2+4N_{G}}_{j}\frac{(P^{\dagger}A^{c}_{0}P)_{jj}}{i\omega_{n}-\mathcal{E}_{j}},$
(S45) $\displaystyle G^{cf}_{0}$
$\displaystyle=\frac{1}{2}\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum^{2+4N_{G}}_{j}\frac{(P^{\dagger}A^{\gamma}_{0}P)_{jj}}{i\omega_{n}-\mathcal{E}_{j}}$
(S46) $\displaystyle G^{fc}_{0}$
$\displaystyle=\frac{1}{2}\frac{T}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum_{\omega_{n}}\sum^{2+4N_{G}}_{j}\frac{(P(A^{\gamma}_{0})^{\dagger}P^{\dagger})_{jj}}{i\omega_{n}-\mathcal{E}_{j}},$
(S47)

where

$\displaystyle
A^{\gamma}_{0}=\begin{bmatrix}\begin{array}[]{ccccc}0&0&\cdots&0&\mathbb{1}^{0}_{4\times
2}\\\ 0&0&\cdots&0&\mathbb{1}^{0}_{4\times 2}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\cdots&0&\mathbb{1}^{0}_{4\times
2}\\\ \mathbb{1}^{0}_{2\times 4}&\mathbb{1}^{0}_{2\times
4}&\cdots&\mathbb{1}^{0}_{2\times 4}&0\\\ \end{array}\end{bmatrix}$ (S48)

and

$\displaystyle\mathbb{1}^{0}_{4\times
2}=\begin{bmatrix}\begin{array}[]{ccccc}1&0\\\ 0&1\\\ 0&0\\\ 0&0\\\
\end{array}\end{bmatrix}$ (S49)

with

$\displaystyle\mathbb{1}^{0}_{2\times
4}=\begin{bmatrix}\begin{array}[]{ccccc}1&0&0&0\\\ 0&1&0&0\\\
\end{array}\end{bmatrix}.$ (S50)

##  IV Numerics

In this section, we give the data details of the numerics. We also give more
numerical details. First of all, the three mean field equations Eq. S51, S52,
and S51 are coupled with each other. As one approaches the Kondo region which
means $\rho\approx 0$, the Eq. S51 will be automatically satisfied. One no
longer needs to consider Eq. S51 in the Kondo region. Second, to solve the
three equations self-consistently, we numerically calculate
$\mathcal{E}_{j}(\bar{\lambda},\bar{\rho},\mu_{f},\mathbb{k})$ by
diagonalizing $A_{0}$ for each value of $\bar{\lambda}$, $\bar{\rho}$, and
$\mu_{f}$. We set a certain value of temperature $T$ and set $\mu_{c}=\mu_{f}$
then go through the parameters regions $\rho\in[0,0.5]$,
$\mu_{f}\in[-100,100]$ meV, and $i\lambda\in[-100,100]$ meV. Since we set the
interaction $U=\infty$, $Q=1$, so we have $q_{0}=1/4$. We find the solutions
to make the $errs\approx 0$. There exist two solutions: one is positive
$\mu_{f}$, and another is negative $\mu_{f}$. We note that momentum
$\mathbb{k}$ belongs to the first mBZ with $\mathbb{k}\in\text{mBZ}$. To check
the convergence, we perform $N_{G}$ from $3$ to $37$ and it turns out
$3~{}\text{mBZs}$ are good enough for the convergence.

$\bar{\rho}\bar{\lambda}=\frac{i}{2N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{\bar{\rho}}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j}),$ (S51)
$q_{0}-\bar{\rho}^{2}=-\frac{i}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{\bar{\lambda}}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j}),$ (S52)
$n_{t}=q_{0}-\bar{\rho}^{2}+\frac{1}{N_{L}}\sum_{\mathbb{k}\in\text{mBZ}}\sum^{2+4N_{G}}_{j}(P^{\dagger}A^{c}_{0}P)_{jj}\cdot
n_{F}(\mathcal{E}_{j})-2N_{G},$ (S53)

The following data are for the self-consistent solutions of $n_{t}=0.8q_{0}$.
As $\rho\approx 0.24$ at $T=24K$, $n_{c}\approx 0$ which means $T=24K$ is the
Kondo temperature for $n_{t}=0.8q_{0}$. The reason why $\rho$ is not zero is
because we set $n_{t}=0.8q_{0}$, this means parts of $n_{c}$ are in higher
energy conduction bands.

Table 1: Temperature Dependence for $n_{t}=0.8q_{0}$. $n_{t}$ | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8
---|---|---|---|---|---|---|---|---|---|---
$q_{0}$ | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4
$T$ | $0.01K$ | $0.4K$ | $0.8K$ | $1.2K$ | $1.6K$ | $2.0K$ | $4.0K$ | $6.0K$ | $8.0K$ | $10.0K$
$\rho$ | 0.4 | 0.4 | 0.41 | 0.4 | 0.4 | 0.4 | 0.4 | 0.39 | 0.4 | 0.39
$\mu_{f}$ | 19 | 19 | 20 | 19 | 20 | 19 | 19 | 18 | 20 | 18
$i\lambda$ | $22$ | $22$ | $23$ | $22$ | $23$ | $22$ | $22$ | $21$ | $24$ | $22$
$n_{f}$ | $0.09$ | $0.09$ | $0.11$ | $0.09$ | $0.09$ | $0.09$ | $0.09$ | $0.09$ | $0.09$ | $0.09$
$n_{c}$ | $0.11$ | $0.11$ | $0.09$ | $0.11$ | $0.11$ | $0.11$ | $0.11$ | $0.11$ | $0.11$ | $0.11$
$errs$ | $0.015$ | $0.015$ | $0.0084$ | $0.011$ | $0.01$ | 0.0087 | 0.013 | 0.0081 | 0.0058 | 0.0052
Table 2: Temperature Dependence for $n_{t}=0.8q_{0}$. $n_{t}$ | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8 | 1/4$\times$0.8
---|---|---|---|---|---|---|---|---|---|---|---
$q_{0}$ | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4 | 1/4
$T$ | $12K$ | $14K$ | $16K$ | $18K$ | $20K$ | $22K$ | $24K$ | $26K$ | $28K$ | $29K$ | $30K$
$\rho$ | 0.4 | 0.39 | 0.37 | 0.35 | 0.32 | 0.28 | 0.24 | 0.19 | 0.13 | 0.09 | 0
$\mu_{f}$ | 20 | 18 | 15 | 12 | 6 | 3 | 0 | -1 | -3 | -3.4 | -4
$i\lambda$ | $25$ | $23$ | $20$ | $17$ | $11$ | 8 | $5$ | $4$ | 2 | 1.6 | 1
$n_{f}$ | $0.09$ | $0.09$ | $0.11$ | $0.13$ | $0.15$ | $0.17$ | $0.19$ | $0.21$ | $0.23$ | $0.24$ | $0.25$
$n_{c}$ | $0.11$ | $0.11$ | $0.09$ | $0.07$ | $0.05$ | $0.03$ | $0.01$ | $-0.01$ | $-0.03$ | $-0.04$ | $-0.05$
$errs$ | $0.0047$ | $0.0071$ | $0.0031$ | $0.0032$ | $0.0043$ | 0.0069 | $0.0059$ | $0.0064$ | 0.0047 | 0.0024 | 0.0034