eqs

$\begin{split}\BODY\end{split}$ (1)

††thanks: These authors contributed equally.††thanks: These authors
contributed equally.††thanks: These authors contributed equally.††thanks:
These authors contributed equally.

# Topological Heavy Fermion Principle For Flat (Narrow) Bands With
Concentrated Quantum Geometry

Jonah Herzog-Arbeitman  Department of Physics, Princeton University,
Princeton, New Jersey 08544, USA  Jiabin Yu  Department of Physics, Princeton
University, Princeton, New Jersey 08544, USA  Dumitru Călugăru∗ Department of
Physics, Princeton University, Princeton, New Jersey 08544, USA  Haoyu Hu
Donostia International Physics Center, P. Manuel de Lardizabal 4, 20018
Donostia-San Sebastian, Spain  Nicolas Regnault  Department of Physics,
Princeton University, Princeton, New Jersey 08544, USA  Laboratoire de
Physique de l’Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne
Université, Université Paris-Diderot, Sorbonne Paris Cité, 75005 Paris, France
Chaoxing Liu  Department of Physics, the Pennsylvania State University,
University Park, Pennsylvania 16802, USA  Sankar Das Sarma  Condensed Matter
Theory Center and Joint Quantum Institute, Department of Physics, University
of Maryland, College Park, Maryland 20742, USA  Oskar Vafek  National High
Magnetic Field Laboratory, Tallahassee, Florida, 32310, USA  Department of
Physics, Florida State University, Tallahassee, Florida 32306, USA  Piers
Coleman  Center for Materials Theory, Department of Physics and Astronomy,
Rutgers University, Piscataway, NJ 08854 USA  Department of Physics, Royal
Holloway University of London, Egham, Surrey TW20 0EX, UK  Alexei Tsvelik
Division of Condensed Matter Physics and Materials Science, Brookhaven
National Laboratory, Upton, NY 11973-5000, USA  Zhi-da Song  International
Center for Quantum Materials, School of Physics, Peking University, Beijing
100871, China  Hefei National Laboratory, Hefei 230088, China  Collaborative
Innovation Center of Quantum Matter, Beijing 100871, China  B. Andrei Bernevig
<EMAIL_ADDRESS>Department of Physics, Princeton University, Princeton,
New Jersey 08544, USA  Donostia International Physics Center, P. Manuel de
Lardizabal 4, 20018 Donostia-San Sebastian, Spain  IKERBASQUE, Basque
Foundation for Science, Bilbao, Spain

###### Abstract

We propose a general principle for the low-energy theory of narrow bands with
concentrated Berry curvature and Fubini-Study metric in the form of a map to
Anderson-“+” models composed of heavy fermions hybridizing and interacting
with semi-metallic modes. This map resolves the obstruction preventing
topological bands from being realized in a local Hamiltonian acting on the
low-energy degrees of freedom. The concentrated quantum geometry is reproduced
through band inversion with a dispersive semi-metal, leaving a nearly flat,
trivial band which becomes the heavy fermion. This representation is natural
when the narrow band is not energetically isolated on the scale of the
interaction and an enlarged Hilbert space is inescapable, but also provides
analytical insight into the projected-interaction limit. First exemplified in
twisted bilayer graphene (TBG), we extend it to (1) the twisted checkerboard,
which we find has a chiral symmetric stable anomaly that forbids a lattice
realization at all energies, and (2) the Lieb lattice with gapless flat bands,
where we show the heavy fermions can be obtained by minimizing a Euclidean
instanton action to saturate its BPS bound. The heavy fermion approach is
widely applicable and physically transparent: heavy electrons carry the strong
correlations and dispersive electrons carry the topology. This simple picture
unifies the dichotomous phenomena observed in TBG and points to connections
between moiré and stoichiometric materials.

##  I Introduction

Flat or narrow bands resulting from quantum interference can appear in any
dimension, in both stoichiometric and engineered (moire) materials, and
engender strongly correlated phases of matter [1, 2, 3]. In stoichiometric
compounds, exactly flat bands arise in local tight binding models on bipartite
lattices and their descendants[4, 5, 6, 7, 8], providing a principle for the
existence of flat or narrow band materials [9, 10, 11, 12, 13, 14]. Moreover,
flat bands arising from strong hoppings are likely to be (fragile) topological
if they are gapped from dispersive bands [5]. Topology, and more generally
quantum geometry, is a crucial feature for flat bands to host exotic
correlated states, since trivial flat bands can be formed from spatially
decoupled atomic limits. In moiré systems, flat bands typically arise in
continuum models away from the tight-binding limit [15, 16], and exhibit
emergent symmetries in the small twist angle limit which frequently render
them topological. In these cases, no exponentially localized Wannier
description which faithfully represents the symmetries of the flat bands is
possible. Since continuum Hamiltonians are infinite, analytic expressions for
the wavefunctions are difficult to find, a notable exception being the
idealized “chiral” limit [17, 18, 19, 20, 21] of twisted bilayer graphene
(TBG) which remains separated from experimentally relevant parameters.

A common characteristic of these flat bands (outside the chiral limit) is the
presence of concentrated, generically non-abelian, quantum geometry in
isolated regions of the Brillouin zone. The question then arises: what is the
interacting physics particular to this type of flat band system? There are
several limits, schematized in Fig. 1. If the interaction scale $U$ is weaker
than the narrow band bandwidth $M$, then a Fermi liquid (with its own
instabilities at Van Hove singularities) is expected and weak coupling
Hartree-Fock methods can be used, though analytic progress is dependent upon
expressions for the eigenstates which may not be available. Second, if the
interaction is much larger than the bandwidth $M$ but smaller than the gap to
dispersive bands $\gamma$, then one can project the interaction to the flat
bands and solve a strong coupling interacting problem with exact Slater
groundstates [22, 23, 24, 25, 26] in the repulsive case at integer filling, or
exact eta-pairing groundstates in the attractive case [27, 28, 29, 30, 31] in
the spin-0 sector. These exact results are also in agreement with Hartree-
Fock, and analytical insight again requires knowledge of the flat eigenstates
to obtain the form factors. Third, if the interaction $U$ is larger than both
the bandwidth $M$ and the gap to the dispersive bands $\gamma$ one cannot
perform any meaningful restriction and the problem is both strongly coupled
and not prone to projection. The $\bf{k}$-space formalism obfuscates any local
description. Realistic models of TBG, as well as gapless flat bands including
the Lieb lattice [4] and twisted trilayer graphene [32], and many
stoichiometric materials fall under this latter category.

Figure 1: Characteristic energy scales of a Fermi liquid $(a)$, a
projected/strong coupling ferromagnet $(b)$, and a heavy fermion model $(c)$.

In twisted bilayer graphene, detailed microscopic calculations[33, 34, 35, 36,
37] have shown that the topological flat bands can be mapped into a set of
trivial and highly localized heavy fermions ($f$-electrons) hybridized with
untwisted bilayer graphene $c$-electrons of mass $\pm M\approx 4$ meV via a
matrix element $\gamma\approx 25$meV. Their band structure matches that of the
Bistritzer-MacDonald (BM) model [15] almost perfectly within the range of the
Coulomb interaction $U\approx 40-50$ meV, obtained assuming a dielectric
constant $\epsilon=6$. $2M$ becomes the bandwidth of the narrow (almost flat)
bands and $\gamma$ is the gap between the flat bands and a set of dispersive
bands in the BM model. This is a faithful, symmetric, unitary transformation
of the low-energy bands of the model, where strain, relaxation, and other
effects can be included [38, 39, 40, 41, 42, 43]. The immediate benefit of
this transformation is that it identifies the heavy carriers of the
interaction $U$ —the $f$-electrons— and the itinerant carriers of the topology
— the $c$-electrons. The interaction can be rewritten in this basis, and one
obtains a generalized Anderson model consisting of a large Hubbard $U\approx
50$ meV on the $f$-electrons, as well as other interactions. This model is
suitable is suitable for DMFT implementation [44, 45, 46, 47, 48, 49, 50, 51],
and can match in a natural way a large amount of the experimental features
observed in a wide array of experiments. Furthermore, it provides for a Kondo
mechanism [52, 53] for screening of the repulsive $U$ as a potential pathway
to understand superconductivity — the core mystery of TBG.

Twisted bilayer and trilayer graphene[54] are not the only systems whose flat
bands can be mapped into a topological heavy fermion model. Recently, a
quadratic touching model that is related to the twisted checkerboard model
[55] has also been mapped to a heavy fermion problem [21]. Here we present a
wider perspective, arguing that a broad class of narrow/flat band models with
concentrated Berry curvature or quantum geometry admit such a heavy fermion
description. This entails a single particle Hamiltonian that separates the
bands into atomic $f$-electrons coupled to topological semi-metals, as well as
periodic Anderson-“+” interacting terms. This broad class includes gapped and
gapless bands in continuum or tight-binding models as we will exemplify.

We first review the emergence of heavy fermions in TBG, and give a simple two-
band example on the lattice of a Chern band hybridizing with a metal. We then
consider the twisted checkerboard lattice in the chiral limit [21], completing
the heavy fermion theory by (i) showing the stable anomaly that obstructs a
chiral-symmetric tight-binding model at all energies and (ii) specifying the
forms of the Anderson-“+” interacting terms. We also consider a family of
gapless flat bands on bipartite lattices whose most famous representative is
the Lieb lattice. There we examine the heavy fermion Wannier state
wavefunction, showing that it can be solved analytically by mapping the
Wannier spread localization problem to a Euclidean instanton action.

Finally, we present arguments that heavy fermion physics arises in interacting
narrow bands with concentrated (even divergent, in the gapless case) quantum
geometry when the interaction strength is larger than the gap and the narrow
band width. To this end, we provide simple estimates for when this regime is
reached, and also characterize flat bands that cannot be mapped to heavy
fermion models.

##  II Gapped Topological Bands

###  II.1 $C_{2z}T$ protected obstruction: Heavy Fermions in Twisted Bilayer
Graphene

Twisted bilayer graphene (TBG) is formed by stacking together two layers of
graphene with a slight twist [15] and exhibits unique electronic properties,
such as correlated insulators [56] and superconducting phases [57]. Its non-
interacting band structure, which can be described by the Bistritzer-MacDonald
(BM) model [15], exhibits eight flat bands with stable/fragile topology (in
the presence/absence of eigenstate particle-hole symmetry) [58, 59, 60] around
charge neutrality. Importantly, the nonabelian Berry curvature associated with
these flat bands is predominantly localized near the $\Gamma$ point [33], as
shown in Fig. 2 (a). This observation enables one to model the TBG as a
topological heavy fermion (THF) system, which consists of strongly correlated,
localized $f$-orbitals and topological, dispersive $c$-electrons [33]. The
concentrated Berry curvature allows over 95% of the states within the flat
bands to be accurately described by the $f$ orbitals as shown in Fig. 2 (b).
The non-trivial topology of the flat bands arises from the hybridization
between $f$ and topological $c$ electrons. Remarkably, this mapping also
enables one to analytically obtain the wavefunctions of the active TBG flat
bands in the zero-bandwidth limit. These flat bands exhibit a chiral symmetry
which differs from the conventional chiral symmetry of the BM model in the
unrealistic $w_{0}/w_{1}=0$ limit [61]. Hence the resulting wavefunctions (see
Section B.1) [33]

$\displaystyle\hat{f}^{\dagger}_{\mathbf{k},1,\eta,s}-\frac{\left(\gamma
e^{-\frac{|\mathbf{k}|^{2}\lambda^{2}}{2}}/v_{\star}\right)\hat{c}^{\dagger}_{\mathbf{k},3,\eta,s}}{\eta
k_{x}+ik_{y}},$ (2)

are a good approximation of the realistic TBG bands. Here
$\hat{d}^{\dagger}_{\mathbf{k},\zeta=+,\eta,s}$ is the electron operator in
the Chern band basis with Chern number $\zeta\cdot\eta$, valley
$\eta\in\\{+,-\\}$, spin $s\in\\{\uparrow,\downarrow\\}$, while
$\hat{f}^{\dagger}_{\mathbf{k},\alpha,\eta,s}$ and
$\hat{c}^{\dagger}_{\mathbf{k}+\mathbf{G},a,\eta,s}$ are the corresponding
electron operators of the $f$\- and $c$-fermions, respectively.
$\alpha\in\\{1,2\\},a\in\\{1,2,3,4\\}$ denote the orbital indices of $f$\- and
$c$-fermions respectively. $\gamma$ characterizes the $f$-$c$ hybridization
strength with a decay factor $\lambda$. $v_{\star}$ is the Fermi velocity of
non-interacting $c$ electrons. $\hat{d}^{\dagger}_{\mathbf{k},\zeta=+,\eta,s}$
electron operators admit analytical expressions in the plane-wave BM basis
[33, 35] using the form of the $f$-electron Wannier states and of the
$c$-electron basis. This allows for computing form factors of the
$\hat{d}^{\dagger}_{\mathbf{k},+,\eta,s}$ operators. Approximate analytic
formulas can be obtained as shown in Section B.1.2 and were used in Ref. [62]
to compute the electron-phonon interaction.

There are two interesting limits of this heavy fermion formalism. First, if
the Hubbard interactions $U_{1}$ of $f$-fermions is larger than the $f$-$c$
hybridization strength $|\gamma|$ (which is also the gap between flat active
bands and remote bands), we should include both flat and dispersive bands and
treat the system as an effective heavy fermion system with both $f$ and $c$
fermions. In the opposite limit, where $|\gamma|>U_{1}$, a flat-band-projected
approach [25, 23] is more suitable for the problem. However, even in this
limit, the THF model can still be employed to derive the analytical
expressions of the flat-band wavefunction and to calculate the excitation
spectrum as we discuss in Appendix. B.

Figure 2:  (a) The non-abelian Berry curvature, expressed in arbitrary units,
of the flat Chern band in TBG at the magic angle. (b) The $f$-electron orbital
weight within the active TBG bands at the magic angle.

###  II.2 Chern bands with localized Berry curvature

A similar heavy fermion construction can also be obtained for gapped Chern
bands with concentrated Berry curvature. To see this, we consider a simple
model with $s$ and $p_{z}$ orbitals located at the $1a$ position of a square
lattice (see Appendix. C) [63]. We use $\hat{f}^{\dagger}_{\mathbf{k},\sigma}$
and $\hat{c}^{\dagger}_{\mathbf{k},\sigma}$ to denote the creation operator
for the $s$ and $p_{z}$ orbitals, respectively, with $\mathbf{k}$ being the
momentum and $\sigma$, the spin index. The $p_{z}$ orbitals host a doubly-
degenerate quadratic band near the $\Gamma$ point with dispersion
$\epsilon_{\mathbf{k},p_{z}}\approx-\Delta\epsilon+\frac{t_{p}\absolutevalue{\mathbf{k}}^{2}}{2}$,
where $t_{p}$ denotes the nearest-neighbor hopping of the $p_{z}$ orbitals,
and $\Delta\epsilon$ controls the energy level of the $p_{z}$ orbitals. The
$s$ orbitals have vanishing kinetic energy, but hybridize with the $p_{z}$
orbitals.

The single-particle Hamiltonian can be written as

$\displaystyle
H_{0}=\sum_{|\mathbf{k}|<\Lambda_{c},\sigma}(\Delta\epsilon-\frac{t_{p}}{2}|k|^{2})\hat{c}^{\dagger}_{\mathbf{k},\sigma}\hat{c}_{\mathbf{k},\sigma}$
$\displaystyle+\frac{t_{sp}}{\sqrt{N}}\sum_{|\mathbf{k}|<\Lambda_{c},\mathbf{R},\sigma}\left[e^{-i\mathbf{k}\cdot\mathbf{R}}(\sigma
k_{y}+ik_{x})\hat{c}^{\dagger}_{\mathbf{k},\sigma}\hat{f}_{\mathbf{R},\sigma}+\text{h.c.}\right]\,.$
(3)

where $N$ is the number of unit cells, $t_{sp}$ characterizes the $f$-$c$
hybridization strength and $\Lambda_{c}$ is the momentum cutoff of $c$
electrons. We consider the limit where $|\Delta\epsilon|,|t_{sp}|\ll|t_{p}|$
and $\Delta\epsilon<0$. Then the system consists of a two-fold degenerate
narrow Chern band with bandwidth $\sim|\Delta\epsilon|$ and Chern numbers $\pm
1$, and dispersive bands that are mostly formed by the dispersive $p_{z}$
orbitals. The non-trivial topology of the narrow band comes from the band
inversion that happens near the $\Gamma$ point, where we can observe that the
Chern bands are mostly formed by $p_{z}$ orbitals near $\Gamma$ and $s$
orbitals away from it. The Berry curvature is also concentrated near the
$\Gamma$ point (see Appendix. C). Therefore, our simple model indeed realizes
a system with both narrow Chern bands and dispersive bands. We can already
observe that, this system is equivalent to a heavy fermion system where the
$s$ orbitals correspond to the localized $f$ electrons, while the $p_{z}$
orbitals represent dispersive $c$ electrons. The flatness of the bands comes
from the localized nature of $s$ orbitals ($f$ fermions) and the non-trivial
topology emerges from the hybridization between $s$ and $p_{z}$ orbitals.
After further introducing the on-site Hubbard interaction of the $s$ orbital
($U$) and $p$ orbital ($V$), along with the Coulomb repulsion between $s$ and
$p$ orbitals ($W$), the full Hamiltonian adapts a formula similar to the THF
Hamiltonian of the TBG ( Appendix. C). Therefore, we have established the
duality between a topological heavy fermion systems and a system that contains
narrow Chern bands with concentrated Berry curvatures. This duality has also
been demonstrated in a moiré system featuring Chern flat bands [64].

###  II.3 Chiral-Symmetry-Protected Wannier Obstruction: Heavy Fermions in the
Twisted Checkerboard

The stable anomaly giving rise to the topology of the TBG narrow bands is
protected by particle-hole symmetry for small twist angles [65]. We now
discuss a stable anomaly protected by another symmetry, chiral symmetry, which
is present in various models, including the twisted checkerboard model which
we now discuss. Given a generic 2D single-particle matrix Hamiltonian
$h(\boldsymbol{k})$ with chiral symmetry $\mathcal{C}$ obeying
$\\{\mathcal{C},h(\boldsymbol{k})\\}=0$ and
$\mathcal{C}^{2}=\mathcal{C}\mathcal{C}^{\dagger}=1$, there may still be an
obstruction to finding exponentially localized Wannier functions for an
isolated set of bands, even if the set of bands have zero total Chern number.
To see this is so, let us consider a generic isolated set of $2N$ chirally
symmetric bands: their energies are symmetric with respect to the zero energy
and their projection matrix $P(\boldsymbol{k})$ commutes with $\mathcal{C}$.
An odd number of chirally symmetric bands would necessarily have an exactly-
flat zero-energy band, and is not our focus here. At a fixed $\boldsymbol{k}$,
we can separate the Hilbert space for $P(\boldsymbol{k})$ into the chiral-even
and chiral-odd subspaces, each of which has the dimension $N$; in other words,
$P(\boldsymbol{k})=P_{+}(\boldsymbol{k})+P_{-}(\boldsymbol{k})$ with
$\mathcal{C}P_{\pm}(\boldsymbol{k})=\pm P_{\pm}(\boldsymbol{k})$. We then
define the Berry curvature and Chern numbers for $P_{\pm}(\boldsymbol{k})$ as
$F_{\pm}(\boldsymbol{k})$ and $\text{Ch}_{\pm}$, respectively, and the total
Chern number of the entire isolated set of $2N$ bands reads
$\text{Ch}=\text{Ch}_{+}+\text{Ch}_{-}$. Even if the total Chern number is
zero, we can still have nonzero chiral Chern number
$\text{Ch}_{+}-\text{Ch}_{-}$, which forbids the construction of the
exponentially-localized Wannier functions with strictly local $\mathcal{C}$ in
the real space. See details in Appendix. D. The nonzero chiral Chern number
$\text{Ch}_{+}-\text{Ch}_{-}$ can occur even if there is $C_{2z}\mathcal{T}$
symmetry, but a stable anomaly cannot be protected by $C_{2z}\mathcal{T}$
alone.

