# Kondo Phase in Twisted Bilayer Graphene - A Unified Theory for Distinct
Experiments

Geng-Dong Zhou  International Center for Quantum Materials, School of Physics,
Peking University, Beijing 100871, China  Zhi-Da Song<EMAIL_ADDRESS>International Center for Quantum Materials, School of Physics, Peking
University, Beijing 100871, China

###### Abstract

A number of interesting physics phenomena have been discovered in magic-angle
twisted bilayer graphene (MATBG), such as superconductivity, correlated gapped
and gapless phases, etc. The gapped phases are believed to be symmetry-
breaking states described by mean-field theories. Gapless phases, however,
exhibit exotic phenomena beyond naive mean-field descriptions. Such phenomena
include (i) zero-energy peaks in spectral density at low temperatures, (ii) a
cascade of transitions as that of a quantum dot at higher temperatures, (iii)
the Pomeranchuk effect where local moment develops upon heating, (iv)
resistance peaks in gapless phases at certain integer fillings, etc. These
phenomena have never been connected or explained on a microscopic level by any
prior theory. In this work, we point out that all these phenomena result from
a simple unified mechanism - the Kondo effect. We applied systematic
analytical and numerical renormalization group analyses to a single-impurity
version of the recently proposed topological heavy fermion (THF) model of
MATBG. We find that the Fermi liquid ground state exhibiting (i) zero-energy
peaks is stabilized by the Kondo screening effect. A higher temperature will
drive the system into a local moment phase that obeys Curie’s law and has
higher entropy. This explains the (ii) transition cascade, (iii) the
Pomeranchuk effect, and (iv) resistance peaks. Remarkably, using realistic
parameters given by the previous THF model work, the computed spectral
densities, entropies, and spin susceptibilities are quantitatively comparable
to experiments. Furthermore, we predict that the ground state in a wide range
of fractional fillings is the heavy Fermi liquid, and compute its dispersion,
quasi-particle weight, and energy surfaces that could be used to verify our
theory if compared to experiments in the future studies. We conjecture that
the heavy Fermi liquid is the parent state of the observed unconventional
superconductivity. It is also possible that the observed “strange metal”
behavior could be connected to the competition between RKKY interaction and
Kondo screening.

##  I Introduction

After the discovery of superconductivity Cao  _et al._ (2018a) and correlated
insulators Cao  _et al._ (2018b) in magic-angle twisted bilayer graphene
(MATBG) Bistritzer and MacDonald (2011), MATBG has become a platform for
studying new correlation effects in flat-band systems and has received
extensive attention. Remarkably rich physics, including the interplay between
superconductivity Lu  _et al._ (2019); Yankowitz  _et al._ (2019); Saito  _et
al._ (2020, 2020); Stepanov  _et al._ (2020); Liu  _et al._ (2021a); Arora
_et al._ (2020); Cao  _et al._ (2021a); Oh  _et al._ (2021) and strong
correlation Lu  _et al._ (2019); Saito  _et al._ (2020); Stepanov  _et al._
(2020); Liu  _et al._ (2021a); Xie  _et al._ (2019); Sharpe  _et al._ (2019);
Choi  _et al._ (2019); Kerelsky  _et al._ (2019); Jiang  _et al._ (2019); Wong
_et al._ (2020); Zondiner  _et al._ (2020); Choi  _et al._ (2021a),
interaction driven Chern insulators Serlin  _et al._ (2020); Nuckolls  _et
al._ (2020); Choi  _et al._ (2021b); Saito  _et al._ (2021a); Das  _et al._
(2021); Wu  _et al._ (2021); Park  _et al._ (2021), strange metal behaviors
Polshyn  _et al._ (2019); Cao  _et al._ (2020); Jaoui  _et al._ (2022), and
the Pomeranchuk effect Rozen  _et al._ (2021); Saito  _et al._ (2021b), etc.,
have been observed in MATBG. Several theoretical understandings of the
correlated gapped states have also been achieved: The strong correlation
arises from the two topological flat bands Bistritzer and MacDonald (2011);
Song  _et al._ (2019); Po  _et al._ (2019); Ahn  _et al._ (2019); Tarnopolsky
_et al._ (2019); Liu  _et al._ (2019); Song  _et al._ (2021), each of which is
four-fold degenerate due to the spin and valley d.o.f. A large U(4) symmetry
group Kang and Vafek (2019); Bultinck  _et al._ (2020); Seo  _et al._ (2019);
Bernevig  _et al._ (2021a); Vafek and Kang (2020) emerges in the flat-band
limit, where the actual bandwidth is counted as negligible. Then the observed
correlated gapped states can be understood as flavor polarized states Kang and
Vafek (2019); Bultinck  _et al._ (2020); Seo  _et al._ (2019); Vafek and Kang
(2020); Lian  _et al._ (2021); Xie  _et al._ (2021); Xie and MacDonald (2020);
Liu and Dai (2021); Cea and Guinea (2020); Venderbos and Fernandes (2018);
Ochi  _et al._ (2018); Zhang  _et al._ (2020); Liu  _et al._ (2021b); Da Liao
_et al._ (2021); Kang and Vafek (2020); Soejima  _et al._ (2020); Hejazi  _et
al._ (2021); Da Liao  _et al._ (2019); Kennes  _et al._ (2018); Classen  _et
al._ (2019) that spontaneously break the U(4) symmetry. The continuous U(4)
degeneracy also leads to Goldstone mode fluctuations Bernevig  _et al._
(2021b); Khalaf  _et al._ (2020) that may destroy the long-range order due to
the Mermin-Wagner theorem.

Less theoretical understandings have been achieved for the gapless states,
which are observed at both fractional and integer fillings. For example, at
$\nu\approx\pm 1$, depending on the experimental setup, both gapped correlated
insulator Lu  _et al._ (2019); Saito  _et al._ (2020); Stepanov  _et al._
(2020) and Fermi liquid Rozen  _et al._ (2021); Wu  _et al._ (2021); Park  _et
al._ (2021); Saito  _et al._ (2021b, 2020); Polshyn  _et al._ (2019); Jaoui
_et al._ (2022) have been reported, suggesting that they are competing ground
states with close energies. (Here $|\nu|$ is the number of electrons ($\nu>0$)
or holes ($\nu<0$) per moiré cell counted from the charge neutrality point
(CNP).) One interesting feature of the Fermi liquids around $\nu=\pm 1$ is
that they usually exhibit resistance peaks Saito  _et al._ (2021b); Park  _et
al._ (2021); Wu  _et al._ (2021); Polshyn  _et al._ (2019); Jaoui  _et al._
(2022) at temperatures above a few kelvins, and the peaks may become
increasingly pronounced as temperature rises. Another interesting feature of
these Fermi liquids is the zero-energy peak. Scanning tunneling microscope
(STM) measurements Nuckolls  _et al._ (2020); Choi  _et al._ (2021b, a); Oh
_et al._ (2021) have constantly seen that, at temperatures of about a few
hundred millikelvins, the conduction (valence) band will be pinned at the
Fermi level and form a sharp resonance peak at the fillings $1\lesssim\nu<2$
($-2<\nu\lesssim-1$). The peaks do not fit the intuition of Stoner instability
given that the interaction is indeed strong. When the temperature increases to
a few or ten kelvins, these peaks develop into a cascade of transitions as
that of a quantum dot model Wong  _et al._ (2020); Choi  _et al._ (2021a).

In this work, based on the recently developed topological heavy fermion (THF)
model Song and Bernevig (2022); Shi and Dai (2022), we find that the gapless
phase at $1\lesssim|\nu|<2$ are the heavy Fermi liquid, and many related
experiments can be naturally explained by the competing Kondo phase and local
moment (LM) phase. The THF model consists of localized $f$-electrons at the
AA-stacking regions of MATBG, which form a triangular lattice, and itinerant
Dirac $c$-electrons. On the one hand, localized $f$-electrons are subject to
strong on-site Coulomb repulsion (50-60meV) and have almost negligible
hoppings ($\sim$0.1meV). Thus, when the valence fluctuation is small,
$f$-electrons behave like LMs. $c$-electrons, on the other hand, form a weakly
interacting Fermi liquid that tends to screen the LM due to the Kondo effect
Nozières and Blandin (1980); Tsvelick and Wiegmann (1983); Coleman (1984);
Bickers (1987); Affleck and Ludwig (1991a, b); Emery and Kivelson (1992);
Tsvelik and Reizer (1993); Furusaki and Nagaosa (1994); Cox and Zawadowski
(1998); Parcollet  _et al._ (1998); Coleman and Schofield (2005); Kotliar  _et
al._ (2006); Werner  _et al._ (2006); Gegenwart  _et al._ (2008); Si and
Steglich (2010); Dzero  _et al._ (2010); Lu  _et al._ (2013). We find that the
zero-energy peak, as well as the Fermi liquids at $1\lesssim|\nu|<2$, are
results of the Kondo resonance. At a higher temperature exceeding the Kondo
energy scale, the LMs become visible in experiment, giving rise to the
transition cascades and the resistance peaks around $\nu=\pm 1$. This is fully
consistent with the Pomeranchuk effect observed around fillings $\nu=\pm 1$
Rozen  _et al._ (2021); Saito  _et al._ (2021b), which shows the presence of
local moments at $T\sim 10$K. Therefore, we conclude that the gapless states
at the fillings $1\lesssim|\nu|<2$ are the strongly correlated heavy Fermi
liquid. Since this filling range overlaps with the superconductivity Cao  _et
al._ (2018a); Lu  _et al._ (2019); Yankowitz  _et al._ (2019); Saito  _et al._
(2020, 2020); Stepanov  _et al._ (2020); Liu  _et al._ (2021a); Arora  _et
al._ (2020); Cao  _et al._ (2021a); Oh  _et al._ (2021) and the strange metal
Polshyn  _et al._ (2019); Cao  _et al._ (2020); Jaoui  _et al._ (2022) around
$\nu=-2+\delta$ (for small $\delta$), the heavy Fermi liquid could be the
parent state for the unconventional superconductivity, as it is in the heavy
fermion materials with $4f$ or $5f$ electrons. This opens a new perspective -
with solid theoretical and experimental basis at the same time - to study the
superconductivity in MATBG.

Here we summarize the primary results of the current work. For this work to be
self-contained, in Sec. II we will review the THF model and its symmetry
shortly. In Sec. III, based on a poor man’s scaling analysis and several
experimental facts, we argue that the Kondo screening effect is irrelevant at
CNP, and hence the ground state at CNP is the previously identified symmetry-
breaking correlated insulator Kang and Vafek (2019); Bultinck  _et al._
(2020); Seo  _et al._ (2019); Vafek and Kang (2020); Lian  _et al._ (2021);
Xie  _et al._ (2021). Then we derive a simpler effective periodic Anderson
model to describe active excitations upon the correlated ground state. We
further simplify the model to a single-impurity version. In Sec. IV, by
applying poor man’s scaling and Wilson’s numerical renormalization group (NRG)
method Krishna-murthy  _et al._ (1980a, b); Bulla  _et al._ (2008) to the
single-impurity problem, we obtain a phase diagram characterized by strong
coupling fixed points and various LM fixed points. The strong coupling phase
consists of a Kondo regime and a frozen impurity regime, and the gapless
$1\lesssim|\nu|<2$ states are found to be in the Kondo regime. We also present
detailed analyses of the RG flows at these fixed points. In Sec. V we
calculate the spectral densities at various fillings and temperatures. At low
temperatures smaller than the Kondo energy scale, the spectral densities
feature sharp Kondo resonances. At higher temperatures that suppress the Kondo
resonances, the upper (lower) Hubbard bands periodically cross the Fermi level
from above (below) as $\nu$ changes from 0 to 4 (-4), matching the cascade of
transitions in STM experiments. In Sec. VI the spin susceptibility and the
entropy are computed. The spin susceptibilities obey Curie’s law at high
temperatures, suggesting the existence of LM, and approach constants at lower
temperatures, suggesting the Fermi liquid behavior. Kondo temperature $T_{K}$
can be estimated as the turning point between the two behaviors of the spin
susceptibility. For $1\lesssim|\nu|<2$, we find $k_{B}T_{K}$ ranges from
0.175meV to 0.680meV. (Here $k_{B}$ is Boltzmann’s constant.) Due to orbital
and spin d.o.f., the LM at $|\nu|=1$ is four-fold degenerate and hence
contributes to an entropy $\ln 4\cdot k_{B}\approx 1.39k_{B}$. A strong in-
plane magnetic field will quench the spin d.o.f. and reduce the entropy to
$\ln 2\cdot k_{B}\approx 0.69k_{B}$. These values will be appreciated if one
notices that the measured entropies at $\nu\approx 1$ and $T\approx 10$K in
the absence and presence of a strong in-plane field are about $1.2k_{B}$ and
$0.6k_{B}$, respectively Rozen  _et al._ (2021). We also calculate the entropy
as a function of the filling at $T\approx 10$K and obtain a curve
quantitatively comparable to the measured curve in Ref. Rozen  _et al._
(2021). In addition, we predict a correlation between the Pomeranchuk effect
and the spectral density among different samples and experiment setups, which
could be used to verify our theory in future experimental researches. In Sec.
VII, we predict the heavy Fermi liquids at $1\lesssim|\nu|<2$ and propose
experiments to verify them. Using spectral information provided by the NRG
calculation, we explicitly calculated the heavy Fermi liquid bands, their
quasi-particle weights, and energy surfaces at $T=0$ and $T>T_{K}$. We find
that a smoking gun signature of the heavy Fermi liquid state is the striking
temperature dependence of the energy surfaces. We also propose a spin-
selective orbital Kondo effect in the presence of an in-plane magnetic field.
In the end, we discuss the possible effects of the RKKY interactions.

It is worth mentioning that all the parameters used in the current work are
extracted from the continuous model, and all the calculations are done in a
first-principle spirit, i.e., no tuning parameter is introduced to fit
experiments.

##  II The THF model

Figure 1: The THF model. (a) Top: red spheres represent the effective
$f$-electrons located at AA-stacking regions of MATBG, and blue spheres
represent the itinerant $c$-electrons. Bottom: the moiré Brillouin zone. (b)
Black bands are given by the free part of the THF model ($\hat{H}_{0}$ in Sec.
III.2). The red and blue bands are the decoupled $f$\- and $c$-bands,
respectively. $M$ is a parameter that determines the bandwidth of the flat
bands. We focus on the $M\to 0$ limit in this work. (c) Excitation spectrum of
the active $c$-electrons upon the symmetry-breaking parent state for $\nu>0$,
described by $\hat{H}_{c}$ in Eq. 13. (d) The hybridization function
$\Delta(\omega)$ contributed by the $c$-bands in (c). A nonzero $M$ will not
change the asymptotic behavior of $\Delta(\omega)$ around the gap.

One theoretical challenge in studying correlation physics in MATBG is the lack
of a fully symmetric lattice model for low energy physics, which is forbidden
by the band topology protected by a $C_{2z}T$ symmetry Song  _et al._ (2019);
Po  _et al._ (2019); Ahn  _et al._ (2019) and an emergent particle-hole
symmetry $P$ Song  _et al._ (2021) \- even though extended Hubbard models Kang
and Vafek (2018); Koshino  _et al._ (2018); Yuan and Fu (2018); Zou  _et al._
(2018); Zang  _et al._ (2022) can be constructed at the sacrifice of either
symmetry or locality. The band topology was thought as fragile Song  _et al._
(2019); Po  _et al._ (2019); Ahn  _et al._ (2019) but was later shown to be a
stable symmetry anomaly jointly protected by $C_{2z}T$ and $P$ Song  _et al._
(2021). The THF model Song and Bernevig (2022); Shi and Dai (2022) resolved
this problem by ascribing the strong correlation to effective $f$-orbitals at
the AA-stacking regions, which form a triangular lattice, and leaving the
remaining low energy states to continuous $c$-bands described by a topological
Dirac Hamiltonian (Fig. 1). The THF model faithfully reproduces the symmetry,
topology, dispersion, and Coulomb interaction of the continuous Bistritzer-
MacDonald model Bistritzer and MacDonald (2011). Its free part is given by

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

Here $\mu$ is the chemical potential, $\hat{N}$ is the particle-number
operator, $c_{{\mathbf{k}}a\eta s}$ is the fermion operator for the
$c$-electron of the momentum ${\mathbf{k}}$, orbital $a$ ($=1,2,3,4$), valley
$\eta$ ($=\pm$), and spin $s$ ($=\uparrow,\downarrow$),
$f_{{\mathbf{k}}\alpha\eta s}$ is the corresponding fermion operator for the
$f$-electron of the orbital $\alpha$ (=1,2). The summation over ${\mathbf{k}}$
for the $c$-bands is in principle limited up to the cutoff $\Lambda_{c}$. But
the theory is well-defined and yields the same low energy physics if we take
the $\Lambda_{c}\to\infty$ limit. Hence, hereafter we will drop the
restriction $|{\mathbf{k}}|<\Lambda_{c}$.
$H^{(c,\eta)}({\mathbf{k}})=v_{\star}(\eta\sigma_{x}\otimes\sigma_{0}k_{x}-\sigma_{y}\otimes\sigma_{z}k_{y})+\mathbb{0}_{2\times
2}\oplus M\sigma_{x}$ is the Dirac Hamiltonian of the $c$-bands. When $M\neq
0$, $c$-bands have a quadratic band touching at the zero energy, whereas when
$M=0$, $c$-bands become linear. The two-by-two block of
$H^{(cf,\eta)}_{a\alpha}({\mathbf{k}})$ for $a=1,2$ is given by
$\gamma\sigma_{0}+v_{\star}^{\prime}(\eta\sigma_{x}k_{x}+\sigma_{y}k_{y})$,
and the two-by-two block of $H^{(cf,\eta)}_{a\alpha}({\mathbf{k}})$ for
$a=3,4$ vanishes. The parameter $\lambda$ in the second line of Eq. 1 is the
spread of the Wannier functions of $f$-electrons, and it truncates the
hybridization at $|{\mathbf{k}}|\gg\lambda^{-1}$. In this work we adopt the
$w_{0}/w_{1}=0.8$ parameters of Ref. Song and Bernevig (2022):
$\gamma=-24.75$meV, $v_{\star}=-4.303\mathrm{eV\cdot\mathring{A}}$,
$v_{\star}^{\prime}=1.623\mathrm{eV\cdot\mathring{A}}$,
$\lambda=1.4131/k_{\theta}$, $k_{\theta}=1.703\mathring{\rm A}^{-1}\cdot
2\sin\frac{\theta_{m}}{2}$ with $\theta_{m}=1.05^{\circ}$ being the first
magic angle. ($w_{0}$ and $w_{1}$ are the interlayer Slater-Koster couplings
of MATBG at the AA-stacking and AB-stacking regions, respectively. Due to the
corrugation effect Uchida  _et al._ (2014); Van Wijk  _et al._ (2015); Dai
_et al._ (2016); Jain  _et al._ (2016), $w_{0}$ is usually smaller than
$w_{1}$. Ref. Koshino  _et al._ (2018) estimates $w_{0}/w_{1}$ as 0.817.) The
resulting band structure with a nonzero $M$ (3.697meV) is shown in Fig. 1(b).
One can see that the topological flat bands result from the hybridization
between $c$\- and $f$-bands. As explained in detail in Ref. Song and Bernevig
(2022) and shown in Fig. 1(b), the parameter $M$ determines the bandwidth
($2|M|$) of the flat-bands.

