# Theory of correlated Chern insulators in twisted bilayer graphene

Xiaoyu Wang<EMAIL_ADDRESS>National High Magnetic Field Lab,
Tallahassee, FL 32310  Oskar Vafek<EMAIL_ADDRESS>National High
Magnetic Field Lab, Tallahassee, FL 32310  Department of Physics, Florida
State University, Tallahassee, FL 32306

###### Abstract

Magic-angle twisted bilayer graphene is the best studied physical platform
featuring moiré potential induced narrow bands with non-trivial topology and
strong electronic correlations. Despite their significance, the Chern
insulating states observed at a finite magnetic field –and extrapolating to a
band filling, $s$, at zero field– remain poorly understood. Unraveling their
nature is among the most important open problems in the province of moiré
materials. Here we present the first comprehensive study of interacting
electrons in finite magnetic field while varying the electron density, twist
angle and heterostrain. Within a panoply of correlated Chern phases emerging
at a range of twist angles, we uncover a unified description for the
ubiquitous sequence of states with the Chern number $t$ for
$(s,t)=\pm(0,4),\pm(1,3),\pm(2,2)$ and $\pm(3,1)$. We also find correlated
Chern insulators at unconventional sequences with $s+t\neq\pm 4$, as well as
with fractional $s$, and elucidate their nature.

The twisted bilayer graphene (TBG) has been a subject of intense theoretical
and experimental investigation, in no small part due to its isolated,
topologicaly non-trivial, narrow bands displaying rich correlated electron
physics when partially occupied [1, 2]. As the twist angle between the two
graphene layers is tuned toward the magic value of $\sim 1.05^{\circ}$ [3],
TBG devices show a plethora of correlated phenomena including
superconductivity, correlated insulating states and (quantized) anomalous Hall
effect [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22,
23, 24]. The non-trivial topology of the pair of narrow bands for a given
valley and spin flavor is protected by the combined two-fold rotation symmetry
about the out-of-plane axis $C_{2z}$ (an emergent symmetry at low twist angle)
and spinless time reversal symmetry $T$ [25, 26]. The narrow band Hilbert
space can thus be decomposed into a Chern $+1$ and a Chern $-1$ band [25, 27,
28]. One way to reveal the non-trivial topology of the narrow bands is to
break $C_{2z}$ via alignment with the hexagonal boron nitride substrate (hBN)
and separate the Chern bands in energy. If in addition, the valley is
spontaneously polarized, thus breaking $T$, the resulting state with one
electron or hole per moiré unit cell becomes a Chern $\pm 1$ insulator [29,
30, 31, 32]. Indeed, experiments have observed anomalous Hall effect (AHE)
near the filling of $3$ electron per moiré unit cell [10, 19] in hBN aligned
samples. Further studies on non-aligned samples [20, 23] have also observed
AHE near 1 electron per moiré unit cell. Theoretically, such zero-field Chern
insulating states (zCIs) have been proposed to be energetically competitive
near magic angle, when the Coulomb interaction exceeds the narrow bandwidth,
even without the hBN alignment [33, 28, 34, 31].

An external magnetic field $B$, which preserves $C_{2z}$ but breaks $T$, has
been argued to be an alternative way to reveal the band topology [11, 18, 21],
as evidenced by the experimental observations of correlated Chern insulating
states (CCIs) with a finite Chern number $t$ and extrapolating to a band
filling $s$ at $B=0$ [6, 7, 8, 11, 13, 14, 15, 16, 17, 18, 20, 21, 22, 23,
24]. Specifically, the most prominent sequence of CCIs has $(s,t)=(0,\pm
4),\pm(1,3),\pm(2,2),\pm(3,1)$, consistent with selective population of the
aforementioned $B=0$ Chern $\pm 1$ bands [11, 20, 21]. These experiments also
report that some CCIs are stable at extremely low $B$, suggesting that they
originate from the zCIs [11, 20].

However, CCIs are also observed in TBG devices away from the magic angle
($\sim 1.27^{\circ}$), where the bandwidth of the narrow bands is expected to
be significantly larger than at $1.05^{\circ}$, without any observation of the
correlated insulators at $B=0$ [9]. Such CCIs appear only above a critical
$B$, below which they transition into nearly compressible states for a fixed
$(s,t)$. Similar phenomenology has also been reported in near-magic angle
devices, leading to an alternative explanation of the CCIs invoking Stoner
ferromagnetism within the magnetic subbands [15, 16, 22, 7], termed Hofstadter
subband ferromagnets (HSFs) [15]. As argued theoretically [35, 36, 37],
realistic heterostrain can also increase the bandwidth dramatically near the
magic angle, likely placing many TBG devices in the intermediate coupling
regime where the zCIs may not be energetically favored.

To date, the nature of these CCIs remains poorly understood. No microscopic
calculation favoring either zCI or HSF has been carried out at $B\neq 0$, nor
tying them to the relevant experiments. Moreover, the interplay of the CCIs
with the competing states at $B=0$ near the magic angle, such as the
intervalley coherent states (IVCs) [38, 27], the incommensurate Kekulé spiral
orders (IKSs) [36], and the striped and nematic states [28, 14, 39, 40],
remains unclear.

Here we report the first comprehensive study of the interacting electrons
within the TBG narrow bands directly at $B\neq 0$, and construct the phase
diagram for a range of twist angles, $B$-fields and electron densities, with
and without heterostrain. Consistent with the experimental observations, we
find CCIs with $(0,\pm 4)$, $\pm(1,3)$, $\pm(2,2)$, and $\pm(3,1)$, see Fig.
1. Such states are found to be stablized at higher $B$ fields for twist angles
as high as $1.38^{\circ}$ (highest twist angle studied in this work), and for
realistic heterostrain; based on analysis of their wavefunctions, we identify
them as correlated Hofstadter ferromagnets (CHFs). Similar to HSFs, they
correspond to selective population of the valley and spin flavors, but of the
interaction-renormalized magnetic subbands (see Fig. 2). Upon lowering $B$ we
find a first order phase transition into nearly compressible states at larger
twist angles, or into incompressible states with intervalley coherence closer
to the magic angle, at non-zero $s$. In the absence of heretostrain and at
larger twist angles we similarly observe CHFs transitioning into nearly
compressible states upon lowering $B$ (see Fig. S2). As we lower the twist
angle toward $1.05^{\circ}$, the incompressible $\pm(3,1)$ state extends to
lower $B$ and crosses over into the finite $B$ analog of the zCI, approaching
maximal sublattice polarization. We refer to this state as the strong coupling
Chern insulator (sCI). The $\pm(2,2)$, $\pm(1,3)$ and $(0,\pm 4)$ states also
cross over into sCIs upon decreasing the twist angle, but they experience a
first order phase transition into IVCs at low twist angles. The details of
this transition depend sensitively on the model parameters, as shown in Fig. 3
and Fig. S4. Because there is no symmetry distinction between them, the CHF
smoothly crosses over into the sCI as the twist angle is decreased toward the
magic angle, see Fig. 3 and Fig. S4. In addition to these prominent CCIs, we
also find Quantum Hall ferromagnetic states (QHFMs) emanating from band edges
and the charge neutrality point, CCIs with $s+t\neq\pm 4$, as well as
fractional $s$ that break the magnetic translation symmetry.

Figure 1:  Single-electron gap $\Delta$ as a function of electron filling and
magnetic flux at twist angles $1.38^{\circ}$ (a), $1.32^{\circ}$ (b),
$1.28^{\circ}$ (c), $1.24^{\circ}$ (d), $1.2^{\circ}$ (e), and $1.05^{\circ}$
(f). The size of the gaps are proportional to the radius of the respective
solid circles, with the circle representation $10$meV shown in (a). The
sequence of $(s,t)=(0,-4)$, $(-1,-3)$, $(-2,-2)$, and $(-3,-1)$ are labeled by
red dashed lines in (c). Chern insulating states with fractional $s$ at
$(s,t)=(-1/2,-3)$, $(-2/3,-3)$, and $(-3/2,-2)$ are labeled by orange dashed
lines in (f).

Model and method: We perform self-consistent Hartree-Fock analysis (B-SCHF) at
$B\neq 0$ using the continuum Bistritzer-MacDonald (BM) Hamiltonian, with
Coulomb interactions projected onto the narrow band Hilbert space. Here we
briefly outline the formalism, additional details are in the Supplementary
Information (SI). Starting with the (strained) BM Hamiltonian at rational
magnetic flux ratios $\phi/\phi_{0}=p/q$ where $p$ and $q$ are coprime
integers, $\phi$ is the magnetic flux per moiré unit cell and $\phi_{0}=h/e$
is the magnetic flux quantum, we solve for the Hofstadter spectra
$\varepsilon_{\eta sr}(\mathbf{k})$ and associated eigenstates
$\ket{\Psi_{\eta sr}(\mathbf{k})}$. Here $\eta=\mathbf{K},\mathbf{K}^{\prime}$
and $s=\uparrow,\downarrow$ are valley and spin quantum numbers, $r=1,\dots
2q$ is the magnetic subband index, and
$\mathbf{k}=k_{1}\mathbf{g}_{1}+k_{2}\mathbf{g}_{2}$ is the magnetic crystal
momentum defined in the magnetic Brillouin zone with $k_{1}\in[0,1)$ and
$k_{2}\in[0,1/q)$, and $\mathbf{g}_{i=1,2}$ are moiré reciprocal lattice
vectors (see SI). Spin Zeeman splitting is also considered in this
calculation. For the model parameters studied in this work, the gap to remote
Hofstadter bands does not close at the magnetic fluxes of interest. We study
interaction effects by projecting the screened Coulomb interaction onto the
narrow band Hilbert space. The Hamiltonian is given by:

$H=\sum_{\eta sr,\mathbf{k}}\varepsilon_{\eta sr}(\mathbf{k})d^{\dagger}_{\eta
sr,\mathbf{k}}d_{\eta
sr,\mathbf{k}}+\frac{1}{2A}\sum_{\mathbf{q}}V_{\mathbf{q}}\delta\hat{\rho}_{\mathbf{q}}\delta\hat{\rho}_{-\mathbf{q}}.$
(1)

Here $A$ is the total area of the system, $d_{\eta sr,\mathbf{k}}$ is the
electron annihilation operator, $\delta\hat{\rho}_{\mathbf{q}}$ is the Fourier
transform of the electron density operator projected onto the narrow bands,
subtracting a background charge density [41]. It is given by:

