# Correlated insulators, semimetals, and superconductivity in twisted trilayer
graphene

Maine Christos  Department of Physics, Harvard University, Cambridge MA 02138,
USA  Subir Sachdev  Department of Physics, Harvard University, Cambridge MA
02138, USA  Mathias S. Scheurer  Institute for Theoretical Physics, University
of Innsbruck, Innsbruck A-6020, Austria

###### Abstract

Motivated by recent experiments indicating strong superconductivity and
intricate correlated insulating and flavor-polarized physics in mirror-
symmetric twisted trilayer graphene, we study the effects of interactions in
this system close to the magic angle, using a combination of analytical and
numerical methods. We identify asymptotically exact correlated many-body
ground states at all integer filling fractions $\nu$ of the flat bands. To
determine their fate when moving away from these fine-tuned points, we apply
self-consistent Hartree-Fock numerics and analytic perturbation theory, with
good agreement between the two approaches. This allows us to construct a phase
diagram for the system as a function of $\nu$ and the displacement field, the
crucial experimental tuning parameter of the system, and study the spectra of
the different phases. The phase diagram is dominated by a correlated
semimetallic intervalley coherent state and an insulating sublattice-polarized
phase around charge neutrality, $\nu=0$, with additional spin-polarization
being present at quarter ($\nu=-2$) or three quarter ($\nu=+2$) fillings of
the quasi-flat bands. We further study the superconducting instabilities
emerging from these correlated states, both in the absence and in the
additional presence of electron-phonon coupling, also taking into account
possible Wess-Zumino-Witten terms. In the experimentally relevant regime, we
find triplet pairing to dominate, possibly explaining the observed violation
of the Pauli limit. Our results have several consequences for experiments as
well as future theoretical work and illustrate the rich physics resulting from
the interplay of almost flat bands and dispersive Dirac cones in twisted
trilayer graphene.

###### Contents

  1. I Introduction
    1. I.1 Connection to experiment and phase diagram
    2. I.2 Organization of the paper
  2. II Model and symmetries
    1. II.1 Continuum model and symmetries
    2. II.2 Adding a Displacement Field
    3. II.3 Projected low-energy model
    4. II.4 Interactions and form factors
  3. III Exact groundstates at $D_{0}=0$
    1. III.1 Hamiltonian and construction of eigenstates
    2. III.2 Discussion of ground states
    3. III.3 Resultant candidate states
  4. IV Hartree-Fock numerics
    1. IV.1 Hartree-Fock Method
    2. IV.2 Band Structures
      1. IV.2.1 Self-Consistent Band Structures at $D_{0}=0$
      2. IV.2.2 Self-Consistent Band Structures for $D_{0}>0$
    3. IV.3 Energies and phase diagram
  5. V Analytical perturbation theory
    1. V.1 Ordering in the graphene-like bands
    2. V.2 Mixing between the bands
    3. V.3 Deforming the form factors
    4. V.4 Finite TBG bandwidth
    5. V.5 Comparison of energetics
    6. V.6 Breaking SU(2)${}_{+}\times$ SU(2)-
  6. VI Numerics for $\nu=\pm 2$
    1. VI.1 Procedure and results
    2. VI.2 Connection to experiment
  7. VII Superconductivity
    1. VII.1 Pairing in the presence of polarization
    2. VII.2 Pairing without polarization
    3. VII.3 Relevance of Dirac cones and topology
  8. VIII Conclusion and discussion
  9. A Hamiltonian for the system
    1. A.1 Continuum model and low-energy degrees of freedom
    2. A.2 Single-particle symmetries
    3. A.3 Interaction and form factors
  10. B Exact statements about the interacting groundstates
    1. B.1 Turning on the coupling between the subsystems
    2. B.2 Behavior of Goldstone modes
  11. C Hartree-Fock functional
    1. C.1 General form
    2. C.2 Numerical procedure
      1. C.2.1 How we fix the phases in numerics
      2. C.2.2 Iterative procedure
      3. C.2.3 Ground State dependence on grid size
    3. C.3 Varying Parameters
  12. D Numerics at $\nu=2$
    1. D.1 Stable Phases at $\nu=2$
    2. D.2 Connection to Experimental Phase Diagram
  13. E Band Structure and solution for $D_{0}\sim w_{1}$
  14. F Different contributions to the energy
    1. F.1 Deformation of the form factors
    2. F.2 Finite bandwidth of the TBG-like bands
      1. F.2.1 First order perturbation theory
      2. F.2.2 Second order perturbation theory: superexchange
    3. F.3 Ordering in the graphene-like bands
    4. F.4 Mixing between the TBG-like and graphene-like bands
    5. F.5 Hund’s Coupling
  15. G Superconducting instabilities
    1. G.1 Fluctuation-induced superconductivity
    2. G.2 Order parameter in the polarized phase
  16. H WZW terms in MSTG

##  I Introduction

Graphene-based moiré superlattice systems have attracted considerable interest
in the last few years, motivated by the strongly correlated physics they
display MacDonald (2019); Andrei and MacDonald (2020); Kennes  _et al._
(2021); Balents  _et al._ (2020); Scheurer (2019); Zaletel (2021). While
originally driven by the experimental realization of near-magic-angle twisted
bilayer graphene (TBG) Cao  _et al._ (2018a, b), related strongly-correlated
moiré superlattices have emerged, such as twisted double-bilayer graphene Shen
_et al._ (2020); Liu  _et al._ (2020); Cao  _et al._ (2020); Burg  _et al._
(2019) and $\mathrm{ABC}$-trilayer graphene on hexagonal boron nitride Chen
_et al._ (2019, 2019, 2020). All of these systems exhibit low-energy bands
which are energetically separated from the rest of the spectrum and can be
tuned to be flat Dos Santos  _et al._ (2007); Bistritzer and MacDonald (2011);
Dos Santos  _et al._ (2012), enhancing the impact of correlations. The
resulting correlated phenomena include interaction-induced insulating states
Cao  _et al._ (2018a); Shen  _et al._ (2020); Liu  _et al._ (2020); Cao  _et
al._ (2020); Chen  _et al._ (2019), superconductivity Cao  _et al._ (2018b);
Shen  _et al._ (2020); Liu  _et al._ (2020); Chen  _et al._ (2019), and
nematic order Kerelsky  _et al._ (2019); Scheurer (2019); Cao  _et al._
(2020); Rubio-Verdú  _et al._ (2020); Samajdar  _et al._ (2021), which are
also integral parts of the phase diagrams of paradigmatic strongly correlated
materials such as the cuprates. Furthermore, TBG has demonstrated that
graphene moiré systems can also exhibit additional cascades of transitions
that “reset” the band structure in an entire range of electron filling
fractions Zondiner  _et al._ (2020); Wong  _et al._ (2020). This is likely
related to the polarization of certain combinations of the internal “flavor”
quantum numbers of the electrons in the quasi-flat bands Christos  _et al._
(2020); Kang  _et al._ (2021), with interesting consequences for
superconductivity Scheurer and Samajdar (2020); Christos  _et al._ (2020).

In addition to the twist angle between the layers as a tuning knob in TBG, the
spectrum of twisted double-bilayer graphene can be efficiently tuned by
applying a perpendicular electric displacement field $D_{0}$. Notwithstanding
the interesting consequences for the correlated physics in the system Shen
_et al._ (2020); Liu  _et al._ (2020); Cao  _et al._ (2020); Burg  _et al._
(2019); Rubio-Verdú  _et al._ (2020); Zhang  _et al._ (2021); Liu  _et al._
(2021a), including the possibility of electrical control of the nematic
director Samajdar  _et al._ (2021), its superconducting properties Shen  _et
al._ (2020); Liu  _et al._ (2020) have been found to be more fragile than in
TBG. Fortunately, mirror symmetric twisted trilayer graphene (MSTG), which
consists of three layers of graphene with alternating relative twist angles,
see Fig. 2(a), has very recently been realized experimentally Park  _et al._
(2021); Hao  _et al._ (2021); Cao  _et al._ (2021) and combines the best of
both worlds: it can be tuned significantly by applying a perpendicular
displacement field $D_{0}$ while exhibiting strong and reproducible
superconductivity. These experiments indicate that MSTG also exhibits both
interaction-induced resistive states around integer filling fractions $\nu$ as
well as the aforementioned spontaneous band resetting in an extended region of
$\nu$ with $2\lesssim|\nu|\lesssim 3$, akin to the cascades in TBG Zondiner
_et al._ (2020); Wong  _et al._ (2020). As a function of $D_{0}$ and $\nu$,
the largest superconducting region is found to emerge out of this
reconstructed normal state. Remarkably, Ref. Cao  _et al._ , 2021 finds that
the superconductor can sustain an in-plane magnetic field much larger than
simple estimates of the Pauli limit. One additional crucial difference between
MSTG and all of the other graphene moiré systems mentioned above is that MSTG
not only exhibits quasi-flat bands but, at the same time, dispersive Dirac
cones Khalaf  _et al._ (2019); Carr  _et al._ (2020); Mora  _et al._ (2019);
Park  _et al._ (2021); Hao  _et al._ (2021); Cǎlugǎru  _et al._ (2021); Shin
_et al._ (2021); Lei  _et al._ (2021). The interplay between these different
types of bands and interactions likely gives rise to rich physics but also
provides theoretical challenges.

So far, only very few theoretical studies of interactions in MSTG exist
Ramires and Lado (2021); Choi and Choi (2021); Lake and Senthil (2021); Qin
and MacDonald (2021); Fischer  _et al._ (2021); Chou  _et al._ (2021) and a
systematic understanding of the nature and origin of possible particle-hole
and superconducting instabilities as a function of $\nu$ and $D_{0}$ is still
missing. The goal of this work is to help fill this gap by providing a
detailed theoretical study of electron-electron interactions in the normal
state of MSTG, allowing for a large class of possible instabilities, which is
then used to analyze the order parameter and origin of the different
superconducting phases of the system. In order to tackle the challenges
associated with the simultaneous presence of flat and dispersive bands, we use
a combination of analytical and numerical approaches.

More specifically, we start from the limit without displacement field,
$D_{0}=0$, where the non-interacting band structure of the system is just that
of TBG and single-layer graphene Khalaf  _et al._ (2019); Cǎlugǎru  _et al._
(2021); Carr  _et al._ (2020); Shin  _et al._ (2021); Mora  _et al._ (2019);
Lei  _et al._ (2021); these two “subsystems” are, however, coupled by the
Coulomb interaction. In the limit where the TBG bands are perfectly flat, we
construct exact eigenstates of the interacting Hamiltonian of MSTG at all
integer $\nu$, which are also shown to be groundstates for a finite range of
the strength of the interaction between the two subsystems. This analysis
shows that, in certain limits, also the interacting groundstates of MSTG are
given by those of TBG and single-layer graphene. This is consistent with
experiment Park  _et al._ (2021); Hao  _et al._ (2021), observing a graphene
Dirac cone at small $D_{0}$.

To be able to address $D_{0}\neq 0$ and more realistic system parameters, we
use Hartree-Fock (HF) numerics and analytic perturbation theory. Our HF
approach is motivated by the success of this approach in TBG Xie and MacDonald
(2020); Bultinck  _et al._ (2019); Liu  _et al._ (2021b); Liao  _et al._
(2021); Bultinck  _et al._ (2020); Zhang  _et al._ (2020) and the fact that
the exact groundstates at $D_{0}=0$ are Slater-determinant states (in the TBG
sector). In the analytic perturbation theory, we start from the exact
groundstates and study the $D_{0}$-induced deformations of the interaction
matrix elements and band structure, the possible ordering in the graphene
bands, the mixing between the graphene and TBG sectors for $D_{0}\neq 0$, as
well as the finite bandwidth of the TBG bands as perturbations. This
complements our HF numerics, as it allows us to pinpoint the novel energetic
contributions in MSTG as compared to TBG and serves as an important validation
of our numerics—in particular, concerning the fate of the graphene Dirac
cones, which provides a challenge to any numerical study as they are
energetically degenerate with the TBG bands only in a small fraction of the
moiré Brillouin zone (MBZ). In our analytics, we take advantage of this
fraction being small and use it as an expansion parameter. Overall, we find
good agreement between the HF numerics and the perturbation theory, both
confirming the form of the exact eigenstates at $D_{0}=0$ (but with realistic
system parameters); discrepancies between the approaches are traced back to
the impact on the energetics coming from additional remote bands only taken
into account in the numerics.

Having established the correlated nature and possible particle-hole
instabilities in the normal state of MSTG as a function of $D_{0}$ and $\nu$,
we analyze the consequences for superconductivity. We study superconductivity
both in the additional presence and absence of flavor polarization. Motivated
by experiments Stepanov  _et al._ (2020); Liu  _et al._ (2021c); Saito  _et
al._ (2020) on TBG, which indicate that electron-phonon coupling is important
for pairing, we follow Ref. Samajdar and Scheurer, 2020 and assume that
electron-phonon coupling stabilizes superconductivity but leaves singlet and
triplet almost degenerate; informed by results for particle-hole instabilities
in MSTG, we can then investigate which of the two will be favored due to
additional particle-hole fluctuations. We also discuss purely electronic
pairing and comment on the relevance of Dirac cones and associated Wess-
Zumino-Witten (WZW) terms, as studied previously in TBG Khalaf  _et al._
(2021); Christos  _et al._ (2020).

Figure 1: Schematic phase diagram for MSTG based on our analysis, as a
function of displacement field $D_{0}$ and filling fraction $\nu$. Due to
particle-hole symmetry of our model, the phase diagram only depends on
$|\nu|$. We show the form of the dominant particle-hole instabilities, the
intervalley coherent (IVC) semimetallic phase and a sublattice polarized (SLP)
insulator at $\nu=0$ (charge neutrality); the IVC and SLP state coexist with
spin-polarization (SP) around $|\nu|=2$. We discuss three different types of
superconducting regimes, labelled SC I–III. We indicate in black their
respective spin structure for the expected sign of the intervalley Hund’s
coupling $J_{H}$ resulting from Coulomb interactions alone [with opposite sign
in gray].

###  I.1 Connection to experiment and phase diagram

To concisely summarize our main findings, we here briefly discuss their
consequences for the phase diagram of MSTG, see Fig. 1, and their relation to
experiment Park  _et al._ (2021); Hao  _et al._ (2021); Cao  _et al._ (2021).
Based on our exact groundstates, perturbation theory, and HF numerics, the
(for $\nu\neq 0$ slightly doped) semimetallic nature of the graphene sector is
expected to be stable against interactions for $D_{0}=0$. At charge
neutrality, $\nu=0$, the leading instability is found to be an intervalley
coherent state (IVC) in the TBG sector. While the two sectors start to mix
when $D_{0}\neq 0$, the semimetallic character of the bands is retained; the
Landau fan sequence of this state, $\pm 2,\pm 6,\pm 10,\dots$, is consistent
with that observed in Ref. Hao  _et al._ , 2021. Once $D_{0}$ increases beyond
a certain critical value $D_{0}^{c}$ [which depends on the relaxation
parameter $w_{0}/w_{1}$ in Eq. (3)], the IVC state transitions into a
sublattice-polarized (SLP) phase, where all Dirac cones are gapped out,
leading to an insulator. For realistic parameters, we expect $D_{0}^{c}$ to be
of the order of or potentially slightly larger than the range of the
experimentally applied fields Park  _et al._ (2021); Hao  _et al._ (2021); Cao
_et al._ (2021).

At $|\nu|=2$, we find self-consistent HF solutions with spin polarization,
which co-exists with additional particle-hole instabilities. At small $D_{0}$,
the leading additional symmetry-breaking state is again an IVC phase. Here,
the TBG bands are completely spin polarized while the graphene cones’
polarization is only partial and approaches zero as $D_{0}\rightarrow 0$. The
spin polarization is our proposed mechanism for the experimentally observed
band reset and reduced flavor number Park  _et al._ (2021); Hao  _et al._
(2021). At larger $D_{0}$, the IVC transitions into an SLP state. In Sec.
VI.2, we discuss how the bandstructures we find can give rise to the
experimentally observed Park  _et al._ (2021) additional sign change of the
Hall density at sufficiently large $D_{0}$. Our obtained behavior at
$|\nu|=0,2$ also agrees with another experimental feature: the conductance is
suppressed Park  _et al._ (2021); Hao  _et al._ (2021) at $\nu=0$ for both
$D_{0}=0$ and $D_{0}>0$ while the suppression at $|\nu|=2$ only sets in above
a finite critical value of $D_{0}$, which is consistent with the IVC (SLP)
being semimetallic (insulating) at $\nu=0$ and metallic (semimetallitc) at
$|\nu|=2$.

As discussed in detail in Sec. V.6, the relative spin-orientation of these
orders in the two valleys depends on the sign of the intervalley Hund’s
coupling $J_{H}$ (see Eq. (30)), which is not part of the Hamiltonian for our
Hartree-Fock computations; the value or sign of $J_{H}$ is not determined in
our computations. For $J_{H}<0$, which is expected for Coulomb interactions,
both the SLP and SP have the same spin polarization in the two valleys, while
the IVC order parameter does not carry any spin.

As indicated in blue in Fig. 1, we investigate three different regimes of
superconducting phases. For superconductivity (SC I) in the presence of spin
polarization with or without additional IVC order, the order parameter will be
a non-unitary triplet for $J_{H}<0$ and an admixed singlet-triplet phase
Scheurer and Samajdar (2020) for $J_{H}>0$. For superconducting phases (SC II)
close to or in the range $2\lesssim|\nu|\lesssim 3$ but not coexisting with
additional particle-hole instabilities, we expect triplet (a singlet-triplet
admixed phase) to dominate for $J_{H}<0$ ($J_{H}>0$). All of these
superconducting states are consistent Scheurer and Samajdar (2020) with the
observed Cao  _et al._ (2021) behavior in low to intermediate magnetic fields.
Finally, our analysis indicates that a superconducting phase (if present)
close charge neutrality (SC III) should be in a singlet state for $J_{H}<0$.

###  I.2 Organization of the paper