In each valley, the Hamiltonian $\hat{H}_{0}$ respects a magnetic space group
$P6^{\prime}2^{\prime}2$ Song  _et al._ (2019) (#177.151 in the BNS setting
Gallego  _et al._ (2012)), generated by $C_{3z}$, $C_{2x}$, $C_{2z}T$, and
translation symmetries. (Single-valley Hamiltonian does not have time-reversal
symmetry $T$ because $T$ will change a valley to another.) The two
$f$-orbitals and the four $c$-bands have effective angular momenta $L=-\eta$,
$\eta$ and $L=-\eta$, $\eta$, $0$, $0$, respectively. As shown in detail in
Song and Bernevig (2022), the six-by-six representations of $C_{3z}$,
$C_{2x}$, $C_{2z}T$ on these orbitals are
$e^{i\eta\frac{2\pi}{3}\sigma_{z}}\oplus
e^{i\eta\frac{2\pi}{3}\sigma_{z}}\oplus\sigma_{0}$, $\mathbb{1}_{3\times
3}\otimes\sigma_{x}$, and $\mathbb{1}_{3\times 3}\otimes\sigma_{x}K$,
respectively, where $\sigma_{0,x,y,z}$ are Pauli matrices and $K$ the complex
conjugation. One can verify that the Hamiltonian matrices $H^{(c,\eta)}$ and
$H^{(cf,\eta)}$ given in the last paragraph respect these crystalline
symmetries. The interaction Hamiltonian given in the following paragraph also
respects these crystalline symmetries.

The interaction Hamiltonian is given by

$\displaystyle\hat{H}_{I}=$
$\displaystyle\frac{U_{1}}{2}\sum_{\mathbf{R}}\delta n^{f}_{\mathbf{R}}\delta
n^{f}_{\mathbf{R}}+\frac{U_{2}}{2}\sum_{\langle\mathbf{R}\mathbf{R}^{\prime}\rangle}\delta
n^{f}_{\mathbf{R}}\delta
n^{f}_{\mathbf{R}^{\prime}}+\frac{1}{2N_{M}}\sum_{\mathbf{q}aa^{\prime}}V(\textbf{q})\delta
n_{-\mathbf{q}a^{\prime}}^{c}\delta
n^{c}_{\mathbf{q}a}+\frac{1}{N_{M}}\sum_{\mathbf{R}\mathbf{q}a}W_{a}e^{-i\mathbf{q}\cdot\mathbf{R}}\delta
n^{f}_{\mathbf{R}}\delta
n^{c}_{\mathbf{q}a}-\frac{J}{2N_{M}}\sum_{\begin{subarray}{c}\eta\eta^{\prime}\alpha\alpha^{\prime}\\\
ss^{\prime}\end{subarray}}\sum_{\begin{subarray}{c}{\mathbf{k}},{\mathbf{k}}^{\prime}\\\
\mathbf{R}\end{subarray}}\bigg{[}$
$\displaystyle(\eta\eta^{\prime}+(-1)^{\alpha+\alpha^{\prime}})e^{-i({\mathbf{k}}-{\mathbf{k}}^{\prime})\cdot\mathbf{R}-\frac{\lambda^{2}({\mathbf{k}}^{2}+{\mathbf{k}}^{\prime
2})}{2}}(f^{\dagger}_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}-\frac{1}{2}\delta_{\eta\eta^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{ss^{\prime}})(c^{\dagger}_{{\mathbf{k}},\alpha+2\eta
s}c_{{\mathbf{k}}^{\prime},\alpha^{\prime}+2\eta^{\prime}s^{\prime}}-\frac{1}{2}\delta_{{\mathbf{k}}{\mathbf{k}}^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{ss^{\prime}})\bigg{]}\
,$ (2)

where $N_{M}$ is the number of moiré cells, $f_{\mathbf{R}\alpha\eta s}$ is
the real space fermion operator for the $f$-electrons, $\mathbf{R}$ are the
triangular lattice shown in Fig. 1(a),
$\langle\mathbf{R}\mathbf{R}^{\prime}\rangle$ represents nearest neighbor
pairs (ordered), $\delta n^{f}_{\mathbf{R}}=\sum_{\alpha\eta
s}(f^{\dagger}_{\mathbf{R}\alpha\eta s}f_{\mathbf{R}\alpha\eta
s}-\frac{1}{2})$ is the total density operator (counted from CNP) of
$f$-electrons at $\mathbf{R}$, $\delta n^{c}_{\mathbf{q}a}=\sum_{\eta
s{\mathbf{k}}}(c^{\dagger}_{{\mathbf{k}}+\mathbf{q}a\eta
s}c_{{\mathbf{k}}a\eta s}-\frac{1}{2}\delta_{\mathbf{q}0})$ is the density
operator (counted from CNP) for $c$-electrons of the orbital $a$. $U_{1,2}$,
$V(\mathbf{q})$, $W_{a}$ are the density-density interaction between $ff$,
$cc$, $cf$ electrons, respectively, and $J$ is an exchange interaction between
$cf$ electrons. We adopt the $w_{0}/w_{1}=0.8$ parameters of Ref. Song and
Bernevig (2022): $U_{1}=57.95$meV, $U_{2}=2.329$meV, $W_{1}=W_{2}=44.03$meV,
$W_{3}=W_{4}=50.20$meV, $J=16.38$meV. We choose $V(\mathbf{q})$ as the double-
gate-screened Coulomb interaction,
$\frac{\pi\xi^{2}U_{\xi}}{\Omega_{0}}\frac{\tanh(\xi|\mathbf{q}|/2)}{\xi|\mathbf{q}|/2}$,
with $\xi=10$nm being distance between MATBG and the gates, $U_{\xi}=24$meV
the Coulomb interaction at the distance $\xi$, and $\Omega_{0}\approx
156\mathrm{nm}^{2}$ the area of the moiré cell.

Hereafter, we mainly focus on the flat-band limit, i.e., $M=0$. In this limit,
an exact U(4) symmetry of $\hat{H}_{0}+\hat{H}_{I}$ between the spin, valley,
and orbital flavors emerge, as previously recognized in the projected Coulomb
Hamiltonian of the continuous model Kang and Vafek (2019); Bultinck  _et al._
(2020); Seo  _et al._ (2019); Bernevig  _et al._ (2021a); Vafek and Kang
(2020). We emphasize that this U(4) symmetry is not related to the so-called
chiral limit Tarnopolsky  _et al._ (2019); Wang  _et al._ (2021), i.e.,
$w_{0}=0$, which leads a distinct U(4) symmetry [Bernevig  _et al._ , 2021a,
Bultinck  _et al._ , 2020]. The U(4) symmetry in the flat-band limit has been
shown as a good approximation for realistic parameters such as
$w_{0}/w_{1}=0.8$ Bernevig  _et al._ (2021a); Lian  _et al._ (2021). The
sixteen U(4) generators acting on $f_{\mathbf{R}\alpha\eta s}$,
$c_{{\mathbf{k}}a\eta s}$ ($a=1,2$), and $c_{{\mathbf{k}}a\eta s}$ ($a=3,4$)
are

$\Sigma^{f}_{\mu\nu}=\\{\sigma_{0}\tau_{0}\varsigma_{\nu},\sigma_{y}\tau_{x}\varsigma_{\nu},\sigma_{y}\tau_{y}\varsigma_{\nu},\sigma_{0}\tau_{z}\varsigma_{\nu}\\}\
,$ (3)
$\Sigma^{c12}_{\mu\nu}=\\{\sigma_{0}\tau_{0}\varsigma_{\nu},\sigma_{y}\tau_{x}\varsigma_{\nu},\sigma_{y}\tau_{y}\varsigma_{\nu},\sigma_{0}\tau_{z}\varsigma_{\nu}\\}\
,$ (4)

and

$\Sigma^{c34}_{\mu\nu}=\\{\sigma_{0}\tau_{0}\varsigma_{\nu},-\sigma_{y}\tau_{x}\varsigma_{\nu},-\sigma_{y}\tau_{y}\varsigma_{\nu},\sigma_{0}\tau_{z}\varsigma_{\nu}\\}\
,$ (5)

respectively, where $\varsigma_{\nu}$ ($\nu=0,x,y,z$) are Pauli matrices
acting in the spin subspace, $\tau_{\mu}$ ($\mu=0,x,y,z$) are Pauli matrices
acting in the valley subspace, and $\sigma_{0,x,y,z}$ are Pauli matrices
acting in the orbital subspace. With the help of U(4) symmetry, the $J$ term
in Eq. 2 can be written as a ferromagnetic coupling between the U(4) LM of
$f$-electrons and the U(4) LM of $c$-electrons Song and Bernevig (2022). When
$M\neq 0$, only the $\mu=0,z$ U(4) generators commute with the Hamiltonian,
leading to a lower U(2)$\times$U(2) symmetry group. The rotation generated by
$\mu=z$, $\nu=0$ is referred to as the valley-U(1) symmetry.

Consistent with previous results Kang and Vafek (2019); Bultinck  _et al._
(2020); Seo  _et al._ (2019); Vafek and Kang (2020); Lian  _et al._ (2021);
Xie  _et al._ (2021), a Hartree-Fock treatment of the THF model has predicted
the ground state at CNP as a U(4) LM state that respects a U(2)$\times$U(2)
subgroup Song and Bernevig (2022). The LM forms a 20-fold multiplet belonging
to the $[2,2]_{4}$ representation of the U(4) group Lian  _et al._ (2021).
These states can be approximately written as

$|\Psi_{0}\rangle=e^{-i\theta_{\mu\nu}\hat{\Sigma}_{\mu\nu}}\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}|{\rm
FS}\rangle\ ,$ (6)

where the $\rm|FS\rangle$ is the Fermi sea state with the half-filled
$c$-bands, $\hat{\Sigma}_{\mu\nu}$ are the U(4) generator operators defined by
the matrices in Eqs. 3, 4 and 5, and $\theta_{\mu\nu}$ are the rotation
parameters, and an implicit summation over repeated $\mu,\nu$ indices is
assumed. When $\theta_{\mu\nu}$’s are zero, $|\Psi_{0}\rangle$ is the valley-
polarized state because all the occupied $f$-electrons are in the $\eta=+$
valley, and the U(2)$\times$U(2) subgroup is generated by
$\hat{\Sigma}_{0,\nu}$ and $\hat{\Sigma}_{z,\nu}$ ($\nu=0,x,y,z$). For nonzero
$\theta_{\mu\nu}$’s, $|\Psi_{0}\rangle$ respects an equivalent
U(2)$\times$U(2) subgroup. The Kramers inter-valley coherent states can be
obtained by choosing nonzero $\theta_{x0}$ and $\theta_{y0}$ satisfying
$\theta_{x0}^{2}+\theta_{y0}^{2}=(\pi/4)^{2}$. When $M\neq 0$, the Kramers
inter-valley coherent states are found to have a lower energy ($\sim$0.1meV)
than the valley polarized states Lian  _et al._ (2021); Song and Bernevig
(2022).

##  III Effective Anderson model for $\nu>0$ states

###  III.1 Irrelevance of Kondo screening at CNP

Here we argue that the Kondo screening effect is irrelevant at CNP; hence, the
U(4) LM state in Eq. 6 is valid as an approximate ground state. We first
examine the energy scale of a fully symmetric Kondo state at CNP. Since the
$f$-sites are almost decoupled from each other, a reasonable approximation is
to view each $f$-site as a single Anderson impurity coupled to a bath of
$c$-electrons. If we only consider the on-site $U_{1}$ interaction and the
hybridization between $f$\- and $c$-electrons ($H^{(cf,\eta)}({\mathbf{k}})$
in Eq. 1), then it is almost a standard Anderson model with eight flavors. The
effect of $c$-bath is described by the hybridization function
$\Delta(\omega)$, defined as the imaginary part of the self-energy of a free
$f$-electron (in the absence of $U_{1}$) coupled to the $c$-bath, i.e.,
$\mathrm{Im}[\Sigma^{(f)}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}(\omega)]=\delta_{\alpha,\alpha^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}\mathrm{sgn}(\omega)\Delta(\omega)$.
The identity matrix form of $\mathrm{Im}[\Sigma_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}(\omega)]$ is guaranteed by the spin-
SU(2) ($\delta_{ss^{\prime}}$), the valley-U(1) and the time-reversal
($\delta_{\eta\eta^{\prime}}$), and crystalline
($\delta_{\alpha\alpha^{\prime}}$) symmetries. Low energy $c$-bands
(${\mathbf{k}}\to 0$) are coupled to the impurity through the constant
$\gamma$ in $H^{(cf,\eta)}({\mathbf{k}})$. Thus, they contribute to a
hybridization function $\Delta(\omega)$ proportional to the density of states
$\rho(\omega)$ of $c$-bands. In the flat-band limit ($M=0$), the linear
dispersion of $c$-bands (Fig. 1(b)) leads to a linear-in-energy density of
states as well as a linear-in-energy hybridization function, i.e.,
$\rho(\omega)\sim|\omega|$ and $\Delta(\omega)\sim|\omega|$. As a consequence,
low-lying states of the impurity will see vanishing bath electrons when the
energy scale is small enough. Both numerical Chen and Jayaprakash (1995);
Gonzalez-Buxton and Ingersent (1998); Ingersent and Si (2002) and analytical
Fritz and Vojta (2004) RG studies have shown that Anderson impurity models
with such a $\Delta(\omega)\sim|\omega|$ hybridization function do not have
the strong coupling fixed point that exhibits Kondo screening. Instead, the
only stable fixed point is the LM phase.

With a finite $M$, the $c$-bands given by $H^{(c,\eta)}({\mathbf{k}})$ in Eq.
1 have a quadratic band touching at the zero energy, i.e.,
$\pm(-M/2+\sqrt{M^{2}/4+v_{\star}^{2}{\mathbf{k}}^{2}})$, leading to finite
$\rho(0)$ and $\Delta(0)$. Nevertheless, as explained in the following, the
Kondo energy scale resulting from realistic parameters is still negligible. In
Sec. B.2 we derived an analytical expression of $\Delta(\omega)$ for the
symmetric state at CNP. For $|\omega|>M$, $\Delta(\omega)$ is almost linear in
$|\omega|$, i.e., $\Delta(\omega)\approx b\cdot|\omega|$. For $|\omega|<M$,
$\Delta(\omega)$ is a constant plus a linear term:
$\Delta(\omega)\approx\Delta(0)(1+|\omega|/M)$. Using the parameters given in
Sec. II and $M=3.697$meV, we have $b\approx 0.129$, $\Delta_{0}\approx
0.239$meV. A rough estimation of the Kondo energy scale can be made by
replacing the $\omega$-dependent $\Delta(\omega)$ with the constant
$\Delta(0)$ at $\omega=0$. Then naively applying the large-$\mathcal{N}$
formula at second order Coleman (2015), i.e., $D_{K}\approx De^{-\frac{\pi
U_{1}}{4\mathcal{N}\Delta(0)}}$ with $\mathcal{N}=8$ being the number of
flavors and $D=U_{1}/2$ the energy scale up to which the perturbation theory
applies, leads to an extremely small $D_{K}\approx U_{1}\cdot 10^{-10}$. A
better estimation is given by a poor man’s scaling that considers the
$\omega$-dependence of $\Delta(\omega)$. Readers may refer to Sec. B.2 for
more details. Here we only present the main results. There are two stages in
the RG process: (i) a stage with energy scale from $U_{1}/2$ \- scale up to
which the perturbation theory applies - to $M$. (ii) a stage with an energy
scale below $M$. RG in the first stage effectively enhances $\Delta(0)$ to
$g_{1}\Delta(0)$ with $g_{1}\approx 2.34$. Then, RG in the second stage gives
the Kondo energy scale

$\displaystyle D_{K}\approx$ $\displaystyle Me^{y}e^{-\frac{\pi
U_{1}}{4\mathcal{N}g_{1}\Delta(0)}}$ $\displaystyle\approx$ $\displaystyle
3.8\times 10^{-4}\mathrm{meV}\qquad(2M=7.394\mathrm{meV}).$ (7)

where $y\approx 1$ is a factor contributed by the $\omega$-dependence of
$\Delta(\omega)$ in the second stage. (In Sec. B.2 we have verified the
validity of the two-stage RG formula by comparing it to NRG results for the
$\mathcal{N}=4$ case.) This value is still much lower than the energy gain of
the symmetry-breaking correlated state Lian  _et al._ (2021); Bernevig  _et
al._ (2021b). The bandwidth of the Goldstone modes at CNP from $\Gamma_{M}$ to
$M_{M}$ is about 8meV. (See Fig. 2 of Ref. Bernevig  _et al._ (2021b)). If we
understand this spectrum as a tight-binding band of the Holstein–Primakoff
bosons on the $f$-sites, which form a triangular lattice, then the nearest
neighbor hopping is about 8meV/8=1meV. This hopping indicates an RKKY
interaction much larger than $k_{B}T_{K}$ evaluated in Sec. III.1.

The actual $D_{K}$ can even be much smaller than the value in Sec. III.1.
First, as $\Delta(0)\to 0$ when $M\to 0$, $T_{K}$ decays exponentially when
$M$ decreases. For example, a bandwidth $2M=5$meV corresponds to

$D_{K}\approx 5.1\times 10^{-6}\mathrm{meV}\quad(2M=5\mathrm{meV})\ .$ (8)

Second, because we only considered the $U_{1}$ interaction and the $cf$
hybridization that gives all flavors of $f$-electrons the same
$\Delta(\omega)$ (guaranteed by symmetries), the single-impurity model has a
U(8) symmetry. This U(8) symmetry must be broken when other interaction terms,
e.g., $J$ in Eq. 2, are taken into account, leading to a multiplet splitting.
When the energy scale in the RG is smaller than the multiplet splitting, the
degeneracy factor $\mathcal{N}$ should be reduced accordingly, and $D_{K}$
will be further suppressed. (One can see section 17.2 of Ref. Coleman (2015)
and Sec. B.3 for examples of how multiplet splitting suppresses $D_{K}$.)

The U(4) LM state at CNP is also supported by various experiments. In contrast
to the Kondo resonance, STM measurements have shown strong suppression of the
density of states at the zero energy at CNP Xie  _et al._ (2019); Wong  _et
al._ (2020); Nuckolls  _et al._ (2020); Choi  _et al._ (2021a, b); Oh  _et
al._ (2021); Choi  _et al._ (2019); Jiang  _et al._ (2019); Kerelsky  _et al._
(2019). Some transport experiments Lu  _et al._ (2019); Polshyn  _et al._
(2019); Stepanov  _et al._ (2020); Saito  _et al._ (2020) also exhibit a gap
behavior at CNP. Although there are also transport experiments showing
semimetal behavior, the gaplessness can be explained if there are fluctuations
of the local moments from site to site, which is possible due to the Goldstone
mode fluctuations Bernevig  _et al._ (2021b); Khalaf  _et al._ (2020) and
possible inhomogeneity of the sample.

###  III.2 Periodic Anderson model for $\nu>0$ states

