# MATBG as Topological Heavy Fermion: I. Exact Mapping and Correlated
Insulators

Zhi-Da Song<EMAIL_ADDRESS>International Center for Quantum Materials,
School of Physics, Peking University, Beijing 100871, China  Department of
Physics, Princeton University, Princeton, New Jersey 08544, USA  B. Andrei
Bernevig  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

Magic-angle ($\theta=1.05^{\circ}$) twisted bilayer graphene (MATBG) has shown
two seemingly contradictory characters: the localization and quantum-dot-like
behavior in STM experiments, and delocalization in transport experiments. We
construct a model, which naturally captures the two aspects, from the
Bistritzer-MacDonald (BM) model in a first principle spirit. A set of local
flat-band orbitals ($f$) centered at the AA-stacking regions are responsible
to the localization. A set of extended topological semi-metallic conduction
bands ($c$), which are at small energetic separation from the local orbitals,
are responsible to the delocalization and transport. The topological flat
bands of the BM model appear as a result of the hybridization of $f$\- and
$c$-electrons. This model then provides a new perspective for the strong
correlation physics, which is now described as strongly correlated
$f$-electrons coupled to nearly free $c$-electrons - we hence name our model
as the topological heavy fermion model. Using this model, we obtain the U(4)
and U(4)$\times$U(4) symmetries of Refs. Kang and Vafek (2019); Bultinck  _et
al._ (2020a); Vafek and Kang (2020); Bernevig  _et al._ (2021a); Lian  _et
al._ (2021) as well as the correlated insulator phases and their energies.
Simple rules for the ground states and their Chern numbers are derived.
Moreover, features such as the large dispersion of the charge $\pm 1$
excitations Bernevig  _et al._ (2021b); Vafek and Kang (2021); Bultinck  _et
al._ (2020a), and the minima of the charge gap at the $\Gamma_{M}$ point can
now, for the first time, be understood both qualitatively and quantitatively
in a simple physical picture. Our mapping opens the prospect of using heavy-
fermion physics machinery to the superconducting physics of MATBG.

Introduction — Since the initial experimental discovery of the correlated
insulator phases Cao  _et al._ (2018a) and superconductivity Cao  _et al._
(2018b) in MATBG Bistritzer and MacDonald (2011), extensive experimental Lu
_et al._ (2019); Yankowitz  _et al._ (2019); Sharpe  _et al._ (2019); Saito
_et al._ (2020); Stepanov  _et al._ (2020); Liu  _et al._ (2021a); Arora  _et
al._ (2020); Serlin  _et al._ (2019); Cao  _et al._ (2020); Polshyn  _et al._
(2019); Xie  _et al._ (2019); Choi  _et al._ (2019); Kerelsky  _et al._
(2019); Jiang  _et al._ (2019); Wong  _et al._ (2020); Zondiner  _et al._
(2020); Nuckolls  _et al._ (2020); Choi  _et al._ (2021); Saito  _et al._
(2021a); Das  _et al._ (2021); Wu  _et al._ (2021); Park  _et al._ (2021);
Saito  _et al._ (2021b); Rozen  _et al._ (2021); Lu  _et al._ (2021) and
theoretical Tarnopolsky  _et al._ (2019); Zou  _et al._ (2018); Liu  _et al._
(2019a); Efimkin and MacDonald (2018); Kang and Vafek (2018); Song  _et al._
(2019); Po  _et al._ (2019); Hejazi  _et al._ (2019a); Padhi  _et al._ (2018);
Lian  _et al._ (2020); Hejazi  _et al._ (2019b); Padhi  _et al._ (2020); Xu
and Balents (2018); Koshino  _et al._ (2018); Ochi  _et al._ (2018); Xu  _et
al._ (2018); Guinea and Walet (2018); Venderbos and Fernandes (2018); You and
Vishwanath (2019); Wu and Das Sarma (2020); Lian  _et al._ (2019); Wu  _et
al._ (2018); Isobe  _et al._ (2018); Liu  _et al._ (2018); Bultinck  _et al._
(2020b); Zhang  _et al._ (2019); Liu  _et al._ (2019b); Wu  _et al._ (2019);
Thomson  _et al._ (2018); Dodaro  _et al._ (2018); Gonzalez and Stauber
(2019); Yuan and Fu (2018); Kang and Vafek (2019); Bultinck  _et al._ (2020a);
Seo  _et al._ (2019); Hejazi  _et al._ (2021); Khalaf  _et al._ (2021); Po
_et al._ (2018a); Xie  _et al._ (2020); Julku  _et al._ (2020); Hu  _et al._
(2019); Kang and Vafek (2020); Soejima  _et al._ (2020); Pixley and Andrei
(2019); König  _et al._ (2020); Christos  _et al._ (2020); Lewandowski  _et
al._ (2021); Kwan  _et al._ (2020, 2021); Xie and MacDonald (2020a); Liu and
Dai (2021); Cea and Guinea (2020); Zhang  _et al._ (2020); Liu  _et al._
(2021b); Da Liao  _et al._ (2019, 2021); Classen  _et al._ (2019); Kennes  _et
al._ (2018); Eugenio and Dağ (2020); Huang  _et al._ (2020, 2019); Guo  _et
al._ (2018); Ledwith  _et al._ (2020); Repellin  _et al._ (2020); Abouelkomsan
_et al._ (2020); Repellin and Senthil (2020); Vafek and Kang (2020); Fernandes
and Venderbos (2020); Wilson  _et al._ (2020); Wang  _et al._ (2021); Bernevig
_et al._ (2021c); Song  _et al._ (2021); Bernevig  _et al._ (2021a); Lian  _et
al._ (2021); Bernevig  _et al._ (2021b); Xie  _et al._ (2021); Zhang  _et al._
(2021); Vafek and Kang (2021); Cha  _et al._ (2021); Chichinadze  _et al._
(2020); Sheffer and Stern (2021); Kang  _et al._ (2021); Hofmann  _et al._
(2021); Calderón and Bascones (2020); Thomson and Alicea (2021) efforts have
been made to understand the nature of these exotic phases. Theoretical
challenges for understanding the correlation physics come from both the strong
interaction compared to relatively small band width as well as from the
topology Liu  _et al._ (2019a); Song  _et al._ (2019); Po  _et al._ (2019);
Song  _et al._ (2021); Tarnopolsky  _et al._ (2019); Hejazi  _et al._ (2019a),
which forbids a symmetric lattice description of the problem. The two flat
bands of MATBG posses strong topology in the presence of $C_{2z}T$ (time-
reversal followed by $C_{2z}$ rotation) and particle-hole ($P$) symmetries
Song  _et al._ (2021), which supersedes the earlier, $C_{2z}T$ symmetry-
protected fragile topology Song  _et al._ (2019); Po  _et al._ (2019). This
strong topology extends to the entire continuum BM model, and implies the
absence of a lattice model for any number of bands. The topology is also
responsible to exotic phases such as quantum anomalous Hall states Bultinck
_et al._ (2020b); Liu and Dai (2021); Lian  _et al._ (2021); Wu and Das Sarma
(2020); Bultinck  _et al._ (2020a); Kwan  _et al._ (2021) and fractional Chern
states Sheffer and Stern (2021); Repellin and Senthil (2020); Abouelkomsan
_et al._ (2020); Ledwith  _et al._ (2020).

Two types of complementary strategies have been proposed to resolve the
problem of the lattice description. One is to construct extended Hubbard
models Po  _et al._ (2018a); Kang and Vafek (2018, 2019); Koshino  _et al._
(2018); Po  _et al._ (2019); Vafek and Kang (2021); Zou  _et al._ (2018); Xu
_et al._ (2018); Yuan and Fu (2018), where either $C_{2z}T$ Kang and Vafek
(2018, 2019); Koshino  _et al._ (2018); Vafek and Kang (2021); Yuan and Fu
(2018) or $P$ Po  _et al._ (2019) becomes non-local in real space. The other
is to adopt a full momentum-space formalism Lian  _et al._ (2021); Bernevig
_et al._ (2021b); Xie  _et al._ (2021); Zhang  _et al._ (2021); Bultinck  _et
al._ (2020a); Cea and Guinea (2020); Zhang  _et al._ (2020); Hofmann  _et al._
(2021), where locality becomes hidden. (Besides the two strategies, some
phenomenological models are also proposed Efimkin and MacDonald (2018); Wu
_et al._ (2019); Xu and Balents (2018); Thomson  _et al._ (2018); Classen  _et
al._ (2019); Eugenio and Dağ (2020); Repellin  _et al._ (2020); Fernandes and
Venderbos (2020).) The real and momentum space strong coupling models
elucidated the nature of the correlated insulator states: they are
ferromagnets - sometimes carrying Chern numbers - in a large U(4) or
U(4)$\times$U(4) symmetry space that contains spin, valley and band quantum
number Kang and Vafek (2019); Bultinck  _et al._ (2020a); Bernevig  _et al._
(2021a). The dispersion of the excitations above the correlated insulators
Bernevig  _et al._ (2021b); Vafek and Kang (2021); Bultinck  _et al._ (2020a)
\- where superconductivity appears upon doping - is, despite being exact - not
physically understood.

In the current manuscript, nevertheless, we find it possible to write down a
fully symmetric model that has a simple real space picture, which, remarkably
and elegantly, solves the aforementioned puzzles. We reformulate and map the
interacting MATBG as an effective topological heavy fermion system, which
consists of local orbitals ($f$) centered at the AA-stacking regions and
delocalized topological conduction bands ($c$). The $f$-electrons are so
localized that they have an almost zero kinetic energy ($\sim 0.1$meV) but a
strong on-site Coulomb repulsion that we compute to be $\sim 60$meV. The
$c$-electrons carry the symmetry anomaly and have unbounded kinetic energies.
The actual flat bands of the BM model are from a hybridization ($\sim$20meV)
between the $f$\- and $c$-bands. The interacting Hamiltonian also couples the
$f$ and $c$ electrons through the presence of several types of interactions.
Using this model, the ground states Po  _et al._ (2018a); Padhi  _et al._
(2018); Bultinck  _et al._ (2020a); Kang and Vafek (2019); Liu and Dai (2021);
Ochi  _et al._ (2018); Venderbos and Fernandes (2018); Dodaro  _et al._
(2018); Seo  _et al._ (2019); Pixley and Andrei (2019); Xie and MacDonald
(2020b); Cea and Guinea (2020); Zhang  _et al._ (2020) and their topologies
can be understood in a simple, physical picture. The quasi-particle excitation
bandwidth can even be analytically determined.

Figure 1: Topological heavy fermion model. (a) A sketch of the moiré unit cell
of MATBG and its heavy fermion analog, where the local moments and itinerant
electrons are formed by the effective $f$-orbitals at the AA-stacking regions
and topological conduction bands ($c$), respectively. (b) The band structure
of the BM model at the magic angle $\theta=1.05^{\circ}$, where the moiré BZ
and high symmetry momenta are illustrated in the upper inset panel. The
overlaps between the Bloch states and the trial WFs are represented by the red
circles. The density profile of the constructed maximally localized WFs
($f$-orbitals) is shown in the lower inset panel. (c) Bands given by the
topological heavy fermion model (black lines) compared to the BM bands (blue
crosses). The $c$\- (blue) and $f$-bands (red) in the decoupled limit, where
$\gamma=v_{\star}^{\prime}=0$, are shown in the inset. Orange dashed lines
indicate evolution of energy levels as $f$-$c$ coupling is turned on.

Topological heavy fermion model — The single-valley BM model has the symmetry
of the magnetic space group $P6^{\prime}2^{\prime}2$, generated by $C_{3z}$,
$C_{2x}$, $C_{2z}T$, and moiré translations. (See Refs. Song  _et al._ (2019);
SM  for this group and its irreducible representations - irreps.) The energy
bands in the valley $\eta=+$ of the BM model are shown in Fig. 1(b), where the
bands are labeled by their irreps. Refs. Song  _et al._ (2019); Po  _et al._
(2019) showed that the irreps formed by the two flat bands, i.e.,
$\Gamma_{1}\oplus\Gamma_{2}$; $M_{1}\oplus M_{2}$; $K_{2}K_{3}$, are not
consistent with any local orbitals (band representations Bradlyn  _et al._
(2017)) and indicate a fragile Po  _et al._ (2018b); Ahn  _et al._ (2019);
Bouhon  _et al._ (2019); Cano  _et al._ (2018) topological obstruction to a
two-band lattice model. Here we resolve the fragile topology by involving
higher energy bands. Suppose we can “borrow” a $\Gamma_{3}$ irrep from higher
($\sim$20meV) energy bands and use it to replace the
$\Gamma_{1}\oplus\Gamma_{2}$ states; then the replaced irreps - $\Gamma_{3}$,
$M_{1}\oplus M_{2}$, $K_{2}K_{3}$ \- are consistent with $p_{x}\pm ip_{y}$
orbitals located at the triangular lattice. We hence introduce two trial
Gaussian-type Wannier functions (WFs) that transform as $p_{x}\pm ip_{y}$
orbitals under the crystalline symmetries. As indicated by the overlaps
between the trial WFs and the Bloch bands (Fig. 1(a)), the trial WFs are
supported by the flat band states at $\mathbf{k}$ away from $\Gamma_{M}$ and
by the lowest higher energy band states around $\Gamma_{M}$. Feeding the
overlaps into the program Wannier90 Marzari and Vanderbilt (1997); Souza  _et
al._ (2001); Pizzi  _et al._ (2020), we obtain the corresponding maximally
localized WFs, density profile of which is shown in Fig. 1(b) SM . (Similar
local states are also discussed using different methods in Refs. Liu  _et al._
(2019a); Calderón and Bascones (2020).) These WFs are extremely localized -
their nearest neighbor hoppings are about 0.1meV - and span 96% percent of the
flat bands.

To recover the irreps and topology of the middle two bands, we have to take
into account the remaining 4% states, without which the localized electrons
could not form a superconductor. To do this, we define the projector into the
WFs as $\mathbb{P}$, the projector into the lowest six bands (per spin valley)
as $\mathbb{I}$, and divide the low energy BM Hamiltonian $H_{BM}$ into four
parts: $H^{(f)}=\mathbb{P}H_{BM}\mathbb{P}$,
$H^{(c)}=\mathbb{Q}H_{BM}\mathbb{Q}$, $H^{(fc)}=\mathbb{P}H_{BM}\mathbb{Q}$,
$H^{(cf)}=H^{(fc)\dagger}$, where $\mathbb{Q}=\mathbb{I}-\mathbb{P}$,
$H^{(c)}$ is the remaining Hamiltonian, and $H^{(fc)}+h.c.$ is the coupling
between WFs and the remaining states. As the couplings between WFs are
extremely weak ($\sim$0.1meV) we find $H^{(f)}\approx 0$. Since the two states
in $\mathbb{P}$ form $\Gamma_{3}$ at $\Gamma_{M}$, the four states in
$\mathbb{Q}$ must form $\Gamma_{3}\oplus\Gamma_{1}\oplus\Gamma_{2}$ at
$\Gamma_{M}$ due to the irrep counting. Due to the crystalline and $P$
symmetries, $H^{(c)}$ in the valley $\eta$ takes the form SM

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

to linear order of $\mathbf{k}$, where the first two-by-two block is spanned
by the $\Gamma_{3}$ states and the second two-by-two block is spanned by the
$\Gamma_{1}\oplus\Gamma_{2}$ states. The $\Gamma_{1}$ and $\Gamma_{2}$ states
are split by the $M$ term (blue bands in Fig. 1(c)), while the $\Gamma_{3}$
states form a quadratic touching at $\mathbf{k}=0$, which is shown in Ref. SM
responsible to the symmetry anomaly Song  _et al._ (2021) jointly protected by
$C_{2z}T$ and $P$. The coupling $H^{(fc)}$ in the valley $\eta$ has the form

$H^{(fc,\eta)}(\mathbf{k})=\begin{pmatrix}\gamma\sigma_{0}+v_{\star}^{\prime}(\eta
k_{x}\sigma_{x}+k_{y}\sigma_{y}),&0_{2\times 2}\end{pmatrix}\ ,$ (2)

where the second block is computed to be extremely small and hence is omitted
and written as $0_{2\times 2}$. $H^{(fc,\eta)}$ will gap $H^{(c,\eta)}$, and
hence provides for both the single particle gap and for the flat band topology
of the BM model. Using a set of usually adopt parameters for MATBG, we find
$v_{\star}=-4.303\mathrm{eV}\cdot\mathrm{\mathring{A}}$, $M=3.697$meV,
$\gamma=-24.75$meV,
$v_{\star}^{\prime}=1.622\mathrm{eV}\cdot\mathrm{\mathring{A}}$.

Since the WFs and the remaining “$c$” degrees of freedom have localized and
plane-wave-like wave functions, respectively, we make the analogy with local
orbitals and conduction bands in heavy fermion systems. We refer to them as
local $f$-orbitals and (topological) conduction $c$-bands, respectively. We
use $f_{\mathbf{R}\alpha\eta s}$ ($\alpha=1,2$, $\eta=\pm$,
$s=\uparrow,\downarrow$) to represent the annihilation operator of the
$\alpha$-th WF of the valley $\eta$ and spin $s$ at the moiré unit cell
$\mathbf{R}$. We use $c_{\mathbf{k}a\eta s}$ ($a=1,2,3,4$) to represent the
annihilation operator of the $a$-th conduction band basis of the valley $\eta$
and spin $s$ at the moiré momentum $\mathbf{k}$. The single-particle
Hamiltonian can be written as

$\displaystyle{\hat{H}}_{0}=$
$\displaystyle\sum_{|\mathbf{k}|<\Lambda_{c}}\sum_{aa^{\prime}\eta
s}H^{(c,\eta)}_{aa^{\prime}}(\mathbf{k})c_{\mathbf{k}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta
s}+\frac{1}{\sqrt{N}}\sum_{\begin{subarray}{c}|\mathbf{k}|<\Lambda_{c}\\\
\mathbf{R}\end{subarray}}\sum_{\alpha a\eta s}\Big{(}$ $\displaystyle
e^{i\mathbf{k}\cdot\mathbf{R}-\frac{|\mathbf{k}|^{2}\lambda^{2}}{2}}H^{(fc,\eta)}_{\alpha
a}(\mathbf{k})f_{\mathbf{R}\alpha\eta s}^{\dagger}c_{\mathbf{k}a\eta
s}+h.c.\Big{)}\ ,$ (3)

where $\Lambda_{c}$ is the momentum cutoff for the $c$-electrons, $a_{M}$ is
the moiré lattice constant, $N$ is the number of moiré unit cell in the
system, and $\lambda$, which is found to be $0.3375a_{M}$, is a damping factor
proportional to the size of WFs. We plot the band structure of ${\hat{H}}_{0}$
in Fig. 1(c), where the splitting of the two $\Gamma_{3}$ states is given by
$2|\gamma|$ and the bandwidth of the two flat bands is given by $2M\approx
7.4\text{meV}$. The spectrum of ${\hat{H}}_{0}$ matches very well with the BM
model (Fig. 1(a)) in the energy range [-70meV, 70meV].