The remainder of the paper is organized as follows. In Sec. II, we introduce
the model for MSTG we study, its symmetries, and establish the basic notation
used in this work. We begin our discussion of correlated physics in the limit
of vanishing displacement field in Sec. III, where we construct exact
groundstates. These build the starting point for our HF numerics, see Sec. IV
for $\nu=0$ and Sec. VI for $|\nu|=2$, and our analytical perturbation theory,
presented in Sec. V. Furthermore, superconducting instabilities will be
analyzed in Sec. VII. Finally, Sec. VIII contains a short summary and
discussion of results. Details of our analysis are provided in a set of
appendices, labelled Appendix A through H.

##  II Model and symmetries

In this section, we introduce the interacting model we consider in this work,
which consists of a continuum model to describe the moiré bands supplemented
by Coulomb repulsion. We further discuss its symmetries and define the
notation used in the remainder of the paper.

###  II.1 Continuum model and symmetries

MSTG is constructed from three parallel sheets of graphene where the top
($l=1$) and bottom ($l=3$) graphene layers are aligned with one another and
the middle ($l=2$) layer of graphene is twisted at a relative angle $\theta$
with respect to the top and bottom layers, as shown in Fig. 2(a). To compute
the non-interacting band structure, we employ a continuum-model description,
which is just a three-layer extension of the frequently applied continuum
model for TBG Dos Santos  _et al._ (2007); Bistritzer and MacDonald (2011);
Dos Santos  _et al._ (2012). Denoting the electronic creation operator for an
electron at position $\boldsymbol{r}\in\mathbbm{R}^{2}$, in sublattice
$\rho=1,2$, layer $l=1,2,3$, valley $\eta=\pm$, and of spin
$s=\uparrow,\downarrow$ by $c^{\dagger}_{\boldsymbol{r};\rho,l,\eta,s}$ the
Hamiltonian reads in the absence of a displacement field as Khalaf  _et al._
(2019); Carr  _et al._ (2020); Mora  _et al._ (2019); Cǎlugǎru  _et al._
(2021)

$\begin{split}H^{\text{Full}}_{0,1}=&\int_{\boldsymbol{r}}c^{\dagger}_{\boldsymbol{r};\rho,l,\eta,s}\left(-iv_{F}(\boldsymbol{\rho}^{\eta}_{\theta_{l}})^{\phantom{\dagger}}_{\rho,\rho^{\prime}}\cdot\boldsymbol{\nabla}\right)c^{\phantom{\dagger}}_{\boldsymbol{r};\rho^{\prime},l,\eta,s}\\\
&+\int_{\boldsymbol{r}}\left[c^{\dagger}_{\boldsymbol{r};\rho,l,\eta,s}(T_{\eta,\boldsymbol{r}})_{\rho,\rho^{\prime}}c^{\phantom{\dagger}}_{\boldsymbol{r};\rho^{\prime},l+1,\eta,s}+\text{H.c.}\right],\end{split}$
(1)

where summation over repeated indices is implied,
$\theta_{l}=(-1)^{l}\theta/2$,
$(\boldsymbol{\rho}^{\eta=+}_{\theta})_{j}=e^{i\theta\rho_{z}/2}\rho_{j}e^{-i\theta\rho_{z}/2}$,
with $\rho_{j}$ representing Pauli matrices in sublattice space (throughout we
will use the same symbol for Pauli matrices and their indices, i.e., $s_{j}$,
$\eta_{j}$ are Pauli matrices in spin and valley space), and
$\boldsymbol{\rho}^{\eta=-}_{\theta}=-(\boldsymbol{\rho}^{\eta=+}_{\theta})^{*}$.
While the first term in Eq. (1) describes the Dirac cones of the two valleys,
$\eta=\pm$, of each individual graphene layer in the continuum expansion
around their respective Kg and K${}^{\prime}_{g}$ points, the second line
captures the tunneling between adjacent graphene layers parameterized by
$T_{\eta,\boldsymbol{r}}$ (the direct hopping process between the outer two
layers is neglected). As a consequence of the moiré superlattice, these
tunneling matrix elements are modulated spatially. As is common, we focus on
the lowest moiré-lattice harmonics for which symmetry allows for only two
independent parameters, $w_{0}$ and $w_{1}$, in $T_{\eta,\boldsymbol{r}}$. It
can be written as Bistritzer and MacDonald (2011)

$T_{+,\boldsymbol{r}}=e^{-i\boldsymbol{q}_{1}\boldsymbol{r}}\left[\mathcal{T}_{1}+\mathcal{T}_{2}e^{-i\boldsymbol{G}_{1}^{\text{M}}\boldsymbol{r}}+\mathcal{T}_{3}e^{-i(\boldsymbol{G}_{1}^{\text{M}}+\boldsymbol{G}_{2}^{\text{M}})\boldsymbol{r}}\right]$
(2)

and $T_{-,\boldsymbol{r}}=(T_{+,\boldsymbol{r}})^{*}$. Here
$\boldsymbol{q}_{1}=k_{\theta}(0,-1)$,
$k_{\theta}=2|\text{K}_{g}|\sin(\theta/2)$, connects the K and K′ points of
the moiré lattice,
$\boldsymbol{G}_{1}^{\text{M}}=-\sqrt{3}k_{\theta}(1,\sqrt{3})^{T}/2$ and
$\boldsymbol{G}_{2}^{\text{M}}=\sqrt{3}k_{\theta}(1,0)^{T}$ are the basis
vectors of the reciprocal lattice (RL) of the moiré lattice,
$\text{RL}:=\\{\sum_{j=1,2}n_{j}\boldsymbol{G}_{j}^{\text{M}},\,n_{j}\in\mathbbm{Z}\\}$.
The matrices $\mathcal{T}_{j}$ in sublattice space in Eq. (2) are given by

$\mathcal{T}_{j}=w_{0}\rho_{0}+w_{1}\begin{pmatrix}0&\omega^{j-1}\\\
\omega^{-(j-1)}&0\end{pmatrix},\quad\omega=e^{-i\frac{2\pi}{3}}.$ (3)

While rigidly rotated graphene layers correspond to $w_{0}=w_{1}$, lattice
relaxation Nam and Koshino (2017); Carr  _et al._ (2019) leads to $w_{0}\neq
w_{1}$. We will here consider $w_{0}$ as a free parameter and study the
physics as a function of it.

As a result of the moiré modulations in Eq. (2), the tunneling matrix elements
between the layers couple momenta related by RL vectors, which reconstructs
the graphene cones of the first line in Eq. (1), leading to a (technically
infinite) set of bands; the band energies will be denoted by
$\epsilon_{n,\eta}(\boldsymbol{k})$, $n\in\mathbbm{Z}$, where $\boldsymbol{k}$
is in the first MBZ and we already used that Eq. (1) is diagonal in the
valley, $\eta$, and trivial in the spin index, $s$; this makes the band
energies independent of spin and allows to label them by their valley index.

MSTG is distinguished from other moiré systems, such as TBG or twisted double-
bilayer graphene, by a reflection symmetry $\sigma_{h}$ under exchange of the
top and bottom layers of MSTG, with action
$\sigma_{h}:\,c_{\boldsymbol{r};\rho,l,\eta,s}\rightarrow
c_{\boldsymbol{r};\rho,\sigma_{h}(l),\eta,s}$ where $\sigma_{h}(2)=2$,
$\sigma_{h}(1)=3$, and $\sigma_{h}(3)=1$. Since it acts trivially in all
internal indices ($\rho,\eta,s$) and $\boldsymbol{r}$, the Hamiltonian (1) can
be decomposed into sectors with different mirror eigenvalues, $\sigma_{h}=\pm
1$, by performing a unitary transformation, $V\in\text{U}(3)$, in layer space
only,
$c_{\boldsymbol{r};\rho,l,\eta,s}=V_{l,\ell}\psi_{\boldsymbol{r};\rho,\ell,\eta,s}$.
As was pointed out before Khalaf  _et al._ (2019) and detailed in Appendix
A.1, the Hamiltonian in the mirror-even sector ($\sigma_{h}=+1$, $\ell=1,2$)
turns out to be that of TBG with an interlayer hopping renormalized by a
factor of $\sqrt{2}$. In the mirror-odd sector ($\sigma_{h}=-1$, $\ell=3$) the
Hamiltonian is that of single-layer graphene (without any moiré modulation).
This can be clearly seen in the spectrum shown in Fig. 2(b), which exhibits
two almost-flat TBG-bands (shown in red) per spin and valley and
unreconstructed Dirac cones at K (K′) associated with the graphene sector of
valley $\eta=+$ ($\eta=-$). As such, there are three Dirac cones at K (and
three at K′), one belonging to the graphene sector of MSTG in a single valley
only and two Dirac crossings belonging to the TBG sector of both valleys. As
can be seen, there are additional remote bands of both subspaces due to back-
folding into the MBZ.

In addition to reflection symmetry, the Hamiltonian in Eq. (1) has more point
symmetries, forming the group $C_{6h}$: the model is invariant under three-
fold rotational symmetry $C_{3z}$, with action,
$C_{3z}:\,c_{\boldsymbol{r}}\rightarrow
e^{i\frac{2\pi}{3}\rho_{z}\eta_{z}}c_{C_{3z}\boldsymbol{r}}$, which is also an
exact symmetry of the moiré lattice in Fig. 2(a). While not an exact lattice
symmetry, two-fold rotation perpendicular to the plane, $C_{2z}$, is a good
approximate symmetry for small twist angles; in fact, it is an exact symmetry
of the continuum in Eq. (1) which can be verified by applying its action,
$C_{2z}:\,c_{\boldsymbol{r}}\rightarrow\eta_{x}\rho_{x}c_{-\boldsymbol{r}}$.

Besides these point symmetries, the model also exhibits the following exact
internal symmetries: since there is no coupling between the valleys, it is
invariant under a valley U(1) transformation,
$\text{U}(1)_{v}:\,c_{\boldsymbol{r}}\rightarrow
e^{i\eta_{z}\varphi}c_{\boldsymbol{r}}$. In combination with the absence of
spin-orbit coupling, it is further invariant under the separate spin-rotation
in each valley, forming the group $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$
with $\text{SU}(2)_{\pm}$ acting as $c_{\boldsymbol{r}}\rightarrow
e^{i\boldsymbol{\varphi}\cdot\boldsymbol{s}(\eta_{0}\pm\eta_{z})/2}c_{\boldsymbol{r}}$.
Furthermore, there is time-reversal symmetry, which is associated with the
anti-unitary operator $\Theta$ with $\Theta
c_{\boldsymbol{r}}\Theta^{\dagger}=\eta_{x}c_{\boldsymbol{r}}$; unless stated
otherwise, we will always refer to this form of (spinless) time-reversal
symmetry throughout the text, although combinations with spin-rotations
(spinful time-reversal), $\Theta_{s}$, and with $\text{U}(1)_{v}$ rotations,
$\widetilde{\Theta}$, will play a role further below. Since the mirror
symmetry protects any mixing between the graphene and TBG subspaces and
$\text{U}(1)_{v}$ any mixing between the two valleys, the aforementioned Dirac
crossings at K and K′ are protected by the combination $C_{2z}\Theta$, exactly
as in graphene and TBG.

Focusing on specific limits, there are also additional internal symmetries,
similar to TBG Tarnopolsky  _et al._ (2019); Bultinck  _et al._ (2019);
Bernevig  _et al._ (2020): when $w_{0}=0$, $H^{\text{Full}}_{0,1}$ in Eq. (1)
changes its sign, $H^{\text{Full}}_{0,1}\rightarrow-H^{\text{Full}}_{0,1}$, if
we apply the chiral symmetry operator
$C:\,c_{\boldsymbol{r}}\rightarrow\rho_{z}L_{C}c_{\boldsymbol{r}}$ where
$L_{C}$ is a matrix in layer space given by
$L_{C}=V\text{diag}(1,1,-1)V^{\dagger}$. Therefore, $w_{0}=0$ will be referred
to as the chiral limit Tarnopolsky  _et al._ (2019). Among other consequences
for the interaction terms to be discussed below, the band structure has to be
symmetric about zero energy at every given momentum $\boldsymbol{k}$ when
$w_{0}=0$, i.e.,
$\text{spec}_{\boldsymbol{k},\eta}=-\text{spec}_{\boldsymbol{k},\eta}$, with
$\text{spec}_{\boldsymbol{k},\eta}=\\{\epsilon_{n,\eta}(\boldsymbol{k}),n\in\mathbb{Z}\\}$.

Finally, when the additional rotation of the sublattice matrices in the first
line of Eq. (1) is neglected,
$\boldsymbol{\rho}^{\eta}_{\theta_{l}}\rightarrow\boldsymbol{\rho}^{\eta}_{\theta=0}$,
the “unitary particle-hole symmetry”, previously discussed in TBG Bultinck
_et al._ (2019); Bernevig  _et al._ (2020), can be extended to MSTG: defining
the unitary operator $P$ which acts as
$P\psi_{\boldsymbol{r};l=1}P^{\dagger}=\eta_{z}\psi_{-\boldsymbol{r};l=2}$,
$P\psi_{\boldsymbol{r};l=2}P^{\dagger}=-\eta_{z}\psi_{-\boldsymbol{r};l=1}$,
and
$P\psi_{\boldsymbol{r};l=3}P^{\dagger}=i\rho_{y}\eta_{y}\psi_{-\boldsymbol{r};l=3}$,
one finds $PH^{\text{Full}}_{0,1}P^{\dagger}=-H^{\text{Full}}_{0,1}$, showing
that the spectrum must obey
$\text{spec}_{\boldsymbol{k}}=-\text{spec}_{-\boldsymbol{k}}$ where
$\text{spec}_{\boldsymbol{k}}=\cup_{\eta}\text{spec}_{\boldsymbol{k},\eta}$.
Since we will be focusing on small twist angles $\theta$ below, we will always
assume that $\boldsymbol{\rho}^{\eta}_{\theta_{l}}$ has been replaced by
$\boldsymbol{\rho}^{\eta}_{\theta=0}$ in Eq. (1) when analyzing symmetries. We
will also later set $\theta=0$ in our numerics.

For the detailed form of the model and the action of the symmetries in
momentum space we refer the reader to Appendix A.1 and A.2, respectively. A
brief list of the symmetries of $H_{0}$ and when they apply can be found in
Table 1.

The primary interaction term we shall consider, in Eq. (8), will preserve the
symmetries of the continuum model. However, we will also consider the
consequences of a Hund’s coupling in Eq. (30), which will break the
$\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ symmetry down to a more realistic
$\text{SU}(2)_{s}$ (the group of simultaneous spin rotations in the two
valleys).

###  II.2 Adding a Displacement Field

Motivated by recent experimental discoveries of electric-field-tunable
correlated physics in MSTG Park  _et al._ (2021); Hao  _et al._ (2021); Cao
_et al._ (2021), we now extend the model for MSTG in Eq. (1) and include a
perpendicular displacement field $D_{0}$. Finite $D_{0}$ results in a
potential difference between the top and bottom layers of graphene.
Suppressing all indices except layer $l$, a displacement field is represented
in the continuum model as

$H^{\text{Full}}_{0,2}=D_{0}\sum_{\boldsymbol{r}}\left[c^{\dagger}_{\boldsymbol{r};l=1}c^{\phantom{\dagger}}_{\boldsymbol{r};l=1}-c^{\dagger}_{\boldsymbol{r};l=3}c^{\phantom{\dagger}}_{\boldsymbol{r};l=3}\right]$
(4)

and the full Hamiltonian becomes
$H^{\text{Full}}_{0}=H^{\text{Full}}_{0,1}+H^{\text{Full}}_{0,2}$. The
displacement field $D_{0}$ breaks the mirror symmetry $\sigma_{h}$, resulting
in hybridization between bands in the TBG and graphene sectors; this
hybridization can only take place between graphene and TBG bands of the same
valley and spin flavor, as $\text{U}_{v}(1)$ and spin-rotation symmetry [in
fact the full $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$] are preserved by
$H^{\text{Full}}_{0,2}$.

(a)

(b)

(c)

Figure 2:  (a) Schematic of the geometry of MSTG. The top and bottom layers
shown in red and green respectively are aligned, while the middle layer (blue)
is twisted at relative angle $\theta$ to the top and bottom layers. (b) Band
structure of MSTG at $D_{0}=0$ (left) and $D_{0}>0$ (right) for $w_{1}=124$
meV and twist angle $\theta=1.53^{\circ}$ along the one-dimensional cut
through the MBZ shown in (c). We represent the spin-degenerate bands in the
two different valleys $\eta=\pm$ with solid ($\eta=+$) and dotted ($\eta=-$)
lines. At $D_{0}=0$, the band structure is given by that of TBG and single-
layer graphene and we have colored the (two per spin and valley) quasi-flat
TBG bands red. For $D_{0}>0$, the bands of these two subsystems mix but we
continue to label the two bands per spin and valley closest to the Fermi as
“TBG-like” bands, indicated in red.

Furthermore, when $D_{0}$ is finite, the system retains $C_{2z}$, $C_{3z}$,
leading to the point group $C_{6}$. Time-reversal symmetry, $\Theta$, and, in
the limit $w_{0}=0$, chiral symmetry, $C$, also persist. However, the unitary
particle hole symmetry $P$ is no longer an exact symmetry for finite $D_{0}$.

At the K (K′) point of the MBZ, the Dirac crossings which belong to the valley
$\eta=+$ ($\eta=-$) will hybridize as $D_{0}$ increases with the graphene
bands of the same valley which also exhibit a Dirac cone at the Fermi level;
this results in the Dirac points, which were pinned at the Fermi level at
$D_{0}=0$, being pushed away from the Fermi level. As $C_{2z}\Theta$ symmetry
is preserved for $D_{0}\neq 0$, these Dirac crossings cannot be gapped out. As
a consequence of the absence of a graphene cone to hybrize with, the Dirac
point of TBG in the other valley, $\eta=-$ ($\eta=+$), remains pinned at the
Fermi level as $D_{0}$ increases. These features can be clearly observed in
the band structure of MSTG with finite displacement field shown in Fig. 2(b).

###  II.3 Projected low-energy model