In practice, the nonzero value of $\text{Ch}_{+}-\text{Ch}_{-}$ can be
indicated by the total winding number of the chiral protected band-touching
points at zero energy. Chiral symmetry can protect stable band touching points
at zero energy in 2D [66], labeled $\boldsymbol{k}_{l}$. The robustness of the
band touching points is indicated by their nonzero chiral-protected winding
numbers, labelled by $W_{l}$. As elaborated in Appendix. D, we prove that

$\sum_{l}W_{l}=\text{Ch}_{+}-\text{Ch}_{-}\ .$ (4)

Therefore, if $\sum_{l}W_{l}\neq 0$, a Wannier obstruction exists for any
isolated chiral-symmetric set of $2N$ bands, indicating a stable anomaly
protected by chiral symmetry. Continuously deforming $h(\boldsymbol{k})$
cannot change $\sum_{l}W_{l}$, unless the continuous deformation involves
making $h(\boldsymbol{k})$ into a vanishing matrix. Using Eq. (4), we can
conclude that the single-valley spinless single-particle Hamiltonian of
twisted bilayer graphene in the chiral limit [65] has chiral-protected
anomaly, since it has two Dirac cones at zero energy whose total winding
number is 2.

The twisted checkerboard model in the chiral limit proposed in Ref. [67] has a
stable anomaly protected by chiral symmetry. The model consists of two layers
of checkerboard lattice with lattice constant $1$; each layer of the
checkerboard system has two atoms (A and B) in one unit cell. For one
individual layer, the model is parameterized by two hoppings $t$ and
$t^{\prime}$ (Fig. 3(a)). Its low-energy physics at half filling happens
around $M$ point with an effective Hamiltonian
$h_{0}(\boldsymbol{p})=tp_{x}p_{y}\sigma_{x}+t^{\prime}(p_{x}^{2}+p_{y}^{2})\sigma_{y}$
in the basis of
$(\widetilde{c}^{\dagger}_{\boldsymbol{p},+},\widetilde{c}^{\dagger}_{\boldsymbol{p},-})=(c^{\dagger}_{\boldsymbol{p},A}-\mathrm{i}c^{\dagger}_{\boldsymbol{p},B},-\mathrm{i}c^{\dagger}_{\boldsymbol{p},A}+c^{\dagger}_{\boldsymbol{p},B})/\sqrt{2}$.
The effective Hamiltonian carries a chiral winding number 2 with chiral
operator $\sigma_{z}$, and the quadratic touching is pinned by $C_{4z}$
symmetry of the system—the chiral winding number 2 only guarantees two linear
touching points without pinning their positions (see Appendix. E.) The moiré
model is built from low-energy models around $M$ from both layers, which reads

$h_{\boldsymbol{r}}=\left(\begin{matrix}h_{0}(-\mathrm{i}\nabla_{\boldsymbol{r}})&W(\boldsymbol{r})\\\
W^{\dagger}(\boldsymbol{r})&h_{0}(-\mathrm{i}\nabla_{\boldsymbol{r}})\end{matrix}\right)\otimes
s_{0}\ ,$ (5)

where the basis is
$\left(\begin{matrix}\psi^{\dagger}_{\boldsymbol{r},t},\psi^{\dagger}_{\boldsymbol{r},b}\end{matrix}\right)$
with

$\psi^{\dagger}_{\boldsymbol{r},l}=(\psi^{\dagger}_{\boldsymbol{r},l,+,\uparrow},\psi^{\dagger}_{\boldsymbol{r},l,+,\downarrow},\psi^{\dagger}_{\boldsymbol{r},l,-,\uparrow},\psi^{\dagger}_{\boldsymbol{r},l,-,\downarrow})\
,$ (6)

$l=t,b$ labels the layer,

$W(\boldsymbol{r})=2\sum_{i=1}^{2}T_{i}\cos(\boldsymbol{r}\cdot\boldsymbol{q}_{i})+2\sum_{i=1}^{4}T_{i}^{\prime}\cos(\boldsymbol{r}\cdot\boldsymbol{g}_{i})\
,$ (7)

$\boldsymbol{q}_{n}=C_{4z}^{n-1}\frac{k_{\theta}}{\sqrt{2}}(1,1)^{T}$ with
$n=1,2,3,4$ (see $\boldsymbol{q}_{1}$ in Fig. 3(c)),
$\boldsymbol{g}_{i}=\boldsymbol{b}_{M,i}+\boldsymbol{q}_{1}$ with $i=1,2$,
$\boldsymbol{g}_{i}=C_{4z}\boldsymbol{g}_{i-2}$ with $i=3,4$,
$k_{\theta}=2\sqrt{2}\pi\sin(\frac{\theta}{2})$,
$\boldsymbol{b}_{M,1}=\boldsymbol{q}_{1}+\boldsymbol{q}_{4}$,
$\boldsymbol{b}_{M,2}=\boldsymbol{q}_{1}+\boldsymbol{q}_{2}$. Just like the
original BM model [15], the intralayer potential in Eq. (5) is neglected for
simplicity. In practice, we choose $T$ and $T^{\prime}$ to have the following
forms which preserve chiral symmetry:
$T_{1}=\sigma_{z}T_{2}\sigma_{z}=w_{1}\sigma_{x}$ and
$T_{1}^{\prime}=T_{2}^{\prime}=w_{1}^{\prime}\sigma_{x}$ with
$T_{3}^{\prime}=T_{4}^{\prime}=\sigma_{z}T_{1}^{\prime}\sigma_{z}$. Besides
chiral symmetry, the moiré model has spin $\mathrm{SU}(2)$, spinless $C_{4z}$,
spinless $C_{2x}$, $\mathcal{T}$, moiré-lattice translations, and an effective
$m_{z}$ symmetry that exchange two layers. (See Appendix. E for detailed
expressions.) We do not have valleys here since the moiré model is constructed
from the monolayer M points which are time-reversal symmetric. If we set
$w_{1}=w_{1}^{\prime}=0$, we have two quadratic touching points (from two
layers) at the zero energies at $\Gamma_{M}$ and $M_{M}$ in the moiré
Brillouin zone, and each of it has chiral winding number 2, adding uo to a
total chiral winding number of $\sum_{l}W_{l}=4$. The nonzero total chiral
winding number indicates a stable anomaly, which holds for nonzero $w_{1}$ and
$w_{1}^{\prime}$ due to the presence of chiral symmetry.

The stable topology of the nearly flat bands with concentrated curvature
present in this model are particularly suited for the construction of
topological heavy fermion model [21]. To show this, we choose
$t=4/k_{\theta}^{2}$, $t^{\prime}=1.26/k_{\theta}^{2}$, $w_{1}=0.66$ and
$w_{1}^{\prime}=-0.4$, resulting in nearly flat bands with nontrivial topology
and concentrated chiral Berry curvature
$F_{+}(\boldsymbol{k})-F_{-}(\boldsymbol{k})$ as shown in Fig. 3(b-c). (The
topology can also be indicated by the symmetry eigenvalues as discussed in
Appendix. E.) The concentrated chiral Berry curvature allows us to construct
the exponentially-localized Wannier modes $f_{\boldsymbol{k},\alpha}$,
$\alpha=1,2$ by combining states of nearly flat bands near zero energy and the
states near $X_{M}$ and $Y_{M}$ from the remote bands (see Fig. 3(b)). The
resulting $f$ modes have the symmetry representations of $p_{x}-p_{y}$
orbitals at $1a$ position in the moiré unit cell. There are 6 low-energy
states at $X_{M}$ (or $Y_{M}$) per spin, and the construction of the $f$ modes
uses 2 of them; the remaining 4 states at $X_{M}$ (or $Y_{M}$) per spin give
rise to the itinerant $c$ modes. The resulting heavy fermion model consists of
three parts: the $f$-mode part with zero Hamiltonian, the $c$-mode part, and
the $f-c$ coupling. The $c$-mode part reads

$H_{c}=\sum_{\boldsymbol{p}}^{\Lambda}c^{\dagger}_{X_{M}+\boldsymbol{p}}h_{cc}^{X_{M}}(\boldsymbol{p})\otimes
s_{0}c_{X_{M}+\boldsymbol{p}}+\text{{\hbox{C_{4z}}} partner}\ ,$ (8)

where $c^{\dagger}_{X_{M}+\boldsymbol{p}}$
$=(c^{\dagger}_{X_{M}+\boldsymbol{p},1}$,
$c^{\dagger}_{X_{M}+\boldsymbol{p},2}$,
$c^{\dagger}_{X_{M}+\boldsymbol{p},3}$,
$c^{\dagger}_{X_{M}+\boldsymbol{p},4})$ $\otimes(\uparrow,\downarrow)$,

$h_{cc}^{X_{M}}(\boldsymbol{p})=\left(\begin{matrix}m\sigma_{x}&\mathrm{i}v_{2}\left[(p_{x}+p_{y})\sigma_{x}+(p_{x}-p_{y})\sigma_{y}\right]\\\
h.c.&0_{2\times 2}\ .\end{matrix}\right)$ (9)

to the linear order of $\boldsymbol{p}$, and the $C_{4z}$ partner term stands
for the $Y_{M}$ Hamiltonian:

$h_{cc}^{Y_{M}}(\boldsymbol{p})=h_{cc}^{X_{M}}(C_{4z}^{-1}\boldsymbol{p})\ .$
(10)

For the spinless part of the basis $c^{\dagger}_{X_{M}+\boldsymbol{p},\beta}$
with $\beta=1,...,4$, $\beta=1,2$ have opposite spinless-$C_{2z}$ eigenvalues
than $\beta=3,4$, explaining the linear-in-$\boldsymbol{p}$ terms that couple
them. Combined with all other symmetries that leave $X_{M}$ unchanged,
$\beta=1,2$ furnish a reducible representation of them while $\beta=3,4$
furnish an irreducible representation, explaining the appearance of the mass
term for $\beta=1,2$ but not for $\beta=3,4$ (see Appendix. E). The $f-c$
coupling reads

$\displaystyle\begin{split}H_{fc}&=\sum_{\boldsymbol{p}}^{\Lambda}f^{\dagger}_{\boldsymbol{p}}h_{fc}^{X_{M}}(\boldsymbol{p})c_{X_{M}+\boldsymbol{p}}e^{-\frac{|\boldsymbol{p}|^{2}\lambda^{2}}{2}}+(C_{4z}\text{
partner})\\\ &\quad+h.c.\ ,\end{split}$ (11)

where $\lambda=0.43a_{M}$ is square root of the Wannier spread of the $f$
modes with $a_{M}$ the moiré lattice constant, and

$h_{fc}^{X_{M}}(\boldsymbol{p})=\left(\begin{matrix}\mathrm{i}v_{1}(p_{x}\sigma_{x}+p_{y}\sigma_{y})&\gamma(\sigma_{x}+\sigma_{y})\end{matrix}\right)\
.$ (12)

From the tight-binding parameters and the Wannier states, we compute
$v_{1}=-0.52/k_{\theta}$, $\gamma=0.45$, $m=0.51$ and $v_{2}=1.53/k_{\theta}$.
The heavy fermion bands in Fig. 3(d) match the continuum model remarkably
well. At the interacting level, we project the Coulomb interaction to obtain
an Anderson-“+” model with nine interaction terms. The larges term is the
density-density Hubbard $U$ for $f$ modes which can be treated with Hartree
Fock and DMFT methods. Other important terms, which can all be well-
approximated in mean field, are the density-density interaction $V/W$ among
$c$/ $f-c$ modes, and the flavor-dependent interactions $J$ and $K$ between
$f$ and $c$. They are listed in Appendix. E.

Figure 3: (a) The monolayer checkerboard model parametrized by two hoppings
proposed in Ref. [67]. (b) The band structure (red) of the twisted
checkerboard model with $t=4$, $t^{\prime}=1.26$, $w_{1}=0.66$ and
$w_{1}^{\prime}=-0.4$. The blue dots mark the probability of the $f$ modes.
(c) The chiral Berry curvature $F_{+}(\boldsymbol{k})-F_{-}(\boldsymbol{k})$
for the isolated set of two bands around zero energy in (b). (d) The
comparison between the continuum-model band structure (red) and the heavy
fermion-model band structure (blue).

We specify that Ref. [21] studied the heavy fermion framework for a quadratic
band touching model with a perfectly flat band protected by chiral symmetry.
We have here complemented their analysis by uncovering the stable chiral-
symmetric anomaly and obtaining the interaction Hamiltonian.

##  III Gapless Flat Band Models

Figure 4: Divergent quantum geometry of the Lieb lattice. The band structure
$(a)$ exhibits an exactly flat band protected by the bipartite lattice
structure shown in $(b)$ where the $L/L^{\prime}$ sublattice is red/blue
respectively. A singular band crossing at the $M=(\pi,\pi)$ point results in
divergent Fubini-Study metric $(c)$ in the flat band, going as
$g(M+\mathbf{k})\to 1/|\mathbf{k}|^{2}$. Blue shading in $(a)$ shows the
weight of the heavy fermion wavefunction on the dispersive bands.

We now turn to a family of lattice models known as bipartite crystalline flat
bands. These models support both gapped topological flat bands as well as flat
bands which are degenerate with dispersive bands at isolated momenta, which
will be our focus here. To introduce this family of models, we study the Lieb
lattice (see Fig. 4a). In momentum space, the Hamiltonian can be written in
bipartite form

$\displaystyle h(\mathbf{k})$
$\displaystyle=2t\begin{pmatrix}0&S^{\dagger}(\mathbf{k})\\\
S(\mathbf{k})&0\end{pmatrix},\quad S(\mathbf{k})=\begin{pmatrix}\cos
k_{y}/2\\\ \cos k_{x}/2\end{pmatrix}\ .$ (13)

There exist $N_{L}-N_{L^{\prime}}$ flat bands in a general bipartite models if
$S(\mathbf{k})$ is a rectangular matrix of dimensions $N_{L}\times
N_{L^{\prime}}$, where $L$ and $L^{\prime}$ are the two sublattices of the
crystal with $N_{L}$ and $N_{L^{\prime}}$ orbitals per unit cell respectively.
The flat band wavefunction $U_{0}(\mathbf{k})=(\psi(\mathbf{k}),0)^{T}$ is
supported only on the $L$ sublattice and has exactly zero energy since there
always exists a vector $\psi(\mathbf{k})$ such that
$S^{\dagger}(\mathbf{k})\psi(\mathbf{k})=0$. In the Lieb lattice,

$\displaystyle\psi(\mathbf{k})=(-\cos\frac{k_{x}}{2},\cos\frac{k_{y}}{2})^{T}/\sqrt{\cos^{2}\frac{k_{x}}{2}+\cos^{2}\frac{k_{y}}{2}}\
.$ (14)

This construction can be generalized to encompass realistic materials [5],
identifying minimal models for flat and nearly bands away from zero energy.

At generic $\mathbf{k}$, the dimension of the null space of
$S^{\dagger}(\mathbf{k})$ is $N_{L}-N_{L^{\prime}}=2-1=1$, but at high-
symmetry points the nullity may be larger, indicating band touching points.
Plainly, Eq. (13) shows a band touching occurs at $\mathbf{k}=(\pi,\pi)$, the
$M$ point at which $S^{\dagger}(\mathbf{k})=0$, as seen in the band structure
in Fig. 4b. Its existence fundamentally obstructs a local, low-energy model of
the flat band alone, even though most of the flat band is clearly separated
from the dispersive bands. From Eq. (14), we see that

$\displaystyle\psi(\mathbf{k})\to(\cos\phi,-\sin\phi)$ (15)

has a singular limit at the gap closing point, where
$\mathbf{k}-M=|k|(\cos\phi,\sin\phi)$ and $|k|\to 0$. While the Berry
curvature in this model is strictly zero due to $C_{2}\mathcal{T}$ symmetry
rendering all eigenstates real, the quantum geometry of the model, measured by
the Fubini-Study metric, is small away from $M$ and is peaked (indeed,
diverges) at the gap closing, as shown in Fig. 4c. This suggests that a
topological heavy fermion model can overcome this obstruction.

The band structure in Fig. 4a shows a flat band with pierced by a dispersive
bands. What is the low-energy theory at partial filling of the flat bands?
Naively, one may expect a theory of Dirac-like electrons and coupled to the a
localized $s$ orbital in the atomic basis. Instead, we will show analytically
that the low-energy degrees of freedom are a quadratic semi-metal with $2\pi$
Berry phase (unlike the linear, $\pi$ Berry phase Dirac) and an obstructed
atomic $d$-orbital (unlike the atomic basis). The (unitary) transformation of
the microscopic degrees of freedom into this basis is accomplished through the
heavy fermion framework.

Symmetries play a key role in identifying the low-energy degrees of freedom of
the heavy fermion model. We will need the rotation symmetry
$D[C_{4z}]=1\oplus\sigma_{1}$ and time-reversal $D[\mathcal{T}]=K$, obeying
$D[g]h(\mathbf{k})D^{\dagger}[g]=h(g\mathbf{k})$, as well as the anti-
commuting chiral symmetry $D[\Sigma]=1\oplus-\sigma_{0}$ obeying
$D[\Sigma]h(\mathbf{k})D^{\dagger}[\Sigma]=-h(\mathbf{k})$. The combination of
$C_{2}$ and $\mathcal{T}$ ensures the Berry curvature is identically zero
everywhere in the Brillouin zone when the bands are gapped. The irrep tables
in this space group, $p41^{\prime}$, are [68]

$\displaystyle\begin{array}[]{c|ccc}41^{\prime}&1&C_{4}\\\
\hline\cr\Gamma_{1},M_{1}&1&1\\\ \Gamma_{2},M_{2}&1&-1\\\
\Gamma_{3}\Gamma_{4},M_{3}M_{4}&2&0\\\ \end{array},\quad\
\begin{array}[]{c|cc}21^{\prime}&1&C_{2}\\\ \hline\cr X_{1}&1&1\\\
X_{2}&1&-1\\\ \end{array}\ .$ (16)

The band touching at the $M$ point is protected by these symmetries. Ref. [5]
has shown that the irreps of the flat band are given by the formal subtraction
of the representations of the two sublattices

$\displaystyle\mathcal{B}_{0}$
$\displaystyle=\mathcal{B}_{L}\ominus\mathcal{B}_{L^{\prime}}=\Gamma_{2}+X_{1}+(M_{3}M_{4}\ominus
M_{2})$ (17)

where $\mathcal{B}_{L}=\Gamma_{1}\oplus\Gamma_{2}+X_{1}\oplus
X_{2}+M_{3}M_{4}$ and $\mathcal{B}_{L^{\prime}}=\Gamma_{1}+X_{2}+M_{2}$ can be
found on the Bilbao crystallographic server[68]. The incompatibility of the
$M$ irrep subtraction is interpreted as a gapless point. In fact, these
representations identify a unique atomic (localized) representation available
at low energy. At the $\Gamma,X$ points, the dispersive bands are far away in
energy from the flat band so the only accessible irreps are
$\Gamma_{2},X_{1}$. There is only a single one-band atomic representation
consistent with these irreps:

$\displaystyle\Gamma_{2}+X_{1}+M_{2}=B_{1a}\uparrow G$ (18)

where $B_{1a}$ is a $d_{x^{2}-y^{2}}$-orbital at the origin — importantly off
the sites of the $L$ and $L^{\prime}$ sublattices. This is the representation
of the low-energy heavy fermion. The $M_{2}$ irrep is induced from the
$L^{\prime}$ sublattice, meaning the wavefunction of the heavy fermion obeys
$U_{f}(\mathbf{k}=M)=(1,0,0)$. This requires mixing into the dispersive band
Hilbert space since the flat band wavefunction (carrying the
$\Gamma_{2},X_{1}$ irreps) is entirely supported on the $L$ sublattice. Thus
the heavy fermion wavefunction takes the simple form

$\displaystyle U_{f}(\mathbf{k})$
$\displaystyle=U_{0}(\mathbf{k})\cos\theta(\mathbf{k})+(1,0,0)^{T}\sin\theta(\mathbf{k})$
(19)