We aim for an effective model describing the active excitations upon the
ground state $|\Psi_{0}\rangle$ (Eq. 6) at CNP. Let us first assume the
valley-polarized state, where $\theta_{\mu\nu}$’s in Eq. 6 are all zero such
that all the occupied $f$-electrons are in the $\eta=+$ valley. As detailed in
the supplementary material of Ref. Song and Bernevig (2022) and in Ref.
Bernevig  _et al._ (2021b), the lowest particle and hole excitations are in
the $\eta=-$ and $\eta=+$ valleys, respectively. Thus, for $\nu>0$ states,
only particle excitations in the $\eta=-$ valley will be involved, and the
electrons in the $\eta=+$ valley can be viewed as a static background. The
effective Hamiltonian can be obtained by replacing operators in the $\eta=+$
valley by their expectation values on $|\Psi_{0}\rangle$, which are $\langle
f^{\dagger}_{\mathbf{R}\alpha+s}f_{\mathbf{R}^{\prime}\alpha^{\prime}+s^{\prime}}\rangle=\delta_{\mathbf{R}\mathbf{R}^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{ss^{\prime}}$,
$\langle
c^{\dagger}_{{\mathbf{k}}a+s}c_{{\mathbf{k}}^{\prime}a^{\prime}+s^{\prime}}\rangle\approx\frac{1}{2}\delta_{{\mathbf{k}}{\mathbf{k}}^{\prime}}\delta_{aa^{\prime}}\delta_{ss^{\prime}}$,
$\langle
c^{\dagger}_{{\mathbf{k}}a+s}f_{\mathbf{R}\alpha+s^{\prime}}\rangle=0$.
Substituting these expectation values into $\hat{H}_{0}+\hat{H}_{I}$, we
obtain the effective free Hamiltonian

$\displaystyle\hat{H}_{0}^{\rm eff}=$
$\displaystyle\sum_{\begin{subarray}{c}{\mathbf{k}}aa^{\prime}s\end{subarray}}(H_{aa^{\prime}}^{(c)}({\mathbf{k}})-\mu\delta_{aa^{\prime}})c_{{\mathbf{k}}as}^{\dagger}c_{{\mathbf{k}}a^{\prime}s}-\mu\sum_{\mathbf{R}}n^{f}_{\mathbf{R}}$
$\displaystyle+\sum_{\begin{subarray}{c}{\mathbf{k}}a\alpha
s\end{subarray}}\left(e^{-\frac{1}{2}\lambda^{2}{\mathbf{k}}^{2}}H^{(cf)}_{a\alpha}({\mathbf{k}})c^{\dagger}_{{\mathbf{k}}as}f_{{\mathbf{k}}\alpha
s}+h.c.\right)\ ,$ (9)

where $n_{\mathbf{R}\alpha s}=\sum_{\alpha s}f_{\mathbf{R}\alpha
s}^{\dagger}f_{\mathbf{R}\alpha s}$ is the density operator of $f$-electrons
at $\mathbf{R}$. Here we have dropped the valley index $\eta$ as they are
limited to $\eta=-$. The $H^{(c)}({\mathbf{k}})$ and $H^{(cf)}({\mathbf{k}})$
matrices are given by the $H^{(c,-)}({\mathbf{k}})$ and
$H^{(cf,-)}({\mathbf{k}})$ matrices defined after Eq. 1. We also obtain the
effective interaction Hamiltonian

$\displaystyle\hat{H}^{\rm
eff}_{I}=\frac{U_{1}}{2}\sum_{\mathbf{R}}n^{f}_{\mathbf{R}}n^{f}_{\mathbf{R}}+\frac{U_{2}}{2}\sum_{\langle\mathbf{R}\mathbf{R}^{\prime}\rangle}n^{f}_{\mathbf{R}}n^{f}_{\mathbf{R}^{\prime}}$
$\displaystyle\quad+\frac{1}{2N_{M}}\sum_{\mathbf{q}aa^{\prime}}V(\mathbf{q})\delta
n^{c}_{-\mathbf{q},a^{\prime}}\delta
n^{c}_{\mathbf{q},a}+\frac{1}{N_{M}}\sum_{\mathbf{R}\mathbf{q}a}W_{a}e^{-i\mathbf{q}\cdot\mathbf{R}}n^{f}_{\mathbf{R}}\delta
n^{c}_{\mathbf{q}a}$
$\displaystyle\quad-\frac{J}{N_{M}}\sum_{\mathbf{R}ss^{\prime}}\sum_{{\mathbf{k}}{\mathbf{k}}^{\prime}\alpha}e^{-i({\mathbf{k}}-{\mathbf{k}}^{\prime})\cdot\mathbf{R}-\frac{\lambda^{2}({\mathbf{k}}^{2}+{\mathbf{k}}^{\prime
2})}{2}}$ $\displaystyle\quad\times(f^{\dagger}_{\mathbf{R}\alpha
s^{\prime}}f_{\mathbf{R}\alpha
s}-\frac{1}{2}\delta_{ss^{\prime}})(c^{\dagger}_{{\mathbf{k}},\alpha+2,s}c_{{\mathbf{k}}^{\prime},\alpha+2,s^{\prime}}-\frac{1}{2}\delta_{ss^{\prime}})\
,$ (10)

where $\delta
n^{c}_{\mathbf{q}a}=\sum_{s{\mathbf{k}}}(c_{{\mathbf{k}}+\mathbf{q}as}^{\dagger}c_{{\mathbf{k}}as}-\frac{1}{2}\delta_{\mathbf{q}0})$.
The $\delta n^{f}_{\mathbf{R}}$ operator in Eq. 2, which represents the
density deviation from the charge background at CNP, is now replaced by
$n^{f}_{\mathbf{R}}$, the total density in the $\eta=-$ valley, because the
charge background is compensated by the occupied $f$-electrons in the $\eta=+$
valley. The $J$ term in Eq. 2 is also simplified: As active excitations are
limited to $\eta=-$, the factor
$\eta\eta^{\prime}+(-1)^{\alpha+\alpha^{\prime}}$ becomes
$2\delta_{\alpha,\alpha^{\prime}}$.

In the flat-band limit ($M=0$), $\hat{H}_{0}^{\rm eff}+\hat{H}_{I}^{\rm eff}$
applies to arbitrary U(4) partners of the valley-polarized state, including
the so-called Kramers intervalley coherent state. To be specific, for a
generic $|\Psi_{0}\rangle$ given in Eq. 6, we can always define rotated
operators $f_{\mathbf{R}\alpha s}=Uf_{\mathbf{R}\alpha-s}U^{\dagger}$,
$c_{{\mathbf{k}}as}=Uc_{{\mathbf{k}}a-s}U^{\dagger}$, where
$U=e^{-i\theta_{\mu\nu}\hat{\Sigma}_{\mu\nu}}$ is the U(4) rotation defining
$|\Psi_{0}\rangle$, such that the effective Hamiltonian on the rotated basis
is same as Secs. III.2 and III.2.

The effective Hamiltonian $\hat{H}_{0}^{\rm eff}+\hat{H}_{I}^{\rm eff}$
respects all the crystalline symmetries discussed in Sec. II. The effective
angular momenta of the active two $f$-orbitals and four $c$-bands are $L=1$,
$-1$ and $L=1$, $-1$, $0$, $0$, respectively. And, the six-by-six
representations of $C_{3z}$, $C_{2x}$, $C_{2z}T$ on these orbitals are
$e^{-i\frac{2\pi}{3}\sigma_{z}}\oplus
e^{-i\frac{2\pi}{3}\sigma_{z}}\oplus\sigma_{0}$, $\mathbb{1}_{3\times
3}\otimes\sigma_{x}$, and $\mathbb{1}_{3\times 3}\otimes\sigma_{x}K$,
respectively, with $K$ being the complex conjugation. In the flat-band limit
($M=0$), as $|\Psi_{0}\rangle$ respects a U(2)$\times$U(2) subgroup of the
U(4) group, e.g., independent spin-charge rotations in the two valleys for the
valley polarized $|\Psi_{0}\rangle$, one may expect a U(2)$\times$U(2)
symmetry of $\hat{H}_{0}^{\rm eff}+\hat{H}_{I}^{\rm eff}$. However, since the
effective Hamiltonian only involves half of the d.o.f., e.g., the active
$\eta=-$ valley for the valley polarized $|\Psi_{0}\rangle$, only a single
U(2) factor is meaningful for $\hat{H}_{0}^{\rm eff}+\hat{H}_{I}^{\rm eff}$.
Therefore, hereafter we will say that $\hat{H}_{0}^{\rm eff}+\hat{H}_{I}^{\rm
eff}$ respects a U(2) symmetry group.

| $\hat{H}_{0}+\hat{H}_{I}$ | $\hat{H}_{0}^{\rm eff}+\hat{H}_{I}^{\rm eff}$ | $\hat{H}_{\rm SI}$
---|---|---|---
$M=0$ | U(4) | U(2) | U(2)$\times$U(2)
$M\neq 0$ | U(2)$\times$U(2) | U(2) | U(2)$\times$U(2)
Table 1: Continuous symmetries of the Hamiltonians. $\hat{H}_{0}+\hat{H}_{I}$
is the original THF model. For $\nu>0$ ($\nu<0$), $\hat{H}_{0}^{\rm
eff}+\hat{H}_{I}^{\rm eff}$ is the effective periodic Anderson model for the
active particle (hole) excitations upon the symmetry broken state at CNP.
$\hat{H}_{\rm SI}$ is a single-impurity version of $\hat{H}_{0}^{\rm
eff}+\hat{H}_{I}^{\rm eff}$.

As discussed at the end of Sec. II, when $M\neq 0$, the U(4) symmetry is
broken, and the ground state is the Kramers intervalley coherent state. As a
consequence, the effective Hamiltonian should have additional terms. In
Appendix A we show that to the order of $M^{2}$, the correction is simply an
energy shift

$\frac{M^{2}}{J}\sum_{{\mathbf{k}}}\sum_{a=3,4}\sum_{s}c_{{\mathbf{k}}as}^{\dagger}c_{{\mathbf{k}}as}+\mathcal{O}(M^{4})\
.$ (11)

Thus, the $M$-term breaks neither the crystalline symmetry nor the U(2)
symmetry, and it will play a minor role in the effective theory. To avoid
confusion, in Table 1 we summarize the continuous symmetries of different
Hamiltonians discussed in this work.

The effective model for $\nu<0$ states, which only involve hole excitations,
can be obtained by applying the particle-hole operation $\mathcal{P}_{c}$
Bernevig  _et al._ (2021a); Song and Bernevig (2022) to $\hat{H}_{0}^{\rm
eff}+\hat{H}_{I}^{\rm eff}$.

###  III.3 Single impurity model for $\nu>0$ states

Hereafter we mainly focus on a single-impurity version of $\hat{H}^{\rm
eff}_{0}+\hat{H}^{\rm eff}_{I}$, where only the correlation effects at the
$\mathbf{R}=0$ $f$-site are considered. Interactions not involving this
$f$-site will be treated at the mean-field level.

The STM spectra Nuckolls  _et al._ (2020); Choi  _et al._ (2021b, a); Oh  _et
al._ (2021) that show the zero-energy peaks also clearly show a continuity
between the gapped CNP state and the gapless states at $1\lesssim|\nu|<2$,
implying that they have the same symmetries. Therefore, in this work, we
assume no additional symmetry breaking at $1\lesssim|\nu|<2$. One should be
aware that, according to experimental observations Cao  _et al._ (2018b); Lu
_et al._ (2019); Saito  _et al._ (2020); Stepanov  _et al._ (2020); Liu  _et
al._ (2021a); Xie  _et al._ (2019); Sharpe  _et al._ (2019); Choi  _et al._
(2019); Kerelsky  _et al._ (2019); Jiang  _et al._ (2019); Wong  _et al._
(2020); Zondiner  _et al._ (2020); Choi  _et al._ (2021a) and previous
calculations Kang and Vafek (2019); Bultinck  _et al._ (2020); Seo  _et al._
(2019); Vafek and Kang (2020); Lian  _et al._ (2021); Xie  _et al._ (2021);
Xie and MacDonald (2020); Liu and Dai (2021); Cea and Guinea (2020); Venderbos
and Fernandes (2018); Ochi  _et al._ (2018); Zhang  _et al._ (2020); Liu  _et
al._ (2021b); Da Liao  _et al._ (2021); Kang and Vafek (2020); Soejima  _et
al._ (2020); Hejazi  _et al._ (2021); Da Liao  _et al._ (2019); Kennes  _et
al._ (2018); Classen  _et al._ (2019), additional symmetry breaking may happen
in $|\nu|\geq 2$ states at the zero temperature. Thus, the assumption that no
additional is broken should be invalid for $|\nu|\geq 2$ states. Nevertheless,
for the completeness of the discussion, we also use this assumption to study
the $|\nu|\geq 2$ states. As will be discussed in Secs. V and VI, these states
- even may not be ground states - will become relevant at higher temperatures.

Based on the assumption in the last paragraph, at a given filling $\nu$, the
mean fields are characterized by only a few parameters: $\nu_{f}=\langle
n^{f}_{\mathbf{R}}\rangle$, $\nu_{c,a}=\frac{1}{N_{M}}\langle\delta
n^{c}_{\mathbf{q}=0,a}\rangle$, where $\nu_{c,1}=\nu_{c,2}$,
$\nu_{c,3}=\nu_{c,4}$ due to the crystalline symmetries. The actual values of
$\nu_{f}$ and $\nu_{c,a}$ can be determined self-consistently. Then, the
considered correlated site at $\mathbf{R}=0$ is described by the Hamiltonian

$\hat{H}_{f}=-\mu_{f}\sum_{\alpha s}n_{\alpha s}^{f}+U_{1}\sum_{(\alpha
s)<(\alpha^{\prime}s^{\prime})}n^{f}_{\alpha
s}n^{f}_{\alpha^{\prime}s^{\prime}}\ ,$ (12)

where the lattice index $\mathbf{R}$ (=0) is omitted for simplicity,
$\mu_{f}=-6\nu_{f}U_{2}-\sum_{a}\nu_{c,a}W_{a}-\frac{1}{2}U_{1}+\frac{1}{2}J\nu_{c,3}+\mu$
is an effective chemical potential for the $f$-site, and $\mu$ is the global
chemical potential determined by the total filling $\nu$. The $U_{2}$,
$W_{a}$, $J$ terms in $\mu_{f}$ are contributed by the Hartree mean fields of
the $U_{2}$, $W_{a}$, and $J$ interactions in Sec. III.2, respectively. The
$U_{1}$ term in $\mu_{f}$ is from the diagonal $U_{1}$ interactions in Sec.
III.2, i.e., $\frac{1}{2}U_{1}\sum_{\alpha s}n_{\alpha s}^{f}n_{\alpha
s}^{f}=\frac{1}{2}U_{1}\sum_{\alpha s}n_{\alpha s}^{f}$. The $U_{1}$
interaction in $\hat{H}_{f}$ only contains the off-diagonal $U_{1}$
interactions of Sec. III.2.

The effective Hamiltonian of $c$-electrons is given by

$\hat{H}_{c}=\sum_{{\mathbf{k}}aa^{\prime}s}(H^{(c)}_{aa^{\prime}}({\mathbf{k}})+\Delta
H^{(c)}_{aa^{\prime}}-\mu_{c}\delta_{aa^{\prime}})c_{{\mathbf{k}}as}^{\dagger}c_{{\mathbf{k}}a^{\prime}s}\
,$ (13)

where
$H^{(c)}({\mathbf{k}})=-v_{\star}(\sigma_{x}\otimes\sigma_{0}k_{x}+\sigma_{y}\otimes\sigma_{z}k_{y})$
is the free Dirac Hamiltonian, $\Delta
H^{(c)}_{aa^{\prime}}=G\cdot\delta_{aa^{\prime}}(\delta_{a3}+\delta_{a4})$ is
a mean-field term that split the $a=1,2$ and $a=3,4$ $c$-electrons,
$\mu_{c}=-\nu_{f}W_{1}-\nu_{c}\frac{V}{\Omega_{0}}+\mu$ is an effective
chemical potential of $c$-electrons. The $W_{1}$ and $V$ terms in $\mu_{c}$
are contributed by the $W$ and $V$ interactions in Sec. III.2, respectively.
The mean field term $G$ arises from two interaction terms: (i) A mean field
treatment of the $J$ interaction in Sec. III.2 yields an energy shift
$J(\frac{1}{2}-\frac{1}{4}\nu_{f})$ for $a=3,4$ $c$-electrons. (ii) The
Hartree mean fields of the $W_{a}$ interactions in Sec. III.2 are
$\nu_{f}W_{1}$ and $\nu_{f}W_{3}$ for $a=1,2$ and $a=3,4$ $c$-electrons,
respectively. As we have absorbed $\nu_{f}W_{1}$ to $\mu_{c}$, $a=3,4$
$c$-electrons have an extra energy shift $(W_{3}-W_{1})\nu_{f}$. Combining the
two effects, the parameter $G$ is given by
$G=\frac{J}{2}+(W_{3}-W_{1}-\frac{J}{4})\nu_{f}$. Since the interaction
$V(\mathbf{q})$ of $c$-electrons is completely treated at the mean-field
level, $\hat{H}_{c}$ is an effective free-fermion system. As detailed in
Appendix A, the band structure of Eq. 13 is given by
$G/2\pm\sqrt{G^{2}/4+v_{\star}^{2}{\mathbf{k}}^{2}}$. In Fig. 1(c) we plot the
$c$-bands with $\mu_{c}=0$ and $G=10$meV. It has a gap opened by $G$, where,
according to the symmetry representations given in Sec. III.2, the lowest
conduction and highest valence band states have angular momenta $0,0$ and
$1,-1$, respectively.

The $f$-site is coupled to $c$-electrons via two terms. One is the
hybridization

$\hat{H}_{\rm hyb}=\frac{1}{\sqrt{N_{M}}}\sum_{{\mathbf{k}}a\alpha
s}\left(e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}H^{(cf)}_{a\alpha}({\mathbf{k}})c_{{\mathbf{k}}as}^{\dagger}f_{\alpha
s}+h.c.\right)\ .$ (14)

The other is the remaining ferromagnetic exchange interaction

$\displaystyle\hat{H}_{J}=$
$\displaystyle-\frac{J}{N_{M}}\sum_{\begin{subarray}{c}{\mathbf{k}},{\mathbf{k}}^{\prime}\\\
\alpha
ss^{\prime}\end{subarray}}e^{-\frac{1}{2}\lambda^{2}({\mathbf{k}}^{2}+{\mathbf{k}}^{\prime
2})}(f_{\alpha s}^{\dagger}f_{\alpha
s^{\prime}}-\frac{1}{4}\delta_{ss^{\prime}}\nu_{f})$
$\displaystyle\times(c_{{\mathbf{k}}^{\prime}\alpha+2s^{\prime}}^{\dagger}c_{{\mathbf{k}}\alpha+2s}-\frac{1}{2}\delta_{{\mathbf{k}},{\mathbf{k}}^{\prime}}\delta_{ss^{\prime}}\nu_{c,\alpha+2})\
.$ (15)

By “remaining” we mean that the mean field backgrounds $\frac{1}{4}\nu_{f}$
and $\frac{1}{2}\nu_{c,a}$ are deducted in $\hat{H}_{J}$. As explained below,
$\hat{H}_{J}$ leads to an effective Hund’s coupling that changes the symmetry
of the single-impurity model. There are also remaining density-density
interactions between $c$\- and $f$-electrons, i.e.,
$W_{a}(n^{f}-\nu_{f})(n^{c}_{\mathbf{q},a}-\nu_{c,a})$. However, these
remaining density-density interactions do not change the essence of the
single-impurity problem as $\hat{H}_{J}$ does.

Thanks to the $C_{3z}$ symmetry, $\hat{H}_{\rm hyb}$ and $\hat{H}_{J}$ couple
the $f$-electrons to two independent baths belonging to different angular
momenta. This allows us to treat the two terms separately. In a polar
coordinate, $\hat{H}_{\rm hyb}$ only couples $f$-electrons to

$\widetilde{c}_{k\alpha s}=\frac{1}{\mathcal{A}}\sum_{a}\int d\phi\cdot
H_{\alpha a}^{(cf)}({\mathbf{k}})c_{{\mathbf{k}}as}\quad(\alpha=1,2)\ ,$ (16)

where ${\mathbf{k}}=k(\cos\phi,\sin\phi)$, and $\mathcal{A}$ is a
normalization factor. Explicitly, there are $\widetilde{c}_{k1s}\sim\int
d\phi\cdot(\gamma
c_{{\mathbf{k}}1s}-v_{\star}^{\prime}ke^{i\phi}c_{{\mathbf{k}}2s})$ and
$\widetilde{c}_{k2s}\sim\int d\phi\cdot(\gamma
c_{{\mathbf{k}}2s}-v_{\star}^{\prime}ke^{-i\phi}c_{{\mathbf{k}}1s})$. Under
the $C_{3z}$ operation (defined in Sec. III.2), $\widetilde{c}_{k1s}$ and
$\widetilde{c}_{k2s}$ have effective angular momenta 1, -1, respectively. On
the other hand, $\hat{H}_{J}$ only couples $f$-electrons to

$\widetilde{c}_{kas}=\frac{1}{\mathcal{A}^{\prime}}\int d\phi\cdot
c_{{\mathbf{k}}as}\quad(a=3,4)\ ,$ (17)

where $\mathcal{A}^{\prime}$ is a normalization factor. Both
$c_{{\mathbf{k}}as}$ ($a=3,4$) have the effective angular momentum 0 under
$C_{3z}$. Because $\widetilde{c}_{kas}$ ($a=1,2$) and $\widetilde{c}_{kas}$
($a=3,4$) form different representations of $C_{3z}$, they do not couple to
each other by $\hat{H}_{c}$, hence the $\hat{H}_{\rm hyb}$-bath and the
$\hat{H}_{J}$-bath are indeed independent.