In order to make the analytical and numerical study of interactions feasible
in the system, we restrict our Hamiltonian to a finite set of bands in the
vicinity of the Fermi level. Going forward, we will denote the two bands
closest to the Fermi level per spin and valley as TBG-like bands and the next
two closest bands per spin and valley as graphene-like bands. All other bands
at energies further away from the Fermi level will be referred to as remote
bands. When $D_{0}=0$, the TBG-like bands are labeled by mirror symmetry
eigenvalue $\sigma_{h}=+1$ and can be identified exactly with the two flat
bands of TBG. The graphene-like bands at $D_{0}=0$ are identified exactly with
the bands of the continuum model of graphene in the vicinity of the K and K′
points. For finite $D_{0}$, the bands with different eigenvalues hybridize but
we will retain our naming convention for the four bands closest to the Fermi
level, i.e., the bands indicated in red in Fig. 2 are TBG-like bands.

Note that away from the K points, the graphene-like bands are technically
identified with the first remote bands of TBG with $\sigma_{h}=+1$ for
$D_{0}=0$ within this convention. In our analytical calculations below,
however, we will restrict the analysis to the TBG-like bands in the full MBZ
and to the graphene-like bands only in the vicinity of the K and K′ points
where they retain their graphene-like identity. Denoting the associated
creation operators in these two sets of bands by
$b^{\dagger}_{\boldsymbol{k};p,\eta,s}$ and
$g^{\dagger}_{\boldsymbol{k};p,\eta,s}$, where $\boldsymbol{k}\in\text{MBZ}$
is the momentum, $p=+$ ($p=-$) labels the upper (lower) band in each sector,
$\eta$ is the valley, and $s$ the spin quantum number of the electrons, the
non-interacting Hamiltonian projected to these bands becomes

$\displaystyle\begin{split}H_{0}&=\sum_{\boldsymbol{k}\in\text{MBZ}}W_{\text{TBG}}\epsilon_{(b,p),\eta}(\boldsymbol{k})b^{\dagger}_{\boldsymbol{k};p,\eta,s}b^{\phantom{\dagger}}_{\boldsymbol{k};p,\eta,s}\\\
&+\sum_{\boldsymbol{k}\in\text{MBZ}^{\prime}}\epsilon_{(g,p),\eta}(\boldsymbol{k})g^{\dagger}_{\boldsymbol{k};p,\eta,s}g^{\phantom{\dagger}}_{\boldsymbol{k};p,\eta,s}.\end{split}$
(5)

Here the extra prime in MBZ′ in the second line indicates that we restrict the
graphene-like degrees of freedom to the vicinity of the K and K′ points in our
analytical calculations. Note that we further introduced the dimensionless
parameter $W_{\text{TBG}}\in[0,1]$ that will allow us to organize the
perturbation theory of Sec. V.4 in the TBG-like bandwidth. While the physical
system corresponds to $W_{\text{TBG}}=1$, we define the flat limit as
$W_{\text{TBG}}=0$, which will play an important role below. Note that the
Hamiltonian in Eq. (1) realizes the flat limit exactly,
$\epsilon_{(b,p),\eta}(\boldsymbol{k})=0$, only when $w_{0}=0$, $\theta$ is at
the magic angle, and when $D_{0}=0$.

In our numerics, we will keep all momenta of both sets of (in total $16$)
bands in Eq. (5) and continue to use $g^{\dagger}_{\boldsymbol{k};p,\eta,s}$
as creation operators for the second lowest set of bands throughout the MBZ.
To check convergence, we will also study the impact of adding additional
remote bands, as described in Appendix C.3.

In order to fix the phases of the wavefunctions of the TBG-like and graphene-
like bands, we specify how they transform under the action of the various
discrete symmetries discussed above or, equivalently how they act on the
electronic operators $b^{\dagger}_{\boldsymbol{k};p,\eta,s}$ and
$g^{\dagger}_{\boldsymbol{k};p,\eta,s}$. While a thorough discussion can be
found in Appendix A.2, these representations are summarized in Table 1.

Table 1: Summary of discrete single-particle symmetries (denoted $S$) of our model and when they apply (last column). For convenience of the reader and future reference, we show a redundant set of symmetries. We also indicate whether they are unitary or anti-unitary symmetries (second column) and whether they commute $+$ or anti-commute $-$ (third column) with the non-interacting Hamiltonian $H^{\text{Full}}_{0}$. Finally, we specify their action on the low-energy field operators in order to fix the phase of the Bloch states, using $\sigma_{j}$ to denote Pauli matrices acting on the upper/lower bands subspace. $S$ | unitary | $[\cdot,\cdot]_{\pm}$ | $Sb_{\boldsymbol{k}}S^{\dagger}$ | $Sg_{\boldsymbol{k}}S^{\dagger}$ | condition
---|---|---|---|---|---
$C_{2z}$ | ✓ | $+$ | $\eta_{1}\sigma_{0}b_{-\boldsymbol{k}}$ | $\eta_{1}\sigma_{0}g_{-\boldsymbol{k}}$ | —
$C_{3z}$ | ✓ | $+$ | $\eta_{0}\sigma_{0}b_{C_{3z}\boldsymbol{k}}$ | $\eta_{0}\sigma_{0}g_{C_{3z}\boldsymbol{k}}$ | —
$\sigma_{h}$ | ✓ | $+$ | $\eta_{0}\sigma_{0}b_{\boldsymbol{k}}$ | $-\eta_{0}\sigma_{0}g_{\boldsymbol{k}}$ | $D_{0}=0$
$\Theta$ | ✗ | $+$ | $\eta_{1}\sigma_{0}b_{-\boldsymbol{k}}$ | $\eta_{1}\sigma_{0}g_{-\boldsymbol{k}}$ | —
$\Theta_{s}$ | ✗ | $+$ | $s_{2}\eta_{1}\sigma_{0}b_{-\boldsymbol{k}}$ | $s_{2}\eta_{1}\sigma_{0}g_{-\boldsymbol{k}}$ | —
$\widetilde{\Theta}$ | ✗ | $+$ | $\eta_{2}\sigma_{0}b_{-\boldsymbol{k}}$ | $\eta_{2}\sigma_{0}g_{-\boldsymbol{k}}$ | —
$C$ | ✓ | $-$ | $\eta_{3}\sigma_{2}b_{\boldsymbol{k}}$ | $\eta_{3}\sigma_{2}g_{\boldsymbol{k}}$ | $w_{0}=0$
$P$ | ✓ | $-$ | $-i\eta_{3}\sigma_{2}b_{-\boldsymbol{k}}$ | $i\eta_{1}\sigma_{2}g_{-\boldsymbol{k}}$ |  $D_{0}=0$, $\boldsymbol{\rho}_{\theta\rightarrow 0}$
$C_{2z}\Theta$ | ✗ | $+$ | $\eta_{0}\sigma_{0}b_{\boldsymbol{k}}$ | $\eta_{0}\sigma_{0}g_{\boldsymbol{k}}$ | —
$C_{2z}P$ | ✓ | $-$ | $\eta_{2}\sigma_{2}b_{\boldsymbol{k}}$ | $i\sigma_{2}g_{\boldsymbol{k}}$ |  $D_{0}=0$, $\boldsymbol{\rho}_{\theta\rightarrow 0}$

Based on these symmetry representations, one finds that the dispersions are
constrained to have the form

$\displaystyle\begin{split}\epsilon_{(b,p),\eta}(\boldsymbol{k})&=p\,\bar{E}^{b}_{0}(\boldsymbol{k})+\eta\,\bar{E}^{b}_{1}(\boldsymbol{k})\,w_{0}\\\
&\quad+p\,\eta\,\bar{E}^{b}_{2}(\boldsymbol{k})\,D_{0}+\bar{E}^{b}_{3}(\boldsymbol{k})\,w_{0}D_{0}\end{split}$
(6)

and

$\epsilon_{(g,p),\eta}(\boldsymbol{k})=p\,\bar{E}^{g}_{0}(\boldsymbol{k};\eta)+\bar{E}^{g}_{1}(\boldsymbol{k};\eta)\,w_{0}D_{0}.$
(7)

To illustrate the changes of the form of the band structures when symmetries
are broken by $w_{0}\neq 0$ and $D_{0}\neq 0$, we introduced
$\bar{E}^{b}_{j}(\boldsymbol{k})$ and $\bar{E}^{g}_{j}(\boldsymbol{k};\eta)$
which are functions of $w_{0}$ and $D_{0}$ that stay finite when
$w_{0},D_{0}\rightarrow 0$.

As can be seen, the two terms in the first row of Eq. (6) just correspond to
the TBG band structure (with, $w_{0}=0$, or without, $w_{0}\neq 0$, chiral
symmetry). The second line shows that the hybridization with the graphene
bands for $D_{0}\neq 0$ distorts the TBG band structure in a way not present
in TBG, which is related to the $D_{0}$-induced breaking of $C_{2z}P$
symmetry. One of these additional terms leads to superexchange processes in
the TBG-like bands that are not present in TBG and will be discussed in Sec.
V.4. In Eq. (7), the first term, $\bar{E}^{g}_{0}(\boldsymbol{k};\eta)$,
simply becomes the graphene Dirac dispersion for $D_{0}\rightarrow 0$. We
further see that the property
$\epsilon_{(g,+),\eta}(\boldsymbol{k})=-\epsilon_{(g,-),\eta}(\boldsymbol{k})$
is only violated if $w_{0}$ and $D_{0}$ are simultaneously non-zero. This is
expected, as the graphene-like bands can only “feel” the broken chiral
symmetry if they hybridize with the TBG-like bands. For the parameters in Fig.
2(b), this is a rather weak effect around the Dirac cones.

###  II.4 Interactions and form factors

As our goal is to study the interacting behavior of MSTG, we next add
interaction terms to the Hamiltonian. We assume that these interactions in the
full continuum model are of the density-density form. Our full interaction
Hamiltonian thus reads as

$H_{1}^{\text{Full}}=\frac{1}{2N}\sum_{\boldsymbol{q}}V(\boldsymbol{q}):\rho^{\text{Full}}_{\boldsymbol{q}}\rho^{\text{Full}}_{-\boldsymbol{q}}:$
(8)

where $:\ldots:$ denotes normal ordering, $N$ is the number of moiré sites,
and $\rho^{\text{Full}}_{\boldsymbol{q}}$ is the density operator of momentum
$\boldsymbol{q}$ of the continuum electron operators $c_{\boldsymbol{r}}$ in
Eq. (1). Note $\boldsymbol{q}\in\mathbbm{R}$ in the sum in Eq. (8) which is
not restricted to the MBZ. For our analytical discussions below, we will not
have to specify the explicit form of $V(\boldsymbol{q})$ but only use
$V(\boldsymbol{q})>0$ and $V(\boldsymbol{q})=V(-\boldsymbol{q})$;
consequently, our analytical results will be valid regardless of the details
of the screening processes at high energies and of nearby gates and/or
substrates, that affect the detailed form $V(\boldsymbol{q})$. For our
numerics, we will use the single-gate-screened Coulomb potential:

