# Particle-hole asymmetric phases in doped twisted bilayer graphene

Run Hou  Department of Physics and Astronomy, Rice University, Houston, TX
77005, USA  Shouvik Sur  Department of Physics and Astronomy, Rice University,
Houston, TX 77005, USA  Lucas K. Wagner  Department of Physics, University of
Illinois, Urbana-Champaign, USA  Andriy H. Nevidomskyy  Department of Physics
and Astronomy, Rice University, Houston, TX 77005, USA

###### Abstract

Twisted bilayer graphene (TBG) has emerged as a paradigmatic platform for
exploring the interplay between strong interactions in a multi-band system
with nearly flat bands, while offering unprecedented control over the filling
fraction of electron/hole carriers. Despite much theoretical work, developing
a comprehensive ab initio model for this system has proven challenging due to
the inherent trade-off between accurately describing the band structure and
incorporating the interactions within the Hamiltonian, particularly given the
topological obstruction – so-called fragile topology – to the description of
the model in terms of localized symmetric Wannier functions within the flat
band manifold. Here, we circumvent this obstruction by using an extended
8-orbital model, for which localized Wannier orbitals have been formulated by
Carr et al. [1]. We constructed an extended multi-orbital Hubbard model, and
performed Hartree-Fock (HF) calculations to explore its phase diagram across
commensurate fillings from -3 to 3. We found several nearly-degenerate
insulating states at charge neutrality, all of which exhibit orbital orders.
Crucially, TBG near magic angle is known to be particle-hole asymmetric, which
is naturally captured by the single-particle band structure of our model and
is reflected in the distinction between the symmetry broken states obtained at
electron and hole dopings away from the charge neutral point. At filling -1
and +2, quantum anomalous hall states are obtained, while for the rest of the
integer fillings away from charge neutrality, we found the system to realize
metallic states with various orbital, valley and spin orderings. We also
observed that most of the Hartree–Fock ground states exhibit a generalized
valley Hund’s-like rule, resulting in valley polarization. Importantly, we
show that the incorporation of the intra-valley and inter-valley exchange
interactions is crucial to properly stabilize the ordered symmetry-broken
states. In agreement with experiments, we find significant particle-hole
asymmetry, which underscores the importance of using particle-hole asymmetric
models.

##  I Introduction

The theoretical prediction of nearly flat bands [2] and subsequent discovery
of correlated insulating and superconducting phenomena [3, 4] in magic-angle
twisted bilayer graphene (TBG) has sparked significant interest in the
interplay between topology and strong correlations, leading to extensive
theoretical [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59,
60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78,
79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97,
98, 99, 100, 101, 102, 103] and experimental studies [104, 105, 106, 107, 108,
109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123,
124, 125, 126, 127, 128, 129, 130, 131, 132] in the field of condensed matter
physics. The construction of effective models and the development of theories
to understand this system have emerged as crucial endeavors in contemporary
research.

In order to establish models that are both simple and capable of capturing the
essential features of magic-angle TBG without sacrificing the ab initio
perspective, various approaches have been explored. These include models based
on maximally localized Wannier functions (MLWF) [6, 7, 1], atomic tight-
binding [5, 93], the Bistritzer–MacDonald (BM) model [2], and topological
heavy fermion (THF) model [96, 95, 97, 101]. However, developing a
comprehensive model has proven challenging due to the inherent trade-off
between accurately describing the band structure and incorporating the
interactions within the Hamiltonian, particularly given the fragile
topological nature [8, 9] of the magic-angle TBG system.

The BM model, despite its simplicity in constructing the non-interacting band
structure of magic-angle TBG, is approximate due to its inherent k$\cdot$p
nature and moreover, encounters complexity when incorporating interactions by
projecting out remote high-energy bands. This renders the traditional methods
for studying strongly correlated systems less applicable. Alternatively, the
THF model for TBG has been developed and has yielded promising results [99,
100, 101]. However, MLWF models offer another viable option as they provide
accurate descriptions of the band structure and localized interactions, albeit
with a higher number of orbitals. Regrettably, the MLWF models as alternative
choices for magic-angle TBG have not been extensively investigated.

Therefore, in this study, we adopt an 8-orbital model of magic-angle TBG from
Carr et al. [1] and employ Hartree–Fock (HF) methods to explore its
properties. Unlike the belief that Wannierized models are challenging to use
in calculations, we found the HF method to be straightforward, as it does not
involve complicated form factors such as those arising from projecting the
interactions onto the BM model. The 8-orbital model also has the potential for
future use in more sophisticated strong correlation methods such as DMFT, to
calculate the states that cannot be captured by the HF self-consistent
calculations. Although previous studies [69] had considered the 8-band model,
those calculations were done at the Hartree level (without the Fock terms),
and did not find any interaction-driven insulating states.

In this work, we first used the Wannier orbitals from the 8-orbital model to
numerically evaluate the strength of the electron repulsion integral matrix
elements, enabling us to construct a suitably extended Hubbard model including
both the direct and exchange interactions. Subsequently, we incorporated the
complete Hartree and Fock terms to study this model. We identified several
converged ordered states for each filling. These symmetry-broken ordered
states originate from the valley, spin, and orbital degrees of freedom. Our
results compare favourably with the experimental findings [105, 106, 107, 108,
112, 114, 115, 122, 125, 104, 113, 131, 132], including the quantum anomalous
Hall (QAH) states found at certain fillings. The particle-hole asymmetry is
naturally present in our 8-orbital model, in contrast to the BM model, and is
consistent with the experimentally observed asymmetry between electron- and
hole-doped side of the phase diagram. A detailed comparison between our
results and experimental findings suggests that the symmetry-breaking and gap
formation mechanisms depend on the filling fraction and the specific device
configurations used in the experiments. Our HF calculations underscore the
importance of incorporating particle-hole asymmetric models into theoretical
considerations.

We briefly discuss the mechanisms of spontaneous symmetry breaking here. Some
previous studies had used (approximate) valley-spin SU(4) symmetry to discuss
various symmetry-broken states in TBG systems [13, 103]. The atomistic model,
such as the one employed here [1], does not actually possess full SU(4)
symmetry – instead the kinetic energy terms have
$\text{U(2)}_{+}\times\text{U(2)}_{-}$ symmetry where $\pm$ label the two
graphene valleys, with the
$\text{U(2)}=\text{U(1)}_{\text{charge}}\times\text{SU(2)}_{\text{spin}}$
symmetry separately in each valley. Nevertheless, there are close parallels
with the spontaneously broken valley and spin symmetries that are also present
in SU(4)-type models [103], or in the generalized Stoner picture [113]. We
summarize the various possible spontaneous symmetry-broken states in Table 1.

Importantly, we found that the atomistic model also contains information about
orbital degrees of freedom ($p_{+}$ and $p_{-}$ orbitals in the 8-band model
[1]), which is absent from the BM-type analysis and from the above symmetry
considerations. These orbitals are energetically degenerate in the non-
interacting limit, and this additional symmetry can be spontaneously broken in
the presence of interactions. Indeed, we found this to be the case in nearly
all of the ordered states we identified, with orbital degeneracy spontaneously
lifted. In fact, we find the ground state to be a combination of orbital
ordering and one of the symmetry-broken patterns summarized in Table 1,
depending on the filling fraction.

States | Symmetry
---|---
Normal state | $\text{U(2)}_{+}\times\text{U(2)}_{-}$
(a)Valley polarized (VP) state | Breaking time-reversal symmetry $\mathcal{T}$
(b)Orbital polarized (OP) state | Breaking discretized symmetries such as $C_{2z}\mathcal{T}$ and $C_{2x}$
(c)Spin polarized (SP) state |  $\text{U(1)}_{+}\times\text{U(2)}_{-}$ or $\text{U(1)}_{+}\times\text{U(1)}_{-}$
(d)Inter-valley coherent (IVC) state | $\text{SU(2)}_{+}\times\text{SU(2)}_{-}\times\text{U(1)}$
Table 1: Examples of simplified symmetry-broken orders. When the interaction
strength is large enough, the system develops orders. Realistic models produce
more complicated orders. (a) The VP state breaks time-reversal symmetry. (b)
The OP state breaks the orbital discretized symmetry. (c) The SP state breaks
the spin SU(2) symmetry. (d) The IVC state breaks the
$\text{U(1)}_{+}\times\text{U(1)}_{-}$ symmetry.

This paper is organized as follows. In Sec. II, we review the non-interacting
8-orbital tight-binding model from previous work and numerically evaluate the
interaction parameters. In Sec. III, we discuss the numerical results, where
we present the order parameters, symmetry breaking, and associated Chern
numbers, if relevant, of the obtained HF solutions. Importantly, we also show
that the role of exchange interactions is crucial to correctly capture the
correlated insulating states of the TBG. We compare our findings to the
experiments in Sec. IV, before summarizing our conclusions.

##  II Minimal interacting tight-binding model for MATBLG

In this section, we review the lattice structure and construct the interacting
8-orbital model of magic-angle TBG.

###  II.1 Lattice Structure and Noninteracting Hamiltonian

To begin with, we review the lattice structure of TBG 8-band model derived in
[1]. There are 8 Wannier orbitals per spin and per valley inside a given
triangular unit cell of the moiré lattice as shown in Fig. 1. The orbital
indices, corresponding Wannier orbital centers, and orbital symmetries are
shown in TABLE 2. This 8-orbital tight binding model has correct crystalline
symmetries, avoiding the Wannier obstruction resulting from fragile topology
if only a subset of bands are included [9].

Figure 1: Lattice structure of the 8-orbital model. The Bravais lattice of the system is a triangular lattice. Orange dots in triangular lattice represent orbital 1/2 with $p_{\pm}$ symmetries and orbital 3 with $s$ symmetry. Green triangles in honeycomb lattice represent $p_{z}$ orbitals labeled by orbital indices 4 and 5. Purple diamonds in the kagome lattice represent orbitals 6,7, and 8 with $s$ symmetry. The primitive vectors are $\vec{a}_{1}$ and $\vec{a}_{2}$. Orbital indices | Region | Wyckoff position | Symmetry
---|---|---|---
1,2 | AA | $1a$ | $p_{+},p_{-}$
3 | AA | $1a$ | $s$
4,5 | AB | $2b$ | $p_{z}$
6,7,8 | DW | $3c$ | $s$
Table 2: Orbital indices and corresponding regions and symmetries. The AA
region is located at crossing of triangular lattice usually denoted as $1a$.
The AB region is located at the center of triangles denoted as $2b$. And the
DW region is located at the centers of triangular edges denoted as $3c$.

The Wannier orbitals and noninteracting Hamiltonian is obtained from the
maximally localized Wannier functions (MLWF) method, and the resulting non-
interacting tight-binding model is taken from Carr et al. [1]:

$H_{K}=\sum_{ij,ab,\tau
s}t^{\tau}_{ab}(\bm{R}_{i}-\bm{R}_{j})c^{\dagger}_{a\tau
s}(\bm{R}_{i})c_{b\tau s}(\bm{R}_{j}),$ (1)

where $c^{\dagger}_{a\tau s}(\bm{R}_{i})$ is the electron creation operator in
the $a^{\text{th}}$ Wannier orbital ($a=1,\ldots 8$) positioned within the
Moiré unit cell labeled with a lattice vector $\bm{R}_{i}$. The subscript
$\tau=\pm$ denotes graphene valley indices, and $s=\uparrow,\downarrow$ are
the electron spin indices. The hopping parameters
$t^{\tau}_{ab}(\bm{R_{i}}-\bm{R_{j}})$ are fitted from ab initio
$\bm{k}\cdot\bm{p}$ model with lattice relaxation. In the present work, we
choose as a basis the non-interacting Hamiltonian at twist angle
$\theta=1.10\degree$.

