# The correlated insulators of magic angle twisted bilayer graphene at zero
and one quantum of magnetic flux: a tight-binding study

Miguel Sánchez Sánchez<EMAIL_ADDRESS>Tobias Stauber
<EMAIL_ADDRESS>Instituto de Ciencia de Materiales de Madrid CSIC,
Madrid (Spain)

###### Abstract

Magic angle twisted bilayer graphene (MATBG) has become one of the prominent
topics in Condensed Matter during the last few years, however, fully atomistic
studies of the interacting physics are missing. In this work, we study the
correlated insulator states of MATBG in the setting of a tight-binding model,
under a perpendicular magnetic field of $0$ and $26.5$ T, corresponding to
zero and one quantum of magnetic flux per unit cell. At zero field and for
dopings of two holes ($\nu=-2$) or two electrons ($\nu=+2$) per unit cell, the
Kramers intervalley coherent (KIVC) order is the ground state at the Hartree-
Fock level, although it is stabilized by a different mechanism to that in
continuum model. At charge neutrality, the spin polarized state is competitive
with the KIVC due to the on-site Hubbard energy. We obtain a strongly
electron-hole asymmetric phase diagram with robust insulators for electron
filling and metals for negative filling. In the presence of magnetic flux, we
predict an insulator with Chern number $-2$ for $\nu=-2$, a spin polarized
state at charge neutrality and competing insulators with Chern numbers $+2$
and $0$ at $\nu=+2$. The stability of the $\nu=+2$ insulators is determined by
the screening environment, allowing for the possibility of observing a
topological phase transition.

††preprint: APS/123-QED

##  I Introduction

Magic angle twisted bilayer graphene (MATBG) is a two dimensional quantum
material[1] that exhibits a plethora of exotic phases ranging from
superconductors[2, 3, 4, 5] to strange metals[6, 7], passing through
integer[8, 9] and fractional[10] Chern Insulators. It constitutes a remarkable
platform for the understanding of the many-body problem in Condensed Matter
and the interplay of strong interactions and topology[11, 12, 13], and
inaugurated the field of moiré materials[14, 15, 16, 17].

The correlated insulators arise when the doping level is such that the number
of electrons per unit cell is an integer number. They were discovered before
the superconductivity and other phases[18], and are one of the most studied
phenomena in twisted bilayer graphene (TBG)[19, 20, 21, 22, 4, 3, 23, 9, 8, 5,
24, 25, 18, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43].

On another hand, crystalline systems under magnetic fields are controlled by
the scale given by the magnetic flux quantum $\Phi_{0}=h/e$[44, 45]. When the
magnetic field is such that the flux per unit cell is comparable to
$\Phi_{0}$, the system is in the ’Hofstadter regime’ and the picture of Landau
levels is replaced by a reentrant band structure[46, 47]. In typical materials
such magnetic fields are of the order of $10^{4}$ T, but in MATBG the large
moiré unit cell allows accessible fields of the order of $30$ T. In MATBG the
Landau level spectrum of the correlated insulators has been studied[48, 49, 3,
4, 8, 9]. Also, at one magnetic flux quantum reentrant correlated insulators
have been predicted and observed[50, 51, 52].

On the theory side, the Bistritzer-McDonald (BM) or continuum model[53, 54] is
a low energy theory that has proven very powerful in understanding the physics
of TBG, revealing the emergent symmetries of the system that have led to the
picture of the ’$U(4)$ ferromagnets’ for the correlated insulators[35, 19, 20,
32, 37, 33, 34]. However, the model, with only a handful of parameters, cannot
capture the finer details of the spectrum and the wave functions. These
differences at low energy scales are relevant in the competition between
states.

In this work we employ a tight-binding model for MATBG. The high computational
cost, which makes atomistic studies scarce in this system[28, 27, 29, 55], is
partially bypassed by a projection onto the subspace of the low energy bands
(the ’flat bands’). The external magnetic flux is tuned to zero and one
magnetic flux quantum per unit cell, we focus on samples without strain and
leave electron-phonon coupling for future work.

The explicit breaking of the $U(4)$ symmetry of MATBG is assessed via the
particle-hole asymmetry of the flat bands[56], the intervalley Hund’s
coupling[21] and the on-site Hubbard interaction. We find that the intervalley
Hund’s coupling, relevant for the superconductivity, has antiferromagnetic
contributions due to the long range Coulomb interaction. On the other hand,
the Hubbard term contributes to a ferromagnetic coupling.

To study the spontaneous symmetry breaking in the correlated states, we
perform self-consistent Hartree-Fock simulations. For $\Phi=0$ we find
agreement with previous results on the nature of the ground state at even
filling[20, 32, 34, 21, 36, 22], but the selection mechanism of the ground
state is different. Instead of kinetic energy superexchange as expected in the
BM model, the appearance of ’inter-Chern’ order near the $\Gamma$ point
reduces the exchange energy and stabilizes the Kramers intervalley coherent
state. In addition, the Hubbard interaction makes the spin polarized state
competitive at charge neutrality. The self-consistent states are insulating at
electron doping and metallic at hole doping, signalling the experimentally
reported many-body electron-hole asymmetry[38, 4, 3].

For $\Phi=\Phi_{0}$, the system departs from the strong coupling $U(4)$
picture due to the increased bandwidth of the kinetic energy bands, in the
same manner as the strained samples at zero flux[22, 39]. Consistently for
different screening environments, we observe an insulator with Chern number
$-2$ at $\nu=-2$ and a spin polarized state at $\nu=0$. For $\nu=+2$ we
observe a topological phase transition from an insulator with Chern number
$+2$ for small screening to an intervalley coherent trivial insulator for
large screening. The Chern $+2$ insulator is compatible with the experimental
data of Ref. [52].

The paper is organized as follows. In section II we describe the tight-binding
model of TBG and the Peierls’ substitution under magnetic field, in section
III we introduce the emergent $U(4)$ symmetry at the magic angle and in
section IV we discuss the explicit breaking of $U(4)$ in the lattice model.
Finally, in section V we report the correlated states obtained in the Hartree-
Fock simulations, and in section VI we draw some conclusions.

##  II The model

In graphene, the primitive vectors are $\boldsymbol{a}_{1}=a(1/2,\sqrt{3}/2)$
and $\boldsymbol{a}_{2}=a(-1/2,\sqrt{3}/2)$, with $a=\sqrt{3}a_{0}$ and
$a_{0}=0.142$ nm the carbon-carbon distance. Atoms at lattice points belong to
sublattice $A$, and their nearest neighbours displaced by
$(\boldsymbol{a}_{1}+\boldsymbol{a}_{2})/3$ to sublattice $B$.

Consider two graphene layers stacked on top of each other, at $z=-d_{0}/2$ and
$z=d_{0}/2$ respectively, being $d_{0}=0.335$ nm the interlayer distance, such
that top and bottom atoms are vertically aligned. The bottom layer is rotated
by an angle $-\theta/2$, and the top layer by $\theta/2$, with the center of
rotation being the center of one of the graphene hexagons. We choose a value
of $\theta$ that makes the twisted structure commensurate[54]. In our case, we
parametrize the angle by an integer $n_{\theta}$ such that
$\cos(\theta)=1-1/2(3n_{\theta}^{2}+3n_{\theta}+1)$. The unit vectors of the
superlattice are

$\displaystyle\boldsymbol{L}_{1}$
$\displaystyle=R_{-\theta/2}\big{(}n_{\theta}\boldsymbol{a}_{1}+(n_{\theta}+1)\boldsymbol{a}_{2}\big{)}=L_{M}(0,1),$
$\displaystyle\boldsymbol{L}_{2}$
$\displaystyle=R_{\pi/3}\boldsymbol{L}_{1}=R_{-\theta/2}\big{(}(-n_{\theta}-1)\boldsymbol{a}_{1}+\big{(}2n_{\theta}+1)\boldsymbol{a}_{2}\big{)},$
(1)

with $R_{\alpha}$ a rotation by angle $\alpha$ and $L_{M}$ the lattice
constant. The reciprocal vectors are given by

$\displaystyle a_{0}\boldsymbol{G}_{1}$
$\displaystyle=G_{\theta}R_{-\theta/2}\big{(}(3n_{\theta}+1)\boldsymbol{a}_{1}+\boldsymbol{a}_{2}\big{)},$
$\displaystyle a_{0}\boldsymbol{G}_{2}$
$\displaystyle=R_{-2\pi/3}(\boldsymbol{G}_{1})$
$\displaystyle=G_{\theta}R_{-\theta/2}\big{(}-(3n_{\theta}+2)\boldsymbol{a}_{1}+(3n_{\theta}+1)\boldsymbol{a}_{2}\big{)},$
(2)

where $G_{\theta}=\frac{4\pi}{3a_{0}}(9n_{\theta}^{2}+9n_{\theta}+3)^{-1}$.
The magic angle is approximately given by $n_{\theta}=31$ ($1.05^{\circ}$),
corresponding to a Moiré lattice constant of $L_{M}=13.4$ nm and $11908$ atoms
in the unit cell.

The point group of this structure is the dihedral group $D_{6}$, generated by
six-fold rotations around the $z$ axis, $C_{6z}$, and two-fold rotations
around the $y$ axis, $C_{2y}$, leaving the origin fixed. The combined
operation $C_{6z}^{3}=C_{2z}$ amounts to a two-fold rotation around the $z$
axis, and $C_{2x}=C_{2y}C_{2z}=C_{2x}$ to a two-fold rotation around the $x$
axis. The spin-orbit coupling being small, spinless time-reversal
$\mathcal{T}$ is also a symmetry.

Figure 1: Top view of twisted bilayer graphene for a twist angle of
$\theta=9.43^{\circ}$. The hexagonal Wigner-Seitz cell is indicated. The
center of the unit cell is locally $AA$ stacked (the two layers are on top of
each other), and the corners are $AB$, or $BA$ stacked (only $A$ atoms are on
top of $B$ atoms, or vice versa).

Figure 2: Brillouin zone of MATBG. Left: Brillouin zone of the decoupled
monolayers in purple and green, and Moiré Brillouin zone of TBG in black. Low
energy TBG states belong to valley $K$ or $K^{\prime}$ of each monolayer. The
interlayer tunneling couples both graphene sheets, with negligible mixing
between valleys. Right: Redefined high-symmetry momenta at one magnetic flux
quantum in the periodic Landau gauge. The BZ gets shifted by
$(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2$ (see Appendix A).

Lattice relaxation is included via in-plane distortions following the model of
Ref.[57]. The effect of relaxation is to enlarge the AB and BA regions and
reduce the AA regions of the Moiré pattern (see Fig. 1), preserving all the
crystallographic symmetries.
We employ the Slater-Koster parametrization of the hopping integral of
Ref.[58], with a $p_{z}$ orbital per carbon atom and spin, giving the
Hamiltonian

$\displaystyle
H_{0}=\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}},s}t(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})c^{\dagger}_{\boldsymbol{i}s}c_{\boldsymbol{i}s},$
(3)

$c^{\dagger}_{\boldsymbol{i},s}$ being the creation operator of an electron
with spin $s$ at position $\boldsymbol{r_{i}}$. The hopping integral is
decomposed into $\sigma$ and $\pi$-bond hoppings,

$\displaystyle t(\boldsymbol{r})=-$ $\displaystyle
V_{pp\pi}(r)\Bigg{(}1-\bigg{(}\frac{\boldsymbol{r}\cdot\boldsymbol{\hat{z}}}{r}\bigg{)}^{2}\Bigg{)}+V_{pp\sigma}(r)\bigg{(}\frac{\boldsymbol{r}\cdot\boldsymbol{\hat{z}}}{r}\bigg{)}^{2},$
$\displaystyle V_{pp\pi}(r)=V_{pp\pi}^{0}e^{-(r-a_{0})/r_{0}},$ $\displaystyle
V_{pp\sigma}(r)=V_{pp\sigma}^{0}e^{-(r-d_{0})/r_{0}},$ (4)

with the parameters $V^{0}_{ppp\pi}=2.7$ eV, $V^{0}_{pp\sigma}=0.48$ eV and
$r_{0}=0.0453$ nm.

The Coulomb interaction is implemented by the double-gated potential

$\displaystyle
V=\frac{1}{2}\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}}s_{i}s_{j}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}}):c^{\dagger}_{\boldsymbol{i},s_{i}}c_{\boldsymbol{i},s_{i}}c^{\dagger}_{\boldsymbol{j},s_{j}}c_{\boldsymbol{j},s_{j}}:,$
$\displaystyle
V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})=\frac{e^{2}}{4\pi\epsilon_{0}\epsilon}\sum_{n}\frac{(-1)^{n}}{||\boldsymbol{r_{i}}-\boldsymbol{r_{j}}+n\xi\boldsymbol{\hat{z}}||},$
(5)

which applies for the experimental setups where two metallic plates are placed
at $z=\pm\xi/2$. Unless stated otherwise, we set $\xi=10$ nm. The dielectric
constant $\epsilon$ accounts for the screening due to the substrate and
internal screening due to the electrons. The interaction is normal ordered[59]
with respect to the ground state of two decoupled graphene layers at charge
neutrality. This choice of normal ordering is also called graphene subtraction
scheme[30, 31]. Under magnetic field, we do not include the Zeeman shift when
calculating the graphene state, so that the spin imbalances come entirely from
the flat band physics. The on-site Hubbard term is also considered

$\displaystyle
H_{U}=U\sum_{\boldsymbol{r_{i}}}:c^{\dagger}_{\boldsymbol{i}\uparrow}c_{\boldsymbol{i}\uparrow}c^{\dagger}_{\boldsymbol{i}\downarrow}c_{\boldsymbol{i}\downarrow}:.$
(6)

The total Hamiltonian is then $H=H_{0}+V+H_{U}$.

### Minimal coupling to an external magnetic field

At nonzero magnetic field, the Peierls’ substitution[60] adds a phase to the
hopping elements,

$\displaystyle t(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})\to
t(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})e^{i\theta_{\boldsymbol{i},\boldsymbol{j}}},$
$\displaystyle\theta_{\boldsymbol{i},\boldsymbol{j}}=\frac{2\pi}{\Phi_{0}}\int_{\boldsymbol{r_{i}}\to\boldsymbol{r_{j}}}\boldsymbol{A}(\boldsymbol{r^{\prime}})\cdot
d\boldsymbol{r^{\prime}},$ (7)

where $\Phi_{0}=h/e$ is the quantum of magnetic flux, and the line integral
goes from $\boldsymbol{r_{i}}$ to $\boldsymbol{r_{j}}$ in a straight line if
the orbitals are well localized[46].
In the presence of magnetic flux, the translation operators pick up an
Aharonov-Bohm phase. They act on the single-particle states as[45]

$\displaystyle\tilde{T}_{1}$ $\displaystyle=\sum_{\boldsymbol{r_{i}}}e^{-2\pi
i\xi_{2\boldsymbol{i}}\phi-i\theta_{\boldsymbol{i},\boldsymbol{i}+\boldsymbol{L_{1}}}}c^{\dagger}_{\boldsymbol{i}+\boldsymbol{L}_{1}}c_{\boldsymbol{i}},$
$\displaystyle\tilde{T}_{2}$ $\displaystyle=\sum_{\boldsymbol{r_{i}}}e^{2\pi
i\xi_{1\boldsymbol{i}}\phi-i\theta_{\boldsymbol{i},\boldsymbol{i}+\boldsymbol{L_{2}}}}c^{\dagger}_{\boldsymbol{i}+\boldsymbol{L_{2}}}c_{\boldsymbol{i}}.$
(8)

$\xi_{\boldsymbol{i}1}$ and $\xi_{\boldsymbol{i}2}$ are defined by
$\boldsymbol{r_{i}}=\xi_{\boldsymbol{i}1}\boldsymbol{L_{1}}+\xi_{\boldsymbol{i}2}\boldsymbol{L_{2}}$,
and $\phi=\Phi/\Phi_{0}=BA_{M}/\Phi_{0}$ is the flux per moiré unit cell in
units of $\Phi_{0}$.

It can be shown that
$[\mathcal{H},\tilde{T_{1}}]=[\mathcal{H},\tilde{T_{2}}]=0$ and
$\tilde{T}_{1}\tilde{T}_{2}=e^{-2\pi i\phi}\tilde{T}_{2}\tilde{T}_{1}$, so the
translational symmetries are broken in general. However, if $\phi$ is a
rational number $p/q$ one can choose the set of commuting operators
$\\{\tilde{T}_{1},\tilde{T}_{2}^{q}\\}$, or
$\\{\tilde{T}_{1}^{q},\tilde{T}_{2}\\}$ and diagonalize them simultaneously
with the Hamiltonian. Translational symmetry is then recovered at rational
fluxes with a unit cell that is $q$ times larger than at zero flux, and the
Bloch waves are generalized to magnetic waves having good $\tilde{T}_{1}$ and
$\tilde{T}_{2}^{q}$ quantum numbers.

In the periodic Landau gauge[61]

$\displaystyle\boldsymbol{A}$
$\displaystyle(\boldsymbol{r})=\frac{\Phi}{2\pi}\Bigg{(}\xi_{1}\boldsymbol{G}_{2}-2\pi\boldsymbol{\nabla}\big{(}\xi_{2}\left\lfloor{\xi_{1}+\epsilon}\right\rfloor\big{)}\Bigg{)}$
$\displaystyle=\frac{\Phi}{2\pi}\Bigg{(}-\xi_{2}\sum_{n}\delta(\xi_{1}-n+\epsilon)\boldsymbol{G}_{1}+(\xi_{1}-\left\lfloor{\xi_{1}+\epsilon}\right\rfloor\boldsymbol{G}_{2}\Bigg{)},$
(9)

the phases of the translation operators $\tilde{T}_{2}^{q},$ $\tilde{T}_{1}$
cancel and the Bloch waves have the same form as in zero flux
($\lfloor...\rfloor$ is the floor function). The infinitesimal $\epsilon$
prevents ambiguities if some atoms lie at integer values of $\xi_{1}$. The
momentum $\boldsymbol{k}$ takes the possible values in the magnetic Brillouin
zone of the dual lattice with lattice vectors $\boldsymbol{G}_{1}$ and
$\boldsymbol{G}_{2}/q$. Under magnetic flux, time reversal $\mathcal{T}$ and
rotations $C_{2y}$, $C_{2x}$ reverse the sign of the external field, but the
rotations around the $z$ axis are preserved[62].

Besides orbital effects, the Zeeman energy $-g\mu_{B}Bs_{z}/\hbar$ ($g=2$ is
the gyromagnetic ratio of the electron and $\mu_{B}$ the Bohr magneton) is
also taken into account. For $26.5$ T it amounts to $\pm 1.535$ meV.

$\boldsymbol{\Phi=0}$ $\boldsymbol{\Phi=\Phi_{0}}$

Figure 3: Band structure of MATBG ($\boldsymbol{\theta=1.05^{\circ}}$) along
the high symmetry $K_{M}\Gamma_{M}M_{M}K_{M}$ line. Left: flat bands at zero
field. Right: flat bands at $B=26.5$ T. The Zeeman energy produces a splitting
of $\sim 3$ meV.

a) b)