$\begin{split}\delta\hat{\rho}_{\mathbf{q}}&=\sum_{\eta
s,rr^{\prime},\mathbf{k},\mathbf{p}}\bra{\Psi_{\eta
sr}(\mathbf{k})}e^{-i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi_{\eta
sr^{\prime}}(\mathbf{p})}\\\ &\times\left(d^{\dagger}_{\eta
sr,\mathbf{k}}d_{\eta
sr^{\prime},\mathbf{p}}-\frac{1}{2}\delta_{r,r^{\prime}}\delta_{\mathbf{k},\mathbf{p}}\right).\end{split}$

We consider a dual-gate screened Coulomb interaction of the form
$V(\mathbf{q})=\frac{2\pi
e^{2}}{\epsilon_{0}\epsilon_{r}|\mathbf{q}|}\tanh\left(\frac{|\mathbf{q}|\xi}{2}\right)$,
with relative dielectric constant $\epsilon_{r}=15$ and screening length
$\xi=4\sqrt{|\mathbf{L}_{1}||\mathbf{L}_{2}|}$, where $\mathbf{L}_{i=1,2}$ are
the primitive moiré unit cell vectors (see SI). These parameters are chosen to
match the overall change of chemical potential from empty to full occupation
of the narrow bands in magic angle devices as extracted from the
compressibility measurements as well as STM [6, 42, 43, 7, 44] (see also Fig.
S7).

The B-SCHF calculation is carried out for a range of twist angles from
$1.38^{\circ}$ to $1.05^{\circ}$, both for unstrained model as well as for
realistic uniaxial heterostrain strength of $\epsilon=0.2\%$ and orientation
$\varphi=0^{\circ}$ (see SI for details). The respective non-interacting
Hofstadter spectra and Wannier diagrams are shown in Fig. S1.

The projected Hamiltonian at $B\neq 0$ is invariant under the following set of
symmetries [45, 41, 46]: $C_{2z}$, valley $U(1)_{v}$ and spin $U(1)_{s}$,
many-body particle-hole $P$ , magnetic translation symmetries generated by
$\hat{t}_{\mathbf{L}_{1}}$ and $\hat{t}_{\mathbf{L}_{2}}$ (see SI). In the
absence of heterostrain $C_{3z}$ and $C_{2y}T$ also leave $H$ invariant. $P$
guarantees symmetry about the charge neutrality point, and therefore we
present our results for the hole filling only.

Results: We first address the finite $B$ phase diagram for TBG subject to
$0.2\%$ of heterostrain. Fig. 1 gives an overview of the calculated single
particle excitation gap $\Delta$ as a function of moiré unit cell filling
($n/n_{s}$) and magnetic flux ratio ($\phi/\phi_{0}$) for six twist angles
$1.38^{\circ}$, $1.32^{\circ}$, $1.28^{\circ}$, $1.24^{\circ}$, $1.20^{\circ}$
and $1.05^{\circ}$. The size of the gap is proportional to the radius of the
solid circle. As seen, there is a rich panoply of correlated insulating
states. We start by focusing on the sequence of CCIs with
$(s,t)=(0,-4),(-1,-3),(-2,-2),(-3,-1)$, which are observed for all the twist
angles studied, and marked in red in Fig. 1(c). At larger twist angles, CCIs
along $(-3,-1),(-2,-2),(-1,-3)$ emerge at high $\phi/\phi_{0}$, and are
replaced by nearly compressible states at lower $\phi/\phi_{0}$ via a first
order phase transition. As the twist angle is lowered, they become more robust
and can persist beyond the lowest flux ratio of $1/12$ studied in this work.

Figure 2: Compilation of B-SCHF results at two representative twist angles,
$\theta=1.32^{\circ}$ (a-f) and $\theta=1.2^{\circ}$ (g-l). The non-
interacting spectra are shown in (a) and (g). The renormalized spectra along
$(s,t)=(-3,-1)$ are shown in (b) and (h). The occupied electronic states are
colored in red and green at $1.32^{\circ}$, signifying the first order phase
transition between CHFs (red) and nearly compressible states (green). The
occupied electronic states are colored in red and blue at $1.2^{\circ}$,
signifying the first order phase transition between CHFs (red) and gapped
states with intervalley coherence (blue). The alternating colors near
$\phi/\phi_{0}=0.2$ are due to near energetic degeneracies ($\sim 0.01$meV)
between competing states, and are not resolved within our calculation. (c) and
(i) show the non-monotonic behaviors of $\Delta(\phi/\phi_{0})$ along
$(s,t)=(0,-4)$, $(-1,-3)$, $(-2,-2)$, and $(-3,-1)$ Density matrix analysis.
(d) and (j) show the absolute value of the non-vanishing matrix elements of
the density matrices at $\phi/\phi_{0}=2/5$ at $\mathbf{k}=\mathbf{0}$. $(e)$
and $(j)$ are the absolute values of the representative density matrices at
$\phi/\phi_{0}$ for a nearly compressible state (e) and an intervalley
coherent state (k). (f) and (l) are replots of
$\Delta(n/n_{s},\phi/\phi_{0})$, colored in the same manner as (b) and (h)
respectively.

To better understand their nature, we compile detailed results for two
representative twist angles $1.38^{\circ}$ [Figs. 2(a) to (f)] and
$1.20^{\circ}$ [Figs. 2(g) to (l)]. Fig. 2(a) shows the non-interacting
spectra of valley $\mathbf{K}$ and for one spin component (neglecting Zeeman
effect) at $\theta=1.32^{\circ}$. The magnetic subbands marked in red denote
the Chern -1 group below the charge neutrality point. The analogous group of
subbands above the charge neutrality point is related to it by particle hole
symmetry, and also carries total Chern number $-1$. The remaining two subbands
emanate from the zeroth Landau levels (zLLs) of the energetically split Dirac
points, and carry Chern number $+1$ each, such that the total Chern number of
all magnetic subbands is zero. Fig. 2(b) shows the single-particle spectra
including Coulomb interactions along $(s,t)=(-3,-1)$. At $\phi/\phi_{0}\gtrsim
1/3$, the spectrum is gapped and the occupied states marked in red have a
large overlap onto the states marked in red in Fig. 2(a). To quantify the
overlap we make use of the Hartree-Fock density matrix, defined as
$\hat{Q}_{\eta sr,\eta^{\prime}s^{\prime}r^{\prime}}(\mathbf{k})\equiv\langle
d^{\dagger}_{\eta
sr,\mathbf{k}}d_{\eta^{\prime}s^{\prime}r^{\prime},\mathbf{k}}\rangle$; for a
state with unbroken valley $U(1)_{v}$ and spin $U(1)_{s}$ symmetries, it has
the spin-valley diagonal form $\mathcal{Q}_{r,r^{\prime}}^{(\eta
s)}(\mathbf{k})\delta_{\eta,\eta^{\prime}}\delta_{s,s^{\prime}}$. At high $B$,
the $(-3,-1)$ CCI is valley and spin polarized. A representative
$|\mathcal{Q}^{\mathbf{K}\uparrow}_{r,r^{\prime}}(\mathbf{0})|$ is shown in
Fig. 2(d). It is predominantly diagonal in the magnetic subband index, mostly
occupying the lower $q-p$ magnetic subbands, i.e., the group states with total
Chern -1 marked red in Fig. 2(a). At lower $B$, the gapped spectrum (red) in
Fig. 2(b) transitions into a nearly gapless spectrum (green) via a first order
phase transition. These low $B$ states do not break the valley $U(1)_{v}$ and
spin $U(1)_{s}$ symmetries, and their illustrative density matrix is shown in
the Fig. 2(e). For $(-2,-2)$ and $(-1,-3)$ the density matrices of the CCIs at
higher $B$ are also shown in Fig. 2(d). These CCIs are also valley and spin
polarized, similarly mostly populating the lower Chern $-1$ group of magnetic
subbands for the specified valley and spin. $|\mathcal{Q}^{(\eta
s)}(\mathbf{0})|$ are identical for all occupied flavors. As shown in Fig.
2(f), CCIs along $(-3,-1)$, $(-2,-2)$ and $(-1,-3)$ display quantitative
differences in the critical $B$.

Although these CCIs are closely related to the HSFs discussed in Ref. [15, 7,
22], the band structure renormalization is apparent in the non-vanishing off-
diagonal matrix elements of $\mathcal{Q}^{(\eta
s)}_{r,r^{\prime}}(\mathbf{k})$, signifying hybridization with the higher
energy subbands (marked by grey in Fig. 2(a)). For this reason we refer to
them as CHFs. As further demonstrated in Fig. S4, the density matrices of the
CHFs assume a much simpler structure when expressed in the eigenbasis of the
valley and spin symmetric $(0,-4)$ Chern insulating state, which limits to an
interaction-renormalized semimetal at $B=0$ (see Fig. S6) [36]. Finally, in
Fig. 2(c) we show the evolution of the single particle gap as a function of
$B$ along $(0,-4)$, $(-1,-3)$, $(-2,-2)$, and $(-3,-1)$. Generically, the gap
shows a broad plateau and non-monotonic dependence on $B$.

At the lower twist angle of $\theta=1.2^{\circ}$ the CHFs at larger $B$ are
very similar to those at higher twist angles, however with stronger
hybridization into the higher energy magnetic subbands, as seen in Fig. 2(j).
In contrast to $1.32^{\circ}$, the CHFs lose via a first order phase
transition to gapped states with intervalley coherence at lower $B$, as is
apparent in the non-zero matrix elements
$|\hat{Q}_{\mathbf{K}sr,\mathbf{K}^{\prime}sr^{\prime}}(\mathbf{0})|$
illustrated in Fig. 2(k). The intervalley coherent state may also break the
magnetic translation symmetry $\hat{t}_{\mathbf{L}_{2}}$. However, our B-SCHF
calculations cannot determine a unique wavevector $\mathbf{q}$ associated with
the translation symmetry breaking, suggesting that there is a manifold of
$\mathbf{q}\neq 0$ intervalley coherent states with similar energies [36].
Interestingly, such intervalley coherent states hybridize subbands from the
lower half of the (renormalized) Hofstadter spectra in each valley (see Fig.
S5), similar to the IKS ground states at $B=0$ discussed in Ref. [36].

Figure 3: Averaged $\sigma_{z}\tau_{z}$ per occupied single electron states at
$(s,t)=(0,-4)$, $(-1,-3)$, $(-2,-2)$ and $(-3,-1)$ for (a) $\phi/\phi_{0}=2/5$
with heterostrain, (b) $\phi/\phi_{0}=1/8$ with heterostrain, (c)
$\phi/\phi_{0}=2/5$ without heterostrain, and (d) $\phi/\phi_{0}=1/8$ without
heterostrain. The dotted line in each panel is the value for sCI, and the
dashed line corresponds to HSF.