We Fourier transform the electron creation/annihilation operators into
momentum space using the periodic gauge, meaning that the Fourier phases do
not depend on specific positions of Wannier centers but rather only on the
Moiré cell coordinate $\bm{R_{i}}$:

$c^{\dagger}_{a\tau
s}(\bm{R}_{i})=\frac{1}{\sqrt{N_{k}}}\sum_{\bm{k}}c^{\dagger}_{a\tau
s}(\bm{k})e^{i\bm{R}_{i}\cdot\bm{k}}.$ (2)

The momentum-space non-interacting Hamiltonian is then written in the form

$H_{K}=\sum_{ab\tau s}t^{\tau}_{ab}(\bm{k})c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau s}(\bm{k}).$ (3)

The periodic gauge has the convenient property that
$t^{\tau}_{ab}(\bm{k}+\bm{G})=t^{\tau}_{ab}(\bm{k})$, where $\bm{G}$ is the
Moiré reciprocal lattice vector.

The noninteracting Hamiltonian has several important symmetries. The time-
reversal symmetry relates the hoppings in the two valleys via

$t^{-\tau}_{ab}(\bm{k})=t^{\tau*}_{ab}(-\bm{k}).$ (4)

The Hamiltonian is also invariant under operations of $C_{2z}\mathcal{T}$,
$C_{3z}$, and $C_{2x}$ symmetries as demonstrated in Appendix A. For the
internal symmetries, it is clear that the system preserves
$\text{U(2)}_{+}\times\text{U(2)}_{-}$ symmetry, which means for each valley
$+/-$ the system has spin SU(2) symmetry and preserves the total particle
number $N_{+}+N_{-}$ and the particle number difference $N_{+}-N_{-}$.

In Fig. 2, we showed the non-interacting band structure of the single valley
8-orbital TBG model along with the projected weights of orbitals that have
different orbital symmetries. It clearly shows that the central bands are
mainly composed by $p_{\pm}$ orbitals except for the points near $\Gamma$
point.

Figure 2: The tight-binding band structure of non-interacting Hamiltonian. (a)
The dark orange fat bands show the projected weights of $AA_{p_{\pm}}$
orbitals. (b) The orange fat bands show the projected weights of $AA_{s}$
orbital. (c) The fat bands for $AB_{p_{z}}$ orbitals. (d) The fat bands for
$DW_{s}$ orbitals. (e) The density of states of the whole band structure.

###  II.2 Interacting Hamiltonian

The full Hamiltonian is composed of the noninteracting Hamiltonian $H_{K}$ and
interactions $H_{I}$:

$H=H_{K}+H_{I}.$ (5)

Now, we elaborate how to define the interactions in our multi-orbital model.
The most general form of the interacting Hamiltonian that preserves spin SU(2)
symmetry (since the spin-orbit coupling is negligible in graphene) is as
follows:

$H_{I}=\frac{1}{2}\sum_{\tau_{j},a_{j},\mathbf{R}_{j},s,s^{\prime}}V^{\tau_{1}\tau_{2}\tau_{3}\tau_{4}}_{a_{1}a_{2}a_{3}a_{4}}(\bm{R}_{1},\bm{R}_{2},\bm{R}_{3},\bm{R}_{4})\times\\\
c^{\dagger}_{a_{1}\tau_{1}s}(\bm{R}_{1})c^{\dagger}_{a_{2}\tau_{2}s^{\prime}}(\bm{R}_{2})c_{a_{3}\tau_{3}s^{\prime}}(\bm{R}_{3})c_{a_{4}\tau_{4}s}(\bm{R}_{4}).$
(6)

The four-center electron repulsion integrals
$V^{\tau_{1}\tau_{2}\tau_{3}\tau_{4}}_{a_{1}a_{2}a_{3}a_{4}}(\bm{R}_{1},\bm{R}_{2},\bm{R}_{3},\bm{R}_{4})$
were evaluated using the real-space Wannier functions, and we found the
oscillating factor $e^{i\bm{q}_{V}\cdot(\bm{r}-\bm{r}^{\prime})}$ suppresses
the absolute value of the integral when $\tau_{2}\neq\tau_{3}$ and
$\tau_{1}\neq\tau_{4}$, because the momentum $\bm{q}_{V}$ between the graphene
two valleys is large as explained in Appendix B. With the benefit of MLWF, it
is reasonable to only evaluate the 2-center interactions that have the forms
like
$V^{\tau\tau^{\prime}\tau^{\prime}\tau}_{a_{1}a_{2}a_{3}a_{4}}(\bm{R}_{1},\bm{R}_{2},\bm{R}_{2},\bm{R}_{1})$
and
$V^{\tau\tau^{\prime}\tau^{\prime}\tau}_{a_{1}a_{2}a_{3}a_{4}}(\bm{R}_{1},\bm{R}_{2},\bm{R}_{1},\bm{R}_{2})$.

Among these interaction channels, we found the strongest interaction channels
are density-density channels $H_{C}$, with the exchange channels $H_{X}$
playing a key role away from half-filling. This motivates the form of the
Hamiltonian

$H=H_{K}+H_{C}+H_{X}.$ (7)

The density-density interacting Hamiltonian takes the form

$\begin{split}H_{C}=\frac{1}{2}\sum_{ij;ab;\tau\tau^{\prime}ss^{\prime}}U_{ab}(\bm{R}_{i}-\bm{R}_{j})\times\\\
c^{\dagger}_{a\tau
s}(\bm{R}_{j})c^{\dagger}_{b\tau^{\prime}s^{\prime}}(\bm{R}_{i})c_{b\tau^{\prime}s^{\prime}}(\bm{R}_{i})c_{a\tau
s}(\bm{R}_{j}),\end{split}$ (8)

which is an extended Hubbard model with long-range interaction. The density-
density repulsive interaction $U_{ab}(\bm{R}_{i}-\bm{R}_{j})$ does not depend
on the choice of valleys and is explicitly evaluated assuming a screened
single-gated Coulomb potential $V_{SG}(r)$ often used in modeling of TBG:

$V_{SG}(r)=\frac{1}{4\pi\epsilon\epsilon_{0}}\left(\frac{1}{r}-\frac{1}{\sqrt{r^{2}+(2d)^{2}}}\right).$
(9)

The relative dielectric constant $\epsilon=12$ was set to capture the
screening effect from remote upper and lower bands [69]. The distance between
the sample and capacitor was set to $d=10\text{nm}$, which is comparable with
the magic-angle TBG Moiré lattice length and is close to realistic
experimental settings.

The resulting density-density interaction is

$\begin{split}U_{ab}(\bm{R}_{i}-\bm{R}_{j})&=\int\mathrm{d}\bm{r}^{2}\mathrm{d}{\bm{r}^{\prime}}^{2}V_{SG}(|\bm{r}-\bm{r}^{\prime}|)\times\\\
&\sum_{Y,Y^{\prime}}|\mathcal{W}^{Y}_{\bm{R}_{i}a}(\bm{r})|^{2}|\mathcal{W}^{Y^{\prime}}_{\bm{R}_{j}b}(\bm{r}^{\prime})|^{2},\end{split}$
(10)

where $\mathcal{W}^{Y}_{\bm{R}_{i}a}(\bm{r})$ is the Wannier function for the
orbital $a$ at the unit cell located at $\bm{R}_{i}$ for graphene-sublattice
$Y$[1].

We represent the general trend in the density-density interactions between
selected orbitals in Fig. 3. Noting that the lattice vector $\bm{R}_{i}$ is
not equivalent to the position of the Wannier center $\bm{R}_{i}+\bm{r}_{a}$,
we elected to plot the interaction integrals in Fig. 3 as a function of the
distance $\Delta r=|\bm{R}_{i}+\bm{r}_{a}-\bm{R}_{j}-\bm{r}_{b}|$ between the
two orbital centers in order to compare with the bare single-gated
interaction: $\tilde{U}_{ab}(\Delta r)=U_{ab}(\bm{R}_{i}-\bm{R}_{j})$. The
strongest onsite density-density interaction is around 35 meV, which is
comparable with the bandwidth of the central bands. In practice, we chose the
density-density repulsive interactions to be extended to next-nearest neighbor
unit cells. At longer distances, $V_{SG}(r)$ becomes negligible and is smaller
than the exchange interactions that we will introduce next.

Figure 3: Single-gated screened density-density channel interaction strengths
as a function of the distance between orbital centers $\Delta r$. The
horizontal axis is rescaled with Moiré lattice length $\lambda$. Data points
with different shapes correspond to different orbital pairs. Figure 4: (a) The
exchange interaction arrangement. The short-range exchange interactions
between $p_{\pm}$ and other orbitals within a unit cell are taken into
account. (b) The exchange interactions in the same valley
$X^{\tau\tau}_{ab}(\Delta\bm{R})$ and in different valleys
$X^{\tau\bar{\tau}}_{ab}(\Delta\bm{R})$. We found
$X^{\tau\bar{\tau}}_{ab}(\Delta\bm{R})$ is negligible in comparison to the
exchange interactions in the same valley.

The second part of the interacting Hamiltonian includes the exchange
interaction matrix elements $X_{ab}^{\tau\tau^{\prime}}$ which are obtained as
integrals over the corresponding Wannier orbitals, as detailed in Appendix B:

$\begin{split}&H_{X}=\frac{1}{2}\sum_{a,b,\tau,\tau^{\prime},s,s^{\prime}}\sum_{ij}X^{\tau\tau^{\prime}}_{ab}(\bm{R}_{i}-\bm{R}_{j})\times\\\
&:c^{\dagger}_{a\tau s}(\bm{R}_{i})c_{b\tau
s}(\bm{R}_{j})c^{\dagger}_{b\tau^{\prime}s^{\prime}}(\bm{R}_{j})c_{a\tau^{\prime}s^{\prime}}(\bm{R}_{i}):.\end{split}$
(11)

Rather than keep the full complexity of the matrix
$X_{ab}^{\tau\tau^{\prime}}$, we observed that the dominant effect at
experimentally accessible densities arises from the $p_{\pm}$ orbitals (which
constitute the flat bands in magic-angle TBG) interacting among themselves and
with the other orbitals, as illustrated in Fig. 4(a). In this study, we
therefore dropped the matrix elements $X_{ab}^{\tau\tau^{\prime}}$ whenever
both indices $a$ and $b$ $\not\in\\{1,2\\}$.

Furthermore, while the exchange interaction $X_{ab}^{\tau\tau^{\prime}}$
depends on the valley indices (see Appendix B for details), the computed
values shown in Fig. 4(b) demonstrate that the exchange interaction between
different valleys $X^{\tau\bar{\tau}}_{ab}$, is significantly smaller in
magnitude compared to that within the same valley $X^{\tau\tau}_{ab}$. For the
sake of simplicity, we assumed $X^{\tau\bar{\tau}}_{ab}=0$. The largest
exchange interaction is approximately $2.7$ meV, corresponding to
$X^{\tau\tau}_{13}$ and $X^{\tau\tau}_{23}$, independent of $\tau$, that is,
the exchange between the $p_{\pm}$ and $s$-orbital centered in the AA region
of the same Moiré supercell. We noted that previous studies [69, 102] of this
model did not consider exchange interactions due to its small magnitude
compared with the on-site density-density interaction. In contrast, in this
study we found that exchange interactions have pronounced effects on the
interacting model and should not be ignored.

##  III Numerical results