Figure 4: Integrated Berry curvature at $\boldsymbol{B=26.5}$ T. We calculate
the integrated Berry curvature on a $18\times 18$ grid. a) Valley $K$ valence
band (Chern number $-1$). b) Valley $K$ conduction band (Chern number $+1$).
$C_{2z}$ symmetry interchanges the valleys of the conduction/valence bands
(see the text for the definition of the valley charge).

### Results of the non interacting model

In Fig. 3 we plot the spectrum of MATBG for $n_{\theta}=31$. At the charge
neutrality point there are 8 quasi flat bands ($2$ layers $\times$ $2$ valleys
$\times$ $2$ spin polarizations) with a bandwidth of around $10$ meV and well
separated from the ’remote’ bands. The filling is parametrized by
$\nu\in(-4,4)$, being $\nu=0$ the neutrality point and $\nu=+4$ $(-4)$ the
filled (empty) flat bands. When the external magnetic field is $26.5$ T, we
have $\Phi=\Phi_{0}$ and the flat bands are reentrant. The symmetry
$C_{2z}\mathcal{T}$ that preserves the gapless Dirac points in zero flux[63]
is broken, and the Dirac cones are gapped with a Chern number of $-2$[47]. In
the flat bands the kinetic energy is small and the interacting physics
dominates, giving rise to the rich phase diagram of MATBG.

In Fig. 4 we plot the Berry curvatures of the bands at $26.5$ T. In the
discretized Brillouin zone, we compute the integrated Berry curvature on the
parallelograms defined by the discretization, following Ref. [64]. We see an
emergent symmetry relating the Berry curvatures of valence and conduction
bands.

##  III The U(4) symmetry

A starting point for the understanding the physics of MATBG is the projected
limit, in which the Fermi sea of the remote bands is ’frozen’ and the
scattering is restricted to states within the flat bands. The correlated
insulators are ’generalized ferromagnets’ that spontaneously break a $U(4)$
(or a larger $U(4)\times U(4)$) global symmetry. For a detailed discussion on
this symmetry, we refer the reader to Refs. [12, 21, 33, 50]. Here, we
describe the most ’physical’ subgroups. The full group is generated from
these.

  * •

$U(1)_{v}$. The non interacting eigenstates belong to valley $K$ (valley
charge +1) or $K^{\prime}$ (valley charge -1) of each of the graphene
monolayers, see Fig. 2. Charge non conserving terms of the screened Coulomb
interaction are suppressed, hence $U(1)_{v}$ is a symmetry of the interacting
system. We write $\eta=+1(-1)$ for valley $K$($K^{\prime}$). The valley charge
operator is denoted by $\tau_{z}$, with eigenvalues
$\langle\tau_{z}\rangle=\eta$. Electric charge $U(1)_{c}$ is also a symmetry.

  * •

$SU(2)_{K}\times SU(2)_{K^{\prime}}$. The exchange integral of pairs of states
with different valley charges is suppressed due to the big momentum transfer,
and only the Hartree term contributes to the energy. This implies a symmetry
of independent spin rotations in each valley, which do not change the total
density. At nonzero perpendicular magnetic field the Zeeman effect will break
the degeneracy and fix the quantization axis to the $z$ axis.

  * •

The particle-hole (p-h) symmetry. In the spirit of the continuum model we can
write a generic wave function of valley $\eta$ as

$\displaystyle\Psi_{\eta}(\boldsymbol{r_{i}})=e^{i\eta
n_{\theta}\boldsymbol{G_{2}}\cdot\boldsymbol{r_{i}}}\Psi_{\eta\sigma
l}(\boldsymbol{r_{i}})$ (10)

where $e^{i\eta n_{\theta}\boldsymbol{G_{2}}\cdot\boldsymbol{r_{i}}}$ is the
rapidly oscillating valley phase
($n_{\theta}\boldsymbol{G}_{2}\approx\boldsymbol{K}$, the $K$ point of
graphene), and $\Psi_{\eta\sigma l}(\boldsymbol{r_{i}})$ is a smooth envelope
that depends on the sublattice $\sigma$ and layer $l$ of the point
$\boldsymbol{r_{i}}$, evaluated at that point. The particle-hole operator
$C_{2z}P$ (actually, in the language of Refs. [65, 50], the combined operator
of $C_{2z}$ and the particle-hole operator $P$) is a hermitian operator that
squares to the identity and interchanges the graphene valleys, the sublattice
and the layer, and is trivial on the spin. It acts on the wave functions as

$\displaystyle\hskip
22.76228ptC_{2z}P\big{(}\Psi_{\eta}\big{)}(\boldsymbol{r_{i}})=\eta
s_{l}e^{-i\eta
n_{\theta}\boldsymbol{G_{2}}\cdot\boldsymbol{r_{i}}}e^{-i\eta\boldsymbol{G_{2}}\cdot\boldsymbol{r_{i}}}\Psi_{\eta\overline{\sigma}\overline{l}}(\boldsymbol{r_{i}}),$
(11)

where $\overline{\sigma}$ and $\overline{l}$ denote the opposite sublattice
and layer to those of $\boldsymbol{r_{i}}$ and $s_{l}=+1$ if $l$ is the top
layer and $-1$ if $l$ is the bottom layer. In the continuum theory, $C_{2z}P$
commutes with the Coulomb interaction and anticommutes with the kinetic
energy, hence the name particle-hole operator. As such, it is the generator of
a $U(1)$ subgroup in the flat limit, i.e. when the kinetic energy is
negligible compared to the Coulomb energy and set to zero.

The $U(1)$ subgroup generated by $C_{z}P$ and $SU(2)_{K}\times
SU(2)_{K^{\prime}}\times U(1)_{v}\times U(1)_{c}\simeq U(2)_{K}\times
U(2)_{K^{\prime}}$ do not commute. They are subgroups of the ubiquitous $U(4)$
symmetry of TBG. This $U(4)$ group can be further enlarged to $U(4)\times
U(4)$ if we include another $U(1)$ generated by

  * •

the sublattice operator,

$\displaystyle C=\sum_{\boldsymbol{i}\in
A,s}c^{\dagger}_{\boldsymbol{i}s}c_{\boldsymbol{i}s}-\sum_{\boldsymbol{i}\in
B,s}c^{\dagger}_{\boldsymbol{i},s}c_{\boldsymbol{i},s}.$ (12)

This operator generates a symmetry in the so-called ’chiral’ limit[66], where
the projected $C$ operator in the flat band manifold has eigenvalues $\pm 1$ ,
i.e. there exists a perfectly sublattice polarized basis of the flat bands. In
the real system the polarization is around $0.6$-$0.8$[56], and the symmetry
is moderately broken. On the other hand, at one magnetic flux quantum the
chiral limit is topologically distinct from the real system[50].

The reader might have noticed that both the valley charge and the $C_{2z}P$
operators are emergent in the continuum theory and do not have a direct
analogue in the lattice. In Appendix B we describe our implementation of the
valley charge and $C_{2z}P$ in the lattice model.

### The irrep basis

The natural basis of the flat bands in the strong coupling analysis is the so-
called ’irrep’ basis, with the defining property

$\displaystyle\langle\boldsymbol{k}\eta\lambda|C_{2z}P|\boldsymbol{k}\eta^{\prime}\lambda^{\prime}\rangle=[\tau_{y}]_{\eta\eta^{\prime}}[\lambda_{z}]_{\lambda\lambda^{\prime}},$
(13)

where $\lambda=\pm 1$ is the irrep number. We write $\tau_{0,x,y,z}$ and
$\lambda_{0,x,y,z}$ for the identity and Pauli matrices in valley and irrep
number space, respectively. In this subsection the spin index is omitted, and
we construct two identical copies of the irrep basis, one for each spin
polarization.

Given the property $\\{C,C_{2z}P\\}=0$, the irrep basis is equivalent to the
sublattice polarized basis that diagonalizes the projection of $C$ onto the
flat bands, which we denote by $\overline{C}(\boldsymbol{k})$ for a given
momentum $\boldsymbol{k}$. The sublattice is labeled by $\sigma=A(+1),B(-1)$,
and the identity and Pauli matrices in sublattice space by $\sigma_{0,x,y,z}$.
$A(B)$ sublattice has eigenvalue $+1(-1)$ under $\sigma_{z}$.

However, in the real system the particle-hole symmetry is broken, meaning that
$\overline{C_{2z}P}(\boldsymbol{k})$ is not unitary, as we will see. We have
to define the irrep basis in a different way.

At zero magnetic field, the sublattice basis is adiabatically connected to the
irrep basis of the p-h symmetric limit. This is, if we compute the ’closest’
unitary matrix to $\overline{C_{2z}P}(\boldsymbol{k})$ in the sublattice
polarized basis,
$[\overline{C_{2z}P}(\boldsymbol{k})]_{\eta\sigma\eta^{\prime}\sigma^{\prime}}=\langle\boldsymbol{k}\eta\sigma|C_{2z}P|\boldsymbol{k}\eta^{\prime}\sigma^{\prime}\rangle$,
we get to a very good accuracy (up to a gauge choice)

$\displaystyle\overline{C_{2z}P}(\boldsymbol{k})\Big{(}\overline{C_{2z}P}(\boldsymbol{k})\overline{C_{2z}P}(\boldsymbol{k})^{\dagger}\Big{)}^{-1/2}=\sigma_{y}\tau_{x}.$
(14)

Hence, if we identify $\lambda=\eta\sigma$ we conclude that indeed both basis
are equivalent also with p-h breaking, and they can be used indistinctly. The
sublattice polarized bands $|\boldsymbol{k}\eta\sigma\rangle$ have Chern
numbers equal to $\eta\sigma$[67].

We impose the following $C_{2z}\mathcal{T}$ gauge fixing,

$\displaystyle\langle\boldsymbol{k}\rho|C_{2z}\mathcal{T}|\boldsymbol{k}\rho^{\prime}\rangle=[\sigma_{x}\tau_{0}]_{\rho\rho^{\prime}},$
(15)

with $\rho$ the multi-index for valley and sublattice. We also constrain the
representation of $C_{2z}$ and $\mathcal{T}$,

$\displaystyle\langle[\boldsymbol{-k}]\rho|C_{2z}|\boldsymbol{k}\rho^{\prime}\rangle=$
$\displaystyle[\sigma_{x}\tau_{x}]_{\rho\rho^{\prime}}\ \text{if}\
\boldsymbol{k}\neq[-\boldsymbol{k}]\ \text{or}\
\boldsymbol{k}=\boldsymbol{\Gamma}$ $\displaystyle=$
$\displaystyle-[\sigma_{x}\tau_{x}]_{\rho\rho^{\prime}}\ \text{otherwise},$
$\displaystyle\langle[\boldsymbol{-k}]\rho|\mathcal{T}|\boldsymbol{k}\rho^{\prime}\rangle=$
$\displaystyle[\sigma_{0}\tau_{x}]_{\rho\rho^{\prime}}\ \text{if}\
\boldsymbol{k}\neq[-\boldsymbol{k}]\ \text{or}\
\boldsymbol{k}=\boldsymbol{\Gamma}$ $\displaystyle=$
$\displaystyle-[\sigma_{0}\tau_{x}]_{\rho\rho^{\prime}}\ \text{otherwise},$
(16)

$[\boldsymbol{k}]$ being the momentum equivalent to $\boldsymbol{k}$ inside
the Brillouin zone. The additional signs at parity invariant momenta are due
to a topological obstruction[33].

Notice that the wave functions are not completely defined, and there exists a
phase ambiguity redefining the states as

$\displaystyle|\boldsymbol{k}\rho\rangle$
$\displaystyle\longrightarrow[e^{i\alpha(\boldsymbol{k})\sigma_{z}\tau_{z}}]_{\rho\rho^{\prime}}|\boldsymbol{k}\rho^{\prime}\rangle,$
$\displaystyle\alpha(\boldsymbol{k})$
$\displaystyle=\alpha([-\boldsymbol{k}]).$ (17)

At one magnetic flux quantum the irrep basis is defined by

$\displaystyle\overline{C_{2z}P}(\boldsymbol{k})\Big{(}\overline{C_{2z}P}(\boldsymbol{k})\overline{C_{2z}P}(\boldsymbol{k})^{\dagger}\Big{)}^{-1/2}=\tau_{x}\lambda_{0}.$
(18)

In the continuum model, the chiral limit is topologically distinct from the
real system at one flux quantum[50], a phenomenon that is reflected in the
tight-binding model. In this case, the irrep basis is not maximally polarized,
hence the sublattice polarization and irrep character cannot be simultaneously
manifested.

We further fix the representation of $C_{2z}$ to

$\displaystyle\langle[\boldsymbol{-k}]\rho|C_{2z}|\boldsymbol{k}\rho^{\prime}\rangle=-[\tau_{y}\lambda_{y}]_{\rho\rho^{\prime}},$
(19)

with $\rho$ now the multi-index for valley and irrep.

Because $C_{2z}\mathcal{T}$ is broken due to the magnetic field, there does
not exit a remaining local symmetry to further constrain the basis. As a
consequence, in flux the irrep basis is only defined up to arbitrary
transformations $V(\boldsymbol{k})$ in both valleys

$\displaystyle|\boldsymbol{k}\eta\lambda\rangle\longrightarrow[V(\boldsymbol{k})]_{\lambda\lambda^{\prime}}|\boldsymbol{k}\eta\lambda^{\prime}\rangle,$
$\displaystyle
V^{\dagger}([-\boldsymbol{k}])\lambda_{y}V(\boldsymbol{k})=\lambda_{y}.$ (20)

##  IV EXPLICIT BREAKING OF $\boldsymbol{U(4)}$

The $U(4)$ (or $U(4)\times U(4)$) symmetry is only approximate, and is broken
in the atomistic model down to the physical $SU(2)$ of spin rotations. In this
section we study the strength of such symmetry breaking.

$\Phi=0$
$\Phi=\Phi_{0}$

Figure 5: Singular values of $\boldsymbol{\overline{C_{2z}P}}$ at
$\boldsymbol{\Phi=0}$ and $\boldsymbol{\Phi_{0}}$. They are close to $1$,
signalling a small particle-hole breaking in the flat bands. For $\Phi=0$ they
were previously computed in Ref. [56]

###  $\boldsymbol{SU(2)_{K}\times SU(2)_{K^{\prime}}}$ and p-h breaking

Assuming that deviations from $[V,C_{2z}P]=0$ are negligible also in the
lattice model, the projected flat limit enjoys the symmetry generated by the
p-h operator if the projected $C_{2z}P$ matrix is unitary. Hence, p-h breaking
can be quantified by the singular values (s.v.) of
$\overline{C_{2z}P}(\boldsymbol{k})$. Deviations from 1 of the s.v. measure
the particle-hole asymmetry of the Hilbert space.

The properties
$\overline{C_{2z}P}(\boldsymbol{k})^{\dagger}=\overline{C_{2z}P}(\boldsymbol{k})$
and $\\{\tau_{z},C_{2z}P\\}=0$ force the s.v. to be degenerate in pairs. We
plot the largest and smallest s.v. in Fig. 5 at zero and one flux quantum. The
deviations are small and similar in both cases, with mean values of around
$0.99$ and a minimal value of about $0.98$ at zero field.

Figure 6: Coulomb scattering between electrons with valley charges
$\eta_{a,b,c,d}$ and spin projections $s_{a,b,c,d}$. In the continuum
approximation $V_{abcd}\sim V(\boldsymbol{q})$, with
$\boldsymbol{q}\sim(\eta_{a}-\eta_{c})\boldsymbol{K}\sim(\eta_{d}-\eta_{b})\boldsymbol{K}$.
Lattice-scale effects distort this approximation.

Regarding $U(1)_{v}$ and $SU(2)_{K}\times SU(2)_{K^{\prime}}$, consider the
matrix elements of the Coulomb interaction depicted in Fig. 6,

$\displaystyle V_{abcd}=\langle ab|V|cd\rangle,$ (21)

where the states $|i\rangle$ have valley charge $\eta_{i}$ and spin $s_{i}$,
$i=a,b,c,d$. In the continuum theory, the matrix element is to a first
approximation equal to the Fourier transform of $V(\boldsymbol{r})$ at
momentum
$\boldsymbol{q}\sim(\eta_{a}-\eta_{c})\boldsymbol{K}\sim(\eta_{d}-\eta_{b})\boldsymbol{K}$,
with $\boldsymbol{K}$ the corner of the graphene BZ, the midpoint of
$\boldsymbol{K_{1}}$ and $\boldsymbol{K_{2}}$ in Fig. 2. In turn, $V_{abcd}$
is nonzero only if $\eta_{a}=\eta_{c}$ and $\eta_{b}=\eta_{d}$ given that
$||\boldsymbol{K}||\sim a^{-1}\gg\xi^{-1}$, where $\xi$ is the gate distance
in Eq. 5. The valley charge is conserved,
$\eta_{a}+\eta_{b}=\eta_{c}+\eta_{d}$, so $U(1)_{v}$ is a symmetry. Moreover,
the structure of the matrix elements
$V_{abcd}\propto\delta_{\eta_{a}\eta_{c}}\delta_{\eta_{b}\eta_{d}}\delta_{s_{a}s_{c}}\delta_{s_{b}s_{d}}$
exhibits the symmetry $SU(2)_{K}\times SU(2)_{K^{\prime}}$ consisting of
independent spin rotations for each valley sector.

For the atomistic model, we show in Appendix C that $U(1)_{v}$ is preserved
but $SU(2)_{K}\times SU(2)_{K^{\prime}}$ is broken. Furthermore, it is shown
that the exchange energy $-\langle ab|V|ba\rangle$ when $\eta_{a}=-\eta_{b}$
is always positive, contributing to an antiferomagnetic Hund’s coupling $J<0$
in the language of Ref. [21].

Moreover, the on-site Hubbard Hamiltonian
$U\sum_{\boldsymbol{i}}:n_{\boldsymbol{i}\uparrow}n_{\boldsymbol{i}\downarrow}:$
also incorporates symmetry breaking, favouring magnetically ordered phases.

In what follows we quantify the symmetry breaking in the manifold of ground
state candidates and discuss its importance.

### The manifold of possible ground states