Crossover to strong coupling regime: In order to clarify the connection
between the CHFs and the sCIs, we first note that the sCIs saturate the
expectation value of $\sigma_{z}\tau_{z}$ for the occupied electronic states,
where $\sigma_{z}$ and $\tau_{z}$ are Pauli matrices acting in the sublattice
and valley subspace respectively [27, 28, 41, 47]. The solid lines in Fig.
3(a) and (b) show calculated $\langle\sigma_{z}\tau_{z}\rangle$ per electron
as a function of twist angle in the presence of heterostrain along $(-3,-1)$,
$(-2,-2)$, $(-1,-3)$ and $(0,-4)$, for $\phi/\phi_{0}=2/5$ and $1/8$
respectively. For comparison, the upper dashed line corresponds to the sCI
limit, and the lower dashed line to the HSF limit. This measure shows that the
CHFs are clearly distinct from both the HSFs and the sCIs in the presence of
heterostrain. Qualitatively, the Coulomb repulsion favors an increased
$\sigma_{z}\tau_{z}$. The abrupt jump as angle is tuned (e.g., $\sim
1.5^{\circ}$ for $(-3,-1)$ in Fig. 3(b) ) marks the first order transition
between the aforementioned nearly compressible state and the CHF.

On the other hand, in the absence of heterostrain, as shown in Fig. 3(c) and
(d), the CHFs smoothly crosses over into the sCIs upon lowering the twist
angle. For our model parameters, there is a collapse of the
$\langle\sigma_{z}\tau_{z}\rangle$ at lower twist angles, when the CHFs become
energetically less favorable than populating the Landau quantized excitation
spectra of the IVC states [47, 38, 31]. As further demonstrated in Fig. S3,
this transition depends sensitively on model parameters, and can be pushed
toward lower $\phi/\phi_{0}$ (e.g., by moving toward the chiral limit, see SI
and Ref. [48]). Refs. [20, 11] report that the $(2,2)$ persist down to
$\phi/\phi_{0}\sim 1/100$, and therefore argue that they correspond to the
zCIs (more precisely, sCIs). Our quantitative calculations clearly demonstrate
that such states can indeed be stablized at weak fields by small changes of
the model parameters.

Figure 4: Occupation number $n_{0}(\mathbf{k})$ of the lower zLL of the non-
interacting spectra in valley $\mathbf{K}$ and for spin $\downarrow$. Results
are obtained at $1.05^{\circ}$ for correlated Chern inslating states with
fractional $s$, marked by the orange dashed lines in Fig. 1(f).

Additional CCIs: In addition to the aforementioned CCIs, we also find
additional correlated insulating states in the phase diagram with heterostrain
(Fig. 1). The most prominent states emanate from the charge neutrality point
($n/n_{s}=0$), and are identified as QHFMs within the zLLs of the
energetically split Dirac cones [49, 50, 51]. QHFMs emanating from the band
bottom ($n/n_{s}=-4$) are distinctive at $1.38^{\circ}$ but become
progressively weaker as the angle decreases. Additionally, at higher
$\phi/\phi_{0}$ a gapped state with $(s,t)=(-2,0)$ is observed in Figs.
1(a-c). We identify it as a Quantum Spin Hall (QSH) insulator due to strong
spin splitting near $\phi/\phi_{0}=1/2$, as demonstrated in the non-
interacting Hofstadter spectra in Fig. S1.

Closer to the magic angle, we also find unusual sequences with
$(s,t)=(-3,0),(-2,-1),(-1,-2)$. These states have the same $s$ as the sequence
$(s,t)=(-3,-1),(-2,-2),(-1,-3)$ discussed previously, but with a Chern number
that is higher by $+1$. They all have intervalley coherence, with structure of
the density matrices resembling that of an IKS ground state. However, careful
analysis shows that they do not correspond to populating the Landau quantized
excitation spectra of the CCIs along $(s,t)=(-3,-1),(-2,-2),(-1,-3)$. Rather,
they are distinct intervalley coherent states with a regorganization of the
intervalley coherence between the magnetic subbands (see Fig. S5).

At $1.05^{\circ}$, we also find CCIs with fractional $s$ along
$(-{2}/{3},-3)$, $(-{1}/{2},-3)$, and $(-{3}/{2},-2)$, see Fig. 1(f). These
states break magnetic translation symmetry. We identify them as striped states
with period $m$ along the $\mathbf{L}_{2}$ direction, such that the density
matrix is invariant under $\hat{t}_{\mathbf{L}_{2}}^{m}$. We use the electron
occupation number of the lower zLL (see Fig. 2(a) and (i) for example) in
valley $\mathbf{K}$ and for spin $\downarrow$ to illustrate the striped
states. We define it as $n_{0}(\mathbf{k})$, and show its momentum dependence
in the magnetic Brillouin zone at $\phi/\phi_{0}=1/6$ for $(-{2}/{3},-3)$,
$(-{1}/{2},-3)$, and $(-{3}/{2},-2)$ in Figs. 4 (a), (b) and (c) respectively.
At $(-1/2,-3)$ and $(-3/2,2)$, the fractional part of $s$ corresponds to half-
filling of a valley and spin flavor, and we identify the period of the striped
state as $m=2$. At $(-2/3,-3)$, the fractional part of $s$ corresponds to two
thirds filling of a flavor, and we identify the stripe period as $m=3$.

Conclusion: In summary, by performing a comprehensive self-consistent Hartree-
Fock study of the continuum Bistritzer-MacDonald model in finite magnetic
fields, we unravel the nature of the prominent correlated Chern insulators
observed in a wide range of TBG experiments. For realistic heterostrain, these
CCIs are stablized at higher magnetic fields, and correspond to valley and
spin polarizations of the interaction-renormalized magnetic subbands (CHFs).
Upon lowering magnetic field, the CHFs become energetically less favored,
losing to nearly compressible states at larger twist angles and other gapped
states with intervalley coherence at smaller twist angles. In absence of
heterostrain and at higher magnetic fields, the CHF crosses over to the strong
coupling Chern insulating states (sCIs) as the twist angle decreases. At lower
fields, competing states with intervalley coherence become more energetically
favored, and the transition is marked by a collapse of the averaged sublattice
polarization per occupied single-electron state.

Our calculations also predict additional gapped correlated insulating states
beyond the $(s,t)=(0,\pm 4),\pm(1,3),\pm(2,2),\pm(3,1)$ sequence, notably the
striped states at fractional $s$. Given that our calculations have direct
access to the interaction renormalized single-electron excitation spectra at a
given filling and magnetic field (see Fig. S7), comparisons with experiments
such as STM can be made to facilitate the characterization of the panoply of
correlated insulating states.

Acknowledgements X.W. and O. V. acknowledge invaluable discussions with B.
Andrei Bernevig, J. Finney, M.-H. He and J. Kang. X.W. acknowledges financial
support from the Gordon and Betty Moore Foundation’s EPiQS Initiative Grant
GBMF11070. O.V. was supported by NSF Grant No. DMR-1916958 and is partially
funded by the Gordon and Betty Moore Foundation’s EPiQS Initiative Grant
GBMF11070, National High Magnetic Field Laboratory through NSF Grant No.
DMR-1157490 and the State of Florida. Most of the computing for this project
was performed on the HPC at the Research Computing Center at the Florida State
University (FSU).