The U(4) symmetry — The projected model of MATBG Kang and Vafek (2019);
Bultinck  _et al._ (2020a); Bernevig  _et al._ (2021a) is found to possess a
U(4) symmetry if the kinetic energy of the flat bands is omitted. In the heavy
fermion basis, this U(4) symmetry can be realized by imposing the flat band
condition, i.e., $M=0$. (Note that $2|M|$ is the bandwidth of the flat bands.)
The U(4) moments of the $f$-electrons, $\Gamma_{3}$ $c$-electrons, and
$\Gamma_{1}\oplus\Gamma_{2}$ $c$-electrons are given by SM

$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}(\mathbf{R})=$
$\displaystyle\frac{\delta_{\xi,(-1)^{\alpha-1}\eta}}{2}A^{\mu\nu}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}$ (4)
$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}(\mathbf{q})=$
$\displaystyle\frac{\delta_{\xi,(-1)^{a-1}\eta}}{2N}A^{\mu\nu}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}+\mathbf{q}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s^{\prime}},\;(a=1,2)$
$\displaystyle{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}(\mathbf{q})=$
$\displaystyle\frac{\delta_{\xi,(-1)^{a-1}\eta}}{2N}B^{\mu\nu}_{a\eta
s,a^{\prime}\eta^{\prime}s^{\prime}}c_{\mathbf{k}+\mathbf{q}a\eta
s}^{\dagger}c_{\mathbf{k}a^{\prime}\eta^{\prime}s^{\prime}},\;(a=3,4)$

respectively, where repeated indices should be summed over and
$A^{\mu\nu},B^{\mu\nu}$ ($\mu,\nu=0,x,y,z$) are eight-by-eight matrices

$\displaystyle A^{\mu\nu}=$
$\displaystyle\\{\sigma_{0}\tau_{0}{\varsigma}_{\nu},\sigma_{y}\tau_{x}{\varsigma}_{\nu},\sigma_{y}\tau_{y}{\varsigma}_{\nu},\sigma_{0}\tau_{z}{\varsigma}_{\nu}\\}$
(5) $\displaystyle B^{\mu\nu}=$
$\displaystyle\\{\sigma_{0}\tau_{0}{\varsigma}_{\nu},-\sigma_{y}\tau_{x}{\varsigma}_{\nu},-\sigma_{y}\tau_{y}{\varsigma}_{\nu},\sigma_{0}\tau_{z}{\varsigma}_{\nu}\\}\
,$

with $\sigma_{0,x,y,z}$, $\tau_{0,x,y,z}$, ${\varsigma}_{0,x,y,z}$ being the
Pauli or identity matrices for the orbital, valley, and spin degrees of
freedom, respectively. The $\pm 1$ valued index $\xi$, equal to
$(-1)^{\alpha-1}\eta$ or $(-1)^{a-1}\eta$ in the moments, labels different
fundamental representations of the U(4) group. The global U(4) rotations are
generated by ${\hat{\Sigma}}_{\mu\nu}=\sum_{\xi=\pm
1}{\hat{\Sigma}}_{\mu\nu}^{(f,\xi)}+{\hat{\Sigma}}_{\mu\nu}^{(c\prime,\xi)}+{\hat{\Sigma}}_{\mu\nu}^{(c\prime\prime,\xi)}$.
Unlike the U(4) rotations found in Refs. Kang and Vafek (2019); Bultinck  _et
al._ (2020a); Bernevig  _et al._ (2021a), which only commute the projected
Hamiltonian into the flat bands, the U(4) rotations here commute with the full
Hamiltonian. (Generators of the U(4) or U(4)$\times$U(4) symmetry in the first
chiral limit Tarnopolsky  _et al._ (2019); Wang  _et al._ (2021) is also given
in Ref. SM .)

Interaction Hamiltonian — To obtain the interaction Hamiltonian in the heavy
fermion basis, we can first express the density operator $\rho(\mathbf{r})$ of
the BM model in terms of $f_{\mathbf{R}\alpha\eta s}$ and $c_{\mathbf{k}a\eta
s}$, and then substitute it into the Coulomb interaction,
$\rho(\mathbf{r})V(\mathbf{r}-\mathbf{r}^{\prime})\rho(\mathbf{r}^{\prime})$.
By evaluating the Coulomb integrals, we obtain the interaction Hamiltonian
resembling a periodic Anderson model with extra $f$-$c$ exchange interactions
SM ,

${\hat{H}}_{I}={\hat{H}}_{U_{1}}+{\hat{H}}_{J}+{\hat{H}}_{U_{2}}+{\hat{H}}_{V}+{\hat{H}}_{W}\
.$ (6)

${\hat{H}}_{U_{1}}=\frac{U_{1}}{2}\sum_{\mathbf{R}}:\rho^{f}_{\mathbf{R}}::\rho^{f}_{\mathbf{R}}:$
are the on-site interactions of $f$-electrons, where
$\rho^{f}_{\mathbf{R}}=\sum_{\alpha\eta s}f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha\eta s}$ is the $f$-electrons density and the
colon symbols represent the normal ordered operator with respect to the normal
state.

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

is a ferromagnetic exchange coupling between U(4) moments of $f$-electrons and
$\Gamma_{1}\oplus\Gamma_{2}$ $c$-electrons. Using practical parameters for
MATBG, we obtain $U_{1}=57.95$meV and $J=16.38$meV. The other three terms in
${\hat{H}}_{I}$ are: $H_{U_{2}}$ \- repulsion ($\sim 2.3$meV) between nearest
neighbor $f$-electrons, $H_{V}$ \- repulsion ($\sim 48$meV) between
$c$-electrons, $H_{W}$ \- repulsion ($\sim 47$meV) between $c$\- and
$f$-electrons.

As a widely adopted approximation in heavy fermion materials,
${\hat{H}}_{U_{2}}+{\hat{H}}_{V}+{\hat{H}}_{W}$ can be decoupled in the
Hartree channel due to the delocalized and localized natures of $c$\- and
$f$-electrons. Hence these terms only effectively shift the band energies of
$f$\- and $c$-bands. Then, $U_{1}$ \- the on-site repulsion of the
$f$-electrons - is by far the largest energy scale of the problem - more than
twice the hybridization ($\gamma$) and three times the exchange ($J$). In
Hartree-Fock (HF) calculations $U_{1}$ is found to be the source of
spontaneous symmetry-breakings.

Figure 2: The self-consistent HF bands upon the ground states at the fillings
$\nu=0,-1,-2,-3$. The color of the bands represent the contributing
components, wherein yellow represents the $f$-electron states and blue
represents the $c$-electron states.

Ground states — Since $U_{1}$ is much larger than the couplings ($\gamma$,
$J$, $v_{\star}^{\prime}(\eta k_{x}\sigma_{x}+k_{y}\sigma_{y})$) between $f$\-
and $c$-electrons, a reasonable guess of the ground states would be product
states of $f$-multiplets and the (gapless point) Fermi liquid state
($|\text{FS}\rangle$) of the half-filled $c$-electrons. We call such product
states “the parent states”. E.g., the parent valley-polarized (VP) state at
the charge neutrality ($\nu=0$) is

$|\text{VP}^{\nu=0}_{0}\rangle=\prod_{\mathbf{R}}\prod_{\alpha=1,2}\prod_{s=\uparrow\downarrow}f_{\mathbf{R},\alpha,+,s}^{\dagger}|\text{FS}\rangle\
.$ (8)

The parent Kramers inter-valley-coherent (K-IVC) state is a U(4)-rotation of
$|\text{VP}^{\nu=0}_{0}\rangle$ along the $\tau_{x}$-direction

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

Parent states at other integer fillings ($\nu=0,\pm 1,\pm 2,\pm 3$) can be
similarly constructed SM . They would be ground states of the Hamiltonian if
$\gamma$, $J$, $v_{\star}^{\prime}$ terms vanished; hybridization of $f$\- and
$c$-electrons will develop, i.e., $\langle f^{\dagger}c\rangle\neq 0$,
otherwise. The determination of ground states by self-consistent HF
calculation with initial states given by the parent states is given in Ref. SM
. The numerically found HF ground states at the integer fillings (Fig. 2) are
fully consistent with those in Ref. Lian  _et al._ (2021).

The parent states are so good initial states for the HF calculations that the
one-shot HF is already qualitatively same as the self-consistent HF (see Fig.
3). Thanks to the simplicity of the heavy fermion model, the one-shot energies
can be analytically estimated and we are able to derive two rules for the
ground states SM . _First_ , in the parent state, $f$-electrons at each site
tend to be symmetric under permutation of U(4) indices to save the Coulomb
energy (Hunds’ rule). Both Eqs. 8 and MATBG as Topological Heavy Fermion: I.
Exact Mapping and Correlated Insulators satisfy the first rule. Second, for
U(4)-related states at a given integer filling $\nu$, the state that minimizes
${\hat{H}}_{M}+{\hat{H}}_{J}$ is the ground state, where ${\hat{H}}_{M}$ is
the U(4)-breaking $M$ term in ${\hat{H}}_{0}$ (Eq. 1). This energy can be
estimated by the lowest $\nu+4$ levels of the mean field Hamiltonian
$H^{(\Gamma_{1}\oplus\Gamma_{2})}$ spanned by the $\Gamma_{1}\oplus\Gamma_{2}$
basis of the $c$-bands at $\mathbf{k}=0$, which reads (up to constants)

$H^{(\Gamma_{1}\oplus\Gamma_{2})}=M\sigma_{x}\tau_{0}{\varsigma}_{0}-\frac{J}{2}(\tau_{z}\overline{O}^{fT}\tau_{z}+\sigma_{z}\overline{O}^{fT}\sigma_{z})\
.$ (10)

Here $\overline{O}^{f}_{\alpha\eta
s,\alpha^{\prime}\eta^{\prime}s^{\prime}}=\langle f_{\mathbf{R}\alpha\eta
s}^{\dagger}f_{\mathbf{R}\alpha^{\prime}\eta^{\prime}s^{\prime}}\rangle-\frac{1}{2}\delta_{\alpha\alpha^{\prime}}\delta_{\eta\eta^{\prime}}\delta_{ss^{\prime}}$
is the density matrix of the local $f$-orbitals with homogeneity assumed. We
have assumed that, for all the integer fillings, the eight lowest energy
levels (closest to the gap) are always contributed by the
$\Gamma_{1}\oplus\Gamma_{2}$ $c$-band basis - which are part of the flat bands
- hence we only need to look at the $\Gamma_{1}\oplus\Gamma_{2}$ subspace.
This assumption is fully justified by previous studies based on two-band
projected Hamiltonian, where only $\Gamma_{1}\oplus\Gamma_{2}$ basis exists,
and will become clearer after we discuss the charge $\pm 1$ excitations.

We now apply the second rule to Eq. 8 and (MATBG as Topological Heavy Fermion:
I. Exact Mapping and Correlated Insulators) to determine the one-shot state of
lowest energy. $\overline{O}^{f}$ matrices given by Eqs. 8 and MATBG as
Topological Heavy Fermion: I. Exact Mapping and Correlated Insulators are
$\frac{1}{2}\sigma_{0}\tau_{z}{\varsigma}_{0}$ and
$-\frac{1}{2}\sigma_{y}\tau_{y}{\varsigma}_{0}$, respectively; the resulted
lowest four levels of $H^{(\Gamma_{1}\oplus\Gamma_{2})}$ are $\pm M-J/2$ (each
2-fold) and $-\sqrt{M^{2}+J^{2}/4}$ (4-fold), respectively. It is direct to
verify that the latter (K-IVC) has a lower energy. Applying the two rules to
parent states at other fillings, we obtain consistent results with the full
numerical calculations in Refs. Lian  _et al._ (2021); Xie  _et al._ (2021).
We also obtain an analytical expression for the Chern numbers of ground states
SM .

Figure 3: The one-shot HF bands of the ground states at the fillings
$\nu=0,-1$. The red solid bands are the quasi-particle bands of the decoupled
Hamiltonian, where $\gamma=v_{\star}^{\prime}=J=0$. The horizontal and
dispersive red bands are of the $f$\- and $c$-electrons, respectively. The
touching point of the dispersive red bands at $\Gamma_{M}$ is quadratic, while
since $M$ is small, it may look like linear. The one-shot bands can be
understood as a result of hybridization between $f$\- and $c$-electrons.

Charge $\pm 1$ excitations — As shown in Figs. 2 and 3 and in Refs. Bernevig
_et al._ (2021b); Bultinck  _et al._ (2020a); Vafek and Kang (2021); Cea and
Guinea (2020); Zhang  _et al._ (2020); Kang  _et al._ (2021), at $\mathbf{k}$
away from $\Gamma_{M}$, the quasi-particle bands have a large gap ($\sim
U_{1}$) and are relatively flat; at $\mathbf{k}$ around $\Gamma_{M}$, the
bands have significant dip. Such features are found related to the topology of
the two flat-bands Bernevig  _et al._ (2021b); Vafek and Kang (2021) but have
not yet been quantitatively understood. The heavy fermion model provides a
natural explanation to these features. We first consider the decoupled limit
($\gamma=v_{\star}^{\prime}=J=0$) at $\nu=0$, where the $f$-electron bands are
flat and have a (charge) gap $U_{1}$, and the $c$-electron bands are given by
$H^{(c,\eta)}$ (Fig. 3(a)). Tuning on $\gamma,v_{\star}^{\prime},J$ then
yields the one-shot quasi-particle bands. At $\mathbf{k}=0$, $\gamma$ gaps out
the $\Gamma_{3}$ $c$-bands, and $J$ further gaps out the
$\Gamma_{1}\oplus\Gamma_{2}$ $c$-bands. As the splitting of
$\Gamma_{1}-\Gamma_{2}$ is smaller than that of the $\Gamma_{3}$, the lowest
excitations will carry $\Gamma_{1},\Gamma_{2}$ representations, matching Refs.
Bernevig  _et al._ (2021b); Vafek and Kang (2021); Bultinck  _et al._ (2020a)
and, according to the discussion after Eq. 10, equals to
$2\sqrt{M^{2}+J^{2}/4}$ and $|J-2M|$ for K-IVC and VP states, respectively. At
$\mathbf{k}\neq 0$, the $v_{\star}^{\prime}$ term hybridizes the flat
$f$-bands and dispersive $c$-bands. For large $\mathbf{k}$, where the
$c$-bands have very high kinetic energies, the hybridization is relatively
weak and the gap is still approximately $U_{1}$. Thus the shape of the quasi-
particle bands is explained, and its bandwidth is approximately given by
$(U_{1}-J)/2$ when $M$ is small. As discussed in Ref. SM , the feature that
the larger ($\sim U_{1}$) and smaller ($\sim J$) gaps are contributed by $f$\-
and $c$-electrons, respectively, is reflected in the STM spectra and Landau
levels at different regions (AA or AB sites) of MATBG.

At nonzero fillings, the quasi-particle bands can also be understood as
hybridized flat $f$-bands and dispersive $c$-bands, except that the $f$\- and
$c$-bands may feel different effective chemical potentials due to the density-
density interactions between them. For example, at $\nu=-1$, the upper branch
of the $f$-bands is shifted to an energy close to the quadratic touching of
the $c$-bands (Fig. 3(b)) SM . Thus one of the hybridized bands is extremely
flat.

Discussion — The coexistence of quantum-dot-like behavior Xie  _et al._
(2019); Wong  _et al._ (2020) and superconductivity Cao  _et al._ (2018b); Lu
_et al._ (2019); Yankowitz  _et al._ (2019); Saito  _et al._ (2020); Stepanov
_et al._ (2020) may now be understood - they come from two different types
($f$ and $c$) of carriers. In fact, inspired by the pomeranchuk effect
experiments Saito  _et al._ (2021b); Rozen  _et al._ (2021) and strange metal
behavior Cao  _et al._ (2020); Polshyn  _et al._ (2019), authors of Refs.
Saito  _et al._ (2021b); Rozen  _et al._ (2021); Cha  _et al._ (2021) also
conjecture the possibility of coexistence of local momenta and itinerant
electrons. (The heavy fermion theory analog may also exist in other twisted
materials Ramires and Lado (2021).) Our paper _derives_ and shows the exact
mapping of MATBG to such a heavy-fermion type model. As such, the machinery of
heavy fermions Si and Steglich (2010); Gegenwart  _et al._ (2008); Coleman
(1984); Dzero  _et al._ (2010); Tsvelik and Wiegmann (1983); Werner  _et al._
(2006); Lu  _et al._ (2013); Weng  _et al._ (2014); Kotliar  _et al._ (2006);
Emery and Kivelson (1992); Freyer  _et al._ (2018); Furusaki and Nagaosa
(1994); Cassanello and Fradkin (1996); Martin  _et al._ (1997); Chakravarty
and Rudnick (1995) can now be applied, for the first time, to MATBG. We
speculate that it will lead to pairing Guinea and Walet (2018); You and
Vishwanath (2019); Isobe  _et al._ (2018); Liu  _et al._ (2018); Gonzalez and
Stauber (2019); Khalaf  _et al._ (2021); Christos  _et al._ (2020); König  _et
al._ (2020); Kennes  _et al._ (2018); Huang  _et al._ (2019); Guo  _et al._
(2018); Chichinadze  _et al._ (2020); Lian  _et al._ (2019); Wu  _et al._
(2018); Lewandowski  _et al._ (2021) in nontrivial gap channels.

###### Acknowledgements.

We thank X. Dai, O Vafek, A. Yazdani, P. Coleman, EA. Kim, Q. Si, R Fernandez,
P. Jarillo-Herrero and D. Efetov for discussions. B. A. B.’s work was
primarily supported by the DOE Grant No. DE-SC0016239, the Schmidt Fund for
Innovative Research, Simons Investigator Grant No. 404513, and the Packard
Foundation. Z.-D. S. was supported by ONR No. N00014-20-1-2303, the National
Key Research and Development Program of China (No. 2021YFA1401900), NSF-MRSEC
No. DMR-1420541, and Gordon and Betty Moore Foundation through Grant GBMF8685
towards the Princeton theory program. B.A.B. also acknowledges support from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (Grant Agreement No. 101020833).

