††thanks<EMAIL_ADDRESS>
# Fermionic Monte Carlo study of a realistic model of twisted bilayer graphene

Johannes S. Hofmann  Department of Condensed Matter Physics, Weizmann
Institute of Science, Rehovot 76100, Israel  Eslam Khalaf  Department of
Physics, Harvard University, Cambridge, MA 02138, USA  Ashvin Vishwanath
Department of Physics, Harvard University, Cambridge, MA 02138, USA  Erez Berg
Department of Condensed Matter Physics, Weizmann Institute of Science, Rehovot
76100, Israel  Jong Yeon Lee  Kavli Institute for Theoretical Physics,
University of California, Santa Barbara, CA 93106-4030, USA

###### Abstract

The rich phenomenology of twisted bilayer graphene (TBG) near the magic angle
is believed to arise from electron correlations in topological flat bands. An
unbiased approach to this problem is highly desirable, but also particularly
challenging, given the multiple electron flavors, the topological obstruction
to defining tight-binding models and the long-ranged Coulomb interactions.
While numerical simulations of realistic models have thus far been confined to
zero temperature, typically excluding some spin or valley species, analytic
progress has relied on fixed point models away from the realistic limit. Here
we present for the first time unbiased Monte Carlo simulations of realistic
models of magic angle TBG at charge-neutrality. We establish the absence of a
sign problem for this model in a momentum space approach, and describe a
computationally tractable formulation that applies even on breaking chiral
symmetry and including band dispersion. Our results include (i) the emergence
of an insulating Kramers inter-valley coherent ground state in competition
with a correlated semi-metal phase, (ii) detailed temperature evolution of
order parameters and electronic spectral functions which reveal a ‘pseudogap’
regime, in which gap features are established at a higher temperature than the
onset of order and (iii) predictions for electronic tunneling spectra and
their evolution with temperature. Our results pave the way towards uncovering
the physics of magic angle graphene through exact simulations of over a
hundred electrons across a wide temperature range.

##  I Introduction

The discovery of interaction-driven insulating and superconducting phases in
twisted bilayer graphene (TBG) [1, 2] has inspired intense experimental and
theoretical effort aiming to understand the nature of these phases and the
relationship between them [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40]. Such an understanding remains extremely challenging due
to the strongly interacting nature of the problem combined with the non-
trivial band topology [6, 17, 41, 42, 43, 44, 45] which obstructs a
conventional lattice description. In addition, the presence of spin and valley
results in a very large manifold of possible competing states which further
complicates the analysis. Finally, remarkable phenomena are further observed
at finite temperature, ranging from a flavor cascade [46, 47] to the
Pomeranchuk effect [48, 49] and unusual ‘strange metal’ scaling of resistivity
with temperature [50, 51] which clamor for explanation.

To tackle this problem, various approaches have been employed which can be
divided into two categories. First, there are real-space approaches based on
deriving Hubbard-like models [52, 10, 53]. Such models have been analyzed in
strong coupling approaches extended with perturbation theory [54, 53] in
addition to numerical methods such as density matrix renormalization group
(DMRG) [55] and QMC [37]. However, all these approaches have to deal with the
inherent problem that the band topology is obscured in the real space model
and the symmetries are not naturally represented. Furthermore, it is not
always straightforward to relate the parameters of these models to the
realistic problem. The other class of approaches employs a momentum space
continuum description based on the Bistritzer-MacDonald (BM) model [56]. These
approaches include strong coupling perturbative approaches [27, 57, 58], self-
consistent Hartree-Fock [11, 59, 27, 60, 61, 62, 63, 64], DMRG [29, 30, 65],
exact diagonalization (ED) [66, 67] and many others [7, 68, 69]. Each of these
approaches comes with its own limitations; strong coupling perturbative
approaches requires an extrapolation to access the intermediate coupling
regime which is likely to be of experimental relevance, mean-field approaches
are biased towards Slater determinant states and overestimate symmetry-
breaking tendencies, DMRG is largely restricted to zero temperature
properties, and ED can only access very small system sizes.

Figure 1:  (a) Electron-hole excitation spectrum at charge neutrality
[eigenvalues of $\hat{H}_{0}$ in Eq. (III.2)]. In the DQMC simulations, we
project the Hamiltonian onto the middle two narrow bands (red and blue). Solid
(Dashed) lines are for + ($-$) valley. (b) Bandwidth of $\hat{H}_{0}$ as
function of $\kappa=w_{0}/w_{1}$ (red), along with the bandwidth of the
excitation spectrum in the interaction-only limit, where $\hat{H}_{0}$ is
neglected (blue). (c-f) DQMC results for $\kappa=0.75$, for which
$W=8\,\textrm{meV}$: (c) K-IVC correlation lengths for two system sizes
$L=3,6$ and $9$ in the unit of moiré lattice vector, (d) single-particle
density of state
$\rho(\omega)=\int\differential^{2}{\bm{k}}\,A(\omega,{\bm{k}})$, (e,f)
spectral function $A(\bm{k},\omega)$ along a cut through the moiré Brillouin
zone for $T=29$K and $1.8$K, respectively. The red crosses denote the single
particle dispersion $\hat{H}_{0}$, while the blue crosses represent the single
particle dispersion in the interaction only limit. (g-j) Same as panels (c-f)
for $\kappa=0.3$ ($W=14\,\textrm{meV}$). The parameters used in the
simulations are given in Sec. V.1. In (g), the empty square represents the
zero temperature data obtained by a projective variant of the QMC method [70,
71].

In this work, we overcome these limitations by using a determinant quantum
Monte Carlo (DQMC) method to study the behavior of the strongly interacting
flat bands in TBG at charge neutrality (CN) ($\nu=0$) where we establish the
absence of the sign problem. Our approach is based on a momentum space
description which enables it to accurately capture non-trivial topology of
flat bands and evade the problems that plague lattice real-space descriptions.
In addition, it allows us to handle the long-range Coulomb interaction and is
ideally suited to address the physics at finite temperatures. Experimentally,
the behavior of TBG at CN (in the absence of substrate alignment) is sample-
dependent with some samples exhibiting semi-metallic transport [1, 72] while
others showing an activated behavior pointing to an insulating gap [73, 74].
Other probes such as scanning tunneling microscopy (STM) report a strong
reconstruction of density of states and enhancement of the separation between
van Hove peaks at CN, but, in the absence of further input, does not
conclusively establish a gap [75, 76, 77, 78]. The different behavior between
samples could arise from substrate alignment, disorder[79] or strain which was
argued to induce competition between a correlated insulator and a nematic
semimetal [59, 65]. In the absence of these effects, Ref. [27] have proven
that the ground state is necessarily an insulator in the limit of small
dispersion. On the other hand, for the physically realistic parameters where
interactions themselves give rise to a significant dispersion, such a
correlated insulator can be destabilized. Therefore, unbiased numerical
simulations are needed even to establish the phase diagram of pristine twisted
bilayer graphene.

Furthermore, access to the physics at finite temperature can shed light on a
host of phenomena mentioned previously. More broadly, TBG represents a new
kind of correlated insulator which occurs in topological bands, and revealing
the detailed formation of such a state, as we outline here, can have far
reaching consequences.

The rest of the paper is organized as follows. In Sec. II, we summarize the
main results of the paper. In Sec. III, we examine the TBG Hamiltonian in
detail and define the interacting Hamiltonian to be used in DQMC. As pointed
out in [11, 59, 27], the bare kinetic energy receives corrections in the
interacting problem which needs to be duly incorporated. In Sec. IV we revisit
the sign-problem in DQMC. We elucidate the role played by a particle-hole
symmetry in sign problem. Using this, we show that the sign-problem is absent
in TBG at charge neutrality under a minor approximation, which amounts to
neglecting the small angle rotation of the Dirac matrices [44, 42]. In Sec. V,
we present DQMC results with various parameter, system sizes, and
temperatures. In Sec. VI, we conclude the paper with an outlook.

##  II Summary of results

Our starting point is the continuum model describing the moiré bands in TBG
[56] at charge neutrality, supplemented by long-range Coulomb interactions.
This model respects all the relevant symmetries and the non-trivial topology
of the low-energy bands. For computational simplicity, we keep only the pair
of bands closest to the Fermi level [see Fig. 1(a)], and project out all the
other bands, assumed to be either empty or full 111Including additional bands
does not introduce a sign problem, and can in principle be incorporated in the
DQMC simulations. However, the computational complexity increases as the
number of bands to the third power.. Crucially, we take into account the
renormalization of the single-particle dispersion as a result of interactions
[11, 59, 27], as described in Sec. III.2. As we shall see, this effect can
significantly impact the nature of the ground state.

Our key observation is that the continuum model for TBG at charge neutrality
is nearly free of the fermion sign problem in quantum Monte Carlo, due to its
approximate particle-hole symmetry (PHS). Under a slight deformation of the
model – neglecting the rotation of the Pauli matrices in the individual Dirac
Hamiltonians of the two graphene sheets by the twist angle $\theta$ – this
particle-hole symmetry becomes exact, and the model is completely sign-problem
free. Given the smallness of the magic angle $\theta\approx 1.05^{\circ}$, the
magnitude of the neglected term is $\sim 2\%$ of the total Hamiltonian. We
note that the PHS enhances the degeneracy of different symmetry-breaking
states in the flat-band limit only, while finite dispersion lifts this
degeneracy[27].

The resulting PHS is represented as an anti-unitary operator acting on the
many-body Hilbert space. In addition to transforming an electron into a hole,
it also flips the valley, sublattice, and layer indices [see Eq. (2)]. This
PHS, combined with a $U(1)$ spin or valley conservation symmetry, allows for a
sign-problem free formulation of the DQMC method (Sec. IV). As a benchmark, we
simulated the system in the interaction-only case where the single-particle
dispersion is artificially turned off. In this case, our simulations reproduce
the exact results for the ground state energy, the nature of the ground state
manifold [27], and the single-particle excitation spectrum.

We then turn to the realistic case, where the single-particle dispersion is
included. Our main results are summarized in Fig. 1(c–j). One of key
parameters in the BM model, which changes both the bandwidth and the
wavefunctions of the active bands, is the ratio between the intra-sublattice
(A–A) and inter-sublattice (A–B) tunneling matrix elements between the two
layers. We denote these matrix elements by $w_{0}$ and $w_{1}$, respectively.
The ratio $\kappa=w_{0}/w_{1}$ is expected to be less than unity due to
lattice relaxation [81, 23]. We have performed DQMC simulations for two values
of $\kappa$, $0.75$ and $0.3$. For $\kappa=0.75$, we find that the ground
state is a gapped Kramers inter-valley coherent (K-IVC) state [82],
characterized by a spontaneous hybridization between the two valleys and
preserving an effective time-reversal symmetry. The correlation function of
the K-IVC order parameter as a function of temperature is shown in Fig. 1(c)
for three system sizes, showing a rapid growth of the correlation length below
$T\sim 10\,\textrm{K}$.

Additional insight into the nature of the K-IVC state is gained by calculating
the single-particle excitation spectrum, using the maximum entropy method
[83]. Fig. 1(d) shows the single-particle density of states, $\rho(\omega)$,
summed over both bands and spin/valley flavors, at different temperatures.
$\rho(\omega)$ is closely related to the spectrum measured in scanning
tunneling spectroscopy (STM) experiments. At the lowest temperature,
$T=1.8\,\textrm{K}$, the spectrum shows a clear gap of about
$10\,\textrm{meV}$ for either particle or hole excitations. Surprisingly,
however, a gap-like structure in $\rho(\omega)$ is apparent already at
temperatures as high as $T\approx 60\,\textrm{K}$, where the correlation
length of the K-IVC order is around one moiré lattice spacing. We refer to
this behavior as the formation of a ‘pseudogap’, i.e., a spectral gap that
appears with no accompanying long-range order.

Panels (e),(f) show the spectral function $A(\bm{k},\omega)$ for $\kappa=0.75$
along a cut in the Brillouin zone at two representative temperatures,
$T=1.8\,\textrm{K}$ and $29\,\textrm{K}$. The dispersion at $1.8\,\textrm{K}$
is remarkably close to the spectrum of the single-particle excitation of the
interaction-only model where $\hat{H}_{0}$ is neglected (blue crosses), for
which exact results are available. The agreement between the interaction-only
and DQMC dispersions is surprising, since the bare dispersion bandwidth $W$ is
significant, $W\approx 8\,\textrm{meV}$ for $\kappa=0.75$ (the single-particle
dispersion used in the DQMC simulation is indicated by the red crosses). It is
also worth noting that the minimum of the dispersion is at the center of the
moiré Brillouin zone (the $\Gamma_{M}$ point) [27], in contrast to the single-
particle dispersion which has Dirac points at corners of the moiré Brillouin
zone, $K_{M}$ and $K^{\prime}_{M}$. At the higher temperature,
$T=29\,\textrm{K}$, the peaks in the dispersion are broader 222This broadening
may be an artifact of the maximum entropy algorithm, whose intrinsic frequency
resolution is proportional to temperature, or may represent the physical
lifetime of the excitations at finite $T$., but the main features – including
the gap at the $\Gamma_{M}$ point – are still clearly visible, despite the
absence of long-range order. Note, the spectral function can be accessed by
ARPES measurements, the measurement of which has been initiated in recent
experiments [39, 40].

The behavior for $\kappa=0.3$ in Fig. 1(g,h) is dramatically different from
that of $\kappa=0.75$. In this case, the correlation length of the K-IVC order
parameter remains short and nearly system size-independent down to the lowest
temperatures. This result is corroborated for $T=0$ using a projective variant
of the QMC algorithm [70, 71] as visualized by the open boxes in Fig. 1(g).
The system is semi-metallic with Dirac points at $K_{M}$ and $K^{\prime}_{M}$,
as can be seen from the spectral function, $A(\bm{k},\omega)$, at
$T=1.8\,\textrm{K}$ (panel j). The small gap of about $0.7\,\textrm{meV}$ at
the $K_{M}$, $K^{\prime}_{M}$ points is due to the finite-size charging
energy; its value decreases as $1/L^{2}$ with the system size. The spectral
function at $T=29\,\textrm{K}$ is much broader, although the main features of
the dispersion are already visible.

The fact that the K-IVC state is suppressed with decreasing $\kappa$ may seem
surprising, since in the BM model, the bandwidth vanishes at the magic angle
for $\kappa=0$ (the ‘chiral limit’ [85]). However, we find that the bare
dispersion bandwdith increases with decreasing $\kappa$ [Fig. 1(b)] at
$\theta=1.05^{\circ}$. In particular, for $\kappa=0.3$, the total width of the
single-particle dispersion [panels (i),(j); red crosses] is $W\approx 14$ meV,
almost twice as large as that of $\kappa=0.75$, explaining why in the former
case the ground state is a semi-metal, whereas in the latter we obtain a
gapped K-IVC state.

The semi-metal state at $\kappa=0.3$ is still strongly correlated, however;
only about half of the single-particle spectral weight is contained in the
coherent quasi-particle dispersion [the sharp dispersing peaks in panel (j)].
The rest of the spectral weight is distributed over a wide energy range. At
$K_{M}$ and $K^{\prime}_{M}$ the spectral function exhibits broad peaks near
the excitation energies of the interaction-only model [panel (j), blue
crosses]. This spectral transfer is clearly seen in the total single-particle
density of states, $\rho(\omega)$, shown in panel (h). At $T=58\,\textrm{K}$,
$\rho(\omega)$ exhibits two maxima centered around $\omega=\pm
12\,\textrm{meV}$. Below $T\sim 30\,\textrm{K}$, these evolve into the broad
‘shoulders’ around $\omega=\pm 15\,\textrm{meV}$. The shoulders are distinct
from the sharp maxima associated with the semi-metallic dispersion at
$|\omega|\leq 10\,\textrm{meV}$.

We note the calculated $\rho(\omega)$ can be accessed in STM experiments, on
averaging spectra over the unit cell. The broad similarity between the
experimentally observed spectra near charge neutrality in low temperature STM
measurements eg. Ref. [78] and Figure 1(d) should be noted. More
significantly, we further predict that the suppression of low energy spectral
weight persists upto temperatures of $T\leq 10K$, beyond which it gradually
fills in, which can be probed in future experiments.

##  III Twisted Bilayer Graphene

Before describing the QMC simulations, we first introduce the single-particle
and many-body Hamiltonians for twisted bilayer graphene.

###  III.1 Single Particle Hamiltonian

In an appropriate basis, the continuum model single-particle Hamiltonian of
TBG [56] can be written as follows (See Appendix B):

$\displaystyle\hat{H}_{\textrm{BM}}$
$\displaystyle=\sum_{{\bm{k}}}c^{\dagger}_{\bm{k}}H_{\bm{k}}c_{\bm{k}}+\sum_{{\bm{k}},{\bm{k}}^{\prime}}c^{\dagger}_{{\bm{k}}^{\prime}}T_{{\bm{k}}^{\prime},{\bm{k}}}c^{\,}_{\bm{k}}$
$\displaystyle H_{\bm{k}}$ $\displaystyle=\hbar
v_{F}\quantity(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y})$ $\displaystyle
T_{{\bm{k}}^{\prime},{\bm{k}}}$
$\displaystyle=\frac{1}{2}\sum_{n}\delta_{{\bm{k}}^{\prime},{\bm{k}}-\tau_{z}{\bm{q}}_{n}}\quantity(\mu_{x}-i\mu_{y})T_{n}+\textrm{h.c.}$
$\displaystyle T_{n}$
$\displaystyle=w_{0}e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}+w_{1}e^{2\pi
in\sigma_{z}\tau_{z}/3}\sigma_{x}e^{-2\pi in\sigma_{z}\tau_{z}/3}.$ (1)