## References

  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Graphene bilayers with a twist, Nature Materials 19, 1265 (2020).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, Nature Physics 16, 725 (2020).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * 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).
  * Tomarken  _et al._ [2019] S. L. Tomarken, Y. Cao, A. Demir, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, and R. C. Ashoori, Electronic compressibility of magic-angle graphene superlattices, Phys. Rev. Lett. 123, 046601 (2019).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Flavour hund’s coupling, chern gaps and charge diffusivity in moiré graphene, Nature 592, 43 (2021).
  * 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).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Strongly correlated chern insulators in magic-angle twisted bilayer graphene, Nature 588, 610 (2020).
  * Nuckolls  _et al._ [2023] K. P. Nuckolls, R. L. Lee, M. Oh, D. Wong, T. Soejima, J. P. Hong, D. Călugăru, J. Herzog-Arbeitman, B. A. Bernevig, K. Watanabe, T. Taniguchi, N. Regnault, M. P. Zaletel, and A. Yazdani, Quantum textures of the many-body wavefunctions in magic-angle graphene, Nature 620, 525 (2023).
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Unconventional sequence of correlated chern insulators in magic-angle twisted bilayer graphene, Nature Physics 17, 1210 (2021).
  * Xie  _et al._ [2021] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Fractional chern insulators in magic-angle twisted bilayer graphene, Nature 600, 439 (2021).
  * Saito  _et al._ [2021a] Y. Saito, J. Ge, L. Rademaker, K. Watanabe, T. Taniguchi, D. A. Abanin, and A. F. Young, Hofstadter subband ferromagnetism and symmetry-broken chern insulators in twisted bilayer graphene, Nature Physics 17, 478 (2021a).
  * Yu  _et al._ [2022] J. Yu, B. A. Foutty, Z. Han, M. E. Barber, Y. Schattner, K. Watanabe, T. Taniguchi, P. Phillips, Z.-X. Shen, S. A. Kivelson, and B. E. Feldman, Correlated hofstadter spectrum and flavour phase diagram in magic-angle twisted bilayer graphene, Nature Physics 18, 825 (2022).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Das  _et al._ [2021] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Symmetry-broken chern insulators and rashba-like landau-level crossings in magic-angle bilayer graphene, Nature Physics 17, 710 (2021).
  * Serlin  _et al._ [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Intrinsic quantized anomalous hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * Stepanov  _et al._ [2021] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Competing zero-field chern insulators in superconducting twisted bilayer graphene, Phys. Rev. Lett. 127, 197701 (2021).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Chern insulators, van hove singularities and topological flat bands in magic-angle twisted bilayer graphene, Nature Materials 20, 488 (2021).
  * Choi  _et al._ [2021] Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Correlation-driven topological phases in magic-angle twisted bilayer graphene, Nature 589, 536 (2021).
  * Grover  _et al._ [2022] S. Grover, M. Bocarsly, A. Uri, P. Stepanov, G. Di Battista, I. Roy, J. Xiao, A. Y. Meltzer, Y. Myasoedov, K. Pareek, K. Watanabe, T. Taniguchi, B. Yan, A. Stern, E. Berg, D. K. Efetov, and E. Zeldov, Chern mosaic and berry-curvature magnetism in magic-angle graphene, Nature Physics 18, 885 (2022).
  * Uri  _et al._ [2020] A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani, D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and E. Zeldov, Mapping the twist-angle disorder and landau levels in magic-angle graphene, Nature 581, 47 (2020).
  * 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).
  * 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).
  * 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).
  * 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).
  * Zhang  _et al._ [2019] Y.-H. Zhang, D. Mao, and T. Senthil, Twisted bilayer graphene aligned with hexagonal boron nitride: Anomalous hall effect and a lattice model, Phys. Rev. Res. 1, 033126 (2019).
  * 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).
  * 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 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).
  * 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).
  * 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).
  * Bi  _et al._ [2019] Z. Bi, N. F. Q. Yuan, and L. Fu, Designing flat bands by strain, Phys. Rev. B 100, 035448 (2019).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Kekulé spiral order at all nonzero integer fillings in twisted bilayer graphene, Phys. Rev. X 11, 041063 (2021).
  * Wang  _et al._ [2023] X. Wang, J. Finney, A. L. Sharpe, L. K. Rodenbach, C. L. Hsueh, K. Watanabe, T. Taniguchi, M. A. Kastner, O. Vafek, and D. Goldhaber-Gordon, Unusual magnetotransport in twisted bilayer graphene from strain-induced open fermi surfaces, Proceedings of the National Academy of Sciences 120, e2307151120 (2023).
  * 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).
  * 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. Res. 3, 013033 (2021).
  * Xie  _et al._ [2023] F. Xie, J. Kang, B. A. Bernevig, O. Vafek, and N. Regnault, Phase diagram of twisted bilayer graphene at filling factor $\nu=\pm{}3$, Phys. Rev. B 107, 075156 (2023).
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. iii. interacting hamiltonian and exact symmetries, Phys. Rev. B 103, 205413 (2021).
  * 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, E. Berg, P. Jarillo-Herrero, and S. Ilani, Cascade of phase transitions and dirac revivals in magic-angle graphene, Nature 582, 203 (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 (2020).
  * Saito  _et al._ [2021b] 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 (2021b).
  * 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, Phys. Rev. B 99, 035111 (2019).
  * Herzog-Arbeitman  _et al._ [2022] J. Herzog-Arbeitman, A. Chew, D. K. Efetov, and B. A. Bernevig, Reentrant correlated insulators in twisted bilayer graphene at 25 t ($2\pi$ flux), Phys. Rev. Lett. 129, 076401 (2022).
  * Wang and Vafek [2022] X. Wang and O. Vafek, Narrow bands in magnetic field and strong-coupling hofstadter spectra, Phys. Rev. B 106, L121111 (2022).
  * 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).
  * Nomura and MacDonald [2006] K. Nomura and A. H. MacDonald, Quantum hall ferromagnetism in graphene, Phys. Rev. Lett. 96, 256602 (2006).
  * Kharitonov [2012] M. Kharitonov, Phase diagram for the $\nu=0$ quantum hall state in monolayer graphene, Phys. Rev. B 85, 155439 (2012).
  * Young  _et al._ [2012] A. F. Young, C. R. Dean, L. Wang, H. Ren, P. Cadden-Zimansky, K. Watanabe, T. Taniguchi, J. Hone, K. L. Shepard, and P. Kim, Spin and valley quantum hall ferromagnetism in graphene, Nature Physics 8, 550 (2012).

Supplemental Material for “Theory of correlated Chern insulators in twisted
bilayer graphene”

###### Contents

  1. I TBG with heterostrain at zero magnetic field
  2. II TBG in magnetic field
    1. II.1 Landau level wavefunctions of monolayer graphene
    2. II.2 Magnetic translation group eigenstates
    3. II.3 Matrix elements for a generic operator $A_{\mathbf{q}}\equiv\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}$
    4. II.4 Matrix elements of strained BM Hamiltonian
      1. II.4.1 Interlayer terms
      2. II.4.2 Narrow band eigenstates
      3. II.4.3 Degeneracy of magnetic subbands generated by $\hat{t}_{\mathbf{L}_{2}}$
    5. II.5 Matrix elements of $A_{\mathbf{q}}=\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}$ with respect to narrow band eigenstates
  3. III Self-consistent Hartree-Fock method in finite magnetic field
    1. III.1 Product state and one-particle density matrix
    2. III.2 Hartree Fock Hamiltonian
  4. IV Extended Figures

##  I TBG with heterostrain at zero magnetic field

We consider the limit of small twist angle $\theta$ and small deformations of
moiré superlattice due to uniaxial heterostrain, such that the moiré
reciprocal lattice vectors are given by:
$\mathbf{g}_{i=1,2}=\mathcal{E}^{T}\mathbf{G}_{i=1,2}$ [35], where
$\mathbf{G}_{i}$ are reciprocal lattice vectors of the untwisted and
undeformed monolayer graphene, and

$\mathcal{E}=\mathcal{T}(\theta)+\mathcal{S}(\epsilon,\varphi),$ (S1)

where

$\mathcal{T}(\theta)=\begin{pmatrix}0&-\theta\\\ \theta&0\end{pmatrix},\
\mathcal{S}(\epsilon,\varphi)=R_{\varphi}^{T}\begin{pmatrix}-\epsilon&0\\\
0&\nu\epsilon\end{pmatrix}R_{\varphi}.$ (S2)

Here $\\{\varepsilon,\varphi\\}$ parameterize the uniaxial heterostrain
strength and orientation, $R_{\varphi}$ is the two-dimensional rotation
matrix, and $\nu\approx 0.16$ is the Poisson ratio. The moiré lattice vectors
$\mathbf{L}_{i=1,2}$ are uniquely defined through the relation
$\mathbf{L}_{i}\cdot\mathbf{g}_{j}=2\pi\delta_{ij}$. Specifically,

$\mathbf{L}_{1}=\frac{2\pi}{(\mathbf{g}_{1}\times\mathbf{g}_{2})\cdot\hat{z}}(\mathbf{g}_{2}\times\hat{z}),\
\mathbf{L}_{2}=\frac{2\pi}{(\mathbf{g}_{1}\times\mathbf{g}_{2})\cdot\hat{z}}(\hat{z}\times\mathbf{g}_{1}).$
(S3)

We consider the strained Bistritzer-MacDonald (BM) Hamiltonian discussed in
our recent work [37], with the continuum Hamiltonian given as:

$H_{\eta}=(\sum_{l=t,b}H_{\eta,l}^{intra})+H_{\eta}^{inter},$ (S4)

where $\eta=\pm 1$ labels $\mathbf{K}\ (\mathbf{K}^{\prime})$ valleys of
monolayer graphene, $l=t,b$ labels the top (bottom) of the two graphene
layers. The interlayer Hamiltonian is given by:

$H_{\eta}^{inter}\approx\int\mathrm{d}^{2}\mathbf{r}\psi^{\dagger}_{\eta,t}(\mathbf{r})\left(\sum_{j=1,2,3}T_{\eta,j}e^{-i\eta\mathbf{q}_{j}\cdot\mathbf{r}}\right)\psi_{\eta,b}(\mathbf{r})+h.c.,$
(S5)

where
$\psi_{\eta,l}(\mathbf{r})\equiv(\psi_{\eta,l,A}(\mathbf{r}),\psi_{\eta,l,B}(\mathbf{r}))^{T}$
is a spinor in the sublattice basis for a given valley and layer. We have
suppressed the spin index for simplicity. $\mathbf{q}_{j=1,2,3}$ are the three
nearest neighbor bonds of the reciprocal honeycomb lattice, and

$T_{\eta,j}=w_{0}\sigma_{0}+w_{1}\left(\cos\frac{2\pi(j-1)}{3}\sigma_{x}+\eta\sin\frac{2\pi(j-1)}{3}\sigma_{y}\right).$
(S6)

Here $(\sigma_{0},\sigma_{x},\sigma_{y})$ are Pauli matrices acting on
sublattice degrees of freedom. $w_{0}$ and $w_{1}$ are intra-sublattice and
inter-sublattice tunneling strengths between the two layers. The “chiral
limit” discussed in Ref. [48] corresponds to $w_{0}/w_{1}=0$.

The intra-layer Hamiltonian is given by:

$\begin{split}&H_{\eta,l}^{intra}=\alpha\int\mathrm{d}^{2}\mathbf{r}\psi^{\dagger}_{\eta,l}(\mathbf{r})(\tr[\mathcal{E}_{l}]\sigma_{0})\psi_{\eta,l}(\mathbf{r})\\\
&+\frac{\hbar
v_{F}}{a}\int\mathrm{d}^{2}\mathbf{r}\psi^{\dagger}_{\eta,l}(\mathbf{r})\left[(-i\nabla-\mathbf{A}_{\eta,l})\cdot(\eta\sigma_{x},\sigma_{y})\right]\psi_{\eta,l}(\mathbf{r}),\end{split}$
(S7)

where the first term is the deformation potential that couples to the electron
density. $\mathbf{A}_{\eta,l}$ is the pseudovector potential that comes from
changes in the inter-sublattice hopping due to deformations, and changes sign
between graphene valleys. It is given as
$\mathbf{A}_{\eta,l}=\frac{\sqrt{3}\beta}{2a}\eta(\epsilon_{l,xx}-\epsilon_{l,yy},-2\epsilon_{l,xy})$.

We use $\alpha\approx-4.1\ \mathrm{eV}$ and $\beta\approx 3.14$. We further
fix $\hbar v_{F}/a=2.68\mathrm{eV}$, $w_{1}=110\mathrm{meV}$ and
$w_{0}/w_{1}=0.7$ in our calculations, and also set constants $\hbar=a=1$ in
the remainder of this document. Other values of $w_{0}/w_{1}$ are also used
and discussed in the manuscript as a way to stablize the strong-coupling Chern
insulating states.

##  II TBG in magnetic field

###  II.1 Landau level wavefunctions of monolayer graphene

We begin with a brief discussion of the Landau level (LL) eigenstates of the
Dirac Hamiltonian of monolayer graphene. In a finite magnetic field $B$, the
Dirac Hamiltonian in valleys $\mathbf{K}$ and $\mathbf{K}^{\prime}$ are given
by (here for negative charged electrons, the minimal coupling is
$\mathbf{p}+e\mathbf{A}$ where $e$ is positive):