## References

  * Kang and Vafek (2019) Jian Kang and Oskar Vafek, “Strong Coupling Phases of Partially Filled Twisted Bilayer Graphene Narrow Bands,” Physical Review Letters 122, 246401 (2019).
  * Bultinck  _et al._ (2020a) Nick Bultinck, Eslam Khalaf, Shang Liu, Shubhayu Chatterjee, Ashvin Vishwanath, and Michael P. Zaletel, “Ground state and hidden symmetry of magic-angle graphene at even integer filling,” Phys. Rev. X 10, 031034 (2020a).
  * Vafek and Kang (2020) Oskar Vafek and Jian Kang, “Renormalization group study of hidden symmetry in twisted bilayer graphene with coulomb interactions,” Phys. Rev. Lett. 125, 257602 (2020).
  * Bernevig  _et al._ (2021a) B. Andrei Bernevig, Zhi-Da Song, Nicolas Regnault, and Biao Lian, “Twisted bilayer graphene. III. Interacting Hamiltonian and exact symmetries,” Physical Review B 103, 205413 (2021a), publisher: American Physical Society.
  * Lian  _et al._ (2021) Biao Lian, Zhi-Da Song, Nicolas Regnault, Dmitri K. Efetov, Ali Yazdani, and B. Andrei Bernevig, “Twisted bilayer graphene. IV. Exact insulator ground states and phase diagram,” Physical Review B 103, 205414 (2021), publisher: American Physical Society.
  * Bernevig  _et al._ (2021b) B. Andrei Bernevig, Biao Lian, Aditya Cowsik, Fang Xie, Nicolas Regnault, and Zhi-Da Song, “Twisted bilayer graphene. V. Exact analytic many-body excitations in Coulomb Hamiltonians: Charge gap, Goldstone modes, and absence of Cooper pairing,” Physical Review B 103, 205415 (2021b), publisher: American Physical Society.
  * Vafek and Kang (2021) Oskar Vafek and Jian Kang, “Lattice model for the coulomb interacting chiral limit of the magic angle twisted bilayer graphene: symmetries, obstructions and excitations,” arXiv preprint arXiv:2106.05670 (2021).
  * Cao  _et al._ (2018a) Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L. Tomarken, Jason Y. Luo, Javier D. Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, Ray C. Ashoori, and Pablo Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature 556, 80–84 (2018a).
  * Cao  _et al._ (2018b) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43–50 (2018b).
  * Bistritzer and MacDonald (2011) Rafi Bistritzer and Allan H. MacDonald, “Moiré bands in twisted double-layer graphene,” Proceedings of the National Academy of Sciences 108, 12233–12237 (2011).
  * Lu  _et al._ (2019) Xiaobo Lu, Petr Stepanov, Wei Yang, Ming Xie, Mohammed Ali Aamir, Ipsita Das, Carles Urgell, Kenji Watanabe, Takashi Taniguchi, Guangyu Zhang, Adrian Bachtold, Allan H. MacDonald, and Dmitri K. Efetov, “Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene,” Nature 574, 653–657 (2019).
  * Yankowitz  _et al._ (2019) Matthew Yankowitz, Shaowen Chen, Hryhoriy Polshyn, Yuxuan Zhang, K Watanabe, T Taniguchi, David Graf, Andrea F Young, and Cory R Dean, “Tuning superconductivity in twisted bilayer graphene,” Science 363, 1059–1064 (2019).
  * Sharpe  _et al._ (2019) Aaron L. Sharpe, Eli J. Fox, Arthur W. Barnard, Joe Finney, Kenji Watanabe, Takashi Taniguchi, M. A. Kastner, and David Goldhaber-Gordon, “Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene,” Science 365, 605–608 (2019).
  * Saito  _et al._ (2020) Yu Saito, Jingyuan Ge, Kenji Watanabe, Takashi Taniguchi, and Andrea F. Young, “Independent superconductors and correlated insulators in twisted bilayer graphene,” Nature Physics 16, 926–930 (2020).
  * Stepanov  _et al._ (2020) Petr Stepanov, Ipsita Das, Xiaobo Lu, Ali Fahimniya, Kenji Watanabe, Takashi Taniguchi, Frank H. L. Koppens, Johannes Lischner, Leonid Levitov, and Dmitri K. Efetov, “Untying the insulating and superconducting orders in magic-angle graphene,” Nature 583, 375–378 (2020).
  * Liu  _et al._ (2021a) Xiaoxue Liu, Zhi Wang, Kenji Watanabe, Takashi Taniguchi, Oskar Vafek, and JIA Li, “Tuning electron correlation in magic-angle twisted bilayer graphene using coulomb screening,” Science 371, 1261–1265 (2021a).
  * Arora  _et al._ (2020) Harpreet Singh Arora, Robert Polski, Yiran Zhang, Alex Thomson, Youngjoon Choi, Hyunjin Kim, Zhong Lin, Ilham Zaky Wilson, Xiaodong Xu, Jiun-Haw Chu, and et al., “Superconductivity in metallic twisted bilayer graphene stabilized by wse2,” Nature 583, 379–384 (2020).
  * Serlin  _et al._ (2019) M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, “Intrinsic quantized anomalous hall effect in a moiré heterostructure,” Science 367, 900–903 (2019).
  * Cao  _et al._ (2020) Yuan Cao, Debanjan Chowdhury, Daniel Rodan-Legrain, Oriol Rubies-Bigorda, Kenji Watanabe, Takashi Taniguchi, T. Senthil, and Pablo Jarillo-Herrero, “Strange metal in magic-angle graphene with near planckian dissipation,” Phys. Rev. Lett. 124, 076801 (2020).
  * Polshyn  _et al._ (2019) Hryhoriy Polshyn, Matthew Yankowitz, Shaowen Chen, Yuxuan Zhang, K. Watanabe, T. Taniguchi, Cory R. Dean, and Andrea F. Young, “Large linear-in-temperature resistivity in twisted bilayer graphene,” Nature Physics 15, 1011–1016 (2019).
  * Xie  _et al._ (2019) Yonglong Xie, Biao Lian, Berthold Jäck, Xiaomeng Liu, Cheng-Li Chiu, Kenji Watanabe, Takashi Taniguchi, B Andrei Bernevig, and Ali Yazdani, “Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene,” Nature 572, 101–105 (2019).
  * Choi  _et al._ (2019) Youngjoon Choi, Jeannette Kemmer, Yang Peng, Alex Thomson, Harpreet Arora, Robert Polski, Yiran Zhang, Hechen Ren, Jason Alicea, Gil Refael, and et al., “Electronic correlations in twisted bilayer graphene near the magic angle,” Nature Physics 15, 1174–1180 (2019).
  * Kerelsky  _et al._ (2019) Alexander Kerelsky, Leo J. McGilly, Dante M. Kennes, Lede Xian, Matthew Yankowitz, Shaowen Chen, K. Watanabe, T. Taniguchi, James Hone, Cory Dean, and et al., “Maximized electron interactions at the magic angle in twisted bilayer graphene,” Nature 572, 95–100 (2019).
  * Jiang  _et al._ (2019) Yuhang Jiang, Xinyuan Lai, Kenji Watanabe, Takashi Taniguchi, Kristjan Haule, Jinhai Mao, and Eva Y. Andrei, “Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene,” Nature 573, 91–95 (2019).
  * Wong  _et al._ (2020) Dillon Wong, Kevin P. Nuckolls, Myungchul Oh, Biao Lian, Yonglong Xie, Sangjun Jeon, Kenji Watanabe, Takashi Taniguchi, B. Andrei Bernevig, and Ali Yazdani, “Cascade of electronic transitions in magic-angle twisted bilayer graphene,” Nature 582, 198–202 (2020).
  * Zondiner  _et al._ (2020) U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, Ady Stern, and et al., “Cascade of phase transitions and dirac revivals in magic-angle graphene,” Nature 582, 203–208 (2020).
  * Nuckolls  _et al._ (2020) Kevin P Nuckolls, Myungchul Oh, Dillon Wong, Biao Lian, Kenji Watanabe, Takashi Taniguchi, B Andrei Bernevig, and Ali Yazdani, “Strongly correlated chern insulators in magic-angle twisted bilayer graphene,” Nature 588, 610–615 (2020).
  * Choi  _et al._ (2021) Youngjoon Choi, Hyunjin Kim, Yang Peng, Alex Thomson, Cyprian Lewandowski, Robert Polski, Yiran Zhang, Harpreet Singh Arora, Kenji Watanabe, Takashi Taniguchi, _et al._ , “Correlation-driven topological phases in magic-angle twisted bilayer graphene,” Nature 589, 536–541 (2021).
  * Saito  _et al._ (2021a) Yu Saito, Jingyuan Ge, Louk Rademaker, Kenji Watanabe, Takashi Taniguchi, Dmitry A Abanin, and Andrea F Young, “Hofstadter subband ferromagnetism and symmetry-broken chern insulators in twisted bilayer graphene,” Nature Physics 17, 478–481 (2021a).
  * Das  _et al._ (2021) Ipsita Das, Xiaobo Lu, Jonah Herzog-Arbeitman, Zhi-Da Song, Kenji Watanabe, Takashi Taniguchi, B Andrei Bernevig, and Dmitri K Efetov, “Symmetry-broken chern insulators and rashba-like landau-level crossings in magic-angle bilayer graphene,” Nature Physics 17, 710–714 (2021).
  * Wu  _et al._ (2021) Shuang Wu, Zhenyuan Zhang, K Watanabe, T Taniguchi, and Eva Y Andrei, “Chern insulators, van hove singularities and topological flat bands in magic-angle twisted bilayer graphene,” Nature materials 20, 488–494 (2021).
  * Park  _et al._ (2021) Jeong Min Park, Yuan Cao, Kenji Watanabe, Takashi Taniguchi, and Pablo Jarillo-Herrero, “Flavour hund’s coupling, chern gaps and charge diffusivity in moiré graphene,” Nature 592, 43–48 (2021).
  * Saito  _et al._ (2021b) Yu Saito, Fangyuan Yang, Jingyuan Ge, Xiaoxue Liu, Takashi Taniguchi, Kenji Watanabe, JIA Li, Erez Berg, and Andrea F Young, “Isospin pomeranchuk effect in twisted bilayer graphene,” Nature 592, 220–224 (2021b).
  * Rozen  _et al._ (2021) Asaf Rozen, Jeong Min Park, Uri Zondiner, Yuan Cao, Daniel Rodan-Legrain, Takashi Taniguchi, Kenji Watanabe, Yuval Oreg, Ady Stern, Erez Berg, _et al._ , “Entropic evidence for a pomeranchuk effect in magic-angle graphene,” Nature 592, 214–219 (2021).
  * Lu  _et al._ (2021) Xiaobo Lu, Biao Lian, Gaurav Chaudhary, Benjamin A. Piot, Giulio Romagnoli, Kenji Watanabe, Takashi Taniguchi, Martino Poggio, Allan H. MacDonald, B. Andrei Bernevig, and Dmitri K. Efetov, “Multiple flat bands and topological hofstadter butterfly in twisted bilayer graphene close to the second magic angle,” Proceedings of the National Academy of Sciences 118 (2021), 10.1073/pnas.2100006118.
  * Tarnopolsky  _et al._ (2019) Grigory Tarnopolsky, Alex Jura Kruchkov, and Ashvin Vishwanath, “Origin of Magic Angles in Twisted Bilayer Graphene,” Physical Review Letters 122, 106405 (2019).
  * Zou  _et al._ (2018) Liujun Zou, Hoi Chun Po, Ashvin Vishwanath, and T. Senthil, “Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions,” Phys. Rev. B 98, 085435 (2018).
  * Liu  _et al._ (2019a) Jianpeng Liu, Junwei Liu, and Xi Dai, “Pseudo landau level representation of twisted bilayer graphene: Band topology and implications on the correlated insulating phase,” Physical Review B 99, 155415 (2019a).
  * Efimkin and MacDonald (2018) Dmitry K. Efimkin and Allan H. MacDonald, “Helical network model for twisted bilayer graphene,” Phys. Rev. B 98, 035404 (2018).
  * Kang and Vafek (2018) Jian Kang and Oskar Vafek, “Symmetry, Maximally Localized Wannier States, and a Low-Energy Model for Twisted Bilayer Graphene Narrow Bands,” Phys. Rev. X 8, 031088 (2018).
  * Song  _et al._ (2019) Zhida Song, Zhijun Wang, Wujun Shi, Gang Li, Chen Fang, and B. Andrei Bernevig, “All Magic Angles in Twisted Bilayer Graphene are Topological,” Physical Review Letters 123, 036401 (2019).
  * Po  _et al._ (2019) Hoi Chun Po, Liujun Zou, T. Senthil, and Ashvin Vishwanath, “Faithful tight-binding models and fragile topology of magic-angle bilayer graphene,” Physical Review B 99, 195455 (2019).
  * Hejazi  _et al._ (2019a) Kasra Hejazi, Chunxiao Liu, Hassan Shapourian, Xiao Chen, and Leon Balents, “Multiple topological transitions in twisted bilayer graphene near the first magic angle,” Phys. Rev. B 99, 035111 (2019a).
  * Padhi  _et al._ (2018) Bikash Padhi, Chandan Setty, and Philip W Phillips, “Doped twisted bilayer graphene near magic angles: proximity to wigner crystallization, not mott insulation,” Nano letters 18, 6175–6180 (2018).
  * Lian  _et al._ (2020) Biao Lian, Fang Xie, and B. Andrei Bernevig, “Landau level of fragile topology,” Phys. Rev. B 102, 041402 (2020).
  * Hejazi  _et al._ (2019b) Kasra Hejazi, Chunxiao Liu, and Leon Balents, “Landau levels in twisted bilayer graphene and semiclassical orbits,” Physical Review B 100 (2019b), 10.1103/physrevb.100.035115.
  * Padhi  _et al._ (2020) Bikash Padhi, Apoorv Tiwari, Titus Neupert, and Shinsei Ryu, “Transport across twist angle domains in moiré graphene,” (2020), arXiv:2005.02406 [cond-mat.mes-hall] .
  * Xu and Balents (2018) Cenke Xu and Leon Balents, “Topological superconductivity in twisted multilayer graphene,” Physical review letters 121, 087001 (2018).
  * Koshino  _et al._ (2018) Mikito Koshino, Noah F. Q. Yuan, Takashi Koretsune, Masayuki Ochi, Kazuhiko Kuroki, and Liang Fu, “Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene,” Phys. Rev. X 8, 031087 (2018).
  * Ochi  _et al._ (2018) Masayuki Ochi, Mikito Koshino, and Kazuhiko Kuroki, “Possible correlated insulating states in magic-angle twisted bilayer graphene under strongly competing interactions,” Phys. Rev. B 98, 081102 (2018).
  * Xu  _et al._ (2018) Xiao Yan Xu, K. T. Law, and Patrick A. Lee, “Kekulé valence bond order in an extended hubbard model on the honeycomb lattice with possible applications to twisted bilayer graphene,” Phys. Rev. B 98, 121406 (2018).
  * Guinea and Walet (2018) Francisco Guinea and Niels R. Walet, “Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers,” Proceedings of the National Academy of Sciences 115, 13174–13179 (2018).
  * Venderbos and Fernandes (2018) Jörn W. F. Venderbos and Rafael M. Fernandes, “Correlations and electronic order in a two-orbital honeycomb lattice model for twisted bilayer graphene,” Phys. Rev. B 98, 245103 (2018).
  * You and Vishwanath (2019) Y.-Z. You and A. Vishwanath, “Superconductivity from Valley Fluctuations and Approximate SO(4) Symmetry in a Weak Coupling Theory of Twisted Bilayer Graphene,” npj Quantum Materials 4, 16 (2019).
  * Wu and Das Sarma (2020) Fengcheng Wu and Sankar Das Sarma, “Collective excitations of quantum anomalous hall ferromagnets in twisted bilayer graphene,” Physical Review Letters 124 (2020), 10.1103/physrevlett.124.046403.
  * Lian  _et al._ (2019) Biao Lian, Zhijun Wang, and B. Andrei Bernevig, “Twisted bilayer graphene: A phonon-driven superconductor,” Phys. Rev. Lett. 122, 257002 (2019).
  * Wu  _et al._ (2018) Fengcheng Wu, A. H. MacDonald, and Ivar Martin, “Theory of phonon-mediated superconductivity in twisted bilayer graphene,” Phys. Rev. Lett. 121, 257001 (2018).
  * Isobe  _et al._ (2018) Hiroki Isobe, Noah F. Q. Yuan, and Liang Fu, “Unconventional superconductivity and density waves in twisted bilayer graphene,” Phys. Rev. X 8, 041041 (2018).
  * Liu  _et al._ (2018) Cheng-Cheng Liu, Li-Da Zhang, Wei-Qiang Chen, and Fan Yang, “Chiral spin density wave and d+ i d superconductivity in the magic-angle-twisted bilayer graphene,” Physical review letters 121, 217001 (2018).
  * Bultinck  _et al._ (2020b) Nick Bultinck, Shubhayu Chatterjee, and Michael P. Zaletel, “Mechanism for anomalous hall ferromagnetism in twisted bilayer graphene,” Phys. Rev. Lett. 124, 166601 (2020b).
  * Zhang  _et al._ (2019) Ya-Hui Zhang, Dan Mao, Yuan Cao, Pablo Jarillo-Herrero, and T Senthil, “Nearly flat chern bands in moiré superlattices,” Physical Review B 99, 075127 (2019).
  * Liu  _et al._ (2019b) Jianpeng Liu, Zhen Ma, Jinhua Gao, and Xi Dai, “Quantum valley hall effect, orbital magnetism, and anomalous hall effect in twisted multilayer graphene systems,” Physical Review X 9, 031021 (2019b).
  * Wu  _et al._ (2019) Xiao-Chuan Wu, Chao-Ming Jian, and Cenke Xu, “Coupled-wire description of the correlated physics in twisted bilayer graphene,” Physical Review B 99 (2019), 10.1103/physrevb.99.161405.
  * Thomson  _et al._ (2018) Alex Thomson, Shubhayu Chatterjee, Subir Sachdev, and Mathias S. Scheurer, “Triangular antiferromagnetism on the honeycomb lattice of twisted bilayer graphene,” Physical Review B 98 (2018), 10.1103/physrevb.98.075109.
  * Dodaro  _et al._ (2018) John F Dodaro, Steven A Kivelson, Yoni Schattner, Xiao-Qi Sun, and Chao Wang, “Phases of a phenomenological model of twisted bilayer graphene,” Physical Review B 98, 075154 (2018).
  * Gonzalez and Stauber (2019) Jose Gonzalez and Tobias Stauber, “Kohn-luttinger superconductivity in twisted bilayer graphene,” Physical review letters 122, 026801 (2019).
  * Yuan and Fu (2018) Noah FQ Yuan and Liang Fu, “Model for the metal-insulator transition in graphene superlattices and beyond,” Physical Review B 98, 045103 (2018).
  * Seo  _et al._ (2019) Kangjun Seo, Valeri N. Kotov, and Bruno Uchoa, “Ferromagnetic mott state in twisted graphene bilayers at the magic angle,” Phys. Rev. Lett. 122, 246402 (2019).
  * Hejazi  _et al._ (2021) Kasra Hejazi, Xiao Chen, and Leon Balents, “Hybrid wannier chern bands in magic angle twisted bilayer graphene and the quantized anomalous hall effect,” Phys. Rev. Research 3, 013242 (2021).
  * Khalaf  _et al._ (2021) Eslam Khalaf, Shubhayu Chatterjee, Nick Bultinck, Michael P Zaletel, and Ashvin Vishwanath, “Charged skyrmions and topological origin of superconductivity in magic-angle graphene,” Science advances 7, eabf5299 (2021).
  * Po  _et al._ (2018a) Hoi Chun Po, Liujun Zou, Ashvin Vishwanath, and T. Senthil, “Origin of Mott Insulating Behavior and Superconductivity in Twisted Bilayer Graphene,” Physical Review X 8, 031089 (2018a).
  * Xie  _et al._ (2020) Fang Xie, Zhida Song, Biao Lian, and B. Andrei Bernevig, “Topology-bounded superfluid weight in twisted bilayer graphene,” Phys. Rev. Lett. 124, 167002 (2020).
  * Julku  _et al._ (2020) A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkilä, and P. Törmä, “Superfluid weight and berezinskii-kosterlitz-thouless transition temperature of twisted bilayer graphene,” Physical Review B 101 (2020), 10.1103/physrevb.101.060505.
  * Hu  _et al._ (2019) Xiang Hu, Timo Hyart, Dmitry I. Pikulin, and Enrico Rossi, “Geometric and conventional contribution to the superfluid weight in twisted bilayer graphene,” Phys. Rev. Lett. 123, 237002 (2019).
  * Kang and Vafek (2020) Jian Kang and Oskar Vafek, “Non-abelian dirac node braiding and near-degeneracy of correlated phases at odd integer filling in magic-angle twisted bilayer graphene,” Phys. Rev. B 102, 035161 (2020).
  * Soejima  _et al._ (2020) Tomohiro Soejima, Daniel E. Parker, Nick Bultinck, Johannes Hauschild, and Michael P. Zaletel, “Efficient simulation of moiré materials using the density matrix renormalization group,” Phys. Rev. B 102, 205111 (2020).
  * Pixley and Andrei (2019) Jed H. Pixley and Eva Y. Andrei, “Ferromagnetism in magic-angle graphene,” Science 365, 543–543 (2019).
  * König  _et al._ (2020) E. J. König, Piers Coleman, and A. M. Tsvelik, “Spin magnetometry as a probe of stripe superconductivity in twisted bilayer graphene,” Phys. Rev. B 102, 104514 (2020).
  * Christos  _et al._ (2020) Maine Christos, Subir Sachdev, and Mathias S. Scheurer, “Superconductivity, correlated insulators, and wess–zumino–witten terms in twisted bilayer graphene,” Proceedings of the National Academy of Sciences 117, 29543–29554 (2020).
  * Lewandowski  _et al._ (2021) Cyprian Lewandowski, Debanjan Chowdhury, and Jonathan Ruhman, “Pairing in magic-angle twisted bilayer graphene: Role of phonon and plasmon umklapp,” Phys. Rev. B 103, 235401 (2021).
  * Kwan  _et al._ (2020) Yves H. Kwan, S. A. Parameswaran, and S. L. Sondhi, “Twisted bilayer graphene in a parallel magnetic field,” Phys. Rev. B 101, 205116 (2020).
  * Kwan  _et al._ (2021) Yves H. Kwan, Yichen Hu, Steven H. Simon, and S. A. Parameswaran, “Exciton band topology in spontaneous quantum anomalous hall insulators: Applications to twisted bilayer graphene,” Phys. Rev. Lett. 126, 137601 (2021).
  * Xie and MacDonald (2020a) Ming Xie and A. H. MacDonald, “Nature of the correlated insulator states in twisted bilayer graphene,” Phys. Rev. Lett. 124, 097601 (2020a).
  * Liu and Dai (2021) Jianpeng Liu and Xi Dai, “Theories for the correlated insulating states and quantum anomalous hall effect phenomena in twisted bilayer graphene,” Phys. Rev. B 103, 035427 (2021).
  * Cea and Guinea (2020) Tommaso Cea and Francisco Guinea, “Band structure and insulating states driven by coulomb interaction in twisted bilayer graphene,” Phys. Rev. B 102, 045107 (2020).
  * Zhang  _et al._ (2020) Yi Zhang, Kun Jiang, Ziqiang Wang, and Fuchun Zhang, “Correlated insulating phases of twisted bilayer graphene at commensurate filling fractions: A hartree-fock study,” Phys. Rev. B 102, 035136 (2020).
  * Liu  _et al._ (2021b) Shang Liu, Eslam Khalaf, Jong Yeon Lee, and Ashvin Vishwanath, “Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality,” Phys. Rev. Research 3, 013033 (2021b).
  * Da Liao  _et al._ (2019) Yuan Da Liao, Zi Yang Meng, and Xiao Yan Xu, “Valence bond orders at charge neutrality in a possible two-orbital extended hubbard model for twisted bilayer graphene,” Phys. Rev. Lett. 123, 157601 (2019).
  * Da Liao  _et al._ (2021) Yuan Da Liao, Jian Kang, Clara N. Breiø, Xiao Yan Xu, Han-Qing Wu, Brian M. Andersen, Rafael M. Fernandes, and Zi Yang Meng, “Correlation-induced insulating topological phases at charge neutrality in twisted bilayer graphene,” Phys. Rev. X 11, 011014 (2021).
  * Classen  _et al._ (2019) Laura Classen, Carsten Honerkamp, and Michael M. Scherer, “Competing phases of interacting electrons on triangular lattices in moiré heterostructures,” Phys. Rev. B 99, 195120 (2019).
  * Kennes  _et al._ (2018) Dante M. Kennes, Johannes Lischner, and Christoph Karrasch, “Strong correlations and $d+\mathit{id}$ superconductivity in twisted bilayer graphene,” Phys. Rev. B 98, 241407 (2018).
  * Eugenio and Dağ (2020) P. Myles Eugenio and Ceren B. Dağ, “DMRG study of strongly interacting $\mathbb{Z}_{2}$ flatbands: a toy model inspired by twisted bilayer graphene,” SciPost Phys. Core 3, 15 (2020).
  * Huang  _et al._ (2020) Yixuan Huang, Pavan Hosur, and Hridis K. Pal, “Quasi-flat-band physics in a two-leg ladder model and its relation to magic-angle twisted bilayer graphene,” Phys. Rev. B 102, 155429 (2020).
  * Huang  _et al._ (2019) Tongyun Huang, Lufeng Zhang, and Tianxing Ma, “Antiferromagnetically ordered mott insulator and d+id superconductivity in twisted bilayer graphene: a quantum monte carlo study,” Science Bulletin 64, 310–314 (2019).
  * Guo  _et al._ (2018) Huaiming Guo, Xingchuan Zhu, Shiping Feng, and Richard T. Scalettar, “Pairing symmetry of interacting fermions on a twisted bilayer graphene superlattice,” Phys. Rev. B 97, 235453 (2018).
  * Ledwith  _et al._ (2020) Patrick J. Ledwith, Grigory Tarnopolsky, Eslam Khalaf, and Ashvin Vishwanath, “Fractional chern insulator states in twisted bilayer graphene: An analytical approach,” Phys. Rev. Research 2, 023237 (2020).
  * Repellin  _et al._ (2020) Cécile Repellin, Zhihuan Dong, Ya-Hui Zhang, and T. Senthil, “Ferromagnetism in narrow bands of moiré superlattices,” Phys. Rev. Lett. 124, 187601 (2020).
  * Abouelkomsan  _et al._ (2020) Ahmed Abouelkomsan, Zhao Liu, and Emil J. Bergholtz, “Particle-hole duality, emergent fermi liquids, and fractional chern insulators in moiré flatbands,” Phys. Rev. Lett. 124, 106803 (2020).
  * Repellin and Senthil (2020) Cécile Repellin and T. Senthil, “Chern bands of twisted bilayer graphene: Fractional chern insulators and spin phase transition,” Phys. Rev. Research 2, 023238 (2020).
  * Fernandes and Venderbos (2020) Rafael M. Fernandes and Jörn W. F. Venderbos, “Nematicity with a twist: Rotational symmetry breaking in a moiré superlattice,” Science Advances 6 (2020), 10.1126/sciadv.aba8834.
  * Wilson  _et al._ (2020) Justin H. Wilson, Yixing Fu, S. Das Sarma, and J. H. Pixley, “Disorder in twisted bilayer graphene,” Phys. Rev. Research 2, 023325 (2020).
  * Wang  _et al._ (2021) Jie Wang, Yunqin Zheng, Andrew J. Millis, and Jennifer Cano, “Chiral approximation to twisted bilayer graphene: Exact intravalley inversion symmetry, nodal structure, and implications for higher magic angles,” Phys. Rev. Research 3, 023155 (2021).
  * Bernevig  _et al._ (2021c) B. Andrei Bernevig, Zhi-Da Song, Nicolas Regnault, and Biao Lian, “Twisted bilayer graphene. I. Matrix elements, approximations, perturbation theory, and a $k\cdot p$ two-band model,” Physical Review B 103, 205411 (2021c), publisher: American Physical Society.
  * Song  _et al._ (2021) Zhi-Da Song, Biao Lian, Nicolas Regnault, and B. Andrei Bernevig, “Twisted bilayer graphene. ii. stable symmetry anomaly,” Phys. Rev. B 103, 205412 (2021).
  * Xie  _et al._ (2021) Fang Xie, Aditya Cowsik, Zhi-Da Song, Biao Lian, B. Andrei Bernevig, and Nicolas Regnault, “Twisted bilayer graphene. VI. An exact diagonalization study at nonzero integer filling,” Physical Review B 103, 205416 (2021), publisher: American Physical Society.
  * Zhang  _et al._ (2021) Xu Zhang, Gaopei Pan, Yi Zhang, Jian Kang, and Zi Yang Meng, “Momentum space quantum monte carlo on twisted bilayer graphene,” Chinese Physics Letters 38, 077305 (2021).
  * Cha  _et al._ (2021) Peter Cha, Aavishkar A Patel, and Eun-Ah Kim, “Strange metals from melting correlated insulators in twisted bilayer graphene,” arXiv preprint arXiv:2105.08069 (2021).
  * Chichinadze  _et al._ (2020) Dmitry V. Chichinadze, Laura Classen, and Andrey V. Chubukov, “Nematic superconductivity in twisted bilayer graphene,” Phys. Rev. B 101, 224513 (2020).
  * Sheffer and Stern (2021) Yarden Sheffer and Ady Stern, “Chiral magic-angle twisted bilayer graphene in a magnetic field: Landau level correspondence, exact wavefunctions and fractional chern insulators,” arXiv preprint arXiv:2106.10650 (2021).
  * Kang  _et al._ (2021) Jian Kang, B Andrei Bernevig, and Oskar Vafek, “Cascades between light and heavy fermions in the normal state of magic angle twisted bilayer graphene,” arXiv preprint arXiv:2104.01145 (2021).
  * Hofmann  _et al._ (2021) Johannes S Hofmann, Eslam Khalaf, Ashvin Vishwanath, Erez Berg, and Jong Yeon Lee, “Fermionic monte carlo study of a realistic model of twisted bilayer graphene,” arXiv preprint arXiv:2105.12112 (2021).
  * Calderón and Bascones (2020) M. J. Calderón and E. Bascones, “Interactions in the 8-orbital model for twisted bilayer graphene,” Phys. Rev. B 102, 155149 (2020).
  * Thomson and Alicea (2021) Alex Thomson and Jason Alicea, “Recovery of massless dirac fermions at charge neutrality in strongly interacting twisted bilayer graphene with disorder,” Phys. Rev. B 103, 125138 (2021).
  * Xie and MacDonald (2020b) Ming Xie and A. H. MacDonald, “Nature of the correlated insulator states in twisted bilayer graphene,” Phys. Rev. Lett. 124, 097601 (2020b).
  * (115) See supplementary material, which includes the references Bistritzer and MacDonald (2011); Song  _et al._ (2019); Tarnopolsky  _et al._ (2019); Song  _et al._ (2021); Bernevig  _et al._ (2021a); Elcoro  _et al._ (2020); Xu  _et al._ (2020); Marzari and Vanderbilt (1997); Souza  _et al._ (2001); Pizzi  _et al._ (2020); Winkler (2003); Wang  _et al._ (2021); Bultinck  _et al._ (2020a); Kang and Vafek (2019); Călugăru  _et al._ (2021); Lian  _et al._ (2021); Bernevig  _et al._ (2021b); Vafek and Kang (2021); Cea and Guinea (2020); Zhang  _et al._ (2020); Kang  _et al._ (2021); Xie  _et al._ (2019); Wong  _et al._ (2020); Ledwith  _et al._ (2020); Girvin  _et al._ (1986); Parameswaran  _et al._ (2013); Roy (2014); Sheffer and Stern (2021); Repellin and Senthil (2020); Abouelkomsan  _et al._ (2020).
  * Bradlyn  _et al._ (2017) Barry Bradlyn, L. Elcoro, Jennifer Cano, M. G. Vergniory, Zhijun Wang, C. Felser, M. I. Aroyo, and B. Andrei Bernevig, “Topological quantum chemistry,” Nature 547, 298–305 (2017).
  * Po  _et al._ (2018b) Hoi Chun Po, Haruki Watanabe, and Ashvin Vishwanath, “Fragile Topology and Wannier Obstructions,” Physical Review Letters 121, 126402 (2018b).
  * Ahn  _et al._ (2019) Junyeong Ahn, Sungjoon Park, and Bohm-Jung Yang, “Failure of Nielsen-Ninomiya Theorem and Fragile Topology in Two-Dimensional Systems with Space-Time Inversion Symmetry: Application to Twisted Bilayer Graphene at Magic Angle,” Physical Review X 9, 021013 (2019).
  * Bouhon  _et al._ (2019) Adrien Bouhon, Annica M. Black-Schaffer, and Robert-Jan Slager, “Wilson loop approach to fragile topology of split elementary band representations and topological crystalline insulators with time-reversal symmetry,” Phys. Rev. B 100, 195135 (2019).
  * Cano  _et al._ (2018) Jennifer Cano, Barry Bradlyn, Zhijun Wang, L. Elcoro, M. G. Vergniory, C. Felser, M. I. Aroyo, and B. Andrei Bernevig, “Topology of disconnected elementary band representations,” Phys. Rev. Lett. 120, 266401 (2018).
  * Marzari and Vanderbilt (1997) Nicola Marzari and David Vanderbilt, “Maximally localized generalized wannier functions for composite energy bands,” Phys. Rev. B 56, 12847–12865 (1997).
  * Souza  _et al._ (2001) Ivo Souza, Nicola Marzari, and David Vanderbilt, “Maximally localized wannier functions for entangled energy bands,” Phys. Rev. B 65, 035109 (2001).
  * Pizzi  _et al._ (2020) Giovanni Pizzi, Valerio Vitale, Ryotaro Arita, Stefan Blügel, Frank Freimuth, Guillaume Géranton, Marco Gibertini, Dominik Gresch, Charles Johnson, Takashi Koretsune, Julen Ibañez-Azpiroz, Hyungjun Lee, Jae-Mo Lihm, Daniel Marchand, Antimo Marrazzo, Yuriy Mokrousov, Jamal I Mustafa, Yoshiro Nohara, Yusuke Nomura, Lorenzo Paulatto, Samuel Poncé, Thomas Ponweiser, Junfeng Qiao, Florian Thöle, Stepan S Tsirkin, Małgorzata Wierzbowska, Nicola Marzari, David Vanderbilt, Ivo Souza, Arash A Mostofi, and Jonathan R Yates, “Wannier90 as a community code: new features and applications,” Journal of Physics: Condensed Matter 32, 165902 (2020).
  * Ramires and Lado (2021) Aline Ramires and Jose L. Lado, “Emulating heavy fermions in twisted trilayer graphene,” Phys. Rev. Lett. 127, 026401 (2021).
  * Si and Steglich (2010) Qimiao Si and Frank Steglich, “Heavy fermions and quantum phase transitions,” Science 329, 1161–1166 (2010).
  * Gegenwart  _et al._ (2008) Philipp Gegenwart, Qimiao Si, and Frank Steglich, “Quantum criticality in heavy-fermion metals,” nature physics 4, 186–197 (2008).
  * Coleman (1984) Piers Coleman, “New approach to the mixed-valence problem,” Phys. Rev. B 29, 3035–3044 (1984).
  * Dzero  _et al._ (2010) Maxim Dzero, Kai Sun, Victor Galitski, and Piers Coleman, “Topological kondo insulators,” Phys. Rev. Lett. 104, 106408 (2010).
  * Tsvelik and Wiegmann (1983) AM Tsvelik and PB Wiegmann, “Exact results in the theory of magnetic alloys,” Advances in Physics 32, 453–713 (1983).
  * Werner  _et al._ (2006) Philipp Werner, Armin Comanac, Luca de’ Medici, Matthias Troyer, and Andrew J. Millis, “Continuous-time solver for quantum impurity models,” Phys. Rev. Lett. 97, 076405 (2006).
  * Lu  _et al._ (2013) Feng Lu, JianZhou Zhao, Hongming Weng, Zhong Fang, and Xi Dai, “Correlated topological insulators with mixed valence,” Phys. Rev. Lett. 110, 096401 (2013).
  * Weng  _et al._ (2014) Hongming Weng, Jianzhou Zhao, Zhijun Wang, Zhong Fang, and Xi Dai, “Topological crystalline kondo insulator in mixed valence ytterbium borides,” Phys. Rev. Lett. 112, 016403 (2014).
  * Kotliar  _et al._ (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, “Electronic structure calculations with dynamical mean-field theory,” Rev. Mod. Phys. 78, 865–951 (2006).
  * Emery and Kivelson (1992) V. J. Emery and S. Kivelson, “Mapping of the two-channel kondo problem to a resonant-level model,” Phys. Rev. B 46, 10812–10817 (1992).
  * Freyer  _et al._ (2018) Frederic Freyer, Jan Attig, SungBin Lee, Arun Paramekanti, Simon Trebst, and Yong Baek Kim, “Two-stage multipolar ordering in $\mathrm{Pr}{T}_{2}{\mathrm{al}}_{20}$ kondo materials,” Phys. Rev. B 97, 115111 (2018).
  * Furusaki and Nagaosa (1994) Akira Furusaki and Naoto Nagaosa, “Kondo effect in a tomonaga-luttinger liquid,” Phys. Rev. Lett. 72, 892–895 (1994).
  * Cassanello and Fradkin (1996) Carlos R. Cassanello and Eduardo Fradkin, “Kondo effect in flux phases,” Phys. Rev. B 53, 15079–15094 (1996).
  * Martin  _et al._ (1997) Ivar Martin, Yi Wan, and Philip Phillips, “Size dependence in the disordered kondo problem,” Phys. Rev. Lett. 78, 114–117 (1997).
  * Chakravarty and Rudnick (1995) Sudip Chakravarty and Joseph Rudnick, “Dissipative dynamics of a two-state system, the kondo problem, and the inverse-square ising model,” Phys. Rev. Lett. 75, 501–504 (1995).
  * Elcoro  _et al._ (2020) Luis Elcoro, Benjamin J Wieder, Zhida Song, Yuanfeng Xu, Barry Bradlyn, and B Andrei Bernevig, “Magnetic topological quantum chemistry,” arXiv preprint arXiv:2010.00598 (2020).
  * Xu  _et al._ (2020) Yuanfeng Xu, Luis Elcoro, Zhi-Da Song, Benjamin J Wieder, MG Vergniory, Nicolas Regnault, Yulin Chen, Claudia Felser, and B Andrei Bernevig, “High-throughput calculations of magnetic topological materials,” Nature 586, 702–707 (2020).
  * Winkler (2003) Roland Winkler,  _Spin-orbit Coupling Effects in Two-Dimensional Electron and Hole Systems_ (Springer Science & Business Media, 2003).
  * Călugăru  _et al._ (2021) Dumitru Călugăru, Aaron Chew, Luis Elcoro, Nicolas Regnault, Zhi-Da Song, and B Andrei Bernevig, “General construction and topological classification of all magnetic and non-magnetic flat bands,” arXiv preprint arXiv:2106.05272 (2021).
  * Girvin  _et al._ (1986) S. M. Girvin, A. H. MacDonald, and P. M. Platzman, “Magneto-roton theory of collective excitations in the fractional quantum hall effect,” Phys. Rev. B 33, 2481–2494 (1986).
  * Parameswaran  _et al._ (2013) Siddharth A. Parameswaran, Rahul Roy, and Shivaji L. Sondhi, “Fractional quantum hall physics in topological flat bands,” Comptes Rendus Physique 14, 816–839 (2013), topological insulators / Isolants topologiques.
  * Roy (2014) Rahul Roy, “Band geometry of fractional topological insulators,” Phys. Rev. B 90, 165139 (2014).

Supplementary Materials

###### Contents

  1. S1 Bistritzer-MacDonald model
    1. S1.1 Basis and single-particle Hamiltonian
    2. S1.2 Discrete symmetries of the BM model
  2. S2 The single-particle Hamiltonian of the topological heavy fermion model
    1. S2.1 Maximal localized Wannier functions for the local orbitals
    2. S2.2 The topological conduction bands
    3. S2.3 Single-particle Hamiltonian of the topological conduction bands
    4. S2.4 Coupling between topological conduction bands and the local orbitals
    5. S2.5 Summary of the single-particle Hamiltonian
    6. S2.6 The (first) chiral limit and the chiral U(4) symmetry
    7. S2.7 The flat limit and the flat U(4) symmetry
    8. S2.8 Flat bands beyond the third chiral symmetry
    9. S2.9 The U(4)$\times$U(4) symmetry in chiral flat limit
    10. S2.10 Fubini-Study metric and Berry’s curvature of the flat bands in the chiral flat limit
  3. S3 The interaction Hamiltonian and its symmetries
    1. S3.1 The projected density operator and Coulomb interaction
    2. S3.2 Density-density terms (the first and fifth terms)
      1. S3.2.1 The density-density terms of local orbitals (the first term)
      2. S3.2.2 The density-density terms between local orbitals and conduction bands (the fifth term)
      3. S3.2.3 Density-density interaction of the conduction bands (the second term)
    3. S3.3 Exchange terms (the tenth term)
      1. S3.3.1 Constraints of the coupling matrix
      2. S3.3.2 Leading order coupling
      3. S3.3.3 Other coupling terms
    4. S3.4 Double hybridization terms (the third and fourth terms)
    5. S3.5 Density-hybridization terms (the sixth, seventh, eighth, and ninth terms)
      1. S3.5.1 Vanishing sixth and seventh terms
      2. S3.5.2 The eighth and ninth terms
    6. S3.6 Summary of the interaction Hamiltonian and its symmetries
      1. S3.6.1 Summary of the interaction Hamiltonian
      2. S3.6.2 Symmetries of the interaction Hamiltonian
    7. S3.7 Failure of the strong coupling picture in the first chiral limit of the BM model ($w_{0}=0$)
  4. S4 Correlated insulator states
    1. S4.1 Mean field Hamiltonian
    2. S4.2 Correlated insulator phases at the filling $\nu=0$
      1. S4.2.1 Numerical results
      2. S4.2.2 Analytical analysis of the one-shot HF
    3. S4.3 General rules for the ground states
    4. S4.4 Correlated insulator phases at the fillings $\nu=-1,-2,-3$
      1. S4.4.1 Filling $\nu=-1$
      2. S4.4.2 Filling $\nu=-2$
      3. S4.4.3 Filling $\nu=-3$
    5. S4.5 Chern numbers of the ground states
  5. S5 More discussions

##  S1 Bistritzer-MacDonald model

Here we briefly review the Bistritzer-MacDonald (BM) model Bistritzer and
MacDonald (2011) of magic-angle twisted bilayer graphene (MATBG). Readers may
refer to the supplementary materials of Ref. Song  _et al._ (2019) for more
details.

###  S1.1 Basis and single-particle Hamiltonian

As shown in Fig. S1(a), the low energy states around the Dirac points at the K
momenta of the two graphene layers form the moiré Brillouin zone (BZ) in the
valley K. Similarly, the low energy states around the Dirac points at the K′
momenta of the two layers form another moiré BZ in the valley K′. Thus in the
BM model there are two independent sets of basis from the two valleys. We use
the index $\eta$ ($=+$ for K and $-$ for K′) to label the two graphene
valleys. We denote the basis as $c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}$,
where where $\mathbf{k}$ is a momentum in the moiré Brillouin zone (BZ),
$\mathbf{Q}$ takes values in the lattice shown in Fig. S1(b), $\alpha=1,2$
represents the graphene sublattice, $\eta=\pm$ represents the graphene valley,
and $s=\uparrow\downarrow$ is the spin index. There are two types of Q
lattices: the blue lattice
$\mathcal{Q}_{+}=\\{\mathbf{q}_{2}+n_{1}\mathbf{b}_{M1}+n_{2}\mathbf{b}_{M2}\
|\ n_{1,2}\in\mathbb{Z}\\}$ and the red lattice
$\mathcal{Q}_{-}=\\{-\mathbf{q}_{2}+n_{1}\mathbf{b}_{M1}+n_{2}\mathbf{b}_{M2}\
|\ n_{1,2}\in\mathbb{Z}\\}$, where $\mathbf{b}_{M1},\mathbf{b}_{M2}$ are Moiré
reciprocal lattice basis