where $\theta$ is a twist angle, $v_{F}$ is the Dirac velocity, and $w_{0,1}$
is interlayer moiré hopping term. $c^{\dagger}_{\bm{k}}$ represents a row
vector of electron creation operators with wavevector ${\bm{k}}$ and with
sublattice, valley, layer, and spin indices, denoted by $(\sigma,\tau,\mu,s)$,
respectively. Here ${\bm{q}}_{n}={\cal R}_{2\pi n/3}\,{\bm{q}}_{0}$, where
${\cal R}_{\varphi}=e^{i\varphi\sigma_{y}}$ is a rotation matrix,
${\bm{q}}_{0}=\frac{8\pi\sin(\theta/2)}{3\sqrt{3}a_{0}}(0,-1)^{T}$, and
$a_{0}=1.42\textrm{ \AA}$ is the intra-layer distance between carbon atoms.

In the basis used here, the twist-angle dependence only appears in the $w_{0}$
term and in the moiré wavevector $\bm{q}_{n}$. Therefore, in the chiral limit
[85] $w_{0}\rightarrow 0$, the explicit twist-angle dependency in Eq. (III.1)
disappears. In this case, the Hamiltonian is invariant under the following
many-body anti-unitary particle-hole transformation:

$\displaystyle{\cal C}:\quad c^{\dagger}_{\bm{k}}$
$\displaystyle\mapsto\tau_{x}\sigma_{x}\mu_{y}c^{\,}_{\bm{k}}$ $\displaystyle
i$ $\displaystyle\mapsto-i$ (2)

such that, for $w_{0}=0$,
${\cal{C}}\hat{H}_{\textrm{BM}}{\cal{C}}^{-1}=\hat{H}_{\textrm{BM}}$. Note
that if we denote the action of ${\cal{C}}$ on the flavor indices as
$U_{\cal{C}}=\tau_{x}\sigma_{x}\mu_{y}$, the first-quantized Hamiltonian
satisfies $U_{\cal{C}}H_{\bm{k}}U_{\cal{C}}^{-1}=-H_{\bm{k}}$ and
$U_{\cal{C}}T_{{\bm{k}},{\bm{k}}^{\prime}}U_{\cal{C}}^{-1}=-T_{{\bm{k}},{\bm{k}}^{\prime}}$.

As we will demonstrate in Sec. IV, the presence of an anti-unitary particle-
hole symmetry is important for the absence of the sign problem in DQMC.
However, away from the chiral limit, $w_{0}\neq 0$, ${\cal{C}}$ is weakly
broken. The following _small angle approximation_ promotes the approximate PHS
back to an exact symmetry by neglecting the rotation of the Pauli matrices in
the individual Dirac Hamiltonians of the two graphene sheets,

$w_{0}e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}\mapsto w_{0}\,.$ (3)

Since $\theta/\pi\ll 1$ near the magic angle, the wavefunctions and the
combined bandwidth for conduction and valence bands barely change for
$\kappa<0.8$, while the dispersion is now perfectly particle-hole symmetric
unlike the original BM model Hamiltonian. (See Appendix B). As
$\kappa\rightarrow 1$, the gap between the conduction band and the remote band
starts closing and hybridizes, decreasing the overlap between flat band
subspaces of the original BM model and approximated Hamiltonians. At this
point, the idea of projecting the system into flat bands would not be valid
anymore, and one has to include further neighboring bands to properly capture
the hybridzation effect. However, in the experimentally relevant regime
$\kappa$ is considered to be smaller than $0.8$ according to the first-
principle calculations [23]. Therefore, the small angle approximation is
justified, and will be used henceforth.

The resulting band structure with the correction from the interaction is shown
in Fig. 1(a) and characterized by a pair of nearly flat bands at the Fermi
level, while all other valence and conduction bands are separated by a large
gap of order $15\,\textrm{meV}$. Thus, the remote bands can be considered
inert at low temperatures where valence (conduction) bands are completely
filled (empty). Foreshadowing the Coulomb interaction introduced in the next
section, we project the charge-density operator, with respect to the
background at charge neutrality, onto the nearly flat bands,

$\bar{\rho}_{\bm{q}}=\sum_{{\bm{k}}}d^{\dagger}_{{\bm{k}}+{\bm{q}}}\Lambda({\bm{k}},{\bm{q}})d^{\,}_{{\bm{k}}}-\frac{1}{2}\sum_{{\bm{k}}}\delta_{{\bm{q}},{\bm{G}}}\tr\Lambda({\bm{k}},{\bm{q}})\,,$
(4)

where $d^{\dagger}_{{\bm{k}},i}$ is a creation operator for a Bloch
eigenstate, $\ket{u_{{\bm{k}},i}}$, of flavour $i=(s,\tau,n)$, i.e., spin $s$,
valley $\tau$, and band $n$. $d^{\dagger}_{{\bm{k}}}$ is a row vector
containing the eight operators $d^{\dagger}_{{\bm{k}},i}$. The momentum
${\bm{k}}$ is restricted to the $1^{\mathrm{st}}$ moiré BZ while the transfer
momentum ${\bm{q}}$ remains unrestricted. The form factor is defined as
$\Lambda_{ij}({\bm{k}},{\bm{q}})\equiv\innerproduct{u_{{\bm{k}}+{\bm{q}},i}}{u_{{\bm{k}},j}}$.

The action the symmetry operators on $d_{{\bm{k}}}$ depends on the choice of
gauge. Particularly important for our purpose is the action of ${\cal{C}}$,
responsible for the absence of the sign problem. In Appendix D, we describe a
gauge fixing procedure such that

${\cal C}:\quad d^{\dagger}_{\bm{k}}\mapsto\tau_{x}n_{x}d^{\,}_{\bm{k}},\qquad
i\mapsto-i$ (5)

where $n_{x,y,z}$ are Pauli matrices acting on the band index.

###  III.2 Many-Body Hamiltonian

The parameters typically used for the BM model stem from first-principle
calculations that include interaction effects, for example, the Dirac velocity
is already renormalized. This implies that the QMC simulation has to avoid
double counting [59, 27]. One approach to this subtle problem is to subtract
the Hartree-Fock mean-field contribution of a properly chosen reference state,
whose corresponding single-particle density matrix is given by $P_{0}$, as the
following:

$\displaystyle\hat{H}$ $\displaystyle=\hat{H}_{0}+\hat{V}$
$\displaystyle\hat{H}_{0}$
$\displaystyle=\hat{H}_{\textrm{BM}}-\big{[}\hat{V}\,\big{]}_{P_{0}}$
$\displaystyle\hat{V}$
$\displaystyle=\frac{1}{2}\sum_{\bm{q}}\bar{\rho}_{\bm{q}}V_{\bm{q}}\bar{\rho}_{-{\bm{q}}},$
(6)

where
$V_{\bm{q}}=V_{0}\tanh(\absolutevalue{{\bm{q}}}d)/\absolutevalue{{\bm{q}}}$
with $V_{0}=\frac{e^{2}}{2\epsilon\epsilon_{0}}\frac{1}{A}$ is the Fourier-
transformed Coulomb potential with gate distance $d$ and relative permittivity
$\epsilon$. Here, we use the density matrix with respect to the reference
decoupled single graphene sheets, as discussed in Appendix C. This is an
intuitive choice since the first principle calculation were preformed for a
monolayer graphene, however, this choice is not unique. In fact, one may also
use the non-interacting ground state of the _coupled_ system as the reference
state, and the resulting bandwidth can vary by a factor of $2$, apparent by
comparing [59] and [27]. Therefore, one should view the resulting corrected
bare dispersion as an educated guess. Note that a particle-hole symmetric
reference state $P_{0}$ induces a bare Hamiltonian, $\hat{H}_{0}$, that
preserves the PHS ${\cal{C}}$ and thus will not introduce a sign-problem.

##  IV DQMC method for twisted bilayer graphene

We now set up our problem for DQMC simulations, and show that the sign problem
is absent at charge neutrality due to the PHS (5).

###  IV.1 Momentum-space algorithm and absence of the sign problem

The algorithm we use to simulate the Hamiltonian (III.2) is a variant of the
Blankenbecler–Scalapino–Sugar algorithm [86], adapted to work in momentum
space [87]. We first use a Trotter decomposition to write the partition
function as

$Z={\rm
Tr}\left(e^{-\hat{H}_{0}\delta}\prod_{\bm{q},\alpha=1,2}e^{-\frac{V_{\bm{q}}}{8}\hat{\rho}_{\bm{q},\alpha}^{2}\delta}\right)^{N_{\tau}}+O(\delta^{2}).$
(7)

Here, the trace is taken over the many-body Hilbert space of fermions in the
active bands, $N_{\tau}$ is the number of imaginary time steps,
$\delta=\beta/N_{\tau}$ is the step size, and we have introduce the real and
imaginary parts of the density operator:

$\displaystyle\hat{\rho}_{\bm{q},1}$
$\displaystyle=\bar{\rho}_{\bm{q}}+\bar{\rho}_{-\bm{q}},$
$\displaystyle\hat{\rho}_{\bm{q},2}$
$\displaystyle=\frac{1}{i}\left(\bar{\rho}_{\bm{q}}-\bar{\rho}_{-\bm{q}}\right).$
(8)

A subtlety that arises in Eq. (7) is the ordering of the terms in the product
over $\bm{q}$ and $\alpha$, since the $\hat{\rho}_{\bm{q},\alpha}$ operators
do not commute with each other (due to the projection to the active bands). A
symmetric ordering can be chosen in order to achieve a Trotter error that
scales as $\delta^{2}$, as described in Appendix F. For notational simplicity,
we leave the ordering of the product implicit here.

We simulate a finite lattice of $L\times L$ moiré unit cells. Both crystal
momentum ${\bm{k}}$ and momentum transfer ${\bm{q}}$ in Eq. (7) are
discretized as $(m_{1}\bm{G}_{1}+m_{2}\bm{G}_{2})/L$, where $G_{1,2}$ are
moiré reciprocal lattice vectors and $m_{1,2}\in\mathbb{Z}$. Unlike
${\bm{k}}$, ${\bm{q}}$ is not restricted to the first Brillouin zone, and in
principle has infinite range. However, for any sensible physical systems the
magnitude of $V_{\bm{q}}\bar{\rho}_{\bm{q}}\bar{\rho}_{-{\bm{q}}}$ decreases
rapidly with ${\bm{q}}$ as $\Lambda({\bm{k}},{\bm{q}})$ decays rapidly with
$\absolutevalue{{\bm{q}}}$ 333The main reason for the decay is the decay of
the form factors contained inside $\bar{\rho}_{\bm{q}}$ from the flat band
projection., and one can set a cutoff for the range of ${\bm{q}}$. In our
simulations, we restricted our ${\bm{q}}$ to be within twice the size of the
first Brillouin zone in each direction.

To facilitate DQMC simulations, we perform a Hubbard-Stratonovich (HS)
transformation to decouple the quartic terms in (7), writing the partition
function as

$Z=\sum_{\\{\eta_{\bm{q},\alpha}(l)\\}}W(\\{\eta_{\bm{q},\alpha}(l)\\})+O(\delta^{2}),$
(9)

where $\eta_{\bm{q},\alpha}(l)$ are real, discrete HS fields, and the weight
$W(\\{\eta_{\bm{q},\alpha}(l)\\})$ is given by

$W={\rm
Tr}\left[\prod_{l=1}^{N_{\tau}}\left(e^{-\hat{H}_{0}\delta}\prod_{\bm{q},\alpha}e^{-\hat{H}_{\rm{i}}(\bm{q},\alpha,l)\delta}f\left[\eta_{\bm{q},\alpha}(l)\right]\right)\right].$
(10)

Here,
$\hat{H}_{\rm{i}}(\bm{q},\alpha,l)=i\left(\frac{V_{\bm{q}}}{8}\delta\right)^{1/2}\eta_{\bm{q},\alpha}(l)\hat{\rho}_{\bm{q},\alpha}$,
$f[\eta_{\bm{q},\alpha}(l)]$ is a real, non-negative weight of the HS field,
given explicitly in Eq. (74).

Proving that the sign problem is absent amounts to showing that
$W(\\{\eta_{\bm{q},\alpha}(l)\\})$ is a real and non-negative for any
configuration of the HS fields. To show this, we note that due to the $U(1)$
spin and valley conservation symmetries, we can write
$\hat{H}_{0}=\sum_{s,\tau}\hat{H}_{0,s,\tau}$, where $\hat{H}_{0,s,\tau}$ is
the part of $\hat{H}_{0}$ that acts on electrons with spin $s$ and valley
$\tau$, and similarly $\hat{H}_{\rm{i}}=\sum_{s,\tau}\hat{H}_{\rm{i},s,\tau}$.
Using this fact, we can decompose $W$ as a product:

$W(\\{\eta_{\bm{q},\alpha}(l)\\})=\prod_{s,\tau}W_{s,\tau}(\\{\eta_{\bm{q},\alpha}(l)\\}),$
(11)

where $W_{s,\tau}$ is of the form of Eq. (10), replacing $\hat{H}_{0}$ by
$\hat{H}_{0,s,\tau}$ (and similarly for $\hat{H}_{\rm{i}}$), and taking the
trace over electrons with spin and valley indices $(s,\tau)$.

We now use the symmetry of the Hamiltonian under $\mathcal{C}$ [Eq. (5)]. In
particular,
$\mathcal{C}\hat{H}_{0,s,\tau}\mathcal{C}^{-1}=\hat{H}_{0,s,-\tau}$, and
similarly for $\hat{H}_{\rm{i},s,\tau}$. Using this transformation law, and
accounting for the anti-unitarity of $\mathcal{C}$, one can show that (see
Appendix A, Eq. (23))

$W_{s,\tau}(\\{\eta_{\bm{q},\alpha}(l)\\})=W_{s,-\tau}^{*}(\\{\eta_{\bm{q},\alpha}(l)\\}).$
(12)

Hence, we find that
$W(\\{\eta_{\bm{q},\alpha}(l)\\})=\prod_{s}\left|W_{s,\tau}\right|^{2}\geq 0$,
and there is no sign problem. As a corollary, we note that the weight of a
single spin flavor is real and non-negative on its own. Therefore, a spin-
polarized system with a single spin flavor is also sign problem free.

The DQMC algorithm proceeds by interpreting $W/Z$ as the probability of the
configuration $\\{\eta_{\bm{q},\alpha}(l)\\}$, and sample the configuration
space using a Markov Chain Monte Carlo method. Since both $\hat{H}_{0}$ and
$\hat{H}_{\rm{i}}$ are quadratic in the fermion operators, the fermionic trace
can be performed exactly for a given configuration of the HS fields (Appendix
F).

| Chiral | Nonchiral | Intervalley
---|---|---|---
$\matrixquantity{\textrm{Required}\\\ \textrm{Multiplicity}}$ | $\matrixquantity{\textrm{2 valleys}\\\ \textrm{OR 2 spins}}$ | 2 valleys | $\matrixquantity{\textrm{4 valleys}\\\ \textrm{AND spins}}$
Table 1:  The minimal number of bands required to guarantee the absence of
sign problem in DQMC. We always assume a conduction/valence band pair, and
note the multiplicities of such pairs to obtain a sign problem free
formulation. For the nonchiral case ($\kappa>0$), we always need at least two
conduction-valence band pairs in opposite valleys, since it is the valley-
exchanging anti-unitary symmetry that guarantees the absence of sign problem.
However, with chiral symmetry ($\kappa=0$), even in the presence of band
dispersion, a single valley is sign problem free, provided both spin-species
are present. Finally, we note that if we include intervalley density-density
interactions, that are typically suppressed at small angle, the HS-decoupled
Hamiltonian is not block-diagonal in valley anymore and we need all four pairs
of bands (both spin and valley) to guarantee the absence of sign problem.

###  IV.2 Sign-problem-free extensions

Employing the identity Eq. (23), we now consider various extensions and
physical perturbations to the TBG Hamiltonian (III.2) that allow for sign
problem-free DQMC simulations. Table 1 summarizes different cases of the TBG
where the sign problem is absent, labeled by whether we are working in the
chiral limit ($\kappa=0$) or not, and whether an inter-valley scattering
interaction is included [see Eq. (13)]. For each case, we specify the number
of conserved flavors (valley and spin) needed to avoid a sign problem.

As we have seen, in the most generic case (away from the chiral limit,
$\kappa>0$, and with $\hat{H}_{0}\neq 0$), the DQMC weight associated with
each spin flavor is real and non-negative, and therefore a problem with a
single spin flavor is sign problem-free [Eq. (12) and the discussion below].
What if we have more symmetries? In the chiral limit ($\kappa=0$), we can
relax the condition for the sign problem further. In this case, we have an
extra anti-unitary particle-hole symmetry ${\cal{C}}_{\textrm{sub}}$ whose
unitary part is simply given as $U_{{\cal{C}}_{\textrm{sub}}}=\tau_{z}n_{x}$
in the band basis (See Tab. 3). The origin of ${\cal{C}}_{\textrm{sub}}$ is
the sublattice symmetry of the original graphene layer. Since it maps each
part of the HS-decoupled Hamiltonian with a given $(s,\tau)$ onto itself,
$W_{s,\tau}$ must be real. For a single valley, if we have both spin flavors,
then the presence of the unitary spin-flip symmetry $s_{x}$ implies that
$W_{s,\tau}=W_{-s,\tau}$ due to Eq. (23). Therefore, $W_{\tau}\equiv
W_{\uparrow,\tau}W_{\downarrow,\tau}\geq 0$. In other words, one can simulate
a single valley of TBG without a sign problem in this case.