where the band mixing is parameterized by $\theta(\mathbf{k})$.
$\theta(M)=\pi/2$ and decays to zero away from $M$. We can choose
$\theta(\mathbf{k})$ such that the heavy fermion is maximally localized[69] in
position space in order to obtain an Anderson-“+” model where Hubbard-like
interactions are predominantly carried by the heavy fermion, and correlations
are local. Usually one first uses numerical routines such as Wannier90 to
accomplish this task, or relies on gaussian approximations for the Wannier
states [33]. Remarkably, we now show it is possible to solve for
$\theta(\mathbf{k})$ analytically. This is likely to generalize to other
bipartite models with band touchings, as the $L$ and $L^{\prime}$ support of
the flat band and touching point as well as ansatz in Eq. 19 generalize.

We form the Wannier state

$\displaystyle f^{\dagger}_{\mathbf{R}}$
$\displaystyle=\sum_{\alpha}\int\frac{d^{2}k}{(2\pi)^{2}}e^{i\mathbf{k}\cdot\mathbf{R}}c^{\dagger}_{\mathbf{k},\alpha}U_{f,\alpha}(\mathbf{k})$
(20)

and compute its spread to be (see Appendix. F)

$\displaystyle\frac{1}{2}\braket{r^{2}}$
$\displaystyle=\int\frac{d^{2}k}{(2\pi)^{2}}\left(\frac{1}{2}(\nabla\theta)^{2}+\frac{1}{2}g(\mathbf{k})\cos^{2}\theta\right)$
(21)

using the fact that a smooth $C_{2}\mathcal{T}$-symmetric band can be chosen
real with vanishing Berry connection. Here

$\displaystyle g(\mathbf{k})$ $\displaystyle=\frac{1}{2}\text{Tr}\,(\nabla
P)^{2}=\frac{1-\cos k_{x}\cos k_{y}}{2(2+\cos k_{x}+\cos k_{y})^{2}}$ (22)

is the Fubini-Study (quantum) metric of the flat band with projector
$P=U_{0}U_{0}^{\dagger}$, which diverges as $1/|k|^{2}$ for
$\mathbf{k}=M+|k|(\cos\phi,\sin\phi)$. Only through mixing into the dispersive
bands parameterized by $\sin^{2}\theta$, which reaches 1 at $M$, can the
Wannier spread of the heavy fermion be made finite. We now impose a cutoff
$\Lambda\ll\pi$ around the $M$ point, beyond which $\theta(\mathbf{k})=0$ and
the heavy fermion wavefunction is supported strictly in the flat bands. This
cutoff is analogous to the “energy window” employed in disentanglement
algorithms, and physically corresponds to allowing band mixing only within the
a range of energy given by the interaction. The cutoff $\Lambda\sim U/t$
restricts the heavy fermion to be supported on low-energy states only.
Minimizing Eq. (21) now corresponds to an Euler-Lagrange problem with
Dirichlet boundary conditions given by the action

$\displaystyle S$
$\displaystyle=\int_{\Lambda}\frac{d^{2}k}{(2\pi)^{2}}\left(\frac{1}{2}(\nabla\theta)^{2}+\frac{1}{2k^{2}}\cos^{2}\theta\right)$
(23)
$\displaystyle=\int_{\Lambda}\frac{dt}{(2\pi)^{2}}\left(\frac{1}{2}\dot{\theta}^{2}+\frac{1}{2}\cos^{2}\theta\right)$

where in the last line we have taken $\theta(\mathbf{k})=\theta(|\mathbf{k}|)$
since higher angular momentum dependence strictly increases the spread, and
written $k=e^{t}$ to bring the action into the famous instanton form. Writing
$\dot{\theta}^{2}+\cos^{2}\theta=(\dot{\theta}\pm\cos\theta)^{2}\mp
2\frac{d}{dt}\sin\theta$ allows us to saturate the
Bogomol‘nyi–Prasad–Sommerfield (BPS) bound [70], and find the absolute minimum
with solution (see Appendix. F)

$\displaystyle\cos\theta(\mathbf{k})=\frac{2|\mathbf{k}|/\Lambda}{(|\mathbf{k}|/\Lambda)^{2}+1}$
(24)

for $|\mathbf{k}|\leq\Lambda$ and $\theta(\mathbf{k})=0$ otherwise. Note that
$\theta(\mathbf{k})$ and its first derivative are continuous at $\Lambda$, as
can be seen in Fig. 5a. The real space Wannier wavefunction and its
$d$-orbital character are shown in Fig. 5b. The theory of instantons is well
developed and solutions are known for higher dimensional, non-abelian actions
(see Ref. [71] and references therein), which may find other analogues in
heavy-fermion problems.

The heavy fermion wavefunction carries the $M_{2}$ irrep, but the $M_{3}M_{4}$
irrep must be carried by another particle, the conduction electrons. Their
wavefunctions are easily obtained by looking at the two-dimensional orthogonal
complement of $U_{f}(\mathbf{k})$, and fixing the gauge freedom with time-
reversal symmetry (see Appendix. F). To leading order, we find

$\displaystyle\gamma^{\dagger}_{M+\mathbf{k},c_{\pm}}$
$\displaystyle=\sum_{\alpha}c^{\dagger}_{M+\mathbf{k},\alpha}\big{(}\mp
i\sqrt{2}\,\frac{k_{x}\pm
ik_{y}}{\Lambda},\mp\frac{i}{\sqrt{2}},-\frac{1}{\sqrt{2}}\big{)}_{\alpha}$
(25)

which are explicitly smooth, and have highest probability on the $L$
sublattice (recall at the heavy fermion wavefunction is supported primarily on
the $L^{\prime}$ sublattice within the cutoff). The Hamiltonian in the heavy
fermion basis can now be written in a $k.p$ expansion as ($k,\bar{k}=k_{x}\pm
ik_{y}$)

$\displaystyle
h_{HF}(\mathbf{k})=t\begin{pmatrix}&-2i\bar{k}^{2}/\Lambda&\bar{k}/\sqrt{2}\\\
2ik^{2}/\Lambda&&k/\sqrt{2}\\\ k/\sqrt{2}&\bar{k}/\sqrt{2}&0\\\ \end{pmatrix}$
(26)

in the basis $c_{\mathbf{k},+},c_{\mathbf{k},-},f_{\mathbf{k}}$. The heavy
fermion description reveals the emergent degrees of freedom at low energies to
be obstructed $d$-orbitals coupled to a double vortex topological semi-metal
(see Fig. 5c).

Figure 5: Lieb lattice Heavy Fermion. (a) $c$-electron instanton-tunneling
probability in the flat band. (b) Real-space Wannier function, whose
$d$-orbital character is apparent from the sign structure (blue/red
corresponds to $\pm$ amplitude). The total probability on the $L^{\prime}$
sublattice is set by the cutoff $\Lambda$. $(c)$ We show the decoupled HF and
conduction electron bands (dashed) within $\Lambda$ around the $M$ point. The
conduction electrons carry $2\pi$ Berry phase and are a topological semi-
metal. Linear $f-c$ coupling restores the Lieb dispersion (black).

We now consider the addition of Hubbard interactions

$\displaystyle H_{int}$
$\displaystyle=U\sum_{\mathbf{R}\alpha}n_{\mathbf{R}\alpha,\uparrow}n_{\mathbf{R}\alpha,\downarrow}$
(27)

acting on the spin-degenerate bands. Since the HF basis is unitarily related
to the band basis, it is simple to rotate between them and compute the
periodic Anderson model, which we show in Appendix. F. Let us recall that for
a highly localized heavy fermion like in TBG, the largest interaction is an
onsite Hubbard term, with nearest- and farther-neighbor interactions
controlled by the Wannier spread. A surprising feature that distinguishes the
Lieb lattice from its continuum heavy fermion counterparts is that there is an
obstruction to taking $\lambda\to 0$. This obstruction is easily explained
through topological quantum chemistry[72]: the representation of the Lieb
lattice $f$-mode is a $d$-orbital at a lattice vacancy, so its wavefunction
must always be supported on the neighboring orbitals (see Fig. 5) which causes
strong overlaps between the densities of neighboring $f$-modes. Hence the
$f$-mode density-density interaction always has strong nearest-neighbor
components, while farther-neighbor interactions will be suppressed due to the
localized Wannier wavefunction. Approximating the Wannier state by keeping
only its largest 4 amplitudes (the central 4 $d$-wave like lobes), we find the
$f$-$f$ interaction takes the form

$\displaystyle H_{U}$
$\displaystyle=U\sum_{\mathbf{R}\alpha}\bar{n}^{f}_{\mathbf{R},\alpha,\uparrow}\bar{n}^{f}_{\mathbf{R},\alpha,\downarrow}\approx\frac{U}{4}\sum_{\braket{\mathbf{R}\mathbf{R}^{\prime}}}n_{\braket{\mathbf{R}\mathbf{R}^{\prime}},\uparrow}n_{\braket{\mathbf{R}\mathbf{R}^{\prime}},\downarrow}$
(28)

where $\bar{n}^{f}_{\mathbf{R},\alpha,\sigma}$ is the projection of
$n_{\mathbf{R},\alpha,\sigma}$ to the $f$-modes,
$f^{\dagger}_{\braket{\mathbf{R}\mathbf{R}^{\prime}},\sigma}=\frac{f^{\dagger}_{\mathbf{R},\sigma}+f^{\dagger}_{\mathbf{R}^{\prime},\sigma}}{\sqrt{2}}$
is a bond-centered heavy fermion,
$n_{\braket{\mathbf{R}\mathbf{R}^{\prime}},\sigma}=f^{\dagger}_{\braket{\mathbf{R}\mathbf{R}^{\prime}},\sigma}f_{\braket{\mathbf{R}\mathbf{R}^{\prime}},\sigma}$
is its density, and $\braket{\mathbf{R}\mathbf{R}^{\prime}}$ denotes all
$f$-mode nearest neighbors on the square lattice. The physical meaning of the
symmetry obstruction to an onsite Hubbard $f$-mode term is the following:
while the tight-binding Hubbard interaction is onsite, it is smeared out in
the low-energy degrees of freedom due to the quantum geometry of the heavy
fermion. This shows that the effective model at low energies is dramatically
altered from the original tight-binding description.

We remark that $f$-mode sector $\bar{H}^{f}=H_{U}$ is an analytically
tractable strong coupling problem. This is because projected Hubbard
Hamiltonians with flat bands have exactly solvable groundstates [27] and low-
energy excitations [29], where the minimal Fubini-Study metric [73] is known
to generate the mass of the collective bosonic modes (the spin wave for
repulsive interactions, and the Cooper pair for attractive interactions) [74,
75, 76, 77, 78]. Projecting the Lieb lattice onto the $f$-mode basis yields
such a projected Hubbard Hamiltonian since the $f$-modes have zero kinetic
energy in Eq. (26). However, the hybridization with the conduction electrons
cannot be neglected, and the full interacting Hamiltonian

$\displaystyle H_{int}$ $\displaystyle=H_{U}+H_{W}+H_{J}+\dots$ (29)

contains additional interactions between the $f$ and $c$ modes, including the
$f$-density/$c$-density interaction $H_{W}\sim f^{\dagger}fc^{\dagger}c$ the
exchange interaction $H_{J}\sim f^{\dagger}cc^{\dagger}f$, which are usually
decoupled in mean-field. (The dots in Eq. 29 refer to odd-parity $f$-terms and
the conduction electron Coulomb repulsion which are typically neglected see
Appendix. F). Together, these terms form an Anderson-“+” model in the low-
energy Hilbert space, which can be approached from various directions. The
importance of quantum geometry in the Lieb lattice has been studied in seminal
work on flat band superconductivity [79], and reveals strongly correlated
behavior consistent with the phenomenology of heavy fermions [80]. Although
our focus is on the partially filled flat band for $U<t$, the heavy fermion
representation is simply a unitary rewriting of the model, and can potentially
shed light on other densities and parameter regimes of the model [81].

##  IV General Argument For Heavy Fermion Representations of Flat Bands With
Concentrated Quantum Geometry

We have shown that flat Bands with concentrated non-abelian Berry curvature in
moiré systems, or concentrated Fubini-Study metric in the Lieb lattice, are
accompanied by dispersive bands whose gap is smallest where the curvature or
metric is largest. These systems admit a heavy fermion representation and
serve to illustrate the broader proposal that narrow band systems with
strongly peaked quantum geometry are candidates for heavy fermion systems.

Let us recall the two obstacles to an exponentially localized Wannier basis
with physical symmetry representations in a given set of bands. The first
obstacle arises from symmetry considerations: if the momentum space irreps of
the band cannot be induced from atomic limits [72], then no such Wannier basis
with physical symmetry representations exists. The second obstacle is quantum
geometric: the Wannier spread is lower bounded by the integrated Fubini-Study
metric [69], becoming infinite if the bands carry a nonzero Chern number [82]
(see Ref. [83, 84]) or are gapless, as in e.g. the Lieb lattice. These
obstacles are intrinsic properties of the Hilbert space of the bands. They
cannot be removed without mixing into nearby bands, necessitating a heavy
fermion description.

Symmetries are crucial for restricting the possible $f$-modes. In all the
cases studied so far, the symmetries are powerful enough to uniquely identify
candidate $f$-modes: the $p_{x}-p_{y}$ orbitals in TBG and twisted
checkerboard, and the $d_{x^{2}-y^{2}}$ orbital in the Lieb lattice. It is
worth noting that a class of line graph lattices [85, 5, 9, 86] such as the
kagome are ruled out as heavy fermion candidates by symmetry.

The existence of irrep-compatible $f$-modes is a necessary but not sufficient
condition, since the quantum geometric obstruction must also be overcome. We
will now present a criterion to estimate the $f$-mode spread and determine
whether a heavy fermion description is appropriate.

Given a projector
$P_{\mathcal{S}}(\mathbf{k})=\sum_{n\in\mathcal{S}}U_{n}(\mathbf{k})U_{n}^{\dagger}(\mathbf{k})$
onto a set of bands indexed by $\mathcal{S}$, their abelian Fubini-Study
metric is [87]

$\displaystyle g_{\mathcal{S}}(\mathbf{k})$
$\displaystyle=\frac{1}{2}\text{Tr}\,\partial_{i}P_{\mathcal{S}}(\mathbf{k})\partial_{i}P_{\mathcal{S}}(\mathbf{k})\geq
0$ (30)

where the momentum index $i$ is implicitly summed over. Since
$P_{\mathcal{S}}(\mathbf{k})$ is gauge-invariant,
$g_{\mathcal{S}}(\mathbf{k})$ can be computed without choosing a gauge.
Picking $\mathcal{S}$ to be the set of narrow bands, a large value of
$g_{\mathcal{S}}(\mathbf{k})$ at its peak prevents a highly localized $f$-mode
by lower-bounding the spread of any possible Wannier functions formed in
$\mathcal{S}$. However, by including a set of low-energy dispersive bands
$\mathcal{D}$, it may be possible to find a more localized Wannier function if
$g_{\mathcal{S}\oplus\mathcal{D}}(\mathbf{k})<g_{\mathcal{S}}(\mathbf{k})$,
noting that
$P_{\mathcal{S}\oplus\mathcal{D}}=P_{\mathcal{S}}+P_{\mathcal{D}}$, but
$g_{\mathcal{S}\oplus\mathcal{D}}\neq g_{\mathcal{S}}+g_{\mathcal{D}}$. We
propose a simple estimate for the total Wannier spread $\lambda^{2}_{f}$ of
the $f$-mode

$\displaystyle\lambda^{2}_{f}\approx\Omega\int\frac{d^{2}k}{(2\pi)^{2}}\min\\{\,g_{\mathcal{S}}(\mathbf{k}),g_{\mathcal{S}\oplus\mathcal{D}}(\mathbf{k})\\}$
(31)

where $\Omega$ is the unit cell area. $\lambda^{2}_{f}$ reflects the
possibility of reducing the Wannier obstruction by mixing into the dispersive
bands if the enlarged Hilbert space $\mathcal{S}\oplus\mathcal{D}$ reduces the
Fubini-Study metric. Fig. 6 shows
$g_{\mathcal{S}},g_{\mathcal{S}\oplus\mathcal{D}}$ for TBG and the twisted
checkerboard model, demonstrating how the inclusion of the nearby dispersive
bands can remove the strong peak in $g_{\mathcal{S}}(\mathbf{k})$. For TBG,
this estimate gives $0.29a_{M}$ compared to the Wannier90 value $0.3a_{M}$,
and for twisted checkerboard, this estimate gives $0.45a_{M}$ compared to
$0.42a_{M}$. The good numerical agreement between these estimates and the
Wannier90 results is evidence that the localization properties of different
Hilbert spaces, captured by Eq. (31), control the existence of a heavy fermion
(assuming there is no symmetry obstruction). Further understanding can be
gleaned from the Lieb lattice, where $g_{\mathcal{S}}$ diverges at the gap
closing point whereas $g_{\mathcal{S}\oplus\mathcal{D}}$ vanishes, since
$\mathcal{S}\oplus\mathcal{D}$ is the entire single-particle Hilbert space
whose projector is the identity. One can verify from the value of the
instanton action at its BPS bound that Eq. (31) is exact up to
$\frac{\Lambda^{2}}{2\pi}\ll 1$. Although Eq. (31) is written in terms of the
Fubini-Study metric, the inequality $g(\mathbf{k})\geq|f(\mathbf{k})|$ [88,
89] ensures that small Fubini-Study metric also constrains the Berry
curvature, $f(\mathbf{k})$.

Armed with an approximate value for the Wannier spread $\lambda^{2}_{f}$, it
is also possible to estimate the onsite Hubbard interaction strength via[35]

$\displaystyle U\approx\frac{e^{2}}{\epsilon\lambda_{f}}$ (32)

assuming a Coulomb interaction (with gate-screening length much larger than
$\lambda_{f}$). Given the gap $\gamma$ to the dispersive bands, we can
determine the regime of the resulting Anderson model based on whether

$\displaystyle\frac{\gamma}{U}$ $\displaystyle=\frac{\gamma}{e^{2}/\epsilon
a}\frac{\lambda_{f}}{a}$ (33)

is less than 1 (heavy fermion regime) or greater than 1 (projected regime).
Here $a$ is the lattice constant, and $e^{2}/(\epsilon a)$ is the Coulomb
energy scale. In TBG, using a value of $\epsilon=6$ for the dielectric
constant gives $\frac{\gamma}{e^{2}/\epsilon a}\sim 1$. Due to the strong
localization of the $f$-modes, $\lambda_{f}/a\sim.3$ ensures
$\frac{\gamma}{U}$ is firmly less than 1, putting TBG into the heavy fermion
regime where local moment behavior is expected.

Figure 6: Resolution of the quantum geometric obstruction to localized Wannier
functions. The 2-band Fubini-Study metric for TBG $(a)$ and twisted
checkerboard $(c)$ are peaked in the regions (dashed red) where the nearest 4
dispersive bands are closest. The 6-band Fubini-Study metric for TBG $(b)$ and
twisted checkerboard $(d)$ shows no peak in these regions, demonstrating that
mixing into the dispersive bands (within the low-energy, unshaded regions) can
yield highly localized Wannier functions.

##  V Conclusions

The low-energy description of a system can be very different from its
microscopic model. Untangling these relevant degrees of freedom from the
microscopic model is often the key to understanding their emergent properties,
both intuitively and computationally. The formation of heavy fermions from the
moiré coupling of two massless Dirac particles in TBG is a prototypical
example of this phenomenon in interacting flat bands. In this work, we have
studied a variety of other examples, from moiré systems to tight-binding
lattices, and observed a simple pattern. While band theoretic considerations
would forbid a local description of topological bands, we find that — if their
quantum geometry is concentrated — a heavy fermion representation is possible
by introducing metallic degrees of freedom to carry the obstruction. This
single-particle description is extended to a periodic Anderson-“+” model
incorporating interactions within this basis. The abundance of heavy fermions
thus far revealed in topological, correlated quantum materials underscores the
value in studying this class of models, and sets a direction for future
research.

##  VI Acknowledgements