As a ferromagnetic coupling always flows to zero and becomes irrelevant in low
energy physics, we can integrate out the $\hat{H}_{J}$-bath in a single
attempt. This leads to an effective Hund’s coupling (Appendix A)

$\hat{H}_{H}=J_{H}\sum_{\alpha}n_{\alpha\uparrow}n_{\alpha\downarrow}\ ,$ (18)

where $J_{H}$, estimated as 0.34meV, is the additional energy that two
electrons will acquire if they occupy the same orbital. A nonzero $M$ does not
change the form of $\hat{H}_{H}$.

Integrating out the $\hat{H}_{\rm hyb}$-bath leads to a self-energy correction
$\Sigma^{(f)}_{\alpha s,\alpha^{\prime}s^{\prime}}(\omega)$ to the
$f$-electrons, the imaginary part of which defines the hybridization function
$\Delta(\omega)$, i.e., $\mathrm{Im}[\Sigma^{(f)}_{\alpha
s,\alpha^{\prime}s^{\prime}}(\omega)]=\delta_{\alpha\alpha^{\prime}}\delta_{ss^{\prime}}\mathrm{sgn}(\omega)\Delta(\omega)$.
The identity matrix structure of the self-energy is guaranteed by SU(2) spin
rotation symmetry and crystalline symmetries. In Appendix A we derived the
following analytical expression of $\Delta(\omega)$

$\displaystyle\Delta(\omega)=$
$\displaystyle\frac{\Omega_{0}}{4v_{\star}^{2}}\left|\omega+\mu_{c}-G\right|\left(\gamma^{2}+v_{\star}^{\prime
2}k_{F}^{2}\right)e^{-k_{F}^{2}\lambda^{2}}$
$\displaystyle\times\left[\theta(\omega+\mu_{c}-G)+\theta(-\omega-\mu_{c})\right]\
,$ (19)

where $k_{F}$ is determined by
$v_{\star}^{2}k_{F}^{2}+G^{2}/4=(\omega+\mu_{c}-G/2)^{2}$ for either
$\omega+\mu_{c}>G$ or $\omega+\mu_{c}<0$. As shown in Fig. 1(d),
$\Delta(\omega)$ has an abnormal $\omega$-dependence compared to those in
usual metals. First, $\Delta(\omega)=0$ for $\omega$ in the gap of $c$-bands
($0<\omega<G$). Second, because $f$-electrons ($L=1,-1$) have different
angular momenta as the lowest conduction $c$-bands ($L=0$), hybridization
between them vanishes as ${\mathbf{k}}\to 0$. As a result, $\Delta(\omega)$
linearly approaches zero at the conduction band edge. Third, as $f$-electrons
have the same angular momenta as the highest valence $c$-bands,
$\Delta(\omega)$ approaches a constant at the valence band edge. In summary,
around the gap $\Delta(\omega)$ has the asymptotic behaviors

$\Delta(\omega)\sim\begin{cases}|\omega+\mu_{c}-G|,\qquad&\omega+\mu_{c}\to
G+0^{+}\\\ 0,&0\leq\omega+\mu_{c}\leq G\\\
const.,&\omega+\mu_{c}\to-0^{+}\end{cases}\ .$ (20)

A nonzero $M$ does not change the asymptotic behaviors as these behaviors are
guaranteed by the $C_{3z}$ symmetry that is also respected by $M$. Due to the
damping factor $e^{-\lambda^{2}k_{F}^{2}}$, $c$-electrons with a momentum
$k_{F}\gg 1/\lambda$ will not contribute to $\Delta(\omega)$. Thus,
$\Delta(\omega)$ decays exponentially when $\omega$ exceeds
$v_{\star}/\lambda\sim 95$meV. In the rest of this work, we will restrict the
hybridization to $|\omega|<D=100$meV.

Baths giving rise to the same $\Delta(\omega)$ are physically equivalent.
Therefore, we can choose a bath that is as simple as possible. We introduce
the following effective single-impurity Hamiltonian

$\displaystyle\hat{H}_{\rm SI}=$
$\displaystyle\hat{H}_{f}+\hat{H}_{H}+\sum_{\alpha
s}\int_{-D}^{D}d\epsilon\cdot\epsilon\cdot d_{\alpha
s}^{\dagger}(\epsilon)d_{\alpha s}(\epsilon)$ $\displaystyle+$
$\displaystyle\sum_{\alpha
s}\int_{-D}^{D}d\epsilon\cdot\sqrt{\frac{\Delta(\epsilon)}{\pi}}(f_{\alpha
s}^{\dagger}d_{\alpha s}(\epsilon)+h.c.)\ ,$ (21)

where $\hat{H}_{f}$ and $\hat{H}_{H}$ are given by Eqs. 12 and 18,
respectively, and $d_{\alpha s}(\epsilon)$, satisfying
$\\{d_{\alpha^{\prime}s^{\prime}}(\epsilon^{\prime}),d_{\alpha
s}^{\dagger}(\epsilon)\\}=\delta_{\alpha^{\prime}\alpha}\delta_{s^{\prime}s}\delta(\epsilon^{\prime}-\epsilon)$,
are the auxiliary bath fermions introduced to reproduce the hybridization
function. $\hat{H}_{\rm SI}$ is completely defined by the four parameters:
$\mu_{f}$ the chemical potential of $f$-electrons, $U_{1}$ the Coulomb
repulsion, $J_{H}$ the Hund’s coupling, $\mu_{c}$ the chemical potential of
$c$-electrons, and $G$ the gap of $c$-bands (Fig. 1(c)). As explained at the
beginning of this subsection, the actual values of $\mu_{f}$, $\mu_{c}$, and
$G$ depend on the occupations $\nu_{f}$, $\nu_{c,a}$, which can be determined
by self-consistent calculations at given total fillings $\nu$. We plot the
self-consistent $\mu_{f}$, $\mu_{c}$, and $G$ as functions of $\nu$ in Fig.
2(a). For $\nu$ changing from 0 to 4, $\mu_{c}$ changes from 0.57 to 65.80meV,
$\mu_{f}$ changes from -24.42meV to 128.36meV , and $G$ changes from 8.47meV
to 14.27meV. We can now regard $\mu_{c},\mu_{f},G$ as given parameters that
define the single-impurity problem.

It is worth mentioning that Sec. III.3 has an emergent U(2)$\times$U(2)
symmetry - independent spin-charge rotations in the $\alpha=1,2$ orbitals -
that is higher than the U(2) symmetry of Secs. III.2 and III.2. It is not
surprising that a single-impurity model has a higher symmetry than its lattice
version. If $J=0$, there would be no Hund’s coupling $J_{H}$, and
$\hat{H}_{\rm SI}$ would have a U(4) symmetry, as expected in a four-flavor
Anderson impurity model without multiplet splitting.

##  IV Phase diagram of the single-impurity model

###  IV.1 Poor man’s scaling

Before going to numerical calculations, we first apply a poor man’s scaling to
the single impurity model Sec. III.3. The scaling theory helps us understand
several features of the data obtained by NRG, as will be discussed in the
following subsections. And, it predicts the Kondo energy scales in the same
order as those given by NRG.

We assume that the ground state of the isolated impurity has $n_{f}$ (integer)
occupied $f$-electrons. One should not confuse $n_{f}$ with $\nu_{f}$ \- the
expectation value of $f$-occupation after the impurity is coupled to the bath.
The chemical potential $\mu_{f}$ must be in the range
$(n_{f}-1)U_{1}<\mu_{f}<n_{f}U_{1}$. We apply a Schrieffer-Wolff
transformation to Sec. III.3 to obtain an effective Coqblin–Schrieff model
where the local Hilbert space of $f$-electrons is restricted to $n_{f}$
particles. The transformation involves virtual particle and hole excitations,
the energies of which are $\Delta E_{+}=n_{f}U_{1}-\mu_{f}$ and $\Delta
E_{-}=\mu_{f}-(n_{f}-1)U_{1}$, respectively. (We have ignored the $J_{H}$ term
in $\Delta E_{\pm}$ as it is small compared to $U_{1}$.) Adding the two
contributions, we have

$\displaystyle\hat{H}=\hat{H}_{H}+\sum_{\alpha
s}\int_{-D}^{D}d\epsilon\epsilon d_{\alpha s}^{\dagger}(\epsilon)d_{\alpha
s}(\epsilon)+\frac{4g}{\pi
U_{1}}\sum_{\alpha\alpha^{\prime}ss^{\prime}}\int_{-D}^{D}d\epsilon
d\epsilon^{\prime}\bigg{[}$
$\displaystyle\;\times\sqrt{{\Delta(\epsilon)\Delta(\epsilon^{\prime})}}(f_{\alpha
s}^{\dagger}f_{\alpha^{\prime}s^{\prime}}-x\delta_{\alpha\alpha^{\prime}}\delta_{ss^{\prime}})d_{\alpha^{\prime}s^{\prime}}^{\dagger}(\epsilon^{\prime})d_{\alpha
s}(\epsilon)\bigg{]}.$ (22)

The parameters $g,x$ are given by

$g=\frac{U_{1}}{4}\left(\frac{1}{\Delta E_{+}}+\frac{1}{\Delta
E_{-}}\right),\qquad x=\frac{\Delta E_{-}}{U_{1}}\ .$ (23)

$g$ is a dimensionless parameter characterizing the anti-ferromagnetic
coupling strength between the impurity and the bath. $x$ appears as a “charge
background” of the $f$-electrons. For $\mu_{f}=(n_{f}-\frac{1}{2})U_{1}$,
there is $g=1$, $x=\frac{1}{2}$. For a generic
$(n_{f}-1)U_{1}<\mu_{f}<n_{f}U_{1}$, $g\geq 1$ and $0<x<1$. Flow equations of
$g,x$ are derived in Sec. B.3, where the divergence of $g$ indicates the
strong coupling fixed point that exhibits the Kondo screening. We notice that
$x$ always flows to $n_{f}/4$, i.e., the occupation fraction of $f$-electrons.

One should be careful about the cutoff $D$ in Sec. IV.1. First, it must be
smaller than $\Delta E_{+}$ and $\Delta E_{-}$ for the Schrieffer-Wollf
transformation to be valid. Second, for analytical conveniences, we only keep
the positive branch of $\Delta(\omega)$ (Eq. 20) at $\omega>G-\mu_{c}$ because
when $\nu>0$ the negative branch is far away from the Fermi level (Fig. 1(d)).
Hence, we also require $D<\mu_{c}-G$. We can choose $D=\min(\mu_{c}-G,\Delta
E_{+},\Delta E_{-})$.

We first consider the case $n_{f}=1$. The flow equation of $g(t)$ as the
cutoff is successively reduced to $De^{-t}$ is given by

$\frac{dg}{dt}=\frac{4\Delta(0)}{\pi U}\mathcal{N}g^{2}+\mathcal{O}(e^{-t})\
,$ (24)

and the initial condition $g(0)$ is given by Eq. 23. Here $\mathcal{N}=4$ is
the number of flavors. The local Hilbert space for $n_{f}=1$ is four-fold. The
Hund’s coupling $J_{H}$ does not split the four-fold degeneracy and hence does
not appear in the flow equation. The $\mathcal{O}(e^{-t})$ terms are
irrelevant at small energy scales but they may affect the coupling constant at
an early stage of the RG process. As shown in Fig. 1(d), the positive branch
of $\Delta(\omega)$ can be well approximated by a linear function, i.e.,
$\Delta(\omega)\approx\Delta(0)(1+\omega/(\mu_{c}-G))$. Using this linear
approximation we obtain the Kondo energy scale (Sec. B.3)

$D_{K}^{(1)}=De^{y_{1}}\cdot e^{-\frac{\pi U_{1}}{4\mathcal{N}\Delta(0)g(0)}}$
(25)

where $y_{1}\approx-0.75\frac{D}{\mu_{c}-G}<0$ is factor contributed by the
irrelevant $\mathcal{O}(e^{-t})$ terms and will slightly suppress the Kondo
screening. The suppression is due to the fact that, for the $n_{f}=1$ states,
virtual processes contributing to the RG equation involve more hole
excitations (with smaller $\Delta(\omega<0)$) than particle excitations (with
larger $\Delta(\omega>0)$) in the bath. As a result, the resulting
$k_{B}T_{K}$ is smaller than the standard case ($y_{1}=0$) where the coupling
is a constant, i.e., $\Delta(\omega)=\Delta(0)$.

We then study the RG equations at $n_{f}=2$. Unlike the $n_{f}=1$ case, Hund’s
coupling $J_{H}$ splits the six-dimensional local Hilbert space. According to
Eq. 18, the four states with
$(n^{f}_{1\uparrow},n^{f}_{1\downarrow};n^{f}_{2\uparrow},n^{f}_{2\downarrow})=$(10;10),
(10;01), (01;10), (01;01) do not feel $J_{H}$ and have the energy
$-2\mu_{f}+U_{1}$, whereas the two states (11;00), (00;11) have the energy
$-2\mu_{f}+U_{1}+J_{H}$. We divide the RG process into two stages: (i) a stage
with an energy scale from $D$ to $J_{H}$, (ii) a stage with an energy scale
below $J_{H}$. In the first stage, $J_{H}$ plays a minor role; hence, a flow
equation similar to Eq. 24 applies. If $g$ diverges before the energy scale
reaches $J_{H}$, then the Kondo scale is given by

$D_{K}^{(2)\prime}=D\cdot e^{-\frac{\pi U_{1}}{4\mathcal{N}\Delta(0)g(0)}}$
(26)

there is no $y$-factor as in Eq. 25 because, for the $n_{f}=2$ states, the
virtual processes that contribute to the RG equations evolve equal holes and
particles in the bath. Otherwise, $g$ will be renormalized to

$g_{1}=\frac{g(0)}{1-g(0)\frac{4\Delta(0)}{\pi
U_{1}}\mathcal{N}\ln\frac{D}{J_{H}}}$ (27)

at the energy scale of $J_{H}$. As detailed in Sec. B.3, RG in the second
stage is similar to Eq. 24 except that, due to the multiplet splitting, the
factor $\mathcal{N}=4$ is replaced by 2. This leads to the Kondo energy scale

$D_{K}^{(2)\prime\prime}=J_{H}\cdot e^{-\frac{\pi
U_{1}}{8\Delta(0)g_{1}}}=D\left(\frac{D}{J_{H}}\right)^{\frac{\mathcal{N}}{2}-1}e^{-\frac{\pi
U_{1}}{8\Delta(0)g(0)}}\ .$ (28)

Considering the $g$ may diverge in either the first or the second stage, the
physical Kondo energy scale can be written as

$D_{K}^{(2)}=\begin{cases}k_{B}D_{K}^{(2)\prime},\qquad&k_{B}D_{K}^{(2)\prime}>J_{H}\\\
k_{B}D_{K}^{(2)\prime\prime},\qquad&\text{otherwise}\end{cases}\ .$ (29)

$\nu$ |  $D_{K}$ (meV) |  $\delta_{K}$(meV) |  $k_{B}T_{K}$ (meV)
---|---|---|---
0.75 | 0.292 | 0.168 | 0.282
1.00 | 0.259 | 0.143 | 0.244
1.25 | 0.362 | 0.236 | 0.407
1.50 | 0.740 | 0.979 | 0.796
1.75 | 1.52 | 1.70 | 1.38
2.00 | - | 1.74 | 1.35
2.25 | 2.66 | 1.68 | 1.28
2.50 | 3.21 | 1.76 | 1.48
2.75 | 3.99 | 2.35 | 1.92
3.00 | 4.36 | 2.43 | 2.60
3.25 | 4.15 | 3.64 | 3.50
3.50 | 3.48 | 2.88 | 4.03
Table 2:  Kondo energy scales at various fillings $\nu$. $D_{K}$,
$\delta_{K}$, and $T_{K}$ are Kondo energy scales estimated by the poor man’s
scaling, the NRG spectral density, and the NRG spin susceptibility,
respectively. $\delta_{K}$ is defined as the half width at half maximum of the
spectral peak. $T_{K}$ is defined as the turning temperature where $\chi(T)$
transitions from the Curie-Weiss behavior to the Fermi liquid behavior. There
is no data of $D_{K}$ at $\nu=2$ because $\nu=2$ is close to a mixed-valence
case where $\mu_{f}\approx U_{1}$ and hence the Schrieffer-Wolff
transformation does not apply.

The theory at $n_{f}=3$ is almost the same as the theory at $n_{f}=1$ except
that the four single-particle states are now replaced by the four single-hole
states. Thus, the local Hilbert space is also four-dimensional and will not be
split by Hund’s coupling $J_{H}$. The Kondo energy scale is given by

$D_{K}^{(3)}=De^{y_{3}}\cdot e^{-\frac{\pi U_{1}}{4\mathcal{N}\Delta(0)g(0)}}$
(30)

where $y_{3}\approx 0.75\frac{D}{\mu_{c}-G}>0$ is a factor contributed by the
irrelevant $\mathcal{O}(e^{-t})$ terms and will slightly enhance the Kondo
screening. The enhancement arises from the fact that, in contrast to the case
of $n_{f}=1$, for $n_{f}=3$ states the virtual processes contributing to the
RG equation evolve more particle excitations (with larger $\Delta(\omega>0)$)
than hole excitations (with smaller $\Delta(\omega<0)$) in the bath.

In Table 2 we tabulate the Kondo energy scales obtained by the poor man’s
scaling at various fillings using the filling-dependent $\mu_{c},\mu_{f},G$
parameters given in Fig. 2(a), and compare them to the values obtained by the
NRG method.

###  IV.2 The NRG method

Figure 2: Phase diagram and fixed points. (a) The self-consistent mean-field
values of $\mu_{c},\mu_{f},G$ as functions of the total filling $\nu$ from
$\nu=0$ to 4, where we have enforced the symmetries of the correlated
insulator state at CNP. The left and right vertical axis coordinates are for
the $\mu_{c}$, $\mu_{f}$ and $G$ parameters, respectively. (b) The phase
diagram in the parameter space of $\mu_{c},\mu_{f}$ for $G=8$meV. The white
lines are phase boundaries between local moment (LM) phases and the strong
coupling phase. The dashed black lines are crossover boundaries between the
frozen impurity (FI) and Kondo regimes of the strong coupling phase. The color
maps the half-width of the spectral density peak, reflecting the Kondo energy
scale if in the Kondo regime. The solid black line indicates the trajectory of
$\mu_{c}$ and $\mu_{f}$ determined from a self-consistent calculation as $\nu$
changes from 0 to 4, where the five arrows from left to right represent
$\nu=0,1,2,3,4$, respectively. (c) is the same as (b) but a different
parameter $G=14$meV is used. (d), (e) The RG flows of the many-body spectra of
the scaled Hamiltonian $\widetilde{H}_{N}$ ($N\in$ odd) towards fixed points
in the Kondo regime. (d) uses mean-field parameters at $\nu=1.25$:
$\mu_{c}=34.4$meV, $\mu_{f}/U_{1}=0.536$, $G=9.86$meV, and (e) uses mean-field
parameters at $\nu=2.5$: $\mu_{c}=49.9$meV, $\mu_{f}/U_{1}=1.31$, $G=11.8$meV.
The levels’ colors represent the symmetry sectors labeled by good quantum
numbers $(Q_{1},Q_{2};S_{1},S_{2})$. Since the levels in sectors
$(Q_{1},Q_{2};S_{1},S_{2})$ and $(Q_{2},Q_{1};S_{2},S_{1})$ are identical,
only $|Q_{1}|\geq|Q_{2}|$ sectors are shown for simplicity. The insets are the
resulting single-particle spectral densities that exhibit a Kondo resonance.
(f), (g), (h) The RG flows of many-body spectra of the scaled Hamiltonian
$\widetilde{H}_{N}$ ($N\in$ even) in the LM1,2,3 phases, where $\mu_{c}=5$meV,
$G=$8meV, and $\mu_{f}/U_{1}=0.5,1.5,2.5$, respectively. The insets are the
resulting single-particle spectral densities that exhibit Hubbard bands.