$\nu$ | Order | Insulator | Order parameter |  $C_{2z}\mathcal{T}$ symmetry | Chern number $|C|$
---|---|---|---|---|---
0 | OP (orbital-polarized) | ✓ | $\left<\Sigma_{3}\right>$ | $\crossproduct$ | 0
0 | VOP (valley-orbital polarized) | ✓ | $\left<\tau_{3}\otimes\Sigma_{3}\right>$ | $\crossproduct$ | 4
0 | VOP+OP | ✓ | $\left<\tau_{3}\otimes\Sigma_{3}\right>_{\uparrow},\left<\Sigma_{3}\right>_{\downarrow}$ | $\crossproduct$ | 2
0 | VP (valley-polarized) | $\crossproduct$ | $\left<\tau_{3}\right>$ | ✓ | N/A
0 | SP (spin-polarized) | $\crossproduct$ | $\left<s_{3}\right>$ | $\crossproduct$ | N/A
0 | IVC (inter-valley coherent) | $\crossproduct$ | $\left<\tau_{1,2}\right>$ | ✓ | N/A
2 | VP | ✓ | $\left<\tau_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | 2
2 | SP | ✓ | $\left<s_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | 0
2 | IVC (inter-valley coherent)+VP | $\crossproduct$ | $\left<\tau_{3}\right>_{\uparrow},\left<\tau_{1,2}\right>_{\downarrow},\left<\Sigma_{3}\right>$ | $\crossproduct$ | N/A
-2 | VP | $\crossproduct$ | $\left<\tau_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | 2
-2 | SP | $\crossproduct$ | $\left<s_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | 0
-2 | OP | $\crossproduct$ | $\left<\Sigma_{3}\right>$ | $\crossproduct$ | 0
1 | VPSP | $\crossproduct$ | $\left<s_{3}\right>,\left<\tau_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | N/A
1 | IVCSP | $\crossproduct$ | $\left<s_{3}\right>,\left<\tau_{1,2}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | N/A
-1 | VPSP | ✓ | $\left<s_{3}\right>,\left<\tau_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | 3
-1 | IVCSP | $\crossproduct$ | $\left<s_{3}\right>,\left<\tau_{1,2}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | N/A
-3 | VPSP | $\crossproduct$ | $\left<s_{3}\right>,\left<\tau_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | N/A
-3 | IVCSP | $\crossproduct$ | $\left<s_{3}\right>,\left<\tau_{1,2}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | N/A
3 | VPSP | $\crossproduct$ | $\left<s_{3}\right>,\left<\tau_{3}\right>,\left<\Sigma_{3}\right>$ | $\crossproduct$ | 1
Table 3: Summary of properties of converged states at different filling $\nu$
obtained from HF calculations. The third column labels whether the state is an
insulator or not. The fourth column shows the order parameters of the
corresponding state. The fifth column shows whether the $C_{2z}\mathcal{T}$
symmetry is broken. The last column shows the absolute value of the Chern
number. We use Pauli matrices $s$, $\tau$, and $\Sigma$ to label spin, valley,
and the $p_{\pm}$ degrees of freedom. The notation
$\left<\hat{O}\right>_{\uparrow/\downarrow}$ represents the order in the spin-
up/down species. For IVC states, $\left<\tau_{1,2}\right>$ refers to any order
at the $U(1)$ circle formed by $\tau_{1}$ and $\tau_{2}$. We calculated the
Chern number for all insulating states, also including some states that are
not fully gapped. This is done for reference purposes, as the bands are still
separable.

###  III.1 Order parameters and symmetry breaking

In this section, we discuss the ground state candidates at the filling
fractions $\nu=0,\pm 1,\pm 2,\pm 3$ relevant to the experiment where the
correlated insulating behavior has been observed [3, 110, 105, 108, 109, 104,
107, 106, 125, 111, 124, 122, 120, 121, 113, 114, 123, 115, 132]. Because the
8-orbital model does not have particle-hole symmetry (unlike the approximate
BM model), we found the resulting states not to be the same for the particle-
doped and the hole-doped sides of the phase diagram. We performed the multi-
band Hartree-Fock (HF) calculations of the Hamiltonian Eq. (7) up to
$15\crossproduct 15$ discretized points in momentum space for a fixed twist
angle $\theta=1.10\degree$. The detailed derivation of the HF theory is
presented in Appendix C.

A closer analysis of the converged results shows that all symmetry-broken
states can be characterized by spin, valley, and orbital degrees of freedom.
The most general order parameter is the mean value of fermion bilinears
$\hat{O}(\bm{k})=\sum_{a\tau s,b\tau^{\prime}s^{\prime}}c^{\dagger}_{a\tau
s}(\bm{k})\mathcal{O}_{a\tau
s,b\tau^{\prime}s^{\prime}}c_{b\tau^{\prime}s^{\prime}}(\bm{k})$:

$\left<\hat{O}\right>=\frac{1}{N_{k}}\sum_{\bm{k}}\Tr[\mathcal{O}P(\bm{k})],$
(12)

where $P_{a\tau
s,b\tau^{\prime}s^{\prime}}(\bm{k})=\left<c^{\dagger}(\bm{k})_{a\tau
s}c(\bm{k})_{b\tau^{\prime}s^{\prime}}-\frac{1}{2}\delta_{ab}\delta_{\tau\tau^{\prime}}\delta_{ss^{\prime}}\right>$
is the one-body reduced density matrix, and $N_{k}$ is the total number of
momentum points. While the orbital symmetry breaking is in principle possible
for all 8 orbitals, we found that all the orbital ordering is limited to the
$p_{\pm}$ orbitals that constitute the nearly flat bands in magic-angle TBG.
We thus limited the orbital order to this subset and use the Pauli matrix
$\Sigma$ to label the orders formed by the $p_{\pm}$ orbitals. The ordering of
valley and spin degrees of freedom is labeled by the Pauli matrices $\tau$ and
$s$, respectively.

While the $C_{2z}\mathcal{T}$ symmetry is preserved by the tight-binding
Hamiltonian, the interactions can result in the symmetry being spontaneously
broken. Denoting the $C_{2z}\mathcal{T}$ operator by $g$ for brevity, the
corresponding symmetry breaking strength is given by (see Appendix D for
derivation):

$\displaystyle\mathcal{D}(\bm{k})=|\left<g^{-1}c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})g\right>-\left<c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})\right>|.$ (13)

When summed over all $\mathbf{k}$-points as in Eq. (12), it measures the
difference between the expectation value of the order parameter and its
$C_{2z}\mathcal{T}$ symmetry transformed value. The $C_{2z}\mathcal{T}$
symmetry-breaking strength is zero if the order parameter is invariant under
the $C_{2z}\mathcal{T}$ symmetry.

Furthermore, we calculated the Chern number $C$ which characterizes quantum
anomalous Hall states in insulators. The Chern number is calculated by
measuring the winding of the wavefunction along the Wilson loop as shown in
Appendix F.

Table 3 summaries the main results. we found that the order parameters of
converged states at all fillings are generally characterized by several broken
symmetries: the spin polarized (SP) state with $\left<s_{3}\right>\neq 0$, the
valley polarized (VP) state with $\left<\tau_{3}\right>\neq 0$, the inter-
valley coherent (IVC) state with $\left<\tau_{1,2}\right>\neq 0$, and the
orbital polarized (OP) state with $\left<\Sigma_{3}\right>\neq 0$. Depending
on the filling, we found both the correlated metallic states and insulators
with broken symmetries. Certain semimetallic states near charge neutrality
preserve the $C_{2z}\mathcal{T}$ symmetry, while all other orders break the
$C_{2z}\mathcal{T}$ symmetry.

###  III.2 Hartree–Fock results at all integer fillings

Figure 5: Relative energies of converged ordered states. and schematic phase
diagram of $\theta=1.10^{\degree}$ 8-orbital model. The inset shows details of
relative energies of competing orders. Figure 6: Ground state candidates at
charge neutrality for TBG 8-orbital model. The color bar shows the valley
polarization strength for each k point. (a) Orbital polarized state (OP). (b)
Valley-orbital polarized state (VOP). (c) Valley-orbital polarized and orbital
polarized state (VOP+OP). (d)Spin polarized metal (SP). (e) Valley polarized
semimetal (VP). (f) Inter-valley coherent state (IVC).

By performing Hartree–Fock calculations on the interacting 8-orbital model in
Eq. (7), we analyzed the solutions at various integer fillings. Our results
are summarized in Table 3 and in Figure 5. In some cases, we were able to find
a well-defined ground state that spontaneously breaks the symmetries of the
model; while in other cases, we identified several candidate states nearly
degenerate in energy. The resulting symmetry-broken states and the
characteristic energy differences are outlined in Figure 5. At certain
fillings, in particular $\nu=-3,-2,+2$, we found several candidate states that
have nearly identical total energies, so the competitions of these symmetry-
broken states are possible. we found that both insulating and (semi)metallic
ground states are stabilized, depending on the filling fraction of the central
bands. QAH solutions are also found. Below, we analyse in detail our results
for all integer fillings from $\nu=-3$ to $\nu=+3$.

Numerical results at filling $\bm{\nu=0}$. At the charge neutral point, our
calculations found insulating states and semimetallic states. We generated
several initial guesses by the symmetry argument as discussed in Sec. I, and
put them into the HF self-consistent cycles. Six different converged ordered
states are obtained, whose band structures are shown in Fig. 6 – the orbital
polarized (OP) state, the valley-orbital polarized (VOP) state, a polarized
state in which the up- [down-] spin sector is orbital [valley-orbital]
polarized (VOP+OP state), the spin-polarized (SP) semimetallic state, the
valley polarized (VP) semimetallic state, and the inter-valley coherent (IVC)
state.

The first three ordered states OP, VOP, and VOP+OP shown in Fig. 6 have
similar band structures, with nearly degenerate total energies (energy
difference within $0.13$ meV as shown in Fig. 5) among which the VOP has the
lowest total energy. However, as shown in Table 3, these three states have
different total Chern number $C=0,4,2$.

The semimetallic states in general have higher total energy compared with
gapped states. One would expect opening a gap when the system is fully spin-
polarized or valley-polarized due to Hund’s rule, but this is not the case for
our VP and SP states at charge neutrality. While it is true that the exchange
interaction is repulsive between two valleys and two spins species, its
magnitude is however too small to open up the gap at the $\Gamma$ point. We
note that the semimetallic VP and IVC symmetry-broken solutions preserve the
$C_{2z}\mathcal{T}$ symmetry, resulting in Dirac points at $K$ points. SP
state has the similar structure. It breaks the $C_{2z}\mathcal{T}$ symmetry
due to spin polarization. However, if only choosing one spin flavour, and
making $\mathcal{T}^{2}=1$, the spinless band has the $C_{2z}\mathcal{T}$
symmetry.

Therefore, our results suggest that orbital orders like
$\tau_{3}\otimes\Sigma_{3}$ (VOP) and $\Sigma_{3}$ (OP) are the reason for
generating insulating states. The mechanism is explained as follows. The
interaction-driven symmetry breaking in orbital-polarized channels can be
understood through a reduced Hamiltonian composed of $p_{\pm}$ orbitals only,
which is obtained by truncating the 8-orbital single-particle Hamiltonian. In
the non-interacting limit, such an effective Hamiltonian has the general form,

$\displaystyle
H_{\text{eff}}(\bm{k})=\sum_{\alpha,\beta=0}^{3}h^{\alpha\beta}(\bm{k})$ (14)

where

$h^{\alpha\beta}(\bm{k})=\sum_{s=\uparrow,\downarrow}t^{\alpha\beta}(\bm{k})\psi_{s}^{\dagger}(\bm{k})~{}\Gamma^{\alpha\beta}~{}\psi_{s}(\bm{k}),$
(15)