In the special case where the bare dispersion $\hat{H}_{0}$ is set to zero in
the chiral limit, we can show that only two bands related by ${\cal{C}}$
symmetry are needed to guarantee the absence of sign problem. Here we can
choose a basis which diagonalizes $n_{x}$, denoted by $\tilde{n}$ (Appendix
D). Then, HS-decoupled Hamiltonian is diagonal in $(s,\tau,\tilde{n})$.
Furthermore, the extra anti-unitary particle-hole symmetry
${\cal{C}}_{\textrm{sub}}$ maps each flavor to itself since it commuts with
$n_{x}$. It implies that the corresponding determinant $M_{s,\tau,\tilde{n}}$
is real. Now, note that the unitary part of ${\cal{C}}$ in this basis is given
as $U_{{\cal{C}}}=\tau_{y}\sigma_{y}$. It implies that
$W_{s,\tau,\tilde{n}}=W_{s,-\tau,-\tilde{n}}^{*}$. Therefore, $W_{s,C}\equiv
W_{s,\tau,\tilde{n}}W_{s,-\tau,-\tilde{n}}$ is real and non-negative, where
$C=\tau_{z}\sigma_{z}$ is the Chern number for this pair of flavors.

Finally, it is worth noting that one can add an inter-valley density
interaction. Such an interaction was not considered in our current
simulations, due to the fact that it is estimated to be smaller than the long-
range Coulomb interaction by a factor of $a_{0}/a_{M}\sim\sin(\theta)$, i.e.
the ratio of the atomic lattice constant $a_{0}$ to the moiré lattice
constant. $a_{M}$. However, this interaction is allowed by symmetry, and
despite its small magnitude it can have significant effects in lifting
degeneracy. The inter-valley interaction is of the form
$\hat{H}_{\rm{iv}}\propto\sum_{\bm{q}}\rho_{\bm{q}}^{\textrm{inter}}\rho_{-{\bm{q}}}^{\textrm{inter}}$,
where the inter-valley density is defined as

$\rho_{\bm{q}}^{\textrm{inter}}=\sum_{\bm{k}}d_{{\bm{k}}+{\bm{q}}}^{\dagger}\Lambda^{\textrm{inter}}({\bm{k}},{\bm{q}})d_{{\bm{k}}}^{\,}.$
(13)

$\Lambda^{\textrm{inter}}({\bm{k}},{\bm{q}})$ is off-diagonal in valley space.
This interaction reduces the manifold for the K-IVC state from U(2) to U(1)
for spin singlet K-IVC case, which allows for a finite-temperature transition.
Once this interaction is introduced, the HS-decoupled Hamiltonian is not
block-diagonal in valley anymore. As a result, using ${\cal{C}}$ symmetry, we
can show that $W_{s}$ is real, but not necessarily non-negative. By further
using spin-flip unitary symmetry $s_{x}$, we can show that
$W_{\uparrow}W_{\downarrow}>0$. Therefore, we need all physical flavors of TBG
to be present in this case to evade the sign problem.

Using our approach, we can further show that the model remains sign-problem
free in the presence of external fields 444First quantized matrix for such
operator should have zero trace, unless it is chemical potential. First,
consider the Zeeman field term

$\hat{H}_{\textrm{Z}}=-h_{z}\sum_{\bm{k}}d^{\dagger}_{\bm{k}}s_{z}d^{\,}_{\bm{k}}.$
(14)

Let $\hat{H}_{s}$ be a HS-decoupled Hamiltonian for spin flavor $s$,
containing both valleys and bands. Then, we can show that
$(s_{x}{\cal{C}})\hat{H}_{s}(s_{x}{\cal{C}})=\hat{H}_{-s}$. Therefore, using
anti-unitary particle-hole symmetry $s_{x}{\cal{C}}$ we obtain that
$W_{\uparrow}^{*}=W_{\downarrow}$ and
$W=W_{\uparrow}W_{\downarrow}=|W_{\uparrow}|^{2}>0$.

We can also include a ‘valley-Zeeman’ field,

$\hat{H}_{\textrm{vZ}}=-h_{vz}\sum_{\bm{k}}d^{\dagger}_{\bm{k}}\tau_{z}d^{\,}_{\bm{k}}.$
(15)

This term is symmetric under ${\cal{C}}$, and does not cause a sign problem
even if we have a single spin flavor. Similarly, we can add the uniform
sublattice polarization ($\sim\tau_{z}n_{x}$ in band basis) and perpendicular
electric field terms, which are also symmetric under ${\cal{C}}$. Therefore,
in principle we can turn on all these four kinds of external fields and
simulate the system without a sign problem.

##  V Simulation Results and Discussion

In this section we report the results of our QMC computations. Our
implementation of the algorithm described above is built on top of the
software package ALF [90, 91]. We investigate the leading instabilities of the
BM model at charge neutrality, and benchmark our simulations by comparing the
results to the analytical solutions in the strong coupling limit [27].
Furthermore, we present a more detailed analysis of the spectral functions
shown in Fig. 1 focusing on their temperature dependence.

###  V.1 Parameters

We used to following set of parameters in the simulations. In the BM model, we
used $\frac{\hbar v_{F}}{3a_{0}/2}=2700\,\textrm{meV}$ and
$w_{1}=105\,\textrm{meV}$. We set $\theta=1.05$ which is close to the magic
angle for the above parameters. For the interaction, we focus on a relative
dielectric constant $\epsilon=10$ and a gate distance of $d=20\,\mathrm{nm}$.
We confirmed the convergence of QMC results for a Trotter step size of down to
$\delta=0.01\,\textrm{meV}^{-1}$ as shown in Fig. 9.

Effectively, two most important parameters in the simulation is the bandwidth
$W$ of bare dispersion $\hat{H}_{0}$ and the interaction energy scale. It is
worth noting that $W$ is a complicated function of the system parameters. Due
to the correction term in Eq. (III.2), it depends on the interaction
parameters, $\epsilon$ and $d$, as well as on the choice of reference state
$P_{0}$ and the hopping ratio $\kappa=w_{0}/w_{1}$. The interaction energy
scale can be represented by the bandwidth of electron and hole excitations in
the strong coupling limit, $\hat{H}_{0}=0$ (Appendix C). In Fig. 1(b), we plot
$W$ and the interaction scale as a function of $\kappa$ at
$\theta=1.05^{\circ}$. In this plot, $W$ decreases with $\kappa$, while the
interaction scale barely changes. However, this behavior is not universal; as
we move away from the magic angle, $W$ can increase with $\kappa$. (Fig. 7).

###  V.2 Observables

In our DQMC simulation, there are two observables we examine: ($i$) the
electron Green’s function $G({\bm{k}},\tau)$, and ($ii$) correlation functions
$S_{a}$ of fermion bilinears ${\cal O}_{a}$,

$\displaystyle S_{a}({\bm{q}})$
$\displaystyle\equiv\frac{1}{L^{2}}\expectationvalue{{\cal
O}_{a}({\bm{q}}){\cal O}_{a}(-{\bm{q}})}$ $\displaystyle{\cal
O}_{a}({\bm{q}})$
$\displaystyle\equiv\sum_{\bm{k}}d^{\dagger}_{{\bm{k}}+{\bm{q}}}\Lambda_{a}({\bm{k}},{\bm{q}})d^{\,}_{\bm{k}}.$
(16)

$N_{M}=L^{2}$ is the number of moiré unit cells in the system and
$\Lambda_{a}$ is the operator-dependent form factor,

$\quantity[\Lambda_{a}({\bm{k}},{\bm{q}})]_{ij}\equiv\bra{u_{{\bm{k}}+{\bm{q}},i}}M_{a}\ket{u_{{\bm{k}},j}},$
(17)

where $i,j$ runs over the flavor indices $(s,\tau,n)$ of the projected flat
bands. The microscopic operator $M_{a}$ for various correlation functions are
defined in Tab. 2, which act on the vector space of the unprojected electron
operator. The normalization of (V.2) implies that a spatially homogeneous
long-range order gives rise to an extensive scaling of $S_{a}({\bm{q}}=0)\sim
L^{2}$. Note that the HS decoupled Hamiltonian is a fermion bilinear such that
higher correlation function can be extracted using Wick’s theorem, as
illustrated in Appendix D 4\. The correlation length, $\xi_{a}$, which we
presented in Fig. 1, is extracted from the momentum dependence of
$S_{a}({\bm{q}})$ at small wavevectors [92],

$\xi^{2}_{a}/a^{2}_{M}=\frac{3}{16\sin^{2}(\pi/L)}\left(\frac{S_{a}({\bm{q}}=0)}{S_{a}({\bm{q}}={\bm{q}}_{nn})}-1\right),$
(18)

where ${\bm{q}}_{nn}$ are the six nearest-neighboring momenta of ${\bm{q}}=0$,
$S_{a}({\bm{q}}={\bm{q}}_{nn})$ is averaged over them.

At each momentum ${\bm{k}}$ we have eight electron flavors, and therefore
there is a large manifold for potential symmetry breaking. In order to examine
the nature of the ground state manifold, we choose six representative order
parameters: valley polarized (VP), spin polarized (SP), valley-Hall (VH), two
types of intervalley coherence (t-IVC/K-IVC), and quantum Hall (QH), whose
correspodning microscopic operators $M_{a}$ are listed in Tab. 2.

Figure 2:  Correlation functions at ${\bm{q}}=0$ and global symmetry $G$ of full-interacting TBG Hamiltonians in four different limits ($\epsilon=10$ and $d_{s}=20$nm). Each limit is characterized by: $(\kappa,W)$ where $\kappa=w_{0}/w_{1}$ and $W$ is the bandwidth of the bare dispersion $\hat{H}_{0}$ in Eq. (III.2). (a) chiral, flat limit: $(0,0\,\textrm{meV})$, (b) chiral, non-flat with artifically reduced bandwidth: $(0,3.8\,\textrm{meV})$, (c) non-chiral, flat: $(0.75,0\,\textrm{meV})$, (d) non-chiral, non-flat: $(0.75,8.3\,\textrm{meV})$, and (e) non-chiral, non-flat: $(0.3,14.3\,\textrm{meV})$. Correlation functions $S_{i}({\bm{q}}=0)$ for K-IVC, VP, VH and QH are shown, which should scale as the system volume if they acquire long-range order. In (a), all channels have equal magnitude of correlation functions. The suppression of $S_{\textrm{QH}}({\bm{q}}=0)$ in this case is due to the difference in the degeneracy of $C_{\textrm{tot}}=0$ and $C_{\textrm{tot}}=2$ groundstate manifold sectors. The deviations from this ideal limit of $(\kappa,W)=(0,0\,\textrm{meV})$ lift the degeneracy among symmetry broken states, as tabulated in Tab. 2. In (b) where we retain a quarter of the bare dispersion to obtain an ordered state, the $U(4)_{\rm chiral}$ relates the K-IVC and VH orders, while the VP is suppressed due to dispersion. In (c), compared to the chiral nonflat case, the system now has a different $\mathrm{U}(4)_{\textrm{flat}}$ symmetry which relates K-IVC and VP states instead of VH. Here, VH is suppressed due to nonchiral part of the interaction. Progressing from (a) to (d), each perturbation reduces the symmetry down to that in the realistic limit, $\mathrm{U}(2)_{\bm{K}}\times\mathrm{U}(2)_{{\bm{K}}^{\prime}}$, where only the K-IVC groundstate survives. In (e), we are in the same symmetry class as in (d), but a larger bandwidth. As a result, the correlation functions saturate to a finite value with increasing system size, resulting in a disordered phase. In (e), although not shown, we emphasize that the magnitudes of both VP and SP correlations are vanishing at low temperature.  | $\matrixquantity{\kappa=0\\\ \,W=0\,}$ | $\matrixquantity{\kappa=0\\\ \,W>0\,}$ | $\matrixquantity{\kappa>0\\\ \,W=0\,}$ | $\matrixquantity{\kappa>0\\\ \,W>0\,}$ | $M_{i}$
---|---|---|---|---|---
VP | ✓ | ✗ | ✓ | ✗ | $\tau_{z}$
SP | ✓ | ✗ | ✓ | ✗ | $s_{z}$
QH | ✓ | ✓ | ✗ | ✗ | $\tau_{z}{\sigma}_{z}$
VH | ✓ | ✓ | ✗ | ✗ | ${\sigma}_{z}$
t-IVC | ✓ | ✗ | ✗ | ✗ | $\tau_{x,y}\sigma_{y}\mu_{y}$
K-IVC | ✓ | ✓ | ✓ | ✓ | $\tau_{x,y}\sigma_{x}\mu_{y}$
Table 2:  Order parameters and regimes where the corresponding symmetry broken
phases are favored for the parameters $(\kappa,W)$ in the strong coupling
limit [27], i.e., small dispersion limit. In the last column, microscopic
operators are listed for each order parameter. In the chiral limit with
dispersionless (flat) band, i.e., $(\kappa,W)=(0,0)$, the system has
$\mathrm{U}(4)\times\mathrm{U}(4)$ global symmetry [27]. Accordingly, symmetry
broken states for the listed order parameters become degenerate. For $W>0$,
the finite dispersion gives rise to an anti-ferromagnetic coupling between
different Chern sectors and disfavors the SP, VP, and t-IVC states. For
$\kappa>0$, the deviation from chiral limit gives rise to inter-sublattice
interactions which disfavors QH, VH, and t-IVC states. Combining both
perturbations, the K-IVC becomes the most favored state among the listed
symmetry broken states. This qualitative behavior agrees with the observations
in Fig. 2

###  V.3 Hierarchy of Symmetries and Degeneracy of Ground-state Manifold

In the following, we expose the hierarchy of scales and symmetries in TBG by
studying the temperature dependence of various correlation functions in
different limits. For chiral and nonchiral limits, we have $\kappa=0$ and
$0.75$ respectively. For flat and nonflat bands, we used $\hat{H}_{0}=0$ and
$\hat{H}_{0}=\alpha(\hat{H}_{\textrm{BM}}-\big{[}\hat{V}\,\big{]}_{P_{0}})$
respectively in Eq. (III.2). For the nonchiral limit, we take $\alpha=1$
corresponding to a bandwidth of $W=8.3\,\textrm{meV}$ at $\kappa=0.75$
(realistic limit). For the chiral limit $\kappa=0$, we suppress the dispersion
with $\alpha=0.25$ so the bandwidth $W^{*}=3.8\,\textrm{meV}$ stays within the
ordered phase and allows a meaningful comparison. $\alpha=1$ leads to a
correlated semimetal, that we study using $\kappa=0.3$ with
$W=14.3\,\textrm{meV}$. In each limit, we have different global symmetries as
summarized in Fig. 2. Therefore, correlation functions must exhibit an exact
degeneracy induced by the global symmetry. For a detailed symmetry analysis,
we refer to Appendix D.2. The results for the correlation function
$S_{a}({\bm{q}}=0)$ are depicted in Fig. 2 as function of temperature for four
different classes: (a) flat (i.e., $\hat{H}_{0}=0$) chiral, (b) nonflat
chiral, (c) flat nonchiral, and (d,e) nonflat nonchiral. Note that (d) and (e)
have the same global symmetry, but (e) has a larger bandwith. For readability,
we omit the correlation functions for SP and t-IVC order parameters. The
magnitude of SP is equivalent to that of VP by symmetry. The magnitude of
t-IVC is equivalent to the others only at chiral flat limit, and is suppressed
in all other cases. Fig. 2(d),(e) demonstrate that all short-ranged
correlation functions converged with system size while
$S_{\mathrm{K-IVC}}(q=0)$ for $\kappa=0.75$ [Fig. 2(a)-(d)] systematically
increases with system size at low temperature. In fact,
$S_{\mathrm{K-IVC}}(q=0)/L^{2}$ determines the K-IVC order parameter at zero
temperature by first taking the limit $T\rightarrow 0$ and then extrapolating
to the thermodynamic limit.

The correlation function saturates as we decrease the temperature to a value
that scales with the system’s volume for the $L=3,\,6$, as shown in Fig. 2.
Based on this observation, we can conclude that at these temperatures and
system sizes, the system appears to be long-range ordered. However, it should
be emphasized that no true long-range or even algebraic order is expected at
finite temperature in thermodynamic limit, as the order parameter break a
continuous non-Abelian symmetry.

The data confirms the theoretical understanding of the hierarchy of symmetries
in TBG. In the flat chiral limit, conduction and valence bands are degenerate
and the form factors for the Coulomb interaction obey an enhanced symmetry
(see Appendix D.2). As a result, the system exhibits a
$\mathrm{U}(4)_{+}\times\mathrm{U}(4)_{-}$ symmetry, where each
$\mathrm{U}(4)_{\pm}$ acts on the respective Chern sector $C=n_{x}$ in the
band basis ($\tau_{z}\sigma_{z}$ in a sublattice basis of [27]). Within each
Chern sector, there are four electron flavors (valley index $\tau=\pm 1$ and
spin index $s=\pm 1$), and the system is invariant under rotations of those
flavor. Due to this large symmetry group, there is a massive ground state
degeneracy at charge neutrality. All the ground states are annihilated by the
Coulomb interaction. The ground states reside in the union of following
disjoint manifolds in the thermodynamic limit (See Appendix E):

$\displaystyle\quantity(\frac{\mathrm{U}(4)}{\mathrm{U}(2)\times\mathrm{U}(2)})^{\otimes
2}\bigcup\quantity(\frac{\mathrm{U}(4)}{\mathrm{U}(3)\times\mathrm{U}(1)})^{\otimes
2}\bigcup\,\mathbb{Z}_{2}$ (19)

where the three manifolds correspond to the states with Chern number 0, $\pm
2$, and $\pm 4$.

Such a disjoint set of ground state manifolds indicates that the system is
fine-tuned to a first order phase transition. This situation is typically
challenging for the DQMC algorithm due to the breakdown of ergodicity. Note
that the PHS symmetry which guarantees the absence of the sign problem is
preserved by each configuration and therefore enforces an equal sampling of
opposite Chern sectors, $\pm C$. This already facilitates the sampling.
Additionally, we have carefully analyzed QMC runs and fine-tuned updates, e.g.
by using both discrete and continuous auxiliary fields for the ${\bm{q}}=0$
contribution of the Coulomb repulsion. The results presented in Fig. 1 are
taken at finite $\kappa$ where this potential issue is absent.