$V(\boldsymbol{q})=\frac{1-e^{-2|\boldsymbol{q}|d_{s}}}{2A_{\text{moir\'{e}}}\epsilon_{0}\epsilon|\boldsymbol{q}|},$
(9)

where we have normalized by the real space area of the moiré unit cell
$A_{\text{moir\'{e}}}$. We vary the screening length $d_{s}$ in our numerical
calculations from $d_{s}=10$ to $80$ nm and find relatively little dependence
of the relative energies between phases and no qualitative change in the
ground states. We will therefore take $d_{s}=40$ nm as our default value. We
also vary the dielectric constant $\epsilon$ from $\epsilon=4$ to
$\epsilon=15$. We find the phase boundaries depend on $\epsilon$, though no
new phases emerge as $\epsilon$ is varied [the effect of varying interaction
parameters is further discussed in Appendix C.3, especially Fig. 9]. Unless
otherwise specified, we will take $\epsilon=7$.

We note that the above interaction term has all the symmetries of the
continuum model introduced in Sec. II, including the
$\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ symmetry. In the realistic system,
there is a Hund’s coupling, $J_{H}\neq 0$, which breaks
$\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ down to $\text{SU}(2)_{s}$. We
discuss the form of the Hund’s coupling in Appendix F.5 and Eq. (30), and will
note its consequences in Sec. III.3.

Neglecting normal ordering and replacing
$\rho^{\text{Full}}_{\boldsymbol{q}}\rightarrow\tilde{\rho}^{\text{Full}}_{\boldsymbol{q}}=\rho^{\text{Full}}_{\boldsymbol{q}}-\text{const.}\times\delta_{\boldsymbol{q},0}$
in Eq. (8) leads to a new form of the interaction,
$\widetilde{H}_{1}^{\text{Full}}$, which, however, can be rewritten as
$H_{1}^{\text{Full}}$ by a redefinition of the chemical potential (and energy
reference point). This is not the case anymore after projecting
$\rho^{\text{Full}}_{\boldsymbol{q}}$ and
$\tilde{\rho}^{\text{Full}}_{\boldsymbol{q}}$ in $H_{1}^{\text{Full}}$ and
$\widetilde{H}_{1}^{\text{Full}}$, respectively, to a finite set of bands of
$H_{0}$. As described in more detail in Appendix A.3, we follow Bultinck  _et
al._ (2019); Bernevig  _et al._ (2020) and rewrite the interaction before
projection such that it exhibits particle-hole symmetry with respect to the
charge neutrality point before and after projection. Denoting the electronic
creation operators for band $n$, of valley $\eta$, spin $s$, and with momentum
$\boldsymbol{k}\in\text{MBZ}$ by $f^{\dagger}_{\boldsymbol{k};n,\eta,s}$, the
projected interaction becomes

$\widetilde{H}_{1}=\frac{1}{2N}\sum_{\boldsymbol{q}}V(\boldsymbol{q})\delta\rho_{\boldsymbol{q}}\delta\rho_{-\boldsymbol{q}},$
(10a) with the symmetrized density operators
$\displaystyle\delta\rho_{\boldsymbol{q}}=\sum_{\boldsymbol{k}\in\text{MBZ}}\left[f^{\dagger}_{\text{MBZ}(\boldsymbol{k}+\boldsymbol{q})}F_{\boldsymbol{k},\boldsymbol{q}}f^{\phantom{\dagger}}_{\boldsymbol{k}}-\frac{1}{2}\sum_{\boldsymbol{G}\in\text{RL}}\delta_{\boldsymbol{q},\boldsymbol{G}}\text{tr}(F_{\boldsymbol{k},\boldsymbol{G}})\right].$
(10b)

Here
$\text{MBZ}(\boldsymbol{k}):=\boldsymbol{k}-\boldsymbol{G}_{\boldsymbol{k}}\in\text{MBZ}$
for the unique reciprocal lattice vector
$\boldsymbol{G}_{\boldsymbol{k}}\in\text{RL}$. We suppressed all indices of
$f_{\boldsymbol{k}}$ and $f^{\dagger}_{\boldsymbol{k}}$, which should be
viewed as column and row vectors in band, valley, and spin space and
introduced the matrix-valued form factors $F_{\boldsymbol{k},\boldsymbol{q}}$,
which contain all the microscopic details of the wavefunctions of the bands of
$H_{0}$ [see Eq. (71) for a formal definition]. Note that Eq. (10) holds for
any subset of bands that we want to keep. As already discussed in Sec. II.3
above, we will only keep the TBG-like bands and graphene-like bands (around
the K/K’ points) in the analytics, while we will allow for more bands in the
HF numerics.

We refer the interested reader to Appendix A.3, where a detailed discussion of
constraints on $F_{\boldsymbol{k},\boldsymbol{q}}$ resulting from symmetries,
Hermiticity, and the structure of $H_{0}$ can be found, and here only state a
few properties of the form factors that we will explicitly refer back to in
the main text. First, as a consequence of $\text{U}(1)_{v}$ and spin rotation
symmetry, the form factors can only have non-diagonal matrix structure in band
space,

$\left(F_{\boldsymbol{k},\boldsymbol{q}}\right)_{(n,\eta,s),(n^{\prime},\eta,s)}=\delta_{s,s^{\prime}}\delta_{\eta,\eta^{\prime}}\left(F^{\eta}_{\boldsymbol{k},\boldsymbol{q}}\right)_{n,n^{\prime}}.$
(11)

In accordance with our notation in Eq. (5), we will use the multi-index
notation $n=(t,p)$ where the “type” $t$ distinguishes between the graphene-
like, $t=g$, and TBG-like, $t=b$, bands and $p=+$ ($p=-$) labels the upper
(lower) band; for instance, for the electronic operators, it holds
$f_{\boldsymbol{k};(b,p),\eta,s}=b_{\boldsymbol{k};p,\eta,s}$ and
$f_{\boldsymbol{k};(g,p),\eta,s}=g_{\boldsymbol{k};p,\eta,s}$. As a
consequence of the $\sigma_{h}$ symmetry for $D_{0}=0$, the form factors
become block-diagonal in the sectors with different $\sigma_{h}$ eigenvalue.
As such, $F_{\boldsymbol{k},\boldsymbol{q}}^{tt^{\prime}}$ defined via

$\left(F^{tt^{\prime}}_{\boldsymbol{k},\boldsymbol{q}}\right)_{(p,\eta,s),(p^{\prime},\eta^{\prime},s^{\prime})}=\delta_{s,s^{\prime}}\delta_{\eta,\eta^{\prime}}\left(F^{\eta}_{\boldsymbol{k},\boldsymbol{q}}\right)_{(t,p)(t^{\prime},p^{\prime})}$
(12)

obeys

$F^{tt^{\prime}}_{\boldsymbol{k},\boldsymbol{q}}\propto\delta_{t,t^{\prime}},\quad\text{for
}D_{0}=0.$ (13)

Consequently, the entire interacting Hamiltonian $H_{0}+\widetilde{H}_{1}$
preserves the charge in the graphene and TBG system separately. This leads to
simplifications, which will be exploited in Sec. III below. For $D_{0}\neq 0$,
this is not the case anymore, as the form factors in Eq. (12) and, hence, the
density operators in the interaction will not be diagonal in the $t$ index and
scatter electrons between the two different types of bands.

In the limit $D_{0}=0$, one can also compute the form factors in the graphene
sector analytically [see Eq. (85) for the full expression]. Most importantly
for our purposes here, one finds that

$F^{gg}_{\boldsymbol{k},\boldsymbol{G}}=\delta_{\boldsymbol{G},0}\sigma_{0}\eta_{0}s_{0},\quad\boldsymbol{G}\in\text{RL},\quad\text{for
}D_{0}=0,$ (14)

where $\sigma_{0}$ are Pauli matrices in the “band space” (with indices
$p=\pm$). Here only the $\boldsymbol{G}=0$ component is finite, which is
related to the fact that the graphene bands do not “feel” the moiré
superlattice for $D_{0}=0$.

##  III Exact groundstates at $D_{0}=0$

We begin our discussion of the interacting physics in the decoupled limit,
defined as $D_{0}=0$, where all bands can be labeled by their mirror
eigenvalue $\sigma_{h}=\pm 1$, and the low-energy bands are those of TBG and
single-layer graphene. While one might intuitively expect that the presence of
the additional graphene Dirac cones, which have a much lower density of states
than the (almost) flat bands close to the magic angle of TBG, is not strong
enough to change the symmetry of the correlated insulating phase of TBG, it is
a priori not clear whether the density-density coupling between the two
subsystems can also induce the same symmetry-breaking order (and potentially
gap out) the graphene Dirac cones. It is further not clear whether exact
interacting ground states of MSTG can be identified in certain limits, similar
to TBG Lian  _et al._ (2020); Kang and Vafek (2019); Bultinck  _et al._
(2019). These aspects will be addressed in this section.

###  III.1 Hamiltonian and construction of eigenstates

To this end, let us focus on the flat limit, $W_{\text{TBG}}=0$, and postpone
the perturbative treatment of the finite TBG bandwidth to Sec. V.4. From our
discussion in Sec. II, we can read off that the low-energy Hamiltonian of MSTG
in the flat-decoupled limit is given by

$H_{\text{FD}}=H^{\text{g}}+H^{\text{b}}+\lambda H^{\text{gb}},$ (15)

consisting of three terms: $H^{\text{g}}$ is the Hamiltonian of (both valleys
of) single-layer graphene with Coulomb repulsion [see Eq. (90)],
$H^{\text{b}}$ the interacting Hamiltonian of TBG in the flat limit, given by

$H^{\text{b}}=\frac{1}{2N}\sum_{\boldsymbol{q}}V(\boldsymbol{q})\delta\rho^{\text{b}}_{\boldsymbol{q}}\delta\rho^{\text{b}}_{-\boldsymbol{q}}$
(16)

where we defined the projected subsystem density operators

$\delta\rho^{t}_{\boldsymbol{q}}=\sum_{\boldsymbol{k}\in\text{MBZ}}\left[t^{\dagger}_{\text{MBZ}(\boldsymbol{k}+\boldsymbol{q})}\left(F^{tt}_{\boldsymbol{k},\boldsymbol{q}}\right)t^{\phantom{\dagger}}_{\boldsymbol{k}}-\frac{1}{2}\sum_{\boldsymbol{G}\in\text{RL}}\delta_{\boldsymbol{q},\boldsymbol{G}}\text{tr}(F^{tt}_{\boldsymbol{k},\boldsymbol{G}})\right]$
(17)

where $t=g,b$. Finally, the last term in Eq. (15) describes the coupling
between the two subsystems via a density-density interaction,

$\lambda
H^{\text{gb}}=\frac{\lambda}{N}\sum_{\boldsymbol{q}}V(\boldsymbol{q})\delta\rho^{\text{g}}_{-\boldsymbol{q}}\delta\rho^{\text{b}}_{\boldsymbol{q}},$
(18)

where the additional prefactor $\lambda\in[0,1]$ has been introduced to
adiabatically turn on the $H^{\text{gb}}$ interaction (the physical system
corresponds to $\lambda=1$).

For $\lambda=0$, the Hamiltonian is just the sum $H^{\text{g}}+H^{\text{b}}$
of the two commuting subsystem Hamiltonians. So its eigenstates are just given
by (all combinations) of the individual eigenstates of graphene and TBG with
the correct particle number: let us fix a certain integer filling
$\nu_{\text{b}}=0,\pm 1,\pm 2,\pm 3,\pm 4$ of the TBG system. Then the
graphene system will be at a corresponding filling, which we write formally as
$\nu_{g}=\nu_{g}(\nu_{\text{b}})$ and its ground state will be a semimetal
with ($\nu_{g}=0$) or without ($\nu_{g}\neq 0$) doping, that exhibits
correlations but does not spontaneously break any symmetries; these properties
are well established (theoretically and experimentally) for graphene Castro
Neto  _et al._ (2009). Let us denote the ground-state of the graphene system
at filling $\nu_{g}=\nu_{g}(\nu_{\text{b}})$ by
$\ket{\Psi^{\text{g}}_{0}(\nu_{\text{b}})}$ and its (gapless) excited states
with the same particle number by $\ket{\Psi^{\text{g}}_{j}(\nu_{\text{b}})}$,
$j>0$. In the hypothetical absence of any correlations in the graphene
subspace, $\ket{\Psi^{\text{g}}_{0}(\nu_{\text{b}})}$
($\ket{\Psi^{\text{g}}_{j}(\nu_{\text{b}})}$) would just be Slater-determinant
state(s) with the Dirac cones filled up to to the chemical potential (and some
additional particle-hole excitations).

Exact ground states of the flat-band TBG Hamiltonian in Eq. (16) have been
discussed previously Lian  _et al._ (2020), which we will very briefly review
here using our notation, in order to set the stage for the extension to MSTG.
Upon defining new operators according to [referred to as “chiral basis” in
Bernevig  _et al._ (2020); in this basis the chiral form factors, see Eq.
(93), are diagonal]

$\tilde{b}^{\phantom{\dagger}}_{\boldsymbol{k};c,\eta,s}=U_{c,p}b^{\phantom{\dagger}}_{\boldsymbol{k};p,\eta,s},\quad
U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-i\\\ 1&i\end{pmatrix},$ (19)

consider the set of states

$\ket{\Psi_{0}^{\text{b}}(\nu_{b})}=\prod_{\boldsymbol{k}\in\text{MBZ}}\prod_{c=\pm}\prod_{j_{c}=1}^{\nu_{c}}\tilde{b}^{\dagger}_{\boldsymbol{k};c,\eta^{c}_{j_{c}},s^{c}_{j_{c}}}\ket{0^{\text{b}}}$
(20)

with arbitrary combinations of occupied flavors
$\\{\eta^{\pm}_{j},s^{\pm}_{j}\\}$ such that
$\nu_{+}+\nu_{-}=4+\nu_{\text{b}}$. It was shown in Lian  _et al._ (2020) that

$\delta\rho^{\text{b}}_{\boldsymbol{q}}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}=\sum_{\boldsymbol{G}\in\text{RL}}\delta_{\boldsymbol{q},\boldsymbol{G}}R_{\boldsymbol{G}}\ket{\Psi_{0}^{\text{b}}(\nu_{b})},$
(21)

with
$R_{\boldsymbol{G}}=\nu_{b}\sum_{\boldsymbol{k}}\text{tr}[F^{bb}(\boldsymbol{k},\boldsymbol{G})]/8$
when $w_{0}=0$, i.e., in the chiral limit. As such, all of these states are
exact eigenstates of $H^{\text{b}}$ in Eq. (16). Ref. Lian  _et al._ , 2020
further showed that these states will always be ground states of
$H^{\text{b}}$ for $\nu_{b}=0$; the same holds for all other integer $\nu_{b}$
as long as the flat-metric condition,

$F^{bb}_{\boldsymbol{k},\boldsymbol{G}}=\mathbbm{1}f(\boldsymbol{G}),\qquad\forall\,\boldsymbol{k},\boldsymbol{G},$
(22)

is not violated by a significantly large amount. Furthermore, when turning on
$w_{0}\neq 0$, the subset of states in Eq. (20) with
$\eta_{j}^{+}=\eta_{j}^{-}$ and $s_{j}^{+}=s_{j}^{-}$, which are necessarily
at even integer $\nu_{b}$, still obey Eq. (21) and remain ground states of
$H_{b}$ [unless $\nu_{b}=\pm 2$ and Eq. (22) is sufficiently violated].

Having established the spectrum of $H_{\text{FD}}$ in Eq. (15) for
$\lambda=0$, let us next discuss what happens once $\lambda$ is turned on.
Using the fact that the graphene form factors obey Eq. (14), we show in
Appendix B.1 that $H^{\text{gb}}$ can be rewritten in the low-energy spectrum
of MSTG as

$H^{\text{gb}}=\frac{1}{N}\sum_{\boldsymbol{q}}V(\boldsymbol{q})\delta\rho^{\text{g}}_{-\boldsymbol{q}}\left(\delta\rho^{\text{b}}_{\boldsymbol{q}}-\sum_{\boldsymbol{G}\in\text{RL}}\delta_{\boldsymbol{q},\boldsymbol{G}}R_{\boldsymbol{G}}\right)+E_{0}(\nu_{b}),$
(23)

where $E_{0}(\nu_{b})$ is just a constant energy. So we immediately see that
the property (21) of all of the exact TBG states
$\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ defined above, implies

$H^{\text{gb}}\ket{\Psi^{\text{g}}_{j}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}=E_{0}(\nu_{b})\ket{\Psi^{\text{g}}_{j}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}.$
(24)

Consequently, all of the states
$\ket{\Psi^{\text{g}}_{j}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ remain
exact eigenstates of the full MSTG Hamiltonian $H_{\text{FD}}$ in the flat-
decoupled limit, at arbitrary $\lambda$.

Whether the states
$\ket{\Psi^{\text{g}}_{0}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ will
also remain the exact ground states is a more subtle question: since the
states $\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ break symmetries,
$H_{\text{FD}}|_{\lambda=0}$ will have a gapless Goldstone spectrum. In
principle, an arbitrarily small $\lambda$ could lower the energy of some of
those states below that of
$\ket{\Psi^{\text{g}}_{0}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$.
However, we show in Appendix B.2 that is not the case if $\nu_{b}=0$ or
$w_{0}=0$ or Eq. (22) holds. Therefore, a finite $\lambda>0$ is required
before $\ket{\Psi^{\text{g}}_{0}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$
cease to be the exact ground states.

###  III.2 Discussion of ground states

Taken together, we have shown that the states
$\ket{\Psi^{\text{g}}_{j}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ where
$\ket{\Psi^{\text{g}}_{j}(\nu_{b})}$ is just the spectrum of single-layer
graphene at filling $\nu_{g}=\nu_{g}(\nu_{\text{b}})$ and
$\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ is any of the states in Eq. (20) are
exact eigenstates of the MSTG Hamiltonian in the chiral-flat-decoupled limit
($w_{0}=W_{\text{TBG}}=D_{0}=0$), $H_{\text{FD}}|_{w_{0}=0}$, for any integer
$\nu_{b}$. Furthermore, there is a finite region of $\lambda$ for which
$\ket{\Psi^{\text{g}}_{0}(\nu_{b})}\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ will
remain a groundstate of MSTG if $\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ is a
groundstate of TBG (recall that $\ket{\Psi_{0}^{\text{b}}(\nu_{b})}$ is
guaranteed to be a groundstate of TBG for $\nu_{b}=0$ without further
assumptions while it requires that the flat-metric condition is not too
strongly violated for $\nu_{b}\neq 0$). Finally, away from the chiral limit
$w_{0}\neq 0$, the subset of states in Eq. (20) with
$\eta_{j}^{+}=\eta_{j}^{-}$ and $s_{j}^{+}=s_{j}^{-}$ are known to be ground
states of TBG in the flat limit for integer $\nu_{b}$ and if the flat metric
condition holds Lian  _et al._ (2020). Our analysis shows that these states
remain exact eigenstates for $\lambda\neq 0$ and also ground states for
$|\lambda|<\lambda_{c}>0$ in the non-chiral-flat-decoupled limit
($W_{\text{TBG}}=D_{0}=0$, $w_{0}\neq 0$).

In all of these limits, we see that the graphene subsystem retains its
(correlated but symmetry-unbroken and, depending on $\nu_{b}$, doped)
semimetallic properties for all integer filling fractions $\nu_{b}$. This is
consistent with experiment, where quantum oscillations indicate a dispersive
Dirac cone at $D_{0}=0$ Park  _et al._ (2021); Hao  _et al._ (2021).
Furthermore, the exact eigenstates established above will be used as our
starting point for further analytical considerations in Sec. V and their
product-state nature motivates our HF numerical study of the problem in Sec.
IV. Both numerics and analytics will complement the discussion presented above
by (i) validating the stability of the Dirac cones at $D_{0}=0$ in schemes
that do not rely on $\lambda$ being small and (ii) by tuning away from the
exactly solvable limits ($W_{\text{TBG}}=D_{0}\neq 0$) and (iii), for the
numerics, including additional remote bands.

Table 2: We list the different candidate phases in the TBG-like subspace, constructed as the discrete set of states that are part of the large manifold of exact ground states in the chiral-flat-decoupled limit, see Sec. III, but transform under the irreducible representations of the symmetries of the real system. Here $\boldsymbol{1}$ ($\boldsymbol{3}$) is the singlet (triplet) representation of $\text{SU}(2)_{s}$ and $0$ ($1$) the one(two)-dimensional representation of $\text{U}(1)_{v}$. For future reference in Sec. VII, we list the behavior ($\pm$ denoting even/odd, and ✗ indicating absence) under both spinful ($\Theta_{s}$), and valley ($\widetilde{\Theta}$) time-reversal symmetry, see Table 1. The last two columns indicate which states are Hund’s partners Scheurer and Samajdar (2020), i.e., transform into each other when reversing the sign of the Hund’s coupling $J_{H}$ while being exactly degenerate in the $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$-symmetric limit, and which sign of $J_{H}$ favors the respective state. Type | Short form | $Q^{b}$ | $\widetilde{Q}^{b}$ | SU(2)s | U(1)v | $C_{2z}$ | $\Theta_{s}/\widetilde{\Theta}$ | Hund’s part. | $J_{H}$
---|---|---|---|---|---|---|---|---|---
spin polarized | SP | $\sigma_{0}\eta_{0}\boldsymbol{s}$ | $\widetilde{\sigma}_{0}\eta_{0}\boldsymbol{s}$ | $\boldsymbol{3}$ | $0$ | ✓ |  $-$/✗ | SVP | $<0$
valley polarized | VP | $\sigma_{0}\eta_{z}s_{0}$ | $\widetilde{\sigma}_{0}\eta_{z}s_{0}$ | $\boldsymbol{1}$ | $0$ | ✗ |  $-$/$-$ | — | $0$
$\Theta$-even IVC | IVC+ | $\sigma_{0}\eta_{x,y}s_{0}$ | $\widetilde{\sigma}_{0}\eta_{x,y}s_{0}$ | $\boldsymbol{1}$ | $1$ | ✓ |  $+$/$-$ | SIVC+ | $>0$
$\Theta$-odd IVC | IVC- | $\sigma_{y}\eta_{x,y}s_{0}$ | $\widetilde{\sigma}_{z}\eta_{x,y}s_{0}$ | $\boldsymbol{1}$ | $1$ | ✓ |  $-$/$+$ | SIVC- | $>0$
$\Theta$-odd, sublattice pol./Hall | SLP- | $\sigma_{y}\eta_{0}s_{0}$ | $\widetilde{\sigma}_{z}\eta_{0}s_{0}$ | $\boldsymbol{1}$ | $0$ | ✓ |  $-$/$-$ | — | $0$
$\Theta$-even, sublattice pol./valley Hall | SLP+ | $\sigma_{y}\eta_{z}s_{0}$ | $\widetilde{\sigma}_{z}\eta_{z}s_{0}$ | $\boldsymbol{1}$ | $0$ | ✗ |  $+$/$+$ | — | $0$
$\Theta$-odd, spin-sublattice-pol./spin Hall | SSLP- | $\sigma_{y}\eta_{0}\boldsymbol{s}$ | $\widetilde{\sigma}_{z}\eta_{0}\boldsymbol{s}$ | $\boldsymbol{3}$ | $0$ | ✓ |  $+$/✗ | SSLP+ | $<0$
spin-valley polarized | SVP | $\sigma_{0}\eta_{z}\boldsymbol{s}$ | $\widetilde{\sigma}_{0}\eta_{z}\boldsymbol{s}$ | $\boldsymbol{3}$ | $0$ | ✗ |  $+$/✗ | SP | $>0$
$\Theta$-even, spin-pol. IVC | SIVC+ | $\sigma_{0}\eta_{x,y}\boldsymbol{s}$ | $\widetilde{\sigma}_{0}\eta_{x,y}\boldsymbol{s}$ | $\boldsymbol{3}$ | $1$ | ✓ |  $-$/✗ | IVC+ | $<0$
$\Theta$-odd, spin-pol. IVC | SIVC- | $\sigma_{y}\eta_{x,y}\boldsymbol{s}$ | $\widetilde{\sigma}_{z}\eta_{x,y}\boldsymbol{s}$ | $\boldsymbol{3}$ | $1$ | ✓ |  $+$/✗ | IVC- | $<0$
$\Theta$-even, spin-subl. pol./spin-valley Hall | SSLP+ | $\sigma_{y}\eta_{z}\boldsymbol{s}$ | $\widetilde{\sigma}_{z}\eta_{z}\boldsymbol{s}$ | $\boldsymbol{3}$ | $0$ | ✗ |  $-$/✗ | SSLP- | $>0$

###  III.3 Resultant candidate states

To build the foundation for these additional analytical and numerical
computations, we will use the exact (and highly degenerate) ground states
established above in the chiral-flat-decoupled limit to construct a finite set
of candidate phases and their respective order parameters.

We first define the correlation matrix $P_{\boldsymbol{k}}$ with elements

$\left(P_{\boldsymbol{k}}\right)_{(t,p,\eta,s),(t^{\prime},p^{\prime},\eta^{\prime},s^{\prime})}:=\braket{\Psi_{0}}{f^{\dagger}_{\boldsymbol{k};(t,p),\eta,s}f^{\phantom{\dagger}}_{\boldsymbol{k};(t^{\prime},p^{\prime}),\eta^{\prime},s^{\prime}}}{\Psi_{0}}$
(25)

to characterize a given ground state $\ket{\Psi_{0}}$ of MSTG. Hermiticity
implies $P^{\dagger}_{\boldsymbol{k}}=P^{\phantom{\dagger}}_{\boldsymbol{k}}$.
As is common, we further write
$P_{\boldsymbol{k}}=\frac{1}{2}(\mathbbm{1}+Q_{\boldsymbol{k}})$ and will use
$Q_{\boldsymbol{k}}$ as our “order parameter” to characterize the (potentially
symmetry-broken) structure of $\ket{\Psi_{0}}$. It must obey

$Q_{\boldsymbol{k}}^{\dagger}=Q^{\phantom{\dagger}}_{\boldsymbol{k}},\quad\frac{1}{N}\sum_{\boldsymbol{k}}\text{tr}\left[Q_{\boldsymbol{k}}\right]=\nu_{b}+\nu_{g}\equiv\nu.$
(26)

As we have seen above, the ground states of MSTG in the flat-decoupled limit
($W_{\text{TBG}}=D_{0}=0$) obey

$\left(Q_{\boldsymbol{k}}\right)_{(t,p,\eta,s),(t^{\prime},p^{\prime},\eta^{\prime},s^{\prime})}=\delta_{t,t^{\prime}}\left(Q^{t}_{\boldsymbol{k}}\right)_{(p,\eta,s),(p^{\prime},\eta^{\prime},s^{\prime})},$
(27)

i.e., do not exhibit any “coherence” between the graphene and TBG sectors.
This is expected as the presence of $\sigma_{h}$ requires any order parameter
to be either even (diagonal in $t$ space) or odd (off-diagonal) under
$\sigma_{h}$; due to the large density of states in the TBG sector, we expect
the former to dominate. Once $D_{0}\neq 0$, mixing is allowed, as we will see
in our numerics below and discuss in detail analytically in Sec. V.2.

Furthermore, the analysis above reveals that the ground state in the TBG
sector will be of the form of Eq. (20). For instance, for $\nu_{b}=0$ with
$\eta_{j}^{+}=\eta_{j}^{-}=(-1)^{j}$ and $s_{j}^{+}=s_{j}^{-}=\,\uparrow$,
$j=1,2$, it holds
$\widetilde{Q}^{b}_{\boldsymbol{k}}=\widetilde{\sigma}_{0}\eta_{0}s_{z}$ where
$\widetilde{Q}^{b}_{\boldsymbol{k}}=U^{*}Q^{b}_{\boldsymbol{k}}U^{T}$ is the
order parameter $Q^{b}_{\boldsymbol{k}}$ in Eq. (27) in the TBG subspace
transformed to the chiral basis of Eq. (19); in the basis of Eq. (25), it
holds $Q^{b}_{\boldsymbol{k}}=\sigma_{0}\eta_{0}s_{z}$, which we will refer to
as spin polarized (SP) state. Here and in the following we will use
$\sigma_{j}$ ($\widetilde{\sigma}_{j}$) to denote Pauli matrices in the band-
space with index $p$ (in the chiral basis with index $c$). Besides the SP
state, Eq. (20) describes many other possible ground states, that are exactly
degenerate in the chiral-flat-decoupled limit. A systematic way of seeing this
proceeds by noting that the $\text{U}(4)\times\text{U}(4)$ identified Bultinck
_et al._ (2019); Kang and Vafek (2019); Bernevig  _et al._ (2020) for TBG also
persists as a symmetry of MSTG in chiral-flat-decoupled limit Cǎlugǎru  _et
al._ (2021); this immediately follows from the structure of $H_{\text{FD}}$ in
Eq. (15). We will here refer to this symmetry group as
$(\text{U}(4)\times\text{U}(4))_{\text{b,cf}}$ and its action is particularly
simple in the chiral basis Bernevig  _et al._ (2020),

$\widetilde{b}_{\boldsymbol{k}}\,\rightarrow\,\mathcal{U}\widetilde{b}_{\boldsymbol{k}},\quad\mathcal{U}=e^{i\sum_{j=0,3}\sum_{\mu,\mu^{\prime}=0}^{3}\varphi_{j,\mu,\mu^{\prime}}\widetilde{\sigma}_{j}\eta_{\mu}s_{\mu^{\prime}}}.$
(28)

The form of these transformations is readily inferred from Eq. (94) which
indicates that all $\mathcal{U}$ with $[\mathcal{U},\widetilde{\sigma}_{3}]=0$
will leave $H_{\text{FD}}$ invariant. Under Eq. (28), the order parameter
defined above transforms as
$\widetilde{Q}^{b}_{\boldsymbol{k}}\rightarrow\mathcal{U}^{*}\widetilde{Q}^{b}_{\boldsymbol{k}}\mathcal{U}^{T}$,
which allows us to generate the entire (continuous) set of exactly degenerate
ground states from one “seed” state, such as the SP state,
$\widetilde{Q}^{b}_{\boldsymbol{k}}=\widetilde{\sigma}_{0}\eta_{0}s_{z}$.
Since this seed state and $\mathcal{U}$ commute with $\widetilde{\sigma}_{3}$,
we know that $[\widetilde{Q}^{b}_{\boldsymbol{k}},\widetilde{\sigma}_{3}]=0$
for all ground states. Further noting that
$(\widetilde{Q}^{b}_{\boldsymbol{k}})^{2}=\mathbbm{1}$ (physically related to
the Slater-determinant nature), we can, thus, summarize the properties in the
original basis as

$Q^{b}_{\boldsymbol{k}}=Q^{b},\quad[Q^{b},\sigma_{2}]=0,\quad(Q^{b})^{2}=\mathbbm{1},\quad\text{tr}\left[Q^{b}\right]=0,$
(29)

at charge neutrality.

The actual Hamiltonian of MSTG is not in the chiral-flat-decoupled limit and
does not exhibit an exact $(\text{U}(4)\times\text{U}(4))_{\text{b,cf}}$
symmetry. Intuitively, this can be thought of as generating an easy axis in
this multi-dimensional space of degenerate states, favoring a specific
(subspace of) state(s) in Eq. (29). While energetics is required to decide
which phase is ultimately preferred by the system—the aim of the subsequent
sections—we can use symmetries to derive the discrete and finite set of
possible “candidate states”: to this end, we impose only $\text{U}(1)_{v}$ and
global spin rotations, $\text{SU}(2)_{s}$, as exact continuous symmetries. We
then know that the candidate order parameters must transform under the
irreducible representations of these symmetry groups (and be even or odd under
the exact discrete symmetries $C_{2z}$ and $\Theta$), leading to the $11$
options listed in Table 2. In order to connect smoothly to the limit
$D_{0}\rightarrow 0$, we take here $Q^{g}_{\boldsymbol{k}}=-\sigma_{z}$ in the
graphene subspace but emphasize that these order parameters are only used to
define the different states and characterize their symmetries; for our
numerical and analytical discussion below, they are only taken to be the
starting point and we will allow for (and also find) mixing between the TBG-
like and graphene-like sectors when $D_{0}\neq 0$ as well as momentum
dependence in $Q_{\boldsymbol{k}}$.

Finally, we point out that the model introduced in Sec. II, and which we study
energetics in below, has an exact $\text{SU}(2)_{+}\times\text{SU}(2)_{-}$
symmetry. Therefore, certain pairs of states, which we call Hund’s partners
following Scheurer and Samajdar (2020), have to be exactly degenerate, see
Table 2. As we noted earlier, in the realistic system, there is a non-zero
intervalley Hund’s coupling, $J_{H}\neq 0$, which will break
$\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ down to $\text{SU}(2)_{s}$, albeit
weakly, and favor one member of each of the pairs over the other. For most of
the following study we will focus on the
$\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ limit and, hence, can, without loss
of generality, restrict the discussion to the first $7$ states above the line
in Table 2. However, one has to keep in mind that the real system will realize
only one state of each Hund’s pair, which will depend on the (unknown) sign of
$J_{H}$ and precise form of the Hund’s coupling [see Appendix F.5 for more
details]. Taking, for concreteness, the intervalley Hund’s coupling to be of
the form

$H_{2}=\frac{J_{H}}{N}\sum_{\boldsymbol{q}}\boldsymbol{S}_{\boldsymbol{q}}^{+}\cdot\boldsymbol{S}_{\boldsymbol{q}}^{-},\quad\boldsymbol{S}_{\boldsymbol{q}}^{\pm}=\frac{1}{2}c^{\dagger}_{\boldsymbol{k}+\boldsymbol{q}}(\mathbbm{1}\pm\eta_{z})c_{\boldsymbol{k}},$
(30)

we can, in the chiral-flat-decoupled limit, uniquely associate a single state
of each Hund’s pair with a given sign of $J_{H}$; this is indicated in the
last column in Table 2.

##  IV Hartree-Fock numerics

As it facilitates the presentation of the results, we will begin the
discussion of correlated phases away from the chiral-flat-decoupled limit with
the HF numerics and postpone the complementary analytics to Sec. V.
Furthermore, we will first focus on the charge-neutrality point, $\nu=0$.

###  IV.1 Hartree-Fock Method

In the HF approximation, one focuses on Slater-determinant states
$\ket{\Psi[P_{\boldsymbol{k}}]}$ characterized by the correlation matrix
$P_{\boldsymbol{k}}$ as defined in Eq. (25). Consequently, it holds
$P^{2}_{\boldsymbol{k}}=P^{\phantom{2}}_{\boldsymbol{k}}$ or, equivalently,
$Q_{\boldsymbol{k}}^{2}=\mathbbm{1}$; this is also true for the exact
candidate ground states constructed in Sec. III.3 for the chiral-flat-
decoupled limit and, hence, the HF approximation is expected to provide
reliable results.

The goal of our HF numerics will be to determine the optimal
$Q_{\boldsymbol{k}}$ that yields the lowest energy expectation value with
respect to the interacting Hamiltonian for MSTG introduced in Sec. II. To be
more specific, we start from the full Hamiltonian,
$H=H^{\text{Full}}_{0,1}+H^{\text{Full}}_{0,2}+H_{1}$, consisting of the
continuum model in Eq. (1), the displacement-field term in Eq. (4),
supplemented by the density-density interaction in Eq. (8), and perform a
mean-field decoupling. Using the same notation as in Eq. (10b), the resulting
HF mean-field Hamiltonian reads as

$\begin{split}H^{\text{MF}}&=\sum_{\boldsymbol{k}\in\text{MBZ}}\epsilon_{n,\eta}(\boldsymbol{k})f^{\dagger}_{\boldsymbol{k};n,\eta,s}f^{\phantom{\dagger}}_{\boldsymbol{k};n,\eta,s}\\\
&+\sum_{\boldsymbol{k}\in\text{MBZ}}f^{\dagger}_{\boldsymbol{k}}\left[h_{\text{H}}[P](\boldsymbol{k})+h_{\text{F}}[P](\boldsymbol{k})\right]f_{\boldsymbol{k}}\\\
&-\frac{1}{2}\sum_{\boldsymbol{k}\in\text{MBZ}}\text{Tr}\left[h_{\text{H}}[P](\boldsymbol{k})P_{\boldsymbol{k}}^{T}+h_{\text{F}}[P](\boldsymbol{k})P_{\boldsymbol{k}}^{T}\right],\end{split}$
(31)

where the Hartree and Fock contributions to the mean-field Hamiltonian can be
written in terms of the projector $P_{\boldsymbol{k}}$ as

$h_{\text{H}}[P](\boldsymbol{k})=\frac{1}{N}\sum_{\boldsymbol{G}\in\text{RL}}V(\boldsymbol{G})F_{\boldsymbol{k},\boldsymbol{G}}\sum_{\boldsymbol{k}^{\prime}\in\text{MBZ}}\text{Tr}\left[F^{*}_{\boldsymbol{k}^{\prime},\boldsymbol{G}}P_{\boldsymbol{k}^{\prime}}\right]$
(32)

and

$h_{\text{F}}[P](\boldsymbol{k})=-\frac{1}{N}\sum_{\boldsymbol{q}}V(\boldsymbol{q})F^{\dagger}_{\boldsymbol{k},\boldsymbol{q}}P^{T}_{\boldsymbol{k}+\boldsymbol{q}}F_{\boldsymbol{k},\boldsymbol{q}},$
(33)

respectively. This form of the HF mean-field Hamiltonian is valid for an
arbitrary number of bands kept. In the numerics presented here, we will focus
on the four bands for each spin and valley flavor that are closest to the
Fermi level, which contains the graphene-like and TBG-like bands we focus on
in the analytics. We verify for representative values of $D_{0}$ and $w_{0}$
that the solutions we obtain are stable against doubling the number of remote
bands in our self-consistent calculation in Appendix C.3.

As pointed out in several HF works on TBG Bultinck  _et al._ (2019); Liu  _et
al._ (2021b); Xie and MacDonald (2020); Liao  _et al._ (2021), it is important
to note that the continuum model $H_{0,1}^{\text{Full}}$ already references
electron-electron interactions in the experimentally determined values for
microscopic model parameters and therefore we must define a reference
subtraction projector $P_{0}$ such that interactions will not be double
counted in our numerics. We here choose $P_{0}$ such that the projected low-
energy Hamiltonian exhibits the manifestly particle-hole symmetric interaction
in Eq. (10). As shown in Bernevig  _et al._ (2020) for TBG, this ansatz has
the natural interpretation of effectively taking into account the HF
contributions from all remote bands that have been projected out.

Figure 3: HF band structures in the limit $D_{0}=0$ (left) and nonzero $D_{0}$
(right) as obtained from Eq. (31) for the optimal $P_{\boldsymbol{k}}$ in the
respective symmetry-breaking channel defined in Table 2. We show only the
states which are either ground states or subleading energy states for some
region of our phase diagram, with a full set of band structures available in
Appendix C.3. Each band of the HF Hamiltonian is colored according to the
expectation value of the mirror symmetry operator within each band as a
function of $\boldsymbol{k}$, with the expectation value ranging from -1 (red)
to +1 (blue).

To determine the optimal $P_{\boldsymbol{k}}$, we start with an initial guess
for it with the symmetries of a given candidate order in Table 2. We then use
the HF Hamiltonian in Eq. (31) to compute a new projector $P_{\boldsymbol{k}}$
and iterate until $P_{\boldsymbol{k}}$ converges. More details on our
iterative HF procedure and subtraction point are given in Appendix C.2.

###  IV.2 Band Structures

In this section, we will discuss the band structures of the self-consistent
solutions we find at $D_{0}=0$ for each of the states in Table 2 and how these
band structures evolve as $D_{0}$ increases.

####  IV.2.1 Self-Consistent Band Structures at $D_{0}=0$

At $D_{0}=0$, we can separately describe the behavior of the graphene-like and
TBG-like bands for each type of state, since, for all states we consider, the
graphene-like bands near the Fermi-level at $D_{0}=0$ do not mix with any
other bands near the K and K′ points. They have a bandwidth larger than the
scale of Coulomb interactions. For this reason, the Dirac cones of the
graphene-like bands prefer a $Q^{g}_{\boldsymbol{k}}$ which equally fills the
lower bands of the continuum model and preserves all point group symmetries
for every class of solution we study. The Dirac cones thus remain semimetallic
and are, in this sense, “spectators” at $D_{0}=0$, in agreement with Sec. III.
On the other hand, the TBG-like bands have a bandwidth (5-10 meV) smaller than
the scale of the Coulomb energy at $D_{0}=0$ and therefore become insulating
as they are polarized for a given symmetry-breaking $P_{\boldsymbol{k}}$.

We find converged solutions for each ansatz in Table 2 and show representative
band structures for those states which have the lowest energy in the leftmost
panels of Fig. 3. Additional band structures for solutions not shown in Fig. 3
can be found in Appendix C.3. We note the similarity of the IVC- band
structure in the TBG-like bands to the band structure of the ground state in
Ref. Bultinck  _et al._ , 2019.

Figure 4: Phases obtained by self-consistent HF method as a function of
$w_{0}$ and $D_{0}$. Both the leading instability (red) and the subleading
ones (black, in paraenthesis) are shown. The phases are identified by the
symmetries they break, see Table 2, and representative band structures can be
found in Fig. 3. We show only one state as a representative of each pair of
Hund’s partners which are degenerate in our Hartree-Fock procedure. We label
the nearly degenerate SLP+, SLP-, SSLP+, and SSLP- phases with "SLP", though
note the preference in our numerics is for the SLP- state, from a Hartree
contribution.

####  IV.2.2 Self-Consistent Band Structures for $D_{0}>0$

All the solutions described for $D_{0}=0$ are insulating in the TBG-like bands
and semimetallic in the graphene-like bands. However, as $D_{0}$ increases and
the TBG-like and graphene-like bands begin to hybridize, the graphene-like
bands begin to play a more important role.

For the spin and valley polarized states (SP and VP) which preserve
$C_{2z}\Theta$ and U(1)v, the hybridized Dirac crossings are protected (and
pinned to the K/K′ points due to $C_{3z}$), meaning if the TBG-like bands
acquire a spin or valley polarization, the Dirac crossings of the graphene-
like bands must connect to these polarized bands as they are pushed away from
the Fermi level. Away from the K points of the MBZ, the graphene-like bands
will likely still prefer to fill the lowest bands of the non-interacting
model. Therefore, the VP and SP are generically expected to be metallic for
$D_{0}>0$. This is indeed what we find, as shown in the right panels for the
SP and VP states in Fig. 3.

The band structures of the IVC± states also retain Dirac crossings at the K/K′
points for nonzero $D_{0}$, as they exhibit $C_{3z}$ and a
$\boldsymbol{k}$-local anti-unitary symmetry that commutes with $C_{3z}$ (for
the IVC+ and IVC- these are $C_{2z}\Theta$ and $C_{2z}\widetilde{\Theta}$,
respectively). However, unlike the SP and VP states, the IVC± states also
preserve SU(2)s and $C_{2z}$, which pin the Dirac crossings at the Fermi level
at $\nu=0$. We therefore expect the intervalley coherent states will remain
semimetallic as $D_{0}$ increases. We observe this to be true for the self
consistent solutions, as can be seen in the right IVC- panel in Fig. 3.

The last class of states are the sublattice-polarized,
$C_{2z}\Theta$-symmetry-breaking states (SLP±, SSLP-) which preserve U(1)v. We
expect these states will generally be insulating for nonzero $D_{0}$ as there
are no protected Dirac crossings and both the TBG-like and graphene-like bands
can be gapped out. This is indeed seen in our numerics, with insulating band
structures for the SLP+, SLP-, and SSLP- states. The band structure of the
SLP- state is shown in Fig. 3.

###  IV.3 Energies and phase diagram

Having established the band structures of the different possible phases, we
next turn to their relative energetics and discuss which states are expected
to be favored energetically.

The evolution of the energies of each of our self-consistent solutions as a
function of $w_{0}$ and $D_{0}$ is shown in Figs. 5 and 6. As mentioned
before, at $D_{0}=0$, the Hamiltonian of the system is given by the sum of the
Hamiltonian of TBG and that of graphene, both with Coulomb interactions, which
are further coupled to each other by a density-density interaction. While we
only constructed exact eigenstates in Sec. III for the chiral-flat-decoupled
limit, we expect a similar picture when $w_{0},W_{\text{TBG}}\neq 0$: given
the bandwidth of the graphene-like bands is large compared to the scale of the
Coulomb interactions, we expect the graphene bands will prefer to fill the
lower bands of the continuum model. As the graphene density of states is small
compared to that of the flat bands of TBG, it should not crucially alter the
ground state in the TBG sector—at least close to the magic angle. Based on
previous work Bultinck  _et al._ (2019); Liao  _et al._ (2021), we thus expect
that the IVC- state has the lowest energy for $D_{0}=0$ (though with a smaller
energy difference than in previous works between our IVC- and spin polarized
phase due to our choice of subtraction point). Both expectations for the
graphene-like and TBG-like bands are confirmed by our numerics which finds the
IVC- state has the lowest energy of all our candidates for all values of
$w_{0}$ studied in the decoupled limit, $D_{0}=0$. We will also recover these
observations analytically in Sec. V.

Figure 5: Energies obtained from self-consistent HF calculations for selected
candidate states of Table 2. All energies shown are relative to the energy of
the self-consistent SP solution and measured in meV per moiré unit cell. We
take $\epsilon=7$, $d_{s}=40$ nm, and $w_{1}=124$ meV and 4 bands per spin and
valley. We find the IVC- is lower in energy than the SP state for all values
of $w_{0}$ and $D_{0}$. The IVC+ state is always higher in energy than the SP
in agreement with Table 3. The VP state is nearly exactly degenerate with the
SP state, even as $D_{0}$ increases. Both the SLP+ and SLP- states have a
region of energy lower than both the SP and IVC- states at small $w_{0}$ and
large $D_{0}$. While they appear nearly degenerate, a small Hartree
contribution favors the time-reversal odd SLP- over the SLP+ state. However,
this difference between the two sublattice polarized states (< 1 meV) is not
visible in the energies we show.

The lowest energy state for $D_{0}\neq 0$ cannot be directly inferred from
knowledge of the physics of TBG as a finite $D_{0}$ induces hybridization
between the TBG-like and graphene-like bands near the K/K′ points of the MBZ;
it further breaks symmetries in the TBG sector and, hence, changes the basic
form of its dispersion and interaction matrix elements (form factors).
Consequently, it is not clear whether the ground state in the TBG-like and/or
graphene-like sector changes with increasing $D_{0}$. As can be seen in Figs.
5 and 6 as well as in the corresponding phase diagram in Fig. 4, we find
within HF that the IVC- remains the ground state for an extended range of
$D_{0}$, which increases with $w_{0}$. For reference, the range of $D_{0}$ in
Figs. 4, 5, and 6 when combined with additional studies at larger $D_{0}$ in
Appendix E corresponds roughly to the range of displacement fields studied
experimentally in Ref. Park  _et al._ , 2021.

Beyond the critical value of $D_{0}$ for the IVC-, the sublattice-polarized
(SLP) group of states (SLP±, SSLP-) dominates. While these latter three states
are almost degenerate for all parameters $D_{0}$, $w_{0}$ studied, there is a
slight preference towards the time-reversal-odd SLP- (quantum Hall),
predominantly associated with the Hartree energy. While a SP or VP phase does
not appear as a ground state in Fig. 4, we find the energetically close SP or
VP states are either the second or third lowest energy state to the IVC- and
SLP group across the phase diagram. In the next subsection, we will recover
many of these features analytically by investigating the aforementioned
energetic contributions perturbatively.

Figure 6: (a) Linecuts for fixed $w_{0}$ showing the energy of our self
consistent solutions relative to a SP state as a function of $D_{0}$. We use
the same parameters as in Fig. 5. Note that the SLP± and SSLP- are not exactly
degenerate, as discussed in the main text, but are plotted as a single line
here since the splitting is too small to be visible on the scale of the plot.
(b) Matrix form of our converged IVC- and SLP- orders $P_{\boldsymbol{k}}$ for
$\boldsymbol{k}$ near the $\Gamma$ point and near the K point. The matrix
structure of $P_{\boldsymbol{k}}$ is organized such that the largest block
outlined in blue denotes spin flavor, the next largest block outlined in red
denotes valley flavor, the green block denotes upper and lower bands in the
continuum model and the final two boxes denote the graphene-like and TBG-like
band in the upper-band box and the TBG-like then graphene-like bands in the
lower-band box.

##  V Analytical perturbation theory

Finally, we complement the HF numerics with an analytical study of the
behavior of the energies and the order parameters of the degenerate ground
states of the chiral-flat-decoupled limit ($w_{0}=W_{\text{TBG}}=D_{0}=0$) of
Sec. III and Table 2 when turning on $W_{\text{TBG}}$, $D_{0}$, and $w_{0}$.
We first investigate the ordering tendencies of the graphene-like bands (Sec.
V.1) and their mixing with the TBG-like bands (Sec. V.2), before addressing
the energetic contributions coming from the $D_{0}$,$w_{0}$-induced distortion
of the TBG-like form factors (Sec. V.3) and band structure (Sec. V.4).

###  V.1 Ordering in the graphene-like bands

Since the graphene and graphene-like bands in Fig. 2(b) are highly dispersive,
starting from a flat-band limit, as is natural for the TBG-like sector, is not
possible for the graphene(-like) bands. Since their bands are coupled, without
any band gap, to the TBG-like bands with a high density of states, treating
the interactions between the two subsystems as a perturbation is also not
necessarily controlled for realistic parameters. Instead, we here use a
different control parameter: based on the band structure, we expect the effect
of the TBG(-like) bands on the graphene(-like) bands to be the strongest
around the K/K′ points and very weak away from it. To formalize this, let us
assume that the graphene-like bands remain in their filled-lower-bands state
away from the K/K′ points, but allow them to be “deformed” in the region
$\mathcal{A}=\mathcal{A}_{+}\cup\mathcal{A}_{-}$ of the MBZ, where
$\mathcal{A}_{+}$ ($\mathcal{A}_{-}$) are simply connected and centered around
the K (K′) point. As explained in detail in Appendix F.3, we study the energy
of symmetry-allowed ordering tendencies of the graphene-like bands for the
different candidate phases in Table 2 in the limit where the area of
$\mathcal{A}$ is small (compared to that of the MBZ).

To illustrate this procedure, let us consider the SP state. Since the TBG-like
bands break spin-rotation symmetry, it is natural to assume that the same
happens to the graphene-like bands in $\mathcal{A}$; postponing the discussion
of mixing between the bands to Sec. V.2, this means that
$(Q^{g}_{\boldsymbol{k}})_{\eta,\eta^{\prime}}=\delta_{\eta,\eta^{\prime}}\sigma_{0}\boldsymbol{s}$,
$\boldsymbol{k}\in\mathcal{A}_{\eta}$, in the notation introduced in Sec.
III.3, while
$(Q^{g}_{\boldsymbol{k}})_{\eta,\eta^{\prime}}=-\delta_{\eta,\eta^{\prime}}\sigma_{3}$
for all other $\boldsymbol{k}$. Here we already anticipated (as is also
readily checked within this formalism) that it is energetically more favorable
if only the graphene valley $\eta=+$ ($\eta=-$) that is at low energies in the
region $\mathcal{A}_{+}$ ($\mathcal{A}_{-}$) exhibits spin polarization.
Denoting the linear size of $\mathcal{A}_{\pm}$ by $\Delta k>0$, the structure
of the energetic change associated with the deformation of the graphene order
is asymptotically given by

$\displaystyle\begin{split}\Delta E&\sim d_{1}|D_{0}|(\Delta
k)^{2}+d_{2}(\Delta k)^{3}\\\ &\qquad+I_{gg}(\Delta
k)^{3}-I_{gb}D_{0}^{2}(\Delta k)^{2}\end{split}$ (34)

for small $\Delta k$ and $D_{0}$; here, $d_{1,2}>0$ are positive constants
(independent of $\Delta k$ and $D_{0}$, but dependent on $w_{0}$) associated
with the graphene-like dispersion, while $I_{gg}$ depends on the graphene-
graphene form factors $F^{gg}$ and $I_{gb}>0$ on $F^{gb}$ and $F^{gg}$. The
explicit form and derivation can be found in Appendix F.3.

From Eq. (34), we can read off the behavior of the graphene-like bands of the
SP state. For $D_{0}=0$, we see that $\Delta E>0$ (deformation is
energetically disfavored) as long as $d_{2}>I_{gg}$. While this is what our
numerics shows, we point out that $d_{2}<I_{gg}$ would imply that single-layer
graphene spontaneously magnetizes, which is known to be not the case. This
agrees with our result in Sec. III based on adibatically turning on the
coupling between the TBG and graphene system and our numerics which displays
unpolarized graphene Dirac cones in Fig. 3 at $D_{0}=0$. When $D_{0}\neq 0$,
we see that the energetic cost coming from the dispersion in the first line of
Eq. (34) will always overcompensate the energetic gain coming from the
scattering between the graphene- and TBG-like bands (at least in a finite
range of $D_{0}\neq 0$). This is also consistent with our HF numerics: as can
be seen in the band structure in Fig. 3, the graphene-like bands are not spin-
polarized for non-zero $D_{0}$ [see also Fig. 17 where $Q_{\boldsymbol{k}}$ as
obtained in HF are shown].

The ordering tendencies in the graphene-like bands for all other states in
Table 2 can be analyzed in a similar way, see Appendix F.3. For instance, for
the VP we find that the graphene-like bands are found to not develop valley
polarization either—both for $D_{0}=0$ and $D_{0}\neq 0$ and again consistent
with the HF numerics. The same holds for the IVC± states, where breaking of
$\text{U}(1)_{v}$ in the graphene-like bands is strongly suppressed due to the
fact that, at the K (K′) point, the $\eta=-$ ($\eta=+$) graphene-like band
energies are large (already at $D_{0}=0$), see Fig. 2.

The situation is different for the SLP group of states (SLP± and SSLP±). To
begin with the SLP± states, we have
$(Q^{g}_{\boldsymbol{k}})_{\eta,\eta^{\prime}}=\delta_{\eta,\eta^{\prime}}(-\sigma_{3}\cos\theta_{\boldsymbol{k}}+\sigma_{2}\sin\theta_{\boldsymbol{k}})$,
$\boldsymbol{k}\in\mathcal{A}_{\eta}$, and we need to determine
$\theta_{\boldsymbol{k}}$, with
$\theta_{\boldsymbol{k}}=\mp\theta_{-\boldsymbol{k}}$ for SLP± due to $\Theta$
and $C_{2z}$, by minimizing the energy. The change of the energy as a
consequence of this “deformation” is of the form

$\Delta
E_{gg}[\theta_{\boldsymbol{k}}]=\sum_{\eta}\sum_{\boldsymbol{k}\in\mathcal{A}_{\eta}}\left[A_{\boldsymbol{k},\eta}(1-\cos\theta_{\boldsymbol{k}})+B_{\boldsymbol{k},\eta}\sin\theta_{\boldsymbol{k}}\right]$
(35)

for both SLP+ and SLP-. The explicit expressions for $A_{\boldsymbol{k},\eta}$
and $B_{\boldsymbol{k},\eta}$, which are functions of $w_{0}$ and $D_{0}$, are
given in Appendix F.3. From these expressions, it follows that
$B_{\boldsymbol{k},\eta}=0$ and $A_{\boldsymbol{k},\eta}>0$ for $D_{0}=0$ such
that $\Delta E[\theta_{\boldsymbol{k}}]$ is minimized when
$\sin\theta_{\boldsymbol{k}}=0$ (with $\Delta E=0$) and there is no order in
the graphene bands—again in agreement with Sec. III and the HF numerics. Once
$D_{0}\neq 0$, we get $B_{\boldsymbol{k},\eta}\neq 0$ and $\Delta E<0$ by
choosing a profile with $\sin\theta_{\boldsymbol{k}}\neq 0$. In other words,
the graphene-like bands will develop SLP± order for any non-zero $D_{0}$, as
is visible in the $Q_{\boldsymbol{k}}$ for the SLP- state close to the K point
shown in Fig. 6(b). This gaps out the graphene cone, as can also be seen in
our HF band structure in Fig. 3. The energetic gain scales as $g_{\pm}(\Delta
k)^{2}|D_{0}|^{3}$, $g_{\pm}>0$, for small $D_{0}$ and $\Delta k$ for the SLP±
state; the prefactors differ, $g_{+}\neq g_{-}$, due to the symmetry-imposed
constraint $\theta_{\boldsymbol{k}}=\mp\theta_{-\boldsymbol{k}}$. For small
$w_{0}$, we can show that $B_{\boldsymbol{k},\eta}>0$ such that $g_{-}>g_{+}$,
i.e., the SLP- state can gain more energy than the SLP+.

For the SSLP- state and its Hund’s partner SSLP+, which we discuss here
explicitly for reasons that will become clear shortly, we have
$(Q^{g}_{\boldsymbol{k}})_{\eta,\eta^{\prime}}=\delta_{\eta,\eta^{\prime}}(-\sigma_{3}\cos\theta_{\boldsymbol{k}}+\sigma_{2}\boldsymbol{s}\sin\theta_{\boldsymbol{k}})$,
$\boldsymbol{k}\in\mathcal{A}_{\eta}$. Symmetry imposes
$\theta_{\boldsymbol{k}}=\mp\theta_{-\boldsymbol{k}}$ for SSLP±. The
deformation-related energy change $\Delta E$ is found to be again of the form
of Eq. (35); however, we here obtain $B_{\boldsymbol{k},\eta}=0$ for any
$D_{0}$ or $w_{0}$. This is consistent with
$\text{SU}(2)_{+}\times\text{SU}(2)_{-}$ which requires that SSLP± have the
same energy. We, hence, have shown that the graphene-like bands do not develop
any direct SSLP± order and gain energy in the process. They will, however,
hybridize with the TBG-like bands, as we discuss in the next subsection, which
also gaps out the Dirac cones.

In Table 3, we summarize our findings that only the SLP± states can benefit
from ordering the graphene-like bands and how the respective energy gain
scales.

###  V.2 Mixing between the bands

As anticipated above, we next look into the mixing between the TBG-like and
graphene-like bands. While it is clear by symmetry that the mixing has to
vanish for $D_{0}=0$ and be generically present for $D_{0}\neq 0$, we here
investigate the associated energetic gain and precise from of the band mixing
for our candidate states in Table 2. As before, we here outline the basic
strategy and discuss the major results of this calculation and defer the
details to Appendix F.4.

Our starting point are product states characterized by a correlator with
$Q_{\boldsymbol{k}}$ as given in Eq. (27), i.e., without any coherence between
the TBG-like and graphene-like bands. In accordance with our analysis of Sec.
III, we will take $Q^{g}_{\boldsymbol{k}}=-\sigma_{z}$ and let
$Q^{b}_{\boldsymbol{k}}$ be any of the candidate orders. To introduce
momentum-dependent coherence between these sets of bands, we “deform”
$Q_{\boldsymbol{k}}$ by a unitary transformation, $U_{\boldsymbol{k}}$, and
take

$Q^{\prime}_{\boldsymbol{k}}=U^{\phantom{\dagger}}_{\boldsymbol{k}}Q_{\boldsymbol{k}}U^{\dagger}_{\boldsymbol{k}},\quad
U_{\boldsymbol{k}}=e^{i\Lambda_{\boldsymbol{k}}},\quad\Lambda_{\boldsymbol{k}}^{\dagger}=\Lambda_{\boldsymbol{k}}^{\phantom{\dagger}}\in\mathbb{C}^{16\times
16},$ (36)

as ansatz for the correlator. Our goal will be to find the optimal momentum-
dependent $\Lambda_{\boldsymbol{k}}$ to minimize the energy. Since we are
interested in band mixing, we will restrict $\Lambda_{\boldsymbol{k}}$ to act
as a superposition of $\zeta_{1}$ and $\zeta_{2}$,
$\Lambda_{\boldsymbol{k}}=\sum_{j=1,2}\mathcal{M}_{\boldsymbol{k},j}\zeta_{j}$,
$\mathcal{M}_{\boldsymbol{k},j}\in\mathbb{C}^{8\times 8}$, with $\zeta_{j}$
denoting Pauli matrices acting in the space of TBG-like and graphene-like
bands (with index $t=b,g$). Furthermore, $\Lambda_{\boldsymbol{k}}$ will be
constrained by the symmetries of the state under consideration.

To illustrate the procedure, let us focus on the IVC- since this state was
found to be dominant in the HF numerics. Choosing
$Q^{b}_{\boldsymbol{k}}=\sigma_{y}\eta_{y}$ for concreteness, this state
preserves the $C_{2z}\Theta$ symmetry of Table 1, which forces
$\Lambda_{\boldsymbol{k}}$ to obey
$(\Lambda_{\boldsymbol{k}})^{*}=-\Lambda_{\boldsymbol{k}}$. Furthermore noting
that the IVC- state does not break the SU(2)s symmetry (while postponing the
consequences of $C_{2z}$ which will relate $\Lambda_{\boldsymbol{k}}$ and
$\Lambda_{-\boldsymbol{k}}$), it follows that $\Lambda_{\boldsymbol{k}}$ has
to be a momentum-dependent superposition of the $16$ generators

$\zeta_{x}\sigma_{y}\eta_{0,x,z},\quad\zeta_{x}\sigma_{0,x,z}\eta_{y},\quad\zeta_{y}\sigma_{0,x,z}\eta_{0,x,z},\quad\zeta_{y}\sigma_{y}\eta_{y}.$
(37)

To simplify further, we can focus on those $8$ linear combinations [see Eq.
(172) for their explicit form] of the terms in Eq. (37) that anti-commute
rather than commute with $Q_{\boldsymbol{k}}$ of the IVC- state.

Table 3: Summary of the different energetic contributions, as discussed in Sec. V, for the candidate orders in Table 2 when tuning away from the chiral-flat-decoupled limit ($w_{0}=W_{\text{TBG}}=D_{0}=0$). By construction, the energies are identical for the respective Hund’s partners and are, hence, omitted. Here $\Delta E_{bb}(W_{\text{TBG}}=0)$ is the change of energy, relative to SP phase, coming from the modifications of the form factors in the TBG-like bands when turning on $w_{0}$ and $D_{0}$ in the flat limit ($W_{\text{TBG}}=0$). As indicated in the column labelled $\Delta E_{gg}$, only the SLP± states gain energy by ordering in the graphene-like bands immediately when $D_{0}\neq 0$. In all cases, $D_{0}\neq 0$ leads to mixing between the bands, with energetic gain as listed in the chiral limit, $w_{0}=0$, in the column $\Delta E_{gb}$; here $\Delta k$ is the linear size of the fraction $\mathcal{A}_{\pm}$ of the MBS where ordering in the graphene-like bands and mixing take place. The four columns with $\bar{E}_{j}^{b}$ indicate which of the four contributions to the TBG dispersion in Eq. (6) can lower the energies to second order in $W_{\text{TBG}}$. Finally, in the last columns, we list the energy change associated with a finite value of the Hund’s coupling in Eq. (30) for the respective state ($\Delta E_{J}$) and its Hund’s partner ($\Delta E^{H}_{J}$), if it exists, in the chiral-flat-decoupled limit. All coefficients obey $c_{j},g_{\pm},g_{j},\beta_{j}>0$ and explicit expressions can be found in Appendix F. Type | $Q^{b}$ | $\widetilde{Q}^{b}$ | $\Delta E_{bb}(W_{\text{TBG}}=0)$ | $\Delta E_{gg}$ | $\Delta E_{gb}$ | $\,\bar{E}^{b}_{0}\,$ | $\bar{E}^{b}_{1}w_{0}$ | $\bar{E}^{b}_{2}D_{0}$ | $\bar{E}^{b}_{3}w_{0}D_{0}$ | $\Delta E^{\phantom{H}}_{J}$ | $\Delta E^{H}_{J}$
---|---|---|---|---|---|---|---|---|---|---|---
SP | $\sigma_{0}\eta_{0}\boldsymbol{s}$ | $\widetilde{\sigma}_{0}\eta_{0}\boldsymbol{s}$ |  $0$ (by definition) | $0$ | $-g_{1}(\Delta k)^{2}D_{0}^{2}$ | ✗ | ✗ | ✗ | ✗ | $2J_{H}\beta_{1}$ | $-2J_{H}\beta_{1}$
VP | $\sigma_{0}\eta_{z}s_{0}$ | $\widetilde{\sigma}_{0}\eta_{z}s_{0}$ | $0$ | $0$ | $-g_{1}(\Delta k)^{2}D_{0}^{2}$ | ✗ | ✗ | ✗ | ✗ | $0$ | —
IVC+ | $\sigma_{0}\eta_{x,y}$ | $\widetilde{\sigma}_{0}\eta_{x,y}$ | $c_{1}D_{0}^{4}+c_{3}w_{0}^{2}$ | $0$ | $-g_{2}(\Delta k)^{2}D_{0}^{2}$ | ✗ | ✓ | ✓ | ✗ | $-3J_{H}\beta_{2}$ | $J_{H}\beta_{2}$
IVC- | $\sigma_{y}\eta_{x,y}$ | $\widetilde{\sigma}_{z}\eta_{x,y}$ | $c_{1}D_{0}^{4}+c_{4}w_{0}^{2}D_{0}^{4}$ | $0$ | $-g_{2}(\Delta k)^{2}D_{0}^{2}$ | ✓ | ✓ | ✗ | ✗ | $-3J_{H}\beta_{2}$ | $J_{H}\beta_{2}$
SLP- | $\sigma_{y}$ | $\widetilde{\sigma}_{z}$ | $c_{3}w_{0}^{2}+c_{4}w_{0}^{2}D_{0}^{4}$ | $-g_{-}(\Delta k)^{2}|D_{0}|^{3}$ | $-g_{1}(\Delta k)^{2}D_{0}^{2}$ | ✓ | ✗ | ✓ | ✗ | $0$ | —
SLP+ | $\sigma_{y}\eta_{z}$ | $\widetilde{\sigma}_{z}\eta_{z}$ | $c_{2}D_{0}^{4}+c_{3}w_{0}^{2}+c_{4}w_{0}^{2}D_{0}^{4}$ | $-g_{+}(\Delta k)^{2}|D_{0}|^{3}$ | $-g_{3}(\Delta k)^{2}D_{0}^{2}$ | ✓ | ✗ | ✓ | ✗ | $0$ | —
SSLP- | $\sigma_{y}s_{z}$ | $\widetilde{\sigma}_{z}s_{z}$ | $c_{3}w_{0}^{2}+c_{4}w_{0}^{2}D_{0}^{4}$ | $0$ | $-g_{1}(\Delta k)^{2}D_{0}^{2}$ | ✓ | ✗ | ✓ | ✗ | $2J_{H}\beta_{3}$ | $-2J_{H}\beta_{3}$

The first energetic constraint we will take into account is related to the
fact that the graphene-like bands of valley $\eta=-$ ($\eta=+$) are far away
from the Fermi level at the K (K′) point, see Fig. 2. This means that all
mixing processes in $\Lambda_{\boldsymbol{k}}$ that induce a finite occupation
of the upper or unoccupied states in the lower graphene-like band of valley
$\eta=-$ ($\eta=+$) at the K (K′) point are suppressed. In Appendix F.4, we
show that this is equivalent to demanding that
$(\mathcal{M}_{\boldsymbol{k}})_{(p^{\prime},\eta^{\prime},s^{\prime}),(p,\mp,s)}=0$
for $\boldsymbol{k}\in\mathcal{A}_{\pm}$. This reduces the number of
generators further from $8$ to only $4$ [given in Eq. (178)]. For each of
these $4$ generators, $\Lambda_{j}$, $j=1,2,3,4$, we compute the change of
energy associated with the deformation in Eq. (36) where
$U_{\boldsymbol{k}}=e^{i\varphi_{\boldsymbol{k}}\Lambda_{j}}$. It is found to
be of the form

$\Delta
E_{gb}[\varphi_{\boldsymbol{k}}]=\sum_{\boldsymbol{k}\in\mathcal{A}_{+}}\alpha_{\boldsymbol{k}}\sin^{2}\varphi_{\boldsymbol{k}}+\beta_{\boldsymbol{k}}\sin\varphi_{\boldsymbol{k}}\cos\varphi_{\boldsymbol{k}},$
(38)

where $\alpha_{\boldsymbol{k}}$ and $\beta_{\boldsymbol{k}}$ are expressions
involving the form factors, the interaction $V(\boldsymbol{q})$, and the band
structure [see Appendix F.4 for explicit form], and thus depend on $D_{0}$ and
$w_{0}$. In Eq. (38), we have already taken into account the $C_{2z}$
constraint that allowed us to write it as a sum over
$\boldsymbol{k}\in\mathcal{A}_{+}$ only.

First, we find $\beta_{\boldsymbol{k}}=0$ for $D_{0}=0$ and, hence, no mixing
between the bands as expected by symmetry. However, even when $D_{0}\neq 0$,
we still obtain $\beta_{\boldsymbol{k}}=0$ for two of the four generators,
while $\beta_{\boldsymbol{k}}\neq 0$ for the other two (say $\Lambda_{1,2}$)
in the chiral limit ($w_{0}=0$). Consequently, at least when $w_{0}$ is not
too large, the generators $\Lambda_{1,2}$ have to dominate. Which of those
remaining two is dominant, cannot be determined purely analytically as it will
depend on non-universal values of the form factors, so let us compare with
numerics. The mixing matrix, defined as

$\left(M^{Q}_{\boldsymbol{k}}\right)_{(p,\eta,s),(p^{\prime},\eta^{\prime},s^{\prime})}:=\left(Q_{\boldsymbol{k}}\right)_{(b,p,\eta,s),(g,p^{\prime},\eta^{\prime},s^{\prime})}$
(39)

and can be computed straightforwardly for the two candidate generators
$\Lambda_{1,2}$. We find

$M^{Q}_{\boldsymbol{k}}=2\sqrt{2}i\left(\sigma_{x}P_{\pm}\pm\sigma_{0}\eta^{\mp}\right)\varphi_{\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{A}_{\pm},$
(40a) and $M^{Q}_{\boldsymbol{k}}=2\sqrt{2}i\left(\sigma_{z}P_{\pm}\pm
i\sigma_{y}\eta^{\mp}\right)\varphi_{\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{A}_{\pm},$
(40b)

for $\Lambda_{1}$ and $\Lambda_{2}$, respectively. Here, we defined
$P_{\pm}=(\eta_{0}\pm\eta_{z})/2$ and $\eta^{\pm}=(\eta_{x}\pm i\eta_{y})/2$.
Intuitively, the first mixing matrix means that the lower (upper) graphene-
like bands of the valley that is at low energies at K or K′ mixes with the
upper (lower) and lower (upper) TBG-like bands of the same and opposite
valley, respectively. The second option in Eq. (40b) describes the “twisted”
situation where the lower (upper) graphene-like band mixes with the lower
(upper) TBG-like bands of the same valley and the upper (lower) band in the
opposite valley. The HF result for $Q_{\boldsymbol{k}}$ close to the K point,
shown in Fig. 6(b), is consistent with Eq. (40a).

Minimizing Eq. (38) and noting that $\beta_{\boldsymbol{k}}$ scales linearly
with $D_{0}$ for small $D_{0}$ while $\alpha_{\boldsymbol{k}}>0$ at $D_{0}=0$,
we find the scaling of the energetic gain to be $\Delta E_{gb}\sim-
g_{3}(\Delta k)^{2}D_{0}^{2}$, where $g_{3}>0$; this is indicated in Table 3.

In the same way, all other candidate states can be studied, see Appendix F.4.
In accordance with our expectation based on symmetry, we find no mixing and
vanishing energetic gain, $\Delta E_{gb}=0$, when $D_{0}=0$, while mixing
immediately sets in and $\Delta E_{gb}<0$ once $D_{0}\neq 0$. In all cases,
$\Delta E_{gb}$ scales as $\Delta E_{gb}\sim-g_{j}(\Delta k)^{2}D_{0}^{2}$,
$g_{j}>0$, as $D_{0}\rightarrow 0$. While $g_{j}=g_{j}(w_{0})$ depend on
$w_{0}$, we list $\Delta E_{gb}$ in the chiral limit, $w_{0}=0$, in Table 3 as
it reveals some structure: we see that the energetic gain coming from
hybridization is identical for the two IVCs and it is the same for the four
states SP, VP, SLP-, and SSLP- in the chiral limit.

The mixing matrices $M^{Q}$ in Eq. (39) for these states provide some
additional consistency checks between analytics and numerics. For instance,
the IVC+’s analysis closely parallels the one outlined above for the IVC-: out
of $16$ symmetry-allowed generators, only two candidate combinations remain.
Interestingly, their associated $\alpha_{\boldsymbol{k}}$ and
$\beta_{\boldsymbol{k}}$ in Eq. (38) are the same (to leading order in
$D_{0}$) as those of the two IVC- candidates. Having identified
$M^{Q}_{\boldsymbol{k}}$ in Eq. (40a) as being dominant by comparison with
numerics, we can read off which of the two analytical candidates for the IVC+
must be realized; indeed, we find the same one in numerics [compare Eq. (185)
and Fig. 17]. To provide a second example, we find a mixing matrix for the
SLP- given by

$M^{Q}_{\boldsymbol{k}}=2i(\sigma_{y}+\sigma_{z})P_{\pm}\varphi_{\boldsymbol{k}},\quad\boldsymbol{k}\in\mathcal{A}_{\pm}.$
(41)

This means that there is only mixing in the valley for which the graphene-like
bands are at low energies and that both upper and lower graphene-like bands
mix with both TBG-like bands. This is consistent with the $Q_{\boldsymbol{k}}$
in Fig. 6(b) we found for the SLP- in HF close to the K point.

###  V.3 Deforming the form factors

Apart from the mixing with and the “proximity-induced” order in the graphene-
like bands, there are also important energetic contributions directly within
the TBG-like bands, when tuning away from the chiral-flat-decoupled limit.
While some of these contributions are analogous to TBG Lian  _et al._ (2020);
Bultinck  _et al._ (2019), some others are not: non-zero $D_{0}$ strongly
breaks $P$ (see Table 1) and, hence, induces terms in the TBG-like band
structure, see second line in Eq. (6), without any analogue in TBG; it also
leads to terms in the form factors [see Eq. (79a) for details] that cannot be
present in TBG and have not been studied in the literature.

We here begin with the impact of the form factors and restrict ourselves, for
now, to the flat limit, $W_{\text{TBG}}=0$. In the column labeled $\Delta
E_{bb}$ in Table 3 we show the energy of the candidate states for non-zero
$w_{0}$ and $D_{0}$ relative to the SP [see Appendix F.1 for explicit
expressions for the prefactors $c_{j}>0$]. First, for $D_{0}=0$, we recover
the previous result Bultinck  _et al._ (2019) that the IVC- and VP are the
only states besides the SP that are not penalized when turning on $w_{0}$.
Second, we see that this changes once $D_{0}\neq 0$: if $w_{0}=0$, it is the
SLP group of states that is not suppressed by breaking $P$ with $D_{0}$
111Note that the term $c_{2}$, which comes from the Hartree term, turns out to
be numerically very small., while both IVCs have increasing energy. When both
$w_{0}$ and $D_{0}$ are simultaneously non-zero, also the SLP group is
suppressed (by the exact same amount as the IVC-) compared to the SP and VP.
Algebraically, this is related to the fact that the SP and VP order parameters
are the only ones that commute with all form factors once $w_{0},D_{0}\neq 0$
[cf. Eq. (79a)].

###  V.4 Finite TBG bandwidth

Finally, we take into account the finite bandwidth of the TBG-like bands by
doing perturbation theory in $W_{\text{TBG}}$ in Eq. (5), starting from the
product states associated with the candidate orders of Table 2. We will only
outline the results here and refer the interested reader to Appendix F.2.

Since the correction to the relative energies of the different candidate
orders vanishes to first order in $W_{\text{TBG}}$, we focus on second order
perturbation theory, which, if non-zero, will always lower the energy of the
states and can be thought of as “superexchange”. We consider the superexchange
processes associated with all four terms, $E^{b}_{j}$, $j=0,1,2,3$, in Eq. (6)
of the TBG-like dispersion; the corresponding energetic gain will scale as
$W_{\text{TBG}}^{2}/U$ for $j=0$, $W_{\text{TBG}}^{2}w_{0}^{2}/U$ for $j=1$,
$W_{\text{TBG}}^{2}D_{0}^{2}/U$ for $j=2$, and
$W_{\text{TBG}}^{2}w_{0}^{2}D_{0}^{2}/U$ for $j=3$ to leading order in $D_{0}$
and $w_{0}$, where $U$ is the energy scale associated with occupying a
$\boldsymbol{k}$-state of an unoccupied flavor in a given ground state. Which
of these four superexchange processes are “active” for the candidate states is
listed in Table 3. We find that, by virtue of being proportional to the
identity in Eq. (6), $E^{b}_{3}$ does not affect the energy of the states.
Furthermore, we see that the displacement-induced superexchange process favors
the IVC+ and the SLP group of states.

###  V.5 Comparison of energetics

Taken together, the energetics obtained analytically as summarized in Table 3
agrees well with the numerics in several aspects: we can see that the energies
of the SLP+, SLP-, and SSLP- are expected to be very close, exactly as seen in
HF, see Fig. 6. There is a very small splitting between the three associated
with the fact that the SSLP- cannot benefit from ordering in the graphene-like
bands, that that effect is weaker for the SLP+ than for the SLP- (recall
$g_{-}>g_{+}$, for small $w_{0}$), and that the SLP+ is slightly suppressed by
the Hartree term. We have seen that there is a small energetic preference
towards the SLP- in the numeris as well. In addition, we can also read off
from Table 3 that the SLP group of states should be preferred for small
$w_{0}$ and large $D_{0}$, which is consistently seen in our HF phase diagram
in Fig. 4.

Notwithstanding the good agreement between the HF and the analytics,
concerning the graphene-ordering, mixing between the bands, and the
energetics, there are also differences. These can be traced back to the
additional presence of remote bands, in particular, the bands with
$\sigma_{h}$ eigenvalues $+1$ just above (below) the almost flat TBG-like band
$p=+$ ($p=-$) away from the K/K′ points. While these are included in our HF
numerics, they are not taken into account in the analytics. Most notably, we
have seen that mixing between these bands and the almost flat TBG-like bands
further lowers the energy of IVC- state relative to the SP and VP states.

###  V.6 Breaking SU(2)${}_{+}\times$ SU(2)-

Finally, we come back to the fact that certain Hund’s partners of states are
degenerate in the model we have focused on so far and defined in Sec. II; this
can be traced back to the presence of the SU(2)${}_{+}\times$ SU(2)- spin
symmetry. Turning on a finite Hund’s coupling, $J_{H}\neq 0$, will lift this
degeneracy and also slightly affect the relative energetics of the candidate
states. For the form of the Hund’s coupling defined in Eq. (30), we have
computed the respective change of energy. With more general and explicit
expressions available in Appendix F.5, we present their impact on the energies
of the states in the last two columns of Table 3 in the chiral-flat-decoupled
limit. Here $\beta_{j}$ are positive constants obeying
$\beta_{2}>\beta_{1},\beta_{3}$. While estimates based on the Coulomb
interaction yield an energy change (per unit cell) of
$2J_{H}\beta_{1}\simeq-0.2\,\text{meV}$ for the SP state in TBG Chatterjee
_et al._ (2020), its actual effective value might be smaller due to screening
processes and additional electron-phonon coupling, which can also change the
sign of $J_{H}$.

However, irrespective of the sign of $J_{H}$ and its precise magnitude, we see
that it will favor the SSLP± states in the otherwise almost degenerate SLP
group of states, since the SLP± cannot lower their energy with $J_{H}\neq 0$.
This is due to the fact that SLP± are their own Hund’s partners. Similarly,
the SP and SVP states will be favored over the VP phase. Finally, we note that
the energetic gain coming from $J_{H}$ is of the same order of magnitude (even
larger) for the SIVC- (IVC-) compared to the SP (SVP) state. Consequently, we
do not expect that $J_{H}$ will change the fact that the IVC- dominates over
the SP in the phase diagram in Fig. 4.

##  VI Numerics for $\nu=\pm 2$

In this section, we will discuss the numerical results for correlated states
at half-filling of the lower or upper TBG-like bands, $\nu=\pm 2$.

###  VI.1 Procedure and results

We find self-consistent HF solutions with energies lower than the symmetry-
unbroken normal state where, at $D_{0}=0$ ($D_{0}\neq 0$), the TBG (TBG-like
and graphene-like) bands are spin polarized and, on top of this, exhibit any
of the candidate orders defined for $\nu=0$ in Table 2. The obtained spin
polarization is our explanation for the observed Park  _et al._ (2021); Hao
_et al._ (2021) reduced flavor degeneracy setting in around $|\nu|=2$.

To be more explicit, our self-consistent solutions are found by starting, say
for $\nu=-2$, from a correlator
$P^{b}_{\nu=-2}=\frac{1}{2}(\mathbbm{1}+s_{z})\frac{1}{2}(\mathbbm{1}+Q^{b})$
in the TBG-like subspace, where $Q^{b}$ is any of the candidate orders in
Table 2; the initial correlator for the graphene-like sector is taken to be
filled lower bands of the continuum model at charge neutrality,
$P^{g}=\frac{1}{2}(\mathbbm{1}-\sigma_{z})$. The corresponding correlator at
$\nu=+2$ in the TBG-like sector is simply given by
$P^{b}_{\nu=+2}=\mathbbm{1}_{8\times 8}-P^{b}_{\nu=-2}$. We apply the same
iterative procedure as for $\nu=0$ to obtain self-consistent solutions.

We show the band structures resulting from our self-consistent calculation for
the spin-polarized IVC- and spin-polarized SLP states in Fig. 7. The energies
we obtain up to $D_{0}/w_{1}=0.5$ are shown in Appendix D, but are
qualitatively similar to the numerical energies we computed at $\nu=0$; we
note though the spin-valley-polarized and spin-polarized IVC- are closer in
energy at $\nu=2$ than the SVP and IVC- state at $\nu=0$. As can be seen in
Fig. 7, the graphene bands at $D_{0}=0$ remain unpolarized, in agreement with
the analytics in Sec. III. Once a finite displacement field is applied, also
the graphene-like bands develop some small polarization. However, exactly as
for $\nu=0$, the Dirac cones of the IVC- state are not gapped out. We note the
IVC-, which was semimetallic at $\nu=0$, is metallic at $\nu=\pm 2$ (albeit
with small Fermi surfaces). This follows from the additional spin polarization
at $|\nu|=2$, breaking SU(2)s which guaranteed its semimetallic character at
$\nu=0$. In Fig. 7(a,b), we see that the SLP state, which was insulating at
$\nu=0$, becomes semimetallic at $|\nu|=2$, since the spin-polarization of the
TBG-like bands only allows for SLP order in one of the spin species of the
graphene-like bands.

###  VI.2 Connection to experiment

An important feature in the experimental data on MSTG Park  _et al._ (2021);
Hao  _et al._ (2021), which may be related to our numerics, is the observation
of enhanced resistivity at $\nu=0$ for both $D_{0}=0$ and $D_{0}>0$ and a
state with high resistivity observed only at a finite value of $D_{0}$ at
$\nu=2$. Taking our results for band structures at $\nu=0$ and $|\nu|=2$
together, we note that at $\nu=0$, the leading phase in our numerics is either
a semimetallic IVC- state or an insulating SLP state, both of which may be
compatible with the observed high resistivity state at $\nu=0$. At $\nu=2$, we
find our leading solution is a metallic IVC- state for small values of $D_{0}$
and a semimetallic SLP state for large enough $D_{0}$—a possible explanation
for why high resistivity at $|\nu|=2$ only sets in above a finite value of
$D_{0}$ in experiment.

Another notable experimental observation is that the reset of the band
structure at $|\nu|=2$ splits into a Dirac-like feature at $|\nu|=2$ and a
van-Hove singulartiy, associated with a rapid increase and sign change in the
Hall density, at $|\nu|=2-\delta(D_{0})$, $\delta(D_{0})\ll 1$ Park  _et al._
(2021). If we assume that the spin polarization we found persists for a finite
range of $|\nu|$ below $2$, the bandstructures in Fig. 7 provide a natural
explanation: lowering the chemical potential, e.g., in Fig. 7(d), until it
hits the lower, almost flat set of bands will lead to a Lifshitz transition
where the hole pockets around the K and K′ points merge. This could explain
the observed behavior of the Hall density. We note that the tendency, visible
in Fig. 7, that increasing $D_{0}$ pushes these almost flat bands away from
the Fermi level at $|\nu|=2$, is consistent with this feature only being
visible at non-zero $D_{0}$ and the associated $\delta(D_{0})$ increasing with
$|D_{0}|$ in experiment Park  _et al._ (2021). Numerical results for larger
values of $D_{0}$ are discussed in Appendix E.

Figure 7: HF band structures at $\nu=2$ for the spin polarized SLP- state (a)
at $D_{0}=0$ and (b) $D_{0}>0$ and spin polarized IVC- state (c) at $D_{0}=0$
and (d) $D_{0}>0$.

##  VII Superconductivity

Having established the nature of the correlated phases in MSTG at various
filling fractions, we next study the consequences for the superconducting
states.

###  VII.1 Pairing in the presence of polarization

Let us begin with the range of electron filling, $2<|\nu|<3$, where
superconductivity is most prominently observed in experiment Park  _et al._
(2021); Hao  _et al._ (2021); Cao  _et al._ (2021). We have seen in Sec. VI
that spin polarization supplemented with IVC- order is favored for realistic
parameters at $|\nu|=2$. Because the associated reduction of the number of
flavors in the normal state is seen in experiment over the entire or most of
the superconducting range of $\nu$, we assume that spin polarization and
superconductivity co-exist at least in part of the phase diagram. As follows
from the analysis in Ref. Scheurer and Samajdar, 2020, where a classification
of pairing in almost SU(2)${}_{+}\times$ SU(2)--symmetric graphene moiré
systems in the presence and absence of flavor polarization can be found, the
superconducting state has to be in a non-unitary triplet phase—irrespective of
the precise pairing mechanism. We emphasize that this also holds if the
additional IVC- ordering found in the HF at $|\nu|=2$ co-exists with
superconductivity in a finite range of $\nu$: the bands above the Fermi level
in the IVC- band structure in Fig. 7 still exhibit Kramers partners at momenta
$\boldsymbol{k}$ and $-\boldsymbol{k}$ with the same spin, due to the
preserved spinless $\widetilde{\Theta}$ symmetry; these degenerate states can
form Cooper pairs with non-unitary triplet vector. Note that this would not be
the case, e.g., for the VP state [shown in Fig. 13], which does not exhibit
exactly degenerate energy levels at $\boldsymbol{k}$ and $-\boldsymbol{k}$
above the Fermi level. Given the strong tendency of MSTG towards
superconductivity, this VP phase is, thus, a less natural candidate order.
Indeed, we find it to be subleading in our HF numerics, see Sec. VI.

It is important to note, however, that the spin polarization (and additional
IVC- state) is only realized for $J_{H}<0$. While this is expected to be the
case if $J_{H}$ stems entirely from the Coulomb interaction, its sign is
unknown. For $J_{H}>0$, we would, instead, obtain the SVP phase with
additional SIVC- order (cf. Table 2). In that case, the associated
superconductor would also be the corresponding Hund’s partner, which is an
admixture between a singlet and unitary triplet Scheurer and Samajdar (2020).
In the presence of a magnetic field there is a crucial difference between the
two scenarios: while, for $J_{H}<0$, the SP will just align with the Zeemann
field and the superconductor will remain a non-unitary triplet, the two
antiparallel spin polarizations in the two valleys of the SVP state will be
canted gradually; a coexisting singlet-unitary-triplet superconductor will
continuously transition into a non-unitary triplet Scheurer and Samajdar
(2020). For completeness, we demonstrate this explicitly in Appendix G.2.
Based on the estimate $J_{H}\lesssim 0.2\,\text{meV}$ Chatterjee  _et al._
(2020), we obtain $3\,\textrm{T}$ as the characteristic magnetic field
strength of the transition.

###  VII.2 Pairing without polarization

In experiment Park  _et al._ (2021); Hao  _et al._ (2021); Cao  _et al._
(2021), there are also regions of $\nu$ and $D_{0}$ with superconductivity but
without any signs of flavor polarization in the corresponding normal state (SC
II in Fig. 1). We here discuss the nature and origin of these superconducting
phases.

In light of recent experiments in TBG Stepanov  _et al._ (2020); Liu  _et al._
(2021c); Saito  _et al._ (2020), which indicate that electron-phonon coupling
plays an important role for pairing in graphene moiré systems, the picture
proposed in Ref. Samajdar and Scheurer, 2020 provides a very promising
microscopic scenario: while electron-phonon coupling is likely important to
stabilize superconductivity and crucially determines the critical temperature,
it might predominantly mediate an SU(2)${}_{+}\times$ SU(2)- symmetric pairing
interaction. This leads to the near degeneracy of singlet and triplet pairing.
In the additional presence of flavor polarization, its structure determines
whether singlet or triplet is realized, as discussed above. In the absence of
polarization, electron-phonon coupling by itself can never favor triplet over
singlet and at most make the two degenerate Brydon  _et al._ (2014); Scheurer
(2016); however, additional Cooper-channel interactions coming from the
fluctuation of particle-hole orders can determine whether the superconductor
will be in a singlet or triplet state Samajdar and Scheurer (2020). To analyze
this, we will consider the action

$\mathcal{S}=\mathcal{S}_{0}+\mathcal{S}_{\phi}+\mathcal{S}_{\phi
f}+\mathcal{S}_{\text{phonon}},$ (42)

which consists of the bare non-interacting action $\mathcal{S}_{0}$ associated
with $H_{0}$ in Eq. (5) and $\mathcal{S}_{\phi}$, given by

$\mathcal{S}_{\phi}=\frac{1}{2}\int_{q}\phi_{q}^{j}\left[\chi^{-1}(i\Omega_{n},\boldsymbol{q})\right]_{j,j^{\prime}}\phi_{-q}^{j^{\prime}},$
(43)

where $\phi_{q}^{j}$ is the set of real bosonic fields (labelled by $j$)
describing the fluctuations of a given candidate particle-hole order in Table
2. We use $q=(i\Omega_{n},\boldsymbol{q})$ labeling bosonic Matsubara
frequencies $\Omega_{n}$ and momentum $\boldsymbol{q}$, introduce the short-
cut $\int_{q}\ldots\equiv T\sum_{\Omega_{n}}\sum_{\boldsymbol{q}}\ldots$, and
denote the (fully renormalized, low-energy) susceptibility in the particle-
hole channel under consideration by $\chi(i\Omega_{n},\boldsymbol{q})$. The
bosons $\phi_{q}^{j}$ are coupled to the low-energy electron fields
$f_{k;n,\eta,s}$, with $k=(i\omega_{n},\boldsymbol{k})$, via

$\mathcal{S}_{\phi
f}=\int_{q}\int_{k}f^{\dagger}_{k+q}\lambda^{j}(\boldsymbol{k}+\boldsymbol{q},\boldsymbol{k})f^{\phantom{\dagger}}_{k}\,\phi_{q}^{j}.$
(44)

Note that the coupling vertex
$\lambda^{j}(\boldsymbol{k}+\boldsymbol{q},\boldsymbol{k})$ is in general a
matrix in valley, spin, and band space. For instance, a minimal description of
fluctuations of the IVC- state is given by the two-component boson
$\phi^{j}_{q}$, $j=x,y$, with
$[\chi(i\Omega_{n},\boldsymbol{q})]_{j,j^{\prime}}\propto\delta_{j,j^{\prime}}/(\Omega_{n}^{2}+c^{2}\boldsymbol{q}^{2}+\xi^{-2})$,
where the coherence length $\xi$ parametrizes the proximity to the critical
point. Furthermore,
$\lambda^{j}(\boldsymbol{k},\boldsymbol{k}^{\prime})=s_{0}\eta_{j}f_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$,
where $f_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$ is a matrix in band space
($p=\pm$, $t=b,g$) only, which obeys
$f^{\phantom{\dagger}}_{\boldsymbol{k},\boldsymbol{k}^{\prime}}=f^{\dagger}_{\boldsymbol{k}^{\prime},\boldsymbol{k}}=-f^{T}_{\boldsymbol{k}^{\prime},\boldsymbol{k}}$
due to Hermiticity and $\Theta C_{2z}$ symmetry. Finally, the last part
$\mathcal{S}_{\text{phonon}}$ in Eq. (42) stands for the action of the phonons
and their coupling to the electrons.

We integrate out the bosonic modes $\phi_{q}^{j}$ in the action of Eq. (42),
leading to an interaction between the fermions $f$. In the saddle-point
equations of the Cooper channel, this interaction can be viewed as a
SU(2)${}_{+}\times$ SU(2)--symmetry breaking correction to the
SU(2)${}_{+}\times$ SU(2)--symmetric interaction coming from the phonons in
$\mathcal{S}_{\text{phonon}}$.

As shown in Ref. Samajdar and Scheurer, 2020, whether this tips the balance
towards singlet or triplet is determined by the behavior of the fluctuating
modes $\phi_{q}^{j}$ under spinful time-reversal $\Theta_{s}$: if the bosonic
mode is even (odd) under $\Theta_{s}$, as indicated by $+$ ($-$) in the column
$\Theta_{s}$ in Table 2, it will generically favor singlet (triplet) over
triplet (singlet) pairing—in the presence of fine-tuning or additional
symmetries, the two might remain degenerate. We emphasize that this conclusion
does not depend on microscopic details such the precise form of $\chi$ in Eq.
(43) or of $\lambda^{j}$ in Eq. (44). Due to U(1)v symmetry, the system also
exhibits the spinless time-reversal symmetry $\widetilde{\Theta}$, see Table
1. We have generalized the analysis of Ref. Samajdar and Scheurer, 2020 in
Appendix G to also include this form of time-reversal symmetry and proved that
any $\phi_{q}^{j}$ that is even (odd) under $\widetilde{\Theta}$ will
generically favor triplet (singlet). This also implies that any $\phi_{q}^{j}$
which has the same behavior under $\Theta_{s}$ and $\widetilde{\Theta}$ will
keep singlet and triplet degenerate. As expected, this is precisely the case
for all order parameters in Table 2 that are their own Hund’s partner.

First, we focus on the superconducting domes (indicated by SC II in Fig. 1)
that emerge in the vicinity of but outside the region with flavor polarization
and reconstructed band structure Park  _et al._ (2021). As mentioned above,
with Coulomb interactions only, the IVC- order with additional spin
polarization is favored over the Hund’s partner, the SIVC- with additional SVP
order. We can read off from Table 2 that triplet will then be favored over
singlet pairing. As shown in Ref. Scheurer and Samajdar, 2020, the resulting
triplet will be unitary within mean-field theory, but can become non-unitary
when spin fluctuations corrections become significant. These two scenarios can
be distinguished experimentally since the unitary triplet will exhibit a BKT
transition (of a charge-$4e$, spin-rotation-invariant order parameter
combination), while the non-unitary state will not. For the other sign of
$J_{H}$, singlet pairing (if mean-field theory applies) or a mixed singlet-
triplet phase (if the SVP fluctuations dominate) will be realized Scheurer and
Samajdar (2020).

Second, we also comment on possible pairing phases close to charge neutrality,
although these have not been seen experimentally. For the same sign,
$J_{H}<0$, that favors SP at $\nu=2$, we obtain the SIVC- at $\nu=0$ as the
dominant, semimetallic instability at small $D_{0}$ and large $w_{0}$, see
Fig. 4. As can be seen from Table 2, it is even under $\Theta_{s}$ and will,
hence, favor singlet pairing. Among the SLP group of states, realized for
larger $D_{0}$, only fluctuations of the SSLP± states will break the
SU(2)${}_{+}\times$ SU(2)--symmetry; the SSLP-, realized for $J_{H}<0$, will
also favor singlet.

Finally, if electron-phonon coupling does not play any role for pairing in
MSTG, i.e., the term $\mathcal{S}_{\text{phonon}}$ in Eq. (42) can be
neglected, we still obtain the same results (see Appendix G.1 for a
derivation): if $\phi_{q}^{j}$ are even/odd (odd/even) under $\Theta_{s}$