with $\Gamma^{\alpha\beta}=\tau_{\alpha}\otimes\Sigma_{\beta}$ are the
generators of $SU(4)$, $\psi_{s}=(c_{1,+,s}\quad c_{2,+,s}\quad c_{1,-,s}\quad
c_{2,-,s})^{\intercal}$, and $t^{\alpha\beta}(\bm{k})$ encodes the respective
hopping amplitudes. We note that, here, the spin degrees of freedom are
trivially summed over because the Hamiltonian is identity in spin space, as
the spin-orbit coupling is negligible in graphene. For the model we study, the
only non-vanishing terms in $H_{\text{eff}}$ correspond to
$(\alpha,\beta)=(0,0),(0,1),(0,2),(3,1),(3,2)$, and $(3,0)$. While the first
five terms obtain contributions from nearest and next-nearest neighbor
hoppings, $t^{30}$ lacks contributions from nearest neighbor hoppings.
Therefore, $t^{30}$ is exponentially suppressed compared to the remaining
hopping amplitudes in $H_{\text{eff}}$. We restrict the following discussion
to nearest-neighbor hoppings only, such that

$\displaystyle H_{\text{eff}}\to
H_{\text{eff},\text{NN}}=h^{00}+h^{01}+h^{02}+h^{31}+h^{32},$ (16)

where we have suppressed the $\bm{k}$-dependence for notational simplicity. We
will comment on the potential impact of further neighbor hoppings at the end
of the present analysis. We note that the four bands produced by diagonalizing
$H_{\text{eff},\text{NN}}$ are twofold spin-degenerate.

The term $t^{00}(\bm{k})$, while being of the same order as the other four
hopping amplitudes, does not directly control the gapping of the Dirac points
in the mini-Brillouin zone. This is because it effectively serves as a local
chemical potential inside each unit-cell, and commutes with all patterns of
internal symmetry breaking. Moreover, in the vicinity of the $K$ points of the
mini-Brillouin zone, $t^{00}(\bm{k})$ is only weakly $\bm{k}$-dependent.

If the order parameter of a symmetry-broken state globally anti-commutes with
$H_{\text{eff}}$, then it follows from the properties of the Clifford algebra
that the ordered state must be gapped.

For our effective model in Eq. (16), we found the following (anti-)commutation
relations:

$\displaystyle\left\\{(H_{\text{eff},\text{NN}}-h^{00}),\Gamma^{03}\right\\}=0,$
$\displaystyle\left\\{(H_{\text{eff},\text{NN}}-h^{00}),\Gamma^{33}\right\\}=0,$
$\displaystyle\left[(H_{\text{eff},\text{NN}}-h^{00}),\Gamma^{30}\right]=0.$
(17)

The generators $\Gamma^{03}$ and $\Gamma^{33}$ correspond to the OP and VOP
order parameters, respectively, and the anticommutation agrees with our
finding of the gap opening at charge neutrality. The semimetallic states VP
(generated by $\Gamma^{30}$) and SP (generated by $s_{3}$ Pauli matrix) have
order parameters that commute with the Hamiltonian, thus remaining gapless.

The size of the gap opened by the OP and VOP order parameters is controlled by
the strength of $H_{C}$, $U_{11}\approx U_{22}$ which is much larger than
$\text{max}\left[t^{30}(\bm{k})\right]$. Therefore, switching-on of the next-
nearest neighbor hoppings acts as weak perturbations, and does not
qualitatively change the above conclusions.

Numerical results at filling $\bm{\nu=\pm 2}$. For the filling $\nu=+2$, the
VP insulator, SP insulator, and IVC+VP metallic solutions are found (see Fig.
8a,b,c), with small energy differences (around 0.5 meV) among them. The two
insulating states have small band gaps around 3 meV. Among the three orders,
the VP state turns out to have the lowest total energy and realizes an
anomalous quantum Hall state, with the Chern number $|C|=2$ (by contrast, the
SP insulating state is topologically trivial with $C=0$).

At the filling $\nu=-2$, VP, SP, and OP ordered metallic states are obtained.
Among the three states, the VP state has the lowest total energy, however all
three states have very small energy differences within about 0.01 meV. As can
be seen from the band structures displayed in Appendix G (Fig. 8d,e,f), all
three states are semimetallic – they are not fully gapped between the hole
band centered at the $\Gamma$ point and the electron band minimum centered on
the $\Gamma-M$ line.

This phenomenon shows that higher-lying (non-flat) bands participate in the
formation of the order, which is the consequence of small band gaps around the
central bands of our model, and the fragile topology of the system. Compared
with the BM models, the 8-orbital model is away from the chiral limit and has
smaller gaps.

Although at hole doping $\nu=-2$ all the states we obtained are metallic, but
it is still possible to manually separate the higher lying conduction bands
and calculate the Chern number of the valence bands as a reference. Defined in
this fashion, the VP and SP states have $|C|=2$ and $C=0$ respectively, same
as the Chern numbers at $\nu=+2$ filling. We thus conclude that the VP
semimetallic state is a ‘failed’ Chern insulator, with non-trivial (non-
quantized) anomalous Hall response.

Unlike the approximate BM model which is particle-hole symmetric, the more
realistic 8-orbital model from Carr et al [1] model is not. Thus generically,
one does not expect to find the same solution at filling $\pm|\nu|$. This is
manifestly the case in our calculations at $|\nu|=2$, where we found the
candidate ground states to all be metallic on the hole-doped $\nu=-2$ side, in
contrast to the aforementioned opening of the gap at $\nu=2$.

Numerical results at filling $\bm{\nu=\pm 1}$. At fillings $\nu=\pm 1$, we
obtained two states on both the electron and hole-doped sides as shown in
Appendix G (Fig. 9): the intervalley-coherent spin-polarized (IVCSP) state and
a valley- and spin-polarized (VPSP) state. However, only one insulating state
is found – the VPSP state at $\nu=-1$ with a direct gap around 5.6 meV. On the
other hand, the IVCSP state (whose energy is 2 meV higher) has a biquadratic
band touching at the $\Gamma$ point where the gap closes. By contrast, the
$\nu=+1$ states are metallic states. Finally, we noted that all the
aforementioned states break the $C_{2z}\mathcal{T}$ symmetry. The insulating
VPSP state at $\nu=-1$ is a QAH state with Chern number $|C|=3$.

Note that the particle-hole asymmetry characteristic of the non-interacting
8-orbital model remains clearly evident in the symmetry-broken states. This
underlies our finding that the resulting Hartree–Fock band structures in the
presence of the interactions appear very asymmetric between $\nu=-1$ (Fig.
9a,b) and $\nu=+1$ (Fig. 9c,d), similar to what we found for $|\nu|=2$ above.

Figure 7: Comparison of insulating band structures with and without exchange
interaction. (a),(b) Valley polarized and spin polarized state at $\nu=-1$.
The exchange interaction enlarges the band gap. (c),(d) the valley polarized
state at $\nu=+2$. The exchange interaction opens a gap around fermi level.

Numerical results at filling $\bm{\nu=\pm 3}$. Fewer states are found at
filling $\nu=\pm 3$. As shown in Fig. 10, only a VPSP ordered state is found
at filling $\nu=+3$. The state is not fully gapped, but if we separate the
bands and calculate the total Chern number of the valence bands, we obtained
$|C|=1$, illustrating that this is an example of a ‘failed’ Chern insulator
similar to the VP state we found at $\nu=-2$.

By contrast, at filling $\nu=-3$, two metallic states IVCSP and VPSP are
obtained. The two states show the competing feature as their total energy
difference is within 0.01 meV.

###  III.3 The significance of Fock terms and exchange interactions

The previous work [69] studied the magic-angle TBG system without Fock terms
and exchange interactions. The resulting state found in the self-consistent
Hartree approximation is always a metallic state for all fillings, which
cannot explain most of the experimental results.

Fock terms are important for two reasons. On the one hand, Fock terms tend to
create symmetry-broken states that are comparable to experimental results. On
the other hand, keeping Fock terms makes the HF state a variational Slater
determinant and obey Wick’s theorem, which means systematic improvements
beyond HF can be made in the future based on HF states found in this work.

We emphasize that in the multi-orbital setting such as realized in TBG, the
“exchange interaction” is not synonymous with the “Fock term” that is
occasionally used, somewhat confusingly, to describe the exchange. Indeed, the
distinction between the direct Coulomb interaction $H_{C}$ in Eq. (8) and the
exchange interaction $H_{X}$ in Eq. (11) is made at the level of the
Hamiltonian, regardless of the approximate method such as Hartree or
Hartree–Fock used to optimize the ground state. Below, we analyze the effect
of the exchange interactions $H_{X}$ on the physics of TBG.

Given the much larger value of direct density-density interactions $H_{C}$ (as
large as 35 meV, see Fig. 3) compared to the strength of exchange interactions
$H_{X}$ ($\lesssim 3$ meV, Fig. 4), some of the previous studies [69] have
neglected the exchange interactions for this reason, performing calculations
with only density-density interactions. It is reasonable to ask whether
exchange interactions play a role in our HF calculations, if any, in
stabilizing the various symmetry-broken states we had identified.

To answer this question, we analyzed the results of our modeling with only the
direct density-density interactions in Eq. (8) vs. the full interacting
Hamiltonian including the exchange terms $H_{X}$ as in Eq. (11). We found that
the exchange interaction is essential for obtaining the correlated insulating
states we reported, except perhaps for the cases at charge neutrality. For
example, comparing the two approaches at filling $\nu=+2$ in Fig. 7(c)(d), we
found that excluding the exchange interaction would result in the gap closing
at the $\Gamma$ point. The same behavior is found at filling $\nu=-1$ as shown
in Fig. 7(a) and (b), where the band gap at $\Gamma$ point is much larger if
the exchange interactions are included.

These numerical observations can be simply explained by the generalized Hund’s
interaction involving the valley and spin degrees of freedom – the
corresponding indices are treated on an equal footing in the exchange
Hamiltonian in Eq. (11). The Hund’s effect originates from the terms that have
equal indices $\tau=\tau^{\prime}$ and $s=s^{\prime}$. In other words, the
fully spin-polarized or valley-polarized bands are preferentially filled to
lower the total energy. Of course the $\mathbf{k}$ dependence of the kinetic
energy makes the picture more complicated, but the above argument still
applies qualitatively. By the same token, the Hund’s rule would suggest that
intervalley-coherent states are less energetically favorable, which is indeed
reflected in our numerical findings summarized in Fig. 5.

##  IV Comparison with experiments

The particle-hole asymmetric nature of the correlated states we found in our
calculations agrees broadly with the absence of particle-hole symmetry in the
experiments. In this section, we provide a detailed comparison with
experimental data, highlighting the agreement with our computational results.

We start from the comparison at the charge-neutral point, where scanning
tunneling microscopy (STM) experiments reported strong local correlations
[112, 109, 110, 118], while transport measurements identified semimetallic
features [3, 4, 105, 119]. Our results suggest several candidate insulating
and semimetallic states at the charge-neutral point. Crucially, our
semimetallic states exhibit symmetry-breaking, distinguishing them from the
strongly correlated symmetry-preserving states obtained by the heavy fermion
model [101, 100, 102].

The experiments, although they depend on twist angles and sample alignment
with hBN substrates, have reported quantized Hall conductance at certain
fillings. For instance, in transport measurements [106, 107], QAH states at
$\nu=+3$ were observed. Notably, our calculation yielded a ‘failed’ Chern
insulator with $|C|=1$ at this filling. Conversely, on the hole side at
$\nu=-3$, our findings indicate a metallic state. Consistent with our results,
there is no observed resistance peak at $\nu=-3$ in prior studies [107, 114].
This discrepancy may stem from the particle-hole asymmetric nature,
accentuated by lattice relaxation.

Our results at $\nu=-1$ and $\nu=+2$ also display particle-hole asymmetric
features, but unlike $|\nu|=3$, the QAH states are actually insulating in our
calculations and possess well-defined Chern numbers. Experimental observations
at $\nu=+1$ show a correlated Chern insulator and a Chern number transition
from $C=1$ to $C=3$ under an external magnetic field [128], whereas our
findings only reveal a $C=3$ QAH state at hole doping $\nu=-1$.