Let us note a few consequences arising from the different dimensions of the
ground state manifolds. In a finite-size system with $N_{M}=L^{2}$ moiré unit
cells, one can calculate the degeneracy of each manifold using Young’s
tableau. The calculation in Appendix E shows that while the degeneracy of
$C=0$ ground states scales as $N_{M}^{8}$, the degeneracy of $|C|=2$ ground
states scales as $N_{M}^{6}$. For $|C|=4$, there are only two ground states.
Hence, the $C=0$ sector (containing e.g. K-IVC, t-IVC, SP, VP, VH states)
dominates the partition function. Indeed, if we look at the QH channel
(corresponding to $|C|=2,\,4$) in Fig. 2(a), the QH correlation function
decreases with the system size while the other channels (K-IVC, VP, VH, etc.)
increase. This aligns with theoretical expectation: the relative scaling of
the degeneracy between $C=0$ and $\absolutevalue{C}=2$ scales as $1/N_{M}^{2}$
(the $\absolutevalue{C}=4$ only contributes 2 states and can be neglected).
Therefore the QH correlation function at ${\bm{q}}=0$ should scale as
$N_{M}\cdot 1/N_{M}^{2}=1/N_{M}$, where the first factor of $N_{M}$ is the
volume scaling. This accounts for the suppression of the QH channel at low
temperature. Nevertheless, the correlation function is strongly peaked at
${\bm{q}}=0$ giving rise to a very large correlation length (not shown).

As we turn on the dispersion, the perturbation induces anti-ferromagnetic like
interaction among local order parameters [27] that disfavors the t-IVC, SP,
and VP ground states, as illustrated in Tab. 2. As a result, the groundstate
manifold shrinks and the saturation value of the order parameter magnitudes
should increase, which agrees with our observation in Fig. 2(b). One important
caveat is that an anti-ferromagnetic interaction favors totally singlet
ground-states. For example, on a finite-size lattice, the Heisenberg
antiferromagnet has a unique groundstate, unlike the ferromagnet whose ground
state degeneracy scales as the system’s volume. However, once we consider the
Anderson tower of states within a small energy window $[0,1/\beta]$ set by
finite temperature, the number of states inside the window rouhgly scales with
$N_{M}$ as if the order parameter belonged in the following manifold:

$\displaystyle\quantity(\frac{\mathrm{U}(4)}{\mathrm{U}(2)\times\mathrm{U}(2)})\bigcup\quantity(\frac{\mathrm{U}(4)}{\mathrm{U}(3)\times\mathrm{U}(1)})\bigcup\,\mathbb{Z}_{2}$
(20)

Accordingly, the quasi-degeneracy of each manifold with a different Chern
number roughly scales as $N_{M}^{4}$, $N_{M}^{3}$, and $2$, respectively. For
a thorough discussion, see Appendix E. This implies that QH correlation
functions should now scale as $N_{M}\cdot 1/N_{M}\sim\textrm{Const}$, which
aligns with our observation in Fig. 2(b) where $S_{\textrm{QH}}({\bm{q}}=0)$
at different system sizes collapse at low temperature. Note that this scaling
behavior becomes asymptotically exact at thermodynamic limit for any finite
temperature $T>0$ since the energy spacing between Anderson tower states
scales as the inverse of system’s volume. At zero temperature or temperature
scale far below the excitation gap of Anderson tower states, $C=0$, $C=\pm 2$,
and $C=\pm 4$ sectors will collapse into unique ground states respectively.
Therefore, the QH correlation function should exhibit a volume-law scaling
behavior. However, a rough estimate of the Anderson tower excitation energy
for $L=6$ shows that this limit is reached when $T\sim 0.2\,\rm{K}$, which is
lower than the temperature in our DQMC simulations.

We can use the results of Fig. 1 to determine the stiffness of the non-linear
sigma model describing the ground state manifold [38, 57]. This value controls
the dispersion of the low-energy charge neutral bosonic modes [57] as well as
skyrmion excitations [38]. As shown in Appendix G, the stiffness $\rho$ can be
extracted using the correlation length data. By fitting the relation $\xi\sim
e^{\pi\rho/2T}$ which is valid above the saturation temperature Fig. 11, we
extract $\rho\sim 1\,\textrm{meV}$. The value of $\rho$ obtained here is in
rough agreement with the value obtained through Hartree-Fock calculation in
Ref. [38].

Figure 3:  Spectral function $A(\omega,k)$, summed over band, valley and spin,
at the Dirac point ${\bm{k}}=\textbf{K}_{M}$ (a,c) and the BZ center
${\bm{k}}=\Gamma_{M}$ (b,d) for $\kappa=0.75$ (top row) and $\kappa=0.3$
(bottow row) and $L=9\times 9$. The dashed vertical lines indicate the non-
interacting excitation energy, shifted by the charging energy, the dotted line
represent interaction-only limit. We find that the spectra for $\kappa=0.75$
are in fairly good agreement with the excitation energies of the interaction-
only limit. For $\kappa=0.3$, instead, sharp peaks arise at low temperatures
at the Dirac point ${\bm{k}}=\textbf{K}_{M}$, their energy agrees with the
finite-size charging energy, indicating the semi-metallic phase. Note the
spectral weight at $\omega\sim\pm 15\,\textrm{meV}$, reminiscent of the purely
interacting mode. One cannot distinquish those two modes at ${\bm{k}}=\Gamma$.

###  V.4 Spectral Function and Charge Gap

The momentum dependence of the spectral function has already been shown in
Fig. 1 and has been discussed in Sec. II. Let us therefore focus on the
frequency and temperature dependence of spectral function at particular
momenta, presented in Fig. 3. The spectral function is reconstructed using
maximum entropy method [83] from the electron Green’s function $G(\tau)$,
shown in Fig. 10 of Appendix F.

Spectral evolution with temperature in the K-IVC insulator: In Fig. 3(a), we
show the spectral function $A(\omega,\bm{k}=\textbf{K}_{M})$ (at the Dirac
point) for $\kappa=0.75$ at different temperatures. At the lowest
temperatures, we observe two clear peaks located around $\pm
19\,\textrm{meV}$, whose values are very close to the exact excitation
energies of the interaction only Hamiltonian ($H_{0}=0$), which also has a
K-IVC ground state [27]. Clearly, since the $H_{0}$ has gapless excitations at
the Dirac points, the observed gap is a result of interactions.

As we increase the temperature, the K-IVC correlation length become smaller
than one lattice spacing around $10\,\textrm{K}$, as shown in Fig. 1(c) and
Fig. 2(d). Surprisingly, this does not immediately lead to the disappearance
of the gap-like features in the spectral function. Rather, the location of the
peak in $A(\omega,\bm{k}=\textbf{K}_{M})$ broadens and moves towards $\pm
15\,\textrm{meV}$. Such peaks persist for an intermediate range
$10\,\textrm{K}\lesssim T\lesssim 100\,\textrm{K}$. The spectral weight at
$\omega=0$ increases with increasing $T$, which may arise as the superposed
tails of two thermally broadened peaks. By $T=116\,\textrm{K}$, the gap
finally fills in, and there is a single, broad peak centered at $\omega=0$.

The spectrum at the $\Gamma_{M}$-point is shown in Fig. 3(b). The picture is
qualitatively similar to that at the $K_{M}$ point, although the magnitude of
the gap is smaller. The gap features fills in by $T=60\,\textrm{K}$.

To summarize the results for $\kappa=0.75$, the spectra show a pseudo-gapped
regime at intermediate temperatures before a fully gapped state with
substantially long-ranged K-IVC correlations develops below $10\,\rm{K}$. The
energy density obtained by the DQMC simulation for $L=6$ at
$T=1.8\,\textrm{K}$ is $E=-0.377\pm 0.006\,\textrm{meV}$ per moiré unit cell,
which is in good agreement with the self-consistent Hartree-Fock calculation
result, $E=-0.384\,\textrm{meV}$ per moiré unit cell. Furthermore, the
excitation spectrum obtained in the DQMC coincide with the HF solution. This
indicates that the Hartree-Fock calculation is accurate in this case, i.e.,
the ground state is well approximated by a Slater determinant.

Figure 4: Characterization of the strongly correlated semi-metal phase at
$\kappa=0.3$. (a) Proxy for the quasiparticle weight $Z$ of gapless Dirac
particle in comparison to the non-interacting Dirac metal. (b) Flavor resolved
double occupancy averaged over all spin, valley, and momenta indices.

Spectral evolution with temperature in the semi-metal: The second row of Fig.
3(c,d) shows the same quantities as (a,b), but for $\kappa=0.3$. In Fig. 1(j),
we saw that the system behaves like a semimetal at the lowest temperature. The
dispersion of the low-energy quasiparticle peaks is close to that of the bare
Hamiltonian $\hat{H}_{0}$. This is in stark contrast with the K-IVC state at
$\kappa=0.75$, where the dispersion is close to that of the interaction-only
Hamiltonian [Fig. 1(d)].

To understand this behavior better, we focus on ${\bm{k}}=\textbf{K}_{M}$
[Fig. 3(c)]. At $T=116\,\rm{K}$, we observe a broad maximum around $\omega=0$.
At lower temperatures, two distinct narrow peaks appear at $\omega\sim\pm
0.7\,\textrm{meV}$. These peaks can be identified with the quasi-particle
excitation of the Dirac semi-metal. This finite gap originats from the finite-
size charging energy $V_{{\bm{q}}=0}$, which is $0.7\,\textrm{meV}$ for $L=9$,
and vanishes in the thermodynamic limit.

Additionally, there are two broad peaks that give rise to the ‘shoulders’ at
$\omega\sim 15\,\textrm{meV}$, close to the excitation energy of interaction-
only Hamiltonian, shown by the thin dashed lines in Fig. 3(c) [see also Fig.
1(j)]. Thus, the low-temperature spectrum at $\textbf{K}_{M}$ contains two
distinct features, one near the excitation energy of $\hat{H}_{0}$ and another
at the interaction energy scale. This behavior is reminiscent of that of a
Hubbard model, where coherent low-energy quasiparticle peaks can coexist with
high-energy ‘Hubbard bands’ at the scale of the local interaction $U$. The
similarity suggests that electron correlations play an important role in the
semi-metal phase at $\kappa=0.3$. Fig. 3(d) shows the spectral function at
${\bm{k}}=\Gamma$, which is qualitatively similar to that at $\kappa=0.75$.

Next, we present two additional analyses to characterize further the
correlated semi-metal state: (a) the spectral weight at low energies for the
Dirac point and (b) averaged double occupancy per spin, valley, and momentum.
The former,
$\rho_{\mathrm{Dirac}}\equiv(2\pi)^{-1}\int_{-\Delta}^{+\Delta}d\omega
A(\omega,K_{M})$, is a proxy for the quasiparticle weight $Z$ of the Dirac
fermion at the Fermi level, where $\Delta$ is set to twice the charging
energy. This quantity converges to the quasiparticle weight in thermodynamic
limit. In Fig. 4(a), we find that $\rho_{\mathrm{Dirac}}$ increases with
decreasing temperature due to the suppression of thermal spectral broadening
and approaches $\rho_{\mathrm{Dirac}}\sim 2.75$ for $L=3$ and
$\rho_{\mathrm{Dirac}}\sim 2$ for $L=6,9$ at low temperature. Note that
$\rho_{\mathrm{Dirac}}=4$ for the Slater state of a non-interacting symmetric
Dirac semi-metal. Hence, roughly half of the spectral weight has been
transferred from the Dirac quasi-particles to the continuum. This phenomena
resembles the reduction of quasiparticle weight in Hubbard model with large
interaction parameter in a metallic phase before it undergoes a transition to
Mott insulator [93].

Fig. 4(b) shows the averaged flavor-resolved double occupancy,
$D\equiv\sum_{i}\frac{1}{4L^{2}}\expectationvalue{(n_{i,\rm{val}}+n_{i,\rm{con}}-1)^{2}}$,
where $i=(s_{z},\tau_{z},{\bm{k}})$ is a collective index for spin, valley,
and momentum, and $n_{i,\rm{con}}$ ($n_{i,\rm{val}}$) are the occupation
numbers of conduction (valence) band at state $i$, respectively. In the
correlated semi-metal state, we observe that $D$ converges to the value around
$0.2$ at low temperature with increasing system size. Note that in the ground
state of a non-interacting semi-metal, $D\rightarrow 0$ as
$L\rightarrow\infty$. The double occupancy $D$ is a direct measure of
fluctuations in $n_{i}=n_{i,\textrm{val}}+n_{i,\textrm{con}}$. Finite values
of $D$ can have multiple origins, e.g., fluctuations between different
valleys, spins, or momenta, even in the absence of long-range order. However,
fluctuations between conduction and valence band, relevant for the QH and VH
channel, do not contribute to $D$. Fig. 4(b) indicates an extensive double
occupancy even in the limit $T\rightarrow 0$ and thus cannot be explained by
the usual excitation of the semi-metal in a small region around the Dirac
points.

Fig. 2(e) indicates the absence of any long-range correlations at low
temperature in the semi-metal. There is an enhancement of short-ranged
correlations for K-IVC, VH and QH channels at low temperature, while VP and SP
channels are highly suppressed. This indicates that the correlated semi-metal
state consists of significant quantum fluctuations between conduction and
valence bands (VH/QH channels) as well as a quantum fluctuations between
different valleys (K-IVC channel).

To summarize the results for $\kappa=0.3$, we find spectra which show features
of the semi-metal, but also a considerable transfer of spectral weight to a
high-energy correlated mode. The observed semi-metal phase at $\kappa=0.3$ is
in stark contrast with results of self-consistent HF, which predicts a K-IVC
groundstate with a gapped single-particle spectrum with a gap of
$16.7\,\textrm{meV}$ at the $\textbf{K}_{M}$ point. The HF calculation with
$L=6$ predicts an energy density $E=-3.861\,\textrm{meV}$ per moiré unit cell,
while the DQMC result at $T=1.8\,\textrm{K}$ gives $E=-5.388\pm
0.002\,\textrm{meV}$ per moiré unit cell. The energy of the semimetal state is
significantly lower than the energy of the HF K-IVC state, which explains why
the semi-metal state is preferred over the K-IVC state at $\kappa=0.3$, and
indicates that the ground state wavefunction of the correlated semi-metal is
far from a Slater determinant state.

##  VI Conclusions and Outlook

In this work, we have demonstrated that the rich physics of magic-angle
twisted bilayer graphene at charge neutrality can be simulated using the
unbiased, numerically-exact DQMC method, free of the fermion sign problem.
This is a rare example where a nearly-realistic model of a strongly-correlated
electron system can be fully solved. Our simulations reveal that the nature of
the ground state depends sensitively on the bare dispersion bandwidth. Varying
the bandwidth, we found either a semi-metal or a K-IVC ordered state. The
K-IVC state exhibits a non-trivial evolution with temperature, including a
‘pseudogap’ feature in single-particle density of states that onsets far above
the temperature at which the correlation length of the K-IVC order begins to
grow rapidly. The semi-metal phase appears in a regime where Hartree-Fock
would predict an insulating state with a large gap, and exhibits signatures of
strong correlations, such as coexistence of low-energy coherent quasiparticles
and high-energy correlation-induced modes.

Our work paves the way for more detailed studies of TBG and direct comparisons
to experiments. The phase diagram as a function of the single-particle
bandwidth (that can be tuned, e.g., by the distance to the gates) and
temperature can be mapped, including any possible intermediate phases between
the semi-metal and the K-IVC. The finite-temperature crossover regime above
the K-IVC can be characterized in more detail. This regime likely hosts strong
bosonic fluctuations, whose entropy (directly calculable in DQMC) can give
rise to a Pomeranchuk-like effect [48, 49]. Our method to eliminate the sign
problem at small twist angles can be extended to include various physical
perturbations, such as strain, a sub-lattice potential, and perpendicular
electric and magnetic fields. Furthermore, novel electrical transport at
elevated temperatures has been experimentally reported [50, 51], and our setup
will help clarify the intrinsic conductivity of TBG, free of phonons and other
non-electronic scattering mechanisms. These investigations will be reported in
a forthcoming series of works.

Our method should also be useful to probe strong coupling superconductivity
that leaves an imprint on excitations in the parent insulator. Signatures of
incipient superconductivity in future calculations could provide strong
support for a purely electronic pairing mechanism. One such proposal is
Skyrmion superconductivity [38, 32, 94]. Probing for such charge–$2e$
excitations above charge neutrality may shed light on this unusual mechanism
for superconductivity.

###### Acknowledgements.