For even values of $\nu$, $\nu=0,\pm 2$, the ground states in strong coupling
are ’$U(4)$ ferromagnets’. Any $U(4)$ rotation $U$ of the valley-spin
polarized states is a possible ground state[34, 20, 21, 37, 51]

$\displaystyle|\text{GS}\rangle=U\bigg{(}\prod_{\boldsymbol{k}}\prod_{j=1,\nu/2+2}d^{\dagger}_{\boldsymbol{k}\eta_{j}+1s_{j}}d^{\dagger}_{\boldsymbol{k}\eta_{j}-1s_{j}}|0\rangle\bigg{)},$
(22)

$d^{\dagger}_{\boldsymbol{k}\eta\lambda s}$ denoting the creation operator of
state $|\boldsymbol{k}\eta\lambda s\rangle$, $|0\rangle$ the state with the
filled remote bands, and $s=\uparrow,\downarrow$ the spin index. The identity
and Pauli matrices in spin will be denoted by $s_{0,x,y,z}$. The valley-spin
flavors $\eta_{j}$, $s_{j}$ can be chosen arbitrarily, as different choices
are related by a $U(4)$ transformation.

The valley polarized states $|\text{VP}\rangle$ correspond to choosing $U=1$
and $(\eta_{1}$, $s_{1})=(K$, $\uparrow)$, $(\eta_{2}$, $s_{2})=(K$,
$\downarrow)$, $(\eta_{3}$, $s_{3})=(K^{\prime}$, $\uparrow)$ above, such that
the total valley charge is maximized.

For $\Phi=0$, the Kramers intervalley coherent (KIVC) state belongs to the
$U(4)$ manifold, and is related to the valley polarized (VP) state by a
$C_{2z}P$ angle of $\pi/4$

$\displaystyle|\text{KIVC}\rangle$
$\displaystyle=\text{exp}\bigg{(}i\frac{\pi}{4}\mathcal{S}\bigg{)}|\text{VP}\rangle,$
(23)

with $\mathcal{S}$ the generator of the particle-hole $U(1)$ in the projected
system,

$\displaystyle\mathcal{S}=$
$\displaystyle\sum_{\boldsymbol{k}}\sum_{\rho\rho^{\prime}}c^{\dagger}_{\boldsymbol{k}\rho}[\tau_{y}\lambda_{z}s_{0}]_{\rho\rho^{\prime}}c_{\boldsymbol{k}\rho^{\prime}},$
(24)

again with the multi-index $\rho$ denoting valley, irrep and spin. The unitary
$\text{exp}\big{(}i\phi\mathcal{S}\big{)}$ transforms the basis as

$\displaystyle
e^{i\phi\mathcal{S}}d^{\dagger}_{\boldsymbol{k}\rho}e^{-i\phi\mathcal{S}}=\bigg{[}\text{exp}\Big{(}i\phi(\tau_{y}\lambda_{z}s_{0})^{T}\Big{)}\bigg{]}_{\rho\rho^{\prime}}d^{\dagger}_{\boldsymbol{k}\rho^{\prime}}.$
(25)

On another hand, he chiral $U(4)\times U(4)$ group contain operations that
rotate each Chern sector independently[41]. The time reversal intervalley
coherent (TIVC) order relates to VP via a $U(4)\times U(4)$ rotation with
$C_{2z}P$ angles of $+\pi/4$ and $-\pi/4$,

$\displaystyle|\text{TIVC}\rangle=\text{exp}\bigg{(}i\frac{\pi}{4}\mathcal{S}_{+1}-i\frac{\pi}{4}\mathcal{S}_{-1}\bigg{)}|\text{VP}\rangle.$
(26)

Here $\mathcal{S}_{\lambda}=P_{\lambda}\mathcal{S}P_{\lambda}$, with
$P_{\lambda}$ the projector onto irrep $\lambda$.

In addition, particularizing to $\nu=0$ we consider the fully spin polarized
(SP) state,

$\displaystyle|\text{SP}\rangle$
$\displaystyle=\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}K+1\uparrow}d^{\dagger}_{\boldsymbol{k}K-1\uparrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}+1\uparrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}-1\uparrow}|0\rangle,$
(27)

that is derived from the valley-spin polarized (VSP) state,

$\displaystyle|\text{VSP}\rangle$
$\displaystyle=\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}K+1\uparrow}d^{\dagger}_{\boldsymbol{k}K-1\uparrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}+1\downarrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}-1\downarrow}|0\rangle,$
(28)

after a spin rotation in valley $K^{\prime}$, belonging to $SU(2)_{K}\times
SU(2)_{K^{\prime}}$.

$\displaystyle|\text{SP}\rangle=\text{exp}\bigg{(}i\frac{\pi}{2}P_{K^{\prime}}s_{y}P_{K^{\prime}}\bigg{)}|\text{VSP}\rangle,$
(29)

with $P_{K^{\prime}}$ the projector onto valley $K^{\prime}$. At zero field,
by $\mathcal{T}$ symmetry the VP and the VSP states have the same energy.

At $\Phi=\Phi_{0}$, the intervalley coherent state corresponds to a $C_{2z}P$
rotation of the VP, of angle $\pi/4$,

$\displaystyle|\text{IVC}\rangle=$
$\displaystyle\text{exp}\bigg{(}i\frac{\pi}{4}\mathcal{S}\bigg{)}|\text{VP}\rangle,$
(30)

with $\mathcal{S}$ taking a different form in accordance with our gauge choice
of Eq. 18,

$\displaystyle\mathcal{S}=\sum_{\boldsymbol{k}}\sum_{\rho\rho^{\prime}}c^{\dagger}_{\boldsymbol{k}\rho}[\lambda_{0}\tau_{x}s_{0}]_{\rho\rho^{\prime}}c_{\boldsymbol{k}\rho^{\prime}}.$
(31)

### Explicit symmetry breaking in the ground state manifold

Figure 7: Breaking of the chiral $\boldsymbol{U(4)\times U(4)}$ symmetry.
Coulomb energy ($\epsilon$ is set to $1$) per unit cell of candidates states
in the $U(4)\times U(4)$ manifold at filling $\nu=-2$, which is explicitly
broken to $U(4)$ with $\phi_{+1}=\phi_{-1}$. The zero point of the energy is
arbitrary. VP, KIVC and TIVC are the valley polarized, Kramers intervalley
coherent and time reversal intervalley coherent states, respectively. Sending
$\phi_{\pm 1}\to-\phi_{\pm 1}$ amounts to a $U(1)_{v}$ transformation of angle
$\pi$. $U(1)_{v}$ is visibly preserved.

The strength of the explicit symmetry breaking processes can be measured by
the difference in Coulomb energy of the different ground state candidates. In
Tables F2 and F1 of Appendix F we tabulate the different contributions to the
total energy (Hartree, Fock, kinetic and Hubbard energy) of the states
discussed above. We note that, remarkably, the kinetic energy is degenerate
for all the states considered at a a given filling. Here we focus on the
states at charge neutrality.

At $\Phi=0$, p-h breaking corresponds to the difference between VP and KIVC
energies,

$\displaystyle\epsilon\big{(}\langle V\rangle_{\text{VP}}-\langle
V\rangle_{\text{KIVC}}\big{)}/N_{M}=2.26\ \text{meV},$ (32)

the breaking of the chiral symmetry by the difference between VP and TIVC
states,

$\displaystyle\epsilon\big{(}\langle V\rangle_{\text{TIVC}}-\langle
V\rangle_{\text{VP}}\big{)}/N_{M}=14.01\ \text{meV},$ (33)

and the breaking of $SU(2)_{K}\times SU(2)_{K^{\prime}}$ by the energies of
VSP (same energy as VP) and SP,

$\displaystyle\epsilon\big{(}\langle V\rangle_{\text{SP}}-\langle
V\rangle_{\text{VP}}\big{)}/N_{M}=10.32\ \text{meV}.$ (34)

Here $N_{M}$ is the number of unit cells, and we have multiplied by $\epsilon$
so that we are comparing energies per unit cell for $\epsilon=1$.

Very similarly for $\Phi=\Phi_{0}$, p-h breaking corresponds to the difference
between VP and IVC energies,

$\displaystyle\epsilon\big{(}\langle V\rangle_{\text{IVC}}-\langle
V\rangle_{\text{VP}}\big{)}/N_{M}=0.74\ \text{meV},$ (35)

and the breaking of $SU(2)_{K}\times SU(2)_{K^{\prime}}$ by the VSP and SP
states,

$\displaystyle\epsilon\big{(}\langle V\rangle_{\text{SP}}-\langle
V\rangle_{\text{VSP}}\big{)}/N_{M}=5.23\ \text{meV}.$ (36)

Additionally, in Fig. 7 we plot the Coulomb energy of states of the form
$\text{exp}\big{(}i\phi_{+1}\mathcal{S}_{+1}+i\phi_{-1}\mathcal{S}_{-1}\big{)}|\text{VP}\rangle$
for $\nu=-2$ and $B=0$ T. The breaking of the chiral $U(4)\times U(4)$
symmetry down to $U(4)$ is evident, and the inversion symmetry of the plot
shows the conservation of $U(1)_{v}$.

In light of the results, we conclude that $C_{2z}P$ breaking is smaller than
the breaking of the chiral symmetry at $\Phi=0$ and of $SU(2)_{K}\times
SU(2)_{K^{\prime}}$.

The $SU(2)_{K}\times SU(2)_{K^{\prime}}$ breaking is stronger and comparable
to the breaking of the chiral symmetry at zero field. This effect was called
intervalley Hund’s interaction in Ref. [21], where it was argued that it is
the smallest energy scale in the hierarchy of symmetry breakings. Our
calculations in the tight-binding model involve lattice-scale interactions,
providing reliable values for the splitting of the $SU(2)_{K}\times
SU(2)_{K^{\prime}}$ multiplets[68].

The Hubbard interaction clearly breaks $SU(2)_{K}\times SU(2)_{K^{\prime}}$,
and favours the states with a net spin polarization. On the other hand, we
showed that the long ranged Coulomb energy is Hund antiferromagnetic,
supporting states with opposite spins in different valleys. Although in the
real system we expect their values to be correlated, the interplay between $U$
and $\epsilon$ (and possibly other effects due to phonons[68]) determines the
sign of the Hund’s coupling $J$. For instance, the splitting between the VSP
and SP states changes sign when $\epsilon U=8.82$ eV($9.02$ eV) at zero
field(one flux quantum). Also, if the Hubbard interaction is strong enough, it
can go beyond selecting the state of the $SU(2)_{K}\times SU(2)_{K^{\prime}}$
multiplets and stabilize spin polarized phases, as we will see.

Finally, notice that the small $C_{2z}P$ breaking favours the KIVC phase as
can be seen in Fig. 7 or Eq. 32. However, the gain in energy of the ’dressed’
self-consistent states is larger than this small splitting between the VP and
KIVC, so the p-h breaking is not the decisive factor in the stability of the
ground states.

##  V SELF-CONSISTENT HARTREE-FOCK

We have carried out self-consistent Hartree-Fock simulations in a system of
$12\times 12$ unit cells, focusing on filling factors $\nu=-2,0,+2$. We
describe the Hartree-Fock formalism and the flat band projection method in
Appendix D.

Typical values for $\epsilon$ found in the literature range from about $7$ to
$12$[26, 21], so we choose $\epsilon=10$ and a realistic value for $U$ of $4$
eV[27, 42]. However, it has been argued that internal screening is large in
these systems and a more appropriate value for $\epsilon$ is several times
larger[27, 69]. This agrees with the fact that lower values of $\epsilon$
overestimate the gap of the insulators, which in transport are found to be
$\lesssim$ 1 meV[4, 3, 5, 52, 38]. We account for both scenarios and report
results also for $\epsilon=50$ and $0.5$ eV.

The self-consistent states are characterized by the $Q$ matrix, defined by

$\displaystyle[Q(\boldsymbol{k})]_{\rho\rho^{\prime}}$
$\displaystyle=2[P(\boldsymbol{k})]_{\rho\rho^{\prime}}-\delta_{\rho\rho^{\prime}},$
$\displaystyle[P(\boldsymbol{k})]_{\rho\rho^{\prime}}$ $\displaystyle=\langle
d^{\dagger}_{\boldsymbol{k}\rho}d_{\boldsymbol{k}\rho^{\prime}}\rangle,$ (37)

with the properties $Q(\boldsymbol{k})=Q(\boldsymbol{k})^{\dagger}$,
$Q(\boldsymbol{k})^{2}=1$ and $\text{tr}(Q(\boldsymbol{k}))=2\nu$. In most
cases, as we discuss below, $Q$ will be diagonal in the spin,
$Q(\boldsymbol{k})=Q_{\uparrow}(\boldsymbol{k})P_{\uparrow}+Q_{\downarrow}(\boldsymbol{k})P_{\downarrow}$
($P_{\uparrow(\downarrow)}$ is the projector onto spin
$\uparrow(\downarrow)$), with each spin polarization either completely empty
($Q_{s}(\boldsymbol{k})=-1$), completely full ($Q_{s}(\boldsymbol{k})=+1$) or
half filled. If half filled, $Q_{s}$ can be expressed as a linear combination
of products of Pauli matrices,

$\displaystyle
Q_{s}(\boldsymbol{k})=\mathop{\sum_{\alpha,\beta=0,x,y,z}}_{(\alpha,\beta)\neq(0,0)}A^{s}_{\alpha\beta}(\boldsymbol{k})\sigma_{\alpha}\tau_{\beta}=\mathop{\sum_{\alpha,\beta=0,x,y,z}}_{(\alpha,\beta)\neq(0,0)}B^{s}_{\alpha\beta}(\boldsymbol{k})\lambda_{\alpha}\tau_{\beta},$
(38)

with $A^{s}_{\alpha,\beta}(\boldsymbol{k})$,
$B^{s}_{\alpha,\beta}(\boldsymbol{k})$ real coefficients and
$\sum_{\alpha\beta}(A^{s}_{\alpha\beta}(\boldsymbol{k}))^{2}=\sum_{\alpha\beta}(B^{s}_{\alpha\beta}(\boldsymbol{k}))^{2}=1$.
In the following we will write $\langle\sigma_{\alpha}\tau_{\beta}\rangle$,
$\langle\lambda_{\alpha}\tau_{\beta}\rangle$ to denote the coefficients
$A^{s}_{\alpha\beta}(\boldsymbol{k})$, $B^{s}_{\alpha\beta}(\boldsymbol{k})$.
The momentum dependence is left implicit, and the spin can be deduced
depending on the context.

The preferred state of the analytical approaches and numerical studies at zero
magnetic field is the KIVC[21, 34, 22, 36]. At $\nu=0$ it can be either spin
singlet, $Q(\boldsymbol{k})=\sigma_{y}\tau_{y}$, or ’spin triplet’,
$Q(\boldsymbol{k})=\sigma_{y}\tau_{y}\boldsymbol{n}\cdot\boldsymbol{s}$, with
$\boldsymbol{n}$ denoting an spontaneous direction. An antiferromagnetic
Hund’s coupling, $J<0$, favours the KIVC singlet whereas $J>0$ prefers the
’triplet’[21]. However, in our Hartree Fock numerics we restrict the state to
be a direct product of spin up and spin down wave functions, so $Q$ will be
diagonal in the spin index and we cannot access the ’triplet’ state. Setting
$\boldsymbol{n}=\boldsymbol{\hat{z}}$ amounts to a valley rotation which does
not change the kinetic or Coulomb energy. By a similar argument to the one
given in Appendix D, the Hubbard energy does not change either.

a) b) phantom

Figure 8: KIVC order at $\boldsymbol{\Phi=0}$. a) $\Delta_{\text{KIVC}}$ at
$\nu=+2$ and $\epsilon=10,U=4$ eV. The inter-Chern parameter completes the
total weight, $\Delta_{\text{IC}}=1-(\Delta_{\text{KIVC}})^{2}$, reaching a
value of $0.95$ at $\Gamma_{M}$. In the case of $\epsilon=50,U=0.5$ eV the
only difference is that $\Delta_{\text{IC}}$ is smaller with a maximal value
of $0.46$ at $\Gamma_{M}$. For $\nu=0$ the order parameter is almost
identical. b) KIVC order parameter at $\nu=-2$ and $\epsilon=50,\ U=0.5$ eV.
Around $\Gamma_{M}$ there are extra holes creating a Fermi pocket, and the
order parameter is ill-defined in that region.

For $\nu=\pm 2$, a ferromagnetic Hund’s coupling leads to
$Q(\boldsymbol{k})=\sigma_{y}\tau_{y}P_{\uparrow}\pm P_{\downarrow}$ where the
spins in both valleys are aligned. Antiferromagnetic coupling on the other
hand promotes the state where the spins of the two valleys are anti-aligned.
Again, we can only detect the spin-diagonal order.

Notice that the $U(1)_{v}$ symmetry allows for an arbitrary global rotation in
the order parameter, $\tau_{y}\to\cos(\theta)\tau_{y}+\sin(\theta)\tau_{x}$.
However, the difference in the intervalley angle for different values of
$\boldsymbol{k}$ cannot be removed and is physical. In any case, we do not
observe textures or windings in the IVC angle.

Contrary to the the BM theory[42], in the tight-binding model the on-site
Hubbard interaction is implemented naturally. At charge neutrality the SP
state is not competitive if we consider only the Coulomb energy, but for
sufficiently large $U$ it will be the lowest energy state. For electron and
hole dopings, the candidate states are spin polarized, see Table F2 in
Apppendix F, and the preferred states (KIVC, VP) are degenerate in $H_{U}$, so
we do not expect different orders for different $U$.

At flux $\Phi_{0}$, analytical studies[51] suggest that the Zeeman energy
drives the system to maximize the spin polarization in the $U(4)$ manifold. At
charge neutrality the SP state with $Q(\boldsymbol{k})=s_{z}$ is the ground
state and at $\nu=-2(+2)$ there are two possible orders, IVC, with
$Q(\boldsymbol{k})=\lambda_{0}\tau_{y}P_{\uparrow}-P_{\downarrow}(\lambda_{0}\tau_{y}P_{\downarrow}+P_{\uparrow})$,
or valley polarization,
$Q(\boldsymbol{k})=\lambda_{0}\tau_{z}P_{\uparrow}-P_{\downarrow}(\lambda_{0}\tau_{z}P_{\downarrow}+P_{\uparrow})$.

$\epsilon=10,\ U=4$ eV

$\epsilon=50,\ U=0.5$ eV

Figure 9: KIVC band structures of MATBG at $\boldsymbol{B=0}$ T. Different
spins are shown in different colors. From left to right, the filling is
$\nu=-2,0,+2$. At $\nu=\pm 2$ these are the true ground states, whereas at
charge neutrality there is competition between the KIVC and the spin polarized
state. The many-body electron-hole asymmetry is apparent.