$\mathbf{b}_{M1}=\mathbf{q}_{2}-\mathbf{q}_{1},\qquad\mathbf{b}_{M_{2}}=\mathbf{q}_{3}-\mathbf{q}_{1}\
.$ (S1)

For $\mathbf{Q}\in\mathcal{Q}_{+}$, the basis is defined as

$(\mathbf{Q}\in\mathcal{Q}_{+})\qquad
c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}=\begin{cases}\displaystyle\frac{1}{\sqrt{N_{\rm
tot}}}\sum_{\mathbf{R}\in\mathrm{top}}e^{i(\mathbf{K}_{+}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}+\mathbf{t}_{\alpha})}c_{\mathbf{R},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=+\\\
\displaystyle\frac{1}{\sqrt{N_{\rm
tot}}}\sum_{\mathbf{R}^{\prime}\in\mathrm{bottom}}e^{i(-\mathbf{K}_{-}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\mathbf{t}_{\alpha}^{\prime})}c_{\mathbf{R}^{\prime},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=-\end{cases}\
,$ (S2)

and for $\mathbf{Q}\in\mathcal{Q}_{-}$, the basis is defined as

$(\mathbf{Q}\in\mathcal{Q}_{-})\qquad
c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}=\begin{cases}\displaystyle\frac{1}{\sqrt{N_{\rm
tot}}}\sum_{\mathbf{R}^{\prime}\in\mathrm{bottom}}e^{i(\mathbf{K}_{-}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\mathbf{t}_{\alpha}^{\prime})}c_{\mathbf{R}^{\prime},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=+\\\
\displaystyle\frac{1}{\sqrt{N_{\rm
tot}}}\sum_{\mathbf{R}\in\mathrm{top}}e^{i(-\mathbf{K}_{+}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}+\mathbf{t}_{\alpha})}c_{\mathbf{R},\alpha,s}^{\dagger}\qquad&\text{if}\;\eta=-\end{cases}\
.$ (S3)

Here $N_{\rm tot}$ is the number of graphene unit cells in each layer,
$\mathbf{R}$ and $\mathbf{R}^{\prime}$ indexes graphene lattices in the top
and bottom layers, respectively, $\mathbf{t}_{\alpha}$ and
$\mathbf{t}_{\alpha}^{\prime}$ are the sublattice vectors of the two layers,
respectively, $\mathbf{K}_{+}$ and $\mathbf{K}_{-}$ are the K momentum of the
top and bottom layers, respectively, and $c_{\mathbf{R},\alpha,s}$
($c_{\mathbf{R}^{\prime},\alpha,s}$) is the fermion operator with spin $s$ at
the atom site $\mathbf{R}+\mathbf{t}_{\alpha}$
($\mathbf{R}^{\prime}+\mathbf{t}_{\alpha}^{\prime}$). The $\mathcal{Q}_{+}$
($\mathcal{Q}_{-}$) lattice is defined in such a way that
$\eta\mathbf{K}_{+}+\mathbf{k}-\mathbf{Q}$
($\eta\mathbf{K}_{-}+\mathbf{k}-\mathbf{Q}$) with
$\mathbf{Q}\in\mathcal{Q}_{+}$ ($\mathbf{Q}\in\mathcal{Q}_{-}$) is the Dirac
point position $\eta\mathbf{K}_{+}$ ($\eta\mathbf{K}_{-}$) when $\mathbf{k}$
equals to the high symmetry point $\eta K_{M}$ ($\eta K_{M}^{\prime}$) of the
moiré BZ, as sown in Fig. S1(b).

Figure S1: The moiré BZ and the $\mathbf{Q}$ lattice. (a) The blue and red
hexagons represent the BZ’s of the top and bottom graphene layers,
respectively. The blue and red dots represent the Diract points (K momentum)
from the two layers. The vectors
$\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3}$ connecting the Dirac points of
the two layers span the moiré BZ in the valley K (shaded with light yellow).
Similarly, the Dirac points at the K′ momentum form another moiré BZ in the
valley K′ (shaded with light orange). (b) We denote the blue and red
$\mathbf{Q}$ lattices as $\mathcal{Q}_{+}$ and $\mathcal{Q}_{-}$ respectively.
In the BM model, the states in the valley K (K′) of the top layer and in the
valley K′ (K) of the bottom layer contribute to the lattice $\mathcal{Q}_{+}$
($\mathcal{Q}_{-}$). (c) The unit cell in real space, where the maximal
Wyckoff positions $1a$, $2c$, $3f$ are marked by the black dots.