At $\nu=\pm 2$, experiments typically indicate correlated insulating features,
although outcomes may vary between devices [114]. Our calculations for $\nu=2$
demonstrate a anomalous Hall insulator with a direct gap and $|C|=2$. On the
hole side, $\nu=-2$, we obtained a metallic state without Chern number.

To conclude, a comparison between our results and experimental findings
suggests that the mechanisms for spontaneous symmetry breaking and gap
formation exhibit particle-hole asymmetry and depend on the experimental
details such as the twist angle and sample alignment with hBN substrate. Our
HF calculations underscore the importance of incorporating particle-hole
asymmetric models into theoretical considerations.

##  V Conclusion

In this work, motivated by the challenge of accurately capturing the
interactions in the magic-angle TBG without sacrificing the ab initio
perspective, we have studied a multi-orbital model that circumvents the
fragile topological obstruction by enlarging the active orbital space. We
first employed the localized Wannier orbitals to numerically evaluate the
matrix elements of the electron repulsion, enabling us to construct a suitably
extended Hubbard model including both the direct and exchange interactions.
Subsequently, we incorporated the complete Hartree and Fock terms to study
this model.

By performing Hartree–Fock calculations, we have obtained a variety of states
as a function of integer fillings from $-3$ to $+3$ of the nearly-flat bands.
These symmetry-broken ordered states originate from the valley, spin, and
orbital degrees of freedom. The orbital polarization, which spontaneously
breaks the discrete $C_{2z}\mathcal{T}$ symmetry, appears promptly at several
fillings and is particularly important, as it is missing from the alternative
treatments such as those based on the Bistritzer–MacDonald model. Our results
are not symmetric with respect to electron and hole doping, consistent with
the experimental observations on TBG near the magic angle (see previous
Section IV for detailed comparison with experiments). This asymmetry is
natural for the model derived from first principles which lacks particle-hole
symmetry, in contrast to approximate models such as the BM model that have
been subject of many previous studies.

We found symmetry-broken insulating states at $\nu=0,-1,+2$. In many cases,
several ground state candidates have very close total energies, less than 0.05
meV apart, suggesting that the precise nature of the ground state depends
delicately on the microscopic model parameters, which in turn depend
sensitively on the twist angle, encapsulation, defects, and strain effects.
Properly including exchange interactions is crucial to obtain insulating
states especially at $\nu=-1,+2$.

We found the insulating states at $\nu=-1$ and $\nu=2$ to be anomalous Hall
insulators characterized by a non-zero Chern number of the occupied bands.
Several semimetallic states with clearly separable but overlapping in energy
conduction and valence bands can also be described as ‘failed’ Chern
insulators, with a well defined Chern number of the valence bands. These
metallic states are predicted to have (non-quantized) anomalous Hall effect
due to non-zero integrated Berry curvature of the occupied bands.

Hartree-Fock is of course a limited theory for interacting systems. In
comparison to Hartree-only approximations, the Fock term tends to create
broken-symmetry mean field states, of which we find a number that roughly
correspond to experimental observations; however, we view it as unlikely that
the Hartree-Fock ground states are sufficiently accurate to determine the
ground state with high confidence. Recently, [102], there have been
DMFT+Hartree studies of the same model, along with the aforementioned Hartree-
only [69] calculations. Those calculations drop both the exchange $H_{X}$ and
the Fock term from consideration, which we show are both important to
stabilizing symmetry-broken states. In future work, it would be interesting to
see if correlated calculations correct towards the symmetry-broken states.

###### Acknowledgements.

We would like to express our gratitude to Fang Xie for engaging in helpful
discussions. This research was supported by the U.S. Department of Energy
Computational Materials Sciences (CMS) program under Award Number DE-
SC0020177. A.H.N. is grateful for the hospitality of the Aspen Center for
Physics, supported by the National Science Foundation grant PHY-2210452, where
a portion of this work was performed.

##  Appendix A Symmetries in 8-orbital TBG model

We verified the symmetries of the non-interacting Hamiltonian. While we
initially described the non-interacting Hamiltonian using a periodic gauge for
the sake of convenience in HF calculations, it is more straightforward to
describe the symmetry transformation using a physical gauge. These two
Hamiltonians are connected through a gauge transformation denoted as
$U_{\text{gauge}}(\bm{k})$. Whenever we require a specific gauge choice, we
can perform the transformation as follows:

$U^{\dagger}_{\text{gauge}}(\bm{k})H^{\text{Periodic}}_{K}(\bm{k})U_{\text{gauge}}(\bm{k})=H^{\text{Physical}}_{K}(\bm{k}).$
(18)

We discuss the symmetry operations under physical gauge. The non-interacting
Hamiltonian written in momentum space is

$H_{K}=\sum_{ab\tau s}t^{\tau}_{ab}(\bm{k})c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau s}(\bm{k}).$ (19)

We simply denote the matrix of $[t^{\tau}_{ab}(\bm{k})]$ as
$h^{\tau}(\bm{k})$. When there is no valley dependence, index $\tau$ will be
omitted. In the following verification, We use the notation $D(g)$ to label
the representation matrix of a symmetry operation $g$. We consider a valley-
independent symmetry operation acting on a creation operator:

$g^{-1}c^{\dagger}_{a\tau s}(\bm{k})g=\sum_{b}D^{*}_{ba}(g)c^{\dagger}_{b\tau
s}(g\bm{k}).$ (20)

This gives the definition of matrix $D(g)$. Using the symmetry condition
$g^{-1}H_{K}g=H_{K}$, we could derive how the matrix $h^{\tau}(\bm{k})$
transforms under a symmetry operation.

The above definition assumes the operation does not mix two valleys. We
consider an example that does mix two valleys—the time-reversal symmetry. The
time-reversal symmetry is anti-unitary and relates two valleys

$\mathcal{T}^{-1}c_{a\tau s}(\bm{k})\mathcal{T}=c_{a-\tau s}(-\bm{k}).$ (21)

Note here we choose the spinless anti-unitary convention $\mathcal{T}^{2}=1$,
but if an order is related to any spin polarization, we choose back to
original definition $\mathcal{T}^{2}=-1$. It is easy to verify that

$h^{\tau*}(-\bm{k})=h^{-\tau}(\bm{k}).$ (22)

The relationship helps to construct the whole non-interacting Hamiltonian from
a single-valley Hamiltonian.

There are three important Symmetries $C_{2x}$, $C_{2z}\mathcal{T}$, and
$C_{3z}$ in the magic-angle TBG system. According to the lattice structure and
the symmetries of Wannier orbitals, we can write down the symmetry
representation matrices easily. The Wannier basis in a single valley is
denoted as
$\left(\ket{1,p_{+}},\ket{2,p_{-}},\ket{3,s},\ket{4,p_{z}},\ket{5,p_{z}},\ket{6,s},\ket{7,s},\ket{8,s}\right)$
where the numbers label the orbital indices as in Table 2. We use
$\mathcal{K}$ to represent the complex conjugation, in case there is a anti-
unitary operator. The representation matrices of these three operators can be
derived as follows:

$D(C_{2z}\mathcal{T})=\begin{pmatrix}0&-1&0&0&0&0&0&0\\\ -1&0&0&0&0&0&0&0\\\
0&0&1&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&0&1&0&0\\\
0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\\\ \end{pmatrix},$ (23)

$D(C_{2x})=\begin{pmatrix}0&-1&0&0&0&0&0&0\\\ -1&0&0&0&0&0&0&0\\\
0&0&1&0&0&0&0&0\\\ 0&0&0&-1&0&0&0&0\\\ 0&0&0&0&-1&0&0&0\\\ 0&0&0&0&0&1&0&0\\\
0&0&0&0&0&0&0&1\\\ 0&0&0&0&0&0&1&0\\\ \end{pmatrix},$ (24)

$D(C_{3z})=\begin{pmatrix}e^{-i2\pi/3}&0&0&0&0&0&0&0\\\
0&e^{i2\pi/3}&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\
0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&0&1\\\ 0&0&0&0&0&1&0&0\\\
\end{pmatrix}.$ (25)

Applying these three matrices to the non-interacting Hamiltonian $h(\bm{k})$,
the Hamiltonian must satisfy the following relationships due to symmetry
constraint:

$D^{-1}(C_{2z}\mathcal{T})h^{*}(\bm{k})D(C_{2z}\mathcal{T})=h(\bm{k}),$ (26)
$D^{-1}(C_{2x})h(\bm{k})D(C_{2x})=h(C_{2x}\bm{k}),$ (27)

and

$D^{-1}(C_{3z})h(\bm{k})D(C_{3z})=h(C_{3z}\bm{k}).$ (28)

We verified that the non-interacting Hamiltonian has the correct symmetry.

##  Appendix B Interacting Hamiltonian

The four-center electron repulsion integrals are evaluated numerically using
Wannier functions:

$V^{\tau_{1}\tau_{2}\tau_{3}\tau_{4}}_{a_{1}a_{2}a_{3}a_{4}}(\bm{R}_{1},\bm{R}_{2},\bm{R}_{3},\bm{R}_{4})=\int\mathrm{d}\bm{r}^{2}\mathrm{d}{\bm{r}^{\prime}}^{2}V_{SG}(|\bm{r}-\bm{r}^{\prime}|)\sum_{XX^{\prime}}\mathcal{W}^{X}_{\bm{R}_{1}a_{1}\tau_{1}}(\bm{r})\mathcal{W}^{X}_{\bm{R}_{4}a_{4}\tau_{4}}(\bm{r})\mathcal{W}^{X^{\prime}}_{\bm{R}_{2}a_{2}\tau_{2}}(\bm{r}^{\prime})\mathcal{W}^{X^{\prime}}_{\bm{R}_{3}a_{3}\tau_{3}}(\bm{r}^{\prime}),$
(29)

where the Wannier functions of two valleys are related by complex conjugation:

$\mathcal{W}^{X}_{\bm{R}a}(\bm{r}):=\mathcal{W}^{X}_{\bm{R}a\tau=+}(\bm{r})=\mathcal{W}^{X}_{\bm{R}a\tau=-}(\bm{r})^{*}.$
(30)

The Wannier function is proportional to an overall valley-dependent
modulation:

$\mathcal{W}^{X}_{\bm{R}a\tau}(\bm{r})\propto\sum^{3}_{j=1}e^{i\tau\bm{K}^{j}_{\tau}\cdot(\bm{r}-\bm{R}_{X})}.$
(31)

$\bm{K}^{j}_{\tau}$ are pristine graphene valley vectors and $\bm{R}_{X}$ are
sublattice vectors in real-space. This modulation structure suppress those
integrals whose $\tau_{1}\neq\tau_{4}$ and $\tau_{2}\neq\tau_{3}$, as
explained in the main text. We only considered two-center integrals assuming
other integrals are small. We can define the direct-channel interaction
integrals
$U_{a_{1}a_{2}}(\bm{R}_{1}-\bm{R}_{2}):=V^{\tau\tau^{\prime}\tau^{\prime}\tau}_{a_{1}a_{2}a_{2}a_{1}}(\bm{R}_{1},\bm{R}_{2},\bm{R}_{2},\bm{R}_{1})$,
and exchange-channel interaction integrals
$X^{\tau\tau^{\prime}}_{a_{1}a_{2}}(\bm{R}_{1}-\bm{R}_{2}):=V^{\tau\tau^{\prime}\tau^{\prime}\tau}_{a_{1}a_{2}a_{1}a_{2}}(\bm{R}_{1},\bm{R}_{2},\bm{R}_{1},\bm{R}_{2})$.