### Results for $\boldsymbol{\Phi=0}$

Examining Table F2 we deduce that the possible orders for $B=0$ T are the VP
or KIVC, or the SP at charge neutrality. We compute self-consistent states
with initial guesses for the KIVC, VP and SP orders, whose energies are
tabulated in Table F4 of Appendix F. We plot the band structures of the KIVC
states in Fig. 9, of the SP states in Fig. 10 and of the VP states in Fig. E.6
of Appendix E. The KIVC, with order parameter

$\displaystyle\Delta_{\text{KIVC}}=\langle\sigma_{y}\tau_{y}\rangle,$ (39)

is the ground state for $\nu=\pm 2$, and the KIVC and SP are competitive at
charge neutrality.

$\epsilon=10,\ U=4$ eV

$\epsilon=50,\ U=0.5$ eV

Figure 10: Band structures of the spin polarized phase at charge neutrality.
Different spins are shown with different colors. The spin polarization is not
full, and the system is a metal.

A general feature of the self-consistent states is that deviations from the
predominant order parameter occur near the $\Gamma_{M}$ point.

In the gapped phases they involve the inter-Chern parameter,

$\displaystyle\Delta_{\text{IC}}=\langle\sigma_{x}\tau_{0}\rangle^{2}+\langle\sigma_{y}\tau_{z}\rangle^{2},$
(40)

signalling coherent superpositions of bands with different Chern numbers.
Notice how the expression for $\Delta_{\text{IC}}$ is invariant with respect
to the gauge ambiguity of Eq. 17. The effect of $\Delta_{\text{IC}}$ near
$\Gamma_{M}$ is to reduce the exchange energy at the expense of kinetic (and
Hartree for $\nu=\pm 2$) energy. In consequence, $\Delta_{\text{IC}}$ becomes
larger for lower $\epsilon$, where the kinetic energy penalty is comparatively
less important.

For $\nu=-2$ the system is metallic with a hole pocket around $\Gamma_{M}$.
Because the states near $\Gamma_{M}$ show a homogeneous density profile, the
Hartree energy is optimized if the excess density relative to charge
neutrality originates from this region of the Brillouin zone. This is
manifested by the characteristic dips at electron doping, or rises at hole
doping in the band structures[70]. At $\nu=-2$ the holes are populated,
optimizing the Hartree energy.

In Fig. 8 we plot the KIVC order parameter for two different cases, where both
phenomena are illustrated. Additionally, as depicted in Fig. 10 the bands of
the spin polarized states show also the Fermi pockets at $\Gamma_{M}$, and the
additional electrons and holes appear in inter-Chern coherence.

The different behaviour at $\nu=+2$ and $\nu=-2$ is a sign of the many-body
electron-hole asymmetry, which is a robust experimental feature[38, 4, 3]. We
argue that the flat band limit shows electron-hole symmetry but the kinetic
energy disturbs it at weaker couplings. Indeed, the $H_{0}$ spectrum is very
asymmetric at zero flux, see Fig. 3. At $\Phi_{0}$ flux the spectrum is
visibly more p-h symmetric, and the electron-hole relation between $\nu=+2$
and $-2$ is more apparent for $\epsilon=10$, $U=4$ eV (Fig. 13).

Let us now discuss the energies of the different self-consistent phases at
charge neutrality. In the tables below we write the energies of the KIVC, SP
and VP self-consistent states at relative to the pure KIVC state, in meV per
unit cell (the value of $U$ is in units of eV).

$\epsilon=10$, $U=4$ eV
$\epsilon\times$Hartree $\epsilon\times$Fock Hubb./$U$ kinetic total KIVC
$-0.49$ $-64.60$ $0.01$ $1.28$ $-5.19$ SP $-0.13$ $-8.19$ $-1.10$ $0.61$
$-4.62$ VP $-0.13$ $-17.97$ $0.00$ $0.61$ $-1.20$

$\epsilon=50$, $U=0.5$ eV
$\epsilon\times$Hartree $\epsilon\times$Fock Hubb./$U$ kinetic total KIVC
$-0.91$ $-31.14$ $0.00$ $0.40$ $-0.24$ SP $-0.67$ $18.61$ $-1.13$ $-0.16$
$-0.37$ VP $-0.67$ $8.48$ $0.00$ $-0.16$ $-0.00$

The KIVC always has lower energy than the VP state, but the competition
between KIVC and SP is more complex. For $\epsilon=10$, $U=4$ eV the KIVC
state is $0.57$ meV per unit cell lower than the SP state, whereas for
$\epsilon=50$, $U=0.5$ eV the SP state is the ground state with $0.13$
meV/cell of difference. Assuming that the self-consistent state does not
change much with $U$ (as a matter of fact, we find that the states for
$\epsilon=10,U=4$ are almost identical to those with $\epsilon=10,U=0.5$ eV),
we expect a transition to the SP state for $U$ grater that a critical value of
$U_{c}=4.51$ eV when $\epsilon=10$, and to the KIVC for $U$ smaller than
$U_{c}=0.39$ when $\epsilon=50$.

### Results for $\boldsymbol{\Phi=\Phi_{0}}$

The energies of self-consistent states for $\Phi=\Phi_{0}$ are tabulated in
Table F3 of Appendix F, and their band structures plotted in Fig. 13. The spin
polarization is maximal for all fillings, however, we observe gapped phases
with a completely different order compared to the previously predicted for
$\nu=\pm 2$.

The order parameter

$\displaystyle\Delta_{xyz,0}=$
$\displaystyle\langle\lambda_{x}\tau_{0}\rangle^{2}+\langle\lambda_{y}\tau_{0}\rangle^{2}+\langle\lambda_{z}\tau_{0}\rangle^{2},$
(41)

which is invariant under the gauge ambiguity of Eq. 20, is predominant at hole
doping. $\Delta_{xyz,0}$ corresponds to choosing a particular irrep basis
allowed by Eq. 20, and filling the bands with the same irrep number of the two
valleys.

At electron doping, the state for $\epsilon=10$, $U=4$ eV is very similar to
the many-body particle-hole partner of the $\nu=-2$ insulator, with the same
order parameter. On the other hand, when $\epsilon$ is larger than some
critical value between $30$ and $50$, see Appendix E, there is a phase
transition that promotes the IVC state. For both $\epsilon=10$, $U=4$ eV and
$\epsilon=50$, $U=0.5$ eV the second largest order parameter is found to be
$\Delta_{xyz,z}$, which is obtained after replacing $\tau_{0}$ by $\tau_{z}$
in Eq. 41. The total ’spectral weight’ is essentially saturated by the two
terms.

Figure 11: The order parameter for $\boldsymbol{\nu=-2}$ at
$\boldsymbol{\Phi=\Phi_{0}}$. The dominant order parameter of the self-
consistent states, $\Delta_{xyz,0}$, is shown. For $\epsilon=10$ it is very
similar to its many-body partner at $\nu=+2$ in Fig. 12. The $\Delta_{xyz,0}$
character is enhanced for weaker couplings.

Figure 12: Phase transition at $\boldsymbol{\nu=+2}$ and
$\boldsymbol{\Phi=\Phi_{0}}$. The order parameters of the self-consistent
states are shown for $\epsilon=10$, $U=4$ eV and $\epsilon=50$, $U=0.5$ eV.
For some value of $\epsilon$ between $50$ and $30$ there is a transition from
the IVC to the $\Delta_{xyz,0}$ state.

Let us discuss the discrepancies between our numerical results and the results
of Ref. [51]. Firstly, the kinetic spectrum at $B=26.5$ T is much broader than
at zero field, with a gap between the valence and conduction bands of about
$4$ meV. The system is in an intermediate coupling regime as the kinetic
energy is comparable to the interaction energy. Secondly, at $B=0$ T the order
is exchange-driven (at least at fillings $\nu=\pm 2$ when the Hubbard term is
not relevant), and the Fock energy is optimized in the ground state. For
$B=26.5$ T we find that the direct term is more influential. In Ref. [51] the
’flat metric condition’ assumes that the Hartree energy is trivial, and the
strong coupling analysis suggests that the kinetic energy is a subdominant
scale. Both premises are at odds with our numerical results in the tight-
binding model.

The ground state at hole doping has Chern number -2 and is adiabatically
connected to the ground state of the non interacting Hamiltonian. The
predominant order parameter, $\Delta_{xyz,0}$, is plotted in Fig. 11. In the
following table we write the energy gain of the $\nu=-2$ ground states with
respect to the candidate IVC state for $\epsilon=10$, $U=4$ eV and
$\epsilon=50$, $U=0.5$ eV, in meV per unit cell ($U$ is expressed in eV).

$\nu=-2$
$\epsilon,\ U$ $\epsilon\times$Hartree $\epsilon\times$Fock Hubb./$U$ kinetic
total $10,\ 4$ $-100.27$ $-41.38$ $-0.13$ $-2.89$ $-17.58$ $50,\ 0.5$
$-105.80$ $19.14$ $-0.15$ $-4.94$ $-6.75$

Clearly, the Hartree and kinetic energies drive the system towards a different
order to the previously proposed.

On the other hand, we write in the table below the energy differences of the
self-consistent states for $\nu=+2$ and the pure IVC state.

$\nu=+2$
$\epsilon,\ U$ $\epsilon\times$Hartree $\epsilon\times$Fock Hubb./$U$ kinetic
total $10,\ 4$ $-155.20$ $-3.40$ $-0.26$ $5.21$ $-11.69$ $50,\ 0.5$ $-1.65$
$-9.29$ $0.00$ $-0.19$ $-0.41$

In this case, we conclude that the interplay between the Hartree and kinetic
energies determines the transition between the two competing orders, with
their order parameters plotted in Fig. 12. The phase transition at
intermediate values of $\epsilon$ is topological, since the state for
$\epsilon=10$, $U=4$ eV is found to have Chern number $+2$, and the IVC state
is trivial with Chern number $0$.

$\epsilon=10,\ U=4$ eV

$\epsilon=50,\ U=0.5$ eV

Figure 13: Quasiparticle band structures of MATBG at $\boldsymbol{B=26.5}$ T.
Spin up bands are shown in blue and spin down bands in red. From left to
right, the filling factor is $\nu=-2,0,+2$. The states are metallic at charge
neutrality and insulating for electron and hole doping. When $\epsilon=10,U=4$
eV the $\nu=\pm 2$ states are related by a particle-hole transformation.

##  VI CONCLUSIONS

In this work, we have studied magic angle twisted bilayer graphene under
external magnetic fluxes of zero and one flux quantum per unit cell, and
dopings of $-2$, $0$ and $+2$ electrons per unit cell. In the atomistic tight-
binding model we capture details of the physics that are not available in
continuum theories.

For instance, we have established numerically and analytically that the long
ranged Coulomb interaction contributes to an antiferromagnetic intervalley
Hund’s coupling. The on-site Hubbard term, which is generally not included in
the continuum model[42], contributes to a ferromagnetic coupling. The
character of these valley exchanging interactions is important for the
superconducting order parameter in Moiré systems[71, 72, 73].

In addition, we studied the particle-hole asymmetry of the Hilbert space of
the flat bands, and concluded that p-h breaking effects represent a small
energy scale in the problem. In particular, the effects of p-h breaking cannot
stabilize the time reversal intervalley coherent (TIVC) state at $\nu=-2$
observed recently[40], and it must rely on other mechanisms like electron-
phonon coupling[36, 74].

In order to elucidate the nature of the ground state we have performed self-
consistent Hartree-Fock simulations. For $\Phi=0$, the Kramers intervalley
coherent state is the ground state at electron and hole doping, and is
competitive with the spin polarized state, supported by the Hubbard energy, at
charge neutrality. In contrast to the BM model, where the self-consistent KIVC
is favoured by kinetic energy superexchange[20, 32, 34, 21, 36, 22], in the
tight-binding model the inter-Chern order at $\Gamma_{M}$ lowers the exchange
at the expense of kinetic and Hartree energies. Also, in the BM model the on-
site Hubbard interaction is generally not included[42], and the SP state is
never a candidate. We predict a phase transition from the KIVC to the SP at a
critical value of $U_{c}=4.51$ eV when $\epsilon=10$. Finally, we also detect
the many-body electron-hole asymmetry found in the experiments[4, 38, 3].

For $\Phi_{0}$ flux, the self-consistent states are spin polarized due to the
Zeeman effect. Also, the Dirac cones are gapped and the bandwidth of the flat
bands is increased, driving the system to an intermediate coupling regime.
There are correlated insulators for $\nu=\pm 2$, and at charge neutrality the
state is (almost) fully spin polarized. For electron doping and small
screening, we observe an insulator with Chern number $+2$, consistent with
experiments[52], and a trivial insulator for large screening. We envision the
possibility of observing a topological phase transition by changing the
screening environment in the experimental setup.

Strain[36, 39] and electron-phonon coupling[22, 43, 74] are relevant to the
physics and their treatment is left for future work. Also, an accurate
understanding of internal screening is needed, and it has been suggested that
it is large in MATBG[75, 27]. Improvements of the Hartree-Fock method such as
the inclusion of more bands in the projection, and specially allowing for
coherence between spins, like in the KIVC ’triplet’ state, will further refine
the outcome.

Finally, our results shed light into the correlated states in magic angle
twisted bilayer graphene with precise atomic-scale modelling, both at $B=0$
and $26.5$ T. Furthermore, our calculations constitute the first self-
consistent study of the correlated phases in MATBG under one magnetic flux
quantum.

##  VII Acknowledgements

This work has been supported by MICINN (Spain) under Grant No.
PID2020-113164GBI00, as well as by the CSIC Research Platform on Quantum
Technologies PTI-001. The access to computational resources of CESGA (Centro
de Supercomputación de Galicia) is also gratefully acknowledged.