We thank Michael Zaletel, Zhenjiu Wang, and Fakher F. Assaad for explanation
of their inspiring work [95]. AV was supported by a Simons Investigator award
and by the Simons Collaboration on Ultra-Quantum Matter, which is a grant from
the Simons Foundation (651440, AV). EK was supported by the German National
Academy of Sciences Leopoldina through grant LPDS 2018-02 Leopoldina
fellowship. EB and JH were supported by the European Research Council (ERC)
under grant HQMAT (Grant Agreement No. 817799), the Israel-US Binational
Science Foundation (BSF), CRC 183 of the Deutsche Forschungsgemeinschaft, and
by a Research grant from Irving and Cherna Moskowitz. JYL was supported by
Gordon and Betty Moore Foundation through Grant GBMF8690 to UCSB and by the
National Science Foundation under Grant No. NSF PHY-1748958. The auxillary
field QMC simulations were carried out using the ALF package available at
https://alf.physik.uni-wuerzburg.de. Note Added: After the completion of this
work, Ref. [96] appeared which notes absence of a sign problem in certain
limits of TBG.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Bistritzer and MacDonald [2011a] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011a).
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Phys. Rev. B 86, 155449 (2012).
  * Xu and Balents [2018] C. Xu and L. Balents, Topological superconductivity in twisted multilayer graphene, Phys. Rev. Lett. 121, 087001 (2018).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of mott insulating behavior and superconductivity in twisted bilayer graphene, Phys. Rev. X 8, 031089 (2018).
  * Isobe  _et al._ [2018] H. Isobe, N. F. Q. Yuan, and L. Fu, Unconventional superconductivity and density waves in twisted bilayer graphene, Phys. Rev. X 8, 041041 (2018).
  * Thomson  _et al._ [2018] A. Thomson, S. Chatterjee, S. Sachdev, and M. S. Scheurer, Triangular antiferromagnetism on the honeycomb lattice of twisted bilayer graphene, Phys. Rev. B 98, 075109 (2018).
  * You and Vishwanath [2019] Y.-Z. You and A. Vishwanath, Superconductivity from valley fluctuations and approximate so(4) symmetry in a weak coupling theory of twisted bilayer graphene, npj Quantum Materials 4, 16 (2019).
  * Kang and Vafek [2018] J. Kang and O. Vafek, Symmetry, maximally localized wannier states, and a low-energy model for twisted bilayer graphene narrow bands, Phys. Rev. X 8, 031088 (2018).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. 124, 097601 (2020).
  * Lin and Nandkishore [2019] Y.-P. Lin and R. M. Nandkishore, Chiral twist on the high-${T}_{c}$ phase diagram in moiré heterostructures, Phys. Rev. B 100, 085136 (2019).
  * Dodaro  _et al._ [2018] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun, and C. Wang, Phases of a phenomenological model of twisted bilayer graphene, Phys. Rev. B 98, 075154 (2018).
  * Padhi  _et al._ [2018] B. Padhi, C. Setty, and P. W. Phillips, Doped twisted bilayer graphene near magic angles: Proximity to wigner crystallization, not mott insulation, Nano letters 18, 6175 (2018).
  * Wu  _et al._ [2018] F. Wu, A. H. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene, Phys. Rev. Lett. 121, 257001 (2018).
  * Lian  _et al._ [2018] B. Lian, Z. Wang, and B. A. Bernevig, Twisted bilayer graphene: A phonon driven superconductor, arXiv preprint arXiv:1807.04382 (2018).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions, Phys. Rev. B 98, 085435 (2018).
  * Roy and Juričić [2019] B. Roy and V. Juričić, Unconventional superconductivity in nearly flat bands in twisted bilayer graphene, Phys. Rev. B 99, 121407 (2019).
  * Zhang  _et al._ [2019] Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Nearly flat chern bands in moiré superlattices, Phys. Rev. B 99, 075127 (2019).
  * Khalaf  _et al._ [2019] E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Magic angle hierarchy in twisted graphene multilayers, Phys. Rev. B 100, 085109 (2019).
  * Mora  _et al._ [2019] C. Mora, N. Regnault, and B. A. Bernevig, Flatbands and perfect metal in trilayer moiré graphene, Phys. Rev. Lett. 123, 026402 (2019).
  * Cea  _et al._ [2019a] T. Cea, N. R. Walet, and F. Guinea, Twists and the electronic structure of graphitic materials,  _Nano Letters_ , Nano Letters 19, 8683 (2019a).
  * Fang  _et al._ [2019] S. Fang, S. Carr, Z. Zhu, D. Massatt, and E. Kaxiras, Angle-dependent ab initio low-energy hamiltonians for a relaxed twisted bilayer graphene heterostructure (2019), arXiv:1908.00058 [cond-mat.mes-hall] .
  * Bi  _et al._ [2019] Z. Bi, N. F. Q. Yuan, and L. Fu, Designing flat bands by strain, Phys. Rev. B 100, 035448 (2019).
  * Lee  _et al._ [2019] J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, and A. Vishwanath, Theory of correlated insulating behaviour and spin-triplet superconductivity in twisted double bilayer graphene, Nature Communications 10, 5333 (2019).
  * Da Liao  _et al._ [2019] Y. Da Liao, Z. Y. Meng, and X. Y. Xu, Valence bond orders at charge neutrality in a possible two-orbital extended hubbard model for twisted bilayer graphene, Phys. Rev. Lett. 123, 157601 (2019).
  * Bultinck  _et al._ [2020a] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Phys. Rev. X 10, 031034 (2020a).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Kekulé spiral order at all nonzero integer fillings in twisted bilayer graphene (2021), arXiv:2105.05857 [cond-mat.str-el] .
  * Kang and Vafek [2020] J. Kang and O. Vafek, Non-abelian dirac node braiding and near-degeneracy of correlated phases at odd integer filling in magic-angle twisted bilayer graphene, Phys. Rev. B 102, 035161 (2020).
  * Soejima  _et al._ [2020] T. Soejima, D. E. Parker, N. Bultinck, J. Hauschild, and M. P. Zaletel, Efficient simulation of moiré materials using the density matrix renormalization group, Phys. Rev. B 102, 205111 (2020).
  * Huang  _et al._ [2020] Y. Huang, P. Hosur, and H. K. Pal, Quasi-flat-band physics in a two-leg ladder model and its relation to magic-angle twisted bilayer graphene, Phys. Rev. B 102, 155429 (2020).
  * Chatterjee  _et al._ [2020] S. Chatterjee, M. Ippoliti, and M. P. Zaletel, Skyrmion superconductivity: Dmrg evidence for a topological route to superconductivity (2020), arXiv:2010.01144 [cond-mat.str-el] .
  * Hejazi  _et al._ [2021] K. Hejazi, X. Chen, and L. Balents, Hybrid wannier chern bands in magic angle twisted bilayer graphene and the quantized anomalous hall effect, Phys. Rev. Research 3, 013242 (2021).
  * Liu  _et al._ [2019] J. Liu, J. Liu, and X. Dai, Pseudo landau level representation of twisted bilayer graphene: Band topology and implications on the correlated insulating phase, Phys. Rev. B 99, 155415 (2019).
  * Vafek and Kang [2020] O. Vafek and J. Kang, Renormalization group study of hidden symmetry in twisted bilayer graphene with coulomb interactions, Phys. Rev. Lett. 125, 257602 (2020).
  * Xie  _et al._ [2021a] F. Xie, A. Cowsik, Z.-D. Song, B. Lian, B. A. Bernevig, and N. Regnault, Twisted bilayer graphene. vi. an exact diagonalization study at nonzero integer filling, Phys. Rev. B 103, 205416 (2021a).
  * Da Liao  _et al._ [2021] Y. Da Liao, J. Kang, C. N. Breiø, X. Y. Xu, H.-Q. Wu, B. M. Andersen, R. M. Fernandes, and Z. Y. Meng, Correlation-induced insulating topological phases at charge neutrality in twisted bilayer graphene, Phys. Rev. X 11, 011014 (2021).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Charged skyrmions and topological origin of superconductivity in magic-angle graphene, Science Advances 7, 10.1126/sciadv.abf5299 (2021).
  * Utama  _et al._ [2021] M. I. B. Utama, R. J. Koch, K. Lee, N. Leconte, H. Li, S. Zhao, L. Jiang, J. Zhu, K. Watanabe, T. Taniguchi, P. D. Ashby, A. Weber-Bargioni, A. Zettl, C. Jozwiak, J. Jung, E. Rotenberg, A. Bostwick, and F. Wang, Visualization of the flat electronic band in twisted bilayer graphene near the magic angle twist, Nature Physics 17, 184 (2021).
  * Lisi  _et al._ [2021] S. Lisi, X. Lu, T. Benschop, T. A. de Jong, P. Stepanov, J. R. Duran, F. Margot, I. Cucchi, E. Cappelli, A. Hunter, A. Tamai, V. Kandyba, A. Giampietri, A. Barinov, J. Jobst, V. Stalman, M. Leeuwenhoek, K. Watanabe, T. Taniguchi, L. Rademaker, S. J. van der Molen, M. P. Allan, D. K. Efetov, and F. Baumberger, Observation of flat bands in twisted bilayer graphene, Nature Physics 17, 189 (2021).
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Failure of nielsen-ninomiya theorem and fragile topology in two-dimensional systems with space-time inversion symmetry: Application to twisted bilayer graphene at magic angle, Phys. Rev. X 9, 021013 (2019).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, All magic angles in twisted bilayer graphene are topological, Phys. Rev. Lett. 123, 036401 (2019).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Faithful tight-binding models and fragile topology of magic-angle bilayer graphene, Phys. Rev. B 99, 195455 (2019).
  * Hejazi  _et al._ [2019] K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Multiple topological transitions in twisted bilayer graphene near the first magic angle, Physical Review B 99, 035111 (2019).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Twisted bilayer graphene. ii. stable symmetry anomaly, Physical Review B 103, 10.1103/physrevb.103.205412 (2021).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, and et al., Cascade of phase transitions and dirac revivals in magic-angle graphene, Nature 582, 203–208 (2020).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Cascade of electronic transitions in magic-angle twisted bilayer graphene, Nature 582, 198–202 (2020).
  * Saito  _et al._ [2021] Y. Saito, F. Yang, J. Ge, X. Liu, T. Taniguchi, K. Watanabe, J. I. A. Li, E. Berg, and A. F. Young, Isospin pomeranchuk effect in twisted bilayer graphene, Nature 592, 220 (2021).
  * Rozen  _et al._ [2021] A. Rozen, J. M. Park, U. Zondiner, Y. Cao, D. Rodan-Legrain, T. Taniguchi, K. Watanabe, Y. Oreg, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, Entropic evidence for a pomeranchuk effect in magic-angle graphene, Nature 592, 214 (2021).
  * Cao  _et al._ [2020] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Strange metal in magic-angle graphene with near planckian dissipation, Physical Review Letters 124, 10.1103/physrevlett.124.076801 (2020).
  * Polshyn  _et al._ [2019] H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Large linear-in-temperature resistivity in twisted bilayer graphene, Nature Physics 15, 1011–1016 (2019).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Phys. Rev. X 8, 031087 (2018).
  * Seo  _et al._ [2019] K. Seo, V. N. Kotov, and B. Uchoa, Ferromagnetic mott state in twisted graphene bilayers at the magic angle, Phys. Rev. Lett. 122, 246402 (2019).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong coupling phases of partially filled twisted bilayer graphene narrow bands, Phys. Rev. Lett. 122, 246401 (2019).
  * Chen  _et al._ [2020] B.-B. Chen, Y. Da Liao, Z. Chen, O. Vafek, J. Kang, W. Li, and Z. Y. Meng, Realization of topological mott insulator in a twisted bilayer graphene lattice model, arXiv preprint arXiv:2011.07602 (2020).
  * Bistritzer and MacDonald [2011b] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011b).
  * Khalaf  _et al._ [2020] E. Khalaf, N. Bultinck, A. Vishwanath, and M. P. Zaletel, Soft modes in magic angle twisted bilayer graphene, arXiv preprint arXiv:2009.14827 (2020).
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Twisted bilayer graphene. iv. exact insulator ground states and phase diagram, Phys. Rev. B 103, 205414 (2021).
  * Liu  _et al._ [2021] S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality, Phys. Rev. Research 3 (2021).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers, Proceedings of the National Academy of Sciences 115, 13174 (2018).
  * Cea  _et al._ [2019b] T. Cea, N. R. Walet, and F. Guinea, Electronic band structure and pinning of fermi energy to van hove singularities in twisted bilayer graphene: A self-consistent approach, Phys. Rev. B 100, 205113 (2019b).
  * Cea and Guinea [2020] T. Cea and F. Guinea, Band structure and insulating states driven by coulomb interaction in twisted bilayer graphene, Phys. Rev. B 102, 045107 (2020).
  * Liu and Dai [2021] J. Liu and X. Dai, Theories for the correlated insulating states and quantum anomalous hall effect phenomena in twisted bilayer graphene, Phys. Rev. B 103, 035427 (2021).
  * Zhang  _et al._ [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Correlated insulating phases of twisted bilayer graphene at commensurate filling fractions: A hartree-fock study, Phys. Rev. B 102, 035136 (2020).
  * Parker  _et al._ [2020] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Strain-induced quantum phase transitions in magic angle graphene, arXiv preprint arXiv:2012.09885 (2020).
  * Xie  _et al._ [2021b] F. Xie, A. Cowsik, Z.-D. Song, B. Lian, B. A. Bernevig, and N. Regnault, Twisted bilayer graphene. vi. an exact diagonalization study at nonzero integer filling, Phys. Rev. B 103, 205416 (2021b).
  * Potasz  _et al._ [2021] P. Potasz, M. Xie, and A. H. MacDonald, Exact diagonalization for magic-angle twisted bilayer graphene, arXiv preprint arXiv:2102.02256 (2021).
  * Chichinadze  _et al._ [2020] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Valley magnetism, nematicity, and density wave orders in twisted bilayer graphene, Phys. Rev. B 102, 125120 (2020).
  * Classen  _et al._ [2020] L. Classen, A. V. Chubukov, C. Honerkamp, and M. M. Scherer, Competing orders at higher-order van hove points, Phys. Rev. B 102, 125141 (2020).
  * Sugiyama and Koonin [1986] G. Sugiyama and S. Koonin, Auxiliary field monte-carlo for quantum many-body ground states, Annals of Physics 168, 1 (1986).
  * Sorella  _et al._ [1989] S. Sorella, S. Baroni, R. Car, and M. Parrinello, A novel technique for the simulation of interacting fermion systems, Europhysics Letters (EPL) 8, 663 (1989).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science , eaav1910 (2019).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Untying the insulating and superconducting orders in magic-angle graphene, Nature 583, 375 (2020).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Electronic correlations in twisted bilayer graphene near the magic angle, Nature Physics 10.1038/s41567-019-0606-5 (2019).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene, Nature 573, 91 (2019).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Maximized electron interactions at the magic angle in twisted bilayer graphene, Nature 572, 95 (2019).
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene, Nature 572, 101 (2019).
  * Thomson and Alicea [2021] A. Thomson and J. Alicea, Recovery of massless dirac fermions at charge neutrality in strongly interacting twisted bilayer graphene with disorder, Phys. Rev. B 103, 125138 (2021).
  * Note [1] Including additional bands does not introduce a sign problem, and can in principle be incorporated in the DQMC simulations. However, the computational complexity increases as the number of bands to the third power.
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017).
  * Bultinck  _et al._ [2020b] N. Bultinck, S. Chatterjee, and M. P. Zaletel, Mechanism for anomalous hall ferromagnetism in twisted bilayer graphene, Phys. Rev. Lett. 124, 166601 (2020b).
  * Beach [2004] K. S. D. Beach, Identifying the maximum entropy method as a special limit of stochastic analytic continuation, arXiv e-prints , cond-mat/0403055 (2004), arXiv:cond-mat/0403055 [cond-mat.str-el] .
  * Note [2] This broadening may be an artifact of the maximum entropy algorithm, whose intrinsic frequency resolution is proportional to temperature, or may represent the physical lifetime of the excitations at finite $T$.
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of magic angles in twisted bilayer graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte carlo calculations of coupled boson-fermion systems. i, Phys. Rev. D 24, 2278 (1981).
  * Ippoliti  _et al._ [2018] M. Ippoliti, R. S. K. Mong, F. F. Assaad, and M. P. Zaletel, Half-filled landau levels: A continuum and sign-free regularization for three-dimensional quantum critical points, Phys. Rev. B 98, 235108 (2018).
  * Note [3] The main reason for the decay is the decay of the form factors contained inside $\bar{\rho}_{\bm{q}}$ from the flat band projection.
  * Note [4] First quantized matrix for such operator should have zero trace, unless it is chemical potential.
  * Bercx  _et al._ [2017] M. Bercx, F. Goth, J. S. Hofmann, and F. Assaad, The alf (algorithms for lattice fermions) project release 1.0. documentation for the auxiliary field quantum monte carlo code, SciPost Physics 3, 013 (2017), arXiv:1704.00131 [cond-mat.str-el] .
  * Collaboration  _et al._ [2021] A. Collaboration, F. F. Assaad, M. Bercx, F. Goth, A. Götz, J. S. Hofmann, E. Huffman, Z. Liu, F. P. Toldin, J. S. E. Portela, and J. Schwab, The alf (algorithms for lattice fermions) project release 2.0. documentation for the auxiliary-field quantum monte carlo code (2021), arXiv:2012.11914 [cond-mat.str-el] .
  * Parisen Toldin  _et al._ [2015] F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Fermionic quantum criticality in honeycomb and $\pi$-flux hubbard models: Finite-size scaling of renormalization-group-invariant observables from quantum monte carlo, Phys. Rev. B 91, 165108 (2015).
  * Georges  _et al._ [2004] A. Georges, S. Florens, and T. A. Costi, The Mott transition: Unconventional transport, spectral weight transfers, and critical behaviour, Journal de Physique IV 114, 165 (2004), arXiv:cond-mat/0311520 [cond-mat.str-el] .
  * Christos  _et al._ [2020] M. Christos, S. Sachdev, and M. S. Scheurer, Superconductivity, correlated insulators, and wess–zumino–witten terms in twisted bilayer graphene, Proceedings of the National Academy of Sciences 117, 29543–29554 (2020).
  * Wang  _et al._ [2021] Z. Wang, M. P. Zaletel, R. S. K. Mong, and F. F. Assaad, Phases of the ($2+1$) dimensional so(5) nonlinear sigma model with topological term, Phys. Rev. Lett. 126, 045701 (2021).
  * Zhang  _et al._ [2021] X. Zhang, G. Pan, Y. Zhang, J. Kang, and Z. Y. Meng, Momentum space quantum Monte Carlo on twisted bilayer Graphene, arXiv e-prints , arXiv:2105.07010 (2021), arXiv:2105.07010 [cond-mat.str-el] .
  * Hirsch [1985] J. E. Hirsch, Two-dimensional hubbard model: Numerical simulation study, Phys. Rev. B 31, 4403 (1985).
  * Wu and Zhang [2005] C. Wu and S.-C. Zhang, Sufficient condition for absence of the sign problem in the fermionic quantum monte carlo algorithm, Phys. Rev. B 71, 155115 (2005).
  * Li  _et al._ [2015a] Z.-X. Li, Y.-F. Jiang, and H. Yao, Solving the fermion sign problem in quantum Monte Carlo simulations by Majorana representation, Phys. Rev. B 91, 241117 (2015a), arXiv:1408.2269 [cond-mat.str-el] .
  * Li  _et al._ [2015b] Z.-X. Li, Y.-F. Jiang, and H. Yao, Fermion-sign-free Majarana-quantum-Monte-Carlo studies of quantum critical phenomena of Dirac fermions in two dimensions, New Journal of Physics 17, 085003 (2015b), arXiv:1411.7383 [cond-mat.str-el] .
  * Li  _et al._ [2016] Z.-X. Li, Y.-F. Jiang, and H. Yao, Majorana-time-reversal symmetries: A fundamental principle for sign-problem-free quantum monte carlo simulations, Phys. Rev. Lett. 117, 267002 (2016).
  * Repellin  _et al._ [2020] C. Repellin, Z. Dong, Y.-H. Zhang, and T. Senthil, Ferromagnetism in narrow bands of moiré superlattices, Phys. Rev. Lett. 124, 187601 (2020).
  * Note [5] Note that the groundstates with the same Chern number can be rotated into each other by the global symmetry.
  * Greiter and Rachel [2007] M. Greiter and S. Rachel, Valence bond solids for $\mathrm{SU}(n)$ spin chains: Exact models, spinon confinement, and the haldane gap, Phys. Rev. B 75, 184441 (2007).
  * Note [6] http://cftp.ist.utl.pt/~gernot.eichmann/2014-hadron-physics/hadron-app-1.pdf.
  * Motome and Imada [1997] Y. Motome and M. Imada, A quantum monte carlo method and its applications to multi-orbital hubbard models, Journal of the Physical Society of Japan 66, 1872 (1997), http://dx.doi.org/10.1143/JPSJ.66.1872 .
  * Assaad  _et al._ [1997] F. F. Assaad, M. Imada, and D. J. Scalapino, Charge and spin structures of a ${d}_{{x}^{2}-{y}^{2}}$ superconductor in the proximity of an antiferromagnetic mott insulator, Phys. Rev. B 56, 15001 (1997).
  * Note [7] We employed the block-diagonal structure in both valley and spin. Otherwise, the linear dimension would increase by a factor of up to 4\.
  * Note [8] For the Hubbard interaction, the matrix actually is only a scalar number or a $2\times 2$ matrix at the most.
  * Polyakov [1975] A. M. Polyakov, Interaction of goldstone particles in two dimensions. applications to ferromagnets and massive yang-mills fields, Physics Letters B 59, 79 (1975).
  * Hikami [1981] S. Hikami, Three-loop ß-functions of non-linear $\sigma$ models on symmetric spaces, Physics Letters B 98, 208 (1981).
  * Wegner [1989] F. Wegner, Four-loop-order $\beta$-function of nonlinear $\sigma$-models in symmetric spaces, Nuclear Physics B 316, 663 (1989).
  * Evers and Mirlin [2008] F. Evers and A. D. Mirlin, Anderson transitions, Reviews of Modern Physics 80, 1355 (2008).
  * Note [9] This is in a sense similar to distinguishing the charge and spin order in a spin-triplet superconductor. While the former is associated with a diverging correlations function at finite $T=T_{\mathrm{BKT}}$, the correlations function of latter only diverges at $T=0$.