The authors are grateful for important conversations and collaborations with
Antoine Georges, Xi Dai, Andy Millis, Francisco Guinea, Roser Valentí, Giorgio
Sangiovanni, Tim Wehling, Qimiao Si, Jian Kang, Dmitri K. Efetov, Keshav
Singh, Liam Lau, Daniel Kaplan, Ryan Lee, Gautam Rai, Lorenzo Crippa, and
Michael Scheer. This work was supported by Office of Basic Energy Sciences,
Material Sciences and Engineering Division, U.S. Department of Energy (DOE)
under Contracts No. DE-SC0016239 (B. A. B.), DE-FG02-99ER45790 (P. C.) and DE-
SC0012704 (A. M. T.). J. H.-A. is supported by a Hertz Fellowship, with
additional support from DOE Grant No. DE-SC0016239. J. Y. is supported by the
Gordon and Betty Moore Foundation through Grant No. GBMF8685 towards the
Princeton theory program and through the Gordon and Betty Moore Foundation’s
EPiQS Initiative (Grant No. GBMF11070). D.C. acknowledges the hospitality of
the Donostia International Physics Center, at which this work was carried out,
and acknowledges support by the Simons Investigator Grant No. 404513, the
Gordon and Betty Moore Foundation through Grant No. GBMF8685 towards the
Princeton theory program, the Gordon and Betty Moore Foundation’s EPiQS
Initiative (Grant No. GBMF11070), Office of Naval Research (ONR Grant No.
N00014-20-1-2303), Global Collaborative Network Grant at Princeton University,
BSF Israel US foundation No. 2018226, NSF-MERSEC (Grant No. MERSEC DMR
2011750). H.H. was supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (Grant Agreement
No. 101020833).

###### Contents

  1. I Introduction
  2. II Gapped Topological Bands
    1. II.1 $C_{2z}T$ protected obstruction: Heavy Fermions in Twisted Bilayer Graphene
    2. II.2 Chern bands with localized Berry curvature
    3. II.3 Chiral-Symmetry-Protected Wannier Obstruction: Heavy Fermions in the Twisted Checkerboard
  3. III Gapless Flat Band Models
  4. IV General Argument For Heavy Fermion Representations of Flat Bands With Concentrated Quantum Geometry
  5. V Conclusions
  6. VI Acknowledgements
  7. A Twisted bilayer graphene
    1. A.1 Review of the Bistritzer-MacDonald model
      1. A.1.1 Single-particle Hamiltonian
      2. A.1.2 Many-body Hamiltonian
    2. A.2 Review of the topological heavy fermion model for TBG
      1. A.2.1 The THF model fermions
      2. A.2.2 The single-particle THF model
      3. A.2.3 The interaction THF Hamiltonian
      4. A.2.4 Berry Curvature of the TBG active bands
  8. B Analytical results within the THF model
    1. B.1 Flat-band wave functions
      1. B.1.1 Analytical expression for the TBG flat-band wave functions
      2. B.1.2 Analytical expressions for the TBG form factors
    2. B.2 The charge-one excitation spectrum of TBG in the $\absolutevalue{\gamma}>U_{1}$ limit
      1. B.2.1 Charge-one excitation matrices
      2. B.2.2 Analytical expressions for the charge-one excitation matrices of TBG
      3. B.2.3 Flat-band projected limit of the THF model
  9. C Chern bands with localized Berry curvature
  10. D Chiral-Symmetry-Protected Stable Anomaly
  11. E Twisted Checkerboard Model
    1. E.1 Chiral-Symmetry Protected Stable Anomaly
    2. E.2 Heavy Fermion Model: Single-Particle
    3. E.3 Heavy Fermion Model: Interaction
      1. E.3.1 $H_{int,1}$
      2. E.3.2 $H_{int,2}$
      3. E.3.3 $H_{int,3}$
      4. E.3.4 $H_{int,4}$
      5. E.3.5 $H_{int,5}$
      6. E.3.6 $H_{int,6}$
      7. E.3.7 $H_{int,7}$
      8. E.3.8 Summary of The Interaction
    4. E.4 Alternative Convention
  12. F Lieb Lattice and Analytical Wannier States from the Instanton Action
    1. F.1 Conduction electrons and Hamiltonian
    2. F.2 Anderson Model

##  Appendix A Twisted bilayer graphene

Twisted bilayer graphene (TBG) consist of two graphene layers rotated by a
small angle $\theta$ relative to each other. Near the magic angle
$\theta=$1.05\text{\,}\mathrm{\SIUnitSymbolDegree}$$, flat bands emerge near
charge neutrality, which gives rise to various exotic quantum phases,
including correlated insulators and superconductivity. Extensive experimental
efforts, including both transport [90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114,
115, 116, 117] and spectroscopic measurements [118, 119, 120, 121, 122, 123,
124, 125, 126, 127, 128, 129, 130, 131, 132], have been taken to investigate
the nature of these quantum phases. Concurrently, significant theoretical
efforts have also been made towards constructing the microscopic models [133,
134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148,
149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163,
164, 59, 165, 166, 167, 168, 169], investigating experimental response [170,
171, 172, 173, 174, 120, 175, 176, 177, 178], understanding the topology [153,
179, 180, 181, 182, 183, 162, 184, 59] and excitations [185, 186, 187, 188,
189, 24, 190, 191], elucidating the mechanism behind the correlated insulating
states [192, 172, 193, 194, 195, 196, 197, 198, 199, 200, 120, 141, 201, 202,
203, 204, 205, 206, 207, 208, 209, 187, 210, 211, 212, 213, 214, 215, 216,
189, 23, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230,
231, 232, 233, 234, 190, 235, 236, 237, 154, 238, 191, 22, 239, 240, 241], and
superconductivity [242, 243, 244, 245, 246, 186, 247, 248, 249, 250, 251, 252,
253, 254, 255, 256, 257, 218, 220, 221, 227, 258, 259, 230, 260, 231, 233,
261, 235, 262, 183, 263, 264, 265, 266, 267, 241, 268]. In this appendix, we
briefly review both the Bistritzer-MacDonald (BM) model [15] and the
topological heavy fermion (THF) model [33] of TBG.

###  A.1 Review of the Bistritzer-MacDonald model

In this section, we briefly review the Bistritzer-MacDonald (BM) model [15] of
twisted bilayer graphene (TBG). Throughout this work, we use the notations
introduced in Refs. [149, 60, 269, 23, 24, 270, 271, 272]. We use
$\hat{a}^{\dagger}_{\mathbf{p},\alpha,s,l}$ to denote the creation operator of
an electron with momentum $\mathbf{p}$, spin $s=\uparrow,\downarrow$,
sublattice index $\alpha=A,B$, and within the graphene layer $l=\pm$. We use
$\mathbf{K}_{l}$ to label the $K$ point of the $l$-layer graphene Brillouin
Zone (BZ), and use $\eta=\pm$ to denote the two valleys of graphene which
corresponds to the momenta $\eta\mathbf{K}_{\pm}$. In addition, we introduce
the following auxiliary vectors

$\mathbf{q}_{1}=k_{\theta}\left(0,-1\right)^{T},\qquad\mathbf{q}_{2}=k_{\theta}\left(\frac{\sqrt{3}}{2},\frac{1}{2}\right)^{T},\qquad\mathbf{q}_{3}=k_{\theta}\left(-\frac{\sqrt{3}}{2},\frac{1}{2}\right)^{T},$
(34)

where
$k_{\theta}=\absolutevalue{\mathbf{K}_{+}-\mathbf{K}_{-}}=2\absolutevalue{\mathbf{K}_{+}}\sin\frac{\theta}{2}$.
The reciprocal lattice of the moiré system is then given by

$\displaystyle\mathcal{Q}_{0}=\mathbb{Z}\mathbf{b}_{M1}+\mathbb{Z}\mathbf{b}_{M2}$
$\displaystyle\mathbf{b}_{M1}=\mathbf{q}_{3}-\mathbf{q}_{1}\
,\qquad\mathbf{b}_{M2}=\mathbf{q}_{3}-\mathbf{q}_{2},$ (35)

For future convenience, we also introduce
$\mathcal{Q}_{\pm}=\mathcal{Q}\pm\mathbf{q}_{1}$. The corresponding moiré
lattice vectors in the real space are defined as
$\mathbf{a}_{M1}=\frac{2\pi}{3k_{\theta}}\left(\sqrt{3},1\right)$ and
$\mathbf{a}_{M2}=\frac{2\pi}{3k_{\theta}}\left(-\sqrt{3},1\right)$.

####  A.1.1 Single-particle Hamiltonian

To define the Hamiltonian of TBG, we consider the following electron operators

$\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}\equiv\hat{a}^{\dagger}_{\eta\mathbf{K}_{l}+\mathbf{k}-\mathbf{Q},\alpha,s,l},\mbox{\quad
for\quad}\mathbf{Q}\in\mathcal{Q}_{\eta l}.$ (36)

The single-particle Hamiltonian can then be written as [149, 60, 269]

$\hat{H}^{\text{TBG}}_{0}=\sum_{\mathbf{k}}\sum_{\eta,\alpha,\beta,s}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{\pm}}\left[h^{\text{TBG},\eta}_{\mathbf{Q},\mathbf{Q}^{\prime}}\left(\mathbf{k}\right)\right]_{\alpha\beta}\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\eta,\alpha,s}\hat{c}_{\mathbf{k},\mathbf{Q}^{\prime},\eta,\beta,s},$
(37)

We can decompose matrix
$h^{\text{TBG},\eta}_{\mathbf{Q},\mathbf{Q}^{\prime}}\left(\mathbf{k}\right)$
into two parts, which correspond, respectively, to the original Dirac cone
graphene Hamiltonian $h^{\text{D},\eta}_{\mathbf{Q}}\left(\mathbf{k}\right)$
and the inter-layer coupling term
$h^{\text{I},\eta}_{\mathbf{Q},\mathbf{Q}^{\prime}}$

$h^{\text{TBG},\eta}_{\mathbf{Q},\mathbf{Q}^{\prime}}\left(\mathbf{k}\right)=h^{\text{D},\eta}_{\mathbf{Q}}\left(\mathbf{k}\right)\delta_{\mathbf{Q},\mathbf{Q}^{\prime}}+h^{\text{I},\eta}_{\mathbf{Q},\mathbf{Q}^{\prime}}.$
(38)

For valley $+$, we have

$\displaystyle h^{\text{D},+}_{\mathbf{Q}}\left(\mathbf{k}\right)$
$\displaystyle=v_{F}\left(\mathbf{k}-\mathbf{Q}\right)\cdot\boldsymbol{\sigma},$
(39) $\displaystyle h^{\text{I},+}_{\mathbf{Q},\mathbf{Q}^{\prime}}$
$\displaystyle=\sum_{j=1}^{3}T_{j}\delta_{\mathbf{Q},\mathbf{Q}^{\prime}\pm\mathbf{q}_{j}}.$
(40)

where $v_{F}$ is the Fermi velocity of graphene,
$\boldsymbol{\sigma}=\left(\sigma_{x},\sigma_{y}\right)$ is the Pauli vector,
and $T_{j}$ characterizes the inter-layer tunneling with

$T_{j}=w_{0}\sigma_{0}+w_{1}\left[\sigma_{x}\cos\frac{2\pi\left(j-1\right)}{3}+\sigma_{y}\sin\frac{2\pi\left(j-1\right)}{3}\right],\quad\text{for}\quad
j=1,2,3.$ (41)

In this work, we take $w_{1}=$110\text{\,}\mathrm{meV}$$,
$v_{F}=$5.944\text{\,}\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$,
$\absolutevalue{\mathbf{K}_{+}}=$1.703\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{-1}$$,
and vary $0\leq w_{0}/w_{1}\leq 1$. The Hamiltonian of valley $\eta=-$ can be
obtained from the time-reversal transformation

$h^{\text{TBG},-}_{\mathbf{Q},\mathbf{Q}^{\prime}}\left(\mathbf{k}\right)=\sigma_{x}h^{\text{TBG},+}_{-\mathbf{Q},-\mathbf{Q}^{\prime}}\left(-\mathbf{k}\right)\sigma_{x}.$
(42)

####  A.1.2 Many-body Hamiltonian

For a double-gated experimental setup, the electron-electron interaction
potential is given by [269, 23]

$V\left(\mathbf{r}\right)=U_{\xi}\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{\sqrt{\left(\absolutevalue{\mathbf{r}}/\xi\right)^{2}+n^{2}}}\mbox{\quad
and\quad}V\left(\mathbf{q}\right)=\left(\pi
U_{\xi}\xi^{2}\right)\frac{\tanh\left(\absolutevalue{\mathbf{q}}\xi/2\right)}{\absolutevalue{\mathbf{q}}\xi/2},$
(43)

where
$U_{\xi}=\frac{e^{2}}{4\pi\epsilon_{0}\epsilon\xi}=$24\text{\,}\mathrm{meV}$$,
$\epsilon=6$ is the dielectric constant, and $\xi=$10\text{\,}\mathrm{nm}$$ is
the distance between the two screening gates. The interaction term takes the
form of [22, 25, 269]

$\hat{H}^{\text{TBG}}_{I}=\frac{1}{2}\frac{1}{N_{0}\Omega_{0}}\sum_{\mathbf{q}}\sum_{\mathbf{G}\in\mathcal{Q}_{0}}V\left(\mathbf{q}+\mathbf{G}\right)\delta\rho\left(-\mathbf{q}-\mathbf{G}\right)\delta\rho\left(\mathbf{q}+\mathbf{G}\right),$
(44)

where the normal-ordered density operator
$\delta\rho\left(\mathbf{q}+\mathbf{G}\right)$ is defined as

$\delta\rho\left(\mathbf{q}+\mathbf{G}\right)=\rho\left(\mathbf{q}+\mathbf{G}\right)-\frac{1}{2}\delta_{\mathbf{q},\mathbf{0}}\delta_{\mathbf{G},\mathbf{0}}=\sum_{\eta,\alpha,s}\sum_{\mathbf{k}}\sum_{\mathbf{Q}\in\mathcal{Q}_{\pm}}\left(\hat{c}^{\dagger}_{\mathbf{k}+\mathbf{q},\mathbf{Q}-\mathbf{G},\alpha,\eta,s}\hat{c}_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}-\frac{1}{2}\delta_{\mathbf{q},\mathbf{0}}\delta_{\mathbf{G},\mathbf{0}}\right),$
(45)

$N_{0}$ is the number of moiré unit cells, and $\Omega_{0}$ is the area of a
single moiré unit cell. The full TBG Hamiltonian is written as
$\hat{H}^{\text{TBG}}=\hat{H}^{\text{TBG}}_{0}+\hat{H}^{\text{TBG}}_{I}$.

###  A.2 Review of the topological heavy fermion model for TBG

We now review the topological heavy fermion (THF) description [33, 273, 35,
274, 275, 276, 277, 50, 45, 49, 278, 171] of the BM low-energy bands. The THF
model resolves the stable topological obstruction of the entire BM model [60,
33, 58, 25, 279, 280] by introducing two fermionic species: “heavy” ($f$) and
“conduction” ($c$) electrons. In contrast to the $f$-fermions, which form
effective $p_{x}\pm ip_{y}$ orbitals localized at the moiré AA sites, the
$c$-fermions are semi-metallic, dispersive, and form anomalous electronic
bands.

####  A.2.1 The THF model fermions

In terms of the low-energy operators of the BM model introduced in Eq. 36, the
$f$-fermions are given in momentum space by

$\hat{f}^{\dagger}_{\mathbf{k},\alpha,\eta,s}=\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;\alpha}\left(\mathbf{k}\right)\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s},$
(46)

where $v^{\eta}_{\mathbf{Q}\beta;\alpha}\left(\mathbf{k}\right)$ is the
$f$-fermion momentum-space wave function in valley $\eta$, with $\alpha=1,2$
denotes the orbital quantum number. Although usually obtained numerically
through a Wannierization and disentanglement procedure, the $f$-fermion wave
function has an analytical approximation in terms of Gaussian profiles [33,
35]. Defining the Fourier transformation for the $f$-fermion wave functions,

$\displaystyle v^{\eta}_{\mathbf{Q}\beta;\alpha}\left(\mathbf{k}\right)$
$\displaystyle=\frac{1}{\sqrt{\Omega_{0}}}\int\text{d}^{2}{r}w^{\eta}_{l\beta;\alpha}\left(\mathbf{r}\right)e^{-i\left(\mathbf{k}-\mathbf{Q}\right)\cdot\mathbf{r}},\mbox{\quad
for\quad}\mathbf{Q}\in\mathcal{Q}_{\eta l},$ (47) $\displaystyle
w^{\eta}_{l\beta;\alpha}\left(\mathbf{r}\right)$
$\displaystyle=\frac{1}{N_{0}\sqrt{\Omega_{0}}}\sum_{\mathbf{k}}\sum_{\mathbf{Q}\in\mathcal{Q}_{\eta
l},\beta}v^{\eta}_{\mathbf{Q}\beta;\alpha}\left(\mathbf{k}\right)e^{i\left(\mathbf{k}-\mathbf{Q}\right)\cdot\mathbf{r}},$
(48)

their real-space Wannier wave functions can be approximated by [33],

$\displaystyle w_{l1;1}^{\eta}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{\alpha_{1}}{\sqrt{2\pi\lambda_{1}^{2}}}e^{i\frac{\pi}{4}l\eta-\frac{r^{2}}{2\lambda_{1}^{2}}},\qquad$
$\displaystyle w_{l2;1}^{\eta}(\mathbf{r})$ $\displaystyle=$
$\displaystyle-l\frac{\alpha_{2}}{\sqrt{2}}\frac{r_{x}+i\eta
r_{y}}{\lambda_{2}^{2}\sqrt{\pi}}e^{i\frac{\pi}{4}l\eta-\frac{r^{2}}{2\lambda_{2}^{2}}},$
(49) $\displaystyle w_{l1;2}^{\eta}(\mathbf{r})$ $\displaystyle=$
$\displaystyle l\frac{\alpha_{2}}{\sqrt{2}}\frac{r_{x}-i\eta
r_{y}}{\lambda_{2}^{2}\sqrt{\pi}}e^{-i\frac{\pi}{4}l\eta-\frac{r^{2}}{2\lambda_{2}^{2}}},\qquad$
$\displaystyle w_{l2;2}^{\eta}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{\alpha_{1}}{\sqrt{2\pi\lambda_{1}^{2}}}e^{-i\frac{\pi}{4}l\eta-\frac{r^{2}}{2\lambda_{1}^{2}}},$
(50)

where the amplitudes $\alpha_{1}$ and $\alpha_{2}$, as well as the spreads
$\lambda_{1}$ and $\lambda_{2}$ have been obtained both numerically and
analytically across a large parameter regime [33, 35]. In momentum-spaceacross
a large parameter regime [4, 12]. In momentum-space the approximations from
Eqs. (A16) and (A17) can be, the approximations from Eqs. 49 and 50 can be
computed using Eq. 48

$\displaystyle v_{\mathbf{Q}1;1}^{(\eta)}(\mathbf{k})=$
$\displaystyle\alpha_{1}\sqrt{\frac{2\pi\lambda_{1}^{2}}{\Omega_{0}\mathcal{N}_{f,\mathbf{k}}}}e^{i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}},$
$\displaystyle v_{\mathbf{Q}2;1}^{(\eta)}(\mathbf{k})=$
$\displaystyle\alpha_{2}\sqrt{\frac{2\pi\lambda_{2}^{4}}{\Omega_{0}\mathcal{N}_{f,\mathbf{k}}}}\zeta_{\mathbf{Q}}\left[i\eta(k_{x}-Q_{x})-(k_{y}-Q_{y})\right]e^{i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}(\mathbf{k}-\mathbf{Q})^{2}\lambda_{2}^{2}},$
(51) $\displaystyle v_{\mathbf{Q}1;2}^{(\eta)}(\mathbf{k})=$
$\displaystyle\alpha_{2}\sqrt{\frac{2\pi\lambda_{2}^{4}}{\Omega_{0}\mathcal{N}_{f,\mathbf{k}}}}\zeta_{\mathbf{Q}}\left[-i\eta(k_{x}-Q_{x})-(k_{y}-Q_{y})\right]e^{-i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}(\mathbf{k}-\mathbf{Q})^{2}\lambda_{2}^{2}},$
$\displaystyle v_{\mathbf{Q}2;2}^{(\eta)}(\mathbf{k})=$
$\displaystyle\alpha_{1}\sqrt{\frac{2\pi\lambda_{1}^{2}}{\Omega_{0}\mathcal{N}_{f,\mathbf{k}}}}e^{-i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}}.$
(52)