$\displaystyle\hat{H}^{\mathbf{K}}_{l}(\mathbf{p}+e\mathbf{A})$
$\displaystyle=v_{F}\vec{\sigma}\cdot(\mathbf{p}-\mathbf{K}_{l}+e\mathbf{A}),$
(S8) $\displaystyle\hat{H}_{l}^{\mathbf{K}^{\prime}}(\mathbf{p}+e\mathbf{A})$
$\displaystyle=-v_{F}\vec{\sigma}^{*}\cdot(\mathbf{p}-\mathbf{K}_{l}^{\prime}+e\mathbf{A}).$
(S9)

Here $l=1,2$ is the layer index, $\mathbf{K}_{l}=(K_{l,x},K_{l,y})$ is the
position of the Dirac cone in the reciprocal space.

We choose the Landau gauge such that the vector potential $\mathbf{A}\equiv
Bx\hat{y}$, where the global $x-y$ coordinate system is defined such that
$\mathbf{L}_{2}$ is along the $+\hat{y}$ direction. When TBG is subject to
heterstrain, this amounts to a small angle rotation $\theta_{0}$ of the global
coordinate system with respect to the case in absence of heterostrain. The
eigenstates of the Dirac Hamiltonian are solved by going to the harmonic
oscillator basis: $x=\frac{\ell}{\sqrt{2}}(a+a^{\dagger})$, and
$p_{x}=\frac{1}{i\sqrt{2}\ell}(a-a^{\dagger})$, where $\ell\equiv 1/\sqrt{eB}$
is the magnetic length.

In valley $\mathbf{K}$, the particle-hole symmetric LL eigenstates are given
as:

$\bra{\mathbf{r}}\ket{\psi_{n\gamma
l}(k_{2})}=e^{iK_{l,x}x}e^{i\frac{2\pi}{L_{2}}k_{2}y}\frac{1}{\sqrt{2}}\begin{pmatrix}-i\gamma
e^{-i\theta_{0}}\phi_{n-1}(x+\tilde{k}_{2,l}\ell^{2})\\\
\phi_{n}(x+\tilde{k}_{2,l}\ell^{2})\end{pmatrix},$ (S10)

where $\epsilon_{n\gamma}=\frac{v_{F}}{\ell}\gamma\sqrt{2n}$ is the energy of
the Dirac Hamiltonian, labeled by $n=1,2,\dots$, and $\gamma=\pm 1$
corresponds to positive and negative energy solutions. We have defined
$L_{2}\equiv|\mathbf{L}_{2}|$. In addition, there is an anomalous zero energy
state given by

$\bra{\mathbf{r}}\ket{\psi_{0l}(k_{2})}=e^{iK_{l,x}x}e^{i\frac{2\pi}{L_{2}}k_{2}y}\begin{pmatrix}0\\\
\phi_{0}(x+\tilde{k}_{2,l}\ell^{2})\end{pmatrix}$ (S11)

which lives on the B sublattice. $\phi_{n}(x)$ is the eigenfunction of
$a^{\dagger}a$, and is given to be:

$\phi_{n}(x)=\frac{1}{\pi^{1/4}}\frac{1}{\sqrt{2^{n}n!}}e^{-x^{2}/2\ell^{2}}H_{n}(x/\ell),$
(S12)

where $H_{n}(x)$ is the Hermite polynomial. The shift in the position for a
given wavevector $k_{2}\mathbf{g}_{2}$ is given by:

$\tilde{k}_{2,l}\ell^{2}=\left(\frac{2\pi}{L_{2}}k_{2}-K_{l,y}\right)\ell^{2}.$
(S13)

In valley $\mathbf{K}^{\prime}$, the Landau level wavefunctions are solved in
a similar manner, and the LL eigenbasis is given by:

$\bra{\mathbf{r}}\ket{\psi_{n\gamma
l}(k_{2})}=e^{iK_{l,x}^{\prime}x}e^{i\frac{2\pi}{L_{2}}k_{2}y}\frac{1}{\sqrt{2}}\begin{pmatrix}\phi_{n}(x+\tilde{k}_{2,l}\ell^{2})\\\
i\gamma e^{-i\theta_{0}}\phi_{n-1}(x+\tilde{k}_{2,l}\ell^{2})\end{pmatrix}$
(S14)

with eigenenergy $\epsilon_{n\gamma}=v_{F}\gamma\frac{\sqrt{2n}}{\ell}$, and
eigenstate for the zeroth LL:

$\bra{\mathbf{r}}\ket{\psi_{0l}(k_{2})}=e^{iK_{l,x}^{\prime}x}e^{i\frac{2\pi}{L_{2}}k_{2}y}\begin{pmatrix}\phi_{0}(x+\tilde{k}_{2,l}\ell^{2})\\\
0\end{pmatrix}.$ (S15)

For notational convenience we define a ket state $\ket{\Phi^{(\eta)}_{n\gamma
l}(k_{2})}$ such that:

$\bra{\mathbf{r}}\ket{\psi^{(\eta)}_{n\gamma l}(k_{2})}=e^{iK_{\eta
l,x}x}e^{i\frac{2\pi}{L_{2}}k_{2}y}\bra{\mathbf{r}}\ket{\Phi^{(\eta)}_{n\gamma
l}(k_{2})}.$ (S16)

The ket state has the following real space properties:

$\bra{\mathbf{r}}\ket{\Phi^{(\eta)}_{n\gamma
l}(k_{2})}\equiv\Phi^{(\eta)}_{n\gamma l}(x+\tilde{k}_{2,l}\ell^{2}),$ (S17)

where the definitions of $\Phi^{(\eta)}_{n\gamma
l}(x+\tilde{k}_{2,l}\ell^{2})$ are inferred from Eqs. (S10), (S11), (S14) and
(S15).

###  II.2 Magnetic translation group eigenstates

At rational magnetic flux ratios $\phi/\phi_{0}=p/q$ where $\phi_{0}=h/e$ is
the magnetic flux quantum, $p$ and $q$ are coprime integers, the strained BM
Hamiltonian is invariant under magnetic translations,
$\hat{t}_{\mathbf{L}_{1}}$ and $\hat{t}_{\mathbf{L}_{2}}$, such that:

$\hat{t}_{\mathbf{L}_{1}}=e^{-i\mathbf{q}_{\phi}\cdot(\mathbf{r}-\mathbf{L}_{1}/2)}\hat{T}_{\mathbf{L}_{1}},\
\hat{t}_{\mathbf{L}_{2}}=\hat{T}_{\mathbf{L}_{2}},\
\mathbf{q}_{\phi}=\frac{2\pi}{|\mathbf{L}_{2}|}\frac{\phi}{\phi_{0}}\hat{y},\
\hat{t}_{\mathbf{L}_{2}}\hat{t}_{\mathbf{L}_{1}}=e^{i2\pi\phi/\phi_{0}}\hat{t}_{\mathbf{L}_{1}}\hat{t}_{\mathbf{L}_{2}}.$
(S18)

Here
$\hat{T}_{\mathbf{L}_{i}}\psi(\mathbf{r})=\psi(\mathbf{r}-\mathbf{L}_{i})$. As
a result, the eigenstates of the strained Bistritzer-MacDonald Hamiltonian are
simultaneous eigenstates of the magnetic translation group (MTG). The MTG
basis states can be generated from the LL basis states discussed in the
previous section:

$\begin{split}\ket{V^{(\eta)}_{n\gamma
l}(\mathbf{k})}&=\frac{1}{\sqrt{\mathcal{N}}}\sum_{s_{1}=-\infty}^{\infty}e^{i2\pi
k_{1}s_{1}}\hat{t}_{\mathbf{L}_{1}}^{s_{1}}\ket{\psi^{(\eta)}_{n\gamma
l}(k_{2})}\\\
&=\frac{1}{\sqrt{\mathcal{N}}}\sum_{s_{1}}e^{i2\pi(k_{1}-k_{2}\frac{L_{1y}}{L_{2}})s_{1}}e^{i\frac{1}{2}s_{1}^{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-is_{1}K_{\eta
l,x}L_{1,x}}\ket{\psi^{(\eta)}_{n\gamma
l}(k_{2}-s_{1}\frac{\phi}{\phi_{0}})},\end{split}$ (S19)

where $\mathcal{N}$ is a normalization factor, and $\mathbf{k}\equiv
k_{1}\mathbf{g}_{1}+k_{2}\mathbf{g}_{2}$. It is straightforward to check that:

$\displaystyle\ket{V_{n\gamma l}^{(\eta)}(k_{1}+1,k_{2})}$
$\displaystyle=\ket{V_{n\gamma l}^{(\eta)}(k_{1},k_{2})},$ (S20)
$\displaystyle\ket{V_{n\gamma l}^{(\eta)}(k_{1},k_{2}+\frac{\phi}{\phi_{0}})}$
$\displaystyle=e^{iK_{l,x}L_{1,x}}e^{-i2\pi(k_{1}-k_{2}\frac{L_{1y}}{L_{2}})}\ket{V_{n\gamma
l}^{(\eta)}(k_{1},k_{2})},$ (S21)
$\displaystyle\hat{t}_{\mathbf{L}_{2}}\ket{V_{n\gamma
l}^{(\eta)}(k_{1},k_{2})}$ $\displaystyle=e^{-i2\pi k_{2}}\ket{V_{n\gamma
l}^{(\eta)}(k_{1}+\frac{\phi}{\phi_{0}},k_{2})},$ (S22)
$\displaystyle\bra{V_{n\gamma
l}^{(\eta)}(k_{1},k_{2})}\ket{V_{n^{\prime}\gamma^{\prime}l^{\prime}}^{(\eta^{\prime})}(p_{1},p_{2})}$
$\displaystyle=\delta_{\eta,\eta^{\prime}}\delta_{l,l^{\prime}}\delta_{k_{1},k_{2}}\delta_{p_{1},p_{2}}\delta_{n,n^{\prime}}\delta_{\gamma,\gamma^{\prime}}.$
(S23)

Therefore, the MTG eigenstates defined in
$(k_{1},k_{2})\in[0,1)\times[0,\frac{\phi}{\phi_{0}})$ form a complete and
orthornomal basis set in a finite magnetic field. This is referred to as the
magnetic Brillouin zone. In order to find the eigenstates of the narrow bands
for strained BM Hamiltonian in $B$, we also need to compute the matrix
elements of the interlayer tunneling terms in the above MTG basis states. From
now on we work with rational flux ratio $\phi/\phi_{0}=p/q$, where $p$ and $q$
are coprime.

###  II.3 Matrix elements for a generic operator
$A_{\mathbf{q}}\equiv\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}$

