Theory of correlated Chern insulators in twisted bilayer graphene

Magic-angle twisted bilayer graphene is the best studied physical platform featuring moire 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\'e 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.


I. INTRODUCTION
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 ∼ 1.05 • [3], TBG devices show a plethora of correlated phenomena including superconductivity, correlated insulating states and (quantized) anomalous Hall effect .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 [27,28].The narrow band Hilbert space can thus be decomposed into a Chern +1 and a Chern −1 band [27,29,30].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 ±1 insulator [31][32][33][34].Indeed, experiments have observed anomalous Hall effect (AHE) near the filling of 3 electrons per moiré unit cell [10,20] in hBN aligned samples.Further studies on non-aligned samples [21,25] 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 [30,33,35,36].
However, CCIs are also observed in TBG devices away from the magic angle (∼ 1.27 • ), where the bandwidth of the narrow bands is expected to be significantly larger than at 1.05 • , 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 [7,15,16,23], termed Hofstadter subband ferromagnets (HSFs) [15].As argued theoretically [37][38][39], 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 zCI, HSF or other states has been carried out at B ̸ = 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 (IVC) states [29,33,[40][41][42], the incommensurate Kekulé spiral (IKS) orders [38], and the stripe and nematic states [14,30,43,44], remains unclear.
Here we report the first comprehensive study of the interacting electrons within the TBG narrow bands directly at B ̸ = 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, ±4), ±(1, 3), ±(2, 2), and ±(3, 1), as shown in the Fig. 1 for the case with heterostrain and Fig. 2 for the case without heterostrain.These figures plot the single particle excitation gap at the Fermi level obtained using the self-consistent Hartree-Fock method for each electron density and each magnetic field that we studied.In either strained or unstrained case, CCIs are found to be stablized at higher B fields for twist angles as high as 1.38 • (highest twist angle studied in this work).Based on an analysis of their wavefunctions, we identify them as correlated Hofstadter ferromagnets (CHFs).Similar to HSFs, the CHFs correspond to selective population of the valley and spin flavors, but of the interaction-renormalized magnetic subbands (see Fig. 3).Unlike HSFs however, CHFs may include -but aren't limited to-spin and/or valley polarized states which correspond to B = 0 Chern insulators whose interaction renormalized bands are Landau quantized at B ̸ = 0.Although only metastable at B = 0, such Chern insulators can be stabilized at B ̸ = 0 in the form of CHFs as we demonstrate in the Fig. 6(d-e).
For realistic heterostrain (Fig. 1), upon lowering B and at a non-zero s we find a phase transition into incompressible states with intervalley coherence.These states break the discrete magnetic translation symmetry but preserve the combination of the discrete magnetic translation and a U v (1) valley transformation (see Sec. III B), and therefore are the finite B analogs of the IKS states [38], albeit carrying a non-zero Chern number t. Upon lowering B and at larger twist angles, the IKS states transition into nearly compressible states; at lower twist angles they remain robust down to the lowest magnetic flux studied in this work.
In the absence of heterostrain (Fig. 2) and at larger twist angles we instead observe CHFs transitioning directly into nearly compressible states upon lowering B without the intermediate IKS states.As we lower the twist angle toward 1.05 • , the incompressible state at ±(3, 1) 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).Because there is no symmetry distinction between them, there is no true phase transition between CHF and the sCI, rather it is a smooth crossover as shown in the Fig. 4(c,d).The ±(2, 2), ±(1, 3) and (0, ±4) states also approach sCIs upon decreasing the twist angle, but they experience a first order phase transition into IVCs at low twist angles [29,33,[40][41][42].IVCs break the U v (1) valley symmetry, but unlike the IKSs, they preserve the magnetic translation symmetry while also strongly hybridize the entirety of the narrow band Hilbert space.The details of this transition depend sensitively on the model parameters, as shown in the Fig. S19.
Conversely, in the presence of heterostrain, the sCIs are absent in the Fig. 1, as evidenced by the fact that the sublattice polarization remains low and far from saturating the sCI bound (Fig. 4 a,b).At the magic angle with heterostrain, we also identify Chern 0 IKS states along (−3, 0) and (−2, 0) (Fig. 1f), consistent with what has been reported previously in B = 0 Hartree-Fock calculations [38].Unlike IKS states at a non-zero t, these states are less robust upon increasing B and lose to nearly compressible states.Similar phenomenon has been reported experimentally in the Ref. [15,21,22,45].
In addition to these prominent CCIs, we also find gapped Chern states emanating from band edges and the charge neutrality point (CNP) as shown in the Figs.1-2.In the presence of heterostrain, they are quantum Hall ferromagnetic states (QHFMs) corresponding to spinvalley symmetry breaking states within a multi-flavored Landau level, with more details discussed in Figs.S6 and S7.Without heterostrain, gapped states emanating from the band edges are also identified as QHFMs (Fig. S13).However, the gapped states emanating from the CNP assume a different character (Fig. S14), developing intervalley coherence and a correlation gap at B = 0 even at the largest twist angle studied.In contrast, QHFMs are field-induced symmetry breaking states, and are absent at B = 0.At larger twist angles, a quantum spin Hall insulating state (QSH) along (−2, 0) is identified (see orange dashed line in Figs.1a and 2a), whose origin can be traced back to the spin-split non-interacting magnetic subbands at half flux quantum per moiré unit cell (Figs.S1 and S11).CCIs extrapolating to a fractional s at B = 0 can also be found at 1.05 • with heterostrain and for all twist angles in the absence of heterostrain.Further examination of their many body wavefunctions reveal that the magnetic translation symmetries are broken (Figs. 7, S9, S16), similar to the symmetry-broken Chern insulators (SBCIs) discussed in relation to the experiments of Ref. [14].

II. MODEL AND METHOD
We perform self-consistent Hartree-Fock analysis (B-SCHF) at B ̸ = 0 using the minimal continuum Hamiltonian (BM) [3,46], with Coulomb interactions projected onto the narrow band Hilbert space.Here we briefly outline the formalism, additional details are in the Sup-plementary Information (SI).Our starting point is the (strained) BM Hamiltonian at rational magnetic flux ratios ϕ/ϕ 0 = p/q where p and q are coprime integers, ϕ is the magnetic flux per moiré unit cell and ϕ 0 = h/e is the magnetic flux quantum.We choose a Landau gauge such that the magnetic vector potential A(r) = eB ŷ, where ŷ is defined along L 2 direction, with L i=1,2 the two (strain-deformed) moiré unit cell vectors.The interacting Hamiltonian is invariant under discrete magnetic translation symmetries tL1 (r) = e −i2π ϕ ϕ 0 y |L 2 | TL1 and tL2 = TL2 , where TLi=1,2 denote discrete moiré translations at B = 0. tL1 and tL2 are non-commuting but satisfy tL1 , tq L2 = 0.This allows us to define a magnetic Brillouin zone described by the magnetic crystal momentum k = k 1 g 1 + k 2 g 2 , with k 1 ∈ [0, 1) and k 2 ∈ [0, 1/q), and g i=1,2 are moiré reciprocal lattice vectors (for more detailed information see SI Sec.I and II).We therefore first solve for the non-interacting Hofstadter spectra ε ηsr (k) and associated eigenstates |Ψ ηsr (k)⟩, with η = K, K ′ and s =↑, ↓ denoting valley and spin quantum numbers, and r = 1, . . .2q is the magnetic subband index.Spin Zeeman splitting is also considered in this calculation.In earlier works [47,48], we used the hybrid Wannier states (hWS) at B = 0 to construct the finite B Hilbert space.Although accurate and numerically efficient at low B, the hWS approach to TBG was shown to break down above moderate flux ratios (e.g., ϕ/ϕ 0 ≳ 0.2) [47].As one of the purposes of this work is to connect the CCIs between low and high B, we instead solve the BM Hamiltonian by expanding in Landau level basis of each graphene layer.By scaling the upper Landau level cutoff as ϕ/ϕ 0 decreases, we ensure an accurate construction of the Hilbert space as well as the non-interacting Hofstadter spectra.While the high B regime is straightforward in the LL basis, the low B is computationally expensive.We were nevertherless able to reach ϕ/ϕ 0 = 1/12 using the LL basis, which corresponds to ≈ 2.2Tesla at the magic angle (for all twist angle studied here we reached B < 4Tesla).This value is therefore sufficiently low to make direct comparison with experiments.The (interacting) results shown in Figs.1-2 are for 1/12 ≤ ϕ/ϕ 0 ≤ 1/2, with the maximum value of q being 12 and 1 < p < q.
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: Here A is the total area of the system, d ηsr,k is the electron annihilation operator, δ ρq is the Fourier transform of the electron density operator projected onto the narrow bands, subtracting a background charge density [47,49,50].It is given by: We consider a dual-gate screened Coulomb interaction of the form V (q) = 2πe 2 ϵ0ϵr|q| tanh |q|ξ

2
, with relative dielectric constant ϵ r = 15 and screening length ξ = 4 |L 1 ||L 2 |.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,7,[51][52][53] (see also Fig. S2).
In the B-SCHF procedure, we minimize the total energy for a fixed particle number using many body wavefunctions expressable as product states: where q 0 is an arbitrary wavevector shift between single electron states in opposite valleys, whose value is constrained by the discrete momentum mesh such that k+q 0 is on the same momentum mesh as k (modulo a reciprocal lattice vector).The constrained product ′ n,k is over all the occupied states labeled by n, k. {α s ′ r ′ ,k+q0 } are variational parameters that minimize the total energy, and they obey sr |α for any {n, k}.The total energy is also optimized with respect to q 0 , allowing us to probe IKS-like states (see SI Sec.III D).An equivalent formulation of the B-SCHF procedure (for details see SI Sec.III) is based on the one-particle density matrix: where for notational convenience we defined d † . Note that Q contains information about q 0 .In the remaining text we discuss the numerical results based on either the oneparticle density matrix or the many body wavefunction whichever is more convenient.
The projected Hamiltonian at B ̸ = 0 is invariant under the following set of symmetries [49,54,55]: C 2z , valley U v (1) and spin U s (1), many-body particle-hole P , magnetic translation symmetries generated by tL1 and tL2 .We fix the gauge such that The above gauge choice is useful in identifying a magnetic translation symmetry breaking from the density matrix (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.The B-SCHF calculation is carried out for a range of twist angles from 1.38 • to 1.05 • , both for unstrained model as well as for realistic uniaxial heterostrain magnitude ϵ = 0.2% and orientation φ = 0 • (see SI and Ref. [39] for the definition of the uniaxial strain orientation).We choose the Fermi velocity such that ℏv F /a = 2482meV with the graphene lattice constant a ≈ 2.46 Å, and interlayer tunneling parameters w 0 = 77meV (intrasublattice) and w 1 = 110meV (inter-sublattice).These parameters place the magic angle near 1.05 • .The respective non-interacting Hofstadter spectra and Wannier diagrams with/without heterostrain are shown in Fig. S1  and S11.
At a given B and electron density, obtaining the true Hartree-Fock ground state is a highly non-trivial task due to competing states of similar energy.We typically run ∼ 6 random initializations of the single-particle density matrix, as well as several educated guesses (e.g., flavor polarized or intervalley coherent states), and report the lowest-energy state as the ground state.Optimal damping algorithm is also used to speed up the numerical convergence [29].This elaborate procedure turns out to be adequate (i.e.random and educated initializations converge to the same state) for establishing incompressible ground states with big excitation gaps, such as the CCIs with (s, t) = (0, ±4) ±(1, 3), ±(2, 2) and ±(3, 1), QHFMs, and CCIs extrapolating to a fractional s.However it may face convergence issues for nearly compressible states, which are abundant in the phase diagrams shown in Figs. 1 and 2. We make no assertions regarding the nature of such nearly compressible states in this paper and mainly use them in order to highlight the contrast with the ground states with large excitation gaps or, at B = 0, to elucidate the physics of the B-induced incompressible states as in the Fig. 6(a-c).We often find competing states close in energy, with ≲ 1meV difference in the Hartree-Fock energy per moiré unit cell.Throughout the paper we try to adopt the philosophy of identifying these competing states [2], and comment on how small variations in model parameters (e.g., dielectric constant, kinetic terms beyond BM Hamiltonian) may tip the balance between ground states and metastable states.

A. Phase diagram and main CCIs
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 (i.e. the charge gap) as a function of moiré unit cell filling (n/n s ) and magnetic flux ratio (ϕ/ϕ 0 ) for six twist angles 1.38 • , 1.32 • , 1.28 • , 1.24 • , 1.20 • and 1.05 • .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 by red dashed lines in Fig. 1(c).At larger twist angles, CCIs along (−3, −1), (−2, −2), (−1, −3) emerge at high ϕ/ϕ 0 , and are replaced by nearly compressible states at lower ϕ/ϕ 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.
To better elucidate their nature, we compile detailed results for a representative twist angle 1.20 • in Fig. 3. Results for other twist angles can be found in the SI.Fig. 3(a) shows the non-interacting spectra of valley K and for one spin component (neglecting Zeeman effect).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 either from the zeroth Landau levels (zLLs) of the energetically split Dirac points (ϕ/ϕ 0 ≳ 0.1) or the ±1 LLs (ϕ/ϕ 0 ≲ 0.1), and carry Chern number +1 each, such that the total Chern number of all magnetic subbands is zero.Fig. 3(b) shows the single-particle spectra including Coulomb interactions along (s, t) = (−3, −1), where the occupied states are colored according to their valley polarization |⟨τ z ⟩|, where τ z is the Pauli matrix acting on valley degrees of freedom.|⟨τ z ⟩| → 1(0) implies maximal valley polarization (maximal intervalley mixing).At ϕ/ϕ 0 ≳ 1/3, the occupied states are maximally valley polarized, and have a large overlap onto the states marked in red in Fig. 3(a).To quantify the overlap we make use of the density matrix defined in Eq. ( 4), which has the spinvalley diagonal form Q (ηs) r,r ′ (k)δ η,η ′ δ s,s ′ for a state with unbroken valley and spin symmetries.A representative |Q K↑ r,r ′ (0)| is shown in Fig. 3(f).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 number -1 marked red in Fig. 3(a).For (−2, −2) and (−1, −3) and at higher B, the CCIs are also valley and spin polarized, similarly mostly populating the lower Chern -1 group of magnetic subbands for the specified valley and spin, with |Q (ηs) (0)| identical for all occupied flavors.Although these CCIs are closely related to the HSFs discussed in Ref. [7,15,23], the band structure renormalization is apparent in the non-vanishing off-diagonal matrix elements of Q (ηs) r,r ′ (k), signifying hybridization with the higher energy subbands (marked by grey in Fig. 3(a)).For this reason we refer to them as CHFs.As further demonstrated in Fig. S3, 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. S2 and also Ref. [38]).
At lower B, the valley and spin polarized CCIs along (−3, −1), (−2, −2) and (−1, −3) all transition into gapped states with strong intervalley mixing, as reflected by |⟨τ z ⟩| → 0 for the occupied states shown in Figs.3(be).Along (−3, −1) and (−1, −3), the intervalley mixing is between the same spin species (Fig. 3g), and the resulting states do not suffer from the Zeeman energy cost compared to CHFs which have the same spin polarization.However, along (−2, −2) the intervalley mixing is between opposite spins (mixing between the same spins is achieved as a metastable state), leading to an extra Zeeman energy cost.This likely explains the lowered critical field for the transition between the CHF and the IKS states for (−2, −2) compared to (−3, −1) and (−1, −3), as seen in Figs.3(b-e) by the valley polarization of occupied single-electron states.Based on a detailed analysis of the density matrix we identify these intervalley coherent CCIs as the finite B analog of IKS states [38] carrying a finite Chern number.We postpone a detailed discussion of identifying IKS states to Sec.III B. Finally in Fig. 3h we show the Hartree-Fock energy difference between the ground state and the CHFs which become metastable at lower B. It shows that the phase transition between CHF and IKS is likely first order along (−3, −1) and (−1, −3), where the IKS order parameter -qualitatively captured by the | QKsr,K ′ sr ′ (k)|, has an abrupt onset upon lowering B. On the other hand, the transition along (−2, −2) is most likely second order, as | QKsr,K ′ sr ′ (k)| gradually increases as B decreases.
The qualitative picture described above is universal across all the twist angles we have studied (see Fig. S4).However at larger twist angles, in addition to the phase transition between CHF and IKS, upon further decreasing B we observe a transition from a gapped IKS into a nearly compressible state, e.g., along (−3, −1) at 1.38 • as shown in Fig. 1(a).This can be attributed to the larger non-interacting bandwidth and comparatively weaker Coulomb interaction.

B. IKS states
In previous Hartree-Fock studies at B = 0, a notable finding is that in the presence of a small amount of uniaxial heterostrain (e.g., 0.2%), the ground state at integer filling fractions are the incommensurate Kekulé spiral ordered states [38].Unlike the IVC states achieved when heterostrain is absent [29,40], IKS states predominantly mix the lower energy bands (instead of the full narrow band Hilbert space) between opposite valleys of the non-interacting BM Hamiltonian, while developing a wavevector Q IKS incommensurate with the underlying moiré lattice.On the moiré scale, a modified discrete translation symmetry is preserved, i.e., discrete moiré translations ( TLj=1,2 ) followed by a U v (1) valley transformation: In our calculations, we find that the gapped states along (−3, −1), (−2, −2) and (−1, −3) at lower B [Fig. 3(b-e,g)] also hybridize electronic states predominantly of the lower Chern −1 group of the non-interacting magnetic subbands, break the magnetic translation symmetries generated by tLj=1,2 , but preserve a modified translation e − i 2 (QIKS•Lj )τz tLj .This further constrains the many body wavefunction in Eq. ( 3) and the density matrix in Eq. ( 4), with details discussed in SI Sec.III.D. Therefore, we can identify them as the finite B analog of the zero B IKS states.At larger twist angles, we only find gapped IKS states with a finite Chern number (i.e.non-zero t), and an absence of Chern 0 IKS states.The latter only occur close to the magic angle [e.g.Fig. 1(e,f)] along (−3, 0) and (−2, 0).At higher B, these Chern 0 IKS states become energetically unfavorable and lose to nearly compressible states.This phenomenon has been reported in various experiments [15,21,22,45].

C. Robustness of the CHF-IKS phase transition against perturbations
Given the small energy differences between the CHF and IKS states along (−3, −1), (−2, −2) and (−1, −3) shown in Fig. 3(h), it is natural to ask how robust is this phase transition to small variations in model parameters.Although this is difficult to answer definitively, it is possible to present some qualitative arguments favoring such a phase transition.Along (−2, −2), as argued above, Zeeman coupling is expected to favor a CHF state, as it is fully spin polarized, while the IKS is not.The Zeeman energy cost of both CHF and IKS can be computed via g s µ B tr{ Qs z }, where g s = 2, µ B is the Bohr magneton, and s z is the z-component of the spin.We find that numerically switching off the Zeeman splitting increases the critical field at which the CHF-IKS transition occurs.Interestingly, the same argument does not hold along (−3, −1) and (−1, −3), as IKS does not suffer from extra Zeeman energy cost compared to CHF.It would appear that the relative strength of the Coulomb interaction and the non-interacting bandwidth of the relevant magnetic subbands (red colored states in Fig. 3a) also plays an important role.We find numerical evidence that decreasing the Coulomb interaction (e.g., changing the relative dielectric constant from 15 to 25) makes IKS more stable and increases the critical field, conversely increasing Coulomb interaction favors a CHF state.Finally, there is also evidence that increasing the strength of the uniaxial heterostrain (e.g., from 0.2% to 0.3%) favors IKS over CHF, and increases the critical field.
In a recent experiment, the spin polarizations of the main CCIs near the magic angle are identified by edge conductance measurements, with (−3, −1) and (2, 2) being spin polarized, and (−2, −2) being spin unpolarized [17].Our theory naturally recovers the spin polarization along (−3, −1).While CHF is fully spin polarized, as mentioned, the IKS is not.Therefore, spin polarization along ±(2, 2) could be used to differentiate between CHF and IKS states as shown in Fig. 3(f,g).The experimental results therefore suggest an IKS state along (−2, −2) and CHF along (2, 2), although to fully explain them, more careful analysis of the edge modes and modeling beyond the BM Hamiltonian [56] which captures the particlehole asymmetry are necessary.

D. 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 ex-pectation value of σ z τ z for the occupied electronic states, where σ z and τ z are Pauli matrices acting in the sublattice and valley subspace respectively [29,30,47,49].The solid lines in Fig. 4(a) and (b) show calculated ⟨σ z τ z ⟩ per electron as a function of twist angle in the presence of heterostrain along (0, −4), (−1, −3), (−2, −2), and (−3, −1), for ϕ/ϕ 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 quantitatively different from both the HSFs and the sCIs in the presence of heterostrain, but closer to HSFs.Intuitively, the increased ⟨σ z τ z ⟩ of CHFs compared to HSFs may originate in the short-ranged part of the Coulomb repulsion disfavoring two electrons sitting close to each other in real space.
On the other hand, in the absence of heterostrain, as shown in Fig. 4(c) and (d), the CHFs along (−3, −1) smoothly cross over into the sCIs upon lowering the twist angle.For our model parameters, there is a collapse of the ⟨σ z τ z ⟩ along (−2, −2) (−1, −3) and (0, −4) at lower twist angles, when the CHFs become energetically less favorable than populating the Landau quantized excitation spectra of the IVC states [33,40,47,57].As further demonstrated in Fig. S19, this transition depends sensitively on model parameters, and can be pushed toward lower ϕ/ϕ 0 (e.g., by moving toward the chiral limit, see SI and Ref. [58]).Refs.[11,21] report that the (2, 2) persist down to ϕ/ϕ 0 ∼ 1/25, 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.
In our earlier work [47] we computed the single particle excitations of IVC insulators in the strong coupling limit at B ̸ = 0, which have been demonstrated to be the ground states at the CNP and n/n s = −2.At the CNP and at low B we find two-fold degenerate Landau levels (LLs) 0, ±2, ±4, . . . .At n/n s = ±2 we find instead 0, ±1, ±2, . . . .These have been corroborated by other works [57,59].Although the results in Ref. [47] assumed adding a single electron or a single hole, they are expected to hold at an asymptotically low B even at a small but finite density along (0, t) or (±2, t).This is because the energy difference between competing many-body states due to addition of a small density of carriers necessary to fill the excited LLs is expected to be proportional to the number of flux quanta, while the energy difference between competing states at B = 0 is extensive and thus proportional to the total particle number.Therefore, a finite critical B-field would be necessary to tip the energy balance in favor of a state such as the sCI, distinct from the one obtained by a naive rigid filling of the LLs.Closer examination of Fig. 2(f) indeed demonstrates this.At n/n s = −2, we find a prominent (−2, −1) which corresponds to emptying one energetically well-separated LL from the spectra of the (−2, 0) IVC state, consistent with our earlier studies.Along (−2, −2) and at ϕ/ϕ 0 < 1/7, a weaker gapped state is observed where two LLs are emptied.However at ϕ/ϕ 0 > 1/7 the sCI is stabilized via a first order phase transition.At n/n s = 0, we similarly find a prominent (0, −2) corresponding to emptying a two-fold degenerate LL of the (0, 0) state (with (0, −1) being QHFM of (0, −2)).Interestingly, gapped states along (0, −4) -expected based on the results in Ref. [47]are not observed down to ϕ/ϕ 0 = 1/12.Given that this corresponds to an electron filling of n/n s = −1/3, we attribute the absence of (0, −4) to band renormalization effects at finite electron densities not captured in Ref. [47].Fig. 2(f) also shows gapped states emanating from n/n s = −3 and −1 that can be characterized either as sCIs [(−3, −1), (−1, ±1)] or via (de)population of the sCIs' Landau quantized excitation spectra [(−3, 0), (−1, −2), (−1, 0)].We refer interested readers to details presented in Fig. S18.We note that in the strong coupling limit there are near degeneracies between competing states.For example, along (−1, −3), a Chern −3 sCI can be found as a metastable state, with a Hartree-Fock energy ≲ 0.01meV higher than the nearly compressible state plotted in Fig. 2(f).

E. Main CCIs in the absence of heterostrain
Though most of the experiments so far are believed to be strongly influenced by uniaxial heterostrain, there are a few experiments which appear to be in the ultra low heterostrain regime [12,18].In addition to the main CCIs along (0, ±4), ±(1, 3), ±(2, 2) and ±(3, 1), these experiments show that the CNP develops a gap at B = 0 at low temperatures without any apparent hBN alignment.This is in contrast to a (gapless) semimetal found at a moderate amount of uniaxial heterostrain (e.g.0.2%) [12,38].As discussed in the previous section, in the strong coupling limit, there are prominent gapped states along (−2, −1), while (−1, −3) and (0, −4) are gapless at reasonable magnetic flux ratios.In contrast, (−2, −1) does not appear to be a prominent gapped state in experiments without hBN alignment, and the main CCIs are ubiquitous.We therefore conclude that the above low heterostrain experiments [12,18] cannot be in the strong coupling regime with a negligible non-interacting bandwidth.
We investigate the impact of a finite non-interacting bandwidth by systematically studying the finite B phase diagram for a range of twist angles from 1.38 • down to the magic angle of 1.05 • , in an analogous fashion compared to studies with 0.2% heterostrain.The single-electron excitation gaps at different twist angles, electron densities, and B field are presented in Fig. 2. Through separate B = 0 Hartree-Fock calculations we demonstrate that at CNP the ground state is a gapped IVC state for all twist angles, albeit with an IVC order parameter localized to the vicinity of the K points of the moiré Brillouin zone at larger twist angles.At the largest twist angle 1.38 • with the largest non-interacting bandwidth, the gap structures are very similar to that in the presence of heterostrain.gapped state consistent with previous works [47,57,59].However, gapped state along (0, −4) is missing, and we attribute it to a finite B and finite electron density regime where populating the rigid excitation spectra of n/n s = 0 IVC insulator is no longer energetically favorable.At the twist angle 1.2 • , as shown in Fig. 2(e), the main CCIs along (−3, −1), (−2, −2) and (−1, −3) remain robust down to the lowest magnetic flux ratio ϕ/ϕ 0 = 1/12 studied.Simultaneously there are strong IVC states along (0, 0) and (−2, 0).We also do not find a strong gapped state along (−2, −1).Moreover, the main gapped states emanating from the CNP are along (0, 0), (0, −2) and (0, −4), consistent with the mentioned experiments.To better understand these states we present their detailed Hartree-Fock spectra and representative density matrices in Fig. 5.At n/n s = 0 and along (0, 0), the IVC state hybridizes {K, ↑ (↓)} with {K ′ , ↓ (↑)}, creating bonding and antibonding states.The occupied bonding states have net Chern number 0, resulting in the gapped state along (0, 0).In contrast, the gapped state along (0, −4) is a CHF, where the occupied group of states in each valley/spin sector contributes a Chern number −1.Along (0, −2), we identify the gapped state as a CHF in the spin ↑ sector which accounts for the net Chern number −2, and an IVC in the spin ↓ sector having zero Chern number.Analogously at n/n s = −2, the gapped state along (−2, 0) forms intervalley coherence between {K, ↓} and {K ′ , ↓}, whereas the gapped state along (−2, −2) is a CHF.Along (−2, −1), we find that the state is described by adding one Landau level worth of electrons on top of a CHF ground state along (−2, −2), rather than adding one Landau level worth of holes to an IVC state along (−2, 0).The latter state is metastable at 1.2 • but is the Hartree-Fock ground state at 1.05 • , as illustrated in Figs.2(e,f).Our results at twist angle 1.2 • therefore capture the most salient features of the experiments in Refs.[12,18].
It is also interesting to note that at finite B, the (0, 0) IVC state mixes opposite spin species.Within our numerics we can find an IVC state which mixes the same spin species, which has a slightly higher energy.Nevertheless, we find that the energetic difference between these two IVC states grows as B increases.For example, at ϕ/ϕ 0 = 1/8 the energy difference between these two states is ∼ 0.02meV per moiré unit cell, but grows to ∼ 0.4meV at ϕ/ϕ 0 = 5/11.Although not strictly in strong coupling, they are similar to the zCIs discussed in the strong coupling limit.The spectra of CHFs smoothly extrapolate to these states as B decreases even though some are metastable at B = 0.At 1.28 • , the B = 0 states are stable ground states within our Hartree-Fock calculation, preserve C2zT and are compressible, but displaying sizable exchange splittings (∼ 12meV) reflected in the relative energy shift of the Dirac cones of different valley and spin flavors.The spectra of CHFs also extrapolate to the spectra of the B = 0 states, albeit with some differences such as the exact locations of the zLLs of the exchange split Dirac cones.In (a) and (d), the {K, ↓} states (red) are hidden behind {K, ↑} (blue) at B = 0, and hidden behind {K ′ , ↓} at B ̸ = 0.In (b) and (e) at B ̸ = 0, the spectrum of {K, ↓} (red) is hidden behind {K ′ , ↓} (green), while that of {K, ↑} (blue) behind {K ′ , ↑} (orange).In (c) and (f) and finite B, {K, ↑} (blue) states are hidden behind {K ′ , ↑} (orange).
As we showed in the previous section, at 1.05 • , at low B, and in the absence of heterostrain, this can be understood from the Landau quantization [47,57] of the strong coupling single-particle excitation spectra of the flavor symmetry breaking ground states at B = 0 [59,60].At integer fillings on the hole-doped side of the CNP, a doped hole (i.e.moving away from CNP) has a light mass, leading to well separated LLs (large cyclotron frequency) in relatively low B. In contrast, a doped electron (i.e.moving toward the CNP) has a heavy mass, and the LLs are much more densely spaced at a comparable B. While this explanation holds at low B for 1.05 • twist angle without heterostrain, Figs. 1 and 2 show that the main Landau fans point away from the CNP even at larger twist angles or when we include moderate heterostrain.Under such conditions the strong coupling limit is not reached as can be seen by the absence of the correlated insulators at B = 0 or the presence of IKS states which do not saturate the sublattice-valley polarization (Fig. 4a,b).
In order to understand why this result persists away from the strong coupling, we start by focusing on 1.28 • in the absence of heterostrain whose phase diagram is presented in Fig. 2(c).In Fig. 6(a-c), the left columns show the Hartree-Fock spectra of the respective ground states at B = 0 and n/n s = −1, −2, −3.These ground states preserve the valley U v (1), spin U s (1), C 3z and C 2z T symmetries.As a result, the Dirac points (DPs) of any given valley/spin flavor are protected and located at the K point of the moiré Brillouin zone.Crucially however, for all three fillings there is an exchange splitting (∼ 12meV) between different flavors, as reflected in the relative energetic shift between the DPs of different valley/spin character.These ground states are all compressible, where the DPs are shifted above the Fermi energy, with the charge being compensated by the Fermi pockets from the exchange split bands.Moreover, there is a band renormalization effect for all flavors, reflected in the narrowing (sharpening) of the bands below (above) the DPs.The band renormalization becomes stronger as the filling is tuned to the band edge (e.g., n/n s = −3).Such band flattening effect has been discussed in the literature as a Hartree effect [24] (our calculation here also shows exchange splitting due to Fock terms).The second and third columns show the B ̸ = 0 Hartree-Fock spectra of the CHFs which share the same valley and spin quantum numbers as the B = 0 states, and color-coded in the same manner.It is evident that spectra of the CHFs can be smoothly extrapolated to the B = 0 Hartree-Fock dispersions, e.g. by matching the energies of the DP zLLs to DPs at B = 0.This demonstrates that the CHFs emerge from the respective compressible ground states at B = 0.
Generally, the heavy-light dichotomy at high twist angles (such as 1.28 • ) can be understood by examining the B = 0 dispersions.We address it using the CHF (or lack thereof) along (−3, ±1) as an example.As shown in Fig. 6(c), the CHF is formed by populating the Chern −1 group of magnetic subbands below the DP zLLs of a given valley/spin flavor (in this case {K, ↓} shown in the second column).However, Landau quantizing the B = 0 dispersions show that CHF is not an energetically favorable state at infinitesimally weak B field, as the DP LLs of the K, ↓ are buried inside the dense LLs associated with the Fermi pockets (i.e., heavier electronic states due to band flattening) of the other three flavors.Only at larger B, can the Chern −1 magnetic subband group be separated from the rest of the spectra due to the wide LL spacings of a linearly dispersing Dirac cone.This ties the finite B CHF along (−3, −1) directly to the Landau quantizations of the B = 0 ground state.Conversely, the zLLs of the DPs cannot be separated from the dense LLs, therefore the state along (−3, 1) is nearly compressible (unless intercepted by gapped states along another (s, t)).The B = 0 and finite B correspondence can also explain why the critical field for the onset of CHFs along (−3, −1) occurs at a lower B than (−2, −2), which in turn occurs at a lower B than (−1, −3).As n/n s increases from −3 to −1, the flavor degeneracy of the DPs (marked by green or red crossings) increases while the flavor degeneracy of the Fermi pockets -originating from the exchange split bands-decreases.As a result, the size of the Fermi pockets grows and the bottom of these pockets sinks deeper below the DPs.A higher B is therefore necessary to separate the aforementioned Chern -1 group of magnetic subbands per flavor from the rest of the spectra.
Next, we focus on 1.20 • in the absence of heterostrain, where interaction effects are stronger than at 1.28 • , with the CHFs along (−3, −1), (−2, −2) and (−1, −3) all persisting down to the lowest magnetic flux ratio studied in this work.The finite B phase diagram is presented in Fig. 2(e), and the comparisons to B = 0 states are presented in Fig. 6(d-f).These B = 0 states are obtained by enforcing the valley U v (1) and spin U s (1) symmetries, and share the same valley/spin quantum numbers as the respective CHFs at finite B. While the state at n/n s = −3 is the Hartree-Fock ground state, states at n/n s = −2, −1 are metastable.Comparing the B = 0 and finite B spectra we conclude that the CHFs originate from the Landau quantization of these B = 0 bands.A notable difference compared to 1.28 • is that the C 2z T symmetry is spontaneously broken at all n/n s = −3, −2, −1, making them (meta-)stable and gapped zCIs (with Chern numbers ±1, ±2, ±3 respectively), analogous to those discussed in the strong coupling limit [32,33].Given that the B = 0 states at n/n s = −2, −1 are metastable, at very low B, these CHFs must lose to populating LLs of the excitation spectra of the respective true ground states (with IVC order, see Fig. S20), likely via a first order phase transition.This is in contrast to higher twist angles (e.g., 1.28 • ) where there are no competing zCIs states nearby.
We finally address the heavy-light dichotomy for the calculations performed with 0.2% of uniaxial heterostrain.Due to competing Chern 0 IKS states (gapped or compressible) which persist to higher twist angles, it is challenging to find the B = 0 metastable states from which the finite B CCIs descend, as we can in the absence of heterostrain.We therefore postpone tying the finite B and B = 0 physics to a future work.Here we instead provide an intuitive picture by examining Hartree-Fock spectra directly at finite B. We use CCIs (or lack thereof) along (−3, ±1) at 1.28 • twist angle as an example.Along (−3, −1), both CHF and IKS predominantly involve the lower Chern -1 group of the noninteracting magnetic subbands, separated from the zLLs and the upper Chern -1 group by a gap (see Fig. 3(a) for 1.2 • ; at 1.28 • the structure of the magnetic subbands is qualitatively similar).The bandwidth of the lower Chern -1 group narrows upon increasing B, creating a favorable condition for symmetry breaking states (either CHF or IKS) driven by Coulomb interactions [15].Along (−3, 1), to get a similar flavor symmetry breaking state, one would involve both the lower Chern -1 group and the two zLLs.However, as B increases, this Chern +1 composite of magnetic subbands has an increasing bandwidth, making it energetically more costly to break the valley/spin flavor symmetry.Conversely, it is more energetically favorable to populate two extra LLs from the nearby, exchange split, group of states on top of the (−3, −1) state.Due to band flattening effect discussed above, an added electron to the bottom of the spectra of a given valley/spin flavor on the hole side of the CNP is heavy, and reflected as nearly compressible along (−3, 1).This is further supported numerically from Fig. S4, where the spectra above the gap at (−3, −1) (gray colored) do not have well-separated LLs.The finite B perspective presented here complements the perspective of tying finite B to B = 0 (meta-)stable states, and relies on the same Hartree-Fock effects of band (magnetic subband) renormalizations and exchange splitting.
At low twist angles in the presence of heterostrain, such as 1.20 • and 1.05 • shown in Figs.1(e) and (f), the gapped IKS states along (−3, −1) (−2, −2) and (−1, −3) remain robust down to the lowest magnetic flux ratio studied in this work, while gapped Chern 0 IKS states along (−3, 0) (both (−3, 0) and (−2, 0) at 1.05 • ) are also found.Moreover, we find that at 1.05 • , the IKS wavevector of these two classes of IKS states are very different (≈ 1  2 g 1 for Chern-0 IKS and ≈ 1 q g 2 for IKS carrying a Chern number).Given the earlier discussions of competing zCI and IVC states in the absence of heterostrain and at lower twist angles, it is tempting to conjecture that the gapped IKS states along (−3, −1), (−2, −2) and (−1, −3) may descend from B = 0 "topological IKS" (tIKS) states which break the C 2z T symmetry, i.e. a distinct IKS state from those reported in the literature [38].Should they exist, the tIKS state must be metastable at B = 0 and energetically unfavorable compared to the IKS state that preserves C 2z T .We leave more elaborate analysis of such a conjecture to future studies.

G. Additional CCIs
Besides the aforementioned CCIs, we also find additional correlated insulating states in the phase diagram both with and without heterostrain, as shown in Fig. 1 and Fig. 2 respectively.
In the presence of heterostrain, the most prominent states emanate from the CNP, and are identified as quantum Hall ferromagnetic states (QHFMs) within the zLLs of the energetically split Dirac cones [61][62][63] (see Fig. S7).QHFMs emanating from the band minimum (n/n s = −4) are well developed at 1.38 • but become progressively weaker as the angle decreases (see Fig. S6).Moreover, at higher ϕ/ϕ 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 ϕ/ϕ 0 = 1/2, as demonstrated in the non-interacting Hofstadter spectra in Fig. S1 and representative density matrices in Fig. S8.At 1.05 • , 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 L 2 direction, such that the density matrix is invariant under ( tL2 ) m .Their respective density matrices are shown in Fig. S9.We use the electron occupation number of the lower zLL (see e.g.Fig. 3(a)) in valley K and for spin ↓ to illustrate the striped states.We define it as n 0 (k), and show its momentum dependence in the magnetic Brillouin zone at ϕ/ϕ 0 = 1/6 for (−2/3, −3), (−1/2, −3), and (−3/2, −2) in Figs.7 (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.While some of the fractional s CCIs show intervalley coherence and others show valley/spin polarization, the energy difference between these two kinds of states are ∼ 0.05meV per moiré unit cell.Within our numerical accuracy, we cannot say with certainty if either will be observed in future experimental works.Furthermore, we are also limited to probe striped states only along L 2 direction, but not striped states along L 1 or checkerboard states modulating along both directions.Nevertheless, it is clear that states with broken magnetic translation symmetries (being striped or checkerboard states) are more energetically favorable than states that preserve them (or IKS states that preserve modified magnetic translation symmetries).More careful studies of fractional s states are beyond the scope of the present work, and will be left for the future.Ref. [14] observed fractional s CCIs in hBN aligned TBG devices, and presented a qualitative argument based on translation symmetry breaking phases at B = 0. Here, our results demonstrate that they can be stablized purely due to interactions and without hBN alignment.
In the absence of heterostrain and large twist angles, we identify gapped states emanating from the band minimum as QHFMs similar to the strained case (see Fig. S13).However the most prominent gapped states emanating from the CNP are no longer QHFMs even at the largest twist angle studied.As alluded to in Sec.III E, the CNP develops a small but finite correlation gap even at 1.38 • , supported by separate B = 0 Hartree-Fock studies.However the energy minimum of electron-like excitations (or maximum of hole-like excitations) remains at the hexagon corners of the moiré Brillouin zone.As the twist angle decreases, QHFMs emanating from the band minimum fade away similar to the cases in the presence of heterostrain, the size of the IVC gap at CNP grows and Coulomb interaction hybridizes electronic states further from the vicinity of the Dirac cones of the noninteracting band.This eventually leads to "inverted" excitation spectra at the magic angle 1.05 • , with the energy minimum of electron-like excitations shifted to the Γ point of the moiré Brillouin zone.Additionally, as in the strained case, fractional s CCIs are also found at the magic angle, but persist to higher twist angles compared to strained case (see Fig. S16).

IV. CONCLUSION
In summary, by performing a comprehensive selfconsistent 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 correlated Chern insulators are stablized at higher magnetic fields, and correspond to valley and spin polarizations of the interaction-renormalized magnetic subbands that we dub correlated Hofstadter ferromagnets (CHFs).Upon lowering magnetic field, the CHFs become energetically less favored, losing to gapped states with intervalley coherence.In the 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, ±4), ±(1, 3), ±(2, 2), ±(3, 1) sequence, notably the striped states at fractional s.Given that our calculations have direct access to the interaction renormalized singleelectron excitation spectra at a given filling and magnetic field (see Fig. S10), comparisons with experiments such as STM can be made to facilitate the characterization of the panoply of correlated insulating states.
Supplemental Material for "Theory of correlated Chern insulators in twisted bilayer graphene"

I. TBG WITH HETEROSTRAIN AT ZERO MAGNETIC FIELD
We consider the limit of small twist angle θ and small deformations of moiré superlattice due to uniaxial heterostrain, such that the moiré reciprocal lattice vectors are given by: g i=1,2 = E T G i=1,2 [37], where G i are reciprocal lattice vectors of the untwisted and undeformed monolayer graphene, and where Here {ε, φ} parameterize the uniaxial heterostrain strength and orientation, R φ is the two-dimensional rotation matrix, and ν ≈ 0.16 is the Poisson ratio.The moiré lattice vectors L i=1,2 are uniquely defined through the relation We consider the strained Bistritzer-MacDonald (BM) Hamiltonian discussed in our recent work [39], with the continuum Hamiltonian given as: where η = ±1 labels K (K ′ ) valleys of monolayer graphene, l = t, b labels the top (bottom) of the two graphene layers.The interlayer Hamiltonian is given by: where ψ η,l (r) ≡ (ψ η,l,A (r), ψ η,l,B (r)) T is a spinor in the sublattice basis for a given valley and layer.We have suppressed the spin index for simplicity.q j=1,2,3 are the three nearest neighbor bonds of the reciprocal honeycomb lattice, and Here (σ 0 , σ x , σ 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. [58] corresponds to w 0 /w 1 = 0.The intra-layer Hamiltonian is given by: where the first term is the deformation potential that couples to the electron density.A η,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 2a η(ϵ l,xx − ϵ l,yy , −2ϵ l,xy ).We use α ≈ −4.1 eV and β ≈ 3.14.We further fix ℏv F /a = 2.68eV, w 1 = 110meV and w 0 /w 1 = 0.7 in our calculations, and also set constants ℏ = 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 A. 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 K and K ′ are given by (here for negative charged electrons, the minimal coupling is p + eA where e is positive): Here l = 1, 2 is the layer index, 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 A ≡ Bxŷ, where the global x − y coordinate system is defined such that L 2 is along the +ŷ direction.When TBG is subject to heterstrain, this amounts to a small angle rotation θ 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: √ eB is the magnetic length.In valley K, the particle-hole symmetric LL eigenstates are given as: where is the energy of the Dirac Hamiltonian, labeled by n = 1, 2, . . ., and γ = ±1 corresponds to positive and negative energy solutions.We have defined L 2 ≡ |L 2 |.In addition, there is an anomalous zero energy state given by which lives on the B sublattice.ϕ n (x) is the eigenfunction of a † a, and is given to be: where H n (x) is the Hermite polynomial.The shift in the position for a given wavevector k 2 g 2 is given by: In valley K ′ , the Landau level wavefunctions are solved in a similar manner, and the LL eigenbasis is given by: ℓ , and eigenstate for the zeroth LL: For notational convenience we define a ket state Φ (η) nγl (k 2 ) such that: The ket state has the following real space properties: where the definitions of Φ (η) nγl (x + k2,l ℓ 2 ) are inferred from Eqs. (S10), (S11), (S14) and (S15).
The notation [b] a represents b modulo a, with a > 0. We have defined: We proceed to work out the 1D integration by noting that: x and as a result: where Âηη ′ nγl,n ′ γ ′ l ′ can be expressed in terms of associated Laguerre polynomials [47], and: Next we consider the implications of the δ function constraints.Firstly we have: Secondly we have: This shows that above matrix elements are non-vanishing only if: Finally, we have the following expression for the matrix elements: )s e i 2 s 2 q ϕ •L1 e −isK η ′ l ′ ,x L1,x e iqx k2,ηl ℓ 2 e i 2 qx qyℓ D. Matrix elements of strained BM Hamiltonian

Interlayer terms
The interlayer terms are a special case of the generic A q discussed in the previous section.For bottom to top layer tunneling we have (see Eq. (S24)): )s e i s 2 2 q ϕ •L1 e −isKηt,xL1,x e iqx k2,ηl ℓ 2 e i 2 qx qyℓ where q j = 0, g 2 , g 1 + g 2 .For notational convenience we have also defined q j ≡ q 1 g 1 + q 2 g 2 where q 1,2 are integers (with implicit dependence on j).Plugging these into the δ-function constraints we have: We can split the momentum p 2 , k 2 into strips of [0, 1/q), . . .[(p − 1)/q, p/q), by redefining: It is clear then that (k 1 , k 2 ) ∈ [0, 1) × [0, 1/q) are good quantum numbers under moire periodic potential, and that:

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 ) ∈ [0, 1) × [0, 1/q), and an additional 2q quantum numbers per valley and spin.Generally we can denote them as: 3. Degeneracy of magnetic subbands generated by tL 2 The non-interacting Hamiltonian H non-int in valley η acting on the narrow band eigenstates is: We apply the magnetic translation tL2 to these energy eigenstates.Note that they are eigenstates of tL1 and tq L2 ; the reason why tL2 is applied q times is due to non-commuting nature of magnetic translation operators tL1 and tL2 .We see that: Furthermore since [H non-int , tL2 ] = 0, we have: namely that tL2 Ψ (η) a (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 → k 1 + 1/q, and that all the distinct energy solutions can be found in the domain of (k 1 , k 2 ) ∈ [0, 1/q) × [0, 1/q).Furthermore, the tL2 translation gives us one way of gauge fixing between k 1 values in different strips of [0, 1/q).Specifically, we can define: as the gauge-fixed eigenstate wavefunction at ([k 1 + p/q] 1 , k 2 ).This amounts to the following definition: nγlr,a (k 1 , k 2 ), r 1 = 0, . . .q − 1. (S41) E. Matrix elements of Aq = Âe −iq•r with respect to narrow band eigenstates This can be computed as follows: )s e i s 2 2 q ϕ •L1 e −isK η ′ l ′ ,x L1,x e iqx k2,ηl ℓ 2 e i 2 qx qyℓ As a result, the Coulomb piece of the Hamiltonian can be written as: ) where G = mg 1 + ng 2 is the moiré reciprocal lattice vector.We note first on the constraint of the δ-functions on the allowed values of momentum.Specifically we have: Since p 1 − k 1 ∈ (−1, 1) and p 2 − k 2 ∈ (−1/q, 1/q), the following constraint applies: A. Product state and one-particle density matrix At a given filling, we consider a ground state |Ω⟩ to be a product state, given by partial fillings of the narrow band states |Ψ ηs a (k)⟩.|Ω⟩ can always be expressed in terms of a product of Bogoliubov quasiparticle creation operators acting on vacuum, where number of γ † i,k appearing in the product is constrained by the filling.Generally the Bogoliubov quasiparticles are related by the original fermions d ηsa,k via a unitary transformation: The Hartree-Fock procedure is to find the optimal set of unitary rotations {U ηsa,i (k)} which minimizes the total energy.Alternatively we can define a one-particle density matrix (note here we used a different definition compared to the main text, and for clairity we also changed the notation from Q in the main text to P here): It is related to {U ηsa,i (k)} via: Note that tL2 takes k to k + ϕ/ϕ 0 g 1 , and as a result, a difference of the density matrix at P (k) and P (k + ϕ/ϕ 0 g 1 ) allows us to probe magnetic translation symmetry breaking of tL2 .