The exchange interaction is important to determine the phases, although it is
small than the density-density interaction $U_{ab}(\bm{R}_{i}-\bm{R}_{j})$. We
denote the exchange integral of two orbitals $a$ and $b$ located at unit cells
$\bm{R}_{i}$ and $\bm{R}_{j}$ as
$X^{\tau\tau^{\prime}}_{ab}(\bm{R}_{i}-\bm{R}_{j})$. Because it has dependence
on valleys, we can write it in two parts: the intra-valley exchange integral
$J_{ab}(\bm{R}_{i}-\bm{R}_{j})$ when $\tau=\tau^{\prime}$ and the inter-valley
exchange integral $P_{ab}(\bm{R}_{i}-\bm{R}_{j})$ when
$\tau\neq\tau^{\prime}$. The matrix form of the integral is written as

$\left[X^{\tau\tau^{\prime}}_{ab}(\bm{R}_{i}-\bm{R}_{j})\right]=\begin{pmatrix}J_{ab}(\bm{R}_{i}-\bm{R}_{j})&P_{ab}(\bm{R}_{i}-\bm{R}_{j})\\\
P_{ab}(\bm{R}_{i}-\bm{R}_{j})^{*}&J_{ab}(\bm{R}_{i}-\bm{R}_{j})\end{pmatrix}.$
(32)

Using the condition that the Wannier functions of two valleys are connected
with a complex conjugation, these two integrals are evaluated numerically
using the expressions as following:

$J_{ab}(\bm{R}_{i}-\bm{R}_{j})=\int\mathrm{d}\bm{r}^{2}\mathrm{d}{\bm{r}^{\prime}}^{2}V_{SG}(|\bm{r}-\bm{r}^{\prime}|)\sum_{XX^{\prime}}\mathcal{W}^{X}_{\bm{R}_{i}a}(\bm{r})^{*}\mathcal{W}^{X}_{\bm{R}_{j}b}(\bm{r})\mathcal{W}^{X^{\prime}}_{\bm{R}_{j}b}(\bm{r}^{\prime})^{*}\mathcal{W}^{X^{\prime}}_{\bm{R}_{i}a}(\bm{r}^{\prime}),$
(33)

$P_{ab}(\bm{R}_{i}-\bm{R}_{j})=\int\mathrm{d}\bm{r}^{2}\mathrm{d}{\bm{r}^{\prime}}^{2}V_{SG}(|\bm{r}-\bm{r}^{\prime}|)\sum_{XX^{\prime}}\mathcal{W}^{X}_{\bm{R_{i}}a}(\bm{r})^{*}\mathcal{W}^{X}_{\bm{R}_{j}b}(\bm{r})\mathcal{W}^{X^{\prime}}_{\bm{R}_{i}a}(\bm{r}^{\prime})^{*}\mathcal{W}^{X^{\prime}}_{\bm{R}_{j}b}(\bm{r}^{\prime}).$
(34)

We set the elements $X^{\tau\tau^{\prime}}_{aa}(\bm{0})=0$ to avoid double
counting of interactions that have already been calculated in density-density
interaction.

Combining all the interaction terms with the tight-binding model, the full
interacting Hamiltonian is

$H_{C}+H_{X}=H_{C}+\frac{1}{2}\sum_{a,b,\tau,\tau^{\prime},s,s^{\prime}}\sum_{ij}X^{\tau\tau^{\prime}}_{ab}(\bm{\Delta}\bm{R}_{j}):c^{\dagger}_{a\tau
s}(\bm{R}_{i})c_{b\tau
s}(\bm{R}_{i}+\bm{\Delta}\bm{R}_{j})c^{\dagger}_{b\tau^{\prime}s^{\prime}}(\bm{R}_{i}+\bm{\Delta}\bm{R}_{j})c_{a\tau^{\prime}s^{\prime}}(\bm{R}_{i}):.$
(35)

##  Appendix C Hartree-Fock mean-field theory

In this section, we review the HF method and basic notations.

The Hartree-Fock order parameter which can be also called the one-body reduced
density matrix (1rdm) is defined in the orbital Bloch basis:

$P_{a\tau s,b\tau^{\prime}s^{\prime}}(\bm{k})=\left<c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})-\frac{1}{2}\delta_{ab}\delta_{\tau\tau^{\prime}}\delta_{ss^{\prime}}\right>.$
(36)

Here the subtraction
$\frac{1}{2}\delta_{ab}\delta_{\tau\tau^{\prime}}\delta_{ss^{\prime}}$ is made
to counter the double counting of the interaction. The Hartree-Fock mean-field
Hamiltonian $\mathcal{H}^{HF}(\bm{k})$ has dependence on 1rdm $P_{a\tau
s,b\tau^{\prime}s^{\prime}}(\bm{k})$ and is derived in the Appendix E. Solving
the self-consistent eigenvalue problem of the HF Hamiltonian

$\mathcal{H}^{HF}(\bm{k})\ket{u_{n}(\bm{k})}=\varepsilon_{n}(\bm{k})\ket{u_{n}(\bm{k})},$
(37)

we obtained the HF band dispersion $\varepsilon_{n}(\bm{k})$ and corresponding
band wave functions $u_{n}(\bm{k})$. The direct inversion of the iterative
subspace (DIIS) [133, 134] and energy-DIIS (EDIIS) [135] algorithms are
implemented to accelerate the convergence.

There are three important symmetries in magic-angle TBG system, in-plane
$120\degree$ rotational symmetry $C_{3z}$, out-of-plane $180\degree$
rotational symmetry $C_{2x}$, and $C_{2z}T$ symmetry. $C_{2z}T$ symmetry
combines in-plane $180\degree$ rotation with time-reversal symmetry. The
representation matrix $D(g)$ for symmetry $g$ can be found in Appendix A.
Order parameters with breaking symmetries can lead to topological non-trivial
states. For example, breaking the $C_{2z}T$ symmetry is necessary for getting
non-zero Chern number state; Therefore, We define the $C_{2z}T$ symmetry-
breaking order parameter to capture the $C_{2z}T$ symmetry-breaking strength.
The commutator-like order parameter denotes:

$\mathcal{D}(\bm{k})=\left|P(\bm{k})D(C_{2z}\mathcal{T})-D(C_{2z}\mathcal{T})P(\bm{k})^{*}\right|$
(38)

where $D(C_{2z}T)$ is the representation matrix of $C_{2z}T$ symmetry.

The insulating states with non-zero Chern number are identified as QAH states.
The total Chern number $C$ is calculated from the determinant of the Wilson
loop matrix of the occupied bands (see Appendix F). Many trial states are
randomly generated to initialize the self-consistent calculation, but only
several relevant converged results come out.

##  Appendix D Derivation of the symmetry breaking strength

The symmetry breaking strength of a HF converged state can be defined using a
commutation-like expression:

$\left|\left<g^{-1}c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})g\right>-\left<c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})\right>\right|,$ (39)

where $g$ is the symmetry operation. Using the representation matrix $D(g)$
defined in Appendix A, we obtain

$\begin{split}&\left|\left<g^{-1}c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})g\right>-\left<c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})\right>\right|\\\
=&\left|D^{*}_{ca}(g)D_{db}(g)\left<c^{\dagger}_{c\tau s}(g\bm{k})c_{d\tau
s}(g\bm{k})\right>-\left<c^{\dagger}_{a\tau s}(\bm{k})c_{b\tau
s}(\bm{k})\right>\right|\\\
=&\left|D^{\dagger}(g)P(g\bm{k})D(g)-P(\bm{k})\right|\end{split}$ (40)

for unitary transformation. Multiplying $D(g)$ on the left side, the symmetry
breaking strength $\mathcal{D}(\bm{k})$ can be expressed as

$\mathcal{D}(\bm{k})=\left|P(g\bm{k})D(g)-D(g)P(\bm{k})\right|.$ (41)

For anti-unitary transformation, the same procedure gives

$\mathcal{D}(\bm{k})=\left|P(g\bm{k})D(\mathcal{K}g)-D(\mathcal{K}g)P(\bm{k})^{*}\right|.$
(42)

We mainly checked the breaking of $C_{2z}\mathcal{T}$ symmetry, because
breaking $C_{2z}\mathcal{T}$ symmetry is the necessary condition for a finite
Chern number.

##  Appendix E Hartree-Fock mean-field Hamiltonian

We used Hartree-Fock mean-field method to study the 8-band extended Hubbard
model plus exchange interactions $H=H_{K}+H_{C}+H_{X}$. In this section, we
derived the HF self-consistent equations.

###  E.1 Hartree-Fock decomposition of density-density interacting Hamiltonian

In the first part, we derive the HF self-consistent equations of the density-
density interaction part. The density-density interacting Hamiltonian is
expressed as

$H_{C}=\sum_{a,b,\tau,\tau^{\prime},s,s^{\prime}}\left(\sum_{i}\frac{1}{2}U_{ab}(\bm{0}):n_{a\tau
s}(\bm{R}_{i})n_{b\tau^{\prime}s^{\prime}}(\bm{R}_{i}):+\sum_{ij}\frac{1}{2}U_{ab}(\bm{\Delta}\bm{R}_{j}):n_{a\tau
s}(\bm{R}_{i})n_{b\tau^{\prime}s^{\prime}}(\bm{R}_{i}+\bm{\Delta}\bm{R}_{j}):\right).$
(43)

Assuming there is no translational symmetry breaking, we do the Fourier
transform and obtain the Hamiltonian in momentum space,

$H_{C}=\sum_{a,b,\tau,\tau^{\prime},s,s^{\prime}}\sum_{\bm{q}}\frac{1}{2N_{k}}\bar{U}_{ab}(\bm{q})\delta\rho_{a\tau
s}(\bm{q})\delta\rho_{b\tau^{\prime}s^{\prime}}(-\bm{q}),$ (44)

in which we define the charge operator $\delta\rho_{a\tau s}(\bm{q})$ as

$\delta\rho_{a\tau s}(\bm{q})=\sum_{\bm{k}}c^{\dagger}_{a\tau
s}(\bm{k})c_{a\tau s}(\bm{q}+\bm{k})-\frac{1}{2}\delta_{\bm{q},\bm{G}}.$ (45)

When doing the Fourier transform, we chose periodic gauge, so
$\delta\rho(\bm{q}+\bm{G})=\delta\rho(\bm{q})$. And $\bar{U}_{ab}(\bm{q})$ is
the Fourier transformation of real-space Coulomb repulsion, which means

$\bar{U}_{ab}(\bm{q})=\left(V_{ab}(\bm{0})+\sum_{j}U_{ab}(\bm{\Delta}\bm{R}_{j})e^{i\bm{\Delta}\bm{R}_{j}\cdot\bm{q}}\right).$
(46)

To perform the HF decomposition, we define the one-body reduced density matrix
as

$P_{a\tau s,b\tau^{\prime}s^{\prime}}(\bm{k})=\left<c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})-\frac{1}{2}\delta_{ab}\delta_{\tau\tau^{\prime}}\delta_{ss^{\prime}}\right>,$
(47)

in which the braket is evaluated using mean-field Slater determinant. The
single-body Hamiltonian from Hartree decomposition is obtained:

$\mathcal{H}^{(H)}=\frac{1}{N_{k}}\sum_{\bm{k}^{\prime}\bm{k}}\sum_{a,b,\tau,\tau^{\prime},s,s^{\prime}}\bar{U}_{ab}(\bm{0})P_{a\tau
s;a\tau
s}(\bm{k}^{\prime})\left(c^{\dagger}_{b\tau^{\prime}s^{\prime}}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})-\frac{1}{2}\right),$
(48)

and the Hamiltonian from Fock decomposition is

$\mathcal{H}^{(F)}=-\frac{1}{2N_{k}}\sum_{\bm{k}\bm{k}^{\prime}}\sum_{ab\tau\tau^{\prime}ss^{\prime}}\bar{U}_{ab}(\bm{k}-\bm{k}^{\prime})P_{b\tau^{\prime}s^{\prime};a\tau
s}(\bm{k}^{\prime})\left(c^{\dagger}_{a\tau
s}(\bm{k})c_{b\tau^{\prime}s^{\prime}}(\bm{k})-\frac{1}{2}\delta_{ab}\delta_{\tau\tau^{\prime}}\delta_{ss^{\prime}}\right)+h.c..$
(49)