In Eqs. 51 and 52, we have introduced the normalization factor
$\mathcal{N}_{f,\mathbf{k}}$, which can be determined to be

$\mathcal{N}_{f,\mathbf{k}}=\alpha_{1}^{2}\frac{2\pi\lambda_{1}^{2}}{\Omega_{0}}\sum_{\mathbf{Q}}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}}+\alpha_{2}^{2}\frac{2\pi\lambda_{2}^{4}}{\Omega_{0}}\sum_{\mathbf{Q}}(\mathbf{k}-\mathbf{Q})^{2}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{2}^{2}}\
.$ (53)

In order to correctly capture the strong topology of the system, Ref. [33]
also introduces four anomalous $c$-electrons around the $\Gamma_{M}$ point.
The $c$-fermion operators are defined according to

$\hat{c}^{\dagger}_{\mathbf{k},a,\eta,s}=\sum_{\mathbf{Q},\beta}\tilde{u}^{\eta}_{\mathbf{Q}\beta;a}\left(\mathbf{k}\right)\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s},\mbox{\quad
for\quad}1\leq a\leq 4.$ (54)

In Eq. 54, $a$ indexes the four $c$-electron states (for each spin and
valley), while $\tilde{u}^{\eta}_{\mathbf{Q}\beta;a}\left(\mathbf{k}\right)$
denotes their wave function. The $\hat{c}^{\dagger}_{\mathbf{k},a,\eta,s}$
fermions have a large kinetic energy away from the $\Gamma_{M}$ point, and as
a result, they are only defined for momenta below a certain cutoff
$\absolutevalue{\mathbf{k}}<\Lambda_{c}$, and their wave function in the
proximity of the $\Gamma_{M}$ point can be approximated by [33]

$\tilde{u}^{\eta}_{\mathbf{Q}\beta;a}\left(\mathbf{k}\right)\approx\tilde{u}^{\eta}_{\mathbf{Q}\beta;a}\left(\mathbf{0}\right).$
(55)

Analytical expressions for $c$-electron wave functions were found in Ref.
[35]. These were not directly derived from the BM model approximations [149,
35]. The approximations for the wave functions at $\mathbf{k}=\mathbf{0}$ are
given by [35]

$\displaystyle\tilde{u}_{\mathbf{Q}1,1}^{(\eta)}(\mathbf{0})$
$\displaystyle=-\alpha_{c1}\lambda_{c1}\sqrt{\frac{2\pi}{\Omega_{0}\mathcal{N}_{c1}}}e^{-i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c1}^{2}},\qquad$
$\displaystyle\tilde{u}_{\mathbf{Q}2,1}^{(\eta)}(\mathbf{0})$
$\displaystyle=\alpha_{c2}\lambda_{c2}^{3}\sqrt{\frac{\pi}{\Omega_{0}\mathcal{N}_{c1}}}(i\eta
Q_{x}+Q_{y})^{2}e^{-i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c2}^{2}},$
(56) $\displaystyle\tilde{u}_{\mathbf{Q}1,2}^{(\eta)}(\mathbf{0})$
$\displaystyle=\alpha_{c2}\lambda_{c2}^{3}\sqrt{\frac{\pi}{\Omega_{0}\mathcal{N}_{c2}}}(-i\eta
Q_{x}+Q_{y})^{2}e^{i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c2}^{2}},\qquad$
$\displaystyle\tilde{u}_{\mathbf{Q}2,2}^{(\eta)}(\mathbf{0})$
$\displaystyle=-\alpha_{c1}\lambda_{c2}\sqrt{\frac{2\pi}{\Omega_{0}\mathcal{N}_{c2}}}e^{i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c1}^{2}},$
(57) $\displaystyle\tilde{u}_{\mathbf{Q}1,3}^{(\eta)}(\mathbf{0})$
$\displaystyle=\alpha_{c3}\lambda_{c3}^{2}\sqrt{\frac{2\pi}{\Omega_{0}\mathcal{N}_{c3}}}\zeta_{\mathbf{Q}}(-i\eta
Q_{x}+Q_{y})e^{-i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c3}^{2}},\qquad$
$\displaystyle\tilde{u}_{\mathbf{Q}2,3}^{(\eta)}(\mathbf{0})$
$\displaystyle=\alpha_{c4}\lambda_{c4}^{3}\sqrt{\frac{\pi}{\Omega_{0}\mathcal{N}_{c3}}}(-i\eta
Q_{x}+Q_{y})^{2}e^{-i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c4}^{2}},$
(58) $\displaystyle\tilde{u}_{\mathbf{Q}1,4}^{(\eta)}(\mathbf{0})$
$\displaystyle=\alpha_{c4}\lambda_{c4}^{3}\sqrt{\frac{\pi}{\Omega_{0}\mathcal{N}_{c4}}}(i\eta
Q_{x}+Q_{y})^{2}e^{i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c4}^{2}},\qquad$
$\displaystyle\tilde{u}_{\mathbf{Q}2,4}^{(\eta)}(\mathbf{0})$
$\displaystyle=\alpha_{c3}\lambda_{c3}^{2}\sqrt{\frac{2\pi}{\Omega_{0}\mathcal{N}_{c4}}}\zeta_{\mathbf{Q}}(i\eta
Q_{x}+Q_{y})e^{i\frac{\pi}{4}\zeta_{\mathbf{Q}}-\frac{1}{2}\mathbf{Q}^{2}\lambda_{c3}^{2}},$
(59)

where we have used the sublattice factor $\zeta_{\mathbf{Q}}=+1$
($\zeta_{\mathbf{Q}}=-1$) for $\mathbf{Q}\in\mathcal{Q}_{+}$
($\mathbf{Q}\in\mathcal{Q}_{+}$). These expressions are also extended for
small momenta around the $\Gamma_{M}$ point using Eq. 55 [33]. For
$\theta=$1.05\text{\,}\mathrm{\SIUnitSymbolDegree}$$ and $w_{0}/w_{1}=0.8$,
the parameters appearing in Eqs. 56, 57, 58 and 59 have been fitted to the
numerical solution and are given by [35]

$\displaystyle\alpha_{c1}=0.3958,\qquad$
$\displaystyle\alpha_{c2}=0.9183,\qquad$
$\displaystyle\mathcal{N}_{c1}=\mathcal{N}_{c2}=1.2905,\qquad$
$\displaystyle\lambda_{c1}=0.2194\absolutevalue{\mathbf{a}_{M1}},\qquad$
$\displaystyle\lambda_{c2}=0.3299\absolutevalue{\mathbf{a}_{M1}},\qquad$ (60)
$\displaystyle\alpha_{c3}=0.9257,\qquad$
$\displaystyle\alpha_{c4}=0.3783,\qquad$
$\displaystyle\mathcal{N}_{c3}=\mathcal{N}_{c4}=1.1102,\qquad$
$\displaystyle\lambda_{c3}=0.2430\absolutevalue{\mathbf{a}_{M1}},\qquad$
$\displaystyle\lambda_{c4}=0.2241\absolutevalue{\mathbf{a}_{M1}}.$ (61)

####  A.2.2 The single-particle THF model

The single-particle THF model Hamiltonian for TBG is obtained by projecting
the low-energy fermions from Eq. 36 into the THF model basis via the
approximation

$\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s}\approx\begin{cases}\sum_{\alpha}v^{*\eta}_{\mathbf{Q}\beta;\alpha}\left(\mathbf{k}\right)\hat{f}^{\dagger}_{\mathbf{k},\alpha,\eta,s}+\sum_{a}\tilde{u}^{*\eta}_{\mathbf{Q}\beta;a}\left(\mathbf{k}\right)\hat{c}^{\dagger}_{\mathbf{k},a,\eta,s}&\mbox{\quad
for\quad}\absolutevalue{\mathbf{k}}\leq\Lambda_{c},\\\
\sum_{\alpha}v^{*\eta}_{\mathbf{Q}\beta;\alpha}\left(\mathbf{k}\right)\hat{f}^{\dagger}_{\mathbf{k},\alpha,\eta,s}&\mbox{\quad
for\quad}\absolutevalue{\mathbf{k}}>\Lambda_{c}.\end{cases}$ (62)

and reads as

$\displaystyle H^{\text{TBG}}_{0}=$
$\displaystyle\sum_{\begin{subarray}{c}\absolutevalue{\mathbf{k}}\leq\Lambda_{c}\\\
\eta,s\end{subarray}}\left[\sum_{a,a^{\prime}}h^{cc,\eta}_{aa^{\prime}}\left(\mathbf{k}\right)\hat{c}^{\dagger}_{\mathbf{k},a,\eta,s}\hat{c}_{\mathbf{k},a,\eta,s}+\left(\sum_{a,\alpha}h^{cf,\eta}_{a\alpha}\left(\mathbf{k}\right)\hat{c}^{\dagger}_{\mathbf{k},a,\eta,s}\hat{f}_{\mathbf{k},\alpha,\eta,s}+\text{h.c.}\right)\right]$
$\displaystyle+$
$\displaystyle\sum_{\begin{subarray}{c}\mathbf{k},\alpha,\alpha^{\prime}\\\
\eta,s\end{subarray}}h^{ff,\eta}_{\alpha\alpha^{\prime}}\left(\mathbf{k}\right)\hat{f}^{\dagger}_{\mathbf{k},\alpha,\eta,s}\hat{f}_{\mathbf{k},\alpha^{\prime},\eta,s}.$
(63)

In Eq. 63, “$+\text{h.c.}$” denotes the addition of the Hermitian conjugate.
In order to differentiate the BM and THF model Hamiltonians, the former are
always denoted with a hat. The matrix blocks appearing in Eq. 63 are given by
[33]

$\displaystyle h^{cc,\eta}\left(\mathbf{k}\right)=$
$\displaystyle\begin{pmatrix}\mathbb{0}&v_{\star}\left(\eta
k_{x}\sigma_{0}+ik_{y}\sigma_{z}\right)\\\ v_{*}\left(\eta
k_{x}\sigma_{0}-ik_{y}\sigma_{z}\right)&M\sigma_{z}\end{pmatrix},$ (64)
$\displaystyle h^{cf,\eta}\left(\mathbf{k}\right)=$
$\displaystyle\begin{pmatrix}\gamma\sigma_{0}+v_{*}^{\prime}\left(\eta
k_{x}\sigma_{x}+k_{y}\sigma_{y}\right)\\\ v_{*}^{\prime\prime}\left(\eta
k_{x}\sigma_{x}-k_{y}\sigma_{y}\right)\end{pmatrix}e^{-\frac{\absolutevalue{\mathbf{k}}^{2}\lambda^{2}}{2}},$
(65) $\displaystyle h^{ff,\eta}\left(\mathbf{k}\right)=$
$\displaystyle\mathbb{0},$ (66)

with $\mathbb{0}$ denoting the zero $2\times 2$ (per spin per valley) matrix.
The single-particle parameters appearing in Eqs. 64 and 65 depend on the
parameters of the BM model of TBG and have been obtained both numerically and
analytically across a vast range of twist angles $\theta$ and tunneling ratios
$w_{0}/w_{1}$ [33, 35]. For
$\theta=$1.05\text{\,}\mathrm{\SIUnitSymbolDegree}$$, $w_{0}/w_{1}=0.8$, and
the BM model parameters chosen around Eq. 38,
$\lambda=0.3375\absolutevalue{\mathbf{a}_{M1}}$,
$\gamma=$-24.75\text{\,}\mathrm{meV}$$, $M=$3.697\text{\,}\mathrm{meV}$$,
$v_{\star}=$-4.303\text{\,}\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$,
$v^{\prime}_{\star}=$1.623\text{\,}\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$,
and
$v^{\prime\prime}_{\star}=$-0.0332\text{\,}\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$.

####  A.2.3 The interaction THF Hamiltonian

The THF interaction Hamiltonian was derived by directly projecting density
operator from Eq. 45 with the aid of Eq. 62. Ref. [33] has shown that the
resulting expression can be simplified to a sum of just seven terms which
depend on only six parameters and the Coulomb interaction potential from Eq.
43

$H^{\text{TBG}}_{I}=H_{U_{1}}+H_{W}+H_{V}+H_{J}+H_{U_{2}}+H_{\tilde{J}}+H_{K}.$
(67)

The seven terms of Eq. 67 are given by

$\displaystyle H_{U_{1}}$
$\displaystyle=\frac{U_{1}}{2}\sum_{\mathbf{R}}\sum_{\begin{subarray}{c}\alpha,\eta,s\\\
\alpha^{\prime},\eta^{\prime},s^{\prime}\end{subarray}}:\mathrel{\hat{f}^{\dagger}_{\mathbf{R},\alpha^{\prime},\eta^{\prime},s^{\prime}}\hat{f}_{\mathbf{R},\alpha^{\prime},\eta^{\prime},s^{\prime}}}::\mathrel{\hat{f}^{\dagger}_{\mathbf{R},\alpha,\eta,s}\hat{f}_{\mathbf{R},\alpha,\eta,s}}:,$
(68) $\displaystyle H_{W}$
$\displaystyle=\frac{1}{N_{0}}\sum_{\begin{subarray}{c}\mathbf{k}_{1}\\\
\absolutevalue{\mathbf{k}_{2}}\leq\Lambda_{c}\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{q}\\\
\absolutevalue{\mathbf{k}_{2}-\mathbf{q}}\leq\Lambda_{c}\end{subarray}}\sum_{\begin{subarray}{c}\alpha,\eta,s\\\
a^{\prime},\eta^{\prime},s^{\prime}\end{subarray}}W_{a^{\prime}}:\mathrel{\hat{f}^{\dagger}_{\mathbf{k}_{1}+\mathbf{q},\alpha,\eta,s}\hat{f}_{\mathbf{k}_{1},\alpha,\eta,s}}::\mathrel{\hat{c}^{\dagger}_{\mathbf{k}_{2}-\mathbf{q},a^{\prime},\eta^{\prime},s^{\prime}}\hat{c}_{\mathbf{k}_{2},a^{\prime},\eta^{\prime},s^{\prime}}}:,$
(69) $\displaystyle H_{V}$
$\displaystyle=\frac{1}{2\Omega_{0}N_{0}}\sum_{\absolutevalue{\mathbf{k}_{1}},\absolutevalue{\mathbf{k}_{1}}\leq\Lambda_{c}}\sum_{\begin{subarray}{c}\mathbf{q}\\\
\absolutevalue{\mathbf{k}_{1}+\mathbf{q}},\absolutevalue{\mathbf{k}_{2}+\mathbf{q}}\leq\Lambda_{c}\end{subarray}}\sum_{\begin{subarray}{c}a,\eta,s\\\
a^{\prime},\eta^{\prime},s^{\prime}\end{subarray}}V\left(\mathbf{q}\right):\mathrel{\hat{c}^{\dagger}_{\mathbf{k}_{1}+\mathbf{q},a^{\prime},\eta^{\prime},s^{\prime}}\hat{c}_{\mathbf{k}_{1},a^{\prime},\eta^{\prime},s^{\prime}}}::\mathrel{\hat{c}^{\dagger}_{\mathbf{k}_{2}-\mathbf{q},a,\eta,s}\hat{c}_{\mathbf{k}_{2},a,\eta,s}}:,$
(70) $\displaystyle H_{J}$
$\displaystyle=-\frac{J}{2N_{0}}\sum_{\begin{subarray}{c}\mathbf{k}_{1}\\\
\absolutevalue{\mathbf{k}_{2}}\leq\Lambda_{c}\end{subarray}}\sum_{\begin{subarray}{c}\mathbf{q}\\\
\absolutevalue{\mathbf{k}_{2}+\mathbf{q}}\leq\Lambda_{c}\end{subarray}}\sum_{\begin{subarray}{c}\alpha,\eta,s\\\
a^{\prime},\eta^{\prime},s^{\prime}\end{subarray}}\left[\eta\eta^{\prime}+(-1)^{\alpha+\alpha^{\prime}}\right]:\mathrel{\hat{f}^{\dagger}_{\mathbf{k}_{1}+\mathbf{q},\alpha^{\prime},\eta^{\prime},s^{\prime}}\hat{f}_{\mathbf{k}_{1},\alpha,\eta,s}}::\mathrel{\hat{c}^{\dagger}_{\mathbf{k}_{2}-\mathbf{q},\alpha+2,\eta,s}\hat{c}_{\mathbf{k}_{2},\alpha^{\prime}+2,\eta^{\prime},s^{\prime}}}:,$
(71) $\displaystyle H_{U_{2}}$
$\displaystyle=\frac{U_{2}}{2}\sum_{\left\langle\mathbf{R},\mathbf{R}^{\prime}\right\rangle}\sum_{\begin{subarray}{c}\alpha,\eta,s\\\
\alpha^{\prime},\eta^{\prime},s^{\prime}\end{subarray}}:\mathrel{\hat{f}^{\dagger}_{\mathbf{R}^{\prime},\alpha^{\prime},\eta^{\prime},s^{\prime}}\hat{f}_{\mathbf{R}^{\prime},\alpha^{\prime},\eta^{\prime},s^{\prime}}}::\mathrel{\hat{f}^{\dagger}_{\mathbf{R},\alpha,\eta,s}\hat{f}_{\mathbf{R},\alpha,\eta,s}}:,$
(72) $\displaystyle H_{\tilde{J}}$
$\displaystyle=-\frac{J}{4N_{0}}\sum_{\absolutevalue{\mathbf{k}_{1}},\absolutevalue{\mathbf{k}_{2}}\leq\Lambda_{c}}\sum_{\mathbf{q}}\sum_{\begin{subarray}{c}\alpha,\eta,s\\\
a^{\prime},\eta^{\prime},s^{\prime}\end{subarray}}\left[\eta\eta^{\prime}-(-1)^{\alpha+\alpha^{\prime}}\right]\hat{f}^{\dagger}_{\mathbf{k}_{1}+\mathbf{q},\alpha^{\prime},\eta^{\prime},s^{\prime}}\hat{f}^{\dagger}_{\mathbf{k}_{2}-\mathbf{q},\alpha,\eta,s}\hat{c}_{\mathbf{k}_{2},\alpha+2,\eta,s}\hat{c}_{\mathbf{k}_{1},\alpha^{\prime}+2,\eta^{\prime},s^{\prime}}+\text{h.c.}$
(73) $\displaystyle H_{K}$
$\displaystyle=\frac{K}{N_{0}\Omega_{0}}\sum_{\absolutevalue{\mathbf{k}_{1}},\absolutevalue{\mathbf{k}_{2}},\absolutevalue{\mathbf{k}_{3}}\leq\Lambda_{c}}\sum_{\begin{subarray}{c}\alpha,\eta,s\\\
\eta^{\prime},s^{\prime}\end{subarray}}\eta\eta^{\prime}\left(\hat{c}^{\dagger}_{\mathbf{k}_{1},\bar{\alpha},\eta,s}\hat{c}_{\mathbf{k}_{3},\alpha+2,\eta,s}\hat{f}^{\dagger}_{\mathbf{k}_{1}-\mathbf{k}_{2}-\mathbf{k}_{3},\alpha,\eta^{\prime},s^{\prime}}\hat{c}_{\mathbf{k}_{2},\alpha+2,\eta^{\prime},s^{\prime}}\right.$
$\displaystyle\left.-\hat{f}^{\dagger}_{\mathbf{k}_{2}+\mathbf{k}_{3}-\mathbf{k}_{1},\alpha,\eta^{\prime},s^{\prime}}\hat{c}_{\mathbf{k}_{2},\alpha+2,\eta^{\prime},s^{\prime}}\hat{c}^{\dagger}_{\mathbf{k}_{1},\bar{\alpha}+2,\eta,s}\hat{c}_{\mathbf{k}_{3},\alpha,\eta,s}\right)+\text{h.c.}.$
(74)