##  Appendix A Sign Problem

In the Appendix, we revisit the sign-problem in DQMC. We clarify some
subtleties that can arise from two different perspectives to look into a sign
problems: first quantized and second-quantized Hamiltonians.

In Hirsch’s original paper [97], the sign problem of fermion quantum Monte
Carlo method was discussed in many-body basis. He showed that in the half-
filled Hubbard model on a bipartite lattice, the partition function weight of
a given auxiliary field configuration is positive by applying a many-body
basis transformation that leaves the trace invariant. However, in the same
paper, he also showed how the trace of a many-body Hamiltonian can be
simplified, i.e., one can calculate the trace as the determinant of the
corresponding first-quantized Hamiltonian. This form was later exploited by Wu
et al. [98], where it was shown that the anti-unitary symmetry ${\cal T}$ of a
_single-particle_ (or first-quantized) quadratic Hamiltonian obtained by
Hubbard-Stratonovich decoupling can guarantee the absence of the sign problem
if ${\cal T}^{2}=-1$. This argument was further extended by Li et al. [99,
100, 101], where they discuss how the set of anti-unitary symmetries of a
first quantized (HS-decoupled) Hamiltonain in the Majorana representation can
guarantee the absence of sign problem. The Majorana representation is nice,
since the unitary transformation of a first-quantized Hamiltonian in the
electron representation cannot accommodate particle-hole symmetry, which is a
valid basis transformation in a many-body level. The problem is that single-
particle Hamiltonian encodes the information about the many-body physics with
an explicit $U(1)$ number conservation symmetry. Therefore, if we want to
discuss a generic particle number changing transformation at a first-quantized
level, it advantageous to use the Majorana representation.

However, without bothering to reformulate the problem in Majorana basis, we
can simply work with a (HS-decoupled) many-body Hamiltonian as Hirsch
originally discussed. There is a subtlety. In his original discussion, he only
considered a bipartite system with a real Hamiltonian and associated many-body
particle-hole transformation ${\cal{C}}$, and it was not important whether
${\cal{C}}$ is unitary or anti-unitary at many-body level. For a more general
approach, we have to understand the role of anti-unitary many-body symmetries,
which will be provided below. First note that
$(\bra{n}M\ket{m})^{*}=\bra{m}M^{\dagger}\ket{n}$. The following identity is
useful in the discussion of anti-unitary transformation:

$A=\Tr M\equiv\sum_{n}\bra{n}M\ket{n},\quad\Rightarrow\quad
A^{*}=\sum_{n}(\bra{n}M\ket{n})^{*}=\sum_{n}\bra{n}M^{\dagger}\ket{n}=\Tr
M^{\dagger}=\Tr M^{*},$ (21)

where the last equality comes from the invariance of the trace under
transpose. This part is important as the trotterized partition function cares
about the order of matrix products. Let ${\cal K}$ be the complex-conjugation.
Now define a general anti-unitary symmetry ${\cal T}={\cal U}{\cal K}$, where
${\cal U}$ is a many-body unitary basis transformation. Note that if ${\cal
U}$ is a many-body _unitary_ basis transformation, it can be shown that

$\displaystyle\Tr_{c}e^{-\beta\hat{H}}=\Tr_{c}{\cal U}^{-1}e^{-\beta{\cal
U}\hat{H}{\cal U}^{-1}}{\cal U}=\Tr_{{\cal U}c}e^{-\beta{\cal U}\hat{H}{\cal
U}^{-1}}=\Tr_{c}e^{-\beta{\cal U}\hat{H}{\cal U}^{-1}}$ (22)

where the last equality comes from the cyclic property of trace. By combining
Eq. (22) and Eq. (21), we can show that

$A\equiv\Tr_{c}\prod_{l}e^{-\hat{H}_{l}}\quad\Rightarrow\quad
A^{*}=\Tr_{c}\prod_{l}e^{-{\cal T}\hat{H}_{l}{\cal T}^{-1}}.$ (23)

To illustrate the point, consider a bipartite Hubbard model with nearest-
neighbor hopping terms and on-site repulsive interaction at charge-neutrality:

$\hat{H}=\sum_{\expectationvalue{ij}}\sum_{s}c^{\dagger}_{i,s}H_{ij}c_{j,s}^{\,}+\frac{U}{2}\sum_{i}(n_{i,\uparrow}+n_{i,\downarrow}-1)^{2}$
(24)

The Hirsch’s original approach performed a HS-decoupling into a spin-channel
when the interaction is repulsive ($U>0$) so that the resulting decoupled
interaction term is real. However, even when $U>0$, we can also decouple the
interaction into a density channel with a price that the interaction term
would have an imaginary coefficient, thus non-Hermitian. Under Hubbard-
Stratonovich transformation, the trotterized on-site repulsive interaction
term would be given as the following:

$\displaystyle e^{-\frac{U\delta}{2}(n_{\uparrow}+n_{\downarrow}-1)^{2}}$
$\displaystyle=\frac{1}{2}\sum_{\eta=\pm
1}e^{i\alpha\eta(n_{\uparrow}+n_{\downarrow}-1)}=\frac{1}{2}\sum_{\eta=\pm
1}e^{i\alpha\eta(n_{\uparrow}-1/2)}\cdot e^{i\alpha\eta(n_{\downarrow}-1/2)}$
(25)

where $\alpha=\cos^{-1}e^{-U\delta/2}$ ($U>0$). Under density-channel
decoupling, both spins have the same HS-decoupled Hamiltonians, implying that
$W_{\uparrow}[\eta(\tau)]=W_{\downarrow}[\eta(\tau)]$ (c.f. Eq. (7)). Now,
consider a many-body basis transformation ${\cal U}$:

${\cal U}:c^{\dagger}_{r,s}\mapsto
e^{i\vec{K}\cdot\vec{r}}c^{\,}_{r,s},\qquad\vec{K}=(\pi,\pi)$ (26)

which exchanges particles and holes. Note that ${\cal U}$ alone is not a
symmetry yet. If we consider an anti-unitary operator ${\cal C}\equiv{\cal
U}{\cal K}$, then we see that

$\displaystyle{\cal
C}\Big{[}\sum_{\expectationvalue{ij}}\sum_{s}c^{\dagger}_{i,s}H_{ij}c_{j,s}^{\,}\Big{]}{\cal
C}^{-1}$
$\displaystyle=-\sum_{\expectationvalue{ij}}\sum_{s}c^{\,}_{i,s}H^{*}_{ij}c_{j,s}^{\dagger}=\sum_{\expectationvalue{ji}}\sum_{s}c^{\dagger}_{j,s}H_{ji}c_{i,s}^{\,}$
$\displaystyle{\cal C}\quantity[i\alpha
s_{r,\tau}\quantity(\hat{n}_{r,\eta}-\frac{1}{2})]{\cal C}^{-1}$
$\displaystyle=(-i)\alpha\eta_{r,\tau}\quantity(\frac{1}{2}-\hat{n}_{r,\eta})=i\alpha\eta_{r,\tau}\quantity(\hat{n}_{r,\eta}-\frac{1}{2}).$
(27)

Therefore, Eq. (23) implies the following:

$W^{*}_{\uparrow}[\eta(\tau)]=W_{\uparrow}[\eta(\tau)],\quad
W^{*}_{\downarrow}[\eta(\tau)]=W_{\downarrow}[\eta(\tau)],$ (28)

where $W$ is the weight for a given auxiliary field configuration (c.f. Eq.
(10)) Since $W_{\uparrow}[\eta(\tau)]=W_{\downarrow}[\eta(\tau)]$, we conclude
that

$W=\Tr_{\eta(\tau_{l})}\Tr_{\eta(\tau_{l})}\absolutevalue{W_{\uparrow}[\eta(\tau)]}^{2}>0.$
(29)

This is particularly nice since spin-channel decoupling is not possible for
systems with higher spins, i.e., more than two flavors of electrons per site.
Even when we have two flavors, this is often a better scheme as the density-
decoupling respects the spin-rotation symmetry, while the $S_{z}$-decoupling
does not.

##  Appendix B Twisted Bilayer Graphene and Symmetry

Figure 5: Real space lattice and Brillouin zone for a single graphene layer,
and Brillouin zone origniated from ${\bm{K}}_{+}$ valley for a twisted bilayer
graphene. The figure illustrates all vectors labeled for the discussion in the
manuscript.

Let $\theta$ be a (counterclockwise) twist angle of top layer with respect to
the bottom layer. Then, the TBG hamiltonian for valley $\tau$ has the
following form:

$\displaystyle H_{\tau}$
$\displaystyle=\sum_{{\bm{k}}}\Big{[}\,c_{{\bm{k}},\tau,t}^{\dagger}h_{\theta/2,\tau}({\bm{k}})c_{{\bm{k}},\tau,t}+c_{{\bm{k}},\tau,b}^{\dagger}h_{-\theta/2,\tau}({\bm{k}})c_{{\bm{k}},\tau,b}$
$\displaystyle\quad+\sum_{n}\quantity(c^{\dagger}_{{\bm{k}},\tau,b}T_{n,\tau}c^{~{}}_{{\bm{k}}+{\bm{q}}_{n},\tau,t}+c^{\dagger}_{{\bm{k}}+{\bm{q}}_{n},\tau,t}T^{\dagger}_{n,\tau}c^{~{}}_{{\bm{k}},\tau,b})\Big{]},$
(30)

where $t,b$ denotes top and bottom layers and $h_{\theta,\tau}({\bm{k}})$ is
the matrix for the Dirac Hamiltonian rotated by $\theta$ in counter-clockwise
direction. In other words,
$h_{\theta,\tau}({\bm{k}})=h_{0,\tau}(R_{-\theta}{\bm{k}})$, and the moiré
interlayer hopping term is defined as $T_{n,\tau}=w_{0}+w_{1}e^{2\pi
in\sigma_{z}\tau/3}\sigma_{x}e^{-2\pi in\sigma_{z}\tau/3}$. Here,
$c^{\dagger}_{{\bm{k}},\tau,t/b}$ is a 4-components (spin and sublattice)
electron creation operator with momentum ${\bm{k}}$ measured with respect to
individual Dirac points at valley $\tau$ for top/bottom layers, respectively.
In this paper, we choose the convention where Dirac Hamiltonian is given as
$h_{0,\tau}({\bm{k}})\equiv\tau
k_{1}\sigma_{1}+k_{2}\sigma_{2}={\bm{k}}_{\tau}\cdot{\bm{\sigma}}$ where
${\bm{k}}_{\tau}\equiv(\tau k_{1},k_{2})$. Now, if we rotate momentum by angle
$\varphi$ in counter-clockwise direction, we get the following:

$\displaystyle(R_{\varphi}{\bm{k}})_{\tau}\cdot{\bm{\sigma}}$
$\displaystyle=\tau(k_{1}\cos\varphi-
k_{2}\sin\varphi)\sigma_{1}+(k_{1}\sin\varphi+k_{2}\cos\varphi)\sigma_{2}$
$\displaystyle=\tau
k_{1}(\cos\varphi\sigma_{1}+\sin\varphi\sigma_{2}\tau)+k_{2}(-\sin\varphi\tau\sigma_{1}+\cos\varphi\sigma_{2}).$
(31)

Using the fact that
$e^{i\varphi\sigma_{3}\tau}=\cos\varphi+i\sin\varphi\sigma_{3}\tau$, we can
show that

$\displaystyle
e^{-i\varphi\sigma_{3}\tau/2}\tau\sigma_{1}e^{i\varphi\sigma_{3}\tau/2}$
$\displaystyle=\tau\sigma_{1}e^{i\varphi\sigma_{3}\tau}=\tau(\cos\varphi\sigma_{1}+\sin\varphi\sigma_{2}\tau)$
$\displaystyle
e^{-i\varphi\sigma_{3}\tau/2}\sigma_{2}e^{i\varphi\sigma_{3}\tau/2}$
$\displaystyle=\sigma_{2}e^{i\varphi\sigma_{3}\tau}=\cos\varphi\sigma_{2}-\sin\varphi\tau\sigma_{1}.$
(32)

Therefore,
$(R_{\varphi}{\bm{k}})_{\tau}\cdot{\bm{\sigma}}=({\bm{k}}_{\tau}\cdot{\bm{\sigma}})e^{i\varphi\sigma_{3}\tau}$.
As $\varphi=-\theta\mu_{z}/2$, we find that
$(R_{\varphi}{\bm{k}})_{\tau}\cdot{\bm{\sigma}}=({\bm{k}}_{\tau}\cdot{\bm{\sigma}})e^{-i\theta\sigma_{3}\tau\mu_{3}/2}$.
Therefore, we can unify the descriptions for different layers and valleys in
the following form:

$\displaystyle\hat{H}$
$\displaystyle=\sum_{{\bm{k}}}c^{\dagger}_{\bm{k}}H({\bm{k}})c_{\bm{k}}+\sum_{{\bm{k}},{\bm{k}}^{\prime}}c^{\dagger}_{{\bm{k}}^{\prime}}T({\bm{k}}^{\prime},{\bm{k}})c^{\,}_{\bm{k}}$
$\displaystyle H({\bm{k}})$ $\displaystyle=\hbar
v_{F}\quantity(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y})e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}$
$\displaystyle T({\bm{k}}^{\prime},{\bm{k}})$
$\displaystyle=\frac{1}{2}\sum_{n}\delta_{{\bm{k}}^{\prime},{\bm{k}}-\tau_{z}{\bm{q}}_{n}}\quantity(\mu_{x}-i\mu_{y})T_{n}+\frac{1}{2}\sum_{n}\delta_{{\bm{k}}^{\prime},{\bm{k}}+\tau_{z}{\bm{q}}_{n}}T^{\dagger}_{n}\quantity(\mu_{x}+i\mu_{y})$
(33)