###  E.2 Hartree-Fock theory for exchange interaction

Here we show the HF decomposition of the exchange interacting Hamiltonian
$H_{X}$ as derived in Eq. (35). We do the Fourier transform using periodic
gauge to obtain the momentum-space exchange interaction:

$H_{X}=\sum_{ab\tau\tau^{\prime}ss^{\prime}}\sum_{\bm{k}_{1}...\bm{k}_{4}}\sum_{j}\frac{1}{2N_{k}}X^{\tau\tau^{\prime}}_{ab}(\bm{\Delta}\bm{R}_{j})e^{i(\bm{k}_{3}-\bm{k}_{2})\cdot\bm{\Delta}\bm{R}_{j}}\delta(\bm{k}_{1}-\bm{k}_{2}+\bm{k}_{3}-\bm{k}_{4}):c^{\dagger}_{a\tau
s}(\bm{k}_{1})c_{b\tau
s}(\bm{k}_{2})c^{\dagger}_{b\tau^{\prime}s^{\prime}}(\bm{k}_{3})c_{a\tau^{\prime}s^{\prime}}(\bm{k}_{4}):.$
(50)

Using the same prescriptions as shown in density-density interaction, the
Hartree and Fock terms are derived as

$\mathcal{H}^{(H)}_{X}=\frac{1}{2N_{k}}\sum_{\bm{k}\bm{k}^{\prime}}\sum_{ab\tau\tau^{\prime}ss^{\prime}}\sum_{j}X^{\tau\tau^{\prime}}_{ab}(\bm{\Delta}\bm{R}_{j})e^{i(\bm{k}^{\prime}-\bm{k})\cdot\bm{\Delta}\bm{R}_{j}}P_{a\tau
s,b\tau
s}(\bm{k})\left(c^{\dagger}_{b\tau^{\prime}s^{\prime}}(\bm{k}^{\prime})c_{a\tau^{\prime}s^{\prime}}(\bm{k}^{\prime})-\frac{1}{2}\delta_{ab}\right)+h.c.,$
(51)

and

$\mathcal{H}^{(F)}_{X}=-\frac{1}{2N_{k}}\sum_{\bm{k}\bm{k}^{\prime}}\sum_{ab\tau\tau^{\prime}ss^{\prime}}\sum_{j}X^{\tau\tau^{\prime}}_{ab}(\bm{\Delta}\bm{R}_{j})P_{a\tau
s,a\tau^{\prime}s^{\prime}}(\bm{k})\left(c^{\dagger}_{b\tau^{\prime}s^{\prime}}(\bm{k}^{\prime})c_{b\tau
s}(\bm{k}^{\prime})-\frac{1}{2}\delta_{\tau\tau^{\prime}}\delta_{ss^{\prime}}\right)+h.c..$
(52)

Combining with tight-binding Hamiltonian, the total Hartree-Fock mean-field
Hamiltonian can be written as

$\mathcal{H}^{HF}=H_{K}+\mathcal{H}^{(H)}+\mathcal{H}^{(F)}+\mathcal{H}^{(H)}_{X}+\mathcal{H}^{(F)}_{X}.$
(53)

The mean-field Hamiltonian $\mathcal{H}^{HF}$ is self-consistently determined
by its own eigenvectors. So the self-consistent mean-field states are obtained
by performing self-consistent cycles. The total energy is evaluated by

$E_{tot}=\left<H_{K}+\frac{1}{2}(\mathcal{H}^{(H)}+\mathcal{H}^{(F)}+\mathcal{H}^{(H)}_{X}+\mathcal{H}^{(F)}_{X})\right>.$
(54)

##  Appendix F Wilson loops and Chern numbers

This section explains how to calculate Wilson loop matrices. The total Chern
number is calculated by counting the windings.

We define the momentum points in the first Brillouin zone
$\bm{k}=\frac{1}{2\pi}(\tilde{k}_{1}\bm{G}_{1},\tilde{k}_{2}\bm{G_{2}})=(\frac{n_{x}}{N_{x}}\bm{G}_{1},\frac{n_{y}}{N_{y}}\bm{G_{2}})$,
where $N_{x}$ and $N_{y}$ are the sampling lengths on each direction. The pair
$(n_{x},n_{y})$ labels the discretized momentum, so
$(\tilde{k}_{1},\tilde{k}_{2})$ can be viewed as reduced k points. Supposing
$u_{an}(\mathbf{k})$ is the eigenvector of band $n$ in periodic gauge, we can
define the Berry connection as $\bm{A}=(A^{1},A^{2})$, in which
$A^{1}_{nm}=\sum_{a}iu^{*}_{an}(\tilde{k}_{1},\tilde{k}_{2})\frac{\partial}{\partial_{\tilde{k}_{1}}}u_{am}(\tilde{k}_{1},\tilde{k}_{2})$.

The non-Abelian Wilson loop matrix can be written as

$\begin{split}W_{nm}(\tilde{k}_{2})&=\mathcal{P}\exp{-i\int^{2\pi}_{0}d\tilde{k}_{1}A^{1}_{nm}(\tilde{k}_{1},\tilde{k}_{2})}\\\
&\approx\prod^{2\pi-\Delta\tilde{k}}_{\tilde{k}_{1}=0}\sum_{a}u^{*}_{an}(\tilde{k}_{1},\tilde{k}_{2})u_{am}(\tilde{k}_{1}-\Delta\tilde{k},\tilde{k}_{2}).\end{split}$
(55)

In practice, to make the numerical matrix productions stable, the singular
value decomposition (SVD) is applied. The total Chern number is obtained by
counting the windings of $\det W_{nm}$ of occupied bands.

##  Appendix G Details of the band structures

Figure 8: Band structures of self-consistent ground state candidates at
$\nu=\pm 2$. The color bar shows the valley polarization strength for each k
point. (a) Valley polarized state (VP). (B) Spin polarized state (SP). (C)
Inter-valley coherent and valley polarized state(IVC+VP). (d) Valley polarized
state (VP). (e) Spin polarized state (SP). (f) Orbital polarized state (OP).
Figure 9: Band structures of self-consistent ground state candidates at
$\nu=\pm 1$. The color bar shows the valley polarization strength for each k
point. (a),(c): Valley-polarized and spin-polarized state (VPSP). (b),(d)
Inter-valley coherent and spin polarized state (IVCSP).  Figure 10: Band
structures of self-consistent ground state candidates at $\nu=\pm 3$. The
color bar shows the valley polarization strength for each k point. (a) Valley-
polarized and spin-polarized state (VPSP) for $\nu=-3$. (b) Inter-valley
coherent and spin-polarized state (IVCSP) for $\nu=-3$. (c) VPSP for $\nu=+3$.