## References

  * Polini  _et al._ [2022] M. Polini, F. Giazotto, K. C. Fong, I. M. Pop, C. Schuck, T. Boccali, G. Signorelli, M. D’Elia, R. H. Hadfield, V. Giovannetti, D. Rossini, A. Tredicucci, D. K. Efetov, F. H. L. Koppens, P. Jarillo-Herrero, A. Grassellino, and D. Pisignano, Materials and devices for fundamental quantum science and quantum technologies (2022).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018a).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019), https://www.science.org/doi/pdf/10.1126/science.aav1910 .
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Liu  _et al._ [2021] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. I. A. Li, Tuning electron correlation in magic-angle twisted bilayer graphene using coulomb screening, Science 371, 1261 (2021), https://www.science.org/doi/pdf/10.1126/science.abb8754 .
  * Jaoui  _et al._ [2022] A. Jaoui, I. Das, G. Di Battista, J. Díez-Mérida, X. Lu, K. Watanabe, T. Taniguchi, H. Ishizuka, L. Levitov, and D. K. Efetov, Quantum critical behaviour in magic-angle twisted bilayer graphene, Nature Physics 18, 633 (2022).
  * Cao  _et al._ [2020] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Strange metal in magic-angle graphene with near planckian dissipation, Phys. Rev. Lett. 124, 076801 (2020).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Chern insulators, van hove singularities and topological flat bands in magic-angle twisted bilayer graphene, Nature Materials 20, 488 (2021).
  * Stepanov  _et al._ [2021] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Competing zero-field chern insulators in superconducting twisted bilayer graphene, Phys. Rev. Lett. 127, 197701 (2021).
  * Xie  _et al._ [2021] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Fractional chern insulators in magic-angle twisted bilayer graphene, Nature 600, 439 (2021).
  * Repellin and Senthil [2019] C. Repellin and T. Senthil, Chern bands of twisted bilayer graphene: fractional chern insulators and spin phase transition (2019), arXiv:1912.11469 [cond-mat.str-el] .
  * Ledwith  _et al._ [2020] P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath, Fractional chern insulator states in twisted bilayer graphene: An analytical approach, Phys. Rev. Res. 2, 023237 (2020).
  * Chew  _et al._ [2023] A. Chew, Y. Wang, B. A. Bernevig, and Z.-D. Song, Higher-order topological superconductivity in twisted bilayer graphene, Phys. Rev. B 107, 094512 (2023).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Correlated electronic phases in twisted bilayer transition metal dichalcogenides, Nature Materials 19, 861 (2020).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249 (2021).
  * Scheer and Lian [2023] M. G. Scheer and B. Lian, Twistronics of kekulé graphene: Honeycomb and kagome flat bands (2023), arXiv:2305.19927 [cond-mat.mes-hall] .
  * Crépel  _et al._ [2023] V. Crépel, A. Dunbrack, D. Guerci, J. Bonini, and J. Cano, Chiral model of twisted bilayer graphene realized in a monolayer (2023), arXiv:2305.14423 [cond-mat.mes-hall] .
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018b).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of mott insulating behavior and superconductivity in twisted bilayer graphene, Phys. Rev. X 8, 031089 (2018).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong coupling phases of partially filled twisted bilayer graphene narrow bands, Phys. Rev. Lett. 122, 246401 (2019).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Phys. Rev. X 10, 031034 (2020).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Kekulé spiral order at all nonzero integer fillings in twisted bilayer graphene, Phys. Rev. X 11, 041063 (2021).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Untying the insulating and superconducting orders in magic-angle graphene, Nature 583, 375 (2020).
  * Sharpe  _et al._ [2021] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Evidence of orbital ferromagnetism in twisted bilayer graphene aligned to hexagonal boron nitride, Nano Letters 21, 4299 (2021), pMID: 33970644, https://doi.org/10.1021/acs.nanolett.1c00696 .
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Zhang  _et al._ [2022] S. Zhang, X. Lu, and J. Liu, Correlated insulators, density wave states, and their nonlinear optical response in magic-angle twisted bilayer graphene, Phys. Rev. Lett. 128, 247402 (2022).
  * González and Stauber [2021] J. González and T. Stauber, Magnetic phases from competing hubbard and extended coulomb interactions in twisted bilayer graphene, Phys. Rev. B 104, 115110 (2021).
  * González and Stauber [2020] J. González and T. Stauber, Time-reversal symmetry breaking versus chiral symmetry breaking in twisted bilayer graphene, Phys. Rev. B 102, 081118 (2020).
  * Klebl  _et al._ [2021] L. Klebl, Z. A. H. Goodwin, A. A. Mostofi, D. M. Kennes, and J. Lischner, Importance of long-ranged electron-electron interactions for the magnetic phase diagram of twisted bilayer graphene, Phys. Rev. B 103, 195127 (2021).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. 124, 097601 (2020).
  * Faulstich  _et al._ [2023] F. M. Faulstich, K. D. Stubbs, Q. Zhu, T. Soejima, R. Dilip, H. Zhai, R. Kim, M. P. Zaletel, G. K.-L. Chan, and L. Lin, Interacting models for twisted bilayer graphene: A quantum chemistry approach, Phys. Rev. B 107, 235123 (2023).
  * Vafek and Kang [2020] O. Vafek and J. Kang, Renormalization group study of hidden symmetry in twisted bilayer graphene with coulomb interactions, Phys. Rev. Lett. 125, 257602 (2020).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. iii. interacting hamiltonian and exact symmetries, Phys. Rev. B 103, 205413 (2021).
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Twisted bilayer graphene. iv. exact insulator ground states and phase diagram, Phys. Rev. B 103, 205414 (2021).
  * Seo  _et al._ [2019] K. Seo, V. N. Kotov, and B. Uchoa, Ferromagnetic mott state in twisted graphene bilayers at the magic angle, Phys. Rev. Lett. 122, 246402 (2019).
  * Kwan  _et al._ [2023] Y. H. Kwan, G. Wagner, N. Bultinck, S. H. Simon, E. Berg, and S. A. Parameswaran, Electron-phonon coupling and competing kekulé orders in twisted bilayer graphene (2023), arXiv:2303.13602 [cond-mat.str-el] .
  * Ledwith  _et al._ [2021] P. J. Ledwith, E. Khalaf, and A. Vishwanath, Strong coupling theory of magic-angle graphene: A pedagogical introduction, Annals of Physics 435, 168646 (2021), special issue on Philip W. Anderson.
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Unconventional sequence of correlated chern insulators in magic-angle twisted bilayer graphene, Nature Physics 17, 1210 (2021).
  * Parker  _et al._ [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Strain-induced quantum phase transitions in magic-angle graphene, Phys. Rev. Lett. 127, 027601 (2021).
  * Nuckolls  _et al._ [2023] K. P. Nuckolls, R. L. Lee, M. Oh, D. Wong, T. Soejima, J. P. Hong, D. Călugăru, J. Herzog-Arbeitman, B. A. Bernevig, K. Watanabe, T. Taniguchi, N. Regnault, M. P. Zaletel, and A. Yazdani, Quantum textures of the many-body wavefunctions in magic-angle graphene (2023), arXiv:2303.00024 [cond-mat.mes-hall] .
  * Călugăru  _et al._ [2022] D. Călugăru, N. Regnault, M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, A. Yazdani, O. Vafek, and B. A. Bernevig, Spectroscopy of twisted bilayer graphene correlated insulators, Phys. Rev. Lett. 129, 117602 (2022).
  * Jimeno-Pozo  _et al._ [2023] A. Jimeno-Pozo, Z. A. H. Goodwin, P. A. Pantaleón, V. Vitale, L. Klebl, D. M. Kennes, A. Mostofi, J. Lischner, and F. Guinea, Short vs. long range exchange interactions in twisted bilayer graphene (2023), arXiv:2303.18025 [cond-mat.mes-hall] .
  * Blason and Fabrizio [2022] A. Blason and M. Fabrizio, Local kekulé distortion turns twisted bilayer graphene into topological mott insulators and superconductors, Phys. Rev. B 106, 235112 (2022).
  * Hofstadter [1976] D. R. Hofstadter, Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields, Phys. Rev. B 14, 2239 (1976).
  * Herzog-Arbeitman  _et al._ [2020] J. Herzog-Arbeitman, Z.-D. Song, N. Regnault, and B. A. Bernevig, Hofstadter topology: Noncrystalline topological materials at high flux, Phys. Rev. Lett. 125, 236804 (2020).
  * Lian  _et al._ [2020] B. Lian, F. Xie, and B. A. Bernevig, Landau level of fragile topology, Phys. Rev. B 102, 041402 (2020).
  * Guan  _et al._ [2022] Y. Guan, O. V. Yazyev, and A. Kruchkov, Reentrant magic-angle phenomena in twisted bilayer graphene in integer magnetic fluxes, Phys. Rev. B 106, L121115 (2022).
  * Singh  _et al._ [2023] K. Singh, A. Chew, J. Herzog-Arbeitman, B. A. Bernevig, and O. Vafek, Topological heavy fermions in magnetic field (2023), arXiv:2305.08171 [cond-mat.str-el] .
  * Wang and Vafek [2022] X. Wang and O. Vafek, Narrow bands in magnetic field and strong-coupling hofstadter spectra, Phys. Rev. B 106, L121111 (2022).
  * Herzog-Arbeitman  _et al._ [2022a] J. Herzog-Arbeitman, A. Chew, and B. A. Bernevig, Magnetic bloch theorem and reentrant flat bands in twisted bilayer graphene at $2\pi$ flux, Phys. Rev. B 106, 085140 (2022a).
  * Herzog-Arbeitman  _et al._ [2022b] J. Herzog-Arbeitman, A. Chew, D. K. Efetov, and B. A. Bernevig, Reentrant correlated insulators in twisted bilayer graphene at 25 t ($2\pi$ flux), Phys. Rev. Lett. 129, 076401 (2022b).
  * Das  _et al._ [2022] I. Das, C. Shen, A. Jaoui, J. Herzog-Arbeitman, A. Chew, C.-W. Cho, K. Watanabe, T. Taniguchi, B. A. Piot, B. A. Bernevig, and D. K. Efetov, Observation of reentrant correlated insulators and interaction-driven fermi-surface reconstructions at one magnetic flux quantum per moiré unit cell in magic-angle twisted bilayer graphene, Phys. Rev. Lett. 128, 217701 (2022).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011), https://www.pnas.org/doi/pdf/10.1073/pnas.1108174108 .
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Phys. Rev. B 86, 155449 (2012).
  * Goodwin  _et al._ [2020] Z. A. H. Goodwin, V. Vitale, X. Liang, A. A. Mostofi, and J. Lischner, Hartree theory calculations of quasiparticle properties in twisted bilayer graphene, Electronic Structure 2, 034001 (2020).
  * Kang and Vafek [2023] J. Kang and O. Vafek, Pseudomagnetic fields, particle-hole asymmetry, and microscopic effective continuum hamiltonians of twisted bilayer graphene, Phys. Rev. B 107, 075408 (2023).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017).
  * Moon and Koshino [2012] P. Moon and M. Koshino, Energy spectrum and quantum hall effect in twisted bilayer graphene, Phys. Rev. B 85, 195458 (2012).
  * Giuliani and Vignale [2005] G. Giuliani and G. Vignale,  _Quantum Theory of the Electron Liquid_ (Cambridge University Press, 2005).
  * Luttinger [1951] J. M. Luttinger, The effect of a magnetic field on electrons in a periodic potential, Phys. Rev. 84, 814 (1951).
  * Nemec and Cuniberti [2007] N. Nemec and G. Cuniberti, Hofstadter butterflies of bilayer graphene, Phys. Rev. B 75, 201404 (2007).
  * Herzog-Arbeitman  _et al._ [2022c] J. Herzog-Arbeitman, Z.-D. Song, L. Elcoro, and B. A. Bernevig, Hofstadter topology with real space invariants and reentrant projective symmetries (2022c).
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. 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, Phys. Rev. X 9, 021013 (2019).
  * Fukui  _et al._ [2005] T. Fukui, Y. Hatsugai, and H. Suzuki, Chern numbers in discretized brillouin zone: Efficient method of computing (spin) hall conductances, Journal of the Physical Society of Japan 74, 1674 (2005), https://doi.org/10.1143/JPSJ.74.1674 .
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Twisted bilayer graphene. ii. stable symmetry anomaly, Phys. Rev. B 103, 205412 (2021).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of magic angles in twisted bilayer graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Liu  _et al._ [2019] J. Liu, J. Liu, and X. Dai, Pseudo landau level representation of twisted bilayer graphene: Band topology and implications on the correlated insulating phase, Phys. Rev. B 99, 155415 (2019).
  * Chatterjee  _et al._ [2020] S. Chatterjee, N. Bultinck, and M. P. Zaletel, Symmetry breaking and skyrmionic transport in twisted bilayer graphene, Phys. Rev. B 101, 165141 (2020).
  * González and Stauber [2023] J. González and T. Stauber, Ising superconductivity induced from spin-selective valley symmetry breaking in twisted trilayer graphene, Nature Communications 14, 2746 (2023).
  * Cea  _et al._ [2022] T. Cea, P. A. Pantaleón, N. R. Walet, and F. Guinea, Electrostatic interactions in twisted bilayer graphene, Nano Materials Science 4, 27 (2022), special issue on Graphene and 2D Alternative Materials.
  * Scheurer and Samajdar [2020] M. S. Scheurer and R. Samajdar, Pairing in graphene-based moiré superlattices, Phys. Rev. Res. 2, 033062 (2020).
  * Chatterjee  _et al._ [2022] S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, Inter-valley coherent order and isospin fluctuation mediated superconductivity in rhombohedral trilayer graphene, Nature Communications 13, 6013 (2022).
  * Lee  _et al._ [2019] J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, and A. Vishwanath, Theory of correlated insulating behaviour and spin-triplet superconductivity in twisted double bilayer graphene, Nature Communications 10, 5333 (2019).
  * Chen  _et al._ [2023] C. Chen, K. P. Nuckolls, S. Ding, W. Miao, D. Wong, M. Oh, R. L. Lee, S. He, C. Peng, D. Pei, Y. Li, S. Zhang, J. Liu, Z. Liu, C. Jozwiak, A. Bostwick, E. Rotenberg, C. Li, X. Han, D. Pan, X. Dai, C. Liu, B. A. Bernevig, Y. Wang, A. Yazdani, and Y. Chen, Strong inter-valley electron-phonon coupling in magic-angle twisted bilayer graphene (2023), arXiv:2303.14903 [cond-mat.mes-hall] .
  * Goodwin  _et al._ [2019] Z. A. H. Goodwin, F. Corsetti, A. A. Mostofi, and J. Lischner, Attractive electron-electron interactions from internal screening in magic-angle twisted bilayer graphene, Phys. Rev. B 100, 235424 (2019).
  * Ramires and Lado [2019] A. Ramires and J. L. Lado, Impurity-induced triple point fermions in twisted bilayer graphene, Phys. Rev. B 99, 245118 (2019).

##  Appendix A The symmetry operations under magnetic fields

We look for unitary operators realizing the $C_{3z}$ and $C_{2z}$ symmetries,
acting on the creation operators as

$\displaystyle
gc^{\dagger}_{\boldsymbol{i}}g^{-1}=\exp(i\chi_{g}(g(\boldsymbol{r_{i}})))c^{\dagger}_{g(\boldsymbol{i})}.$
(42)

Here we use indistinctly $g$ for the unitary operators and for the linear
transformations acting on points of the lattice. These can always be
distinguished by the context. As in the main text,
$c^{\dagger}_{\boldsymbol{i}}(c^{\dagger}_{g(\boldsymbol{i})})$ is the
creation operator at position $\boldsymbol{r_{i}}$ $(g(\boldsymbol{r_{i}}))$.
The action on the Hamiltonian is