The BM model is given by

${\hat{H}}_{\rm BM}=\sum_{\eta s}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\alpha\alpha^{\prime}}\sum_{\mathbf{Q},\mathbf{Q}^{\prime}}h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(\eta)}(\mathbf{k})c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}c_{\mathbf{k},\mathbf{Q}^{\prime},\alpha^{\prime},\eta,s}$
(S4)

where the single particle Hamiltonian reads

$h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(+)}(\mathbf{k})=v_{F}(\mathbf{k}-\mathbf{Q})\cdot\boldsymbol{\sigma}\delta_{\mathbf{Q},\mathbf{Q}^{\prime}}+\sum_{j=1}^{3}[T_{j}]_{\alpha\alpha^{\prime}}\delta_{\mathbf{Q},\mathbf{Q}^{\prime}\pm\mathbf{q}_{j}},\qquad
h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(-)}(\mathbf{k})=h_{-\mathbf{Q}\alpha,-\mathbf{Q}^{\prime}\alpha^{\prime}}^{(+)*}(-\mathbf{k}),$
(S5)
$T_{j}=w_{0}\sigma_{0}+w_{1}\sigma_{x}\cos\frac{2\pi(j-1)}{3}+w_{1}\sigma_{y}\sin\frac{2\pi(j-1)}{3}\
.$ (S6)

Here $\mathbf{q}_{j}$ are shown in Fig. S1, and $w_{0}$ and $w_{1}$ are the
interlayer couplings in the AA-stacking and AB-stacking regions, respectively.
The length of $\mathbf{q}_{j}$ is determined by the twist angle $\theta$,
i.e., $|\mathbf{q}_{j}|=k_{\theta}=2|\mathbf{K}|\sin\frac{\theta}{2}$. In this
work, we adopt the parameters $v_{F}=5.944\mathrm{eV\cdot\mathring{A}}$,
$|\mathbf{K}|=1.703\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$,
$\theta=1.05^{\circ}$. The relation
$h_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\alpha^{\prime}}^{(-)}(\mathbf{k})=h_{-\mathbf{Q}\alpha,-\mathbf{Q}^{\prime}\alpha^{\prime}}^{(+)*}(-\mathbf{k})$
is due to the time-reversal symmetry that transform the two valleys to each
other. The single-particle Hamiltonian (upon to a unitary transformation Song
_et al._ (2019)) is periodic with respect to the reciprocal lattice vectors
$\mathbf{b}_{M1}=\mathbf{q}_{2}-\mathbf{q}_{1}=k_{\theta}(\frac{\sqrt{3}}{2},\frac{3}{2})$,
$\mathbf{b}_{M2}=\mathbf{q}_{3}-\mathbf{q}_{1}=k_{\theta}(-\frac{\sqrt{3}}{2},\frac{3}{2})$.
The moiré BZ and high symmetry momenta are defined in Fig. S1(b). The real
space unit cell is generated by
$\mathbf{a}_{M1}=\frac{2\pi}{3k_{\theta}}(\sqrt{3},1)$,
$\mathbf{a}_{M2}=\frac{2\pi}{3k_{\theta}}(-\sqrt{3},1)$, which satisfy
$\mathbf{a}_{Mi}\cdot\mathbf{b}_{Mj}=2\pi\delta_{ij}$. Area of the moiré unit
cell is given by
$\Omega_{0}=\mathbf{a}_{M1}\times\mathbf{a}_{M2}=\frac{8\pi^{2}}{3\sqrt{3}k_{\theta}^{2}}$.
Maximal Wyckoff positions in the unit cell are shown in Fig. S1(c). The $1a$
and $2c$ positions correspond to the AA-stacking and AB-stacking regions,
respectively.

We plot the band structures in the valley $\eta=+$ for different $w_{0}$ in
Fig. S2.

The fact that the band structure is labeled by $\mathbf{k}$ (in the moiré BZ)
implies that $\mathbf{k}$ labels the eigenvalues of translation operators of
the moiré lattice. We hence define the translation operator $T_{\mathbf{R}}$
as

$T_{\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}T_{\mathbf{R}}^{-1}=e^{-i\mathbf{k}\cdot\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}\
,$ (S7)

where $\mathbf{R}=n_{1}\mathbf{a}_{M1}+n_{2}\mathbf{a}_{M2}$, with
$n_{1,2}\in\mathbb{Z}$, is a moiré lattice vector. One should not confuse
$T_{\mathbf{R}}$ with the time-reversal symmetry ($T$) defined in Section
S1.2. We now verify Eq. S7 at the commensurate twist angles, where
$\mathbf{a}_{M1}$ and $\mathbf{a}_{M2}$ are integer linear combinations of the
microscopic graphene lattice vectors $\mathbf{a}_{1},\mathbf{a}_{2}$.
$T_{\mathbf{R}}$ can be defined as translation acting on the atomic orbitals:
$\mathbf{R}^{\prime}\to\mathbf{R}^{\prime}+\mathbf{R}$. Then there is

$\displaystyle
T_{\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}T_{\mathbf{R}}^{-1}=\frac{1}{\sqrt{N_{\rm
tot}}}\sum_{\mathbf{R}^{\prime}\in
l}e^{i(\eta\mathbf{K}_{l}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\mathbf{t}_{\alpha})}c_{\mathbf{R}^{\prime}+\mathbf{R},\alpha,s}^{\dagger}=\frac{1}{\sqrt{N_{\rm
tot}}}\sum_{\mathbf{R}^{\prime}\in
l}e^{-i(\eta\mathbf{K}_{l}+\mathbf{k}-\mathbf{Q})\cdot\mathbf{R}}e^{i(\eta\mathbf{K}_{l}+\mathbf{k}-\mathbf{Q})\cdot(\mathbf{R}^{\prime}+\mathbf{t}_{\alpha})}c_{\mathbf{R}^{\prime},\alpha,s}^{\dagger}\
,$ (S8)

where $l=\eta$ ($-\eta$) for $\mathbf{Q}\in\mathcal{Q}_{+}$
($\mathcal{Q}_{-}$). For commensurate angle, $\eta\mathbf{K}_{l}-\mathbf{Q}$
is a moiré reciprocal lattice and hence
$e^{-i(\eta\mathbf{K}_{l}+\mathbf{k}-\mathbf{Q})\cdot\mathbf{R}}=e^{-i\mathbf{k}\cdot\mathbf{R}}$.
Then Eq. S7 is justified at commensurate angles. We emphasize that Eq. S7 is
also well-defined even at non-commensurate angles.

Figure S2: The band structure in the valley $\eta=+$ of MATBG. In (a)-(f) the
band structures with the inter-layer couplings $w_{0}/w_{1}$=1, 0.9, 0.8, 0.7,
0.6, 0.5 are plotted. The red crossings are energy eigenvalues of the BM
model, and the black lines are bands of the effective model. The other
parameters of the BM model are given by
$v_{F}=5.944\mathrm{eV\cdot\mathring{A}}$,
$|\mathbf{K}|=1.703\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$,
$\theta=1.05^{\circ}$. The parameters of the effective model for different
$w_{0}/w_{1}$ are given in Table S4. We have set
$v_{\star}^{\prime\prime}=B=C=0$ for simplicity.

For later convenience, here we define our convention for Fourier
transformation. Let $N$, $\Omega_{0}$, $\Omega_{\rm tot}=N\Omega_{0}$ be the
number of moiré unit cells, area of each moiré unit cell, and the total area
of the sample. Then, for $\mathbf{k},\mathbf{k}^{\prime}$ in the moiré BZ and
$\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta}$ with $l\eta=+$ or $-$,
integral over $\mathbf{r}$ leads to the momentum conservation

$\text{for}\quad\mathbf{Q},\mathbf{Q}^{\prime}\in\mathcal{Q}_{l\eta}\quad\int
d^{2}\mathbf{r}e^{(\mathbf{k}-\mathbf{k}^{\prime}-\mathbf{Q}+\mathbf{Q}^{\prime})\cdot\mathbf{r}}=\Omega_{\rm
tot}\delta_{\mathbf{k}\mathbf{k}^{\prime}}\delta_{\mathbf{Q},\mathbf{Q}^{\prime}}\
.$ (S9)

The completeness relation of the plane waves is

$\frac{1}{\Omega_{\rm
tot}}\sum_{\mathbf{k}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}=\delta(\mathbf{r}),\qquad\text{for}\;l\eta=\pm.$
(S10)

We introduce a continuous real space basis for the BM model

$c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})=\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}\
,$ (S11)

where $\Omega_{\rm tot}=N\Omega_{0}$ is the total area of the sample in
consideration, $N$ is the number of Moiré unit cells, $\Omega_{0}$ is the area
of the Moiré unit cell, $l=\pm$ is the layer index. We can understand the
continuous variable $\mathbf{r}$ as the coarse-grained graphene lattice
vector. The graphene sublattice ($\alpha$) and valley ($\eta$) are now
regarded as inner degrees of freedom. At each $\mathbf{r}$
$c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})$ has sixteen components from two
layers, two Graphene sublattice, two valleys, and two spins. Unlike the usual
real space basis, which just changes its position under a translation
operation, $c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})$ also gains a phase
factor under translation. According to Eq. S7, there is

$\displaystyle
T_{\mathbf{R}}c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})T_{\mathbf{R}}^{-1}=\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}T_{\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}T_{\mathbf{R}}^{-1}=\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}e^{-i\mathbf{k}\cdot\mathbf{R}}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}\
.$ (S12)

For convenience, we define the momenta

$\Delta\mathbf{K}_{l}=\begin{cases}\mathbf{q}_{2},\qquad&l=+\\\
-\mathbf{q}_{3},\qquad&l=-\end{cases}\ .$ (S13)

Notice that $\eta\Delta\mathbf{K}_{l}-\mathbf{Q}$
($\mathbf{Q}\in\mathcal{Q}_{\eta l}$) is a moiré reciprocal lattice and hence
$e^{-i(\eta\Delta\mathbf{K}_{l}-\mathbf{Q})\cdot\mathbf{R}}=1$. We hence have

$T_{\mathbf{R}}c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})T_{\mathbf{R}}^{-1}=\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot(\mathbf{r}+\mathbf{R})}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}=e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r}+\mathbf{R})\
.$ (S14)