We first work out the matrix elements for a generic operator in the MTG
eigenstates, $A_{\mathbf{q}}=\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}$, where
$\hat{A}$ is a space-independent matrix in the layer, sublattice, and valley
basis. For example, the Fourier transform of the electron density operator,
$\hat{\rho}_{\mathbf{q}}\equiv\hat{\rho}e^{-i\mathbf{q}\cdot\mathbf{r}}$, has
$\hat{\rho}=\hat{I}$, where $\hat{I}$ is the identity matrix. We denote
$\mathbf{q}=q_{1}\mathbf{g}_{1}+q_{2}\mathbf{g}_{2}$ and
$\mathbf{p}=p_{1}\mathbf{g}_{1}+p_{2}\mathbf{g}_{2}$. We therefore have:

$\begin{split}&\bra{V^{(\eta)}_{n\gamma
l}(\mathbf{k})}A_{\mathbf{q}}\ket{V^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(\mathbf{p})}\\\
=&\frac{1}{\mathcal{N}}\sum_{s_{1}s_{1}^{\prime}}e^{-i2\pi k_{1}s_{1}}e^{i2\pi
p_{1}s_{1}^{\prime}}\bra{\psi_{n\gamma
l}^{(\eta)}(k_{2})}\hat{t}_{\mathbf{L}_{1}}^{-s_{1}}\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}\hat{t}_{\mathbf{L}_{1}}^{s_{1}^{\prime}}\ket{\psi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(p_{2})}\\\
=&\frac{1}{\mathcal{N}}\sum_{s_{1}s_{1}^{\prime}}e^{-i2\pi(k_{1}+q_{1})s_{1}}e^{i2\pi
p_{1}s_{1}^{\prime}}\bra{\psi_{n\gamma
l}^{(\eta)}(k_{2})}\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}\hat{t}_{\mathbf{L}_{1}}^{s_{1}^{\prime}-s_{1}}\ket{\psi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(p_{2})}\\\
=&\delta_{[p_{1}-k_{1}-q_{1}]_{1},0}\sum_{s}e^{i2\pi p_{1}s}\bra{\psi_{n\gamma
l}^{(\eta)}(k_{2})}\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}\hat{t}_{\mathbf{L}_{1}}^{s}\ket{\psi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(p_{2})}\\\
=&\delta_{[p_{1}-k_{1}-q_{1}]_{1},0}\sum_{s}e^{i2\pi(p_{1}-p_{2}\frac{L_{1y}}{L_{2}})s}e^{\frac{i}{2}s^{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-isK_{\eta^{\prime}l^{\prime},x}L_{1,x}}\bra{\psi_{n\gamma
l}^{(\eta)}(k_{2})}\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}\ket{\psi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(p_{2}-\frac{sp}{q})}\\\
=&\delta_{[p_{1}-k_{1}-q_{1}]_{1},0}\sum_{s}\delta_{p_{2}-sp/q,k_{2}+q_{2}}e^{i2\pi(p_{1}-p_{2}\frac{L_{1y}}{L_{2}})s}e^{\frac{i}{2}s^{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-isK_{\eta^{\prime}l^{\prime},x}L_{1,x}}\bra{\Phi_{n\gamma
l}^{(\eta)}(k_{2})}\hat{A}e^{-i\tilde{q}_{x}x}\ket{\Phi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(k_{2}+q_{2})}\\\
\equiv&\delta_{[p_{1}-k_{1}-q_{1}]_{1},0}\sum_{s}\delta_{p_{2}-sp/q,k_{2}+q_{2}}e^{i2\pi(p_{1}-p_{2}\frac{L_{1y}}{L_{2}})s}e^{\frac{i}{2}s^{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-isK_{\eta^{\prime}l^{\prime},x}L_{1,x}}\mathcal{A}^{\eta\eta^{\prime}}_{n\gamma
l,n^{\prime}\gamma^{\prime}l^{\prime}}(k_{2},k_{2}+q_{2}).\end{split}$ (S24)

The notation $[b]_{a}$ represents $b$ modulo $a$, with $a>0$. We have defined:

$\tilde{q}_{x}\equiv q_{x}+K_{\eta l,x}-K_{\eta^{\prime}l^{\prime},x}.$ (S25)

We proceed to work out the 1D integration by noting that:

$\bra{x}\ket{\Phi^{(\eta)}_{n\gamma l}(k_{2})}=\Phi^{(\eta)}_{n\gamma
l}(x+\tilde{k}_{2,\eta l}\ell^{2})=e^{ip_{x}\tilde{k}_{2,\eta
l}\ell^{2}}\Phi^{(\eta)}_{n\gamma l}(x)\equiv\hat{T}_{x}(\tilde{k}_{2,\eta
l}\ell^{2})\Phi^{(\eta)}_{n\gamma l}(x),$ (S26)

and as a result:

$\begin{split}\mathcal{A}^{\eta\eta^{\prime}}_{n\gamma
l,n^{\prime}\gamma^{\prime}l^{\prime}}(k_{2},k_{2}+q_{2})&=\bra{\Phi_{n\gamma
l}^{(\eta)}(0)}\hat{T}_{x}(-\tilde{k}_{2,\eta
l}\ell^{2})\hat{A}e^{-i\tilde{q}_{x}x}\hat{T}_{x}(\widetilde{(k+q)}_{2,\eta^{\prime}l^{\prime}}\ell^{2})\ket{\Phi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(0)}\\\
&=e^{i\tilde{q}_{x}\tilde{k}_{2,\eta l}\ell^{2}}\bra{\Phi_{n\gamma
l}^{(\eta)}(0)}\hat{A}e^{-i\tilde{q}_{x}x}e^{i\tilde{q}_{y}\ell^{2}p_{x}}\ket{\Phi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(0)}\\\
&=e^{i\tilde{q}_{x}\tilde{k}_{2,\eta
l}\ell^{2}}e^{\frac{i}{2}\tilde{q}_{x}\tilde{q}_{y}\ell^{2}}\bra{\Phi_{n\gamma
l}^{(\eta)}(0)}\hat{A}e^{c_{-}a+c_{+}a^{\dagger}}\ket{\Phi^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(0)}\\\
&\equiv e^{i\tilde{q}_{x}\tilde{k}_{2,\eta
l}\ell^{2}}e^{\frac{i}{2}\tilde{q}_{x}\tilde{q}_{y}\ell^{2}}\hat{\mathcal{A}}_{n\gamma
l,n^{\prime}\gamma^{\prime}l^{\prime}}^{\eta\eta^{\prime}},\end{split}$ (S27)

where $\hat{\mathcal{A}}^{\eta\eta^{\prime}}_{n\gamma
l,n^{\prime}\gamma^{\prime}l^{\prime}}$ can be expressed in terms of
associated Laguerre polynomials [47], and:

$\tilde{q}_{y}\equiv q_{y}+K_{\eta l,y}-K_{\eta^{\prime}l^{\prime},y},\
c_{\pm}=-i\frac{\ell}{\sqrt{2}}(\tilde{q}_{x}\mp i\tilde{q}_{y}).$ (S28)

Next we consider the implications of the $\delta$ function constraints.
Firstly we have:

$q_{1}=p_{1}-k_{1}+r,\ r\in\mathbb{Z}.$ (S29)

Secondly we have:

$q_{2}=p_{2}-k_{2}-\frac{sp}{q},s\in\mathbb{Z}.$ (S30)

This shows that above matrix elements are non-vanishing only if:

$\mathbf{q}\equiv\mathbf{p}-\mathbf{k}+r\mathbf{g}_{1}-\frac{sp}{q}\mathbf{g}_{2},\
r,s\in\mathbb{Z}.$ (S31)

Finally, we have the following expression for the matrix elements:

$\displaystyle\bra{V^{(\eta)}_{n\gamma
l}(\mathbf{k})}A_{\mathbf{q}}\ket{V^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(\mathbf{p})}$
$\displaystyle=$
$\displaystyle\delta_{[p_{1}-k_{1}-q_{1}]_{1},0}\sum_{s}\delta_{p_{2}-sp/q,k_{2}+q_{2}}e^{i2\pi(p_{1}-p_{2}\frac{L_{1y}}{L_{2}})s}e^{\frac{i}{2}s^{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-isK_{\eta^{\prime}l^{\prime},x}L_{1,x}}e^{i\tilde{q}_{x}\tilde{k}_{2,\eta
l}\ell^{2}}e^{\frac{i}{2}\tilde{q}_{x}\tilde{q}_{y}\ell^{2}}\hat{\mathcal{A}}^{\eta\eta^{\prime}}_{n\gamma
l,n^{\prime}\gamma^{\prime}l^{\prime}}.$

with $\tilde{\mathbf{q}}\equiv\mathbf{q}+\mathbf{K}_{\eta
l}-\mathbf{K}_{\eta^{\prime}l^{\prime}}$.

###  II.4 Matrix elements of strained BM Hamiltonian

####  II.4.1 Interlayer terms

The interlayer terms are a special case of the generic $A_{\mathbf{q}}$
discussed in the previous section. For bottom to top layer tunneling we have
(see Eq. (S24)):

$\begin{split}&\bra{V^{(\eta)}_{n\gamma
b}(\mathbf{k})}T_{\eta,j}e^{-i\eta\mathbf{q}_{j}\cdot\mathbf{r}}\ket{V^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}t}(\mathbf{p})}\\\
=&\delta_{\eta\eta^{\prime}}\delta_{[p_{1}-k_{1}-q_{1}]_{1},0}\sum_{s}\delta_{p_{2}-sp/q,k_{2}+q_{2}}e^{i2\pi(p_{1}-p_{2}\frac{L_{1y}}{L_{2}})s}e^{i\frac{s^{2}}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-isK_{\eta
t,x}L_{1,x}}e^{i\tilde{q}_{x}\tilde{k}_{2,\eta
l}\ell^{2}}e^{\frac{i}{2}\tilde{q}_{x}\tilde{q}_{y}\ell^{2}}\hat{\mathcal{T}}^{(\eta,j)}_{n\gamma
b,n^{\prime}\gamma^{\prime}t},\end{split}$ (S32)

where
$\mathbf{q}_{j}=\mathbf{0},\mathbf{g}_{2},\mathbf{g}_{1}+\mathbf{g}_{2}$. For
notational convenience we have also defined $\mathbf{q}_{j}\equiv
q_{1}\mathbf{g}_{1}+q_{2}\mathbf{g}_{2}$ where $q_{1,2}$ are integers (with
implicit dependence on $j$). Plugging these into the $\delta$-function
constraints we have:

$p_{1}=k_{1},\ p_{2}-k_{2}-\frac{sp}{q}\in\mathbb{Z}.$ (S33)

We can split the momentum $p_{2},k_{2}$ into strips of
$[0,1/q),\dots[(p-1)/q,p/q)$, by redefining:

$p_{2}\rightarrow p_{2}+r_{2}^{\prime}/q,\ k_{2}\rightarrow k_{2}+r_{2}/q,\
k_{2},p_{2}\in[0,1/q).$ (S34)

It is clear then that $(k_{1},k_{2})\in[0,1)\times[0,1/q)$ are good quantum
numbers under moiré periodic potential, and that:

$r_{2}^{\prime}=r_{2}+qq_{2}+sp\Rightarrow r_{2}=[r_{2}+qq_{2}]_{p},\
s=-(r_{2}+qq_{2})/p.$ (S35)

####  II.4.2 Narrow band eigenstates

From above analysis, it is clear that the narrow band eigenstates in the
Landau gauge are labeled by $(k_{1},k_{2})\in[0,1)\times[0,1/q)$, and an
additional $2q$ quantum numbers per valley and spin. Generally we can denote
them as:

$\ket{\Psi^{(\eta)}_{a}(\mathbf{k})}\equiv\sum_{n\gamma
l}\sum_{r_{2}=0}^{p-1}U_{n\gamma
lr_{2},a}^{(\eta)}(\mathbf{k})\ket{V^{(\eta)}_{n\gamma
l}(k_{1},k_{2}+r_{2}/q)},\ a=1,\dots,2q.$ (S36)

####  II.4.3 Degeneracy of magnetic subbands generated by
$\hat{t}_{\mathbf{L}_{2}}$

The non-interacting Hamiltonian $H_{\text{non-int}}$ in valley $\eta$ acting
on the narrow band eigenstates is:

$H_{\text{non-
int}}\ket{\Psi^{(\eta)}_{a}(\mathbf{k})}=\varepsilon^{(\eta)}_{a}(\mathbf{k})\ket{\Psi^{(\eta)}_{a}(\mathbf{k})}.$
(S37)

We apply the magnetic translation $\hat{t}_{\mathbf{L}_{2}}$ to these energy
eigenstates. Note that they are eigenstates of $\hat{t}_{\mathbf{L}_{1}}$ and
$\hat{t}^{q}_{\mathbf{L}_{2}}$; the reason why $\hat{t}_{\mathbf{L}_{2}}$ is
applied $q$ times is due to non-commuting nature of magnetic translation
operators $\hat{t}_{\mathbf{L}_{1}}$ and $\hat{t}_{\mathbf{L}_{2}}$. We see
that:

$\begin{split}\hat{t}_{\mathbf{L}_{2}}\ket{\Psi^{(\eta)}_{a}(\mathbf{k})}&=\sum_{n\gamma
l}\sum_{r_{2}=0}^{p-1}U^{(\eta)}_{n\gamma
lr_{2},a}(\mathbf{k})\hat{t}_{\mathbf{L}_{2}}\ket{V_{n\gamma
l}^{(\eta)}(k_{1},k_{2}+\frac{r_{2}}{q})}\\\ &=\sum_{n\gamma
l}\sum_{r_{2}=0}^{p-1}U^{(\eta)}_{n\gamma
lr_{2},a}(\mathbf{k})e^{-i2\pi(k_{2}+r_{2}/q)}\ket{V_{n\gamma
l}^{(\eta)}(k_{1}+p/q,k_{2}+\frac{r_{2}}{q})}.\end{split}$ (S38)

Furthermore since $[H_{\text{non-int}},\hat{t}_{\mathbf{L}_{2}}]=0$, we have:

$H_{\text{non-
int}}\hat{t}_{\mathbf{L}_{2}}\ket{\Psi^{(\eta)}_{a}(\mathbf{k})}=\hat{t}_{\mathbf{L}_{2}}H_{\text{non-
int}}\ket{\Psi^{(\eta)}_{a}(\mathbf{k})}=\epsilon^{(\eta)}_{a}(\mathbf{k})\hat{t}_{\mathbf{L}_{2}}\ket{\Psi^{(\eta)}_{a}(\mathbf{k})},$
(S39)

namely that $\hat{t}_{\mathbf{L}_{2}}\ket{\Psi^{(\eta)}_{a}(\mathbf{k})}$ is
an energy eigenstate at $([k_{1}+p/q]_{1},k_{2})$.

This means that the tower of states are periodic with respect to
$k_{1}\rightarrow k_{1}+1/q$, and that all the distinct energy solutions can
be found in the domain of $(k_{1},k_{2})\in[0,1/q)\times[0,1/q)$. Furthermore,
the $\hat{t}_{\mathbf{L}_{2}}$ translation gives us one way of gauge fixing
between $k_{1}$ values in different strips of $[0,1/q)$. Specifically, we can
define:

$\ket{\Psi^{(\eta)}_{a}([k_{1}+p/q]_{1},k_{2})}\equiv e^{i2\pi
k_{2}}\hat{t}_{\mathbf{L}_{2}}\ket{V_{a}(k_{1},k_{2})},$ (S40)

as the gauge-fixed eigenstate wavefunction at $([k_{1}+p/q]_{1},k_{2})$. This
amounts to the following definition:

$U^{(\eta)}_{n\gamma lr,a}([k_{1}+r_{1}p/q]_{1},k_{2})=e^{-i2\pi
r_{1}(r_{2}/q)}U^{(\eta)}_{n\gamma lr,a}(k_{1},k_{2}),\ r_{1}=0,\dots q-1.$
(S41)

###  II.5 Matrix elements of
$A_{\mathbf{q}}=\hat{A}e^{-i\mathbf{q}\cdot\mathbf{r}}$ with respect to narrow
band eigenstates

This can be computed as follows:

$\begin{split}&\bra{\Psi^{(\eta)}_{a}(\mathbf{k})}A_{\mathbf{q}}\ket{\Psi^{(\eta^{\prime})}_{b}(\mathbf{p})}\\\
=&\sum U^{(\eta)*}_{n\gamma
lr_{2},a}(\mathbf{k})U^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}r_{2}^{\prime},b}(\mathbf{p})\bra{V^{(\eta)}_{n\gamma
l}(k_{1},k_{2}+r_{2}/q)}A_{\mathbf{q}}\ket{V^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}}(p_{1},p_{2}+r_{2}^{\prime}/q)}\\\
=&\sum U^{(\eta)*}_{n\gamma
lr_{2},a}(\mathbf{k})U^{(\eta^{\prime})}_{n^{\prime}\gamma^{\prime}l^{\prime}r_{2}^{\prime},b}(\mathbf{p})\delta_{[p_{1}-k_{1}-q_{1}]_{1},0}\\\
&\times\sum_{s}\delta_{p_{2}-sp/q,k_{2}+q_{2}}e^{i2\pi(p_{1}-p_{2}\frac{L_{1y}}{L_{2}})s}e^{i\frac{s^{2}}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-isK_{\eta^{\prime}l^{\prime},x}L_{1,x}}e^{i\tilde{q}_{x}\tilde{k}_{2,\eta
l}\ell^{2}}e^{\frac{i}{2}\tilde{q}_{x}\tilde{q}_{y}\ell^{2}}\hat{\mathcal{A}}^{\eta\eta^{\prime}}_{n\gamma
l,n^{\prime}\gamma^{\prime}l^{\prime}}.\end{split}$ (S42)