$\displaystyle
gH_{0}g^{-1}=\sum_{\boldsymbol{i},\boldsymbol{j}}t(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})\exp(i\theta_{\boldsymbol{i},\boldsymbol{j}})\exp(i(\chi_{g}(g(\boldsymbol{r_{i}})-\chi_{g}(g(\boldsymbol{r_{j}})))c_{g(\boldsymbol{i})}^{\dagger}c_{g(\boldsymbol{j})}.$
(43)

We are dealing with symmetries at zero flux, so
$t(g(\boldsymbol{r_{i}})-g(\boldsymbol{r_{j}}))=t(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})$.
Then to realize the symmetry (this is, for $gH_{0}g^{-1}=H_{0}$)
$\chi_{g}(\boldsymbol{r})$ must obey

$\displaystyle\theta_{g^{-1}(\boldsymbol{i}),g^{-1}(\boldsymbol{j})}+\chi_{g}(\boldsymbol{r_{i}})-\chi_{g}(\boldsymbol{j})=\theta_{\boldsymbol{i},\boldsymbol{j}}$
$\displaystyle\frac{2\pi}{\Phi_{0}}\int_{g^{-1}(\boldsymbol{r_{i}})}^{g^{-1}(\boldsymbol{r_{j}})}\boldsymbol{A}(\boldsymbol{r^{\prime}})\cdot
d\boldsymbol{r^{\prime}}-\frac{2\pi}{\Phi_{0}}\int_{\boldsymbol{r_{i}}}^{\boldsymbol{r_{j}}}\boldsymbol{A}(\boldsymbol{r^{\prime}})\cdot
d\boldsymbol{r^{\prime}}=\int_{\boldsymbol{r_{i}}}^{\boldsymbol{r_{j}}}\boldsymbol{\nabla}\chi_{g}(\boldsymbol{r^{\prime}})\cdot
d\boldsymbol{r^{\prime}}$
$\displaystyle\frac{2\pi}{\Phi_{0}}\Bigg{(}g\big{(}\boldsymbol{A}(g^{-1}(\boldsymbol{r}))\big{)}-\boldsymbol{A}(\boldsymbol{r})\Bigg{)}=\boldsymbol{\nabla}\chi_{g}(\boldsymbol{r}).$
(44)

In the periodic Landau gauge
$\boldsymbol{A}(\boldsymbol{r})=\frac{\Phi}{2\pi}\Bigg{(}\xi_{1}\boldsymbol{G}_{2}-2\pi\boldsymbol{\nabla}\big{(}\xi_{2}\left\lfloor{\xi_{1}+\epsilon}\right\rfloor\big{)}\Bigg{)}$
where $\xi_{1}$ and $\xi_{2}$ are defined by
$\boldsymbol{r}=\xi_{1}\boldsymbol{L}_{1}+\xi_{2}\boldsymbol{L}_{2}$.

We have for $C_{3z}$

$\displaystyle
C_{3z}\big{(}\boldsymbol{A}(C_{3z}^{-1}(\boldsymbol{r}))\big{)}=\frac{\Phi}{2\pi}\Bigg{(}-\xi_{2}(\boldsymbol{G}_{2}-\boldsymbol{G}_{1})+2\pi\boldsymbol{\nabla}\bigg{(}(\xi_{2}+\xi_{1})\left\lfloor{\xi_{2}+\epsilon}\right\rfloor\bigg{)}\Bigg{)},$
(45)

and hence

$\displaystyle\chi_{C_{3z}}(\boldsymbol{r})=\frac{2\pi
p}{q}\Bigg{(}(\xi_{1}+\xi_{2})\lfloor\xi_{2}+\epsilon\rfloor+\xi_{2}\lfloor\xi_{1}+\epsilon\rfloor-\xi_{1}\xi_{2}-\frac{\xi_{2}^{2}}{2}\Bigg{)}.$
(46)

Similarly for $C_{2z}$ we get

$\displaystyle
C_{2z}\Big{(}\boldsymbol{A}\big{(}C_{2z}^{-1}(\boldsymbol{r})\big{)}\Big{)}$
$\displaystyle=\frac{\Phi}{2\pi}\Bigg{(}\xi_{1}\boldsymbol{G}_{2}+2\pi\boldsymbol{\nabla}\Big{(}\xi_{2}\left\lfloor{-\xi_{1}+\epsilon}\right\rfloor\Big{)}\Bigg{)},$
$\displaystyle\chi_{C_{2z}}(\boldsymbol{r})$ $\displaystyle=\frac{2\pi
p}{q}\Bigg{(}\xi_{2}\lfloor\xi_{1}+\epsilon\rfloor+\xi_{2}\lfloor-\xi_{1}+\epsilon\rfloor\Bigg{)}.$
(47)

Above we have used the facts that for orthogonal transformations $g$ and
scalar functions $f(\boldsymbol{r})$ and
$h(\boldsymbol{r})=f(g^{-1}(\boldsymbol{r}))$, we have
$\boldsymbol{\nabla}h|_{\boldsymbol{r}}=g\big{(}\boldsymbol{\nabla}f|_{g^{-1}(\boldsymbol{r})}\big{)}$,
and that for a function of $\xi_{1}$ and $\xi_{2}$ we have
$2\pi\boldsymbol{\nabla}f=\frac{\partial
f}{\partial\xi_{1}}\boldsymbol{G}_{1}+\frac{\partial
f}{\partial\xi_{2}}\boldsymbol{G}_{2}$. The functions $\chi_{C_{3z}}$ and
$\chi_{C_{2z}}$ have the following periodicity properties

$\displaystyle\chi_{C_{3z}}(\boldsymbol{r}+q\boldsymbol{L}_{2})=\chi_{C_{3z}}(\boldsymbol{r})+\pi
pq\ \text{mod}\ 2\pi,$
$\displaystyle\quad\chi_{C_{3z}}(\boldsymbol{r}+\boldsymbol{L}_{1})=\chi_{C_{3z}}(\boldsymbol{r})+\frac{2\pi
p}{q}\lfloor\xi_{2}+\epsilon\rfloor$
$\displaystyle\chi_{C_{2z}}(\boldsymbol{r}+q\boldsymbol{L}_{2})=\chi_{C_{2z}}(\boldsymbol{r})\
\text{mod}\ 2\pi,$
$\displaystyle\quad\chi_{C_{2z}}(\boldsymbol{r}+\boldsymbol{L}_{1})=\chi_{C_{2z}}(\boldsymbol{r}).$
(48)

We are interested in $p=q=1$, so we can write

$\displaystyle e^{i\chi_{C_{3z}}(\boldsymbol{r})}=$ $\displaystyle
e^{-i\boldsymbol{G}_{2}\cdot\boldsymbol{r}/2}e^{i\overline{\chi}_{C_{3z}}(\boldsymbol{r})}$
$\displaystyle e^{i\chi_{C_{2z}}(\boldsymbol{r})}=$ $\displaystyle
e^{i\overline{\chi}_{C_{2z}}(\boldsymbol{r})},$ (49)

where barred phases are periodic in the Moiré unit cell. As we will see now,
the phases $e^{i\chi_{C_{3z}}(\boldsymbol{r})}$ and
$e^{i\chi_{C_{2z}}(\boldsymbol{r})}$ modify the transformations of the Bloch
waves, redefining the high symmetry points in flux.

The Bloch waves are written

$\displaystyle
c^{\dagger}_{\boldsymbol{k},\boldsymbol{i}}=\frac{1}{\sqrt{N_{M}}}\sum_{\boldsymbol{l}}e^{i\boldsymbol{k}\cdot(\boldsymbol{R_{l}}+\boldsymbol{\delta_{i}})}c^{\dagger}_{\boldsymbol{l},\boldsymbol{i}},$
(50)

with $\boldsymbol{k}$ belonging to the Moiré Brillouin zone, and here
$c^{\dagger}_{\boldsymbol{l},\boldsymbol{i}}$ creates an electron at position
$\boldsymbol{R_{l}}+\boldsymbol{\delta_{i}}$ where $\boldsymbol{R_{l}}$ is a
lattice vector and $\boldsymbol{\delta_{i}}$ belongs to the Wigner-Seitz cell.

Under $C_{3z}$, $c^{\dagger}_{\boldsymbol{k},\boldsymbol{i}}$ transforms as

$\displaystyle
C_{3z}c^{\dagger}_{\boldsymbol{k},\boldsymbol{i}}(C_{3z})^{-1}=\frac{1}{\sqrt{N_{M}}}\sum_{\boldsymbol{l}}e^{i(C_{3z}(\boldsymbol{k})-\boldsymbol{G}_{2}/2)\cdot
C_{3z}(\boldsymbol{R_{l}}+\boldsymbol{\delta_{i}})}e^{i\overline{\chi}_{C_{3z}}(C_{3z}(\boldsymbol{\delta_{i}}))}c^{\dagger}_{C_{3z}(\boldsymbol{l},\boldsymbol{i})}.$
(51)

Here, $c^{\dagger}_{C_{3z}(\boldsymbol{l},\boldsymbol{i})}$ creates an
electron at position $C_{3z}(\boldsymbol{R_{l}}+\boldsymbol{\delta_{i}})$. We
see that $C_{3z}$ sends momentum $\boldsymbol{k}$ to
$C_{3z}(\boldsymbol{k})-\boldsymbol{G}_{2}/2$. Via the embedding relation
$c^{\dagger}_{\boldsymbol{k}+\boldsymbol{G},\boldsymbol{i}}=e^{i\boldsymbol{G}\cdot\boldsymbol{\delta_{i}}}c^{\dagger}_{\boldsymbol{k},\boldsymbol{i}}$
for $\boldsymbol{G}$ a reciprocal lattice vector, the three-fold rotation in
flux acts in the momenta as follows,

$\displaystyle\boldsymbol{k}\xrightarrow{C_{3z}}C_{3z}(\boldsymbol{k})-\boldsymbol{G}_{2}/2\sim
C_{3z}\Big{(}\boldsymbol{k}-(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2\Big{)}+(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2.$
(52)

Also, given that $\chi_{C_{2z}}(\boldsymbol{r})$ is periodic mod $2\pi$ on the
unit cell, the momentum transforms like in zero flux,

$\displaystyle\boldsymbol{k}\xrightarrow{C_{2z}}C_{2z}(\boldsymbol{k})\sim
C_{2z}\Big{(}\boldsymbol{k}-(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2\Big{)}+(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2.$
(53)

The center of rotations has shifted from $\boldsymbol{\Gamma_{M}}$ to one of
the other parity invariant points,
$(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2$, at one magnetic flux quantum.

Now we look for the operator realizing $C_{2y}$. The procedure is the same,
but in this case $C_{2y}H_{0}(C_{2y})^{-1}$ should be equal to the Hamiltonian
with the sign of the magnetic field reversed. Hence,
$\chi_{C_{2y}}(\boldsymbol{r})$ must obey

$\displaystyle
C_{2y}\big{(}\boldsymbol{A}(C_{2y}^{-1}(\boldsymbol{r}))\big{)}+\boldsymbol{A}(\boldsymbol{r})=\boldsymbol{\nabla}\chi_{C_{2y}}(\boldsymbol{r}).$
(54)

We obtain for $\chi_{C_{2y}}(\boldsymbol{r})$

$\displaystyle\chi_{C_{2y}}(\boldsymbol{r})=\frac{2\pi
p}{q}\Bigg{(}-\xi_{2}\lfloor\xi_{1}+\epsilon\rfloor+\xi_{2}\lfloor\xi_{1}+\xi_{2}+\epsilon\rfloor-\frac{\xi_{2}^{2}}{2}\Bigg{)},$
(55)

which obeys the properties

$\displaystyle\chi_{C_{2y}}(\boldsymbol{r}+q\boldsymbol{L}_{2})=\chi_{C_{2y}}(\boldsymbol{r})-\pi
pq\ \text{mod}\ 2\pi,$
$\displaystyle\quad\chi_{C_{2y}}(\boldsymbol{r}+\boldsymbol{L}_{1})=\chi_{C_{2y}}(\boldsymbol{r}).$
(56)

The same analysis as before applies, and under $C_{2y}$ the momentum transform
as

$\displaystyle\boldsymbol{k}\xrightarrow{C_{2y}}C_{2y}(\boldsymbol{k})-\boldsymbol{G}_{2}/2\sim
C_{2y}\Big{(}\boldsymbol{k}-(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2\Big{)}+(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2.$
(57)

For the time reversal operator $\mathcal{T}$, the magnetic field should be
reversed also, and it is trivial to see that the action is the same as for the
zero flux case. $\mathcal{T}$ is an antiunitary operator satisfying
$\mathcal{T}c^{\dagger}_{\boldsymbol{i}}\mathcal{T}^{-1}=c_{\boldsymbol{i}}^{\dagger}$,
and transforming the momentum as

$\displaystyle\boldsymbol{k}\xrightarrow{\mathcal{T}}-\boldsymbol{k}\sim-\Big{(}\boldsymbol{k}-(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2\Big{)}+(\boldsymbol{G}_{1}+\boldsymbol{G}_{2})/2.$
(58)

In conclusion, the action of symmetry operators under one magnetic flux
quantum effectively shift the Brillouin zone by
$(\boldsymbol{G}_{1}+\boldsymbol{G_{2}})/2$, redefining the high symmetry
points.

##  Appendix B Valley charge and $C_{2z}P$ operator on the lattice

We wish to find an operator $\tau_{z}$ implementing the valley charge on the
lattice, such that $\langle\tau_{z}\rangle=+1$ on states nearby the $K$ point
of graphene and $-1$ near the $K^{\prime}$ point. We adopt a slight
generalization of the valley operator of Ref. [76]

$\displaystyle\tau_{z}=\frac{i}{3\sqrt{3}}\sum_{l}\bigg{(}$
$\displaystyle\sum_{\bigtriangledown}e^{-i\theta_{\bigtriangledown(1),\bigtriangledown(2)}}c^{\dagger}_{\bigtriangledown(1)}c_{\bigtriangledown(2)}+e^{-i\theta_{\bigtriangledown(2),\bigtriangledown(3)}}c^{\dagger}_{\bigtriangledown(2)}c_{\bigtriangledown(3)}+e^{-i\theta_{\bigtriangledown(3),\bigtriangledown(1)}}c^{\dagger}_{\bigtriangledown(3)}c_{\bigtriangledown(1)}$
$\displaystyle-\sum_{\triangle}e^{-i\theta_{\triangle(1),\triangle(2)}}c^{\dagger}_{\triangle(1)}c_{\triangle(2)}+e^{-i\theta_{\triangle(2),\triangle(3)}}c^{\dagger}_{\triangle(2)}c_{\triangle(3)}+e^{-i\theta_{\triangle(3),\triangle(1)}}c^{\dagger}_{\triangle(3)}c_{\triangle(1)}\bigg{)}+\text{h.c.}.$
(59)

The sums are over triangles upside down of sublattice $A$ atoms, and triangles
of sublattice $B$, and $l$ denotes the sum over the two layers. We draw an
example of each kind of triangle in Figure C.3. The phases are the Peierls’
phases defined in the main text. It can be shown that valley $K$ states have
$\langle\tau_{z}\rangle=+1+O(a/L_{M})$ and valley $K^{\prime}$ states have
$\langle\tau_{z}\rangle=-1+O(a/L_{M})$.

For the particle-hole operator $C_{2z}P$ (Eq.11) we proceed as follows. Once
we have valley polarized states, we obtain the envelope functions
$\Psi_{\eta\sigma l}$ by multiplying $\Psi_{\eta}$ with the valley factor
$e^{-i\eta n_{\theta}\boldsymbol{G_{2}}\cdot\boldsymbol{r}}$. Afterwards, we
perform a smooth interpolation of the data $\Psi_{\eta\sigma
l}(\boldsymbol{r_{i}})$ being $\boldsymbol{r_{i}}$ the positions of the atoms
at sublattice $\sigma$ and layer $l$. Finally, the smooth functions are
sampled at the points of the opposite sublattice and layer, and the factor
$\eta s_{l}e^{-i\eta
n_{\theta}\boldsymbol{G_{2}}\cdot\boldsymbol{r}}e^{-i\eta\boldsymbol{G_{2}}\cdot\boldsymbol{r}}$
is also added. As a note, When $\Phi=\Phi_{0}$ the envelopes have a
discontinuity at $\xi_{1}=$ an integer in the periodic Landau gauge, and
special care is needed when performing the interpolation.

The projected operator in the flat bands
$[\overline{C_{2z}P}(\boldsymbol{k})]_{\rho\rho^{\prime}}=\langle\boldsymbol{k}\rho|C_{2z}P|\boldsymbol{k}\rho^{\prime}\rangle=\langle\boldsymbol{k}\rho|C_{2z}P(\boldsymbol{k}\rho^{\prime})\rangle$
is then constructed. We have checked that the particular choice of basis for
the interpolation is irrelevant, and the matrix elements of
$\overline{C_{2z}P}(\boldsymbol{k})$ in a different basis computed via
interpolation or unitary conjugation of the original matrix are essentially
identical.

In order to check the validity of our definitions, we have computed several
benchmarks. The mean of the absolute value of the matrix elements of
$\overline{C_{2z}P}(\boldsymbol{k})-\overline{C_{2z}P}(\boldsymbol{k})^{\dagger}$
are always less than $10^{-5}$, so we capture well the hermiticity of
$C_{2z}P$. Also, the mean of the absolute value of the valley off-diagonal
matrix elements of $\overline{C_{2z}P}(\boldsymbol{k})$ are always less than
$10^{-5}$ for $\Phi_{0}$ flux and $10^{-7}$ at zero field. We conclude that
our implementations of the valley charge and the p-h operator are trustworthy.

Figure B.1: $\boldsymbol{C_{2z}P}$ operator at zero flux. Left: for the state
$|\boldsymbol{M_{M}}KA\rangle$ at zero field we plot the smooth envelope at
sublattice $A$ bottom layer, $\Psi_{KAb}(\boldsymbol{r})$. Right: for the
$C_{2z}P$ transformed state, we plot
$e^{in_{\theta}\boldsymbol{G}_{2}\cdot\boldsymbol{r}}e^{i\boldsymbol{G}_{2}\cdot\boldsymbol{r}}C_{2z}P(\Psi_{K})_{Bt}(\boldsymbol{r})$.

Figure B.2: $\boldsymbol{C_{2z}P}$ operator at one flux quantum. In this case
we plot the envelope of $|\boldsymbol{M_{M}}K^{\prime}+1\rangle$ at $26.5$ T
in the right, and
$e^{-in_{\theta}\boldsymbol{G}_{2}\cdot\boldsymbol{r}}e^{-i\boldsymbol{G}_{2}\cdot\boldsymbol{r}}C_{2z}P(\Psi_{K}^{\prime})_{Bt}(\boldsymbol{r})$
in the left plot. Notice the discontinuities of the wave function on the unit
cell edge due to the periodic Landau gauge, and the extra minus sign $\eta=-1$
in the transformed wave fucntion, see Eq. 11.

##  Appendix C $\boldsymbol{U(1)_{v}}$ AND $\boldsymbol{SU(2)_{K}\times
SU(2)_{K^{\prime}}}$ symmetries in the lattice model

In the BM model, a general wave function can be written similarly to Eq. 10,
in a $4$ component notation for sublattice and layer ($t$=top, $b$=bottom),

$\displaystyle
u(\boldsymbol{r})=e^{i\boldsymbol{K_{u}}\cdot\boldsymbol{r}}(u_{At}(\boldsymbol{r}),u_{Bt}(\boldsymbol{r}),u_{Ab}(\boldsymbol{r}),u_{Bb}(\boldsymbol{r}))^{T},$
(60)

with $\boldsymbol{K_{u}}=\pm\boldsymbol{K}$, the graphene Dirac points, and
the envelope functions are smooth on the graphene scale. The matrix elements
of the Coulomb interaction (actually, the four fermion part of the normal
ordered operator) read

$\displaystyle\langle
ab|V|cd\rangle=\frac{1}{2}\sum_{\\{\sigma\\},\\{l\\}}\int
d\boldsymbol{r_{i}}\int
d\boldsymbol{r_{j}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})e^{i(\boldsymbol{K_{c}}-\boldsymbol{K_{a}})\cdot\boldsymbol{r_{i}}}e^{i(\boldsymbol{K_{d}}-\boldsymbol{K_{b}})\cdot\boldsymbol{r_{j}}}a^{*}_{\sigma_{a}l_{a}}(\boldsymbol{r_{i}})b^{*}_{\sigma_{b}l_{b}}(\boldsymbol{r_{j}})c_{\sigma_{c}l_{c}}(\boldsymbol{r_{i}})d_{\sigma_{d}l_{d}}(\boldsymbol{r_{j}}).$
(61)

This integral computes the Fourier transform of the potential at momenta
$\sim\boldsymbol{K_{a}}-\boldsymbol{K_{c}}$, or equivalently
$\sim\boldsymbol{K_{d}}-\boldsymbol{K_{b}}$. Hence, unless
$\boldsymbol{K_{a}}=\boldsymbol{K_{c}}$ and
$\boldsymbol{K_{b}}=\boldsymbol{K_{d}}$ the form factor is negligible given
$||\boldsymbol{K}||\sim a^{-1}\gg\xi^{-1}$. This is the origin of $U(1)_{v}$
and $SU(2)_{K}\times SU(2)_{K^{\prime}}$ in TBG.

Figure C.3: Left: Triangular loops that calculate the valley charge. Right:
Decomposition of the potential. $V_{sm}(r)$ is a smoothed version of $V(r)$

In the tight-binding, we write states $|u\rangle$, with wavefunctions

$\displaystyle\langle\boldsymbol{r_{i}}|u\rangle=e^{i\boldsymbol{K_{u}}\cdot\boldsymbol{r_{i}}}u(\boldsymbol{r_{i}}).$
(62)

$u(\boldsymbol{r_{i}})$ is further decomposed depending on the sublattice and
layer of the point $\boldsymbol{r_{i}}$, $u(\boldsymbol{r_{i}})=\sum_{\sigma
l}u_{\sigma l}(\boldsymbol{r_{i}})\delta_{\boldsymbol{r_{i}}\in\sigma l}$. The
functions $u_{\sigma l}(\boldsymbol{r_{i}})$ vary slowly on the graphene
scale, such that $u_{\sigma l}(\boldsymbol{r_{i}}+\boldsymbol{a}_{1,2})\approx
u_{\sigma l}(\boldsymbol{r_{i}})$. Consider the matrix element of the
interaction (actually, the four fermion part of the normal ordered operator)

$\displaystyle\langle
ab|V|cd\rangle=\frac{1}{2}\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})e^{-i(\boldsymbol{K_{a}}-\boldsymbol{K_{c}})\cdot\boldsymbol{r_{i}}}e^{-i(\boldsymbol{K_{b}}-\boldsymbol{K_{d}})\cdot\boldsymbol{r_{j}}}a^{*}(\boldsymbol{r_{i}})b^{*}(\boldsymbol{r_{j}})c(\boldsymbol{r_{i}})d(\boldsymbol{r_{j}})$
(63)

We will see that if the valley charge is not conserved, i.e.
$\boldsymbol{K_{a}}+\boldsymbol{K_{b}}-\boldsymbol{K_{c}}-\boldsymbol{K_{d}}\neq\boldsymbol{0}$
the matrix element vanishes to a first approximation. Let us take for example
$\boldsymbol{K_{a}}=\boldsymbol{K_{c}}=\boldsymbol{K_{d}}=-\boldsymbol{K_{b}}=\boldsymbol{K}$,
then

$\displaystyle\langle
ab|V|cd\rangle=\frac{1}{2}\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})e^{-2i\boldsymbol{K}\cdot\boldsymbol{r_{i}}}a^{*}(\boldsymbol{r_{i}})b^{*}(\boldsymbol{r_{j}})c(\boldsymbol{r_{i}})d(\boldsymbol{r_{j}}).$
(64)

Shifting the integration variables
$\boldsymbol{r_{i,j}}\to\boldsymbol{r_{i,j}}+\boldsymbol{a}_{1}$ and using the
property of the envelope functions
$u(\boldsymbol{r_{i}}+\boldsymbol{a}_{1,2})\approx u(\boldsymbol{r_{i}})$ we
get

$\displaystyle\langle
ab|V|cd\rangle=e^{-2i\boldsymbol{K}\cdot\boldsymbol{a}_{1}}\langle
ab|V|cd\rangle=e^{-2\pi i/3}\langle ab|V|cd\rangle=0.$ (65)

The argument for other charge non conserving combinations of
$\boldsymbol{K_{a,b,c,d}}$ is the same and we conclude that $U(1)_{v}$ is a
symmetry in the atomistic model also.

If
$\boldsymbol{K_{a}}=\boldsymbol{K_{d}}=-\boldsymbol{K_{b}}=-\boldsymbol{K_{c}}$,
say $\boldsymbol{K_{a}}=\boldsymbol{K}$, we have

$\displaystyle\langle
ab|V|cd\rangle=\frac{1}{2}\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})e^{-2i\boldsymbol{K}\cdot(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})}a^{*}(\boldsymbol{r_{i}})b^{*}(\boldsymbol{r_{j}})c(\boldsymbol{r_{i}})d(\boldsymbol{r_{j}}),$
(66)

and the shifting argument does not work because the phases for
$\boldsymbol{r_{i}}$ and $\boldsymbol{r_{j}}$ cancel. In general these matrix
elements are non zero and break the $SU(2)_{K}\times SU(2)_{K^{\prime}}$
symmetry of independent spin rotations in each valley.

Furthermore, the exchange energy of states with different valley charge, say
$|a\rangle=|d\rangle$ and $|b\rangle=|c\rangle$, is always positive. For that,
decompose the potential into a smooth part $V_{sm}(\boldsymbol{r})$ and a
short-range part $V_{sr}(\boldsymbol{r})$. $V_{sm}(\boldsymbol{r})$ is slowly
varying on the graphene scale, is equal to $V(\boldsymbol{r})$ at long
distances (say $||\boldsymbol{r}||>2a$) and $V_{sm}(\boldsymbol{0})$ equals a
positive constant. On the other hand,
$V_{sr}(\boldsymbol{r})=V(\boldsymbol{r})-V_{sm}(\boldsymbol{r})$ is equal to
$-V_{sm}(\boldsymbol{0})$ at $\boldsymbol{r}=0$ (remember the on-site
interaction is treated separately by the Hubbard term, so
$V(\boldsymbol{0})=0$), exhibits the steep Coulomb repulsion at short
distances and vanishes at long distances. In Fig. C.3 we sketch this
decompositon.