where $\bar{\alpha}=3-\alpha$. These terms the represent various interactions
within the THF model: onsite and nearest-neighbor repulsions between
$f$-electrons ($H_{U_{1}}$ and $H_{U_{2}}$), Coulomb interactions among
$c$-electrons ($H_{V}$), density-density and exchange interactions between $f$
and $c$ electrons ($H_{W}$ and $H_{J}$), and interactions involving
hybridization ($H_{\tilde{J}}$ and $H_{K}$). $H_{K}$ is usually neglected as
it does not conserve the $f$ number, which bears a large energy cost when
$U>|\gamma|$. When $U<|\gamma|$, the formalism changes, and the Coulomb
interaction must be projected into the flat bands, giving rise to the
projected strong interaction picture [269, 24, 281, 25, 282, 23]. The normal
ordering of an operator $\mathcal{O}$ subtracts its expectation value in the
charge neutral ground state, $\ket{\mathrm{G}_{0}}$, which is defined such
that

$\displaystyle\expectationvalue{\hat{f}^{\dagger}_{\mathbf{k},\alpha,\eta,s}\hat{f}_{\mathbf{k}^{\prime},\alpha^{\prime},\eta^{\prime},s^{\prime}}}{\mathrm{G}_{0}}$
$\displaystyle=\frac{1}{2}\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{\alpha\alpha^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}},$
$\displaystyle\expectationvalue{\hat{c}^{\dagger}_{\mathbf{k},a,\eta,s}\hat{c}_{\mathbf{k}^{\prime},a^{\prime},\eta^{\prime},s^{\prime}}}{\mathrm{G}_{0}}$
$\displaystyle=\frac{1}{2}\delta_{\mathbf{k},\mathbf{k}^{\prime}}\delta_{aa^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}},$
$\displaystyle\expectationvalue{\hat{f}^{\dagger}_{\mathbf{k},\alpha,\eta,s}\hat{c}_{\mathbf{k}^{\prime},a^{\prime},\eta^{\prime},s^{\prime}}}{\mathrm{G}_{0}}$
$\displaystyle=\expectationvalue{\hat{c}^{\dagger}_{\mathbf{k},a,\eta,s}\hat{f}_{\mathbf{k}^{\prime},\alpha^{\prime},\eta^{\prime},s^{\prime}}}{\mathrm{G}_{0}}=0.$
(75)

Additionally, in Eq. 72,
$\left\langle\mathbf{R},\mathbf{R}^{\prime}\right\rangle$ denotes nearest-
neighbor lattice sites. For
$\theta=$1.05\text{\,}\mathrm{\SIUnitSymbolDegree}$$, $w_{0}/w_{1}=0.8$, as
well as the relative permittivity and screening length chosen around Eq. 43,
the six interaction parameters ($W_{1}=W_{2}$ and $W_{3}=W_{4}$ follows by
symmetry [33]) are given by [33]

$\displaystyle U_{1}$ $\displaystyle=$57.95\text{\,}\mathrm{meV}$,$
$\displaystyle J$ $\displaystyle=$16.38\text{\,}\mathrm{meV}$,$ $\displaystyle
W_{1}$ $\displaystyle=W_{2}=$44.03\text{\,}\mathrm{meV}$,$ $\displaystyle
U_{2}$ $\displaystyle=$2.239\text{\,}\mathrm{meV}$,$ $\displaystyle K$
$\displaystyle=$4.887\text{\,}\mathrm{meV}$,$ $\displaystyle W_{3}$
$\displaystyle=W_{4}=$50.2\text{\,}\mathrm{meV}$,$ (76)

while the $c$-$c$ interaction potential is the same as in Eq. 43.

####  A.2.4 Berry Curvature of the TBG active bands

For each spin $s$ and valley $\eta$, the active TBG bands near charge
neutrality can be recombined into Chern-$(\pm 1)$ bands [279, 269, 58], whose
operators are denoted by $\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}$.
Specifically, $\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}$ creates a state of
spin $s$ in valley $\eta$ within the Chern band with Chern number
$\zeta\cdot\eta$. In terms of the TBG low-energy operators from Eq. 36, the
Chern-band operators are defined according to

$\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}=\sum_{\mathbf{Q},\beta}U^{\zeta}_{\mathbf{Q}\beta;\eta}(\mathbf{k})\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s},$
(77)

with $U^{\zeta}_{\mathbf{Q}\beta;\eta}(\mathbf{k})$ being the corresponding
Chern-band wave function.

The Berry curvature of the TBG active Chern bands can be defined in terms of
the corresponding projector. Letting

$\left[P^{\zeta,\eta}(\mathbf{k})\right]_{\alpha\mathbf{Q};\beta\mathbf{Q}^{\prime}}=U^{\zeta}_{\mathbf{Q}\alpha;\eta}(\mathbf{k})U^{*\zeta}_{\mathbf{Q}^{\prime}\beta;\eta}(\mathbf{k})$
(78)

be the projector into the Chern band corresponding to
$\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}$, the associated Berry curvature
is given by

$\Omega^{\zeta,\eta}(\mathbf{k})=-i\sum_{\mu,\nu}\epsilon^{\mu\nu}\text{Tr}\left[P^{\zeta,\eta}(\mathbf{k})\partial_{k^{\mu}}P^{\zeta,\eta}(\mathbf{k})\partial_{k^{\nu}}P^{\zeta,\eta}(\mathbf{k})\right],$
(79)

where $\epsilon^{\mu\nu}$ is the two-dimensional Levi-Civita symbol with
$\epsilon^{xx}=\epsilon^{yy}=0$ and $\epsilon^{xy}=-\epsilon^{yx}=0$.

In Fig. 7. we plot the distribution of the Berry curvature for the TBG Chern
bands at the magic angle, across a range of tunneling ratios $0\leq
w_{0}/w_{1}\leq 0.8$. For the realistic value $w_{0}/w_{1}=0.8$, the Berry
curvature is concentrated near the $\Gamma_{M}$ point. Conversely, as we
approach the chiral limit ($w_{0}/w_{1}=0$) [61], the Berry curvature becomes
more evenly distributed across the moiré BZ. This uniform spread of the Berry
curvature can be directly linked to the largely increased $f$-$c$
hybridization in the chiral limit [33, 35].

Figure 7: Berry curvature $\Omega^{+,+}(\mathbf{k})$ of the Chern-$(+1)$ band
of TBG within valley $\eta=+$ at the magic angle
$\theta=$1.05\text{\,}\mathrm{\SIUnitSymbolDegree}$$. The BM model parameters
we employ are listed below Eq. 41. (a)-(e) illustrate the effects of varying
the tunneling amplitude ratio (indicated below each panel) from the
unrealistic chiral limit $w_{0}/w_{1}=0$ [61] to the more realistic value
$w_{0}/w_{1}=0.8$. The Berry curvature is plotted in arbitrary units, as
indicated by the accompanying colorbar.

##  Appendix B Analytical results within the THF model

In Ref. [35], it was shown that around the magic angle, the flat and closest
remote bands of TBG within the single-particle BM model are extremely well-
fitted by the THF model for a wide range of angles and tunneling amplitude
ratios. It was, however, argued that the direct utility of the THF model is
dependent on a separation of scales between the $f$-$c$ hybridization
($\gamma$) and the onsite Hubbard repulsion of the $f$-electrons ($U_{1}$)
[35], where the magnitude of the latter should be larger than the former. In
this case, a _real-space_ description of TBG is possible, in which the
$f$-electrons form local moments [33, 273, 274, 45, 276, 277, 50, 275, 49,
278, 171].

Conversely, when $\absolutevalue{\gamma}>U_{1}$, a real-space description of
the interaction terms is no longer useful, as the large $f$-$c$ hybridization
does not allow the polarization of the $f$-fermions to minimize the onsite
Hubbard repulsion $H_{U_{1}}$. Because in the $\absolutevalue{\gamma}>U_{1}$
limit, the energy scale of the interaction is smaller than the gap between the
active and remote TBG, a flat-band-projected approach constitutes a more
realistic formulation of the problem [282, 23]. Such a description involves
the so-called TBG form factors – overlaps between the active TBG band wave
functions at different momenta [269]. The TBG form factors can easily be
computed numerically, but in order to obtain them analytically, the
expressions of the active TBG wave functions are necessary. Obtaining the
latter is an impossible task in the original BM model with a large number of
bands, but becomes tractable within the THF formalism.

In this appendix, we show that in the $\absolutevalue{\gamma}>U_{1}$ limit,
the single-particle THF model can still be employed to derive analytical
expressions for the TBG flat-band wave functions. We then show how the latter
can be used to compute the TBG form factors analytically. In turn, this
enables us to analytically obtain the charge-one excitations of TBG using the
method introduced by Ref. [24]. Finally, we also show that for
$\absolutevalue{\gamma}>U_{1}$, a flat-band projected limit of the THF model
can be defined, wherein the charge-one excitations of TBG around charge
neutrality can also be obtained analytically. This confirms the conceptual
utility of the THF model not just in the $\absolutevalue{\gamma}<U_{1}$ limit,
but also in the opposite limit $\absolutevalue{\gamma}>U_{1}$.

###  B.1 Flat-band wave functions

In the THF model, the wave functions of the flat bands can be obtained
analytically within the zero-bandwidth limit, which is defined when
$v_{\star}^{\prime}=v_{\star}^{\prime\prime}=0$, and $M=0$ [33]. In this
limit, the model exhibits an additional symmetry, which Ref. [33] dubs “chiral
symmetry”. Notably, the chiral symmetry of the THF model within the zero-
bandwidth limit is _different_ from the chiral sublattice symmetry of the TBG
BM model in the $w_{0}/w_{1}=0$ limit [61]. This can be seen by noting that in
the $w_{0}/w_{1}=0$ limit, the active TBG Chern bands are eigenstates of the
chiral sublattice operator and are, therefore, graphene sublattice-polarized.
However, the active TBG Chern bands, obtained analytically from the THF model
in its zero-bandwidth limit [33] are not polarized on the single-layer
graphene sublattices. As a result, the analytical expression for the TBG flat
bands obtained with the THF model in zero-bandwidth limit are still a good
approximation of the realistic TBG bands.

####  B.1.1 Analytical expression for the TBG flat-band wave functions

Working in the zero-bandwidth limit, the single-particle Hamiltonian from Eq.
63 can be diagonalized analytically. The fermion operators
$\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}$ ($\zeta=\pm$) and
$\hat{b}^{\dagger}_{\mathbf{k},j,\eta,s}$ ($j=\pm 1,\pm 2$), corresponding to
the Chern band eigenstates of the active bands introduced in Eq. 77 and to the
remote-band eigenstates, respectively, are given in the THF model basis by

$\displaystyle\hat{d}^{\dagger}_{\mathbf{k},+,\eta,s}$
$\displaystyle=g_{\mathbf{k}}e^{i\eta\theta_{\mathbf{k}}}\hat{f}^{\dagger}_{\mathbf{k},1,\eta,s}-\eta
h_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k},3,\eta,s},$
$\displaystyle\hat{b}^{\dagger}_{\mathbf{k},+1,\eta,s}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(h_{\mathbf{k}}\hat{f}^{\dagger}_{\mathbf{k},1,\eta,s}-\hat{c}^{\dagger}_{\mathbf{k},1,\eta,s}+\eta
g_{\mathbf{k}}e^{-i\eta\theta_{\mathbf{k}}}\hat{c}^{\dagger}_{\mathbf{k},3,\eta,s}\right),$
$\displaystyle\hat{b}^{\dagger}_{\mathbf{k},+2,\eta,s}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(h_{\mathbf{k}}\hat{f}^{\dagger}_{\mathbf{k},1,\eta,s}+\hat{c}^{\dagger}_{\mathbf{k},1,\eta,s}+\eta
g_{\mathbf{k}}e^{-i\eta\theta_{\mathbf{k}}}\hat{c}^{\dagger}_{\mathbf{k},3,\eta,s}\right),$
$\displaystyle\hat{d}^{\dagger}_{\mathbf{k},-,\eta,s}$
$\displaystyle=g_{\mathbf{k}}e^{-i\eta\theta_{\mathbf{k}}}\hat{f}^{\dagger}_{\mathbf{k},2,\eta,s}-\eta
h_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k},4,\eta,s},$
$\displaystyle\hat{b}^{\dagger}_{\mathbf{k},-1,\eta,s}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(h_{\mathbf{k}}\hat{f}^{\dagger}_{\mathbf{k},2,\eta,s}-\hat{c}^{\dagger}_{\mathbf{k},2,\eta,s}+\eta
g_{\mathbf{k}}e^{i\eta\theta_{\mathbf{k}}}\hat{c}^{\dagger}_{\mathbf{k},4,\eta,s}\right),$
$\displaystyle\hat{b}^{\dagger}_{\mathbf{k},-2,\eta,s}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(h_{\mathbf{k}}\hat{f}^{\dagger}_{\mathbf{k},2,\eta,s}+\hat{c}^{\dagger}_{\mathbf{k},2,\eta,s}+\eta
g_{\mathbf{k}}e^{i\eta\theta_{\mathbf{k}}}\hat{c}^{\dagger}_{\mathbf{k},4,\eta,s}\right),$
(80)

where

$\displaystyle
g_{\mathbf{k}}=\frac{v_{\star}\absolutevalue{\mathbf{k}}}{\sqrt{\absolutevalue{v_{\star}\mathbf{k}}^{2}+\gamma^{2}e^{-|\mathbf{k}|^{2}\lambda^{2}}}},\quad
h_{\mathbf{k}}=\frac{\gamma
e^{-|\mathbf{k}|^{2}\lambda^{2}/2}}{\sqrt{\absolutevalue{v_{\star}\mathbf{k}}^{2}+\gamma^{2}e^{-|\mathbf{k}|^{2}\lambda^{2}}}},\mbox{\quad
and\quad}e^{i\theta_{\mathbf{k}}}=\frac{k_{x}+ik_{y}}{\absolutevalue{\mathbf{k}}}.$
(81)

In the resulting band-basis, the single-particle Hamiltonian can then be
written simply as

$H_{0}^{\text{TBG}}=\sum_{\mathbf{k},j,\eta,s}\sqrt{\gamma^{2}e^{-\absolutevalue{\mathbf{k}}^{2}\lambda^{2}}+v_{\star}^{2}\absolutevalue{\mathbf{k}}^{2}}\left(2\absolutevalue{j}-3\right)\hat{b}^{\dagger}_{\mathbf{k},j,\eta,s}\hat{b}_{\mathbf{k},j,\eta,s}.$
(82)

For simplicity, in Eq. 80, we restricted the $c$-electrons’ momenta to the
first moiré BZ. In the generic case that includes all the $c$-electrons with
momenta $\absolutevalue{\mathbf{k}}<\Lambda_{c}$, the
$\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}$ operators are given by

$\displaystyle\hat{d}^{\dagger}_{\mathbf{k},+,\eta,s}$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\hat{f}^{\dagger}_{\mathbf{k},1,\eta,s}-\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\sum_{\begin{subarray}{c}\mathbf{G}\\\
|\mathbf{k}+\mathbf{G}|<\Lambda_{c}\end{subarray}}\frac{\gamma
e^{-\frac{1}{2}\lambda^{2}\absolutevalue{\mathbf{k}+\mathbf{G}}^{2}}/v_{\star}}{\eta(k_{x}+G_{x})+i(k_{y}+G_{y})}\hat{c}^{\dagger}_{\mathbf{k}+\mathbf{G},3,\eta,s},$
(83) $\displaystyle\hat{d}^{\dagger}_{\mathbf{k},-,\eta,s}$
$\displaystyle=\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\hat{f}^{\dagger}_{\mathbf{k},2,\eta,s}-\frac{1}{\sqrt{\mathcal{N}_{\mathbf{k}}}}\sum_{\begin{subarray}{c}\mathbf{G}\\\
|\mathbf{k}+\mathbf{G}|<\Lambda_{c}\end{subarray}}\frac{\gamma
e^{-\frac{1}{2}\lambda^{2}\absolutevalue{\mathbf{k}+\mathbf{G}}^{2}}/v_{\star}}{\eta(k_{x}+G_{x})-i(k_{y}+G_{y})}\hat{c}^{\dagger}_{\mathbf{k}+\mathbf{G},4,\eta,s},$
(84)

where the normalization prefactor reads as

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

We remind the reader that the operator $d_{\mathbf{k},\zeta,\eta,s}^{\dagger}$
corresponds to a Chern band with Chern number $\zeta\cdot\eta$. It is worth
noting that, for the valley $\eta=+$, the
$\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}$ operator is the same as the
Chern-band operator introduced by Ref. [269]. However, in valley $\eta=-$, the
$\hat{d}^{\dagger}_{\mathbf{k},\zeta,-,s}$ operator introduced here
corresponds to the $\hat{d}^{\dagger}_{\mathbf{k},-\zeta,-,s}$ one of Ref.
[269].

The wave function of the TBG flat bands can, therefore, be expressed
_analytically_ in terms of the $f$\- and $c$-electron operators. The latter
also admit analytical expressions given by Eqs. 49, 50, 56, 57, 58 and 59. In
turn, this enables one to obtain the active band TBG wave functions
analytically within the plane wave basis from Eq. 36, a result which is not
otherwise directly available from the BM model.

By combining Eqs. 46, 54 and 80, we can express the electron operators of the
flat TBG Chern bands in terms of the the low-energy TBG operators from Eq. 36
as

$\displaystyle\hat{d}^{\dagger}_{\mathbf{k},+,\eta,s}=g_{\mathbf{k}}e^{i\eta\theta_{\mathbf{k}}}\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;1}(\mathbf{k})\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s}-\eta
h_{\mathbf{k}}\sum_{\mathbf{Q},\beta}\tilde{u}_{\mathbf{Q}\beta;3}^{\eta}(\mathbf{0})\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s},$
$\displaystyle\hat{d}^{\dagger}_{\mathbf{k},-,\eta,s}=g_{\mathbf{k}}e^{-i\eta\theta_{\mathbf{k}}}\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;2}(\mathbf{k})\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s}-\eta
h_{\mathbf{k}}\sum_{\mathbf{Q},\beta}\tilde{u}_{\mathbf{Q}\beta;4}^{\eta}(\mathbf{0})\hat{c}^{\dagger}_{\mathbf{k},\mathbf{Q},\beta,\eta,s}.$
(86)

Comparing with Eq. 77, we can directly identify the Chern band wave functions

$U_{\mathbf{Q}\beta;\eta}^{\zeta}(\mathbf{k})=e^{i\eta\zeta\theta_{\mathbf{k}}}g_{\mathbf{k}}v^{\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k})-\eta
h_{\mathbf{k}}\tilde{u}^{\eta}_{\mathbf{Q}\eta;i_{\zeta}}(\mathbf{0}),$ (87)

where, for simplicity, we have denoted $i_{\zeta=+}=3$, $i_{\zeta=-}=4$.

####  B.1.2 Analytical expressions for the TBG form factors

Another useful quantity which can be obtained analytically using Eq. 87 are
the TBG form factors corresponding to the active bands. The latter arise when
projecting the TBG density operator into the Chern-band basis. Specifically,
the projected TBG density operator is given by

$\delta\rho\left(\mathbf{q}\right)=\sum_{\mathbf{k},\eta,s}\sum_{\zeta,\zeta^{\prime}}M_{\zeta\zeta^{\prime}}^{(\eta)}(\mathbf{k},\mathbf{k}+\mathbf{q})\left(\hat{d}^{\dagger}_{\mathbf{k}+\mathbf{q},\zeta,\eta,s}\hat{d}_{\mathbf{k},\zeta^{\prime},\eta,s}-\frac{1}{2}\delta_{\mathbf{q},\mathbf{0}}\right),$
(88)

where the corresponding form factors are defined by

$M_{\zeta\zeta^{\prime}}^{(\eta)}(\mathbf{k},\mathbf{k}^{\prime})=\sum_{\mathbf{Q},\beta}U^{*\zeta}_{\mathbf{Q}\beta;\eta}(\mathbf{k})U^{\zeta^{\prime}}_{\mathbf{Q}\beta;\eta}(\mathbf{k}^{\prime}).$
(89)

The form factors can be calculated via Eq. 87. Approximately, we find that