## References

  * Carr et al. [2019] S. Carr, S. Fang, H. C. Po, A. Vishwanath, and E. Kaxiras, Phys. Rev. Res. 1, 033072 (2019).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Cao et al. [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018a).
  * Cao et al. [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Moon and Koshino [2012] P. Moon and M. Koshino, Physical Review B 85, 195458 (2012).
  * Kang and Vafek [2018] J. Kang and O. Vafek, Phys. Rev. X 8, 031088 (2018).
  * Koshino et al. [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * Zou et al. [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Phys. Rev. B 98, 085435 (2018).
  * Po et al. [2018a] H. C. Po, H. Watanabe, and A. Vishwanath, Phys. Rev. Lett. 121, 126402 (2018a).
  * Efimkin and MacDonald [2018] D. K. Efimkin and A. H. MacDonald, Phys. Rev. B 98, 035404 (2018).
  * Isobe et al. [2018] H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X 8, 041041 (2018).
  * Liu et al. [2018] C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Physical review letters 121, 217001 (2018).
  * Xu and Balents [2018] C. Xu and L. Balents, Physical review letters 121, 087001 (2018).
  * Ochi et al. [2018] M. Ochi, M. Koshino, and K. Kuroki, Phys. Rev. B 98, 081102 (2018).
  * Xu et al. [2018] X. Y. Xu, K. T. Law, and P. A. Lee, Phys. Rev. B 98, 121406 (2018).
  * Po et al. [2018b] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Physical Review X 8, 031089 (2018b).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Proceedings of the National Academy of Sciences 115, 13174 (2018).
  * Guo et al. [2018] H. Guo, X. Zhu, S. Feng, and R. T. Scalettar, Phys. Rev. B 97, 235453 (2018).
  * Wu et al. [2018] F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121, 257001 (2018).
  * Yuan and Fu [2018] N. F. Yuan and L. Fu, Physical Review B 98, 045103 (2018).
  * Venderbos and Fernandes [2018] J. W. F. Venderbos and R. M. Fernandes, Phys. Rev. B 98, 245103 (2018).
  * Padhi et al. [2018] B. Padhi, C. Setty, and P. W. Phillips, Nano letters 18, 6175 (2018).
  * Wu et al. [2019] X.-C. Wu, C.-M. Jian, and C. Xu, Physical Review B 99 (2019), 10.1103/physrevb.99.161405.
  * Thomson et al. [2018] A. Thomson, S. Chatterjee, S. Sachdev, and M. S. Scheurer, Physical Review B 98 (2018), 10.1103/physrevb.98.075109.
  * Kennes et al. [2018] D. M. Kennes, J. Lischner, and C. Karrasch, Phys. Rev. B 98, 241407 (2018).
  * Dodaro et al. [2018] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X.-Q. Sun, and C. Wang, Physical Review B 98, 075154 (2018).
  * Song et al. [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Physical Review Letters 123, 036401 (2019).
  * Hejazi et al. [2019a] K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Phys. Rev. B 99, 035111 (2019a).
  * Tarnopolsky et al. [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Physical Review Letters 122, 106405 (2019).
  * Da Liao et al. [2019] Y. Da Liao, Z. Y. Meng, and X. Y. Xu, Phys. Rev. Lett. 123, 157601 (2019).
  * Liu et al. [2019a] J. Liu, J. Liu, and X. Dai, Physical Review B 99, 155415 (2019a).
  * Hejazi et al. [2019b] K. Hejazi, C. Liu, and L. Balents, Physical Review B 100 (2019b), 10.1103/physrevb.100.035115.
  * You and Vishwanath [2019] Y.-Z. You and A. Vishwanath, npj Quantum Materials 4, 16 (2019).
  * Lian et al. [2019] B. Lian, Z. Wang, and B. A. Bernevig, Phys. Rev. Lett. 122, 257002 (2019).
  * Classen et al. [2019] L. Classen, C. Honerkamp, and M. M. Scherer, Phys. Rev. B 99, 195120 (2019).
  * Zhang et al. [2019] Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Physical Review B 99, 075127 (2019).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Physical Review Letters 122, 246401 (2019).
  * Seo et al. [2019] K. Seo, V. N. Kotov, and B. Uchoa, Phys. Rev. Lett. 122, 246402 (2019).
  * Hu et al. [2019] X. Hu, T. Hyart, D. I. Pikulin, and E. Rossi, Phys. Rev. Lett. 123, 237002 (2019).
  * Pixley and Andrei [2019] J. H. Pixley and E. Y. Andrei, Science 365, 543 (2019).
  * Huang et al. [2019] T. Huang, L. Zhang, and T. Ma, Science Bulletin 64, 310 (2019).
  * Liu et al. [2019b] J. Liu, Z. Ma, J. Gao, and X. Dai, Physical Review X 9, 031021 (2019b).
  * Gonzalez and Stauber [2019] J. Gonzalez and T. Stauber, Physical review letters 122, 026801 (2019).
  * Lian et al. [2020] B. Lian, F. Xie, and B. A. Bernevig, Phys. Rev. B 102, 041402 (2020).
  * Padhi et al. [2020] B. Padhi, A. Tiwari, T. Neupert, and S. Ryu, Phys. Rev. Res. 2, 033458 (2020).
  * Wu and Das Sarma [2020] F. Wu and S. Das Sarma, Physical Review Letters 124 (2020), 10.1103/physrevlett.124.046403.
  * Bultinck et al. [2020a] N. Bultinck, S. Chatterjee, and M. P. Zaletel, Phys. Rev. Lett. 124, 166601 (2020a).
  * Bultinck et al. [2020b] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Phys. Rev. X 10, 031034 (2020b).
  * Hejazi et al. [2021] K. Hejazi, X. Chen, and L. Balents, Phys. Rev. Research 3, 013242 (2021).
  * Khalaf et al. [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Science advances 7, eabf5299 (2021).
  * Xie et al. [2020] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Phys. Rev. Lett. 124, 167002 (2020).
  * Julku et al. [2020] A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkilä, and P. Törmä, Physical Review B 101 (2020), 10.1103/physrevb.101.060505.
  * Kang and Vafek [2020] J. Kang and O. Vafek, Phys. Rev. B 102, 035161 (2020).
  * Chichinadze et al. [2020] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Phys. Rev. B 101, 224513 (2020).
  * Soejima et al. [2020] T. Soejima, D. E. Parker, N. Bultinck, J. Hauschild, and M. P. Zaletel, Phys. Rev. B 102, 205111 (2020).
  * König et al. [2020] E. J. König, P. Coleman, and A. M. Tsvelik, Phys. Rev. B 102, 104514 (2020).
  * Christos et al. [2020] M. Christos, S. Sachdev, and M. S. Scheurer, Proceedings of the National Academy of Sciences 117, 29543 (2020).
  * Lewandowski et al. [2021] C. Lewandowski, D. Chowdhury, and J. Ruhman, Phys. Rev. B 103, 235401 (2021).
  * Kwan et al. [2020] Y. H. Kwan, S. A. Parameswaran, and S. L. Sondhi, Phys. Rev. B 101, 205116 (2020).
  * Kwan et al. [2021a] Y. H. Kwan, Y. Hu, S. H. Simon, and S. A. Parameswaran, Phys. Rev. Lett. 126, 137601 (2021a).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Phys. Rev. Lett. 124, 097601 (2020).
  * Liu and Dai [2021] J. Liu and X. Dai, Phys. Rev. B 103, 035427 (2021).
  * Cea and Guinea [2020] T. Cea and F. Guinea, Phys. Rev. B 102, 045107 (2020).
  * Zhang et al. [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Phys. Rev. B 102, 035136 (2020).
  * Liu et al. [2021a] S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, Phys. Rev. Research 3, 013033 (2021a).
  * Da Liao et al. [2021] Y. Da Liao, J. Kang, C. N. Breiø, X. Y. Xu, H.-Q. Wu, B. M. Andersen, R. M. Fernandes, and Z. Y. Meng, Phys. Rev. X 11, 011014 (2021).
  * Eugenio and Dağ [2020] P. M. Eugenio and C. B. Dağ, SciPost Phys. Core 3, 15 (2020).
  * Huang et al. [2020] Y. Huang, P. Hosur, and H. K. Pal, Phys. Rev. B 102, 155429 (2020).
  * Calderón and Bascones [2020] M. J. Calderón and E. Bascones, Phys. Rev. B 102, 155149 (2020).
  * Ledwith et al. [2020] P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath, Phys. Rev. Research 2, 023237 (2020).
  * Repellin et al. [2020] C. Repellin, Z. Dong, Y.-H. Zhang, and T. Senthil, Phys. Rev. Lett. 124, 187601 (2020).
  * Abouelkomsan et al. [2020] A. Abouelkomsan, Z. Liu, and E. J. Bergholtz, Phys. Rev. Lett. 124, 106803 (2020).
  * Repellin and Senthil [2020] C. Repellin and T. Senthil, Phys. Rev. Research 2, 023238 (2020).
  * Vafek and Kang [2020] O. Vafek and J. Kang, Phys. Rev. Lett. 125, 257602 (2020).
  * Fernandes and Venderbos [2020] R. M. Fernandes and J. W. F. Venderbos, Science Advances 6 (2020), 10.1126/sciadv.aba8834.
  * Wilson et al. [2020] J. H. Wilson, Y. Fu, S. Das Sarma, and J. H. Pixley, Phys. Rev. Research 2, 023325 (2020).
  * Wang et al. [2021] J. Wang, Y. Zheng, A. J. Millis, and J. Cano, Phys. Rev. Research 3, 023155 (2021).
  * Bernevig et al. [2021a] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Physical Review B 103, 205411 (2021a), publisher: American Physical Society.
  * Song et al. [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Phys. Rev. B 103, 205412 (2021).
  * Bernevig et al. [2021b] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Physical Review B 103, 205413 (2021b), publisher: American Physical Society.
  * Lian et al. [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Physical Review B 103, 205414 (2021), publisher: American Physical Society.
  * Bernevig et al. [2021c] B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Physical Review B 103, 205415 (2021c), publisher: American Physical Society.
  * Xie et al. [2021] F. Xie, A. Cowsik, Z.-D. Song, B. Lian, B. A. Bernevig, and N. Regnault, Physical Review B 103, 205416 (2021), publisher: American Physical Society.
  * Zhang et al. [2021] X. Zhang, G. Pan, Y. Zhang, J. Kang, and Z. Y. Meng, Chinese Physics Letters 38, 077305 (2021).
  * Vafek and Kang [2021] O. Vafek and J. Kang, Phys. Rev. B 104, 075143 (2021).
  * Cha et al. [2021] P. Cha, A. A. Patel, and E.-A. Kim, Phys. Rev. Lett. 127, 266601 (2021), arXiv:2105.08069 [cond-mat.str-el] .
  * Sheffer and Stern [2021] Y. Sheffer and A. Stern, Phys. Rev. B 104, L121405 (2021), arXiv:2106.10650 [cond-mat.str-el] .
  * Kang et al. [2021] J. Kang, B. A. Bernevig, and O. Vafek, Phys. Rev. Lett. 127, 266402 (2021).
  * Hofmann et al. [2022] J. S. Hofmann, E. Khalaf, A. Vishwanath, E. Berg, and J. Y. Lee, Phys. Rev. X 12, 011061 (2022).
  * Thomson and Alicea [2021] A. Thomson and J. Alicea, Phys. Rev. B 103, 125138 (2021).
  * Kwan et al. [2021b] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Phys. Rev. X 11, 041063 (2021b).
  * Parker et al. [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Phys. Rev. Lett. 127, 027601 (2021).
  * Pathak et al. [2022] S. Pathak, T. Rakib, R. Hou, A. Nevidomskyy, E. Ertekin, H. T. Johnson, and L. K. Wagner, Phys. Rev. B 105, 115141 (2022).
  * Wagner et al. [2022] G. Wagner, Y. H. Kwan, N. Bultinck, S. H. Simon, and S. A. Parameswaran, Phys. Rev. Lett. 128, 156401 (2022).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Phys. Rev. Lett. 129, 047601 (2022).
  * Shi and Dai [2022] H. Shi and X. Dai, Phys. Rev. B 106, 245129 (2022).
  * Călugăru et al. [2023] D. Călugăru, M. Borovkov, L. L. H. Lau, P. Coleman, Z.-D. Song, and B. A. Bernevig, Low Temperature Physics 49, 640 (2023).
  * Xie et al. [2023] F. Xie, J. Kang, B. A. Bernevig, O. Vafek, and N. Regnault, Phys. Rev. B 107, 075156 (2023).
  * Hu et al. [2023a] H. Hu, B. A. Bernevig, and A. M. Tsvelik, Phys. Rev. Lett. 131, 026502 (2023a).
  * Hu et al. [2023b] H. Hu, G. Rai, L. Crippa, J. Herzog-Arbeitman, D. Călugăru, T. Wehling, G. Sangiovanni, R. Valentí, A. M. Tsvelik, and B. A. Bernevig, Phys. Rev. Lett. 131, 166501 (2023b).
  * Chou and Das Sarma [2023] Y.-Z. Chou and S. Das Sarma, Phys. Rev. Lett. 131, 026501 (2023).
  * Datta et al. [2023] A. Datta, M. J. Calderón, A. Camjayi, and E. Bascones, Nature Communications 14, 5036 (2023).
  * Chichinadze et al. [2022] D. V. Chichinadze, L. Classen, Y. Wang, and A. V. Chubukov, Phys. Rev. Lett. 128, 227601 (2022).
  * 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, Nature 574, 653 (2019).
  * Yankowitz et al. [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Sharpe et al. [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605–608 (2019).
  * Serlin et al. [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Science 367, 900 (2020).
  * Polshyn et al. [2019] H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Nature Physics 15, 1011–1016 (2019).
  * Xie et al. [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 572, 101 (2019).
  * Choi et al. [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, and et al., Nature Physics 15, 1174–1180 (2019).
  * Kerelsky et al. [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, and et al., Nature 572, 95–100 (2019).
  * Jiang et al. [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Nature 573, 91–95 (2019).
  * Saito et al. [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Nature Physics 16, 926–930 (2020).
  * 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, Nature 583, 375–378 (2020).
  * Liu et al. [2021b] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. Li, Science 371, 1261 (2021b).
  * Arora et al. [2020] H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, and et al., Nature 583, 379–384 (2020).
  * Cao et al. [2020] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020).
  * Wong et al. [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 582, 198–202 (2020).
  * Zondiner et al. [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, and et al., Nature 582, 203–208 (2020).
  * Nuckolls et al. [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 588, 610 (2020).
  * Choi et al. [2021] Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, et al., Nature 589, 536 (2021).
  * Saito et al. [2021a] Y. Saito, J. Ge, L. Rademaker, K. Watanabe, T. Taniguchi, D. A. Abanin, and A. F. Young, Nature Physics 17, 478 (2021a).
  * Das et al. [2021] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Nature Physics 17, 710 (2021).
  * Wu et al. [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Nature materials 20, 488 (2021).
  * Park et al. [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 592, 43 (2021).
  * Saito et al. [2021b] Y. Saito, F. Yang, J. Ge, X. Liu, T. Taniguchi, K. Watanabe, J. Li, E. Berg, and A. F. Young, Nature 592, 220 (2021b).
  * Rozen et al. [2021] A. Rozen, J. M. Park, U. Zondiner, Y. Cao, D. Rodan-Legrain, T. Taniguchi, K. Watanabe, Y. Oreg, A. Stern, E. Berg, et al., Nature 592, 214 (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, Phys. Rev. Lett. 127, 197701 (2021).
  * Lu et al. [2021] X. Lu, B. Lian, G. Chaudhary, B. A. Piot, G. Romagnoli, K. Watanabe, T. Taniguchi, M. Poggio, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Proceedings of the National Academy of Sciences 118 (2021), 10.1073/pnas.2100006118.
  * Grover et al. [2022] S. Grover, M. Bocarsly, A. Uri, P. Stepanov, G. Di Battista, I. Roy, J. Xiao, A. Y. Meltzer, Y. Myasoedov, K. Pareek, K. Watanabe, T. Taniguchi, B. Yan, A. Stern, E. Berg, D. K. Efetov, and E. Zeldov, Nature Physics 18, 885 (2022).
  * 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, Nature Physics 18, 633 (2022).
  * 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, Phys. Rev. Lett. 128, 217701 (2022).
  * Pulay [1980] P. Pulay, Chemical Physics Letters 73, 393 (1980).
  * Pulay [1982] P. Pulay, Journal of Computational Chemistry 3, 556 (1982).
  * Kudin et al. [2002] K. N. Kudin, G. E. Scuseria, and E. Cances, The Journal of chemical physics 116, 8255 (2002).