Given that $V_{sm}(\boldsymbol{r}+\boldsymbol{a}_{1,2})\approx
V_{sm}(\boldsymbol{r})$, one can apply the shifting only to the fisrt variable
$\boldsymbol{r_{i}}$, and get

$\displaystyle\langle
ab|V_{sm}|ba\rangle=e^{-2i\boldsymbol{K}\cdot\boldsymbol{a}_{1}}\langle
ab|V_{sm}|ba\rangle=e^{-2\pi i/3}\langle ab|V_{sm}|ba\rangle=0.$ (67)

For $V_{sr}(\boldsymbol{r})$ let us consider only the dominant contributions
of intralayer terms with $||\boldsymbol{r_{i}}-\boldsymbol{r_{j}}||=0,a_{0}$
and $a$.

$\displaystyle\langle ab|V_{sr}|ba\rangle=\frac{1}{2}$
$\displaystyle\sum_{\boldsymbol{r_{j}}}V_{0}|a_{\sigma_{j}l_{j}}(\boldsymbol{r_{j}})|^{2}|b_{\sigma_{j}l_{j}}(\boldsymbol{r_{j}})|^{2}+3V_{a}(e^{2\pi
i/3}+e^{-2\pi
i/3})|a_{\sigma_{j}l_{j}}(\boldsymbol{r_{j}})|^{2}|b_{\sigma_{j}l_{j}}(\boldsymbol{r_{j}})|^{2}$
$\displaystyle\hskip 14.22636pt+V_{a_{0}}(1+e^{2\pi i/3}+e^{-2\pi
i/3})a^{*}_{\bar{\sigma}_{j}l_{j}}(\boldsymbol{r_{j}})a_{\sigma_{j}l_{j}}(\boldsymbol{r_{j}})b^{*}_{\sigma_{j}l_{j}}(\boldsymbol{r_{j}})b_{\bar{\sigma}_{j}l_{j}}(\boldsymbol{r_{j}})$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\boldsymbol{r_{j}}}(V_{0}-3V_{a})|a(\boldsymbol{r_{j}})|^{2}|b(\boldsymbol{r_{j}})|^{2}<0,$
(68)

$\sigma_{j}$ and $l_{j}$ denoting the sublattice and layer of the point
$\boldsymbol{r_{j}}$ and $\bar{\sigma}_{j}$ the opposite sublattice to
$\sigma_{j}$. $V_{r}$ encodes the value of $V_{sr}(\boldsymbol{r})$ when
$||\boldsymbol{r}||=r$. Notice that $V_{0}<0$ and $V_{a}>0$. In turn, the
exchange energy $-\langle ab|V|ba\rangle$ is always positive

##  Appendix D The Hartree-Fock method and flat band projection

Consider the normal ordered interaction of Eqs. 5 and 6,

$\displaystyle
V+H_{U}=\frac{1}{2}\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}}s_{i}s_{j}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}}):c^{\dagger}_{\boldsymbol{i},s_{i}}c_{\boldsymbol{i},s_{i}}c^{\dagger}_{\boldsymbol{j},s_{j}}c_{\boldsymbol{j},s_{j}}:+U\sum_{\boldsymbol{i}}:c_{\boldsymbol{i}\uparrow}^{\dagger}c_{\boldsymbol{r_{i}}\uparrow}c_{\boldsymbol{i}\downarrow}^{\dagger}c_{\boldsymbol{i}\downarrow}:$
(69)

The choice of the normal ordering with respect to the ground state of graphene
at charge neutrality is necessary to avoid double counting the interaction[30,
21]. This is, we assume that the hopping integrals $t(\boldsymbol{r})$ are
already renormalized by the interactions with the deep Fermi sea of graphene.
After expanding the normal ordered product[59] and performing the Hartree-Fock
decoupling, the Hamiltonian reads

$\displaystyle V_{\text{HF}}+H_{U\text{HF}}=$
$\displaystyle\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}},s_{i},s_{j}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})c_{\boldsymbol{i}s_{i}}^{\dagger}c_{\boldsymbol{i}s_{i}}\Big{(}\langle
c_{\boldsymbol{j}s_{j}}^{\dagger}c_{\boldsymbol{j}s_{j}}\rangle-\langle
c^{\dagger}_{\boldsymbol{j}s_{j}}c_{\boldsymbol{j}s_{j}}\rangle_{0}\Big{)}-\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}},s}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})c_{\boldsymbol{i}s}^{\dagger}c_{\boldsymbol{j}s}\Big{(}\langle
c_{\boldsymbol{i}s}^{\dagger}c_{\boldsymbol{j}s}\rangle-\langle
c^{\dagger}_{\boldsymbol{i}s}c_{\boldsymbol{j}s}\rangle_{0}\Big{)}^{*}$
$\displaystyle+U\sum_{\boldsymbol{r_{i}}}c_{\boldsymbol{i}\uparrow}^{\dagger}c_{\boldsymbol{i}\uparrow}\Big{(}\langle
c_{\boldsymbol{i}\downarrow}^{\dagger}c_{\boldsymbol{i}\downarrow}\rangle-\langle
c_{\boldsymbol{i}\downarrow}^{\dagger}c_{\boldsymbol{i}\downarrow}\rangle_{0}\Big{)}+U\sum_{\boldsymbol{r_{i}}}c_{\boldsymbol{i}\downarrow}^{\dagger}c_{\boldsymbol{i}\downarrow}\Big{(}\langle
c_{\boldsymbol{i}\uparrow}^{\dagger}c_{\boldsymbol{i}\uparrow}\rangle-\langle
c_{\boldsymbol{i}\uparrow}^{\dagger}c_{\boldsymbol{i}\uparrow}\rangle_{0}\Big{)}+\text{constant},$
(70)

with $\langle...\rangle_{0}$ denoting the expectation value in the ground
state of graphene at charge neutrality, and $\langle...\rangle$ the
expectation value in the particular state of our Hartree-Fock decoupling. In
our implementation we restrict the wave function to be a direct product of
spin up and down electrons, such that $\langle
c^{\dagger}_{\boldsymbol{i}\uparrow}c_{\boldsymbol{j}\downarrow}\rangle=0$ for
all $\boldsymbol{r_{i}},\boldsymbol{r_{j}}$.

In the projected limit we assume that the remote bands are filled and the
relevant physics takes place in the flat bands. In this spirit we compute mean
field interaction restricted to the subspace of the flat bands,

$\displaystyle[V_{\text{HF,p}}(\boldsymbol{k},\boldsymbol{k^{\prime}})+H_{U\text{HF,p}}(\boldsymbol{k},\boldsymbol{k^{\prime}})]_{\rho\rho^{\prime}}=\Big{(}\langle\text{FS}|\otimes\langle\boldsymbol{k}\rho|\Big{)}\Big{(}V_{\text{HF}}+H_{U\text{HF}}\Big{)}\Big{(}|\text{FS}\rangle\otimes|\boldsymbol{k^{\prime}}\rho^{\prime}\rangle\Big{)},$
(71)

with $|\text{FS}\rangle\otimes|\boldsymbol{k}\rho\rangle$ denoting the direct
product of the state with the filled remote bands and the state with momentum
$\boldsymbol{k}$ and multi-index $\rho$. We further assume translational
symmetry that makes the mean field Hamiltonian block-diagonal in momentum
space,
$V_{\text{HF,p}}(\boldsymbol{k},\boldsymbol{k^{\prime}})+H_{U\text{HF,p}}(\boldsymbol{k},\boldsymbol{k^{\prime}})=\Big{(}V_{\text{HF,p}}(\boldsymbol{k})+H_{U\text{HF,p}}(\boldsymbol{k})\Big{)}\delta_{\boldsymbol{k},\boldsymbol{k^{\prime}}}$.

The self-consistent method starts by proposing an ansatz for the ground state
at any given filling, computing the mean field Hamiltonian and performing the
flat band projection. Next, we solve the projected mean filed Hamiltonian

$\displaystyle
H_{\text{HF,p}}(\boldsymbol{k},\boldsymbol{k^{\prime}})=\Big{(}H_{0\text{,p}}(\boldsymbol{k})+V_{\text{HF,p}}(\boldsymbol{k})+H_{U\text{HF,p}}(\boldsymbol{k})\Big{)}\delta_{\boldsymbol{k},\boldsymbol{k^{\prime}}},$
(72)

with
$H_{0\text{,p}}(\boldsymbol{k})\delta_{\boldsymbol{k},\boldsymbol{k^{\prime}}}$
the projected kinetic energy operator. The ground state of this Hamiltonian is
then a new ansatz for the self-consistent ground state and the process is
repeated until convergence is reached.

The energy of the self-consistent state is

$\displaystyle\langle H\rangle=$ $\displaystyle\langle V\rangle+\langle
H_{U}\rangle+\langle H_{0}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}},s_{i},s_{j}}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})\Big{(}\langle
c_{\boldsymbol{i}s_{i}}^{\dagger}c_{\boldsymbol{i}s_{i}}\rangle-\langle
c^{\dagger}_{\boldsymbol{i}s_{i}}c_{\boldsymbol{i}s_{i}}\rangle_{0}\Big{)}\Big{(}\langle
c_{\boldsymbol{j}s_{j}}^{\dagger}c_{\boldsymbol{j}s_{j}}\rangle-\langle
c^{\dagger}_{\boldsymbol{j}s_{j}}c_{\boldsymbol{j}s_{j}}\rangle_{0}\Big{)}-\frac{1}{2}\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}},s}V(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})\Big{\lvert}\Big{\lvert}\langle
c_{\boldsymbol{i}s}^{\dagger}c_{\boldsymbol{j}s}\rangle-\langle
c^{\dagger}_{\boldsymbol{i}s}c_{\boldsymbol{j}s}\rangle_{0}\Big{\rvert}\Big{\rvert}^{2}$
$\displaystyle+U\sum_{\boldsymbol{r_{i}}}\Big{(}\langle
c_{\boldsymbol{i}\uparrow}^{\dagger}c_{\boldsymbol{i}\uparrow}\rangle-\langle
c^{\dagger}_{\boldsymbol{i}\uparrow}c_{\boldsymbol{i}\uparrow}\rangle_{0}\Big{)}\Big{(}\langle
c_{\boldsymbol{i}\downarrow}^{\dagger}c_{\boldsymbol{i}\downarrow}\rangle-\langle
c^{\dagger}_{\boldsymbol{i}\downarrow}c_{\boldsymbol{i}\downarrow}\rangle_{0}\Big{)}+\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}}s}t(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})e^{i\theta_{\boldsymbol{i},\boldsymbol{j}}}\langle
c_{\boldsymbol{i}s}^{\dagger}c_{\boldsymbol{j}s}\rangle.$ (73)

The Coulomb interaction is split into the Hartree or direct and Fock or
exchange terms, with the plus and minus signs in front respectively. In our
algorithm, we always work with the Fock matrix $\langle
c_{\boldsymbol{i}s}^{\dagger}c_{\boldsymbol{i}s}\rangle-\langle
c^{\dagger}_{\boldsymbol{i}s}c_{\boldsymbol{i}s}\rangle_{0}$, so the values
reported for the kinetic energy have a constant offset of
$\sum_{\boldsymbol{r_{i}},\boldsymbol{r_{j}}s}t(\boldsymbol{r_{i}}-\boldsymbol{r_{j}})e^{i\theta_{\boldsymbol{i},\boldsymbol{j}}}\langle
c_{\boldsymbol{i}s}^{\dagger}c_{\boldsymbol{j}s}\rangle_{0}$.

##  Appendix E Additional Hartree-Fock results

a)
$\Phi=0$

$\Phi=\Phi_{0}$

b) c) d) e) f)

Figure E.4: Hartree-Fock results for $\boldsymbol{\epsilon=30,\ U=4}$ eV. a)
The KIVC band structures, for $\Phi=0$, and the ground state bands, for
$\Phi_{0}$. Spin up bands are in blue and spin down bands in red. From left to
right the doping levels are $\nu=-2,0,+2$. The critical value for the SP-KIVC
transition is found to be $U_{c}=1.1$ eV. In b) we show $\Delta_{\text{KIVC}}$
for $\nu=-2$ at zero flux. There is a hole pocket around $\Gamma_{M}$ and two
electron pockets at its sides, where the order parameter is ill-defined. In
Fig. 8b) the extra electrons come from the opposite spin, so the electron
pockets are not visible. In this phase $C_{3z}$ is broken but $C_{2z}$ and
$C_{2x}$ are preserved. $C_{3z}$ breaking is not evident in the plot of a)
because the $C_{3z}$ (or $C_{6z}$) related lines are related by $C_{2x}$ (or
$C_{2y}=C_{2z}C_{2x}$) also. In c) we plot $\Delta_{xyz,0}$ of the hole doped
state under one flux quantum, and in d) and e) the two dominant order
parameters for the state at $\nu=+2$ and $\Phi_{0}$ flux. Unlike for
$\epsilon=10$, $U=4$ eV, both parameters contribute equitably. The phase
transition to the the IVC state takes place at some value of $\epsilon$
between $30$ and $50$. f) The path of the band structure plots.

a) b)
$\Phi=0$ $\Phi=\Phi_{0}$

Figure E.5: Correlated insulators at $\boldsymbol{\nu=+3}$. a) Chern
insulators with Chern number $1$ that have been reported in samples with
explicit $C_{2z}$ breaking due to the substrate[25, 24]. We obtain incipient
valley polarized anomalous Hall states without the need for $C_{2z}$ breaking
at zero field, both for $\epsilon=10$, $U=4$ eV, plotted to the left, and
$\epsilon=50$, $U=4$ eV to the right. b) Gapped insulators at $B=26.5$ T. In
both cases we get valley polarized states with Chern number $+1$.
Interestingly, the band structures are very dissimilar for $\epsilon=10$,
$U=4$ eV plotted to the left, and $\epsilon=50$, $U=4$ eV to the right.

$\epsilon=10$, $U=4$ eV

$\epsilon=50$, $U=0.5$ eV

Figure E.6: The self-consistent valley polarized states. Spin up bands are
shown in blue and spin down bands in red. From left to right the doping levels
are $\nu=-2,0,+2$. The order parameter is the valley polarization
$\langle\sigma_{0}\tau_{z}\rangle$ predominantly. Only for $\epsilon=10$,
$U=4$ eV the inter-Chern parameter appears around $\Gamma_{M}$.

##  Appendix F Tables with the energies of the candidate and self-consistent
states

$\nu=-2$, $\xi=10$ nm, $\Phi=\Phi_{0}$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
IVC |  $\prod_{\boldsymbol{k}}\prod_{\lambda}\frac{1}{\sqrt{2}}(d^{\dagger}_{\boldsymbol{k}K\lambda\uparrow}+id^{\dagger}_{\boldsymbol{k}K^{\prime}\lambda\uparrow})|0\rangle$ |  $-48144.42$ |  $328.32$ |  $-38862.09$ |  $-38533.77$ |  $11.52$
VP |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}K+1\uparrow}d^{\dagger}_{\boldsymbol{k}K-1\uparrow}|0\rangle$ |  $-48144.42$ |  $327.78$ |  $-38861.38$ |  $-38533.60$ |  $11.52$
non int. GS |  $-$ |  $-48150.18$ |  $225.32$ |  $-38753.07$ |  $-38527.76$ |  $11.37$

$\nu=-2$, $\xi=20$ nm, $\Phi=\Phi_{0}$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
IVC |  $\prod_{\boldsymbol{k}}\prod_{\lambda}\frac{1}{\sqrt{2}}(d^{\dagger}_{\boldsymbol{k}K\lambda\uparrow}+id^{\dagger}_{\boldsymbol{k}K^{\prime}\lambda\uparrow})|0\rangle$ |  $-48144.42$ |  $919.39$ |  $-40665.49$ |  $-39746.10$ |  $11.52$
VP |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}K+1\uparrow}d^{\dagger}_{\boldsymbol{k}K-1\uparrow}|0\rangle$ |  $-48144.42$ |  $918.85$ |  $-40664.80$ |  $-39745.95$ |  $11.52$
non int. GS |  $-$ |  $-48150.18$ |  $816.56$ |  $-40547.45$ |  $-39730.89$ |  $11.37$

$\nu=0$, $\xi=10$ nm, $\Phi=\Phi_{0}$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
IVC |  $\prod_{\boldsymbol{k}}\prod_{s\lambda}\frac{1}{\sqrt{2}}(d^{\dagger}_{\boldsymbol{k}K\lambda s}+id^{\dagger}_{\boldsymbol{k}K^{\prime}\lambda s})|0\rangle$ |  $-46554.33$ |  $-343.05$ |  $-38855.98$ |  $-39199.03$ |  $11.45$
VP |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}K+1s}d^{\dagger}_{\boldsymbol{k}K-1s}|0\rangle$ |  $-46554.33$ |  $-345.21$ |  $-38854.56$ |  $-39199.77$ |  $11.52$
VSP |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}K+1\uparrow}d^{\dagger}_{\boldsymbol{k}K-1\uparrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}+1\downarrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}-1\downarrow}|0\rangle$ |  $-46554.33$ |  $-342.64$ |  $-38854.56$ |  $-39197.20$ |  $11.33$
SP |  $\prod_{\boldsymbol{k}}\prod_{\eta\lambda}d^{\dagger}_{\boldsymbol{k}\eta\lambda\uparrow}|0\rangle$ |  $-46554.33$ |  $-342.64$ |  $-38849.33$ |  $-39191.97$ |  $10.75$
non int. GS |  $-$ |  $-46565.86$ |  $8.86$ |  $-38637.94$ |  $-38629.08$ |  $11.77$

$\nu=0$, $\xi=20$ nm, $\Phi=\Phi_{0}$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
IVC |  $\prod_{\boldsymbol{k}}\prod_{s\lambda}\frac{1}{\sqrt{2}}(d^{\dagger}_{\boldsymbol{k}K\lambda s}+id^{\dagger}_{\boldsymbol{k}K^{\prime}\lambda s})|0\rangle$ |  $-46554.33$ |  $-346.62$ |  $-40657.36$ |  $-41003.99$ |  $11.45$
VP |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}K+1s}d^{\dagger}_{\boldsymbol{k}K-1s}|0\rangle$ |  $-46554.34$ |  $-348.80$ |  $-40655.99$ |  $-41004.79$ |  $11.52$
VSP |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}K+1\uparrow}d^{\dagger}_{\boldsymbol{k}K-1\uparrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}+1\downarrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}-1\downarrow}|0\rangle$ |  $-46554.33$ |  $-346.21$ |  $-40655.99$ |  $-41002.20$ |  $11.33$
SP |  $\prod_{\boldsymbol{k}}\prod_{\eta\lambda}d^{\dagger}_{\boldsymbol{k}\eta\lambda\uparrow}|0\rangle$ |  $-46554.33$ |  $-346.21$ |  $-40650.70$ |  $-40996.92$ |  $10.75$
non int. GS |  $-$ |  $-46565.86$ |  $4.22$ |  $-40421.29$ |  $-40417.07$ |  $11.77$