B. Mean field energy
Here we calculate the mean field energy E Ω with respect to the density matrix P defined above.For simplicity and generality we derive E Ω with the following Hamiltonian instead of the notation used in Eq.S45: Firstly the contribution to the energy from the kinetic term is given as: Next we address the Coulomb piece E Ω .We define , and: The last term can be evaluate as follows: Collecting all the terms we get: Therefore, the total mean field energy is given as: where c Ω is a constant number independent of the density matrix.

C. Energy optimization and Hartree-Fock Hamiltonian
Eq. S61 needs to be minimized subject to the constraint that P + 1 2 I = P + 1 2 I 2 is a projector onto the narrow band Hilbert space.This is equivalent to minimizing the following "free energy" [38], Here { X} are the Lagrange multipliers.If Pc is the desired density matrix, then for any small perturbation P1 satisfying tr{ P1 } = 0, the coefficient to linear order term in P1 must vanish.We obtain: where on the last line we defined the Hartree Fock mean field Hamiltonian as: Requiring Eq.S63 to vanish leads to the following mean field equation: Observing the mean field equation, we first see that a solution must satisfy ĤMF ( Pc ) P T c + 1 2 This is equivalent to the following commutation relations: ĤMF ( Pc ), Pc = 0, X, Pc = 0. (S67) Finally, with respect to the optimized Pc , the total energy is given by: Our B-SCHF results are obtained by iteratively solving the self-consistency equations (Eq.S67).In practice, we run both random initializations and educated guesses for the density matrix P , such as a strong coupling Chern insulating state, flavor polarized state, intervalley coherent state, etc.Typically, every point on the phase diagram in Fig. 1 and Fig. 2 of the main text is a result of ∼ 20 initializations of the density matrix.
For each magnetic flux ratio ϕ/ϕ 0 , we study momentum space mesh consisting of n q q points along g 1 direction in the range of [0, 1), and n q points along g 2 direction in the range of [0, 1/q), where n q is chosen as the maximum integer such that n q q ≤ 14.Our calculations presented in the main text correspond to 23 magnetic flux ratios between 1/12 and 1/2, corresponding to all {p, q} satisfying q ≤ 12 and p < q.