One may remove the factor $e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}$ by
redefining the basis as $c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})\to
e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{r}}c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})$.
However, this redefining will complicate the representation of crystalline
symmetries (Section S1.2). In this work we will stick to the definition Eq.
S11 of the real space basis.

###  S1.2 Discrete symmetries of the BM model

The single-particle Hamiltonian in each valley has the symmetry of the
magnetic space group $P6^{\prime}2^{\prime}2$ (# 177.151 in the BNS setting),
which is generated by $C_{2z}T$, $C_{3z}$, $C_{2x}$ and translation
symmetries. Since the two valleys are related by time-reversal symmetry $T$,
the total system also has the $C_{2z}$ symmetry (product of $C_{2z}T$ and
$T$). The full crystalline symmetries of the two valleys form the space group
$P622$ (#177), which is generated by $C_{6z}$, $C_{2x}$, and translation
symmetries. We write the symmetry action on the fermion operators as

$\hat{g}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}\hat{g}^{-1}=\sum_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime}}c_{g\mathbf{k},\mathbf{Q}^{\prime},\alpha^{\prime},\eta^{\prime},s}^{\dagger}D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(g)\
.$ (S15)

The $D$ matrices of the time-reversal symmetry and the single-valley
symmetries are given by

$D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(T)=\delta_{\mathbf{Q}^{\prime},-\mathbf{Q}}[\sigma_{0}]_{\alpha^{\prime}\alpha}[\tau_{x}]_{\eta^{\prime}\eta},\qquad
D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(C_{3z})=\delta_{\mathbf{Q}^{\prime},C_{3z}\mathbf{Q}}[e^{i\frac{2\pi}{3}\sigma_{z}\tau_{z}}]_{\alpha^{\prime}\eta^{\prime},\alpha\eta}$
(S16)
$D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(C_{2x})=\delta_{\mathbf{Q}^{\prime},C_{2x}\mathbf{Q}}[\sigma_{x}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta},\qquad
D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(C_{2z}T)=\delta_{\mathbf{Q}^{\prime},\mathbf{Q}}[\sigma_{x}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta},$
(S17)

where $\sigma_{x,y,z}$ ($\sigma_{0}$), $\tau_{x,y,z}$ ($\tau_{0}$) are Pauli
(identity) matrices in the sublattice and valley spaces, and
$\zeta_{\mathbf{Q}}=\pm 1$ for $\mathbf{Q}\in\mathcal{Q}_{\pm}$.

Besides the crystalline symmetries, the BM model also has a unitary particle-
hole symmetry $P$ that anti-commutes with the single-particle Hamiltonian. $P$
transforms $\mathbf{k}$ to $-\mathbf{k}$ and the corresponding $D$ matrix
reads

$D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(P)=\zeta_{\mathbf{Q}}\delta_{\mathbf{Q}^{\prime},-\mathbf{Q}}[\sigma_{0}]_{\alpha^{\prime}\alpha}[\tau_{z}]_{\eta^{\prime}\eta}\
.$ (S18)

(Unlike the $P$ symmetry defined in Refs. Bernevig  _et al._ (2021a); Lian
_et al._ (2021); Bernevig  _et al._ (2021b), where $P$ is valley-independent,
here $P$ is chosen have opposite signs in the two valleys. The two definitions
are equivalent because the Hamiltonian respects the valley-U(1) and one can
redefine the particle-hole symmetry as $P\to\tau_{z}P$.) In the first-
quantized formalism the algebra between $P$ and other symmetries are given by

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

The particle-hole symmetry will be broken if the $\theta$-dependence of the
single-layer Hamiltonian or quadratic terms in $\mathbf{k}$ of the single-
layer Hamiltonian are taken into account.

If, furthermore, $w_{0}=0$, the model acquires an effective chiral symmetry
$C$ that anti-commute with the single-particle Hamiltonian Tarnopolsky  _et
al._ (2019). $C$ is referred to as the first chiral symmetry in Refs. Song
_et al._ (2021); Bernevig  _et al._ (2021a). $C$ leaves $\mathbf{k}$ invariant
and its $D$ matrix reads

$D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(C)=\delta_{\mathbf{Q}^{\prime},\mathbf{Q}}[\sigma_{z}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta}\
.$ (S20)

In the first-quantized formalism the algebra between $C$ and other symmetries
are given by

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

In the end we derive the symmetry actions on the continuous real space basis.
For unitary or anti-unitary operators, there are

$\hat{g}c_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})\hat{g}^{-1}=\sum_{l^{\prime}\alpha^{\prime}\eta^{\prime}}c_{l^{\prime},\alpha^{\prime},\eta^{\prime},s}^{\dagger}(g\mathbf{r})D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(g)\
.$ (S22)

Notice that $P\mathbf{r}=-\mathbf{r}$. ($P$ is an inversion-like operation
that interchanges the two layers but leave the sublattice invariant.) Using
the $D$ matrices in momentum space and the definition of
$c_{l,\alpha,\eta,s}(\mathbf{r})$ (Eq. S11), we can derive the $D$ matrices
for the real space continuous basis

$D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(T)=[\varrho_{0}]_{l^{\prime}l}[\sigma_{0}]_{\alpha^{\prime}\alpha}[\tau_{x}]_{\eta^{\prime}\eta},\qquad
D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(C_{3z})=[\varrho_{0}]_{l^{\prime}l}[e^{i\frac{2\pi}{3}\sigma_{z}\tau_{z}}]_{\alpha^{\prime}\eta^{\prime},\alpha\eta}\
,$ (S23)
$D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(C_{2x})=[\varrho_{x}]_{l^{\prime}l}[\sigma_{x}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta},\qquad
D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(C_{2z}T)=[\varrho_{0}]_{l^{\prime}l}[\sigma_{x}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta}\
,$ (S24)
$D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(P)=[-i\varrho_{y}]_{l^{\prime},l}[\sigma_{0}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta},\qquad
D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(C)=[\varrho_{0}]_{l^{\prime},l}[\sigma_{z}]_{\alpha^{\prime}\alpha}[\tau_{0}]_{\eta^{\prime}\eta}\
.$ (S25)

Here $\sigma_{x,y,z}$ ($\sigma_{0}$), $\tau_{x,y,z}$ ($\tau_{0}$),
$\varrho_{x,y,z}$ ($\varrho_{0}$) are Pauli (identity) matrices in the space
of sublattice, valley, and layer, respectively. We emphasize that, unlike the
momentum space basis where $D(P)$ is proportional to $\tau_{z}$, the real
space basis has $D$ proportional to $\tau_{0}$. Here we take the particle-hole
symmetry as an example to show how the $D$ matrices are derived. By
definition, there is

$\displaystyle
Pc_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})P^{-1}=\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}Pc_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}P^{-1}$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}c_{-\mathbf{k},-\mathbf{Q},\alpha,\eta,s}^{\dagger}\
\eta\zeta_{\mathbf{Q}}$ (S26)

We can replace $-\mathbf{Q}$ by $\mathbf{Q}^{\prime}$ with
$\mathbf{Q}^{\prime}\in\mathcal{Q}_{-l\eta}$. Relabeling $\mathbf{Q}^{\prime}$
and $\mathbf{k}$ as $\mathbf{Q}$ and $-\mathbf{k}$, respectively, we obtain

$\displaystyle
Pc_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})P^{-1}=\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{-l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot(-\mathbf{r})}c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}\eta\zeta_{-\mathbf{Q}}$
(S27)

Notice that, according to Eqs. S2 and S3, for
$\mathbf{Q}\in\mathcal{Q}_{-l\eta}$, the state
$c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}$ is in the layer $-l$, and
the factor $\eta\zeta_{-\mathbf{Q}}$ equals to $l$. We hence obtain

$Pc_{l,\alpha,\eta,s}^{\dagger}(\mathbf{r})P^{-1}=lc_{-l,\alpha,\eta,s}^{\dagger}(-\mathbf{r})\
,$ (S28)

consistent with Eq. S25.

For later convenience, we also derive the $D$ matrices of $C_{2z}P$ for the
momentum space basis and real space basis as

$D_{\mathbf{Q}^{\prime}\alpha^{\prime}\eta^{\prime},\mathbf{Q}\alpha\eta}(C_{2z}P)=-i\zeta_{\mathbf{Q}}\delta_{\mathbf{Q}^{\prime},\mathbf{Q}}[\sigma_{x}]_{\alpha^{\prime}\alpha}[\tau_{y}]_{\eta^{\prime}\eta},\qquad
D_{l^{\prime}\alpha^{\prime}\eta^{\prime},l\alpha\eta}(C_{2z}P)=[-i\varrho_{y}]_{l^{\prime},l}[\sigma_{x}]_{\alpha^{\prime}\alpha}[\tau_{x}]_{\eta^{\prime}\eta}\
,$ (S29)

respectively.

##  S2 The single-particle Hamiltonian of the topological heavy fermion model

###  S2.1 Maximal localized Wannier functions for the local orbitals

We first summarize the representation analysis of the Wannier obstruction
given in Ref. Song  _et al._ (2019). As proved in Ref. Song  _et al._ (2019),
due to the particle-hole symmetry, the middle two bands in each valley must
form the irreducible co-representations (dubbed as irreps)

$\Gamma_{1}\oplus\Gamma_{2};\ M_{1}\oplus M_{2};\ K_{2}K_{3}$ (S30)

of the magnetic space group $P6^{\prime}2^{\prime}2$ as long as they are
gapped from higher and lower energy bands. One can find the definitions of the
irreps in Table S1 or using the Corepresentations toolkit on the Bilbao
Crystallographic Server Elcoro  _et al._ (2020); Xu  _et al._ (2020). In Table
S2 we tabulate all the elementary band representations (EBRs), i.e., the
minimal local orbitals, and their corresponding irreps in momentum space. One
can also find the definition of EBRs using the MBANDREP toolkit on the Bilbao
Crystallographic Server Elcoro  _et al._ (2020); Xu  _et al._ (2020). There
are only six allowed two-band band representations, i.e.,
$2[A_{1}]_{a}\uparrow G$, $2[A_{2}]_{a}\uparrow G$, $[A_{1}\oplus
A_{2}]_{a}\uparrow G$, $[E]_{a}\uparrow G$, $[A_{1}]_{c}\uparrow G$, and
$[A_{2}]_{c}\uparrow G$, while none of them has the irreps in Eq. S30. (The
notation $[\rho]_{w}\uparrow G$ represents the EBR of the magnetic space group
$G$ induced from the representation $\rho$ at the Wyckoff possible $w$.) Thus
the middle two bands are not consistent with any local orbitals and must be
topological.

| $\Gamma_{1}$ | $\Gamma_{2}$ | $\Gamma_{3}$ |  | $M_{1}$ | $M_{2}$ |  | $K_{1}$ | $K_{2}K_{3}$
---|---|---|---|---|---|---|---|---|---
$E$ | 1 | 1 | 2 | $E$ | 1 | 1 | $E$ | 1 | 2
$2C_{3}$ | 1 | 1 | -1 | $C_{2}^{\prime}$ | 1 | -1 | $C_{3}$ | 1 | -1
$3C_{2}^{\prime}$ | 1 | -1 | 0 |  |  |  | $C_{3}^{-1}$ | 1 | -1
Table S1:  Character table of irreps at high symmetry momenta in magnetic space group $P6^{\prime}2^{\prime}2$ (#177.151 in BNS settings). For the little group of $\Gamma$, $E$, $2C_{3}$, and $3C_{2}^{\prime}$ represent the conjugation classes generated from identity, $C_{3z}$, and $C_{2x}$, respectively. Symbols for conjugate class at $M$ and $K$ are defined in similar ways.  Wyckoff pos. |  $1a$ $\left(000\right)$ |  $2c$ $\left(\frac{1}{3}\frac{2}{3}0\right)$, $\left(\frac{2}{3}\frac{1}{3}0\right)$ |  $3f$ $(\frac{1}{2}00)$, $(0\frac{1}{2}0)$, $(\frac{1}{2}\frac{1}{2}0)$
---|---|---|---
Site sym. |  $6^{\prime}22^{\prime}$, $32$ |  $32$, $32$ |  $2^{\prime}2^{\prime}2$, $2$
EBR | $[A_{1}]_{a}\uparrow G$ | $[A_{2}]_{a}\uparrow G$ | $[E]_{a}\uparrow G$ | $[A_{1}]_{c}\uparrow G$ | $[A_{2}]_{c}\uparrow G$ | $[E]_{c}\uparrow G$ | $[A]_{f}\uparrow G$ | $[B]_{f}\uparrow G$
Orbitals | $s$ | $p_{z}$ | $p_{x},p_{y}$ | $s$ | $p_{z}$ | $p_{x},p_{y}$ | $s$ | $p_{z}$
$\Gamma\left(000\right)$ | $\Gamma_{1}$ | $\Gamma_{2}$ | $\Gamma_{3}$ | $2\Gamma_{1}$ | $2\Gamma_{2}$ | $2\Gamma_{3}$ | $\Gamma_{1}\oplus\Gamma_{3}$ | $\Gamma_{2}\oplus\Gamma_{3}$
$K\left(\frac{1}{3}\frac{1}{3}0\right)$ | $K_{1}$ | $K_{1}$ | $K_{2}K_{3}$ | $K_{2}K_{3}$ | $K_{2}K_{3}$ | $2K_{1}\oplus K_{2}K_{3}$ | $K_{1}\oplus K_{2}K_{3}$ | $K_{1}\oplus K_{2}K_{3}$
$M\left(\frac{1}{2}00\right)$ | $M_{1}$ | $M_{2}$ | $M_{1}\oplus M_{2}$ | $2M_{1}$ | $2M_{2}$ | $2M_{1}\oplus 2M_{2}$ | $2M_{1}\oplus M_{2}$ | $M_{1}\oplus 2M_{2}$
Table S2:  EBRs of the magnetic space group $P6^{\prime}2^{\prime}2$ (#177.151
in the BNS setting). In the first row are the symbols and coordinates of
Wyckoff positions. In the second row are the corresponding site symmetry
groups (magnetic point groups) and their unitary subgroups (point groups). In
the third row are the names of EBRs. The notation $[\rho]_{w}\uparrow G$
represents the EBR of the magnetic space group $G$ induced from the irrep
$\rho$ at the Wyckoff possible $w$. The orbitals that form the representation
$[\rho]_{w}$ are given in the fourth row. The fifth row to seventh row give
the irreps at high symmetry momenta of each EBR.

We now resolve the Wannier obstruction by adding higher energy bands in MATBG.
Both the first “passive” bands above the middle two bands and the first
“passive” bands below the middle two bands form the degenerate irrep
$\Gamma_{3}$ at the $\Gamma_{M}$ point. If we hybridize the middle two bands
with these higher energy bands, we can extract two topologically trivial bands
that form the band representation (in terms of irreps)

$[E]_{a}\uparrow G:\qquad\Gamma_{3};\ M_{1}\oplus M_{2};\ K_{2}K_{3}$ (S31)

where $[E]_{a}$ is formed by a $p_{x}$-like orbital and a $p_{y}$-like orbital
at the $1a$ position (triangular lattice). In the remaining of this subsection
we will discuss an explicit construction of the Wannier functions forming the
representation $[E]_{a}$. It has large overlaps (Fig. S3) with the lowest two
bands for $\mathbf{k}$ away from $\Gamma_{M}$ — in most area of the moiré BZ
excluding a small neighborhood of $\Gamma_{M}$, the lowest two bands of MATBG
are almost completely spanned by the constructed Wannier functions. While at
(or around) $\Gamma_{M}$, the Wannier functions have zero (or small) overlap
with the lowest two bands of MATBG. Another set of basis, which form the
representation $\Gamma_{1}\oplus\Gamma_{2}\oplus\Gamma_{3}$ at $\Gamma_{M}$,
will be introduced around $\Gamma_{M}$. In Sections S2.2, S2.3, S2.4 and S2.5
we will write the single-particle Hamiltonian on the basis of the constructed
Wannier functions and the remaining low-energy bands, which together form the
representation $\Gamma_{1}\oplus\Gamma_{2}\oplus 2\Gamma_{3}$ at the
$\Gamma_{M}$ point.

For simplicity in this subsection we will use the first-quantized formalism.
We define the basis

$|\mathbf{k},\mathbf{Q},\alpha,\eta,s\rangle=c_{\mathbf{k},\mathbf{Q},\alpha,\eta,s}^{\dagger}|0\rangle\
,$ (S32) $|\mathbf{r},l,\alpha,\eta,s\rangle=c_{l\alpha\eta
s}^{\dagger}(\mathbf{r})|0\rangle=\frac{1}{\sqrt{\Omega_{\rm
tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}|\mathbf{k},\mathbf{Q},\alpha,\eta,s\rangle\
.$ (S33)

We write the trial Wannier functions as

$|W^{\prime}_{\mathbf{R},\alpha=1,\eta,s}\rangle=\frac{1}{\lambda_{0}\sqrt{2\pi}}\sum_{l}\int
d^{2}\mathbf{r}\
e^{i\frac{\pi}{4}l\eta-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}-(\mathbf{r}-\mathbf{R})^{2}/(2\lambda_{0}^{2})}|\mathbf{r},l,1,\eta,s\rangle$
(S34)
$|W^{\prime}_{\mathbf{R},\alpha=2,\eta,s}\rangle=\frac{1}{\lambda_{0}\sqrt{2\pi}}\sum_{l}\int
d^{2}\mathbf{r}\
e^{-i\frac{\pi}{4}l\eta-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}-(\mathbf{r}-\mathbf{R})^{2}/(2\lambda_{0}^{2})}|\mathbf{r},l,2,\eta,s\rangle\
,$ (S35)

where $\lambda_{0}$ represents the size of the Gaussian function and
$\mathbf{R}$ takes value in the triangular lattice sites, such that they form
the representations

$\displaystyle D^{f}(T)=\sigma_{0}\tau_{x},\quad
D^{f}(C_{3z})=e^{i\frac{2\pi}{3}\sigma_{z}\tau_{z}},\quad
D^{f}(C_{2x})=\sigma_{x}\tau_{0},\quad D^{f}(C_{2z}T)=\sigma_{x}\tau_{0},$
$\displaystyle D^{f}(P)=i\sigma_{z}\tau_{z},\quad
D^{f}(C_{2z}P)=-i\sigma_{y}\tau_{y}\ ,$ (S36)

same as $p_{x}\pm ip_{y}$ at the $1a$ Wyckoff position (per spin per valley).
Here $\sigma_{0,x,y,z}$, $\tau_{0,x,y,z}$ are the identity or Pauli matrices
in for the orbital ($\alpha=1,2$) and valley ($\eta=\pm$) degrees of freedom,
respectively. Since the Wannier functions transform diagonally (i.e.,
$D^{f}(P)$ in Eq. S36) under the $P$ symmetry and $P$ flips the layer (Eq.
S25), each Wannier function should be a linear combination of states from the
two layers. The phase factors $e^{\pm i\frac{\pi}{4}l\eta}$ in the Wannier
functions are introduced to fix the gauge of the representation matrices
$D^{f}$’s. The phase factor $e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}$ is
required by the translation symmetry (Eq. S14). One can verify that such
defined trial Wannier functions satisfy

$T_{\Delta\mathbf{R}}|W^{\prime}_{\mathbf{R},\alpha,\eta,s}\rangle=|W^{\prime}_{\mathbf{R}+\Delta\mathbf{R},\alpha,\eta,s}\rangle\
,$ (S37)