In the NRG method Krishna-murthy  _et al._ (1980a, b); Bulla  _et al._ (2008),
the bath is alternatively realized by a Wilson chain constructed in the way
described below. First, the energy window $[-D,D]$ is discretized on a
logarithmic scale, i.e., $\omega_{n}=D/\Lambda^{n-1}$ ($n=1,2\cdots$), where
$\Lambda>1$ is a scaling factor (chosen as 3 in this work). Then for each
energy shell $\omega_{n}\leq|\omega|<\omega_{n+1}$ two auxiliary bath
electrons corresponding to the positive and negative energies are introduced
to reproduce the corresponding $\Delta(\omega)$. These auxiliary electrons are
further recombined into a Wilson chain, $d_{n\alpha s}$, such that (i) only
the first site $d_{1\alpha s}$ couples to the impurity, (ii) the chain has
only on-site and nearest neighbor hopping terms. The Hamiltonian Sec. III.3 is
now mapped to an impurity plus a Wilson chain

$\displaystyle\hat{H}_{N}=\hat{H}_{f}+\hat{H}_{H}+\sum_{\alpha
s}t_{0}(f_{\alpha s}^{\dagger}d_{1\alpha s}+h.c.)$
$\displaystyle+\sum_{n=1}^{N}\sum_{\alpha s}\epsilon_{n}d_{n\alpha
s}^{\dagger}d_{n\alpha s}+\sum_{n=1}^{N-1}\sum_{\alpha s}(t_{n}d_{n+1\alpha
s}^{\dagger}d_{n\alpha s}+h.c.)\ ,$ (31)

where $N$ is an arbitrarily large number. The parameters $\epsilon_{n}$ and
$t_{n}$ can be computed from $\Delta(\omega)$ using a standard iterative
algorithm Bulla  _et al._ (2008). For $n\to\infty$,
$t_{n}\to\frac{1}{2}(1+\Lambda^{-1})\Lambda^{-\frac{1}{2}(n-1)}$ and
$\epsilon_{n}\sim\Lambda^{-n}$ is negligible compared to $t_{n}$. Thus, the
right-most sites represent the lowest-lying bath states.

One can define the $N$th scaled Hamiltonian as
$\widetilde{H}_{N}=(\Lambda)^{\frac{1}{2}N-1}\hat{H}_{N}$. They can be
constructed iteratively

$\displaystyle\widetilde{H}_{N+1}=$
$\displaystyle\Lambda^{\frac{1}{2}}\widetilde{H}_{N}+\Lambda^{\frac{1}{2}(N-1)}\sum_{\alpha
s}\big{(}\epsilon_{N+1}d_{N+1,\alpha s}^{\dagger}d_{N+1,\alpha s}$
$\displaystyle+t_{N}d_{N+1,\alpha s}^{\dagger}d_{N,\alpha s}+t_{N}d_{N,\alpha
s}^{\dagger}d_{N+1,\alpha s}\big{)}\ .$ (32)

The Hilbert space dimension increases exponentially in this iterative process.
The NRG algorithm truncates the Hilbert space by keeping a fixed number
(chosen to be $\sim$1200 in this work) of the lowest-lying states at each
step. In order to keep the symmetry in the truncated Hilbert space, in
practice we keep all the states up to the lowest gap above the 1200th state.

###  IV.3 Phase diagram and fixed points

Two successive transformations (Sec. IV.2) that take $\widetilde{H}_{N}$ to
$\widetilde{H}_{N+2}$ can be thought as a renormalization group operation
Krishna-murthy  _et al._ (1980a, b). The system is said to achieve a fixed
point when $\widetilde{H}_{N}$ and $\widetilde{H}_{N+2}$ have the same low-
lying many-body spectrum. We can obtain a zero temperature phase diagram in
the parameter space of $\mu_{c},\mu_{f},G$ by analyzing the fixed points. For
the completeness of discussions, here we let $\mu_{f}$ take value in
$[-0.5U_{1},3.5U_{1}]$ such that the corresponding impurity occupation (in the
decoupled limit) $n_{f}=\mu_{f}/U_{1}+1/2$ takes value in $[0,4]$. In Fig.
2(b) and (c) we show the obtained phase diagrams in the parameter space of
$\mu_{c},\mu_{f}$ for $G=8$meV and $G=14$meV, respectively. 8meV and 14meV are
chosen to be close to the minimal (8.47meV) and maximal (14.27meV) values of
$G$ (Fig. 2(a)), respectively.

Due to the U(2)$\times$U(2) symmetry, all the many-body levels can be
classified into symmetry sectors labeled by the good quantum numbers
$(Q_{1},Q_{2};S_{1},S_{2})$, where $Q_{i}$ and $S_{i}$ are the charge and spin
of the $i$th U(2) symmetry, respectively. Here we take the convention that
$Q_{1}+Q_{2}=0$ corresponds to a total occupation $2N+2$ ($2N$) for odd (even)
$N$. A fixed point is characterized by low energy many-body levels and the
associated quantum numbers. In the whole phase diagram, we find two distinct
types of stable fixed points: (i) the strong coupling fixed point exhibiting a
Fermi liquid behavior and (ii) the LM fixed points exhibiting nonzero spin
momenta. At a strong coupling fixed point, as exampled in Fig. 2(d), (e), for
either even or odd $N$, the ground state is a singlet and has
$(Q_{1},Q_{2};S_{1},S_{2})=(2k,2k;0,0)$ for some order one integer $k$, which
in most cases equals to 0. The low-lying many-body spectrum is identical to
the one of a free-fermion chain defined by $\epsilon_{n}$ and $t_{n}$ with an
additional chemical potential term. In other words, the impurity acts as if it
was nonexistent. The underlying mechanism is either the Kondo screening, where
the impurity is an effective LM screened by the bath, or the impurity
freezing, where the impurity occupation $\nu_{f}$ is effectively empty or
full. We refer to the two cases as the Kondo regime and the frozen impurity
regime, respectively. The fixed points shown in Fig. 2(d), (e) are in the
Kondo regime because, if we continuously change $\mu_{c}$ to 0, they evolve to
the LM states with one and two particles, respectively. There is a crossover
between Kondo and frozen impurity regimes as one changes $\mu_{f}$, as
indicated by the dashed lines in Fig. 2(b), (c). Later we will determine the
crossover boundary using the spectral density.

At an LM fixed point, as exampled in Fig. 2(f), the low-lying many-body
spectrum is identical to a free-fermion chain plus a detached LM. Depending on
the representation of the ground state, the LM fixed points can be further
classified into LMn, where $n=1,2,3$ is the effective impurity occupation. The
flows of the spectra towards these fixed points are shown in Fig. 2(f), (g),
(h), respectively. LMn ground states have the same SU(2)$\times$SU(2)
representations as ground states of $\hat{H}_{f}+\hat{H}_{H}$ (Eqs. 12 and 18)
with $n$ impurity electrons, where the Hubbard interaction freezes charge
excitations and the Hund’s coupling prefers states with electrons lying in
different orbitals. For $n=1$, the ground states are four-fold degenerate and
belong to the symmetry sectors
$(Q_{1},Q_{2};S_{1},S_{2})=(2k+1,2k;\frac{1}{2},0)$ and
$(2k,2k+1;0,\frac{1}{2})$, corresponding to the spin-$\frac{1}{2}$ states of
the two SU(2)’s, respectively. The SU(2) representations are the same as those
of the four single-particle states of $\hat{H}_{f}+\hat{H}_{H}$:
$(n^{f}_{1\uparrow},n^{f}_{1\downarrow};n^{f}_{2\uparrow},n^{f}_{2\downarrow})=$
(10;00), (01;00), (00;10), (00;01). For $n=2$, the ground states are also
four-fold degenerate but belong to a different symmetry sector
$(2k+1,2k+1;\frac{1}{2},\frac{1}{2})$. One can understand them as the product
states of two spin-$\frac{1}{2}$ states of the two SU(2)’s. They have the same
SU(2) representations as the four two-particle states of
$\hat{H}_{f}+\hat{H}_{H}$: (10;10), (10;01), (01;10), (01;01). Without Hund’s
coupling $J_{H}$, the LM2 ground states would be $\binom{4}{2}=6$-fold
degenerate. A finite $J_{H}$ raises the energy of the two many-body states
with two electrons occupying the same orbital, i.e., (11;00), (00;11). For
$n=3$, the ground states are still four-fold degenerate but belong to the
symmetry sectors $(2k-1,2k;\frac{1}{2},0)$ and $(2k,2k-1;0,\frac{1}{2})$.
Their SU(2) representations are same as the four single-hole states
$\hat{H}_{f}+\hat{H}_{H}$: (01;11), (10;11), (11;01), (11;10).

It is also helpful to look at the global U(2) symmetry representations, where
the two U(2) rotations are the same. The total charge and spin of LM1,2,3 are
1, 2, 3 (mod 4) and $\frac{1}{2}$, $\frac{1}{2}\pm\frac{1}{2}$, $\frac{1}{2}$,
respectively.

Phase boundaries between different LM phases and the strong coupling phase are
described by unstable fixed points where different ground states cross with
each other. The phase boundaries are shown by the white lines in Fig. 2(b),
(c). Starting from an LMn phase, increasing $\mu_{c}$ will eventually drive it
into a strong coupling phase due to the enhancement of hybridization. The
critical $\mu_{c}$, as expected, is close to $G$, the conduction band edge
(Fig. 1(c), (d)).

Figure 3: Spectral densities $A(\omega,T)$ at various fillings $\nu$ and
temperatures. (a) Spectral densities at the zero temperature for $\nu=0$, 0.1,
$0.2\cdots 4$. The curves are offset by $0.1\nu$ (meV-1) for clarity. (b) is
the same as (a) but is shown with a smaller vertical scale for clarity of
Hubbard bands, marked by inverted triangles. The curves are offset by
$0.01\nu$ (meV-1). (c) and (d) are spectral densities at finite temperatures,
where the curves are offset by $0.01\nu$ (meV-1).

##  V Spectral density

We calculate the spectral density of the $f$-electrons,
$A(\omega,T)=-\frac{1}{\pi}\sum_{\alpha s}G_{\alpha s}(\omega,T)$, with
$G_{\alpha s}(\omega,T)$ being the retarded Green’s function of $f_{\alpha s}$
at the temperature $T$. We use the method described in Ref. Bulla  _et al._
(2001) to collect the many-body levels at different RG steps to compute
$A(\omega,T)$. The fixed points in Fig. 2(d), (e) are in the Kondo regime and
hence have sharp zero-energy peaks due to the Kondo resonance. The fixed
points in Fig. 2(f), (g), (h) are in the LM1,2,3 phases, respectively, thus,
their spectral density are dominated by the upper and lower Hubbard bands. We
compute the spectral densities for all the points in the phase diagrams in
Fig. 2(b), (c). For every point we identify a central peak and measure its
half-width $\delta_{K}$ at half maximum. (If there is no central peak, e.g.,
Fig. 2(f), $\delta_{K}=0$.) $\delta_{K}$ is indicated by the color in Fig.
2(b), (c). We can distinguish the Kondo and frozen impurity regimes in the
strong coupling phase through spectral density. Intuitively, a state in the
Kondo regime should have a Kondo resonance. By contrast, a state in the frozen
impurity regime should have its main spectral weight away from zero energy
because the impurity occupation is either empty or full. Thus, we identify a
phase point in the Kondo regime if $\delta_{K}$ covers the zero energy and in
the frozen impurity regime otherwise. The crossover boundary between the two
regimes is indicated by dashed lines in Fig. 2(b), (c).

Several features of $\delta_{K}$ in Fig. 2(b), (c) can be understood using the
poor man’s scaling developed in Sec. IV.1. First, there are three domes around
$\mu_{f}=\frac{1}{2}U_{1},\frac{3}{2}U_{1},\frac{5}{2}U_{1}$ where
$\delta_{K}$ is relatively small. They correspond to the $n_{f}=1,2,3$ cases
discussed in Sec. IV.1. From the poor man’s scaling perspective, these three
$\mu_{f}$’s correspond to the minimal initial value of the coupling constant
$g$ (Eq. 23), which then leads to smaller $D_{K}$’s. Second, when $\mu_{c}$ is
small ($\lesssim 30$meV), $\delta_{K}$ in the middle dome is smallest. The
reason is that the Kondo energy scale $D_{K}$ for $n_{f}=2$ will be strongly
suppressed due to the multiplet splitting if $D_{K}$ is smaller than $J_{H}$.
One can see that the $\mathcal{N}$ factor in the exponential function of Eq.
28 is replaced by 2. Third, for the same $\mu_{c}$, the first dome has lower
$\delta_{K}$ than the third dome. This difference is a result of the
suppression factor $y_{1}$ for $n_{f}=1$ (Eq. 25) and the enhancement factor
$y_{3}$ for $n_{f}=3$ (Eq. 30) due to the particle-hole asymmetry of
$\Delta(\omega)$, as discussed in detail in Sec. IV.1.

In order to compare our results with STM measurements, we need to adopt
physical $\mu_{c}$, $\mu_{f}$, $G$ parameters. As discussed in Sec. III.3,
mean-field parameters $\mu_{c}$, $\mu_{f}$, $G$ can be determined as functions
of the filling $\nu$ via a self-consistent calculation of Sec. III.2. Results
are shown in Fig. 2(a). The resulting spectral densities at the zero
temperature and fillings from 0 to 4 are shown in Fig. 3(a), (b). For $\nu=0$,
the state is in the frozen impurity regime with an (almost) zero occupation;
hence the spectral weight is mainly distributed at positive energy. As $\nu$
increases, the spectral peak moves to the zero energy and is eventually pinned
at the zero energy to form a Kondo resonance. This behavior precisely matches
STM experiments at low temperatures ($T<1$K) Nuckolls  _et al._ (2020); Choi
_et al._ (2021b, a); Oh  _et al._ (2021). One can also observe the evolution
of Hubbard bands at $T=0$ as $\nu$ changes (Fig. 3(b)), but they are
relatively weak compared to the Kondo resonance peaks. At finite temperatures
(Fig. 3(c), (d)), the Kondo resonance peaks are smeared by thermal
fluctuations and the evolution of Hubbard bands becomes clearer. As $\nu$
increases from 0 to 4, the Hubbard bands periodically pass through the zero
energy, matching the cascade of transitions seen in STM experiments at higher
temperatures Wong  _et al._ (2020); Choi  _et al._ (2021a).

As explained in Sec. III.3, we have assumed in the current calculations that
all the $\nu\neq 0$ states have the same symmetries as the correlated
insulator at CNP. Thus, our current calculation does not reproduce the
correlated gaps around $|\nu|\approx 2$ Nuckolls  _et al._ (2020); Choi  _et
al._ (2021b, a); Oh  _et al._ (2021) that arise from further symmetry
breaking, which are estimated as 0.35-0.37meV Lu  _et al._ (2019).
Nevertheless, our results should apply to higher temperature states where the
symmetry is recovered by thermal fluctuations.

One can use $\delta_{K}$ to estimate the Kondo energy scale. Using the
$\nu$-dependent $\mu_{c},\mu_{f},G$ parameters given in Fig. 2(a), we tabulate
the $\delta_{K}$’s at different fillings in Table 2. Comparing it to the Kondo
energy scale ($D_{K}$) estimated by the poor man’s scaling, we find $D_{K}$ is
about $\frac{1}{2}\delta_{K}\sim 2\delta_{K}$.

##  VI Local moments and the Pomeranchuk effect

At a temperature exceeding the Kondo energy scale, the LM will become
effectively decoupled from the bath and visible in experimental measurements.
This is the mechanism of the Pomeranchuk effect Rozen  _et al._ (2021); Saito
_et al._ (2021b). Refs. Rozen  _et al._ (2021) observed a higher entropy
($\sim 1k_{B}$ per moiré cell with $k_{B}$ being the Boltzmann’s constant)
state at $\nu\approx 1$ at the temperature $T\approx 10$K. As this entropy can
be quenched by an in-plane magnetic field, it is ascribed to a free local
moment. Ref. Saito  _et al._ (2021b) observed a similar effect at
$\nu\approx-1$ and showed that an additional resistivity peak that is absent
at $T=0$ develops in the higher entropy state at $T\approx 10$K. These
observations can be naturally explained by the transition from the Fermi
liquid phase to the LM phase as the temperature increases.

Figure 4: Spin susceptibility and entropy contributed by the impurity. (a)
$\chi_{\rm loc}(T)/\chi_{\rm loc}(0)$ as a function of filling $\nu$ and
temperature $T$. The black curve indicates $T_{K}$ defined by $\chi_{\rm
loc}(T_{K})/\chi_{\rm loc}(0)=0.5$. (b) The local spin susceptibilities
$\chi_{\rm loc}(T)$ at fillings $\nu=0,0.5,1\cdots 4$. (c) The entropy
contributed by the impurity as a function of $\nu$ and $T$. (d) The entropy
contributed by the impurity $S_{\rm imp}(T)/(k_{B}\ln 2)$ at fillings
$\nu=0,0.5,1\cdots 4$. (e) The spectral densities at $\nu=1$ at various
temperatures. (f) The entropy contributed by the impurity as a function of
$\nu$ at $B=0,5,10,15,20$ T and temperature $k_{B}T=0.82$ meV.

The LM phase around $\nu\approx 1$ already manifests itself in the RG flow in
Fig. 2(d). At an early stage of the NRG calculation ($N<7$), the ground states
belong to the symmetry sectors $(S_{1},S_{2})=(\frac{1}{2},0)$,
$(0,\frac{1}{2})$, which, according to discussions in Sec. IV.3, form the LM1
phase. Here $S_{i}$ is the total spin of the $i$th SU(2) symmetry group. The
ground state changes into the singlet state (0,0) for $N\geq 7$. If we regard
the energy scales in NRG as an effective temperature, then the RG flow implies
that the ground state will change from the strong coupling phase to the LM1
phase as temperature increases.

To further demonstrate the LM phase, we calculate the local spin
susceptibilities $\chi_{\rm loc}(T)$ using the filling-dependent
$\mu_{c},\mu_{f},G$ parameters given in Fig. 2(a). $\chi_{\rm loc}(T)$ is
defined as $\frac{dM_{\rm loc}}{dB_{\rm loc}}$ Bulla  _et al._ (2005); Fang
_et al._ (2015), with $M_{\rm loc}$ being the spin momenta contributed by the
impurity and $B_{\rm loc}$ a local magnetic field that only acts on the
impurity but not the bath. (For $B<20$T, the Zeeman splitting
$\mu_{B}B<1.2$meV is much smaller than $\mu_{c}$, hence its action on bath is
negligible.) As shown in Fig. 4(a), (b), for $\nu\geq 0.5$, $\chi_{\rm
loc}(T)$ approaches a constant as $T\to 0$, and obeys the Curie’s law, i.e.,
$\chi_{\rm loc}(T)\sim T^{-1}$, when $T$ is larger than the Kondo energy
scale. Obeying Curie’s law is a clear indication of a free LM. (One may notice
that the $T$-dependences of $\chi_{\rm loc}(T)$’s for $\nu<0.5$ are non-
monotonous. The $\nu<0.5$ states lie in the frozen impurity regime, thus the
spin susceptibilities are extremely small when $T\to 0$, and will start to
increase when $T$ is able to excite the LM1 states.) For $\nu>0.5$, we can
define the Kondo temperature $T_{K}$ as the turning point between Curie’s
behavior and Fermi liquid behavior. Specifically, it can be obtained as the
temperature where $\chi_{\rm loc}(T)=0.5\chi_{\rm loc}(0)$. We tabulate the
resulting $T_{K}$ in Table 2. As shown in the table, such defined $T_{K}$ is
about $\frac{1}{2}\delta_{K}\sim 2\delta_{K}$ with $\delta_{K}$ being the
half-width of the spectral density discussed in Sec. V.