Note that here the $\delta$-constraint depends on $r_{2},r_{2}^{\prime}$, and
matrix elements are only non-vanishing for:

$\mathbf{q}=\mathbf{p}-\mathbf{k}+\frac{r_{2}^{\prime}-r_{2}}{q}\mathbf{g}_{2}+r\mathbf{g}_{1}-\frac{sp}{q}\mathbf{g}_{2}.$
(S43)

Computing above bruteforce is costly. However we could make use of
$\hat{t}_{\mathbf{L}_{2}}$ to reduce the amount of calculations by a factor of
$q^{2}/(2q-1)$. This is seen by noting that:

$\begin{split}&\bra{V^{(\eta)}_{n\gamma
l}([k_{1}+r_{1}p/q]_{1},k_{2})}A_{\mathbf{q}}\ket{V^{(\eta)}_{n\gamma
l}([p_{1}+r_{1}p/q]_{1},p_{2})}\\\ =&e^{i2\pi
r_{1}(p_{2}-k_{2})}\bra{V^{(\eta)}_{n\gamma
l}(\mathbf{k})}\hat{t}_{\mathbf{L}_{2}}^{-r_{1}}A_{\mathbf{q}}\hat{t}_{\mathbf{L}_{2}}^{r_{1}}\ket{V^{(\eta)}_{n^{\prime}\gamma^{\prime}l^{\prime}}(\mathbf{p})}\\\
=&e^{i2\pi r_{1}(p_{2}-k_{2}-q_{2})}\bra{V^{(\eta)}_{n\gamma
l}(\mathbf{k})}A_{\mathbf{q}}\ket{V^{(\eta)}_{{n^{\prime}\gamma^{\prime}l^{\prime}}}(\mathbf{p})}\\\
=&e^{-i2\pi r_{1}n/q}\bra{V^{(\eta)}_{n\gamma
l}(\mathbf{k})}A_{\mathbf{q}}\ket{V^{(\eta)}_{n^{\prime}\gamma^{\prime}l^{\prime}}(\mathbf{p})},\
n\in\mathbb{Z}.\end{split}$ (S44)

In the last step we made use of the constraint $q_{2}=p_{2}-k_{2}+n/q$.
Therefore, rather than calculating $q^{2}$ blocks of matrices corresponding to
differences in $k_{1},p_{1}$, we only need to compute $2q-1$ blocks, and
generate all remaining ones by multiplying a global phase factor.

##  III Self-consistent Hartree-Fock method in finite magnetic field

We proceed to discuss the Hartree-Fock procedure in a finite magnetic field
(B-SCHF). We consider the strained BM Hamiltonian in the presence of screened
Coulomb interactions, and project it onto the finite field Hilbert space. The
Hamiltonian is given by:

$H=\sum_{\eta
sa,\mathbf{k}}\varepsilon^{(\eta)}_{sa}(\mathbf{k})d^{\dagger}_{\eta
sa,\mathbf{k}}d_{\eta
sa,\mathbf{k}}+\frac{1}{2A}\sum_{\mathbf{q}}V_{\mathbf{q}}\delta\hat{\rho}_{\mathbf{q}}\delta\hat{\rho}_{-\mathbf{q}},$
(S45)

where $A$ is the area of the system, $\varepsilon^{(\eta)}_{sa}(\mathbf{k})$
is the Hofstadter spectra including spin Zeeman splitting, and
$V_{\mathbf{q}}$ is the Fourier transform of the screened Coulomb interaction.
The projected density operator is given by:

$\delta\hat{\rho}_{\mathbf{q}}=\sum_{\eta
s}\bra{\Psi^{(\eta)}_{a}(\mathbf{k})}\hat{\rho}e^{-i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi^{(\eta)}_{b}(\mathbf{p})}\left(d^{\dagger}_{\eta
sa,\mathbf{k}}d_{\eta
sb,\mathbf{p}}-\frac{1}{2}\delta_{a,b}\delta_{\mathbf{k},\mathbf{p}}\right)$
(S46)

where $\hat{\rho}=\hat{I}$. We emphasize again that the matrix elements are
non-zero only for wavevectors satisfying:

$\mathbf{q}=\mathbf{p}-\mathbf{k}+\frac{r_{2}^{\prime}-r_{2}}{q}\mathbf{g}_{2}+r\mathbf{g}_{1}-\frac{sp}{q}\mathbf{g}_{2}.$
(S47)

Here on we introduce the matrix notation for the structure factor:

$\hat{\Lambda}^{\eta
s}_{\mathbf{q}}(a\mathbf{k},b\mathbf{p})\equiv\bra{\Psi^{(\eta)}_{a}(\mathbf{k})}\hat{\rho}e^{-i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi^{(\eta)}_{b}(\mathbf{p})}.$
(S48)

where we use $s$ to label spin. Pay attention that in prior contexts we have
used $s$ to denote an integer.

As a result, the Coulomb piece of the Hamiltonian can be written as:

$H_{U}=\frac{1}{2A}\sum_{\mathbf{q}}V_{\mathbf{q}}\hat{\Lambda}_{\mathbf{q}}^{\eta
s}(a\mathbf{k},b\mathbf{p})\hat{\Lambda}_{\mathbf{q}}^{\eta^{\prime}s^{\prime}*}(\beta\mathbf{p}^{\prime},\alpha\mathbf{k}^{\prime})\left(d^{\dagger}_{\eta
sa,\mathbf{k}}d_{\eta
sb,\mathbf{p}}-\frac{1}{2}\delta_{a,b}\delta_{\mathbf{k},\mathbf{p}}\delta_{\mathbf{q},\mathbf{G}}\right)\left(d^{\dagger}_{\eta^{\prime}s^{\prime}\alpha,\mathbf{k}^{\prime}}d_{\eta^{\prime}s^{\prime}\beta,\mathbf{p}^{\prime}}-\frac{1}{2}\delta_{\alpha,\beta}\delta_{\mathbf{k}^{\prime},\mathbf{p}^{\prime}}\delta_{\mathbf{q},\mathbf{G}}\right),$
(S49)

where $\mathbf{G}=m\mathbf{g}_{1}+n\mathbf{g}_{2}$ is the moiré reciprocal
lattice vector. We note first on the constraint of the $\delta$-functions on
the allowed values of momentum. Specifically we have:

$\mathbf{q}=\mathbf{p}-\mathbf{k}+m\mathbf{g}_{1}+\frac{n}{q}\mathbf{g}_{2}=\mathbf{k}^{\prime}-\mathbf{p}^{\prime}+m^{\prime}\mathbf{g}_{1}+\frac{n^{\prime}}{q}\mathbf{g}_{2}.$
(S50)

Since $p_{1}-k_{1}\in(-1,1)$ and $p_{2}-k_{2}\in(-1/q,1/q)$, the following
constraint applies:

$\mathbf{p}-\mathbf{k}=\mathbf{k}^{\prime}-\mathbf{p}^{\prime},\
m=m^{\prime},\ n=n^{\prime}.$ (S51)

###  III.1 Product state and one-particle density matrix

At a given filling, we consider a ground state $\ket{\Omega}$ to be a product
state, given by partial fillings of the narrow band states $\ket{\Psi^{\eta
s}_{a}(\mathbf{k})}$. This is equivalent to defining a one-particle density
matrix:

$\hat{Q}^{\eta
s,\eta^{\prime}s^{\prime}}_{a,b}(\mathbf{k})\equiv\bra{\Omega}d^{\dagger}_{\eta
sa,\mathbf{k}}d_{\eta^{\prime}s^{\prime}b,\mathbf{k}}\ket{\Omega}-\frac{1}{2}\delta_{\eta,\eta^{\prime}}\delta_{s,s^{\prime}}\delta_{a,b}.$
(S52)

Note that $\hat{t}_{\mathbf{L}_{2}}$ takes $\mathbf{k}$ to
$\mathbf{k}+\phi/\phi_{0}\mathbf{g}_{1}$, and as a result, a difference of the
density matrix at $\hat{Q}(\mathbf{k})$ and
$\hat{Q}(\mathbf{k}+\phi/\phi_{0}\mathbf{g}_{1})$ allows us to probe magnetic
translation symmetry breaking of $\hat{t}_{\mathbf{L}_{2}}$.

###  III.2 Hartree Fock Hamiltonian

The Hartree-Fock Hamiltonian is given as:

$\begin{split}H^{HF}&=\frac{1}{A}\sum_{\mathbf{q}}V_{\mathbf{q}}\hat{\Lambda}_{\mathbf{q}}^{\eta
s}(a\mathbf{k},b\mathbf{k})\hat{\Lambda}_{\mathbf{q}}^{\eta^{\prime}s^{\prime}*}(\beta\mathbf{k}^{\prime},\alpha\mathbf{k}^{\prime})\hat{Q}^{\eta^{\prime}s^{\prime},\eta^{\prime}s^{\prime}}_{\alpha,\beta}(\mathbf{k}^{\prime})\left(d^{\dagger}_{\eta
sa,\mathbf{k}}d_{\eta sb,\mathbf{k}}-\frac{1}{2}\delta_{a,b}\right)\\\
&-\frac{1}{A}\sum_{\mathbf{q}}V_{\mathbf{q}}\hat{\Lambda}_{\mathbf{q}}^{\eta
s}(a\mathbf{k},b\mathbf{p})(\hat{Q}^{T})^{\eta
s,\eta^{\prime}s^{\prime}}_{b,\alpha}(\mathbf{p})\hat{\Lambda}_{\mathbf{q}}^{\eta^{\prime}s^{\prime}\dagger}(\alpha\mathbf{p},\beta\mathbf{k})\left(d^{\dagger}_{\eta
sa,\mathbf{k}}d_{\eta^{\prime}s^{\prime}\beta,\mathbf{k}}-\frac{1}{2}\delta_{a,\beta}\delta_{\eta,\eta^{\prime}}\delta_{s,s^{\prime}}\right).\end{split}$
(S53)

Our B-SCHF results are obtained by minimizing Eq. S53 with respect to product
state $\ket{\Omega}$. In practice, we achieve minimization by running both
random initializations for the density matrix $\hat{Q}$ as well as for
educated guesses, such as a strong coupling Chern insulating state, flavor
polarized state, intervalley coherent state, etc. Typically, the phase diagram
obtained in Fig. 1 of the main text is a collection of up to 10
initializations schemes of the density matrix.

For each magnetic flux ratio $\phi/\phi_{0}$, we study momentum space mesh
consisting of $n_{q}q$ points along $\mathbf{g}_{1}$ direction in the range of
$[0,1)$, and $n_{q}$ points along $\mathbf{g}_{2}$ direction in the range of
$[0,1/q)$, where $n_{q}$ is chosen as the maximum integer such that
$n_{q}q\leq 14$. Our calculations presented in the main text correspond to 23
magnetic flux ratios between $1/12$ and $1/2$, where we take $q\leq 12$ and
$p<q$.

##  IV Extended Figures

Figure S1:  Non-interacting Hofstadter spectra and Wannier plots for various
twist angles with and without heterostrain. In the Wannier plots, bright
colors correspond to low density of states (gapped) and darker blue colors
correponsd to high density of states (compressible). In obtaining the Wannier
plots we used an energy broadening factor of $\gamma=0.1$meV except for the
magic angle ($\theta=1.05^{\circ}$) without strain, where we used
$\gamma=0.01$meV.  Figure S2:  B-SCHF results for the single-electron
excitation gap as a function of filling ($n/n_{s}$) and magnetic flux ratio
($\phi/\phi_{0}$) for various twist angles in the absence of heterostrain.
Figure S3:  Averaged $\langle\sigma_{z}\tau_{z}\rangle$ for the occupied
electronic states at $\phi/\phi_{0}$, twist angle $\theta=1.05^{\circ}$ in the
absence of heterostrain, for different ratios of $w_{0}/w_{1}$. As
$w_{0}/w_{1}$ is decreased, the correlated Hofstadter ferromagnets persist
toward lower angles, before losing to other gapped states with intervalley
coherence. The phase transition is first order, and marked by a collapse of
$\langle\sigma_{z}\tau_{z}\rangle$. Figure S4:  Density matrices of an
occupied flavor (valley/spin) for the correlated Hofstadter ferromgnets in the
eigenbasis of $(0,-4)$ at flux ratio $\phi/\phi_{0}=2/5$. Here we show that
the CHFs are predominantly flavor symmetry breakings within the renormalized
magnetic subbands obtained at $(0,-4)$. Figure S5:  Full density matrices for
the intervalley coherent states at $(-3,-1)$, $(-2,-2)$ and $(-3,-1)$ (upper
panels), $(-3,0)$, $(-2,-1)$ and $(-1,-2)$ (lower panels). Results for
$\theta=1.05^{\circ}$, $\phi/\phi_{0}=1/8$, and with heterostrain. Distinct
intervalley hybridizations are established for the upper versus lower panels.
Figure S6:  (a) Chemical potential $\mu$ versus filling $n/n_{s}$ for
$\theta=1.05^{\circ}$ and with heterostrain. Results obtained for $16\times
16$ mesh.(b) 3D plot of Hartree-Fock renormalized energy dispersions at
$n/n_{s}=0$ for the same parameterization, showing semimetallic behavior.
Dispersion obtained for $25\times 25$ mesh. Figure S7:  Calculated local
density of states $\mathcal{N}(E)$ at $\phi/\phi_{0}=1/8$ for various twist
angles and with heterostrain. $\mathcal{N}(E)$ is normalized with respect to
the maximum value. $E=0$ corresponds to the position of the chemical potential
$\mu$ at a given filling $n/n_{s}$. Here we used an energy broadening factor
of $\gamma=1\mathrm{meV}$. Dark arrows in (e) and (f) mark $(s,t)=(0,\pm
4),\pm(1,3),\pm(2,2),\pm(3,1)$.