under a translation of the moiré lattice $\Delta\mathbf{R}$. We emphasize
that, for given $D^{f}(C_{3z})$, the sign of $D^{f}(P)$ is not a gauge choice.
Since $[C_{3z},P]=0$, $P$ and $C_{3z}$ can be diagonalized at the same time
and the $P$-eigenvalues in each $C_{3z}$ eigen-space are determined. For the
trial Wannier functions, the $P$ eigenvalue of the $C_{3z}$ eigenstates with
the eigenvalue $e^{i\frac{2\pi}{3}}$ ($e^{-i\frac{2\pi}{3}}$) is $i$ ($-i$).

We can express the trial Wannier functions on the plane-wave basis as

$\displaystyle|W^{\prime}_{\mathbf{R},\alpha=1,\eta,s}\rangle=$
$\displaystyle\frac{1}{\lambda_{0}\sqrt{2\pi\Omega_{\rm
tot}}}\sum_{l=\pm}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}\int d^{2}\mathbf{r}\
e^{i\frac{\pi}{4}l\eta-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}-(\mathbf{r}-\mathbf{R})^{2}/(2\lambda_{0}^{2})-i(\mathbf{k}-\mathbf{Q})\cdot\mathbf{r}}|\mathbf{k},\mathbf{Q},1,\eta,s\rangle$
$\displaystyle=$ $\displaystyle\sqrt{\frac{2\pi\lambda_{0}^{2}}{\Omega_{\rm
tot}}}\sum_{l=\pm}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{i\frac{\pi}{4}l\eta-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}-\mathbf{Q})^{2}}|\mathbf{k},\mathbf{Q},1,\eta,s\rangle$
(S38)
$|W^{\prime}_{\mathbf{R},\alpha=2,\eta,s}\rangle=\sqrt{\frac{2\pi\lambda_{0}^{2}}{\Omega_{\rm
tot}}}\sum_{l=\pm}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i\frac{\pi}{4}l\eta-i\mathbf{k}\cdot\mathbf{R}-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}-\mathbf{Q})^{2}}|\mathbf{k},\mathbf{Q},2,\eta,s\rangle\
.$ (S39)

We denote the energy eigenstates of the continuous BM model as
$|\psi_{\mathbf{k},n,\eta,s}\rangle$. The overlap matrix between the energy
eigenstates and the trial Wannier functions are

$A_{n,1}^{(\eta)}(\mathbf{k})=\langle\psi_{\mathbf{k},n,\eta,s}|W^{\prime}_{0,1,\eta,s}\rangle=\sqrt{\frac{2\pi\lambda_{0}^{2}}{\Omega_{\rm
tot}}}\sum_{l=\pm}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{i\frac{\pi}{4}l\eta-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}-\mathbf{Q})^{2}}\langle\psi_{\mathbf{k},n,\eta,s}|\mathbf{k},\mathbf{Q},1,\eta,s\rangle\
,$ (S40)
$A_{n,2}^{(\eta)}(\mathbf{k})=\langle\psi_{\mathbf{k},n,\eta,s}|W^{\prime}_{0,2,\eta,s}\rangle=\sqrt{\frac{2\pi\lambda_{0}^{2}}{\Omega_{\rm
tot}}}\sum_{l=\pm}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{-i\frac{\pi}{4}l\eta-\frac{1}{2}\lambda_{0}^{2}(\mathbf{k}-\mathbf{Q})^{2}}\langle\psi_{\mathbf{k},n,\eta,s}|\mathbf{k},\mathbf{Q},2,\eta,s\rangle\
.$ (S41)

We plot the $\sum_{\alpha}|A_{n,\alpha}(\mathbf{k})|^{2}$ for each band
$\epsilon_{\mathbf{k},n}$ ($n=\pm 1,\pm 2,\pm 3$) in Fig. S3(a). One can see
that, as expected, at the $\Gamma_{M}$ point, the overlap between the trial
Wannier functions and the $\Gamma_{1}$, $\Gamma_{2}$ states in the middle two
bands are zero.

Figure S3: Wannier functions for $w_{0}/w_{1}=0.8$. In (a) the sizes of the
red dots represent the overlaps between the trial Wannier functions and the
Bloch states, i.e., $\sum_{\alpha}|A_{n,\alpha}(\mathbf{k})|^{2}$, for the
bands $\epsilon_{\mathbf{k},n}$ with $n=\pm 1,\pm 2,\pm 3$. (b) and (c) are
absolute values of the two components $w_{+\beta,1}(\mathbf{r})$ ($\beta=1,2$)
of the constructed maximally localized Wannier functions. The hexagons
represent the moiré unit cell. (d) The density of Wannier functions i.e.,
$\sum_{l\beta}|w^{(\eta)}_{l\beta,\alpha}(\mathbf{r})|^{2}$, which does not
depend on $\eta,\alpha$. (e), (f), (g) are $|w_{+1,1}(\mathbf{r})|$,
$|w_{+2,1}(\mathbf{r})|$, and the density of the approximate Wannier functions
of the Gaussian form (Eqs. S53 and S54), respectively. The used parameters for
the BM model are $v_{F}=5.944\mathrm{eV\cdot\mathring{A}}$,
$|\Delta\mathbf{K}|=1.703\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$,
$\theta=1.05^{\circ}$.

We feed the overlap matrix $A_{n,\alpha}^{(\eta)}(\mathbf{k})$ ($n=\pm 1,\pm
2,\pm 3$, $\alpha=1,2$, $\eta=\pm$) into the machinery of Wannier90 Marzari
and Vanderbilt (1997); Souza  _et al._ (2001); Pizzi  _et al._ (2020) to
construct the maximally localized Wannier functions. We have set
$\lambda_{0}=0.1a_{M}$ in practical calculations, with $a_{M}$ being the moiré
lattice constant. We use a $18\times 18$ momentum mesh and the energy window
$\rm[-80meV,80meV]$ for the disentanglement and Wannierization procedures.
Wannier90 returns $|W_{\mathbf{R},\alpha,\eta,s}\rangle$ on the plane-wave
basis $|\mathbf{k},\mathbf{Q},\beta,\eta,s\rangle$ as

$|W_{\mathbf{R},\alpha,\eta,s}\rangle=\frac{1}{\sqrt{N}}\sum_{l=\pm}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\beta}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}|\mathbf{k},\mathbf{Q},\beta,\eta,s\rangle
e^{-i\mathbf{k}\cdot\mathbf{R}}\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})$
(S42)

Here $N$ is the number of Moiré unit cells.
$\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})$ returned by
Wannier90 is originally stored on the basis of the Bloch wave functions of the
energy bands. Only the lowest six energy bands are involved in the Wannier
functions because we have chosen the energy window $\rm[-80meV,80meV]$ for the
Wannierization, which only covers the lowest six bands (Fig. S3(a)). We then
transform $\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})$ onto
the basis $|\mathbf{k},\mathbf{Q},\beta,\eta,s\rangle$. Since
$|\mathbf{r},l,\beta,\eta,s\rangle$ gains a phase factor
$e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}$ under the translation
$\mathbf{R}$ (Eq. S14), we can write
$\langle\mathbf{r},l,\beta,\eta,s|W_{\mathbf{R}\alpha\eta s}\rangle$ as

$\langle\mathbf{r},l,\beta,\eta,s|W_{\mathbf{R}\alpha\eta
s}\rangle=\langle\mathbf{r},l,\beta,\eta,s|T_{\mathbf{R}}|W_{0\alpha\eta
s}\rangle=e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}\langle\mathbf{r}-\mathbf{R},l,\beta,\eta,s|W_{0\alpha\eta
s}\rangle=e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}w_{l\beta,\alpha}^{(\eta)}(\mathbf{r}-\mathbf{R})$
(S43)

such that
$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r}-\mathbf{R})=\langle\mathbf{r}-\mathbf{R},l,\beta,\eta,s|W_{0\alpha\eta
s}\rangle$ is a function of $\mathbf{r}-\mathbf{R}$. The real space shape of
the Wannier states can be calculated as

$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r}-\mathbf{R})=e^{i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}\langle\mathbf{r},l,\beta,\eta,s|W_{\mathbf{R}\alpha\eta
s}\rangle=\frac{1}{\sqrt{N\Omega_{\rm tot}}}\sum_{\mathbf{k}\in{\rm
MBZ}}\sum_{\mathbf{Q}\in\mathcal{Q}_{l\eta}}e^{i(\mathbf{k}-\mathbf{Q})\cdot(\mathbf{r}-\mathbf{R})}\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})\
.$ (S44)

Making use of Eqs. S9 and S10, we obtain the inverse transformation

$\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})=\frac{1}{\sqrt{\Omega_{0}}}\int
d^{2}\mathbf{r}\
w_{l_{\mathbf{Q},\eta}\beta,\alpha}^{(\eta)}(\mathbf{r}-\mathbf{R})e^{-i(\mathbf{k}-\mathbf{Q})\cdot(\mathbf{r}-\mathbf{R})}\
,$ (S45)

where $l_{\mathbf{Q},\eta}=\zeta_{\mathbf{Q}}\eta$ is the layer the basis
$|\mathbf{k},\mathbf{Q},\beta,\eta,s\rangle$ belonging to and
$\Omega_{0}=\Omega_{\rm tot}/N$ is the area of the moiré unit cell.

We now determine the independent components of the Wannier function. Due to
the spin-SU(2) symmetry, the real space wave functions
$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})$ do not depend on spin. Due to the
time-reversal symmetry, the $w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})$ in the
two valleys are just complex conjugation of each other, i.e.,

$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})=w_{l\beta,\alpha}^{(-\eta)*}(\mathbf{r})\
.$ (S46)

Due to the $P$ symmetry ($\mathbf{r}\to-\mathbf{r}$), which has
$D(P)=-i\varrho_{y}\sigma_{0}\tau_{0}$ for the continuous real space basis
(Eq. S25) and $D^{f}(P)=i\sigma_{z}\tau_{z}$ for the Wannier functions (Eq.
S36), the Wannier functions must satisfy the constraint
$\sum_{l^{\prime}\alpha^{\prime}\eta^{\prime}}D_{l\alpha\eta,l^{\prime}\alpha^{\prime}\eta}(P)w_{l^{\prime}\alpha^{\prime}\eta^{\prime},\alpha}^{(\eta)}(\mathbf{r})=\sum_{\alpha^{\prime}}w_{l\alpha\eta,\alpha^{\prime}}^{(\eta)}(-\mathbf{r})D^{f}_{\eta\alpha^{\prime},\eta\alpha}(P)$,
i.e.,

$w_{-l\beta,\alpha}^{(\eta)}(-\mathbf{r})=il\eta(-1)^{\alpha}w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})$
(S47)

Due to the $C_{2z}T$ symmetry, which has
$D(C_{2z}T)=\varrho_{0}\sigma_{x}\tau_{0}$ for the continuous basis (Eq. S24)
and $D^{f}(C_{2z}T)=\sigma_{x}\tau_{0}$ for the Wannier functions (Eq. S36),
the Wannier functions must satisfy the constraint

$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})=w_{l\overline{\beta},\overline{\alpha}}^{(\eta)*}(-\mathbf{r})$
(S48)

where $\overline{\alpha}=2,1$ for $\alpha=1,2$. The above two conditions
together imply a constraint (given by $PC_{2z}T$) that holds at each point
$\mathbf{r}$

$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})=-il\eta(-1)^{\alpha}w_{-l\overline{\beta},\overline{\alpha}}^{(\eta)*}(\mathbf{r})\
.$ (S49)

Due to the above three independent constraints, only two of the sixteen
components in the real space functions
$w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})$ are independent. (Notice that each
constraint reduce to half the number independent components.) Without loss of
generality, we will focus on the two independent components
$w_{+\beta,1}^{(+)}(\mathbf{r})$ ($\beta=1,2$). We plot the absolute values of
the two components $w_{+\beta,1}^{(+)}(\mathbf{r})$ ($\beta=1,2$) in Fig.
S3(b), (c), respectively. Due to Eqs. S49 and S46, the densities of the
Wannier functions, i.e.,
$\sum_{l\beta}|w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})|^{2}$, do not depend on
$\eta$ and $\alpha$. We plot the density of the Wannier functions in Fig.
S3(d). From the density plot we can see that the Wannier functions are very
well localized. To characterize the localization, we also calculate the square
root of the spread of the Wannier functions, i.e.,

$\lambda=\sqrt{\sum_{l\beta}\int
d^{2}\mathbf{r}|w_{l\beta,\alpha}^{(\eta)}(\mathbf{r})|^{2}\mathbf{r}^{2}}\ .$
(S50)

$\lambda$ does not depend on $\alpha$ and $\eta$ according to Eqs. S49 and
S46. As tabulated in Table S3, $\lambda$ ranges from $0.2965a_{M}$ to
$0.5330a_{M}$ for $w_{0}/w_{1}$ changing from 1.0 to 0.5, where $a_{M}$ is the
lattice constant of moiré unit cell. The magnitude of the hopping between the
constructed Wannier functions at nearest neighbor sites, i.e., $t_{0}=|\langle
W_{\mathbf{R},\alpha,\eta,s}|H_{\rm
BM}|W_{\mathbf{R}^{\prime},-\alpha,\eta,s}\rangle|$ for
$\mathbf{R}\mathbf{R}^{\prime}$ being nearest neighbor pairs, is also
extremely small. ($\langle W_{\mathbf{R},\alpha,\eta,s}|H_{\rm
BM}|W_{\mathbf{R}^{\prime},\alpha,\eta,s}\rangle=0$ due to the crystalline and
$P$ symmetries.) As shown in Table S4, the $t_{0}$ ranges from 0.01512meV to
0.3027meV for $w_{0}/w_{1}$ changing from 1.0 to 0.5. In the rest of this work
we will omit the nearest neighbor hoppings.

$w_{0}/w_{1}$ | $\alpha_{1}$ | $\alpha_{2}$ |  $\lambda_{1}$ ($a_{M}$) |  $\lambda_{2}$ ($a_{M}$) |  $\lambda$ ($a_{M}$)
---|---|---|---|---|---
1.0 | 0.7632 | 0.6461 | 0.1552 | 0.1731 | 0.2965
0.9 | 0.7882 | 0.6155 | 0.1671 | 0.1850 | 0.3029
0.8 | 0.8193 | 0.5734 | 0.1791 | 0.1910 | 0.3375
0.7 | 0.8531 | 0.5217 | 0.1970 | 0.1970 | 0.3792
0.6 | 0.8878 | 0.4603 | 0.2089 | 0.2029 | 0.4380
0.5 | 0.9208 | 0.3901 | 0.2268 | 0.2089 | 0.5330
Table S3: Parameters of Wannier functions for different $w_{0}/w_{1}$.
$\alpha_{1,2}$ and $\lambda_{1,2}$ are parameters for the analytical form
(Eqs. S53 and S54) of the Wannier functions and they are fitted from the
numerical Wannier functions. $\lambda$ is the square root of the spread of the
numerical Wannier functions. The used parameters for the BM model are
$v_{F}=5.944\mathrm{eV\cdot\mathring{A}}$,
$|\Delta\mathbf{K}|=1.703\mathrm{\mathring{A}^{-1}}$, $w_{1}=110\mathrm{meV}$,
$\theta=1.05^{\circ}$.

We find the Wannier functions shown in Fig. S3(b), (c) can be approximated by
Gaussian functions of the form

$w_{+1,1}^{(+)}(\mathbf{r})=\frac{\alpha_{1}}{\sqrt{2}}\frac{1}{\sqrt{\pi\lambda_{1}^{2}}}e^{i\frac{\pi}{4}-\mathbf{r}^{2}/(2\lambda_{1}^{2})},\quad
w_{+2,1}^{(+)}(\mathbf{r})=-\frac{\alpha_{2}}{\sqrt{2}}\frac{x+iy}{\lambda_{2}^{2}\sqrt{\pi}}e^{i\frac{\pi}{4}-\mathbf{r}^{2}/(2\lambda_{2}^{2})}\
.$ (S51)

Even though the Wannier states transform as $p_{x}\pm ip_{y}$ orbitals, which
form the representation $[E]_{a}$ at the $1a$ position (Tables S1 and S2), we
see that the component $w_{+1,1}^{(+)}(\mathbf{r})$ is $s$-orbital-like and
the component $w_{+2,1}^{(+)}(\mathbf{r})$ is ($p_{x}+ip_{y}$)-orbital-like.
The reason is that the inner degrees of freedom, i.e., layer, sublattice, and
valley, transform non-trivially under the crystalline symmetries. For example,
the layer, sublattice, and valley degrees of freedom transform under $C_{3z}$
as $D(C_{3z})=\varrho_{0}e^{i\frac{2\pi}{3}\sigma_{z}\tau_{z}}$ (Eq. S23), in
order for the Wannier functions $w_{+\beta,1}^{(+)}(\mathbf{r})$ to have the
$C_{3z}$ eigenvalue $e^{i\frac{2\pi}{3}}$ (Eq. S36), there must be
$e^{i\frac{2\pi}{3}[\sigma_{z}]_{\beta\beta}}w_{+\beta,1}^{(+)}(C_{3z}^{-1}\mathbf{r})=e^{i\frac{2\pi}{3}}w_{+\beta,1}^{(+)}(\mathbf{r})$.
The two components satisfy
$w_{+1,1}^{(+)}(C_{3z}\mathbf{r})=w_{+1,1}^{(+)}(\mathbf{r})$ and
$w_{+2,1}^{(+)}(C_{3z}\mathbf{r})=e^{i\frac{2\pi}{3}}w_{+2,1}^{(+)}(\mathbf{r})$,
respectively. Thus they can be realized by the polynomials $1$ (constant) and
$x+iy$, respectively. The phase factor $e^{i\frac{\pi}{4}}$ is required by the
symmetry $C_{2x}PC_{2z}T=C_{2y}PT$, whose action action on the continuous
basis and representation on the Wannier states are
$D(C_{2y}PT)=\varrho_{z}\sigma_{0}\tau_{0}$ (Eqs. S24 and S25) and
$D^{f}(C_{2y}PT)=-i\sigma_{z}\tau_{z}$ (Eq. S36), respectively. Given
$w_{+1,1}^{(+)}(\mathbf{r})\propto e^{-\mathbf{r}^{2}/(2\lambda_{1}^{2})}$,
$w_{+1,1}^{(+)}(\mathbf{r})\propto(x+iy)e^{-\mathbf{r}^{2}/(2\lambda_{1}^{2})}$,
the phases of the two components are constrained by the $C_{2y}PT$ symmetry as

$w_{+\beta,1}^{(+)*}(x,-y)=-iw_{+\beta,1}^{(+)}(x,y)\qquad\Rightarrow\qquad
w_{+\beta,1}^{(+)}(\mathbf{r})=\pm
e^{i\frac{\pi}{4}}|w_{+\beta,1}^{(+)}(\mathbf{r})|\ .$ (S52)