We also calculate the impurity entropy $S_{\rm imp}(T)$ for comparison with
experiments. $S_{\rm imp}(T)$ is defined as the difference of the entropy of
$\widetilde{H}_{N}$ and that of a reference free-fermion chain (without the
impurity site) defined by the same $\epsilon_{n},t_{n}$ parameters as in
$\widetilde{H}_{N}$. As shown in Fig. 4(c) and (d), $S_{\rm imp}(T)$ is zero
in the Fermi liquid phase at sufficiently low $T$ and starts to increase when
$T$ reaches the Kondo energy scale. For $\nu=1$, $S_{\rm imp}(T)$ climbs to
about $\ln 4\cdot k_{B}$ at about $k_{B}T\approx 1$meV and (approximately)
stays at this value until $k_{B}T$ reaches 10meV. The entropy $\ln 4\cdot
k_{B}\approx 1.39k_{B}$ is close to the measured value ($\sim 1.2k_{B}$) in
Refs. Rozen  _et al._ (2021); Saito  _et al._ (2021b) and can be understood as
contributed by the four degenerate states in the LM1 phase. Higher excited
states will be involved when $k_{B}T$ is larger than 10meV, and the entropy
continues to increase for larger $k_{B}T$. We also show the spectral density
at $\nu=1$ in Fig. 4(e), one can see that the Kondo resonance peak becomes
weak for $k_{B}T>1$meV, consistent with the entropy and spin susceptibility
results.

In Fig. 4(f), we plot the impurity entropy as a function of the filling at
various magnetic fields $B_{\rm loc}$ and $k_{B}T=0.82$meV ($T\approx 9.5$K).
For $B_{\rm loc}=0$, the entropy has two peaks around $\nu\approx 1$ and
$2.3$, respectively, and a dip around $\nu\approx 1.8$. This entropy curve
quantitatively matches the curve in Fig. 2(e) of Ref. Rozen  _et al._ (2021)
measured at $T\approx 10$K. Our calculation provides a simple understanding of
the two-peak feature. Looking at the phase diagram and the $\mu_{c},\mu_{f}$
trajectories in Fig. 2(b), (c), we find that the two peaks correspond to the
two domes of LM1 and LM2, respectively, where the Kondo temperature is
relatively lower, and the dip corresponds to the mixed valence state between
LM1 and LM2, where the Kondo temperature is higher due to valence fluctuation.
One can also observe the non-monotonous $\nu$-dependence of $T_{K}$ in Fig.
4(a), where $T_{K}$ as a function of $\nu$ is indicated by the black curve. A
finite $B_{\rm loc}$ will polarize the spin and hence suppress the entropy.
However, as discussed in Sec. IV.3, the four-fold degenerate LM1 states
consist of two spin-$\frac{1}{2}$ states due to the orbital degeneracy, thus,
a strong field will not completely quench the entropy. Instead, due to the
orbital degeneracy, the remaining entropy will be $\ln 2\cdot k_{B}\approx
0.69k_{B}$. At $\nu=1$, the entropy saturates to a constant about $\ln 2\cdot
k_{B}\approx 0.69k_{B}$ under a strong field, consistent with the measured
value ($\sim$0.6$k_{B}$) in Ref. Rozen  _et al._ (2021).

The spectral density and the Pomeranchuk effect were not connected in previous
works. Now our theory claims that both of them arise from the Kondo effect of
a local moment. The quantum-dot-behavior in spectral density corresponds to
the higher entropy LM1 state in the compressibility experiment, and the zero-
energy peak corresponds to the lower entropy singlet state. Suppose one can do
the STM and compressibility measurement at the same time, then she/he will
observe correlations between the two phenomena among different samples and
experimental conditions if our theory is correct, i.e., when the spectral
density exhibits the zero-energy peak (quantum-dot-behavior), the
compressibility measurement will show a lower (higher) entropy. This
correlation could be used to verify our theory in future experimental
researches.

In addition, the agreement between our calculation and the experimental data
in Ref. Rozen  _et al._ (2021) at $T\approx 10$K and a wide range of fillings
(including $\nu=\pm 2$) suggests that, even though our treatment ignores the
symmetry-breaking correlated gaps at $\nu=\pm 2$ at the zero temperature, our
results should apply to higher temperature states where the symmetries of the
CNP state are recovered due to thermal fluctuations.

##  VII Discussions

Figure 5: Spectral densities, electron bands, and energy surfaces at zero and
finite temperatures for the filling $\nu=1.5$. (a) The spectral densities at
$\nu=1.5$ at various temperatures. The Kondo temperature estimated from spin
susceptibility is about $k_{B}T_{K}=0.796$meV. (b) Heavy Fermi liquid band at
$\nu=1.5$, $T=0$. As explained in the text, the effective hybridization
between $c$\- and $f$-bands is suppressed by a factor $z^{\frac{1}{2}}$ with
$z$ being the quasi-particle weight of $f$-electrons, which are estimated as
0.188 here. The color of the bands represents the total quasi-particle weight,
which is always larger than $z$. (c) Band in the LM1 phase at $T>T_{K}$.
$f$-band is not shown here since it has zero quasi-particle weight $z$. (d),
(e) The energy surfaces at various energies for $T=0$ and $T>T_{K}$.
$\Gamma_{M},M_{M},K_{M}$ labels the high symmetry point in the moiré Brillouin
zone. The colors of the energy surfaces represent the total quasi-particle
weight. The numbers between energy surfaces indicate the number of bands below
the associated energy (counted from CNP).

Based on the NRG calculations, the poor man’s scaling, and various
experimental observations, we have shown that the gapless states at
$1\lesssim|\nu|<2$ are in the Kondo regime. Considering the translation
invariance of the actual system, these states must be the heavy Fermi liquid.
Here we estimate the heavy Fermi liquid bands from the information provided by
the NRG calculation. At the zero temperature, a spectral density in the Kondo
regime possesses a Lorentz peak around the zero energy, i.e.,
$A(\omega)\approx\frac{4z}{\pi}\frac{\delta_{K}}{\omega^{2}+\delta_{K}^{2}}$,
where $z$ is the quasi-particle weight, $\delta_{K}$ is the half-width defined
in Sec. V, and the factor 4 is from orbital and spin degeneracies. Then the
quasi-particle weight can be estimated as $z=\pi A(0)\delta_{K}/4$.
Substituting the quasi-particle part of $\mathcal{G}_{\alpha s}(i\omega)$,
i.e., $\frac{z}{i\omega}$, into Dyson’s equation of $c$-electrons

$\displaystyle\mathcal{G}^{(c)}(i\omega,{\mathbf{k}})=\mathcal{G}^{(c,0)}(i\omega,{\mathbf{k}})$
$\displaystyle\quad+\mathcal{G}^{(c,0)}(i\omega,{\mathbf{k}})H^{(cf)}({\mathbf{k}})\frac{z}{i\omega}H^{(cf)\dagger}({\mathbf{k}})\mathcal{G}^{(c)}(i\omega,{\mathbf{k}})\
,$ (33)

one can see that the $cf$ hybridization is effectively suppressed by a factor
of $z^{\frac{1}{2}}$. Here $\omega$ is the Matsubara frequency,
$\mathcal{G}^{(c,0)}(i\omega,{\mathbf{k}})=(i\omega-H^{(c)}({\mathbf{k}}))^{-1}$
is the free propagator of $c$-electrons, and $H^{(c)}({\mathbf{k}})$,
$H^{(cf)}({\mathbf{k}})$ are Hamiltonian matrices in Sec. III.2. One should be
aware that the Hubbard band information is completely neglected in this
method.

We now apply the above method to the $\nu=1.5$ state, where $k_{B}T_{K}$ is
estimated as 0.796meV (Table 2). Using the zero-temperature spectral density
in Fig. 5(a), we obtain $z\approx 0.188$ and the heavy Femri liquid bands
shown in Fig. 5(b). For $T>T_{K}$, the Lorentz peak at zero energy is smeared
by thermal fluctuations, corresponding to a LM phase. There is $z=0$ in the LM
phase due to the absence of Lorentz peak, hence the electron band structure is
solely contributed by the $c$-electrons, as shown in Fig. 5(c). (From Fig.
2(a) one can read that $\mu_{f}$ at $\nu=1.5$ is around 40meV and falls in the
range $0<\mu_{f}<U_{1}$, implying that the LM phase is the four-fold LM1
phase.) In Fig. 5(d) and (e) we plot energy surfaces around the Fermi level at
$T=0$ and $T>T_{K}$, respectively. The heavy Fermi liquid phase at $T=0$ has
two shells of electron-type Fermi surfaces with significantly different Fermi
wavevectors; the LM1 phase at $T>T_{K}$ has two degenerate shells of electron-
type Fermi surfaces. Here each shell of the Fermi surfaces is spin-degenerate.
At $T=0$, the two shells of Fermi surfaces enclose a total area corresponding
to the filling $\nu=1.5$; whereas at $T>T_{K}$ the Fermi surfaces enclose a
total area corresponding to the filling $0.5$ because the $f$-electron is not
conductive and does not contribute to electron bands in the LM1 phase. The
temperature dependence of the Fermi surface is a noticeable signature of the
heavy Fermi liquid phase. The most striking signature of energy surfaces
appear at higher energies. At $T=0$, because the Fermi level is close to the
top of the heavy band (Fig. 5(b)), hole-type surfaces will appear as the
energy slightly increases (Fig. 5(d)). In contrast, the energy surface at
$T>T_{K}$ is almost energy-independent (Fig. 5(e)). Similar energy dependent
behaviors of the energy surfaces occur at other fillings between 1 and 2. They
are a smoking gun signature of the heavy Fermi liquid phase at
$1\lesssim|\nu|<2$, that could be observed in spectral experiments such as the
quasi-particle interference.

Figure 6: Spin-selective orbital Kondo effect. (a) and (b) are the total and
spin-selective spectral densities at $T=0$ and various in-plane magnetic
fields, respectively.

Another observable signature of the heavy Fermi liquid phase is the field-
dependence of the spectrum. In Sec. VI we have pointed out that, due to the
orbital degeneracy, the four-fold LM1 state will reduce to a two-fold LM state
under a strong in-plane magnetic field. We find that, when cooling down, the
remaining orbital degeneracy under a strong field leads to a spin-selective
orbital Kondo effect. In Fig. 6(a), we demonstrate the spectral densities
under in-plane magnetic fields at $\nu=1.5$. (Here we use local field $B_{\rm
loc}$ that only acts on the spin d.o.f. of the impurity to approximate the
physical field $B$ that also acts on $c$-electrons and has orbital effects.
For $B<20$T, the Zeeman splitting $\mu_{B}B<1.2$meV is much smaller than
$\mu_{c}\approx 37$meV at $\nu=1.5$, hence its action on $c$-electrons is
negligible. As for the in-plane orbital effect Cao  _et al._ (2021a), it can
be avoided in the magic-angle mirror symmetric twisted trilayer graphene Cao
_et al._ (2021b), the low energy Hamiltonian of which is similar to MATBG.)
When the Zeeman energy is much smaller than $k_{B}T_{K}$, $B_{\rm loc}$ does
not qualitatively affect the spectral density. When the Zeeman energy is
comparable to the Kondo scale, i.e., $B_{\rm loc}\gtrsim 10$T such that
$\mu_{B}B_{\rm locl}\gtrsim 0.58$meV, the Kondo peak splits into a thinner
peak at zero energy and a broadened peak at positive energy away from zero. As
shown in the spin-selective spectral densities in Fig. 6(b), the thinner
resonance peak is contributed by the spin-up (parallel to $B$) electrons, and
the broadened peak is contributed by the spin-down electrons. This behavior
can be understood as a spin-selective orbital Kondo effect. $f$-electrons will
be transferred from the two spin-down $f$-orbitals to the two spin-up
$f$-orbitals by a strong field. Then, the spin-down sector will become an
effective frozen impurity state. The two spin-up $f$-orbitals only contain one
electron - that forms LM1 in the absence of field - and subject to an orbital
Kondo effect, leading to the thinner resonance peak. This spin-selective
orbital Kondo effect does not depend on the details of our calculations.
Instead, it will take place as long as (i) the orbital degeneracy is not
broken, (ii) the ground state of the two orbitals has $n_{f}=1$, because with
these two conditions the spin-up orbitals form a standard $\mathcal{N}=2$
Anderson model. The spin-selective orbital Kondo effect could be observed in
STM experiments.

In this work we have ignored RKKY interactions between LMs at different
$f$-sites that will be induced by $c$-bands. These RKKY interactions can lead
to further symmetry breaking and should be crucial to stabilize the observed
correlated insulator states at $|\nu|=2$. Nevertheless, as we have explained
in Sec. III.3 and demonstrated in Secs. VI and V, where the computed spectral
density and entropy match experimental data at finite temperatures, our
treatment should have yielded correct $|\nu|\geq 2$ states at finite
temperatures where symmetries are recovered due to thermal fluctuations. In
order to have the correlated gaps at the zero temperature, a full self-
consistent treatment including both symmetry-breaking and Kondo screening
effects should be used. It may result in more complicated $\nu$-dependencies
of $\mu_{c},\mu_{f}$ than those shown Fig. 2(a), (b), (c). For example, at
$\nu=2$, $\mu_{c}$ given by a self-consistent symmetry-breaking Hartree-Fock
calculation Song and Bernevig (2022) is about 26meV, which is significantly
lower than the one ($\sim$40meV) in Fig. 2(a). Thus, a possible mechanism for
the correlated insulators at $|\nu|=2$ to win the Kondo screening is that
$\mu_{c}$ drops to a small value such that the Kondo energy scale becomes
irrelevant. (See Fig. 2(b), (c).) The competition between symmetry-breaking
and Kondo screening is also a potential mechanism for the observed strange
metal behaviors Polshyn  _et al._ (2019); Cao  _et al._ (2020); Jaoui  _et
al._ (2022) and could play an important role in the unconventional
superconductivity Cao  _et al._ (2018a); Lu  _et al._ (2019); Yankowitz  _et
al._ (2019); Saito  _et al._ (2020, 2020); Stepanov  _et al._ (2020); Liu  _et
al._ (2021a); Arora  _et al._ (2020); Cao  _et al._ (2021a); Oh  _et al._
(2021). We leave this for future studies.

Note added. During the preparation of the current work, a related work Chou
and Sarma (2022) appeared. This work studied the symmetric Kondo state using a
slave-fermion mean field in a Kondo lattice model derived from the THF model.
Our theory is based on the symmetry-breaking correlated state at CNP. We are
also aware of related works on the Kondo problem in MATBG by A. M. Tsvelik’s
and B. A. Bernevig’s group Hu  _et al._ (2023a, b) and P. Coleman’s group Lau
and Coleman (2023), and a generalization of the THF model to the magic-angle
twisted trilayer graphene Yu  _et al._ (2023). Ref. Hu  _et al._ (2023a) also
obtains a Kondo temperature about $1\sim 2$K around $|\nu|\approx 1$.

###### Acknowledgements.

We are grateful to B. Andrei Bernevig, Ning-Hua Tong, Xi Dai, Jia-Bin Yu,
Xiao-Bo Lu, Yong-Long Xie, Yi-Lin Wang, and Chang-Ming Yue for helpful
discussions. Z.-D. S. and G.-D. Z. were supported by National Natural Science
Foundation of China (General Program No. 12274005), National Key Research and
Development Program of China (No. 2021YFA1401900).

##  Appendix A More details about the effective Hamiltonian

###  A.1 Nonzero $M$ term

A generic trial ground state at CNP is given by (Eq. 6)

$|\Psi_{0}\rangle=U\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}|{\rm
FS}\rangle\ ,$ (34)

where $U=\exp(-i\theta_{\mu\nu}\hat{\Sigma}_{\mu\nu})$ is a U(4) rotation
operator and an implicit summation over repeated $\mu,\nu$ indices is assumed.
We can always define the rotated fermion operators
$\widetilde{c}_{{\mathbf{k}}a\eta s}=Uc_{{\mathbf{k}}a\eta s}U^{\dagger}$,
$\widetilde{f}_{\mathbf{R}a\eta s}=Uf_{\mathbf{R}a\eta s}U^{\dagger}$ such
that $\widetilde{f}_{\mathbf{R}\alpha+s}$’s are occupied in $|\Psi_{0}\rangle$
and $\widetilde{f}_{\mathbf{R}\alpha-s}$’s are empty in $|\Psi_{0}\rangle$.
According to the discussions in the supplementary material section S4B of Ref.
Song and Bernevig (2022), in the flat-band limit ($M=0$), the lowest particle
(hole) excitations only involve $\widetilde{c}_{{\mathbf{k}}a-s}$ and
$\widetilde{f}_{\mathbf{R}a-s}$ ($\widetilde{c}_{{\mathbf{k}}a+s}$ and
$\widetilde{f}_{\mathbf{R}a+s}$). Thus, the effective periodic Anderson model
for $\nu>0$ derived in Sec. III.2 is written in terms of
$c_{{\mathbf{k}}as}=\widetilde{c}_{{\mathbf{k}}a-s}$ and $f_{\mathbf{R}\alpha
s}=\widetilde{f}_{\mathbf{R}\alpha-s}$. Here we give the explicit forms of the
rotated operators

$\widetilde{f}_{\mathbf{R}\alpha\eta
s}=\sum_{\alpha^{\prime}\eta^{\prime}s^{\prime}}\left[e^{i\theta_{\mu\nu}\Sigma^{f}_{\mu\nu}}\right]_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}\
,$ (35)

and

$\widetilde{c}_{{\mathbf{k}}a\eta
s}=\sum_{\begin{subarray}{c}a^{\prime}=1,2\\\
\eta^{\prime}s^{\prime}\end{subarray}}\left[e^{i\theta_{\mu\nu}\Sigma^{c12}_{\mu\nu}}\right]_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{{\mathbf{k}}a^{\prime}\eta^{\prime}s^{\prime}}\quad(a=1,2),$
(36) $\widetilde{c}_{{\mathbf{k}}a\eta
s}=\sum_{\begin{subarray}{c}a^{\prime}=3,4\\\
\eta^{\prime}s^{\prime}\end{subarray}}\left[e^{i\theta_{\mu\nu}\Sigma^{c34}_{\mu\nu}}\right]_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{{\mathbf{k}}a^{\prime}\eta^{\prime}s^{\prime}}\quad(a=3,4),$
(37)

where the eight-by-eight matrices $\Sigma^{f}_{\mu\nu}$,
$\Sigma^{c12}_{\mu\nu}$, $\Sigma^{c34}_{\mu\nu}$ are defined in Eqs. 3, 4 and
5.

The $M$-term in the original basis of the THF model (Eq. 1) is

$M\sum_{aa^{\prime}=3,4}\sum_{{\mathbf{k}}}\sum_{\eta
s}[\sigma_{x}]_{aa^{\prime}}c_{{\mathbf{k}}a\eta
s}^{\dagger}c_{{\mathbf{k}}a^{\prime}\eta s}\ .$ (38)

It favors the Kramers inter-valley coherent state discussed at the end of Sec.
II, where $\theta_{x0}$ and $\theta_{y0}$ are nonzero and satisfy
$\theta_{x0}^{2}+\theta_{y0}^{2}=(\pi/4)^{2}$. Without loss of generality, we
assume $U=\exp(-i\frac{\pi}{4}\hat{\Sigma}_{x0})$ for the Kramers inter-valley
coherent state. Writing this $M$-term in terms of the rotated operators, we
obtain

$\displaystyle M$
$\displaystyle\sum_{{\mathbf{k}}}\sum_{a,a^{\prime}=3,4}\sum_{\eta\eta^{\prime}ss^{\prime}}\widetilde{c}_{{\mathbf{k}}a\eta
s}^{\dagger}O_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}\widetilde{c}_{{\mathbf{k}}a^{\prime}\eta^{\prime}s^{\prime}}\
,$ (39)