$\displaystyle
M_{\zeta\zeta^{\prime}}^{(\eta)}(\mathbf{k},\mathbf{k}^{\prime})\approx\sum_{\mathbf{Q},\beta}\left(g_{\mathbf{k}}g_{\mathbf{k}^{\prime}}e^{-i\theta_{\mathbf{k}}\zeta-i\theta_{\mathbf{k}^{\prime}}\zeta^{\prime}}v^{*\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k})v^{\eta}_{\mathbf{Q}\beta;\zeta^{\prime}}(\mathbf{k}^{\prime})+h_{\mathbf{k}}h_{\mathbf{k}^{\prime}}\tilde{u}^{*\eta}_{\mathbf{Q}\beta;i_{\zeta}}(\mathbf{0})\tilde{u}^{\eta}_{\mathbf{Q}\beta;i_{\zeta^{\prime}}}(\mathbf{0})\right),$
(90)

where we have dropped the cross-terms
$v^{*\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k})\tilde{u}^{\eta}_{\mathbf{Q}\eta;i_{\zeta^{\prime}}}(\mathbf{0})$
and
$v^{\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k}^{\prime})\tilde{u}^{*\eta}_{\mathbf{Q}\eta;i_{\zeta^{\prime}}}(\mathbf{0})$.
In general, the $f$\- and $c$-fermions wave functions are orthogonal at any
momentum $\mathbf{k}$ in the moiré BZ

$\sum_{\mathbf{Q}\beta}v^{*\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k})\tilde{u}^{\eta}_{\mathbf{Q}\beta;i_{\zeta}}(\mathbf{k})=0\,,$
(91)

Coupled with the approximation
$\tilde{u}_{\mathbf{Q}\beta,i}(\mathbf{k})\approx\tilde{u}_{\mathbf{Q}\beta,i}(\mathbf{0})$
from Eq. 55, this implies that the cross-terms should approximately vanish.
Using also the orthonormal condition for the $c$-electron wave functions at
the $\Gamma_{M}$ point,
$\sum_{\mathbf{Q},\beta}\tilde{u}^{*\eta}_{\mathbf{Q}\beta;i_{\zeta}}(\mathbf{0})\tilde{u}^{\eta}_{\mathbf{Q}\beta;i_{\zeta^{\prime}}}(\mathbf{0})=\delta_{\zeta,\zeta^{\prime}}$,
we find that

$M_{\zeta\zeta^{\prime}}^{(\eta)}(\mathbf{k},\mathbf{k}^{\prime})\approx\sum_{\mathbf{Q},\beta}\left(g_{\mathbf{k}}g_{\mathbf{k}^{\prime}}e^{-i(\theta_{\mathbf{k}}-\theta_{\mathbf{k}^{\prime}})\zeta}v^{*\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k})v^{\eta}_{\mathbf{Q}\beta;\zeta^{\prime}}(\mathbf{k}^{\prime})+h_{\mathbf{k}}h_{\mathbf{k}^{\prime}}\delta_{\zeta,\zeta^{\prime}}\right).$
(92)

To further simplify the expression from Eq. 92, we now compute the sum
$\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k})v^{*\eta}_{\mathbf{Q}\beta;\zeta^{\prime}}(\mathbf{k}^{\prime})$
using the analytical formulae for the $f$-electron wave functions from Eqs. 51
and 52. Additionally, we approximate the summation over $\mathbf{Q}$ by an
integral according to

$\sum_{\mathbf{Q}}\rightarrow\frac{\Omega_{0}}{2\pi^{2}}\int\text{d}^{2}Q.$
(93)

This allows us to obtain

$\displaystyle\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;+}(\mathbf{k})v^{*\eta}_{\mathbf{Q}\beta;+}(\mathbf{k}^{\prime})$
$\displaystyle\approx$
$\displaystyle\frac{\Omega_{0}}{2\pi^{2}}\frac{2\pi}{\Omega_{0}\sqrt{\mathcal{N}_{f,\mathbf{k}}\mathcal{N}_{f,\mathbf{k}^{\prime}}}}\int\text{d}^{2}Q\left\\{\lambda_{1}^{2}\alpha_{1}^{2}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}/2-(\mathbf{k}^{\prime}-\mathbf{Q})^{2}\lambda_{1}^{2}/2}\right.$
$\displaystyle\left.+\lambda_{2}^{4}\alpha_{2}^{2}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{2}^{2}/2-(\mathbf{k}^{\prime}-\mathbf{Q})^{2}\lambda_{2}^{2}/2}\left[i\eta(k_{x}-Q_{x})-(k_{y}-Q_{y})\right]\left[-i\eta(k_{x}^{\prime}-Q_{x})-(k_{y}^{\prime}-Q_{y})\right]\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\mathcal{N}_{f,\mathbf{k}}\mathcal{N}_{f,\mathbf{k}^{\prime}}}}\left[\alpha_{1}^{2}e^{-\frac{\absolutevalue{\mathbf{k}_{1}-\mathbf{k}_{2}}^{2}\lambda_{1}^{2}}{4}}+\alpha_{2}^{2}\left(1-\frac{\absolutevalue{\mathbf{k}_{1}-\mathbf{k}_{2}}^{2}\lambda_{1}^{2}}{4}\right)e^{-\frac{\absolutevalue{\mathbf{k}_{1}-\mathbf{k}_{2}}^{2}\lambda_{2}^{2}}{4}}\right],$
(94)
$\displaystyle\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;-}(\mathbf{k})v^{*\eta}_{\mathbf{Q}\beta;-}(\mathbf{k}^{\prime})$
$\displaystyle\approx$
$\displaystyle\frac{\Omega_{0}}{2\pi^{2}}\frac{2\pi}{\Omega_{0}\sqrt{\mathcal{N}_{f,\mathbf{k}}\mathcal{N}_{f,\mathbf{k}^{\prime}}}}\int\text{d}^{2}Q\left\\{\lambda_{1}^{2}\alpha_{1}^{2}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}/2-(\mathbf{k}^{\prime}-\mathbf{Q})^{2}\lambda_{1}^{2}/2}\right.$
$\displaystyle\left.+\lambda_{2}^{4}\alpha_{2}^{2}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}/2-(\mathbf{k}^{\prime}-\mathbf{Q})^{2}\lambda_{1}^{2}/2}\left[-i\eta(k_{x}-Q_{x})-(k_{y}-Q_{y})\right]\left[i\eta(k_{x}^{\prime}-Q_{x})-(k_{y}^{\prime}-Q_{y})\right]\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\mathcal{N}_{f,\mathbf{k}}\mathcal{N}_{f,\mathbf{k}^{\prime}}}}\left[\alpha_{1}^{2}e^{-\frac{\absolutevalue{\mathbf{k}_{1}-\mathbf{k}_{2}}^{2}\lambda_{1}^{2}}{4}}+\alpha_{2}^{2}\left(1-\frac{\absolutevalue{\mathbf{k}_{1}-\mathbf{k}_{2}}^{2}\lambda_{1}^{2}}{4}\right)e^{-\frac{\absolutevalue{\mathbf{k}_{1}-\mathbf{k}_{2}}^{2}\lambda_{2}^{2}}{4}}\right],$
(95)

and

$\displaystyle\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;+}(\mathbf{k})v^{*\eta}_{\mathbf{Q}\beta;-}(\mathbf{k}^{\prime})$
$\displaystyle\approx$
$\displaystyle\frac{\Omega_{0}}{2\pi^{2}}\frac{2\pi}{\Omega_{0}\sqrt{\mathcal{N}_{f,\mathbf{k}}\mathcal{N}_{f,\mathbf{k}^{\prime}}}}\int\text{d}^{2}Q\alpha_{1}\alpha_{2}\lambda_{1}\lambda_{2}^{2}\left[-\eta\left(k_{x}^{\prime}-Q_{x}\right)-i\left(k_{y}^{\prime}-Q_{y}^{\prime}\right)\right]e^{-\frac{1}{2}(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}-\frac{1}{2}(\mathbf{k}^{\prime}-\mathbf{Q})^{2}\lambda_{2}^{2}}$
$\displaystyle+\frac{\Omega_{0}}{2\pi^{2}}\frac{2\pi}{\Omega_{0}\sqrt{\mathcal{N}_{f,\mathbf{k}}\mathcal{N}_{f,\mathbf{k}^{\prime}}}}\int\text{d}^{2}Q\alpha_{1}\alpha_{2}\lambda_{1}\lambda_{2}^{2}\left[-\eta(k_{x}-Q_{x})-i(k_{y}-Q_{y})\right]e^{-\frac{1}{2}(\mathbf{k}-\mathbf{Q})^{2}\lambda_{2}^{2}-\frac{1}{2}(\mathbf{k}^{\prime}-\mathbf{Q})^{2}\lambda_{1}^{2}}=0,$
(96)
$\displaystyle\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;-}(\mathbf{k})v^{*\eta}_{\mathbf{Q}\beta;+}(\mathbf{k}^{\prime})=0.$
(97)

Moreover, the continuum limit from Eq. 93 enables us to approximate the
normalization constant from Eq. 53 and obtain

$\mathcal{N}_{f,\mathbf{k}}\approx\frac{\Omega_{0}}{2\pi^{2}}\int\text{d}^{2}Q\left[\alpha_{1}^{2}\frac{2\pi\lambda_{1}^{2}}{\Omega_{0}}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{1}^{2}}+\alpha_{2}^{2}\frac{2\pi\lambda_{2}^{4}}{\Omega_{0}}(\mathbf{k}-\mathbf{Q})^{2}e^{-(\mathbf{k}-\mathbf{Q})^{2}\lambda_{2}^{2}}\right]=1.$
(98)

As a result, we find that

$\sum_{\mathbf{Q},\beta}v^{\eta}_{\mathbf{Q}\beta;\zeta}(\mathbf{k})v^{*\eta}_{\mathbf{Q}\beta;\zeta^{\prime}}(\mathbf{k}^{\prime})\approx\delta_{\zeta,\zeta^{\prime}}\frac{1}{\alpha_{1}^{2}+\alpha_{2}^{2}}\left[\alpha_{1}^{2}e^{-\frac{\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{1}^{2}}{4}}+\alpha_{2}^{2}\left(1-\frac{\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{1}^{2}}{4}\right)e^{-\frac{\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{2}^{2}}{4}}\right].$
(99)

From Eq. 92 and Eq. 99, the TBG form factors can then be written as

$M_{\zeta\zeta^{\prime}}^{(\eta)}(\mathbf{k},\mathbf{k}^{\prime})=\delta_{\zeta,\zeta^{\prime}}g_{\mathbf{k}}g_{\mathbf{k}^{\prime}}e^{-i(\theta_{\mathbf{k}}-\theta_{\mathbf{k}^{\prime}})\zeta}\left[\alpha_{1}^{2}e^{-\frac{\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{1}^{2}}{4}}+\alpha_{2}^{2}\left(1-\frac{\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{1}^{2}}{4}\right)e^{-\frac{\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{2}^{2}}{4}}\right]+\delta_{\zeta,\zeta^{\prime}}h_{\mathbf{k}}h_{\mathbf{k}^{\prime}}.$
(100)

In the limit $\lambda_{1}\approx\lambda_{2}$ (which holds approximately for
TBG near magic angle), the form factors can be simplified to

$M_{\zeta\zeta^{\prime}}^{(\eta)}(\mathbf{k},\mathbf{k}^{\prime})\approx\delta_{\zeta,\zeta^{\prime}}g_{\mathbf{k}}g_{\mathbf{k}^{\prime}}e^{-i(\theta_{\mathbf{k}}-\theta_{\mathbf{k}^{\prime}})\zeta}\left(1-\frac{\alpha_{2}^{2}\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{1}^{2}}{4}\right)e^{-\frac{\absolutevalue{\mathbf{k}-\mathbf{k}^{\prime}}^{2}\lambda_{1}^{2}}{4}}+\delta_{\zeta,\zeta^{\prime}}h_{\mathbf{k}}h_{\mathbf{k}^{\prime}}\,.$
(101)

###  B.2 The charge-one excitation spectrum of TBG in the
$\absolutevalue{\gamma}>U_{1}$ limit

Armed with analytical expressions of the active TBG bands in both the THF
model basis and the TBG plane wave basis, we now show that the charge-one
excitations of TBG can be obtained analytically in the
$\absolutevalue{\gamma}>U_{1}$ limit. We illustrate two different approaches
to do so: one which leverages the results of Ref. [24] and employs the
analytically derived TBG form factors from Eq. 111, and another one in which
we project the THF interaction Hamiltonian within the active TBG bands using
the analytical eigenstates of the THF single-particle Hamiltonian obtained in
Eq. 80, within the zero-bandwidth limit.

####  B.2.1 Charge-one excitation matrices

As shown in Ref. [24], for the integer-filled exact ground states of TBG
derived by Ref. [23], the charge-one excitation spectrum of TBG can be
calculated directly, as we will review below. For a given exact eigenstate
$\ket{\varphi}$ of the projected interaction $H^{\text{TBG}}_{I}$ [23], the
charge-one excitation is obtained by evaluating the following commutators [24]

$\displaystyle\left[H^{\text{TBG}}_{I}-\mu\hat{N},\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}\right]\ket{\varphi}$
$\displaystyle=\sum_{\zeta^{\prime}}R^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)\hat{d}^{\dagger}_{\mathbf{k},\zeta^{\prime},\eta,s}\ket{\varphi},$
(102)
$\displaystyle\left[H^{\text{TBG}}_{I}-\mu\hat{N},\hat{d}_{\mathbf{k},\zeta,\eta,s}\right]\ket{\varphi}$
$\displaystyle=\sum_{\zeta^{\prime}}\tilde{R}^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)\hat{d}_{\mathbf{k},\zeta^{\prime},\eta,s}\ket{\varphi},$
(103)

where $\mu$ is the chemical potential and $\hat{N}$ is defined as

$\hat{N}=\sum_{\mathbf{k},\zeta,\eta,s}\hat{d}^{\dagger}_{\mathbf{k},\zeta,\eta,s}\hat{d}_{\mathbf{k},\zeta,\eta,s}.$
(104)

The matrices $R^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)$ and
$\tilde{R}^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)$ in Eqs. 102
and 103 are determined by the filling $\nu$ of the ground states
$\ket{\varphi}$ and by the form factors of Eq. 89. The explicit formulas are
given below [24]

$\displaystyle R^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)$
$\displaystyle=\frac{1}{2\Omega_{0}N_{0}}\sum_{\mathbf{G}\in\mathcal{Q}_{0}}\left[\left(\sum_{\mathbf{q}\in\text{MBZ},\zeta^{\prime\prime}}V\left(\mathbf{q}+\mathbf{G}\right)M^{*(\eta)}_{\zeta^{\prime\prime}\zeta^{\prime}}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)M^{(\eta)}_{\zeta^{\prime\prime}\zeta}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)\right)\right.$
$\displaystyle\left.+2A_{-\mathbf{G}}\sqrt{V\left(\mathbf{G}\right)}M^{(\eta)}_{\zeta^{\prime}\zeta}\left(\mathbf{k}+\mathbf{G},\mathbf{k}\right)\right]-\mu\delta_{\zeta^{\prime},\zeta},$
(105)
$\displaystyle\tilde{R}^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)$
$\displaystyle=\frac{1}{2\Omega_{0}N_{0}}\sum_{\mathbf{G}\in\mathcal{Q}_{0}}\left[\left(\sum_{\mathbf{q}\in\text{MBZ},\zeta^{\prime\prime}}V\left(\mathbf{q}+\mathbf{G}\right)M^{(\eta)}_{\zeta^{\prime\prime}\zeta^{\prime}}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)M^{*(\eta)}_{\zeta^{\prime\prime}\zeta}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)\right)\right.$
$\displaystyle\left.-2A_{-\mathbf{G}}\sqrt{V\left(\mathbf{G}\right)}M^{*(\eta)}_{\zeta^{\prime}\zeta}\left(\mathbf{k}+\mathbf{G},\mathbf{k}\right)\right]+\mu\delta_{\zeta^{\prime},\zeta},$
(106)

where MBZ denotes the moiré BZ and $A_{\mathbf{G}}$ is defined as [24]

$A_{\mathbf{G}}=\sqrt{V\left(\mathbf{G}\right)}\sum_{\begin{subarray}{c}\mathbf{k}\in\text{MBZ}\\\
\zeta,\eta\end{subarray}}\frac{\nu}{4}M^{(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{G},\mathbf{k}\right)=\sqrt{V\left(\mathbf{G}\right)}\sum_{\mathbf{k}\in\text{MBZ}}\frac{\nu}{4}\text{Tr}\left[M\left(\mathbf{k}+\mathbf{G},\mathbf{k}\right)\right].$
(107)

Using the fact that the approximate TBG form factors from Eq. 111 are diagonal
in the Chern band basis, we find that the corresponding charge-one excitations
are also diagonal in the same basis, _i.e._
$R^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)=\delta_{\zeta^{\prime},\zeta}R^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)$
and
$\tilde{R}^{\eta}_{\zeta^{\prime}\zeta}\left(\mathbf{k}\right)=\delta_{\zeta^{\prime},\zeta}\tilde{R}^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)$
with