D. Probing IKS like states
In our B-SCHF calculations, the density matrix defined in Eq.S54 can probe states with intervalley coherence, where QKs,K ′ s ′ a,b (k) becomes nonzero.Numerically, we also look for the finite B analog of the incommensurate Kekulé spiral ordered states (IKS) discussed in Ref. [38].In zero magnetic field, IKS has intervalley coherence, and breaks the moiré translation symmetry TLi , i = 1, 2. However the state preserves the discrete moiré translation symmetry followed by U(1) valley gauge transformations: where Q IKS is the IKS wavevector, τz is the Pauli matrix acting on the valley space (τ z |K⟩ = |K⟩ , τz |K ′ ⟩ = − |K ′ ⟩).
In momentum space, this corresponds to a non-vanishing density matrix element: In finite B, the discrete moiré translations are replaced by the non-commuting magnetic translations tLi .We seek analogs of the zero field IKS states by looking for intervalley coherent states which preserves the magnetic translation symmetry followed by U(1) valley gauge transformations: The many body wavefunction for such state can be generally written as follows: Here, n is the band index of the Bogoliubov quasiparticles, k = k 1 g 1 + k 2 g 2 , with k 1 , k 2 ∈ [0, 1/q).The constrained ( ′ ) product is over all occupied Bogoliubov quasiparticle states.We have split the magnetic Brillouin zone in to q strips along the g 1 axis, and label each strip by r = 0, . . .q − 1. q 0 = q 0,1 g 1 + q 0,2 g 2 is the relative wavevector between electrons in opposite valleys.{α, β} are coefficients that we shall constrain below to make above an IKS state respecting the tLi symmetries.Note that from the gauge fixing procedure discussed previously, the following relation applies, Note here that tL2 boosts the magnetic momentum k by the magnetic flux ratio ϕ/ϕ 0 along the g 1 axis.
We begin by noting that |Ψ IKS ⟩ is an eigenstate of e −i 1 2 q0•L1 τz tL1 , as: where θ nk is independent on the spin and magnetic subband quantum numbers, and φ is independent of all quantum numbers, we have: The IKS wavevector is given by: Note due to the cyclic nature of the state under tL2 , φ 2π = n ϕ ϕ0 where n ∈ Z. IKS-like states can be studied using the density matrix as defined in Eq.S54.On the technical side, for any q 0 , we can displace the momentum space mesh of valley K ′ from valley K by q 0 , before carrying out the (magnetic) translationally invariant B-SCHF calculations.This is how we probed these states in the paper, by optimizing the ground state energy with respect to q 0 .We look at the constraints imposed on the density matrix (Eq.S54) by an IKS ground state.Without confusion, we simplify the notations by redefining fermionic operators in valley K ′ as d K ′ sb,k+q0 ≡ d K ′ sb,k in case q 0 ̸ = 0, and keep the definition of fermionic operators the same in valley K.We get: where η, η ′ = +(−)1 denoting valley K (K ′ ).This property of the density matrix for an IKS-like state can be used to reduce the calculation of the Hartree-Fock energy (Eq.S61) by a factor of q.