where
$O=e^{i\frac{\pi}{4}\Sigma^{c34}_{x0}}\sigma_{x}\tau_{0}\varsigma_{0}e^{-i\frac{\pi}{4}\Sigma^{c34}_{x0}}=-\sigma_{z}\tau_{x}\varsigma_{0}$.
The $\tau_{x}$ matrix in $O$ represents couplings between the empty and
occupied single-particle states. If we simply project this $M$-term onto the
empty states, it vanishes, i.e., $[O_{a-s,a^{\prime}-s^{\prime}}]=0$. A better
approximation is applying a Schrieffer-Wolff transformation to decouple the
$\eta=\pm$ states, leading to a second-order correction to the effective
Hamiltonian. As $\langle\Psi_{0}|\widetilde{f}_{\alpha\eta
s}^{\dagger}\widetilde{f}_{\alpha\eta s}|\Psi_{0}\rangle=(1+\eta)/2$, the $J$
term in Eq. 2 yields the following mean field term (see also the supplementary
material section S4B of Ref. Song and Bernevig (2022))

$-\frac{J}{2}\sum_{a=3,4}\sum_{\eta s}\eta\cdot\widetilde{c}_{a\eta
s}^{\dagger}\widetilde{c}_{a\eta s}$ (40)

Then, regarding the $\frac{J}{2}$ term as the zeroth order Hamiltonian and $M$
as a perturbation, a Schrieffer-Wolff transformation leads to the correction

$-\frac{M^{2}}{J}\sum_{|{\mathbf{k}}|<\Lambda_{c}}\sum_{a=3,4}\sum_{\eta
s}\eta\cdot\widetilde{c}_{{\mathbf{k}}a\eta
s}^{\dagger}\widetilde{c}_{{\mathbf{k}}a\eta s}+\mathcal{O}(M^{4})\ .$ (41)

The resulting energy levels $\pm(J/2+M/J^{2})$ at ${\mathbf{k}}=0$ is fully
consistent with a Taylor expansion of the one-shot energy levels
$\pm\sqrt{J^{2}/4+M^{2}}$ derived in Ref. Song and Bernevig (2022). Projecting
the correction to the active d.o.f., i.e.,
$c_{{\mathbf{k}}as}=\widetilde{c}_{{\mathbf{k}}a-s}$, we obtain the correction
to the effective Hamiltonian

$\frac{M^{2}}{J}\sum_{|{\mathbf{k}}|<\Lambda_{c}}\sum_{a=3,4}\sum_{s}{c}_{{\mathbf{k}}as}^{\dagger}{c}_{{\mathbf{k}}as}+\mathcal{O}(M^{4})\
.$ (42)

###  A.2 Hund’s coupling

The four-by-four Hamiltonian matrix $H^{(c)}({\mathbf{k}})+\Delta H^{(c)}$ in
Eq. 13, i.e.,
$-v_{\star}(\sigma_{x}\otimes\sigma_{0}k_{x}+\sigma_{y}\otimes\sigma_{z}k_{y})+\mathbb{0}_{2\times
2}\oplus G\sigma_{0}$, can be diagonalized analytically. As discussed at the
end of the last subsection, to $\mathcal{O}(M^{2})$, the $M$ term simply
shifts the energy of $a=3,4$ electrons by $M^{2}/J$. Thus, all the analysis
below applies to the $M\neq 0$ after $G$ is replaced by $G+M^{2}/J$. We find
the energy eigenvalues and wave-functions of the $H^{(c)}({\mathbf{k}})+\Delta
H^{(c)}$ as

$\begin{aligned}
\epsilon_{1}({\mathbf{k}})=&\epsilon_{+}({\mathbf{k}})=\frac{G}{2}+\sqrt{\frac{G^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{1}({\mathbf{k}})=&\begin{pmatrix}\sin\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}&0&-\cos\frac{\theta_{\mathbf{k}}}{2}&0\end{pmatrix}^{T}\end{aligned}\
,$ (43) $\begin{aligned}
\epsilon_{2}({\mathbf{k}})=&\epsilon_{+}({\mathbf{k}})=\frac{G}{2}+\sqrt{\frac{G^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{2}({\mathbf{k}})=&\begin{pmatrix}0&\sin\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}&0&-\cos\frac{\theta_{\mathbf{k}}}{2}\end{pmatrix}^{T}\end{aligned}\
,$ (44) $\begin{aligned}
\epsilon_{3}({\mathbf{k}})=&\epsilon_{-}({\mathbf{k}})=\frac{G}{2}-\sqrt{\frac{G^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{3}({\mathbf{k}})=&\begin{pmatrix}\cos\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}&0&\sin\frac{\theta_{\mathbf{k}}}{2}&0\end{pmatrix}^{T}\end{aligned}\
,$ (45) $\begin{aligned}
\epsilon_{4}({\mathbf{k}})=&\epsilon_{-}({\mathbf{k}})=\frac{G}{2}-\sqrt{\frac{G^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{4}({\mathbf{k}})=&\begin{pmatrix}0&\cos\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}&0&\sin\frac{\theta_{\mathbf{k}}}{2}\end{pmatrix}^{T}\end{aligned}\
.$ (46)

where

$\theta_{\mathbf{k}}=\arccos\frac{G/2}{\sqrt{G^{2}/4+v_{\star}^{2}{\mathbf{k}}^{2}}}$
(47)

and $\phi_{\mathbf{k}}=\arg(k_{x}+ik_{y})$.

We now derive the effective Hund’s coupling $\hat{H}_{H}$. Applying a second-
order perturbation in terms of $\hat{H}_{J}$, we obtained the correction to
the Hamiltonian

$\displaystyle\Delta\hat{H}=-\frac{J^{2}}{N_{M}^{2}}\sum_{I}\sum_{\begin{subarray}{c}\alpha_{1}\alpha_{2}\\\
s_{1}s_{1}^{\prime}s_{2}s_{2}^{\prime}\end{subarray}}\sum_{\begin{subarray}{c}{\mathbf{k}}_{1},{\mathbf{k}}_{1}^{\prime}\\\
{\mathbf{k}}_{2},{\mathbf{k}}_{2}^{\prime}\end{subarray}}(f_{\alpha_{1}s_{1}}^{\dagger}f_{\alpha_{1}s_{1}^{\prime}}-\frac{\nu_{f}}{4}\delta_{s_{1}s_{1}^{\prime}})$
$\displaystyle\times(f_{\alpha_{2}s_{2}^{\prime}}^{\dagger}f_{\alpha_{2}s_{2}}-\frac{\nu_{f}}{4}\delta_{s_{2}s_{2}^{\prime}})\cdot
e^{-\frac{\lambda^{2}}{2}({\mathbf{k}}_{1}^{2}+{\mathbf{k}}_{1}^{\prime
2}+{\mathbf{k}}_{2}^{2}+{\mathbf{k}}_{2}^{\prime 2})}$
$\displaystyle\times\frac{\langle\Psi_{0}|c_{{\mathbf{k}}_{1}^{\prime}\alpha_{1}+2s_{1}^{\prime}}^{\dagger}c_{{\mathbf{k}}_{1}\alpha_{1}+2s_{1}}|\Psi_{I}\rangle\langle\Psi_{I}|c_{{\mathbf{k}}_{2}\alpha_{2}+2s_{2}}^{\dagger}c_{{\mathbf{k}}_{2}^{\prime}\alpha_{2}+2s_{2}^{\prime}}|\Psi_{0}\rangle}{E_{I}-E_{0}}\
,$ (48)

where $|\Psi_{I}\rangle$ are excited states with a single particle-hole pair
and $E_{I}$ are the energies of the excited states. ${\mathbf{k}}_{1,2}$,
${\mathbf{k}}_{1,2}^{\prime}$ are limited up to the cutoff $\Lambda_{c}$. Due
to the momentum and spin conservation, for the matrix element to be nonzero,
there must be ${\mathbf{k}}_{1}={\mathbf{k}}_{2}$, $s_{1}=s_{2}$,
${\mathbf{k}}_{1}^{\prime}={\mathbf{k}}_{2}^{\prime}$,
$s_{1}^{\prime}=s_{2}^{\prime}$. For simplicity, we rewrite
${\mathbf{k}}_{1}$, ${\mathbf{k}}_{1}^{\prime}$, $s_{1}$, and $s_{1}^{\prime}$
as ${\mathbf{k}}$, ${\mathbf{k}}^{\prime}$, $s$, and $s^{\prime}$,
respectively. $({\mathbf{k}},s)$ and $({\mathbf{k}}^{\prime},s^{\prime})$
label the particle and the hole excitations, respectively. Then the matrix
element in the third line of Sec. A.2 can be written as

$\displaystyle
n_{F}(\epsilon_{-}({\mathbf{k}}^{\prime})-\mu_{c})(1-n_{F}(\epsilon_{+}({\mathbf{k}})-\mu_{c}))$
$\displaystyle\times\langle\Psi_{0}|\wick{\c{1}c_{{\mathbf{k}}^{\prime}\alpha_{1}+2s^{\prime}}^{\dagger}\c{2}c_{{\mathbf{k}}\alpha_{1}+2s}\c{2}c^{\dagger}_{{\mathbf{k}}\alpha_{2}+2s}\c{1}c_{{\mathbf{k}}^{\prime}\alpha_{2}+2s^{\prime}}}|\Psi_{0}$

According to the wave functions given in Eqs. (43)-(46), there are
$\langle\Psi_{0}|c_{{\mathbf{k}}^{\prime}\alpha_{1}+2s^{\prime}}^{\dagger}c_{{\mathbf{k}}^{\prime}\alpha_{2}+2s^{\prime}}|\Psi_{0}\rangle=\delta_{\alpha_{1}\alpha_{2}}\sin^{2}\frac{\theta_{{\mathbf{k}}^{\prime}}}{2}$,
$\langle\Psi_{0}|c_{{\mathbf{k}}\alpha_{1}+2s}c^{\dagger}_{{\mathbf{k}}\alpha_{2}+2s}|\Psi_{0}\rangle=\delta_{\alpha_{1}\alpha_{2}}\cos^{2}\frac{\theta_{{\mathbf{k}}}}{2}$.
The excitation energy $E_{I}-E_{0}$ is given by
$\epsilon_{+}({\mathbf{k}})-\epsilon_{-}({\mathbf{k}}^{\prime})$. Thus,
$\Delta\hat{H}$ is simplified to

$\displaystyle\Delta\hat{H}=-\frac{J^{2}}{N_{M}^{2}}\sum_{\begin{subarray}{c}\alpha
ss^{\prime}\\\ {\mathbf{k}}{\mathbf{k}}^{\prime}\end{subarray}}(f_{\alpha
s}^{\dagger}f_{\alpha
s^{\prime}}-\frac{\nu_{f}}{4}\delta_{ss^{\prime}})(f_{\alpha
s^{\prime}}^{\dagger}f_{\alpha s}-\frac{\nu_{f}}{4}\delta_{ss^{\prime}})$
$\displaystyle\times\frac{n_{F}(\epsilon_{-}({\mathbf{k}}^{\prime})-\mu_{c})(1-n_{F}(\epsilon_{+}({\mathbf{k}})-\mu_{c}))\sin^{2}\frac{\theta_{{\mathbf{k}}^{\prime}}}{2}\cos^{2}\frac{\theta_{{\mathbf{k}}}}{2}}{\epsilon_{+}({\mathbf{k}})-\epsilon_{-}({\mathbf{k}}^{\prime})}$
$\displaystyle\times e^{-\lambda^{2}({\mathbf{k}}^{2}+{\mathbf{k}}^{\prime
2})}$ (49)

The $s=s^{\prime}$ contribution is an effective chemical potential shift,
estimated as 0.17meV at CNP, of the $f$-electrons. As it is much smaller than
$U_{1}$, we will omit the $s=s^{\prime}$ contribution. The $s\neq s^{\prime}$
contribution can be written as

$\hat{H}_{H}=J_{H}\sum_{\alpha}n_{\alpha\uparrow}^{f}n_{\alpha\downarrow}^{f}$
(50)

with $J_{H}$ given by

$\displaystyle J_{H}=$ $\displaystyle
2J^{2}\left(\frac{\Omega_{0}}{2\pi}\right)^{2}\int_{0}^{\Lambda_{c}}dk^{\prime}\cdot
k^{\prime}\int_{k_{0}}^{\Lambda_{c}}dk\cdot k\cdot
e^{-\lambda^{2}(k^{2}+k^{\prime 2})}$
$\displaystyle\times\frac{\sin^{2}\frac{\theta_{k^{\prime}}}{2}\cos^{2}\frac{\theta_{k}}{2}}{\epsilon_{+}(k)-\epsilon_{-}(k^{\prime})}\
,$ (51)

where $k_{0}$ is determined by $\epsilon_{+}(k_{0})=\mu_{c}$. Here we have
made use of the fact that $\epsilon_{\pm}({\mathbf{k}})$ and
$\theta_{{\mathbf{k}}}$ only depends on $|{\mathbf{k}}|$ but not
$\phi_{\mathbf{k}}$. At CNP, $\mu_{c}=0$ and $G=J/2=8.19$meV, taking the limit
$\Lambda_{c}\to\infty$, we obtain $J_{H}\approx 0.34$meV. Using the self-
consistent values of $\mu_{c}$ and $G$ shown in Fig. 2(a), we find $J_{H}$ at
$\nu=1,2,3,4$ are given by 0.29meV, 0.26meV, 0.21meV, 0.19meV, respectively.
As $J_{H}$ is small and does not change significantly with $\nu$, in this
work, we simply set $J_{H}=0.34$meV for simplicity.

###  A.3 Hybridization function

By definition, the hybridization function $\Delta(\omega)$ is given by

$\displaystyle\Delta(\omega)=\frac{\pi}{N}\sum_{{\mathbf{k}}}\sum_{n}|V_{n\alpha}({\mathbf{k}})|^{2}\delta(\omega-\epsilon_{n}({\mathbf{k}}))$
(52)

where
$V_{n\alpha}({\mathbf{k}})=\sum_{a}u^{*}_{an}({\mathbf{k}})H^{(cf)}_{a\alpha}({\mathbf{k}})e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}$
is the hybridization between $f_{\alpha s}$ and the $n$-th energy band of
$c$-electrons. $\Delta(\omega)$ does not depend on $\alpha$ because of the
$C_{2z}T$ or $C_{2x}$ symmetry that flips the $\alpha$ index. Substituting
$\epsilon_{n}({\mathbf{k}})$ and $u_{an}({\mathbf{k}})$ in Eqs. (43)-(46) into
the above equation, we obtain $V_{n\alpha}({\mathbf{k}})$ for $\alpha=1$ as

$\displaystyle V_{11}({\mathbf{k}})=$
$\displaystyle\gamma\sin\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}$
$\displaystyle V_{21}({\mathbf{k}})=$ $\displaystyle
v_{\star}^{\prime}(-k_{x}+ik_{y})\sin\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}$
$\displaystyle V_{31}({\mathbf{k}})=$
$\displaystyle\gamma\cos\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}$
$\displaystyle V_{41}({\mathbf{k}})=$ $\displaystyle
v_{\star}^{\prime}(-k_{x}+ik_{y})\cos\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}\
.$ (53)

Using the energy eigenvalues in Eqs. (43)-(46) and the
$V_{n\alpha}({\mathbf{k}})$ matrix elements given above, it is direct to
obtain

$\displaystyle\Delta(\omega)=\frac{\Omega_{0}}{2v_{\star}^{2}}\left|\omega+\mu_{c}-\frac{G}{2}\right|\left(\gamma^{2}+v_{\star}^{\prime
2}k_{F}^{2}\right)e^{-k_{F}^{2}\lambda^{2}}$
$\displaystyle\times\left[\theta(\omega+\mu_{c}-G)\sin^{2}\frac{\theta_{k_{F}}}{2}+\theta(-\omega-\mu_{c})\cos^{2}\frac{\theta_{k_{F}}}{2}\right]$
(54)

where $k_{F}$ is given by

$\displaystyle
k_{F}=\frac{1}{v_{\star}}\sqrt{\left(\omega+\mu_{c}-{G}/{2}\right)^{2}-\left(G/{2}\right)^{2}}\
.$ (55)

Making use of Eq. 47, Sec. A.3 can be further simplified to

$\displaystyle\Delta(\omega)=$
$\displaystyle\frac{\Omega_{0}}{4v_{\star}^{2}}\left|\omega+\mu_{c}-G\right|\left(\gamma^{2}+v_{\star}^{\prime
2}k_{F}^{2}\right)e^{-k_{F}^{2}\lambda^{2}}$
$\displaystyle\times\left[\theta(\omega+\mu_{c}-G)+\theta(-\omega-\mu_{c})\right]\
.$ (56)

We now verify the asymptotic behaviors of $\Delta(\omega)$. When
$\omega+\mu_{c}\to G^{+}$, $k_{F}\to 0$ and only the first term in the second
line of Sec. A.3 contributes to $\Delta(\omega)$. Then we obtain the
asymptotic behavior of $\Delta(\omega)$ as $\omega+\mu_{c}\to G^{+}$

$\Delta(\omega)=\frac{\Omega_{0}}{4v_{\star}^{2}}\gamma^{2}\cdot(\omega+\mu_{c}-G)+\mathcal{O}((\omega+\mu_{c}-G)^{2})\
.$ (57)

When $\omega+\mu_{c}\to-0^{+}$, $k_{F}\to 0$ and only the second term in the
second line of Sec. A.3 contributes to $\Delta(\omega)$. Then we obtain the
asymptotic behavior of $\Delta(\omega)$ as $\omega+\mu_{c}\to-0^{+}$

$\Delta(\omega)=\frac{\Omega_{0}}{4v_{\star}^{2}}G\gamma^{2}+\mathcal{O}((\omega+\mu_{c}))\
.$ (58)

##  Appendix B Poor man’s scaling of Anderson models with energy-dependent
couplings

###  B.1 Generic theory for U($\mathcal{N}$) models

We consider the Anderson impurity model with $\mathcal{N}$ symmetric flavors

$\displaystyle\hat{H}=$
$\displaystyle-\mu_{f}\hat{N}_{f}+\frac{U}{2}\hat{N}_{f}(\hat{N}_{f}-1)+\sum_{\mu=1}^{\mathcal{N}}\int_{-D}^{D}d\epsilon\cdot\epsilon\cdot
d_{\mu}^{\dagger}(\epsilon)d_{\mu}(\epsilon)$
$\displaystyle+\sum_{\mu=1}^{\mathcal{N}}\int_{-D}^{D}d\epsilon\sqrt{\frac{\Delta(\epsilon)}{\pi}}(f_{\mu}^{\dagger}d_{\mu}(\epsilon)+h.c.)\
,$ (59)

where $\mu$ is the flavor index and
$N_{f}=\sum_{\mu=1}^{\mathcal{N}}f_{\mu}^{\dagger}f_{\mu}$. We assume the
ground state of the isolated impurity has $n_{f}$ $f$-electrons, which can
take the values $1,2\cdots(N-1)$. (We do not consider the empty case
($n_{f}=0$), the full case ($n_{f}=N$), and the mixed valence case where
ground states with different $n_{f}$ are exactly degenerate.) We further
assume the charge gaps to $n_{f}-1$ and $n_{f}+1$ electrons are $\Delta E_{-}$
and $\Delta E_{+}=U-\Delta E_{-}$, respectively. We then apply a Schrieffer-
Wolff transformation to obtain an effective Coqblin–Schrieffer model for the
Hilbert space restricted to $\hat{N}_{f}=n_{f}$

$\displaystyle\hat{H}=$
$\displaystyle\sum_{\mu=1}^{\mathcal{N}}\int_{-D}^{D}d\epsilon\cdot\epsilon\cdot
d_{\mu}^{\dagger}(\epsilon)d_{\mu}(\epsilon)+\frac{4g}{\pi
U}\sum_{\mu,\mu^{\prime}=1}^{\mathcal{N}}\int_{-D}^{D}d\epsilon
d\epsilon^{\prime}\bigg{[}$
$\displaystyle\qquad\sqrt{{\Delta(\epsilon)\Delta(\epsilon^{\prime})}}(f_{\mu}^{\dagger}f_{\mu^{\prime}}-x\delta_{\mu\mu^{\prime}})d_{\mu^{\prime}}^{\dagger}(\epsilon^{\prime})d_{\mu}(\epsilon)\bigg{]}\
.$ (60)