However, the symmetry representation cannot fix the relative sign of the two
components: Suppose $w_{+\beta,1}^{(+)}(\mathbf{r})$ forms the correct
representation of $C_{2y}PT$, then
$(-1)^{\beta-1}w_{+\beta,1}^{(+)}(\mathbf{r})$ also forms the correct
representation. Thus there are two symmetry-allowed choices for the analytical
expressions of $w_{+\beta,1}^{(+)}(\mathbf{r})$. We have computed the overlaps
between the two choices and the numerical Wannier functions and pick the
choice (Eq. S51) that gives the larger overlap. $\lambda_{1}$ can be fitted
using the condition
$|w_{+1,1}^{(+)}(\lambda_{1},0)|=e^{-\frac{1}{2}}|w_{+1,1}^{(+)}(0)|$, and
$\lambda_{2}$ can be fitted by the maximum of $|w_{+2,1}^{(+)}(x_{\max},0)|$,
where $x_{\max}=\sqrt{2}\lambda_{2}$. $\alpha_{1}$ and $\alpha_{2}$ are real
coefficients satisfying $\alpha_{1}^{2}+\alpha_{2}^{2}=1$. For
$w_{0}/w_{1}=0.8$, we find $\lambda_{1}=0.1791a_{M}$,
$\lambda_{2}=0.1850a_{M}$, $\alpha_{1}=0.7882$, $\alpha_{2}=0.6155$. These
parameters for other $w_{0}/w_{1}$ values are tabulated in Table S3. Using the
constraints Eqs. S46, S47 and S48, we obtain the other components as

$w_{l1,1}^{(\eta)}(\mathbf{r})=\frac{\alpha_{1}}{\sqrt{2}}\frac{1}{\sqrt{\pi\lambda_{1}^{2}}}e^{i\frac{\pi}{4}l\eta-\mathbf{r}^{2}/(2\lambda_{1}^{2})},\quad
w_{l2,1}^{(\eta)}(\mathbf{r})=-l\frac{\alpha_{2}}{\sqrt{2}}\frac{x+i\eta
y}{\lambda_{2}^{2}\sqrt{\pi}}e^{i\frac{\pi}{4}l\eta-\mathbf{r}^{2}/(2\lambda_{2}^{2})},$
(S53) $w_{l1,2}^{(\eta)}(\mathbf{r})=l\frac{\alpha_{2}}{\sqrt{2}}\frac{x-i\eta
y}{\lambda_{2}^{2}\sqrt{\pi}}e^{-i\frac{\pi}{4}l\eta-\mathbf{r}^{2}/(2\lambda_{2}^{2})},\qquad
w_{l2,2}^{(\eta)}(\mathbf{r})=\frac{\alpha_{1}}{\sqrt{2}}\frac{1}{\sqrt{\pi\lambda_{1}^{2}}}e^{-i\frac{\pi}{4}l\eta-\mathbf{r}^{2}/(2\lambda_{1}^{2})}\
.$ (S54)

One can verify that the Wannier functions form the representations in Eq. S36.

We introduce the creation operators for the Wannier states as

$f_{\mathbf{R}\alpha\eta s}^{\dagger}=\sum_{l\beta}\int d^{2}\mathbf{r}\
\langle\mathbf{r},l,\beta,s,\eta,s|W_{\mathbf{R},\alpha,\eta,s}\rangle
c_{l,\beta,\eta,s}^{\dagger}(\mathbf{r})=\sum_{l\beta}e^{-i\eta\Delta\mathbf{K}_{l}\cdot\mathbf{R}}\int
d^{2}\mathbf{r}\
w_{l\beta,\alpha}^{(\eta)}(\mathbf{r}-\mathbf{R})c_{l,\beta,\eta,s}^{\dagger}(\mathbf{r})\
.$ (S55)

Under the translation $T_{\Delta\mathbf{R}}$ the operator
$f_{\mathbf{R}\alpha\eta s}^{\dagger}$ transforms as

$T_{\Delta\mathbf{R}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}T_{\Delta\mathbf{R}}^{-1}=f_{\mathbf{R}+\Delta\mathbf{R}\alpha\eta
s}^{\dagger}\ ,$ (S56)

which follows Eq. S37 directly. Correspondingly, the momentum space operators
can be defined as

$f_{\mathbf{k}\alpha\eta
s}^{\dagger}=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{i\mathbf{k}\cdot\mathbf{R}}f_{\mathbf{R}\alpha\eta
s}^{\dagger}=\sum_{\mathbf{Q}\beta}c_{\mathbf{k},\mathbf{Q},\beta,\eta,s}^{\dagger}\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})\
,$ (S57)

with $\widetilde{v}^{(\eta)}_{\mathbf{Q}\beta,\alpha}(\mathbf{k})$ from Eq.
S45. Since $f_{\mathbf{R}\alpha\eta s}^{\dagger}$ creates the Wannier states,
they follows the same representation (Eq. S36) under symmetry operations

$\hat{g}f_{\mathbf{k}\alpha\eta
s}^{\dagger}\hat{g}^{-1}=\sum_{\alpha^{\prime}\eta^{\prime}}f_{g\mathbf{k}\alpha^{\prime}\eta^{\prime}s}^{\dagger}D^{f}_{\alpha^{\prime}\eta^{\prime},\alpha\eta}(g),\qquad\hat{g}f_{\mathbf{R}\alpha\eta
s}^{\dagger}\hat{g}^{-1}=\sum_{\alpha^{\prime}\eta^{\prime}}f_{g\mathbf{R}\alpha^{\prime}\eta^{\prime}s}^{\dagger}D^{f}_{\alpha^{\prime}\eta^{\prime},\alpha\eta}(g),$
(S58)

with $D^{f}(g)$ given by Eq. S36.

Hereafter, we refer to $f_{\mathbf{R},\alpha,\eta,s}^{\dagger}$ as local
orbitals. They transform as $p_{x}\pm ip_{y}$ orbitals, which form the 2D
irreducible representation $[E]_{a}$ at the $1a$ position (Tables S1 and S2),
under the crystalline symmetries.

###  S2.2 The topological conduction bands

As discussed in the last subsection, in the BM model the six low energy states
per valley per spin form the representation $\Gamma_{1}\oplus\Gamma_{2}\oplus
2\Gamma_{3}$ around the $\Gamma_{M}$ point (Fig. S2), where
$\Gamma_{1}\oplus\Gamma_{2}$ are contributed by the middle two bands and
$2\Gamma_{3}$ are contributed by the higher energy bands. The Wannier
functions (in each valley-spin sector) constructed in the last subsection can
span the middle two bands for momenta away from $\Gamma_{M}$, where the
overlap with the flat bands is finite (Fig. S3(a)). At the $\Gamma_{M}$ point,
the Wannier functions span a $\Gamma_{3}$ representation. Thus, in order to
reproduce the correct band structure around the $\Gamma_{M}$ point, we add
four additional degree of freedom that form the representation
$\Gamma_{1}\oplus\Gamma_{2}\oplus\Gamma_{3}$ at $\Gamma_{M}$. We denote the
four states as $c_{\mathbf{k},a,\eta,s}^{\dagger}$ ($a=1,2,3,4$). Formally
they can be written as

$c_{\mathbf{k},a,\eta,s}^{\dagger}=\sum_{\mathbf{Q}\beta}\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})c_{\mathbf{k},\mathbf{Q},\beta,\eta,s}^{\dagger}\
,$ (S59)

where $\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$ is to be
determined. As shown in next subsection, $c_{\mathbf{k},a,\eta,s}^{\dagger}$
has a huge kinetic energy for large $\mathbf{k}$. Therefore, in this work we
will only consider $c_{\mathbf{k},a,\eta,s}^{\dagger}$ within the cutoff
$|\mathbf{k}|<\Lambda_{c}$. Now we determine the coefficients
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$. We define the
projector to the six lowest bands in the two valleys as

$P_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\beta}^{(\eta)}(\mathbf{k})=\sum_{n=\pm
1,\pm 2,\pm
3}u_{\mathbf{Q}\alpha,n}^{(\eta)}(\mathbf{k})u_{\mathbf{Q}^{\prime}\beta,n}^{(\eta)*}(\mathbf{k}),$
(S60)

where $u_{\mathbf{Q}\alpha,n}^{(\eta)}(\mathbf{k})$ is the eigenstate of the
BM Hamiltonian (Eq. S5) for the $|n|$-th band above ($n>0$) or below ($n<0$)
the charge neutrality point. We then define the projector to the two Wannier
functions of the local orbitals as

$Q_{\mathbf{Q}\alpha,\mathbf{Q}^{\prime}\beta}^{(\eta)}(\mathbf{k})=\sum_{\alpha=1,2}\widetilde{v}_{\mathbf{Q}\alpha,\alpha}^{(\eta)}(\mathbf{k})\widetilde{v}_{\mathbf{Q}^{\prime}\beta,\alpha}^{(\eta)*}(\mathbf{k})\
.$ (S61)

Since the two Wannier functions are constructed from the lowest six bands,
there must be
$P^{(\eta)}(\mathbf{k})Q^{(\eta)}(\mathbf{k})P^{(\eta)}(\mathbf{k})=Q^{(\eta)}(\mathbf{k})$.
This is guaranteed by the definition of our Wannier functions: The Bloch
states of the Wannier functions
$\widetilde{v}_{\mathbf{Q}\alpha,\alpha}^{(\eta)}(\mathbf{k})$ are linear
combinations of $u_{\mathbf{Q}\alpha,n}^{(\eta)}(\mathbf{k})$ ($n=\pm 1,\pm
2,\pm 3$), as explained after Eq. S42. The projector to the four low-energy
states in the remaining bands is given by
$P^{(\eta)}(\mathbf{k})-Q^{(\eta)}(\mathbf{k})$. The states
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$ can be obtained as
eigenstates of $P^{(\eta)}(\mathbf{k})-Q^{(\eta)}(\mathbf{k})$ with the
eigenvalue 1.

We emphasize that even though we have formally introduced
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$, in practical
calculations, e.g., Section S2.3, we will always approximate
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$ by
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(0)$ for low energy states, as the
$c$-electrons are only relevant for a small $\mathbf{k}$ region around
$\Gamma_{M}$.

There is a gauge degree of freedom of
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$. We fix the gauge by
requiring the representation matrices under the symmetry transformations

$\hat{g}c_{\mathbf{k},a\eta
s}^{\dagger}\hat{g}^{-1}=\sum_{a^{\prime}\eta^{\prime}}c_{g\mathbf{k},a^{\prime}\eta^{\prime}s}^{\dagger}D^{c}_{a^{\prime}\eta^{\prime},a\eta}(g)$
(S62)

to be

$\displaystyle D^{c}(T)=\sigma_{0}\tau_{x}\oplus\sigma_{0}\tau_{x},\quad
D^{c}(C_{3z})=e^{i\frac{2\pi}{3}\sigma_{z}\tau_{z}}\oplus\sigma_{0}\tau_{0},\quad
D^{c}(C_{2x})=\sigma_{x}\tau_{0}\oplus\sigma_{x}\tau_{0},\quad
D^{c}(C_{2z}T)=\sigma_{x}\tau_{0}\oplus\sigma_{x}\tau_{0}\ .$

such that the $a=1,2$ and $a=3,4$ states form the $\Gamma_{3}$ and
$\Gamma_{1}\oplus\Gamma_{2}$ representations at $\mathbf{k}=0$, respectively.
For example, the $C_{3z}$ representation matrices of the $\Gamma_{3}$ and
$\Gamma_{1}\oplus\Gamma_{2}$ states are given by
$e^{i\frac{2\pi}{3}\sigma_{z}\tau_{z}}$ and $\sigma_{0}\tau_{0}$,
respectively. Since the four states form three different representations at
$\Gamma_{M}$, the states
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$ can be uniquely
determined at $\mathbf{k}=0$. However, for $\mathbf{k}\neq 0$, there is no
easy method to uniquely fix the gauge of the four states. Nevertheless, in
practical calculations, we only need
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(0)$ for the k$\cdot$p expansion.
Using the gauge fixed by Section S2.2 and the algebra between $P$ and the
other symmetries (Eq. S19)

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

we can further determine the representation of $P$. To satisfy $P^{2}=-1$ and
$[P,C_{2z}T]=0$, the representation of $P$ in either the $\Gamma_{3}$ or
$\Gamma_{1}\oplus\Gamma_{2}$ space must have the form $\pm
i\sigma_{z}\tau_{0,z}$. Representations involving $\tau_{x,y}$ are not
considered because $P$ preserves the valley index. The considered
representations, $\pm i\sigma_{z}\tau_{0,z}$, already satisfy
$\\{P,C_{2x}\\}=0$. To commute with $T$, the only allowed form is $\pm
i\sigma_{z}\tau_{z}$. Numerically, we find the representation of $P$ (and
$C_{2z}P$ for later use) as

$D^{c}(P)=(-i\sigma_{z}\tau_{z})\oplus(-i\sigma_{z}\tau_{z}),\quad
D^{c}(C_{2z}P)=(i\sigma_{y}\tau_{y})\oplus(i\sigma_{y}\tau_{y})\ .$ (S65)

As explained after Eq. S36, for fixed $D^{c}(C_{3z})$, the sign of $D^{c}(P)$
is not a gauge choice. One may notice that the first two-by-two block of
$D^{c}(P)$, which spans the $\Gamma_{3}$ irrep, has an opposite sign with
$D^{f}(P)$ of the local $f$-orbitals (Eq. S36), which also span a $\Gamma_{3}$
irrep. Now we explain that they have to have opposite signs. We denote the
$\Gamma_{3}$ irreps formed by the positive energy bands ($n=2,3$) and the
negative energy bands ($n=-2,-3$) of the BM model as $\Gamma_{3+}$ and
$\Gamma_{3-}$, respectively. $\Gamma_{3+}$ and $\Gamma_{3-}$ transform into
each other under $P$. Thus, in the basis of $\Gamma_{3+}$ and $\Gamma_{3-}$,
in each eigen-subspace of $C_{3z}$, the $P$ operator is off-diagonal. To be
explicit, we can write the representation matrices of $C_{3z}$ and $P$

$\begin{pmatrix}e^{i\frac{2\pi}{3}}\sigma_{0}&0\\\
0&e^{-i\frac{2\pi}{3}}\sigma_{0}\end{pmatrix},\qquad\begin{pmatrix}X&0\\\
0&X^{\prime}\end{pmatrix}\ ,$ (S66)

respectively. The basis is ordered as: the first component of $\Gamma_{3+}$,
the first component of $\Gamma_{3-}$, the second component of $\Gamma_{3+}$,
the second component of $\Gamma_{3-}$, such that the first (second) two bases
have the $C_{3z}$ eigenvalue $e^{i\frac{2\pi}{3}}$ ($e^{-i\frac{2\pi}{3}}$).
$X$ and $X^{\prime}$ are off-diagonal such that $P$ interchanges the two
$\Gamma_{3}$ representations. After we diagonalize the $P$ operator, the two
$P$-eigenvalues in each eigen-subspace of $C_{3z}$ in each valley must sum to
zero. Hence the sign of $D^{c}(P)$, which spans a diagonalized $\Gamma_{3}$,
must be opposite to the sign of $D^{f}(P)$, which spans the other
$\Gamma_{3}$.

Hereafter, we refer to $c_{\mathbf{k},a,\eta,s}^{\dagger}$ ($a=1,2,3,4$) as
conduction bands.

###  S2.3 Single-particle Hamiltonian of the topological conduction bands

In this subsection we use the first-quantized formalism for simplicity. We can
divide the single-particle Hamiltonian in the low energy space into four
parts:

$H^{(f,\eta)}(\mathbf{k})=Q^{(\eta)}(\mathbf{k})h^{(\eta)}(\mathbf{k})Q^{(\eta)}(\mathbf{k}),\qquad
H^{(c,\eta)}(\mathbf{k})=(P^{(\eta)}(\mathbf{k})-Q^{(\eta)}(\mathbf{k}))h^{(\eta)}(\mathbf{k})(P^{(\eta)}(\mathbf{k})-Q^{(\eta)}(\mathbf{k}))\
,$ (S67)
$H^{(fc,\eta)}(\mathbf{k})=Q^{(\eta)}(\mathbf{k})h^{(\eta)}(\mathbf{k})(P^{(\eta)}(\mathbf{k})-Q^{(\eta)}(\mathbf{k})),\qquad
H^{(cf,\eta)}(\mathbf{k})=(P^{(\eta)}(\mathbf{k})-Q^{(\eta)}(\mathbf{k}))h^{(\eta)}(\mathbf{k})Q^{(\eta)}(\mathbf{k})\
,$ (S68)

where $h^{(\eta)}(\mathbf{k})$ is the BM model (Eq. S5),
$P^{(\eta)}(\mathbf{k})$ (Eq. S60) and $Q^{(\eta)}(\mathbf{k})$ (Eq. S61) are
projectors to the lowest six bands (per spin valley) and the Wannier
functions, respectively. The first term is the single-particle Hamiltonian of
the local orbitals, which almost vanishes ($\sim 0.15$meV for
$w_{0}/w_{1}=0.8$) as discussed in Section S2.1 and shown in Table S4. The
second term is the Hamiltonian of the conduction bands. The third and fourth
terms are the couplings between the local orbitals and the conduction bands.

We consider the k$\cdot$p expansion of the projected BM model
$H^{(c,\eta)}(\mathbf{k})$ around $\mathbf{k}=0$ ($a,a^{\prime}=1\cdots 4$)

$H^{(c,\eta)}_{aa^{\prime}}(\mathbf{k})=\langle\widetilde{u}_{a}^{(\eta,)}(0)|h^{(\eta)}(\mathbf{k})|\widetilde{u}_{a^{\prime}}^{(\eta)}(0)\rangle\
.$ (S69)

We emphasize that the usual k$\cdot$p expansion Winkler (2003) is made on the
basis $\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(0)$ rather than
$\widetilde{u}^{(\eta)}_{\mathbf{Q}\beta,a}(\mathbf{k})$. This basis choice
has two advantages: (i) It is exact if we keep an infinite number of bands
($a=1\cdots\infty$), (ii) It has simple $\mathbf{k}$-dependence. We find that
in our problem keeping four conduction bands is a good enough approximation.
Hereafter, we alway choose the Bloch states at $\mathbf{k}=0$ as the basis
whenever a k$\cdot$p expansion is needed. Since $h^{(\eta)}(\mathbf{k})$ has
only zeroth order and linear order terms in $\mathbf{k}$,
$H^{(c,\eta)}(\mathbf{k})$ also only has constant terms and linear terms in
$\mathbf{k}$. As proved in the following paragraphs, to linear order terms of
$\mathbf{k}$ the effective Hamiltonian under the gauge Sections S2.2 and S65
has the form

$H^{(c,+)}(\mathbf{k})=\begin{pmatrix}0&v_{\star}(k_{x}\sigma_{0}+ik_{y}\sigma_{z})\\\
v_{\star}(k_{x}\sigma_{0}-ik_{y}\sigma_{z})&M\sigma_{x}\end{pmatrix}\ ,$ (S70)
$H^{(c,-)}(\mathbf{k})=H^{(c,+)*}(-\mathbf{k})\ ,$ (S71)

where the $C_{3z}$ eigenvalues of the four basis in the valley $\eta=+$ are
$e^{i\frac{2\pi}{3}}$, $e^{-i\frac{2\pi}{3}}$, 1, 1, respectively. The
parameters for different $w_{0}/w_{1}$ obtained from Eq. S69 are tabulated in
Table S4. Physical meaning of $M$ is the splitting between the $\Gamma_{1}$
and $\Gamma_{2}$ states.