$\displaystyle R^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)$
$\displaystyle=\frac{1}{2\Omega_{0}N_{0}}\sum_{\mathbf{G}\in\mathcal{Q}_{0}}\left[\left(\sum_{\mathbf{q}\in\text{MBZ}}V\left(\mathbf{q}+\mathbf{G}\right)M^{*(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)M^{(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)\right)\right.$
$\displaystyle\left.+2A_{-\mathbf{G}}\sqrt{V\left(\mathbf{G}\right)}M^{(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{G},\mathbf{k}\right)\right]-\mu,$
(108) $\displaystyle\tilde{R}^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)$
$\displaystyle=\frac{1}{2\Omega_{0}N_{0}}\sum_{\mathbf{G}\in\mathcal{Q}_{0}}\left[\left(\sum_{\mathbf{q}\in\text{MBZ}}V\left(\mathbf{q}+\mathbf{G}\right)M^{(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)M^{*(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)\right)\right.$
$\displaystyle\left.-2A_{-\mathbf{G}}\sqrt{V\left(\mathbf{G}\right)}M^{*(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{G},\mathbf{k}\right)\right]+\mu.$
(109)

At the same time, it is easy to see that
$M^{*(\eta)}_{(-\zeta)(-\zeta)}\left(\mathbf{k},\mathbf{k}^{\prime}\right)=M^{(\eta)}_{\zeta\zeta}\left(\mathbf{k},\mathbf{k}^{\prime}\right)$,
which implies that $A_{\mathbf{G}}$ is real and that
$R^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)=R^{*\eta}_{(-\zeta)(-\zeta)}\left(\mathbf{k}\right)$
and
$\tilde{R}^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)=\tilde{R}^{*\eta}_{(-\zeta)(-\zeta)}\left(\mathbf{k}\right)$.
However, while the first term in each of Eqs. 108 and 109 is real, the second
term is not manifestly so. Generically, the TBG form factors obey the
hermicitiy property
$M^{(\eta)}_{\zeta\zeta^{\prime}}\left(\mathbf{k},\mathbf{k}^{\prime}\right)=M^{*(\eta)}_{\zeta^{\prime}\zeta}\left(\mathbf{k}^{\prime},\mathbf{k}\right)$,
which follows trivially from their definition in Eq. 89. Together with the
even symmetries of $A_{\mathbf{G}}=A_{-\mathbf{G}}$ and
$V(\mathbf{G})=V(-\mathbf{G})$, the hermiticity property of the form factors
coupled with the periodicity in reciprocal space
$M^{(\eta)}_{\zeta\zeta^{\prime}}\left(\mathbf{k}+\mathbf{G},\mathbf{k}^{\prime}+\mathbf{G}\right)=M^{(\eta)}_{\zeta\zeta^{\prime}}\left(\mathbf{k}+\mathbf{G},\mathbf{k}^{\prime}+\mathbf{G}\right)$
is enough to guarantee the reality of the second term of each of Eqs. 108 and
109. However, the approximation in Eq. 111 does not obey the periodicity
property, which makes the charge-one excitation matrices nonhermitian at
$\nu\neq 0$. For $\nu=0$ (or alternatively when the so-called flat metric
condition is imposed [24] at generic $\nu\neq 0$), $A_{\mathbf{G}}=0$, which
implies that the charge-one excitation matrices are proportional to the
identity matrix (and are therefore bona-fide hermitian matrices).

####  B.2.2 Analytical expressions for the charge-one excitation matrices of
TBG

We now provide asymptotic analytical expression for the charge-one excitation
matrices at $\nu=0$. The latter can be rewritten from Eqs. 108 and 109 as

$\displaystyle
R^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)=\tilde{R}^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)=$
$\displaystyle\frac{1}{2\Omega_{0}N_{0}}\frac{\Omega_{0}N_{0}}{(2\pi)^{2}}\sum_{\mathbf{G}\in\mathcal{Q}_{0}}\int_{\text{MBZ}}V\left(\mathbf{q}+\mathbf{G}\right)\absolutevalue{M^{(\eta)}_{\zeta\zeta}\left(\mathbf{k}+\mathbf{q}+\mathbf{G},\mathbf{k}\right)}^{2}\text{d}^{2}q$
$\displaystyle=$
$\displaystyle\frac{1}{8\pi^{2}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}V\left(\mathbf{q}-\mathbf{k}\right)\absolutevalue{M^{(\eta)}_{\zeta\zeta}\left(\mathbf{q},\mathbf{k}\right)}^{2}\text{d}q_{x}\text{d}q_{y}.$
(110)

To move forward, we need to make two simplifying assumptions. First, we
consider the
$\sqrt{\absolutevalue{v_{\star}\mathbf{k}}^{2}+\gamma^{2}e^{-|\mathbf{k}|^{2}\lambda^{2}}}$
denominator of the $g_{\mathbf{k}}$ and $h_{\mathbf{k}}$ functions from Eq.
81, which appears in the expression of the form factors from Eq. 111. For
small values of $\absolutevalue{\mathbf{k}}$ the exponential
$e^{-\absolutevalue{\mathbf{k}}^{2}\lambda^{2}}$ is approximately 1. At large
$\absolutevalue{\mathbf{k}}$, $g_{\mathbf{k}}$ is already multiplied by
decaying exponentials in Eq. 111, while $h_{\mathbf{k}}$ itself contains a
decaying exponential factor in the numerator. In other words, for large
$\absolutevalue{\mathbf{q}}$, the exact values of the $g_{\mathbf{q}}$ and
$h_{\mathbf{q}}$ functions is unimportant, as they become exponentially
suppressed in the form factors. As a result, we will ignore the exponential
decay factor appearing in the denominators of the $g_{\mathbf{q}}$ and
$h_{\mathbf{q}}$ functions in Eq. 110 and approximate the form factor as

$M_{\zeta\zeta}^{(\eta)}(\mathbf{k},\mathbf{q})=g_{\mathbf{k}}\tilde{g}_{\mathbf{q}}e^{-i(\theta_{\mathbf{k}}-\theta_{\mathbf{q}})\zeta}\left[\alpha_{1}^{2}e^{-\frac{\absolutevalue{\mathbf{k}-\mathbf{q}}^{2}\lambda_{1}^{2}}{4}}+\alpha_{2}^{2}\left(1-\frac{\absolutevalue{\mathbf{k}-\mathbf{q}}^{2}\lambda_{1}^{2}}{4}\right)e^{-\frac{\absolutevalue{\mathbf{k}-\mathbf{q}}^{2}\lambda_{2}^{2}}{4}}\right]+h_{\mathbf{k}}\tilde{h}_{\mathbf{q}},$
(111)

where we have defined

$\tilde{g}_{\mathbf{q}}\approx\frac{v_{\star}\absolutevalue{\mathbf{q}}}{\sqrt{\absolutevalue{v_{\star}\mathbf{q}}^{2}+\gamma^{2}}}\mbox{\quad
and\quad}\tilde{h}_{\mathbf{q}}=\frac{\gamma
e^{-|\mathbf{q}|^{2}\lambda^{2}/2}}{\sqrt{\absolutevalue{v_{\star}\mathbf{q}}^{2}+\gamma^{2}}}.$
(112)

Additionally, we note that, while the form factor decays exponentially in
$\absolutevalue{\mathbf{q}}$, the screened Coulomb potential only decays
polynomially over the same momentum scale. As a simple approximation, we can
can therefore approximate the Coulomb potential to a constant

$V\left(\mathbf{q}-\mathbf{k}\right)\approx U_{\xi}\pi\xi^{2}.$ (113)

Finally, we note that the integrand in Eq. 110 only depends on the the
relative angle between $\mathbf{k}$ and $\mathbf{q}$ and their magnitudes.
Letting $\mathbf{k}=\left(k,0\right)$ and
$\mathbf{q}=\left(q\cos\theta,q\sin\theta\right)$ without loss of generality,
the charge-one excitation matrices become

$\displaystyle R^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)=$
$\displaystyle\frac{U_{\xi}\xi^{2}}{8\pi}\int_{0}^{\infty}q\text{d}q\int_{0}^{2\pi}\text{d}\theta\left(g^{2}_{k}\tilde{g}^{2}_{q}f^{2}\left(k,q,\theta\right)+h^{2}_{k}\tilde{h}^{2}_{q}+2g_{k}h_{k}\tilde{g}_{q}\tilde{h}_{q}f\left(k,q,\theta\right)\cos\theta\right),$
(114)

with

$f(k,q,\theta)=\alpha_{1}^{2}e^{-\frac{\lambda_{1}^{2}\left(k^{2}+q^{2}-2kq\cos\theta\right)}{4}}+\alpha_{2}^{2}\left(1-\frac{\lambda_{1}^{2}\left(k^{2}+q^{2}-2kq\cos\theta\right)}{4}\right)e^{-\frac{\lambda_{2}^{2}\left(k^{2}+q^{2}-2kq\cos\theta\right)}{4}}.$
(115)

The integral over the angle in Eq. 114 can be readily performed giving

$\displaystyle
R^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)=\frac{U_{\xi}\xi^{2}}{8\pi}\frac{1}{k^{2}v_{\star}^{2}+\gamma^{2}e^{-\frac{k^{2}\lambda^{2}}{2}}}\int_{0}^{\infty}\frac{q\text{d}q}{q^{2}v_{\star}^{2}+\gamma^{2}}F(k,q),$
(116)

where we have defined

$\displaystyle F(k,q)=$ $\displaystyle
2\pi\gamma^{4}e^{-\lambda^{2}\left(k^{2}+q^{2}\right)}$ $\displaystyle+$
$\displaystyle
2\pi\alpha_{1}^{4}k^{2}q^{2}v_{*}^{4}e^{-\frac{1}{2}\lambda_{1}^{2}\left(k^{2}+q^{2}\right)}I_{0}\left(k\lambda_{1}^{2}q\right)$
$\displaystyle+$
$\displaystyle\frac{1}{8}\pi\alpha_{2}^{4}k^{2}q^{2}v_{*}^{4}e^{-\frac{1}{2}\lambda_{2}^{2}\left(k^{2}+q^{2}\right)}\left(k^{2}\lambda_{1}^{2}+\lambda_{1}^{2}q^{2}-4\right){}^{2}I_{0}\left(k\lambda_{2}^{2}q\right)$
$\displaystyle+$
$\displaystyle\pi\alpha_{1}^{2}\alpha_{2}^{2}k^{2}q^{2}v_{*}^{4}\left(-e^{-\frac{1}{4}\left(\lambda_{1}^{2}+\lambda_{2}^{2}\right)\left(k^{2}+q^{2}\right)}\right)\left(k^{2}\lambda_{1}^{2}+\lambda_{1}^{2}q^{2}-4\right)I_{0}\left(\frac{1}{2}k\left(\lambda_{1}^{2}+\lambda_{2}^{2}\right)q\right)$
$\displaystyle+$ $\displaystyle
4\pi\alpha_{1}^{2}\gamma^{2}kqv_{*}^{2}e^{-\frac{1}{4}\left(2\lambda^{2}+\lambda_{1}^{2}\right)\left(k^{2}+q^{2}\right)}I_{1}\left(\frac{1}{2}k\lambda_{1}^{2}q\right)$
$\displaystyle-$
$\displaystyle\frac{\pi\alpha_{2}^{2}\gamma^{2}kqv_{*}^{2}e^{-\frac{1}{4}\left(2\lambda^{2}+\lambda_{2}^{2}\right)\left(k^{2}+q^{2}\right)}\left(\lambda_{1}^{2}\left(k^{2}\lambda_{2}^{2}+\lambda_{2}^{2}q^{2}-4\right)-4\lambda_{2}^{2}\right)}{\lambda_{2}^{2}}I_{1}\left(\frac{1}{2}k\lambda_{2}^{2}q\right)$
$\displaystyle-$
$\displaystyle\frac{\pi\alpha_{2}^{4}k^{3}\lambda_{1}^{2}q^{3}v_{*}^{4}e^{-\frac{1}{2}\lambda_{2}^{2}\left(k^{2}+q^{2}\right)}\left(\lambda_{1}^{2}\left(k^{2}\lambda_{2}^{2}+\lambda_{2}^{2}q^{2}-1\right)-4\lambda_{2}^{2}\right)}{2\lambda_{2}^{2}}I_{1}\left(k\lambda_{2}^{2}q\right)$
$\displaystyle+$ $\displaystyle
2\pi\alpha_{1}^{2}\alpha_{2}^{2}k^{3}\lambda_{1}^{2}q^{3}v_{*}^{4}e^{-\frac{1}{4}\left(\lambda_{1}^{2}+\lambda_{2}^{2}\right)\left(k^{2}+q^{2}\right)}I_{1}\left(\frac{1}{2}k\left(\lambda_{1}^{2}+\lambda_{2}^{2}\right)q\right)$
$\displaystyle+$ $\displaystyle
2\pi\alpha_{2}^{2}\gamma^{2}k^{2}\lambda_{1}^{2}q^{2}v_{*}^{2}e^{-\frac{1}{4}\left(2\lambda^{2}+\lambda_{2}^{2}\right)\left(k^{2}+q^{2}\right)}I_{2}\left(\frac{1}{2}k\lambda_{2}^{2}q\right)$
$\displaystyle+$
$\displaystyle\frac{1}{2}\pi\alpha_{2}^{4}k^{4}\lambda_{1}^{4}q^{4}v_{*}^{4}e^{-\frac{1}{2}\lambda_{2}^{2}\left(k^{2}+q^{2}\right)}I_{2}\left(k\lambda_{2}^{2}q\right),$
(117)

where $I_{n}\left(z\right)$ is the $n$-th order modified Bessel function of
the first kind. For the first term of $F(k,q)$, the integral over $q$ can be
directly performed

$\int_{0}^{\infty}2\pi\gamma^{4}e^{-\lambda^{2}\left(k^{2}+q^{2}\right)}\frac{q\text{d}q}{q^{2}v_{\star}^{2}+\gamma^{2}}=\frac{\pi\gamma^{4}e^{-\lambda^{2}\left(k^{2}-\frac{\gamma^{2}}{v_{*}^{2}}\right)}\Gamma\left(0,\frac{\gamma^{2}\lambda^{2}}{v_{*}^{2}}\right)}{v_{*}^{2}}.$
(118)

For the other terms, we apply the following strategy. First, we expand the
Bessel functions according to

$I_{n}(z)=\sum_{m=0}^{\infty}\frac{1}{m!(m+n)!}\left(\frac{z}{2}\right)^{2m+n},\mbox{\quad
for\quad}n\in\mathbb{N},$ (119)

which is equivalent to a power series expansion in $q$ of $F(k,q)$. We find it
easier to perform this step for each term of $F(k,q)$ separately. Consider the
typical term in Eq. 117, which has the form
$\mathcal{C}(k)q^{n_{1}}e^{-\frac{\Lambda_{1}^{2}}{2}q^{2}}I_{n_{2}}\left(\frac{kq\Lambda_{2}^{2}}{2}\right)$,
where $\mathcal{C}(k)$ is a function of $k$, $n_{1},n_{2}\in\mathbb{N}$, while
$\Lambda_{1}$ and $\Lambda_{2}$ are combinations of $\lambda$, $\lambda_{1}$,
or $\lambda_{2}$ having dimension of length. For example, the second-to-last
term of Eq. 117 has this form with

$\displaystyle\mathcal{C}(k)$
$\displaystyle=2\pi\alpha_{2}^{2}\gamma^{2}k^{2}\lambda_{1}^{2}v_{*}^{2}e^{-\frac{1}{4}\left(2\lambda^{2}+\lambda_{2}^{2}\right)k^{2}},$
(120) $\displaystyle n_{1}=2,\quad n_{2}$
$\displaystyle=2,\quad\Lambda_{1}=\sqrt{2\lambda^{2}+\lambda^{2}_{2}},\quad\Lambda_{2}=\lambda^{2}_{2}.$
(121)

Such a term can be integrated by first expanding the Bessel function according
to Eq. 119 and then integrating the resulting series in a term-by-term fashion

$\displaystyle\int_{0}^{\infty}\frac{q\text{d}q}{q^{2}v_{\star}^{2}+\gamma^{2}}\mathcal{C}(k)q^{n_{1}}e^{-\frac{\Lambda_{1}^{2}}{2}q^{2}}I_{n_{2}}\left(\frac{kq\Lambda_{2}^{2}}{2}\right)=$
$\displaystyle=$
$\displaystyle\frac{1}{v^{2}_{\star}}\mathcal{C}(k)\absolutevalue{\frac{\gamma}{v_{\star}}}^{n_{1}}\int_{0}^{\infty}\frac{z\text{d}z}{1+z^{2}}z^{n_{1}}e^{-\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}z^{2}}I_{n_{2}}\left(kz\Lambda_{2}^{2}\absolutevalue{\frac{\gamma}{2v_{\star}}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{v^{2}_{\star}}\mathcal{C}(k)\absolutevalue{\frac{\gamma}{v_{\star}}}^{n_{1}}\sum_{m=0}^{\infty}\frac{1}{m!(m+n_{2})!}\left(k\Lambda_{2}^{2}\absolutevalue{\frac{\gamma}{4v_{\star}}}\right)^{2m+n_{2}}\int_{0}^{\infty}\frac{z^{n_{1}+n_{2}+2m+1}\text{d}z}{1+z^{2}}e^{-\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}z^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{v^{2}_{\star}}\mathcal{C}(k)\absolutevalue{\frac{\gamma}{v_{\star}}}^{n_{1}}\sum_{m=0}^{\infty}\frac{\Gamma\left(\frac{n_{1}+n_{2}}{2}+m+1\right)}{m!(m+n_{2})!}\left(k\Lambda_{2}^{2}\absolutevalue{\frac{\gamma}{4v_{\star}}}\right)^{2m+n_{2}}\frac{1}{2}e^{\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}}\Gamma\left(-\frac{n_{1}+n_{2}}{2}-m,\frac{\Lambda_{1}^{2}\gamma^{2}}{4v^{2}_{\star}}\right).$
(122)

We now note that for all the terms in Eq. 117,
$\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\ll 1$, suggesting that the
incomplete Gamma function can be expanded in the second argument. Moreover, we
find that for all the terms in Eq. 117 proportional to a Bessel function,
$n_{1}+n_{2}\geq 2$ is even. As such, we can employ the following
representation of the incomplete Gamma function

$\Gamma(-m,z)=\frac{(-1)^{m}}{m!}\left[\Gamma(0,z)-e^{-z}\sum_{j=0}^{m-1}\frac{j!(-1)^{j}}{z^{j+1}}\right],\mbox{\quad
for\quad}m\in\mathbb{N}.$ (123)

For small $z$, the incomplete Gamma function has the following asymptotic form

$\Gamma(-m,z)=\frac{1}{mz^{m}}+\begin{cases}\mathcal{O}\left(z^{-m+1}\right),&\quad
m\geq 1\\\ \mathcal{O}\left(\log z\right),&\quad m=1\\\ \end{cases}$ (124)

which enables us to rewrite the integral in Eq. 122 as

$\displaystyle\int_{0}^{\infty}\frac{q\text{d}q}{q^{2}v_{\star}^{2}+\gamma^{2}}\mathcal{C}(k)q^{n_{1}}e^{-\frac{\Lambda_{1}^{2}}{2}q^{2}}I_{n_{2}}\left(\frac{kq\Lambda_{2}^{2}}{2}\right)=$
$\displaystyle=$
$\displaystyle\frac{1}{v^{2}_{\star}}\mathcal{C}(k)\absolutevalue{\frac{\gamma}{v_{\star}}}^{n_{1}}\sum_{m=0}^{\infty}\frac{\Gamma\left(\frac{n_{1}+n_{2}}{2}+m+1\right)}{m!(m+n_{2})!}\left(k\Lambda_{2}^{2}\absolutevalue{\frac{\gamma}{4v_{\star}}}\right)^{2m+n_{2}}$
$\displaystyle\times\left[\frac{\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)^{-\frac{n_{1}+n_{2}}{2}-m}}{n_{1}+n_{2}+2m}+\mathcal{O}\left(\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)^{-\frac{n_{1}+n_{2}}{2}-m+1}+\log\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)\right)\right]$
$\displaystyle=$
$\displaystyle\frac{1}{v^{2}_{\star}}\mathcal{C}(k)\absolutevalue{\frac{\gamma}{v_{\star}}}^{n_{1}}\sum_{m=0}^{\infty}\frac{\Gamma\left(\frac{n_{1}+n_{2}}{2}+m+1\right)}{m!(m+n_{2})!}\left(k\Lambda_{2}^{2}\absolutevalue{\frac{\gamma}{4v_{\star}}}\right)^{n_{2}}\left(\frac{k^{2}\Lambda_{2}^{2}}{2\sqrt{2}}\right)^{m}\left(\frac{\Lambda_{2}}{\Lambda_{1}}\right)^{m}\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)^{m}$
$\displaystyle\times\left[\frac{\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)^{-\frac{n_{1}+n_{2}}{2}-m}}{n_{1}+n_{2}+2m}+\mathcal{O}\left(\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)^{-\frac{n_{1}+n_{2}}{2}-m+1}+\log\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)\right)\right].$
(125)

Keeping only the leading contribution in the small parameter
$\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}$, we can approximate the
integral from Eq. 122 as

$\displaystyle\int_{0}^{\infty}\frac{q\text{d}q}{q^{2}v_{\star}^{2}+\gamma^{2}}\mathcal{C}(k)q^{n_{1}}e^{-\frac{\Lambda_{1}^{2}}{2}q^{2}}I_{n_{2}}\left(\frac{kq\Lambda_{2}^{2}}{2}\right)\approx$
$\displaystyle\approx\frac{1}{v^{2}_{\star}}\mathcal{C}(k)\absolutevalue{\frac{\gamma}{v_{\star}}}^{n_{1}}\sum_{m=0}^{\infty}\frac{\Gamma\left(\frac{n_{1}+n_{2}}{2}+m+1\right)}{m!(m+n_{2})!}\left(k\Lambda_{2}^{2}\absolutevalue{\frac{\gamma}{4v_{\star}}}\right)^{n_{2}}\left(\frac{k^{2}\Lambda_{2}^{2}}{2\sqrt{2}}\right)^{m}\left(\frac{\Lambda_{2}}{\Lambda_{1}}\right)^{m}\frac{\left(\frac{\Lambda_{1}^{2}\gamma^{2}}{2v^{2}_{\star}}\right)^{-\frac{n_{1}+n_{2}}{2}}}{n_{1}+n_{2}+2m}.$
(126)

Finally, the series in Eq. 126 can be re-summed into an elementary function.
Performing this approximation for each term in Eq. 117, we arrive at the
following result

$\displaystyle\frac{8\pi}{U_{\xi}\xi^{2}}\left(k^{2}v_{\star}^{2}+\gamma^{2}e^{-\frac{k^{2}\lambda^{2}}{2}}\right)R^{\eta}_{\zeta\zeta}\left(\mathbf{k}\right)\approx$
$\displaystyle\approx$
$\displaystyle\frac{\pi\gamma^{4}e^{-\lambda^{2}\left(k^{2}-\frac{\gamma^{2}}{v_{*}^{2}}\right)}\Gamma\left(0,\frac{\gamma^{2}\lambda^{2}}{v_{*}^{2}}\right)}{v_{*}^{2}}$
$\displaystyle+$