E. Stripe states
Unlike IKS states, stripe states break the magnetic translation symmetries, and lead to density/spin modulations on the moiré scale.In our code, we can probe stripe states commensurate with the enlarged unit cell.Specifically in the Landau gauge, the unit cell is chosen to be 1 moiré unit cell along L 1 direction and q moiré unit cells along L 2 direction.This allows us to probe stripe states with periods divisible by q.For example, at ϕ/ϕ 0 = 1/8, we could probe period 2, 4 stripes along L 2 , but not period 3, 6 and so on.We could not probe checkerboard states, nor stripe states with modulations along L 1 .This is important to keep in mind when examining the phase diagrams in the paper.

IV. EXTENDED FIGURES A. Extended results in the presence of heterostrain
Here we present B-SCHF results obtained for various twist angles in the presence of heterostrain.

FIG. 6 .
FIG. 6. Comparative study of B = 0 Hartree-Fock spectra at integer fillings n/ns = −1, −2, −3 of states that preserve valley Uv(1) and spin Us(1) symmetries, and the B ̸ = 0 Hartree-Fock ground state spectra of CHFs along (−1, −3), (−2, −2) and (−3, −1).The B = 0 states share the same valley and spin quantum numbers with the CHFs.Panels (a-c) are at 1.20 • twist angle and panels (d-f) are at 1.28 • , all in the absence of heterostrain.At 1.20 • the B = 0 states at n/ns = −1, −2 are metastable, while the state at n/ns = −3 is a stable ground state.They all break the C2zT symmetry and are incompressible Chern states with Chern numbers −3, −2, −1 respectively.Although not strictly in strong coupling, they are similar to the zCIs discussed in the strong coupling limit.The spectra of CHFs smoothly extrapolate to these states as B decreases even though some are metastable at B = 0.At 1.28 • , the B = 0 states are stable ground states within our Hartree-Fock calculation, preserve C2zT and are compressible, but displaying sizable exchange splittings (∼ 12meV) reflected in the relative energy shift of the Dirac cones of different valley and spin flavors.The spectra of CHFs also extrapolate to the spectra of the B = 0 states, albeit with some differences such as the exact locations of the zLLs of the exchange split Dirac cones.In (a) and (d), the {K, ↓} states (red) are hidden behind {K, ↑} (blue) at B = 0, and hidden behind {K ′ , ↓} at B ̸ = 0.In (b) and (e) at B ̸ = 0, the spectrum of {K, ↓} (red) is hidden behind {K ′ , ↓} (green), while that of {K, ↑} (blue) behind {K ′ , ↑} (orange).In (c) and (f) and finite B, {K, ↑} (blue) states are hidden behind {K ′ , ↑} (orange).

FIG. 7 .
FIG. 7. Occupation number n0(k) of the lower zLL of the non-interacting spectra in valley K and for spin ↓. Results are obtained at 1.05 • for correlated Chern inslating states with fractional s, marked by the dashed lines in Fig. 1(f).