Terms that only involve $\hat{N}_{f}$ are omitted because they only contribute
to an energy shift. The bandwidth $D$ should be smaller than $\min(\Delta
E_{+},\Delta E_{-})$, otherwise, the Schrieffer-Wolff transformation is
invalid. The parameters $g$ and $x$ are given by

$g=\frac{U}{4}\left(\frac{1}{\Delta E_{+}}+\frac{1}{\Delta
E_{-}}\right),\qquad x=\frac{\Delta E_{-}}{U}\ ,$ (61)

respectively. If $\mu_{f}=(n_{f}-\frac{1}{2})U$, there is $\Delta E_{+}=\Delta
E_{-}=\frac{1}{2}U$ and $g=1$, $x=\frac{1}{2}$.

We now truncate the bandwidth at $D-dD=D(1-dt)$ ($dt\ll 1$) and consider
second order (in $g$) corrections form the virtual particle
($D-dD<\epsilon<D$) and hole ($-D<\epsilon<-D+dD$) excitations. The particle
excitation contributes to the correction

$\displaystyle-\frac{(4g)^{2}}{(\pi
U)^{2}}\frac{1}{D}\sum_{\mu_{1}\mu_{2}\mu_{1}^{\prime}\mu_{2}^{\prime}}\int_{-D+dD}^{D-dD}d\epsilon_{1}d\epsilon_{2}d\int_{D-dD}^{D}d\epsilon_{1}^{\prime}d\epsilon_{2}^{\prime}$
$\displaystyle\times\sqrt{\Delta(\epsilon_{1})\Delta(\epsilon_{2})\Delta(\epsilon_{1}^{\prime})\Delta(\epsilon_{2}^{\prime})}d_{\mu_{1}}^{\dagger}(\epsilon_{1})\langle
d_{\mu_{1}^{\prime}}(\epsilon_{1}^{\prime})d_{\mu_{2}^{\prime}}^{\dagger}(\epsilon_{2}^{\prime})\rangle
d_{\mu_{2}}(\epsilon_{2})$
$\displaystyle\times(f_{\mu_{1}^{\prime}}^{\dagger}f_{\mu_{1}}-x\delta_{\mu_{1}\mu_{1}^{\prime}})\mathbb{P}(f_{\mu_{2}}^{\dagger}f_{\mu_{2}^{\prime}}-x\delta_{\mu_{2}\mu_{2}^{\prime}})\
.$ (62)

The denominator $D$ in the factor is the excitation energy of a virtual
particle. $\mathbb{P}$ is a projector to the restricted Hilbert space, where
$\hat{N}_{f}=n_{f}$. The expectation $\langle
d_{\mu_{1}^{\prime}}(\epsilon_{1}^{\prime})d_{\mu_{2}^{\prime}}^{\dagger}(\epsilon_{2}^{\prime})\rangle$
evaluated on the ground state is
$\delta(\epsilon_{1}^{\prime}-\epsilon_{2}^{\prime})\delta_{\mu_{1}^{\prime}\mu_{2}^{\prime}}$.
Then we have

$\displaystyle-\frac{(4g)^{2}}{(\pi
U)^{2}}\frac{dD}{D}\Delta(D)\sum_{\mu_{1}\mu_{2}\mu^{\prime}}\int_{-D+dD}^{D-dD}d\epsilon_{1}d\epsilon_{2}\sqrt{\Delta(\epsilon_{1})\Delta(\epsilon_{2})}$
$\displaystyle\times
d_{\mu_{1}}^{\dagger}(\epsilon_{1})d_{\mu_{2}}(\epsilon_{2})(f_{\mu^{\prime}}^{\dagger}f_{\mu_{1}}-x\delta_{\mu_{1}\mu^{\prime}})(f_{\mu_{2}}^{\dagger}f_{\mu^{\prime}}-x\delta_{\mu_{2}\mu^{\prime}})\
,$ (63)

where $\mathbb{P}$ is omitted as it commutes with
$f_{\mu_{2}}^{\dagger}f_{\mu^{\prime}}$ and
$f_{\mu^{\prime}}^{\dagger}f_{\mu_{1}}$. After a few steps of algebra, the
four-fermion operator
$\sum_{\mu^{\prime}}f_{\mu^{\prime}}^{\dagger}f_{\mu_{1}}f_{\mu_{2}}^{\dagger}f_{\mu^{\prime}}$
can be simplified to

$f_{\mu_{2}}^{\dagger}f_{\mu_{1}}+\sum_{\mu^{\prime}}f_{\mu^{\prime}}^{\dagger}f_{\mu^{\prime}}f_{\mu_{1}}f_{\mu_{2}}^{\dagger}=f_{\mu_{2}}^{\dagger}f_{\mu_{1}}(1-n_{f})+n_{f}\delta_{\mu_{1}\mu_{2}}\
,$ (64)

where we have made use of the fact that the Hilbert space is restricted to
$\hat{N}_{f}=n_{f}$. Substituting this into Sec. B.1, we obtain the
corrections to parameters $g$ and $xg$ as

$\frac{dg}{dt}\bigg{|}_{p}=\frac{4\Delta(D(t))}{\pi
U}\left((n_{f}-1)+2x\right)g^{2}\ ,$ (65)
$\frac{d(xg)}{dt}\bigg{|}_{p}=\frac{4\Delta(D(t))}{\pi
U}\left(x^{2}+n_{f}\right)g^{2}\ .$ (66)

Here $t$ is the RG parameter and $D(t)=De^{-t}$ is the reduced bandwidth after
successive $t/dt$ RG steps.

We then calculate the contribution from virtual hole excitation. Following the
same process as in the last paragraph, we obtain

$\displaystyle-\frac{(4g)^{2}}{(\pi
U)^{2}}\frac{dD}{D}\Delta(-D)\sum_{\mu_{1}\mu_{2}\mu^{\prime}}\int_{-D+dD}^{D-dD}d\epsilon_{1}d\epsilon_{2}\sqrt{\Delta(\epsilon_{1})\Delta(\epsilon_{2})}$
$\displaystyle\times
d_{\mu_{1}}(\epsilon_{1})d^{\dagger}_{\mu_{2}}(\epsilon_{2})(f_{\mu_{1}}^{\dagger}f_{\mu^{\prime}}-x\delta_{\mu_{1}\mu^{\prime}})\mathbb{P}(f_{\mu^{\prime}}^{\dagger}f_{\mu_{2}}-x\delta_{\mu_{2}\mu^{\prime}})$
$\displaystyle=$ $\displaystyle\frac{(4g)^{2}}{(\pi
U)^{2}}dt\Delta(-D)\sum_{\mu_{1}\mu_{2}\mu^{\prime}}\int_{-D+dD}^{D-dD}d\epsilon_{1}d\epsilon_{2}\sqrt{\Delta(\epsilon_{1})\Delta(\epsilon_{2})}$
$\displaystyle\times
d^{\dagger}_{\mu_{2}}(\epsilon_{2})d_{\mu_{1}}(\epsilon_{1})(f_{\mu_{1}}^{\dagger}f_{\mu^{\prime}}-x\delta_{\mu_{1}\mu^{\prime}})(f_{\mu^{\prime}}^{\dagger}f_{\mu_{2}}-x\delta_{\mu_{2}\mu^{\prime}})\
.$ (67)

In the second equation, we have omitted an energy constant term from the anti-
commutator
$\\{d^{\dagger}_{\mu_{2}}(\epsilon_{2}),d_{\mu_{1}}(\epsilon_{1})\\}$.
$\mathbb{P}$ is the projector to the restricted Hilbert space, where
$\hat{N}_{f}=n_{f}$. It is omitted in the second equation because it commutes
with $f_{\mu^{\prime}}^{\dagger}f_{\mu_{2}}$ and
$f_{\mu_{1}}^{\dagger}f_{\mu^{\prime}}$. The four-fermion operator
$\sum_{\mu^{\prime}}f_{\mu_{1}}^{\dagger}f_{\mu^{\prime}}f_{\mu^{\prime}}^{\dagger}f_{\mu_{2}}$
can be simplified to $(\mathcal{N}-n_{f}+1)f_{\mu_{1}}^{\dagger}f_{\mu_{2}}$
as the Hilbert space is restricted to $\hat{N}_{f}=n_{f}$. Then the
corrections to $g$, $xg$ from Sec. B.1 can be read out as

$\frac{dg}{dt}\bigg{|}_{h}=\frac{4\Delta(-D(t))}{\pi
U}\left(\mathcal{N}-n_{f}+1-2x\right)g^{2}\ ,$ (68)
$\frac{d(xg)}{dt}\bigg{|}_{h}=\frac{4\Delta(-D(t))}{\pi
U}\left(-x^{2}\right)g^{2}\ .$ (69)

Adding up the particle and the hole contributions we can obtain the RG
equations for $g$ and $(xg)$. The Kondo energy scale $T_{K}$ can be estimated
as the reduced bandwidth $D(t)$ where $g$ diverges. For a constant
$\Delta(\omega)=\Delta_{0}$, we obtain

$\frac{dg}{dt}=\frac{4\Delta_{0}}{\pi
U}\mathcal{N}g^{2},\qquad\frac{d(xg)}{dt}=\frac{4\Delta_{0}}{\pi U}n_{f}g^{2}\
$ (70)

and the solution

$g(t)=\frac{g(0)}{1-g(0)\frac{4\Delta_{0}}{\pi U}\mathcal{N}\cdot t}\ ,$ (71)
$x(t)=x(0)\frac{g(0)}{g(t)}+\frac{n_{f}}{\mathcal{N}}\cdot\frac{g(t)-g(0)}{g(t)}\
.$ (72)

where $g(0)$ is the initial condition given in Eq. 61. $g(t)$ diverges at
$t_{K}=\frac{\pi U}{4\mathcal{N}g(0)\Delta_{0}}$, corresponding the Kondo
energy scale $De^{-t_{K}}=De^{-\frac{\pi U}{4\mathcal{N}g(0)\Delta_{0}}}$. As
$g(t)$ diverges as $t\to t_{K}$, the second term in $x(t)$ dominates and there
must be $x\to\frac{n_{f}}{\mathcal{N}}$. In other words, $x$ flows to the
occupation fraction.

###  B.2 Application to the symmetric state at CNP

We assume a symmetric state of the THF model at CNP and examine its Kondo
energy scale. The Hamiltonian for the bath and the impurity-bath hybridization
here are given by $H^{(c,\eta)}$ and $H^{(cf,\eta)}$ in Eq. 1, respectively.
The energy eigenvalues and wave-functions of $H^{(c,-)}$ are

$\begin{aligned}
\epsilon_{1}({\mathbf{k}})=&\frac{M}{2}+\sqrt{\frac{M^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{1}({\mathbf{k}})=&\frac{1}{\sqrt{2}}\begin{pmatrix}\sin\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}&\sin\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}&-\cos\frac{\theta_{\mathbf{k}}}{2}&-\cos\frac{\theta_{\mathbf{k}}}{2}\end{pmatrix}^{T}\end{aligned}\
,$ (73) $\begin{aligned}
\epsilon_{2}({\mathbf{k}})=&-\frac{M}{2}+\sqrt{\frac{M^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{2}({\mathbf{k}})=&\frac{1}{\sqrt{2}}\begin{pmatrix}\cos\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}&-\cos\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}&-\sin\frac{\theta_{\mathbf{k}}}{2}&\sin\frac{\theta_{\mathbf{k}}}{2}\end{pmatrix}^{T}\end{aligned}\
,$ (74) $\begin{aligned}
\epsilon_{3}({\mathbf{k}})=&\frac{M}{2}-\sqrt{\frac{M^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{3}({\mathbf{k}})=&\frac{1}{\sqrt{2}}\begin{pmatrix}\cos\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}&\cos\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}&\sin\frac{\theta_{\mathbf{k}}}{2}&\sin\frac{\theta_{\mathbf{k}}}{2}\end{pmatrix}^{T}\end{aligned}\
,$ (75) $\begin{aligned}
\epsilon_{4}({\mathbf{k}})=&\frac{M}{2}-\sqrt{\frac{M^{2}}{4}+v_{\star}^{2}{\mathbf{k}}^{2}}\\\
u_{4}({\mathbf{k}})=&\frac{1}{\sqrt{2}}\begin{pmatrix}\sin\frac{\theta_{\mathbf{k}}}{2}e^{-i\phi_{\mathbf{k}}}&\sin\frac{\theta_{\mathbf{k}}}{2}e^{i\phi_{\mathbf{k}}}&\cos\frac{\theta_{\mathbf{k}}}{2}&-\cos\frac{\theta_{\mathbf{k}}}{2}\end{pmatrix}^{T}\end{aligned}\
.$ (76)

Following the calculations in Sec. A.3, we obtain
$V^{(\eta)}_{n\alpha}({\mathbf{k}})$ for $\alpha=1,\eta=-$ as

$\displaystyle V_{11}({\mathbf{k}})=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\gamma
e^{i\phi_{\mathbf{k}}}-v_{\star}^{\prime}ke^{-2i\phi_{\mathbf{k}}}\right)\sin\frac{\theta_{\mathbf{k}}}{2}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}$
$\displaystyle V_{21}({\mathbf{k}})=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\gamma
e^{i\phi_{\mathbf{k}}}+v_{\star}^{\prime}ke^{-2i\phi_{\mathbf{k}}}\right)\cos\frac{\theta_{\mathbf{k}}}{2}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}$
$\displaystyle V_{31}({\mathbf{k}})=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\gamma
e^{i\phi_{\mathbf{k}}}-v_{\star}^{\prime}ke^{-2i\phi_{\mathbf{k}}}\right)\cos\frac{\theta_{\mathbf{k}}}{2}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}$
$\displaystyle V_{41}({\mathbf{k}})=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\gamma
e^{i\phi_{\mathbf{k}}}+v_{\star}^{\prime}ke^{-2i\phi_{\mathbf{k}}}\right)\sin\frac{\theta_{\mathbf{k}}}{2}e^{-\frac{\lambda^{2}{\mathbf{k}}^{2}}{2}}\
.$ (77)

where

$\theta_{k}=\arccos\frac{M/2}{\sqrt{M^{2}/4+v_{\star}^{2}k^{2}}}\ .$ (78)

Using Eq. 52 we obtain the hybridization function contributed by the fully
symmetric $c$-bands (Fig. 1(b)) for $\alpha=1$, $\eta=-$ is

$\displaystyle\Delta(\omega)=\frac{\Omega_{0}}{4v_{\star}^{2}}\bigg{[}\left||\omega|-\frac{M}{2}\right|\theta(|\omega|-M)\left(\gamma^{2}+v_{\star}^{\prime
2}k_{F1}^{2}\right)\sin^{2}\frac{\theta_{k_{F1}}}{2}$ $\displaystyle\times
e^{-k_{F1}^{2}\lambda^{2}}+\left||\omega|+\frac{M}{2}\right|\left(\gamma^{2}+v_{\star}^{\prime
2}k_{F2}^{2}\right)\cos^{2}\frac{\theta_{k_{F2}}}{2}e^{-k_{F2}^{2}\lambda^{2}}\bigg{]}\
,$ (79)

where

$k_{F1}=\frac{1}{v_{\star}}\sqrt{\left(|\omega|-M/{2}\right)^{2}-\left(M/{2}\right)^{2}},$
(80)
$k_{F2}=\frac{1}{v_{\star}}\sqrt{\left(|\omega|+M/{2}\right)^{2}-\left(M/{2}\right)^{2}},$
(81)

Due to the time-reversal symmetry and crystalline symmetries, as explained in
Sec. III.1, the hybridization functions for other $\alpha,\eta$ must be the
same. We should choose the initial cutoff $D=\frac{1}{2}U_{1}$ beyond which
the Schrieffer-Wolff transformation is invalid. For these states
$k_{F1,2}\lesssim\frac{U_{1}}{2v_{\star}}$ and hence $v_{\star}^{\prime
2}k_{F1,2}^{2}\lesssim 119.4\mathrm{meV}^{2}$, which is significantly smaller
than $\gamma^{2}\approx 612.6\mathrm{meV}^{2}$. The damping factors
$e^{-\lambda^{2}k_{F}{1,2}^{2}}\gtrsim 0.74$ are also large. Thus, in the
following we approximate $\Delta(\omega)$ ($|\omega|<U_{1}/2$) as

$\displaystyle\Delta(\omega)\approx$
$\displaystyle\frac{\Omega_{0}}{4v_{\star}^{2}}\bigg{[}\left||\omega|-\frac{M}{2}\right|\theta(|\omega|-M)\gamma^{2}\sin^{2}\frac{\theta_{k_{F1}}}{2}$
$\displaystyle\qquad+\left||\omega|+\frac{M}{2}\right|\gamma^{2}\cos^{2}\frac{\theta_{k_{F2}}}{2}\bigg{]}\
.$ (82)

Sec. B.2 could be further simplified making use of Eqs. 80, 81 and 78, which
yields

$\Delta(\omega)=\left\\{\begin{array}[]{cc}b|\omega|,&|\omega|>M\\\
\Delta_{0}\left(1+\frac{|\omega|}{M}\right),&|\omega|<M\end{array}\right.$
(83)

where $b=\frac{\Omega_{0}\gamma^{2}}{4v_{\star}^{2}}\approx
0.1290,\Delta_{0}=\frac{\Omega_{0}M\gamma^{2}}{8v_{\star}^{2}}=bM/2$.

We first consider the flat-band limit ($M=0$), where
$\Delta(\omega)=b|\omega|$. We also assume that there is no multiplet
splitting in the symmetric state such that the effective Anderson model should
be a U(8) theory with $n_{f}=4$. Naively applying the RG equations derived in
the last subsection gives

$\frac{d\widetilde{g}}{dt}=-\widetilde{g}+\frac{4bD}{\pi
U_{1}}\mathcal{N}\widetilde{g}^{2},$ (84)

where $\mathcal{N}=8$, $D=U_{1}/2$, $\widetilde{g}=ge^{-t}$. Due to the
particle-hole symmetry at CNP, the initial condition (Eq. 61) is
$\widetilde{g}(0)=1$. It seems that there would be an unstable fixed point
$\widetilde{g}_{*}=\frac{2\pi}{4\mathcal{N}b}$, the initial $\widetilde{g}$
below (above) which flows to zero (infinity). Using the actual parameters we
find $g_{*}\approx 1.52$, hence the system would not be in the Kondo phase.
This result differs from the standard case with a constant $\Delta(\omega)$,
where a positive $g$ always flows to infinity. Furthermore, a more careful RG
analysis Fritz and Vojta (2004) shows that the fixed point $\widetilde{g}_{*}$
does not really exist. It is a false result of the weak coupling expansion,
which fails for $\Delta(\omega)\sim|\omega|^{r}$ with $r>\frac{1}{2}$. Thus, a
$\Delta(\omega)\sim|\omega|$ bath does not have a strong coupling phase. This
conclusion is also consistent with numerical studies Chen and Jayaprakash
(1995); Gonzalez-Buxton and Ingersent (1998); Ingersent and Si (2002).

We then consider the case with $M\neq 0$. We use the value $M=3.697$meV, which
leads to $\Delta(0)\approx 0.239$meV. The RG process can be divided into two
stages: $D(t)=\frac{1}{2}U_{1}e^{-t}>M$ and $D(t)<M$. The boundary between the
two stages is $t_{1}=\ln\frac{U_{1}}{2M}$. The RG equation in the first stage
reads

$\frac{dg}{dt}=\frac{2\mathcal{N}b}{\pi}g^{2}e^{-t}\;\Rightarrow\;g(t)=\frac{1}{1-\frac{2\mathcal{N}b}{\pi}(1-e^{-t})}\
.$ (85)

Due to the particle-hole symmetry, the initial condition given by Eq. 61 is
$g(0)=1$. We have $g_{1}=g(t_{1})\approx 2.34$. The RG equation in the second
stage is given by

$\displaystyle\frac{dg}{dt}=$ $\displaystyle\frac{4\mathcal{N}\Delta_{0}}{\pi