$\nu=+2$, $\xi=10$ nm, $\Phi=\Phi_{0}$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
IVC |  $\prod_{\boldsymbol{k}}\prod_{\lambda}\frac{1}{\sqrt{2}}(d_{\boldsymbol{k}K\lambda\downarrow}+id_{\boldsymbol{k}K^{\prime}\lambda\downarrow})|+4\rangle$ |  $-44964.25$ |  $474.60$ |  $-38843.22$ |  $-38368.62$ |  $11.34$
VP |  $\prod_{\boldsymbol{k}}d_{\boldsymbol{k}K+1\downarrow}d_{\boldsymbol{k}K-1\downarrow}|+4\rangle$ |  $-44964.25$ |  $474.05$ |  $-38842.51$ |  $-38368.45$ |  $11.34$
non int. GS |  $-$ |  $-44971.37$ |  $821.80$ |  $-38664.38$ |  $-37842.58$ |  $11.73$

$\nu=+2$, $\xi=20$ nm, $\Phi=\Phi_{0}$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
IVC |  $\prod_{\boldsymbol{k}}\prod_{\lambda}\frac{1}{\sqrt{2}}(d_{\boldsymbol{k}K\lambda\downarrow}+id_{\boldsymbol{k}K^{\prime}\lambda\downarrow})|+4\rangle$ |  $-44964.25$ |  $1065.69$ |  $-40642.57$ |  $-39576.89$ |  $11.34$
VP |  $\prod_{\boldsymbol{k}}d_{\boldsymbol{k}K+1\downarrow}d_{\boldsymbol{k}K-1\downarrow}|+4\rangle$ |  $-44964.25$ |  $1065.14$ |  $-40641.89$ |  $-39576.74$ |  $11.34$
non int. GS |  $-$ |  $-44971.37$ |  $1412.04$ |  $-40439.88$ |  $-39027.84$ |  $11.73$

Table F1:  Expectation values of the energy of candidates states for the
correlated insulators of MATBG at one magnetic flux quantum. We list the
states for fillings $\nu=0,\pm 2$ setting $\xi$ to $10$ and $20$ nm. The
values are in units of meV per unit cell, and $U$ is given in eV. The state
non int. GS is the ground state of the non interacting system. Coincidentally
the kinetic ground state at hole doping is given by the filled valence band of
the spin $\uparrow$, and at charge neutrality by the valence bands of both
spins, as can be seen in Fig. 3. $|0\rangle$ denotes the state with filled
remote bands and $|+4\rangle$ the $\nu=+4$ insulator. The Zeeman energy is
$\pm 1.535$ meV per electron and is not explicitly included.

$\nu=-2$, $\xi=10$ nm, $\Phi=0$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KA\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}B\uparrow})(d^{\dagger}_{\boldsymbol{k}KB\uparrow}-d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow})|0\rangle$ |  $-47981.11$ |  $642.52$ |  $-38976.71$ |  $-38334.19$ |  $12.07$
VP |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}KA\uparrow}d^{\dagger}_{\boldsymbol{k}KB\uparrow}|0\rangle$ |  $-47981.11$ |  $642.52$ |  $-38975.58$ |  $-38333.06$ |  $12.07$
TIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KA\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}B\uparrow})(d^{\dagger}_{\boldsymbol{k}KB\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow})|0\rangle$ |  $-47981.11$ |  $641.62$ |  $-38966.78$ |  $-38325.16$ |  $12.07$
QAH-IVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{\sqrt{2}}(d^{\dagger}_{\boldsymbol{k}KAs}+d^{\dagger}_{\boldsymbol{k}K^{\prime}Bs)}|0\rangle$ |  $-47981.11$ |  $641.62$ |  $-38969.92$ |  $-38328.30$ |  $12.41$
QAH-VP |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}KAs}|0\rangle$ |  $-47981.11$ |  $637.34$ |  $-38967.80$ |  $-38330.46$ |  $12.53$
VH |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}KA\uparrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow}|0\rangle$ |  $-47981.11$ |  $637.34$ |  $-38963.67$ |  $-38326.34$ |  $12.07$
non int. GS |  $-$ |  $-47982.83$ |  $581.99$ |  $-38598.70$ |  $-38016.71$ |  $12.34$

$\nu=-2$, $\xi=20$ nm, $\Phi=0$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KA\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}B\uparrow})(d^{\dagger}_{\boldsymbol{k}KB\uparrow}-d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow})|0\rangle$ |  $-47981.11$ |  $1232.78$ |  $-40791.19$ |  $-39558.42$ |  $12.07$
VP |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}KA\uparrow}d^{\dagger}_{\boldsymbol{k}KB\uparrow}|0\rangle$ |  $-47981.11$ |  $1232.78$ |  $-40790.08$ |  $-39557.30$ |  $12.07$
TIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KA\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}B\uparrow})(d^{\dagger}_{\boldsymbol{k}KB\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow})|0\rangle$ |  $-47981.11$ |  $1231.88$ |  $-40780.57$ |  $-39548.70$ |  $12.07$
QAH-IVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{\sqrt{2}}(d^{\dagger}_{\boldsymbol{k}KAs}+d^{\dagger}_{\boldsymbol{k}K^{\prime}Bs)}|0\rangle$ |  $-47981.11$ |  $1231.88$ |  $-40783.73$ |  $-39551.85$ |  $12.41$
QAH-VP |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}KAs}|0\rangle$ |  $-47981.11$ |  $1227.56$ |  $-40781.60$ |  $-39554.04$ |  $12.53$
VH |  $\prod_{\boldsymbol{k}}d^{\dagger}_{\boldsymbol{k}KA\uparrow}d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow}|0\rangle$ |  $-47981.11$ |  $1227.56$ |  $-40777.43$ |  $-39549.87$ |  $12.07$
non int. GS |  $-$ |  $-47982.83$ |  $1172.38$ |  $-40353.58$ |  $-39181.20$ |  $12.34$

$\nu=0$, $\xi=10$ nm, $\Phi=0$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KAs}+d^{\dagger}_{\boldsymbol{k}K^{\prime}Bs})(d^{\dagger}_{\boldsymbol{k}KBs}-d^{\dagger}_{\boldsymbol{k}K^{\prime}As})|0\rangle$ |  $-46395.22$ |  $-364.23$ |  $-38971.96$ |  $-39336.20$ |  $11.98$
VP |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}KAs}d^{\dagger}_{\boldsymbol{k}KBs}|0\rangle$ |  $-46395.22$ |  $-364.23$ |  $-38969.71$ |  $-39333.94$ |  $11.98$
TIVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KA\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}B\uparrow})(d^{\dagger}_{\boldsymbol{k}KB\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow})|0\rangle$ |  $-46395.22$ |  $-367.82$ |  $-38952.11$ |  $-39319.93$ |  $12.18$
SP |  $\prod_{\boldsymbol{k}}\prod_{\eta\sigma}d^{\dagger}_{\boldsymbol{k}\eta\sigma\uparrow}|0\rangle$ |  $-46395.22$ |  $-364.23$ |  $-38959.38$ |  $-39323.62$ |  $10.81$
VH |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}KAs}d^{\dagger}_{\boldsymbol{k}K^{\prime}As}|0\rangle$ |  $-46395.22$ |  $-384.95$ |  $-38945.89$ |  $-39330.85$ |  $12.67$
non int. GS |  $-$ |  $-46398.20$ |  $-364.92$ |  $-38439.73$ |  $-38804.65$ |  $12.00$

$\nu=0$, $\xi=20$ nm, $\Phi=0$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KAs}+d^{\dagger}_{\boldsymbol{k}K^{\prime}Bs})(d^{\dagger}_{\boldsymbol{k}KBs}-d^{\dagger}_{\boldsymbol{k}K^{\prime}As})|0\rangle$ |  $-46395.22$ |  $-368.00$ |  $-40784.72$ |  $-41152.72$ |  $11.98$
VP |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}KAs}d^{\dagger}_{\boldsymbol{k}KBs}|0\rangle$ |  $-46395.22$ |  $-368.00$ |  $-40782.49$ |  $-41150.49$ |  $11.98$
TIVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{2}(d^{\dagger}_{\boldsymbol{k}KA\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}B\uparrow})(d^{\dagger}_{\boldsymbol{k}KB\uparrow}+d^{\dagger}_{\boldsymbol{k}K^{\prime}A\uparrow})|0\rangle$ |  $-46395.22$ |  $-371.62$ |  $-40763.47$ |  $-41135.10$ |  $12.18$
SP |  $\prod_{\boldsymbol{k}}\prod_{\eta\sigma}d^{\dagger}_{\boldsymbol{k}\eta\sigma\uparrow}|0\rangle$ |  $-46395.22$ |  $-368.00$ |  $-40772.06$ |  $-41140.16$ |  $10.81$
VH |  $\prod_{\boldsymbol{k}}\prod_{s}d^{\dagger}_{\boldsymbol{k}KAs}d^{\dagger}_{\boldsymbol{k}K^{\prime}As}|0\rangle$ |  $-46395.22$ |  $-388.88$ |  $-40757.19$ |  $-41146.07$ |  $12.67$
non int. GS |  $-$ |  $-46398.20$ |  $-368.72$ |  $-40176.25$ |  $-40545.97$ |  $12.00$

$\nu=+2$, $\xi=10$ nm, $\Phi=0$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d_{\boldsymbol{k}KA\uparrow}+d_{\boldsymbol{k}K^{\prime}B\uparrow})(d_{\boldsymbol{k}KB\uparrow}-d_{\boldsymbol{k}K^{\prime}A\uparrow})|+4\rangle$ |  $-44809.34$ |  $543.00$ |  $-38954.64$ |  $-38411.64$ |  $11.88$
VP |  $\prod_{\boldsymbol{k}}d_{\boldsymbol{k}KA\uparrow}d_{\boldsymbol{k}KB\uparrow}|+4\rangle$ |  $-44809.34$ |  $543.00$ |  $-38953.51$ |  $-38410.51$ |  $11.88$
TIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d_{\boldsymbol{k}KA\uparrow}+d_{\boldsymbol{k}K^{\prime}B\uparrow})(d_{\boldsymbol{k}KB\uparrow}+d_{\boldsymbol{k}K^{\prime}A\uparrow})|+4\rangle$ |  $-44809.34$ |  $542.11$ |  $-38944.71$ |  $-38402.61$ |  $11.88$
QAH-IVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{\sqrt{2}}(d_{\boldsymbol{k}KAs}+d_{\boldsymbol{k}K^{\prime}Bs)}|+4\rangle$ |  $-44809.34$ |  $542.11$ |  $-38947.85$ |  $-38405.74$ |  $12.22$
QAH-VP |  $\prod_{\boldsymbol{k}}\prod_{s}d_{\boldsymbol{k}KAs}|+4\rangle$ |  $-44809.34$ |  $537.82$ |  $-38945.73$ |  $-38407.91$ |  $12.35$
VH |  $\prod_{\boldsymbol{k}}d_{\boldsymbol{k}KA\uparrow}d_{\boldsymbol{k}K^{\prime}A}|+4\rangle$ |  $-44809.34$ |  $537.82$ |  $-38941.61$ |  $-38403.78$ |  $11.88$
non int. GS |  $-$ |  $-44812.85$ |  $683.22$ |  $-38576.60$ |  $-37893.38$ |  $12.29$

$\nu=+2$, $\xi=20$ nm, $\Phi=0$
---
state |  wave function |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d_{\boldsymbol{k}KA\uparrow}+d_{\boldsymbol{k}K^{\prime}B\uparrow})(d_{\boldsymbol{k}KB\uparrow}-d_{\boldsymbol{k}K^{\prime}A\uparrow})|+4\rangle$ |  $-44809.34$ |  $1133.46$ |  $-40765.58$ |  $-39632.12$ |  $11.88$
VP |  $\prod_{\boldsymbol{k}}d_{\boldsymbol{k}KA\uparrow}d_{\boldsymbol{k}KB\uparrow}|+4\rangle$ |  $-44809.34$ |  $1133.46$ |  $-40764.47$ |  $-39631.01$ |  $11.88$
TIVC |  $\prod_{\boldsymbol{k}}\frac{1}{2}(d_{\boldsymbol{k}KA\uparrow}+d_{\boldsymbol{k}K^{\prime}B\uparrow})(d_{\boldsymbol{k}KB\uparrow}+d_{\boldsymbol{k}K^{\prime}A\uparrow})|+4\rangle$ |  $-44809.34$ |  $1132.55$ |  $-40754.96$ |  $-39622.41$ |  $11.88$
QAH-IVC |  $\prod_{\boldsymbol{k}}\prod_{s}\frac{1}{\sqrt{2}}(d_{\boldsymbol{k}KAs}+d_{\boldsymbol{k}K^{\prime}Bs)}|+4\rangle$ |  $-44809.34$ |  $1132.55$ |  $-40758.11$ |  $-39625.56$ |  $12.22$
QAH-VP |  $\prod_{\boldsymbol{k}}\prod_{s}d_{\boldsymbol{k}KAs}|+4\rangle$ |  $-44809.34$ |  $1128.23$ |  $-40755.98$ |  $-39627.75$ |  $12.35$
VH |  $\prod_{\boldsymbol{k}}d_{\boldsymbol{k}KA\uparrow}d_{\boldsymbol{k}K^{\prime}A}|+4\rangle$ |  $-44809.34$ |  $1128.23$ |  $-40751.81$ |  $-39623.58$ |  $11.88$
non int. GS |  $-$ |  $-44812.85$ |  $1273.27$ |  $-4032.06$ |  $-39058.79$ |  $12.29$

Table F2: Expectation values of the energy of candidates states for the correlated insulators of MATBG at zero magnetic field. We list the states for fillings $\nu=-2,0,+2$ setting $\xi$ to $10$ and $20$ nm. The energies are in units of meV per unit cell, and the value of $U$ is given in eV. The state non int. GS is the ground state of the non interacting system. QAH (Quantum anomalous Hall) states have Chern number 2 and are stabilized by small out of plane magnetic fields[8, 9]. $|0\rangle$ denotes the state with the filled remote bands and $|+4\rangle$ the $\nu=+4$ insulator. $\nu$, ($\epsilon$, $U$) |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
---|---|---|---|---|---
$-2$, $(10,4$ eV) |  $-48147.31$ |  $228.05$ |  $-38903.47$ |  $-38675.42$ |  $11.39$
$-2$, $(50,0.5$ eV) |  $-48149.36$ |  $222.52$ |  $-38842.95$ |  $-38620.43$ |  $11.37$
$0$, $(10,4$ eV) |  $-46552.57$ |  $-340.91$ |  $-38890.94$ |  $-39231.85$ |  $10.90$
$0$, $(50,0.5$ eV) |  $-46554.48$ |  $-341.91$ |  $-38845.54$ |  $-39187.45$ |  $10.76$
$+2$, $(10,4$ eV) |  $-44959.04$ |  $319.40$ |  $-38846.62$ |  $-38527.22$ |  $11.08$
$+2$, $(50,0.5$ eV) |  $-44964.44$ |  $476.25$ |  $-38852.51$ |  $-38376.26$ |  $11.34$
Table F3:  Expectation values of the energy of the ground states for
$\boldsymbol{B=26.5}$ T. The fillings are $\nu=0,\pm 2$ and the interaction
parameters $\epsilon=10$, $U=4$ eV and $\epsilon=10$, $U=0.5$ eV. The gate
distance is set to $\xi=10$ nm. The values are in units of meV per unit cell,
and $U$ is given in eV. The Zeeman energy is $\pm 1.535$ meV per electron and
is not explicitly included.

$\nu=-2$. $\epsilon=10$, $U=4$ eV
---
order |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $-47980.57$ |  $631.44$ |  $-39002.21$ |  $-38370.76$ |  $12.07$
VP |  $-47981.24$ |  $613.44$ |  $-38958.61$ |  $-38345.17$ |  $12.05$

$\nu=-2$. $\epsilon=50$, $U=0.5$ eV
---
order |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $-47981.42$ |  $602.60$ |  $-38951.76$ |  $-38349.17$ |  $12.12$
VP |  $-47981.58$ |  $600.96$ |  $-38936.46$ |  $-38335.50$ |  $12.03$

$\nu=0$. $\epsilon=10$, $U=4$ eV
---
order |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $-46393.94$ |  $-364.72$ |  $-39036.56$ |  $-39401.29$ |  $11.99$
SP |  $-46394.61$ |  $-364.36$ |  $-38980.15$ |  $-39344.50$ |  $10.88$
VP |  $-46394.61$ |  $-364.36$ |  $-38989.93$ |  $-39354.29$ |  $11.98$

$\nu=0$. $\epsilon=50$, $U=0.5$ eV
---
order |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $-46394.82$ |  $-365.14$ |  $-39003.10$ |  $-39368.25$ |  $11.98$
SP |  $-46395.38$ |  $-364.90$ |  $-38953.35$ |  $-39318.25$ |  $10.84$
VP |  $-46395.38$ |  $-364.90$ |  $-38963.48$ |  $-39328.38$ |  $11.98$

$\nu=+2$. $\epsilon=10$, $U=4$ eV
---
order |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $-44808.71$ |  $545.40$ |  $-38986.29$ |  $-38440.89$ |  $11.89$
VP |  $-44809.02$ |  $541.18$ |  $-38961.14$ |  $-38419.95$ |  $11.88$

$\nu=+2$. $\epsilon=50$, $U=0.5$ eV
---
order |  kinetic |  $\epsilon\times$Hartree |  $\epsilon\times$Fock |  $\epsilon\times$Coulomb |  Hubbard$/U$
KIVC |  $-44809.07$ |  $544.39$ |  $-38973.52$ |  $-38429.13$ |  $11.89$
VP |  $-44809.34$ |  $543.00$ |  $-38953.51$ |  $-38410.51$ |  $11.88$

Table F4:  Expectation values of the energy of the self-consistent states for
$\boldsymbol{B=0}$ T. The fillings are $\nu=0,\pm 2$ and the interaction
parameters $\epsilon=10$, $U=4$ eV and $\epsilon=10$, $U=0.5$ eV. The gate
distance is set to $\xi=10$ nm. The states are labeled by the dominant order
parameter. The values are in units of meV per unit cell, and $U$ is given in
eV.