where $\sigma$, $\tau$ and $\mu$ act on sublattice, valley and layer
respectively. Here ${\bm{q}}_{n}={\cal R}_{2\pi n/3}\,{\bm{q}}_{0}$, where
${\cal R}_{\varphi}=e^{i\varphi\sigma_{y}}$ is a rotation matrix,
${\bm{q}}_{0}=\frac{8\pi\sin(\theta/2)}{3\sqrt{3}a_{0}}(0,-1)^{T}$, and
$a_{0}=1.42\textrm{ \AA}$ is the intra-layer distance between carbon atoms.
Also, note that $T_{n}^{\dagger}=T_{n}$. $c^{\dagger}_{\bm{k}}$ is a
16-components creation operator, whose flavors contains spin, valley, layer,
and sublattice.

Figure 6:  Here we plot comparisons between the original BM Hamiltonian and
its small-angle approximated Hamiltonian with varying $\kappa=w_{0}/w_{1}$ at
$\theta=1.05$, $v_{0}=2700\,\textrm{meV}$, and $w_{1}=105\,\textrm{meV}$. (a)
The discrepancy of bandwidth between the original TBG model and approximated
chiral TBG model. Here we also plot the gap between the conduction band and
second conduction (remote) band. Comparing left and right y-axes, we notice
that the flat bands are well-separated from remote bands in the wide range of
$\kappa$. (b) Wavefunction overlaps between the original BM model and small-
angle approximated BM model for valence bands, valence/conduction bands, and
flatband subspaces. As a whole, flat-band subspaces between two Hamiltonians
are almost identical until the gap to remote bands close around $\kappa=1$.

To examine the symmetry properties more closely, we apply the following basis
transformation:

$U_{g}=e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/4},\qquad c_{\bm{k}}\mapsto
U_{g}c_{\bm{k}},\qquad H({\bm{k}})\mapsto
U_{g}H({\bm{k}})U_{g}^{\dagger}=\hbar
v_{F}\quantity(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y}).$ (34)

Under this transformation, we can remove the explicitly $\theta$-dependent
phase factor in $H({\bm{k}})$, while $T({\bm{k}}^{\prime},{\bm{k}})$ is
slightly modified, resulting in

$\displaystyle\hat{H}$
$\displaystyle=\sum_{{\bm{k}}}c^{\dagger}_{\bm{k}}H({\bm{k}})c_{\bm{k}}+\sum_{{\bm{k}},{\bm{k}}^{\prime}}c^{\dagger}_{{\bm{k}}^{\prime}}T({\bm{k}}^{\prime},{\bm{k}})c^{\,}_{\bm{k}}$
$\displaystyle H({\bm{k}})$ $\displaystyle=\hbar
v_{F}\quantity(k_{x}\sigma_{x}\tau_{z}+k_{y}\sigma_{y})$ $\displaystyle
T({\bm{k}}^{\prime},{\bm{k}})$
$\displaystyle=\frac{1}{2}\sum_{n}\delta_{{\bm{k}}^{\prime},{\bm{k}}-\tau_{z}{\bm{q}}_{n}}\quantity(\mu_{x}-i\mu_{y})\tilde{T}_{n}+\textrm{h.c.}$
(35)

where

$\tilde{T}_{n}=w_{0}e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}+w_{1}e^{2\pi
in\sigma_{z}\tau_{z}/3}\sigma_{x}e^{-2\pi in\sigma_{z}\tau_{z}/3}.$ (36)

Now, consider a many-body basis transformation

$\displaystyle{\cal C}:\quad c^{\dagger}_{\bm{k}}$
$\displaystyle\mapsto\tau_{x}\sigma_{x}\mu_{y}c^{\,}_{\bm{k}}$ $\displaystyle
i$ $\displaystyle\mapsto-i$ (37)

which is a proper symmetry when $w_{0}=0$. The operator ${\cal{C}}$ maps
$w_{0}e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}$ to
$w_{0}e^{i\theta\sigma_{z}\mu_{z}\tau_{z}/2}$. For a small $\theta$, the
discrepancy is indeed very small. If we ignore $\theta$-dependence in the
transformed Hamiltonian, i.e., setting $\tilde{T}_{n}$ as

$\tilde{T}_{n}=w_{0}+w_{1}e^{2\pi in\sigma_{z}\tau_{z}/3}\sigma_{x}e^{-2\pi
in\sigma_{z}\tau_{z}/3},$ (38)

enhances the approximate symmetrz ${\cal{C}}$ to an exact one. To test the
validity of this approximation, we examined the wavefunction overlap between
eigenstates of the actual Hamiltonian and approximated Hamiltonian as function
of $w_{0}$ and ${\bm{k}}$ in Fig. 6. As a whole, the flat-band subspace of the
two Hamiltonians are almost identical until the gap to remote bands close
around $\kappa=1$. Individually, the valence band overlap between two
Hamilonians drops around $\kappa\sim 0.25$ and $\kappa\sim 0.85$. The
discrepancy at $\kappa\sim 0.25$ is physically not very important since the
bandwidth illustrated in (a) is extremely flat, implying that valence and
conduction bands are almost degenerate. Therefore, compared to the interaction
energy scale, this discrepancy can be ignored. The discrepancy at $\kappa\sim
0.85$ might have a physical significance, since this implies that the valence
and conduction bands hybridize at certain ${\bm{k}}$-points other than at K
and $\textbf{K}^{\prime}$. Therefore, we confirm that the approximation works
well for $w_{0}<0.8$.

##  Appendix C Hartree-Fock calculation of the bare dispersion

In this section, we review the Hartree-Fock calculation, which is used to
calculate the bare dispersion $\hat{H}_{0}$ In the BM model Hamiltonian, we
obtain the single-particle dispersion under the assumption that the
interaction does not affect the single-particle physics. However, note that
its tight-binding parameters are extracted from first principle calculations
and scattering experiments, where interaction effects are taken into account.
Therefore, the obtained BM model dispersion can be considered as an effective
bandwidth resulting from a bare quadratic terms $\hat{H}_{0}$ and the Coulomb
interaction $\hat{V}$ [59, 27]. In other words,

$\hat{H}_{\textrm{BM}}=\quantity[\hat{H}_{0}+\hat{V}\,]_{\Psi},\qquad\hat{V}=\frac{1}{2}\sum_{\bm{q}}\hat{\rho}_{\bm{q}}V_{\bm{q}}\hat{\rho}_{-{\bm{q}}}$
(39)

for some ground-state wavefunction $\Psi$. If we want to simulate the full-
interacting Hamiltonian, we should use the bare dispersion $\hat{H}_{0}$
instead of naively using $\hat{H}_{\textrm{BM}}$, to avoid double-counting. By
choosing a proper reference state, one can calculate the bare dispersion as

$\hat{H}_{0}=\hat{H}_{\textrm{BM}}-\big{[}\hat{V}\,\big{]}_{\Psi}.$ (40)

When $\Psi$ is a Slater state, its electronic correlation is fully captured by
a single-particle density matrix $P_{0}$ and we can calculate the mean-field
contribution of the interaction $\hat{V}$ with respect to $P_{0}$ using the
Hartree-Fock method. In the following, we elaborate on the details of the
Hartree-Fock calculation and how the bare dispersion $\hat{H}_{0}$ and
interaction-only spectrum $[\hat{V}]$ can be obtained.

Figure 7:  The bandwidths of electron-hole excitation spectra for bare
dispersion $\hat{H}_{0}$ and the interaction-only model at (a)
$\theta=1.05^{\circ}$, (b) $\theta=1.1^{\circ}$, and (c) $\theta=1.2^{\circ}$
under the subtraction scheme in Eq. (C). Here we use $\frac{\hbar
v_{f}}{3a_{0}/2}=2700\,\textrm{meV}$, $w_{0}=105\,\textrm{meV}$, and $6\times
6$ system size. The result shows that the bare dispersion near the BM-model
magic angle $\theta=1.05^{\circ}$ is not minimal. Rather, the bare dispersion
is suppressed at a slightly higher angle of $1.1^{\circ}$, while at
$1.2^{\circ}$, it increases back. Therefore, the so-called ‘magic angle’s for
the BM dispersion and bare dispersion are somewhat different. Around the twist
angle $1.2^{\circ}$, the bare dispersion is significantly larger than that of
$\theta=1.05^{\circ}$ throughout the whole range of $\kappa$. Based on this
figure, we infer that the ground-state at this larger angle would be semi-
metallic even at $\kappa=0.75$.

Let $P_{0}$ be the one-body projector for the slater reference state
$\Psi_{0}$. Following the standard Hartree-Fock calculations and projection to
the flat bands [102], the bare Hamiltonian should be given as the following:

$\displaystyle\big{[}\hat{V}\big{]}_{P_{0}}$
$\displaystyle=\sum_{\bm{k}}\sum_{mn\in{\cal
I}}c^{\dagger}_{{\bm{k}},m}\Big{(}\,\xi^{H}_{{\bm{k}},mn}[P_{0}]+\xi^{F}_{{\bm{k}},mn}[P_{0}]+\xi^{\textrm{Norm}}_{{\bm{k}},mn}-\xi^{H}_{{\bm{k}},mn}[\Delta]-\xi^{F}_{{\bm{k}},mn}[\Delta]\Big{)}c^{\,}_{{\bm{k}},n}$
$\displaystyle\hat{H}_{0}$
$\displaystyle=\hat{H}_{\textrm{BM}}-\sum_{\bm{k}}\sum_{mn\in{\cal
I}}c^{\dagger}_{{\bm{k}},m}\Big{(}\,\xi^{H}_{{\bm{k}},mn}[P_{0}-\Delta]+\xi^{F}_{{\bm{k}},mn}[P_{0}-\Delta]+\xi^{\textrm{Norm}}_{{\bm{k}},mn}\Big{)}c^{\,}_{{\bm{k}},n}$
(41)

where $\cal I$ is the set of indices of the flat-bands. The second equality
comes from the fact that Hartree and Fock terms are linear in the projector
matrix, as in Eq. (C). Here, the first two terms are from the HF contributions
of the reference Slater state represented by $P_{0}$, and the third term is
from the normal-ordering. Note that, for the Hartree-Fock calculation, one
should first normal order the term
$\frac{1}{2}\sum_{\bm{q}}V_{\bm{q}}\rho_{\bm{q}}\rho_{-{\bm{q}}}=\frac{1}{2}\sum_{\bm{q}}V_{\bm{q}}:\rho_{\bm{q}}\rho_{-{\bm{q}}}:+\xi^{\textrm{Norm}}$.
The last two terms stem from the projected-out occupied bands, whose one-body
projector is defined as

$\Delta_{ij}=\begin{cases}\delta_{ij}&\quad\textrm{if }i,j\in\bar{\cal
I}_{o}\\\ 0&\quad\textrm{otherwise}\end{cases}$ (42)

where $\bar{\cal I}$ is the set of indices of the projected-out remote bands,
and the subscript $o$ (or $u$) are for occupied (or unoccupied) bands. Each
contribution is defined as the following:

$\displaystyle\xi^{H}_{{\bm{k}},mn}[M]$
$\displaystyle=\sum_{{\bm{G}}}[\Lambda^{\dagger}({\bm{k}},{\bm{G}})]_{mn}\sum_{{\bm{k}}^{\prime}}\tr(M^{T}({\bm{k}}^{\prime})\cdot\Lambda({\bm{k}}^{\prime},{\bm{G}}))$
$\displaystyle\xi^{F}_{{\bm{k}},mn}[M]$
$\displaystyle=-\sum_{\bm{q}}V_{\bm{q}}\quantity[\Lambda^{\dagger}({\bm{k}},{\bm{q}})M^{T}({\bm{k}}+{\bm{q}})\Lambda({\bm{k}},{\bm{q}})]_{mn}$
$\displaystyle\xi^{\textrm{Norm}}_{{\bm{k}},mn}$
$\displaystyle=\frac{1}{2}\sum_{\bm{q}}V_{{\bm{q}}}\quantity[\Lambda^{\dagger}({\bm{k}},{\bm{q}})\Lambda({\bm{k}},{\bm{q}})]_{mn}=-\frac{1}{2}\xi^{F}_{{\bm{k}},mn}[\mathbb{I}_{\cal
I}]$ (43)

where we used the fact that $V_{\bm{q}}=V_{-{\bm{q}}}$, and the indices for
$\Lambda$ run over ${\cal I}\cup\bar{\cal I}$. Here, $\mathbb{I}_{\cal I}$ is
the matrix which is the identity on the flat-band block and zero on the other.

As discussed in [102], the form of the bare Hamiltonian depends on what forms
of interaction we use. In Eq. (C), the interaction is given as the
$\sim\rho_{\bm{q}}\rho_{-{\bm{q}}}$, which is the raw density-density form,
where the density is measured from the empty filling, not the charge
neutrality. For our purposes, let us use the density $\bar{\rho}$ which is
measured with respect to the charge neutrality and choose the interaction
accordingly, i.e.,
$\hat{V}\sim(\rho_{\bm{q}}-\rho^{0}_{\bm{q}})(\rho_{-{\bm{q}}}-\rho^{0}_{-{\bm{q}}})=\bar{\rho}_{\bm{q}}\bar{\rho}_{-{\bm{q}}}$.
With this modification, the bare Hamiltonian is calculated as the following:

$\displaystyle\hat{H}_{0}$
$\displaystyle=\hat{H}_{\textrm{BM}}-\sum_{\bm{k}}\sum_{mn\in{\cal
I}}c^{\dagger}_{{\bm{k}},m}\Big{(}\,\xi^{H}_{{\bm{k}},mn}[P_{0}-\Delta]+\xi^{F}_{{\bm{k}},mn}[P_{0}-\Delta]+\xi^{\textrm{Norm}}_{{\bm{k}},mn}-\xi^{\textrm{back}}_{{\bm{k}},mn}\Big{)}c^{\,}_{{\bm{k}},n}$
$\displaystyle=\hat{H}_{\textrm{BM}}-\sum_{\bm{k}}\sum_{mn\in{\cal
I}}c^{\dagger}_{{\bm{k}},m}\quantity(\,\xi^{H}_{{\bm{k}},mn}\quantity[P_{0}-\Delta-\frac{1}{2}\mathbb{I}_{\cal
I}]+\xi^{F}_{{\bm{k}},mn}\quantity[P_{0}-\Delta-\frac{1}{2}\mathbb{I}_{\cal
I}])c^{\,}_{{\bm{k}},n}$ $\displaystyle\xi^{\textrm{back}}$
$\displaystyle=\frac{1}{2}\sum_{\bm{q}}V_{\bm{q}}\quantity(\rho_{{\bm{q}}}\rho^{0}_{-{\bm{q}}}+\rho^{0}_{{\bm{q}}}\rho_{-{\bm{q}}})=\frac{1}{2}\xi^{H}[\mathbb{I}_{\cal
I}]$ (44)

Interestingly, quadratic terms originating from the normal ordering and
background charge density can be incorporated into the Hartree and Fock terms
as above. Furthermore, in the case where
$P_{0}-\Delta=\frac{1}{2}(\mathbb{I}+Q)$ is particle-hole symmetric, i.e., for
first quantized expression of PH symmetry ${\cal{C}}$ in Eq. (37),
$\textrm{(PH)}\,Q\,\textrm{(PH)}^{-1}=-Q$, and where
$\textrm{(PH)}\,\Lambda({\bm{k}},{\bm{q}})\,\textrm{(PH)}^{-1}=\Lambda({\bm{k}},{\bm{q}})$,
one can prove that $\xi^{H}$ contribution vanishes. Since this is the case
under small-angle approximation, we are left with only the Fock contribution.
Therefore, in this case our total Hamiltonian $\hat{H}=\hat{H}_{0}+\hat{V}$ is
defined as the following:

$\displaystyle\hat{H}_{0}$
$\displaystyle=\hat{H}_{\textrm{BM}}-\sum_{\bm{k}}\sum_{mn\in{\cal
I}}\xi^{F}_{{\bm{k}},mn}\Big{[}P_{0}-\Delta-\mathbb{I}_{\cal I}/2\Big{]}\cdot
c^{\dagger}_{{\bm{k}},m}c^{\,}_{{\bm{k}},n}$ $\displaystyle\hat{V}$
$\displaystyle=\frac{1}{2}\sum_{\bm{q}}\bar{\rho}_{\bm{q}}V_{\bm{q}}\bar{\rho}_{-{\bm{q}}}-\frac{1}{2}\sum_{\bm{G}}V_{\bm{G}}\rho^{0}_{\bm{G}}\rho^{0}_{-{\bm{G}}},\qquad\bar{\rho}_{\bm{q}}\equiv\rho_{\bm{q}}-\frac{1}{2}\sum_{{\bm{k}},{\bm{G}}}\delta_{{\bm{q}},{\bm{G}}}\tr_{\cal
I}\Lambda({\bm{k}},{\bm{G}})$ (45)

This formulation implies that, if we naively use $\hat{H}_{\textrm{BM}}$ for
the bare dispersion, it is equivalent to take
$P_{0}-\Delta\equiv\mathbb{I}_{\cal I}/2$. Physically, such a choice
corresponds to the case where the reference state is an infinite temperature
ensemble.

Now, we are ready to calculate the bare dispersion following the Eq. (C) with
a proper choice of $P_{0}$. In [59], one takes $P_{0}$ to be the one-body
projector for the groundstate of the BM Hamiltonian. In [27, 30], the
reference state $\Psi_{0}$ is taken to be the groundstate of decoupled bilayer
graphenes. In the second approach, the underlying assumption is that tight-
binding parameters including Dirac velocity are extracted for isolated
graphene layer, which will be renormalized under the interlayer coupling. Here
we take the second approach, where we only considered the flat band subspace
and ignored the contribution from remote bands. The results are illustrated in
Fig. 7, where we plot the bandwidth $W$ of $\hat{H}_{0}$ as function of
$\kappa=w_{0}/w_{1}$ for three different twist angles
$\theta=1.05,\,1.1,\,1.2$ (Red crosses). We also show the excitation spectrum
of the strong coupling limit (interaction-only with $\hat{H}=\hat{V}$) at
charge neutrality. Since any polarized Slater state listed in Tab. 2 is the
exact groundstate in this limit, we can chose the valley-polarized state and
calculated the exact electronic excitation spectrum using the following
equation:

$\displaystyle\big{[}\,\hat{V}\,\big{]}_{P}$
$\displaystyle=\sum_{\bm{k}}\sum_{mn\in{\cal
I}}\xi^{F}_{{\bm{k}},mn}\Big{[}P-\Delta-\mathbb{I}_{\cal I}/2\Big{]}\cdot
c^{\dagger}_{{\bm{k}},m}c^{\,}_{{\bm{k}},n}$ (46)

where we use $P_{ij}=1$ if $i=j$ are in $+$ valley, and 0 otherwise.

##  Appendix D Details of the Numerical Steps

###  D.1 Quantum Monte Carlo Hamiltonian

In the DQMC, we simulate the following Hamiltonian in momentum space:

$\displaystyle\hat{H}$ $\displaystyle=\hat{H}_{0}+\hat{V}$
$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{s,\tau}\quantity[\sum_{n,m}H^{0}_{\tau,{\bm{k}},nm}c^{\dagger}_{{\bm{k}},s,\tau,n}c^{\,}_{{\bm{k}},s,\tau,m}]$
$\displaystyle\hat{V}$
$\displaystyle=\frac{1}{8}\sum_{\bm{q}}V_{\bm{q}}\quantity[(\bar{\rho}_{\bm{q}}+\bar{\rho}_{-{\bm{q}}})^{2}-(\bar{\rho}_{\bm{q}}-\bar{\rho}_{-{\bm{q}}})^{2}]$
(47)

where $\hat{H}_{0}$ is defined in Eq. (40), which is block-diagonal in spin
and valley flavors. Here the Fourier transformation of the interaction
$V_{\bm{q}}$ is defined as the following for a naive screened Coulomb
interaction with screening length $d$:

$V_{\bm{q}}=\frac{1}{A}\int\differential^{2}r\,e^{i{\bm{q}}r}\quantity[\frac{e^{2}}{4\pi\epsilon_{0}\epsilon
r}e^{-\absolutevalue{r}/d_{s}}]=\frac{V_{0}d}{\sqrt{1+\absolutevalue{{\bm{q}}}^{2}d^{2}}},\qquad
V_{0}=\frac{e^{2}}{2\epsilon\epsilon_{0}}\frac{1}{A}=\frac{e^{2}}{\sqrt{3}\epsilon\epsilon_{0}a_{M}^{2}}\frac{1}{N_{M}}$
(48)

where $A=N_{M}A_{0}$ is the total area of the system, $A_{0}$ is the area of a
moiré unit cell, and $N_{M}$ is the number of moiré unit cells in the systems.
$A_{0}$ is given as $A_{0}=\sqrt{3}a_{M}^{2}/2$, where
$\absolutevalue{a_{M}}=\sqrt{3}a/2\sin\frac{\theta}{2}=13.4\textrm{ nm}$ at
$\theta=1.05^{\circ}$ with $a=1.42\,\textrm{\AA}$. For realistic model, one
can consider a more accurate description that involves gate screening, but
they are almost the same in their functional forms. For the doubly gated case
with gate distance $d$, we have the following screened Coulomb interaction in
the momentum space:

$V_{\bm{q}}=\frac{V_{0}\tanh(\absolutevalue{{\bm{q}}}d)}{\absolutevalue{{\bm{q}}}},\qquad
V_{0}=\frac{e^{2}}{2\epsilon\epsilon_{0}}\frac{1}{A}=\frac{e^{2}}{\sqrt{3}\epsilon\epsilon_{0}a_{M}^{2}}\frac{1}{N_{M}}$
(49)

At ${\bm{q}}=0$, we obtain the charging energy
$V_{\textrm{CE}}=\frac{1}{2}V_{0}d$. This is important as the value would
contribute to the spectral function even for a gapless excitations. In the
parameter we use, note that
$V_{\textrm{CE}}\sim\frac{58}{L^{2}}\,\textrm{meV}$. At $L=9$, the charging
energy is 0.7 meV.

###  D.2 Symmetry

Once we project onto the flat-band subspace, there are total 8 flavors of
electrons: spin $s_{z}=\pm 1$, valley $\tau_{z}=\pm 1$, and band $n_{z}=\pm
1$. Under the absence of dispersion and interaction, the system is invariant
under $\mathrm{U}(8)$ symmetry, i.e., for ${\cal U}\in\mathrm{U}(8)$,
$\hat{H}$ is invariant under $d_{\bm{k}}\mapsto{\cal U}d_{\bm{k}}$. There are
64 generators that are given as $s_{\mu}\tau_{\nu}n_{\lambda}$, with
$\mu,\nu,\lambda=0,x,y,z$. Of course, in realstic system this should reduce
into
$\mathrm{U}(2)_{K}\times\mathrm{U}(2)_{K^{\prime}}\sim\mathrm{U}(1)_{\textrm{el}}\times\mathrm{U}(1)_{\textrm{valley}}\times\mathrm{SU}(2)_{K}\times\mathrm{SU}(2)_{K^{\prime}}$.
As we switch the dispersion and interaction back on, the symmetry is reduced,
as elaborated in [27].

Symmetries | Microscopic | Sublattice Basis | Band Basis I | Band Basis II
---|---|---|---|---
${\cal T}$ (${\bm{k}}\mapsto-{\bm{k}}$) | $\tau_{x}\mathcal{K}$ | $\tau_{x}\mathcal{K}$ | $\tau_{x}n_{z}\mathcal{K}$ | $\tau_{x}\mathcal{K}$
$C_{2}$ (${\bm{k}}\mapsto-{\bm{k}}$) | $\tau_{x}\sigma_{x}$ | $\tau_{x}\sigma_{x}e^{i\theta({\bm{k}})}$ | $\tau_{x}$ | $\tau_{x}$
$C_{2}\mathcal{T}$ | $\sigma_{x}\mathcal{K}$ | $e^{i\theta({\bm{k}})}\sigma_{x}\mathcal{K}$ | $n_{z}\mathcal{K}$ | $\mathcal{K}$
Small-angle Approx. ${\cal{C}}$ | $\tau_{x}\sigma_{x}\mu_{y}$ | $\tau_{y}\sigma_{y}$ | $\tau_{x}n_{x}$ | $\tau_{x}n_{y}$
Sublattice ${\cal{C}}_{\textrm{sub}}$ (at $\kappa=0$) | $\sigma_{z}$ | $\sigma_{z}$ | $\tau_{z}n_{x}$ | $\tau_{z}n_{y}$
Table 3:  Representations of symmetry operators in different basis. Sublattice
basis is constructed in [27], and the band basis I is explicitly constructed
in the discussion. The band basis II can be obtained for a different gauge
choice. Among these operators, $\mathcal{T}$ and $C$ flip the momentum
${\bm{k}}$. Note that the above choice is not unique. They should satisfy the
following relations:
$[\mathcal{T},C_{2}\mathcal{T}]=[\mathcal{T},{\cal{C}}_{\textrm{sub}}]=\\{\mathcal{T},{\cal{C}}\\}=0$,
$\\{C_{2}\mathcal{T},{\cal{C}}_{\textrm{sub}}\\}=\\{C_{2}\mathcal{T},{\cal{C}}\\}=0$,
$\\{{\cal{C}},{\cal{C}}_{\textrm{sub}}\\}=0$. Furthermore, in the band basis,
chiral symmetries ${\cal{C}}_{\textrm{sub}}$ and ${\cal{C}}$ should
anticommute with $n_{z}$. Depending on the detailed setting of the simulation,
we may use one basis over the other. Note that ${\cal{C}}_{\textrm{sub}}$ does
not have a proper representation on the flat-band subspace if we are away from
the chiral limit. Therefore, its representation would be only valid at
$\kappa=0$ limit.

0\. _Preliminary_ Before going into detail, let us elaborate on the symmetry
representation of a given basis. For simplicity, we ignore spin for now, which
reduces the number of flavors from 8 to 4. In [27], the sublattice basis was
introduced, but here we choose the band basis. In this discussion, we focus on
the following symmetries: (1) anti-unitary particle-hole symmetry (sublattice
symmetry) ${\cal{C}}_{\textrm{sub}}$, which exists only at $\kappa=0$, (2)
anti-unitary particle-hole symmetry ${\cal{C}}$ from small-angle
approximation, whose definition is given in Eq. (2), (3) time-reversal
$\mathcal{T}$, and (4) composition of $C_{2}$ rotation and time-reversal
$C_{2}\mathcal{T}$. It is important to emphasize that here we examine bloch
wavefunctions, and symmetry operators should be discussed in the first-
quantized language. As a result, ${\cal{C}}$ and ${\cal{C}}_{\textrm{sub}}$
are _unitary_ operators in this context.

For a given representation of symmetries in the band basis, the microscopic
action (_right action_) of a unitary symmetry operator ${\cal O}$ becomes the
basis action (_left action_) as the following:

${\cal O}u_{i}({\bm{k}})={\cal O}\matrixquantity(C^{i}_{\bm{k}}\\\
C^{i}_{{\bm{k}}+{\bm{G}}_{1}}\\\ C^{i}_{{\bm{k}}+{\bm{G}}_{2}}\\\
\vdots)=u_{j}({\bm{k}})\bar{{\cal O}}_{ji}$ (50)

where $u_{i}({\bm{k}})$ is $i$-th eigenvector at momentum ${\bm{k}}$. With
this understanding, the form factor $\Lambda$ and operator-inserted form
factor $\Lambda_{i}$ should satisfy

$\displaystyle\Lambda({\bm{k}},{\bm{q}})$
$\displaystyle=U^{\dagger}({\bm{k}})U({\bm{k}}+{\bm{q}})=U^{\dagger}({\bm{k}})\quantity[{\cal
O}^{-1}{\cal O}]U({\bm{k}}+{\bm{q}})=\bar{{\cal
O}}^{\dagger}\Lambda({\bm{k}},{\bm{q}})\bar{{\cal O}}$
$\displaystyle\Lambda_{i}({\bm{k}},{\bm{q}})$
$\displaystyle=U^{\dagger}({\bm{k}})M_{i}U({\bm{k}}+{\bm{q}})=U^{\dagger}({\bm{k}})\quantity[{\cal
O}^{-1}M_{j}{\cal O}]U({\bm{k}}+{\bm{q}})=\bar{{\cal
O}}^{\dagger}\Lambda_{j}({\bm{k}},{\bm{q}})\bar{{\cal O}}$ (51)

where $U({\bm{k}})$ is now a horizontally stacked eigenvectors, i.e., $N\times
4$ matrix, where $N$ is the number of microscopic degrees of freedom (layer,
valley, sublattice, and moiré reciprocal lattice vectors) and there are valley
and sublattice degrees of freedom. Here we assume that $M_{i}={\cal
O}^{-1}M_{j}{\cal O}$. Now, assume ${\cal O}$ is an anti-unitary operator,
given as ${\cal O}={\cal U}\mathcal{K}$, where ${\cal U}$ is the unitary part.
Assume that this anti-unitary operator has a unitary representation on basis,
i.e., ${\cal O}u_{i}({\bm{k}})=u_{j}({\bm{k}})\bar{{\cal O}}_{ji}$. Since
$\cal O$ is anti-unitary, ${\cal O}\quantity[\alpha
u_{i}({\bm{k}})]=\alpha^{*}u_{j}({\bm{k}})\bar{{\cal O}}_{ji}$. Then,

$\displaystyle\Lambda^{*}({\bm{k}},{\bm{q}})$
$\displaystyle=U^{T}({\bm{k}})U^{*}({\bm{k}}+{\bm{q}})=U^{T}({\bm{k}})\cdot{\cal
U}^{\dagger}{\cal U}\cdot U^{*}({\bm{k}}+{\bm{q}})$
$\displaystyle=\quantity[\,{\cal
U}\mathcal{K}U({\bm{k}})]^{\dagger}\quantity[\,{\cal
U}\mathcal{K}U({\bm{k}}+{\bm{q}})]=\bar{{\cal
O}}^{\dagger}\Lambda({\bm{k}},{\bm{q}})\bar{{\cal O}}$ (52)

Under this band basis, we can choose the gauge such that these have the
following representations: ($i$)
$\mathcal{T}U({\bm{k}})=U(-{\bm{k}})\tau_{x}n_{z}$, ($ii$)
${\cal{C}}U({\bm{k}})=U({\bm{k}})n_{x}\tau_{x}$, and ($iii$)
$(C_{2}\mathcal{T})U({\bm{k}})=U({\bm{k}})n_{z}$. Two conditions that have to
be enforced for the band basis are the followings: (i) $[n_{z},H]=0$ since $H$
is diagonal in the band basis, and (ii) $n_{z}$ anti-commutes with particle-
hole symmetries ${\cal{C}}$ and ${\cal{C}}_{\textrm{sub}}$. To be a valid
symmetry representation, symmetries should satisfy the following commutation
relations originated from the microscopic representation: ${\cal{C}}{\cal
T}=-{\cal T}{\cal{C}}$, ${\cal{C}}(C_{2}{\cal T})=-(C_{2}{\cal T}){\cal{C}}$
and ${\cal T}(C_{2}{\cal T})=(C_{2}{\cal T}){\cal T}$. In Tab. 3, we listed a
few different representations of symmetry operators which satisfy this set of
commutation relations. In this paper, we exclusively use the band basis I in
Tab. 3, but note that this is not unique; for comparison, we show another band
basis representation in Tab. 3.

In the following, let us fix the basis by choosing the representation in Band
Basis I. Using these symmetries, one can deduce the matrix structure of the
projected Hamiltonian and the form factors. The Hamiltonian is diagonal in
valley such that only terms of the form $\tau_{0,z}n_{0,x,y,z}$ are allowed.
Then, $\quantity{{\cal{C}},H}=0$ restricts them further to $\tau_{0}n_{y,z}$
and $\tau_{z}n_{0,x}$. Out of these four terms, $C_{2}\mathcal{T}$ forbids
$\tau_{z}n_{x}$. Therefore, we get

$h({\bm{k}})=h_{0}({\bm{k}})\tau_{z}+h_{y}({\bm{k}})n_{y}+h_{z}({\bm{k}})n_{z}$
(53)

Obtaining form factor is a bit different. In this case, we use Eq. (D.2) and
Eq. (D.2). Also, the coefficients in front of the Pauli matrices can be
complex. We get

$\Lambda_{0}({\bm{k}},{\bm{q}})\sim
C_{0}+C_{1}n_{x}+C_{2}\tau_{z}n_{y}+C_{3}\tau_{z}n_{z}$ (54)

where $C_{0},C_{2},C_{3}$ are real while $C_{1}$ is imaginary. In the chiral
limit, $\kappa=0$, the sublattice chiral symmetry
${\cal{C}}_{\textrm{sub}}=\tau_{z}n_{x}$ further constrains $h({\bm{k}})$ and
$\Lambda_{0}({\bm{k}},{\bm{q}})$, to $h_{0}({\bm{k}})=0$ and $C_{2}=C_{3}=0$.
Now that we know the structure of the full interacting Hamiltonian depending
on the presence of chiral symmetry ${\cal{C}}_{\textrm{sub}}$, we can
understand the global symmetry of the system as follows.

1\. _Chiral limit_ ($\kappa=0$): $\hat{H}_{\textrm{kin}}=0$ and
$\hat{H}_{\textrm{int}}\neq 0$: In this limit, as $\Lambda_{\bm{q}}$ is given
as the linear combination of $\mathbb{I}$ and $n_{x}$, generators of
$\mathrm{U}(8)$ that commute with $n_{x}$ would still be valid symmetry
generators, which halves the number of valid generators from 64 to 32. Those
32 generators can be divided into two sets of 16, each of which generate a
$\mathrm{U}(4)$ symmetry. Therefore, the system has
$\mathrm{U}(4)_{+}\times\mathrm{U}(4)_{-}$ symmetry, where the subscript $\pm$
corresponds to the eigenvalue of $n_{x}$. In the sublattice basis [27], the
operator $n_{x}$ corresponds to $\tau_{z}\sigma_{z}$ which is the Chern
number.

2\. _Chiral limit_ ($\kappa=0$): $\hat{H}_{\textrm{kin}}\neq 0$ and
$\hat{H}_{\textrm{int}}\neq 0$: In this case, a non-zero dispersion is given
by $h({\bm{k}})=h_{y}({\bm{k}})n_{y}+h_{z}({\bm{k}})n_{z}$. Therefore,
symmetry generators must commute with $n_{y}$ in addition to $n_{x}$. This
will again halve the number of symmetry generators from 32 to 16. As a result,
the system has a single $\mathrm{U}(4)_{\textrm{chiral}}$ symmetry with 16
generators, which basically acts trivially on the band index and non-trivially
on the spin-valley indices.

3\. _Non-Chiral limit_ ($\kappa\neq 0$): $\hat{H}_{\textrm{kin}}=0$ and
$\hat{H}_{\textrm{int}}\neq 0$: In this case, the form factor take the
following general form of Eq. (54). As a result, the interaction Hamiltonian
is invariant only under generators that commute with $n_{x}$ and
$\tau_{z}n_{y}$, which decreases the number of valid generators from 64 to 16.
Therefore, the system has a $\mathrm{U}(4)_{\textrm{flat}}$ symmetry. Note
that this $\mathrm{U}(4)_{\textrm{flat}}$ is different from
$\mathrm{U}(4)_{\textrm{chiral}}$ since while the generators of flat one
commute with $n_{y}$, those of the chiral one commute with $\tau_{z}n_{y}$.