Kekul\'e spiral order at all nonzero integer fillings in twisted bilayer graphene

We study magic angle graphene in the presence of both strain and particle-hole symmetry breaking due to non-local inter-layer tunneling. We perform a self-consistent Hartree-Fock study that incorporates these effects alongside realistic interaction and substrate potentials, and explore a comprehensive set of competing orders including those that break translational symmetry at arbitrary wavevectors. We find that at all non-zero integer fillings very small strains, comparable to those measured in scanning tunneling experiments, stabilize a fundamentally new type of time-reversal symmetric and spatially non-uniform order. This order, which we dub the 'incommensurate Kekul\'e spiral' (IKS) order, spontaneously breaks both the emergent valley-charge conservation and moir\'e translation symmetries, but preserves a modified translation symmetry $\hat{T}'$ -- which simultaneously shifts the spatial coordinates and rotates the $U(1)$ angle which characterizes the spontaneous inter-valley coherence. We discuss the phenomenological and microscopic properties of this order. We argue that our findings are consistent with all experimental observations reported so far, suggesting a unified explanation of the global phase diagram in terms of the IKS order.


I. INTRODUCTION
The discovery of superconductivity proximate to correlated insulating behaviour in a variety of graphene moiré heterostructures [1][2][3][4] has triggered intensive efforts to explore the phase structure of these highly tunable two-dimensional materials. In the best-studied example, twisted bilayer graphene (TBG) tuned to a 'magic angle' of approximately 1°, the enhancement of correlations is associated to the formation of extremely narrow bands due to the reconstruction of the electronic dispersion by the moiré superlattice. As noted by Bistritzer and MacDonald (BM) [5], this can be elegantly captured within a continuum model [6][7][8] where Dirac cones contributed by isolated graphene layers are coupled by interlayer tunneling modulated at the moiré scale. The BM model serves as a starting point for most theoretical studies of TBG. On combining the degeneracies corresponding to spin and microscopic valley indices with the Dirac structure enforced by exact and approximate symmetries of the moiré band structure, the model reveals that the striking effects reported in experiments occur when the Fermi energy is tuned to lie within an octet of nearly-flat bands that straddle charge neutrality. The phase structure of TBG then turns on the question of how electron correlations and other perturbations such as strain and substrate effects lift the approximate degeneracy within this subspace to select between a variety of competing ground states.
At first sight, the phase diagram of TBG resembles that of the cuprate high-temperature superconductors, with electrostatic gating playing the role of chemical doping. This prompted initial attempts to model correlation effects within a single-band Hubbard model for electronic states localized to a triangular moiré superlattice. Although this approach has proven fruitful in studying moiré heterostructures of MoS 2 and other transition-metal dichalcogenides [9], its applicability to TBG is limited by the 'fragile topology' of the BM bands [10,11]. The latter requires that the simplest tight-binding model that faithfully captures the symmetries of TBG involves a pair of crystallographically-inequivalent Wannier orbitals centered on sites of a honeycomb lattice, but with their charge densities peaked on a triangular lattice formed by the centers of its hexagons [10,12,13]. A corollary of this Wannier representation is that it implies a higher degree of itineracy than can be captured via a minimal honeycomb lattice Hubbard model with on-site repulsion and nearest-neighbor hopping.
The utility of a Hubbard description is further challenged by the early experimental observation of a quantized anomalous Hall (QAH) resistance in TBG samples aligned with a hexagonal boron nitride (hBN) substrate, at an electron density corresponding to filling 7 of the 8 bands [14]. [In the convention where ν = 0 represents the filling at neutrality, this corresponds to ν = +3.] Since this occurs in the absence of an external magnetic field, it indicates the spontaneous breaking of time reversal symmetry (TRS). One explanation of this phenomenon invokes a compelling analogy to the other paradigmatic setting for strong correlations: the celebrated Landau levels (LLs) of an electron in a magnetic field. The TBG flat bands are endowed with non-trivial topology encoded in the winding of their Bloch functions across the moiré Brillouin zone (mBZ); the inclusion of substrate potential triggers the opening of gaps between the moiré Dirac points (by breaking one of the protecting symmetries), and assigns nonzero Chern numbers to the bands, making them topologically equivalent to LLs. The absence of explicit TRS breaking is reflected in the assignment of equal and opposite Chern numbers to different valleys, which are exchanged by TRS. At ν = 3, electrons in the two remaining unfilled bands spontaneously polarize into one of the two valleys. This allows them to min-2 imize their interaction energy by virtue of Pauli exclusion, leading to an insulator that spontaneously breaks TRS, with a quantized Hall resistivity of ρ xy = h/e 2 protected by the charge gap. Recently, a similar QAH state was also observed at ν = 1 [15]. While the formation of such orbital Chern insulator states can in principle be captured within a Hubbard description [16], its close parallels to quantum Hall ferromagnetism [17][18][19] (QHFM) has motivated a distinct perspective, where TBG is viewed as a generalized multicomponent quantum Hall system. This naturally explains both the observed QAH response as well the propensity for insulating states at commensurate filling, and motivates a sigma model description based on a hierarchy of perturbations around a 'hidden' limit with U (4) × U (4) symmetry [20][21][22]. The QHFM picture receives further experimental support by the observed stabilization of QAH insulators 1 with Chern numbers C = ±3, ±2, ±1 at ν = ±1, ±2, ±3 on applying a small out-of-plane magnetic field, even in the absence of substrate alignment [15,[23][24][25][26]. However, the TBG bands nevertheless retain features absent in LLs. For instance, their dispersion (though small) remains nonzero, and is enhanced when particle-hole symmetry breaking effects are incorporated -especially in the electron-doped regime -or upon inclusion of strain. Such effects are likely important in giving an accurate description of experimental samples. As a case in point, even at commensurate fillings some experiments 2 report gapless states or insulators with Chern numbers distinct from those of the noninteracting bands [15,27,28]. This suggests that departures from the flat band/QHFM limit are non-negligible and that the competition between itineracy and localization characteristic of Hubbard physics remains relevant to TBG.
Given its enticing position at the intersection of two dominant themes of strong correlations, it is natural to conjecture that orders that are 'natural' from both perspectives could be particularly robust candidate ground states in TBG. One example (and our focus here) is furnished by states with broken translation symmetry, which emerge in relatively well-understood limits of both the Hubbard and quantum Hall settings. The formation of charge and/or spin stripe order is believed to be a nearuniversal consequence of hole-doping the cuprates away from commensurate filling: while purely on-site Hubbard repulsion favors phase separation, the inevitably present further-neighbour interactions frustrate this in favor of spatially-ordered phases [29][30][31][32][33][34][35][36][37][38][39][40][41]. For similar energetic reasons a variety of stripe and bubble phases are 1 Strictly speaking, these are not anomalous since they only emerge in a finite magnetic field. However the small fields required (significantly smaller than ∼ 10 T corresponding to one flux per unit cell) suggest that these Chern states remain competitive at zero field. 2 Note that while Ref. [27] reported an anomalous Hall effect at ν = +3, the Hall conductance was not quantized, which is consistent with their samples hosting a gapless state at this filling.
known to be competitive ground states in high Landau levels [42][43][44][45]: phase separation is driven by exchange physics and frustrated by Hartree contributions. As noted above, any Hubbard description of TBG must involve substantial further-neighbour interactions. Meanwhile, corrections to the flat band limit -particularly from strain and particle-hole symmetry breaking -can significantly modify the Hartree potential, penalizing full occupation of the mBZ. Since exchange interactions still favor insulating behaviour, one resolution is to reconstruct the bands via finite-wavevector ordering. Thus, from both points of view, it appears that conditions in TBG might favor translational symmetry breaking states over their competitors. Despite this, relatively little work has focused on this possibility, with rare exceptions [16,22,28,[46][47][48][49][50][51], of which we highlight a few. Refs. [16,49] identified a unidirectional charge density wave order that doubles the moiré unit cell in an interaction-only model with spin and valley degeneracy suppressed, but this does not appear to be energetically competitive with translationallyinvariant states in more realistic situations. More recently, a different period-doubling stripe order stabilized by Hartree effects was invoked to explain the occurrence of commensurate-filling insulators whose Chern numbers depart from those expected from the naive QHFM picture [28] (although we suggest an alternative and possibly more natural translation-breaking order at these fillings below). Ref. [48] focuses on the Dirac cones and classifies the possible flavor-breaking orders that can be connected to superconductivity via Wess-Zumino-Witten terms. Among the candidate normal-state orders are moiré density waves which couple the different mini-valleys of TBG (see also Ref. [52] for similar discussion in twisted trilayer graphene) and hence break translation symmetry, but no microscopic studies have yet been performed to determine their energetic competitiveness. The nature of translation symmetry-breaking differs from that in the state described in the subsequent sections, which retains the modified translation symme-tryT . Hence it is not associated with charge/spin modulations between moiré AA regions, unlike the moiré density waves. Ref. [53] studies a valley spiral state in magnetically encapsulated TBG. This is a similar state to our proposed one, albeit the physical origin and the parameter space within which that state exists are quite different 3 . For completeness, we note that translational symmetry breaking has recently been observed in closely related twisted monolayer-bilayer graphene moiré heterostructures [54], and was proposed theoretically to explain insulating states observed in twisted bilayer WSe 2 [55]. But despite these previous works, to date there has been no systematic analysis of translational symmetry breaking in realistic TBG systems, and so the extent to which such symmetry breaking is a common phenomenon across the wide range of parameters relevant to experimental samples remains unclear.
In this work, we explore translational symmetry breaking order at commensurate integer fillings in TBG. Our analysis incorporates three experimentally important deviations from the Bistritzer-Macdonald limit -particlehole symmetry breaking from non-local tunneling perturbations, a substrate potential, and uniaxial strain 4and also studies different interaction strengths and twist angles. In the balance of this introduction, we provide a digest of our main results, which also serves to signpost the organization of the remainder of this paper.

A. Summary of Results
Our central finding, obtained from extensive Hartree-Fock simulations (discussed in Sec II), is that for modest (as little as 0.2%) uniaxial strain and largely independent of all other parameters, the ground state of TBG at all non-zero integer fillings ν = ±1, ±2, ±3 is a timereversal invariant state that breaks superlattice translational symmetry 5 by modulating intervalley coherence at an incommensurate wavevector q. Since the intervalley coherence corresponds to a Kekulé pattern at the microscopic graphene scale, one can view this order as a Kekulé pattern that rotates at the moiré scale with period 2π/|q|. We therefore dub this the 'incommensurate Kekulé spiral' (IKS) order. The IKS state is insulating at ν = ±2, ±3, but does not show a charge gap at ν = ±1. In contrast to its ubiquity at nonzero integer fillings, the IKS order is absent at charge neutrality, where we instead find the ground state for comparable strain to be a nematic semimetal, as identified in previous work [56]. An overview of all the phases found in Hartree-Fock is given in Fig. 2.
Modulations in the valley coherence are fundamental to the IKS state, which hence relies on the interplay between the moiré pattern and graphene-scale physics. This makes its properties and phenomenology distinct both from the previously-studied period-2 stripe states in TBG [28] and from various stripe orders in other correlated systems. It also differs in a few important ways from other proposed states with spatially-modulated valley coherence in either TBG, twisted bilayer WSe 2 or twisted monolayer-bilayer graphene [46,47,55,57]. First, IKS order generally occurs at an incommensurate wavevector, unlike the moiré density waves in Ref. [48]. Second, it does not rely on the presence of higher-order van Hove singularities or Fermi-surface nesting and thus has valley coherence over almost the entire mBZ. Third, the IKS state apparently requires a small, but non-zero, amount of strain to be stabilized against competing q = 0 orders. And finally, the IKS state is time-reversal symmetric, unlike the state discussed for twisted monolayerbilayer graphene in Ref. [57].
Although the IKS order parameter breaks the moiré translationT ai and U (1) V valley symmetries, it preserveŝ T ai =T ai e iq·aiτz/2 (where a i are moiré lattice vectors and τ µ denote a set of Pauli matrices acting in valley space). Performing a valley-dependent gauge transformation therefore yields eigenstates that satisfy a generalized Bloch theorem. This transformation, which amounts to shifting the dispersion in the valleys by ±q/2 in the mBZ, allows us to label electronic states in the IKS state by an analog of the crystal momentum associated withT ai (rather thanT ai ), without folding the mBZ. This shifted-mBZ perspective yields a simple yet quantitatively accurate ansatz for the HF projector in the IKS state, which in turn allows us to link the origin of the order to features of the interaction-renormalized BM bands. We can also use the preservation ofT ai to derive a modified Lieb-Schulz-Mattis theorem that forces gapped IKS order to only appear at integer fillings unless it triggers subsidiary topological or symmetry-breaking order. This explains why, despite involving a modulation that is incommensurate with the moiré lattice, the IKS insulator is tied to densities that are commensurate with it.
As we show below, the precise magnitude and direction of the spiral wavevector q are controlled by details of the dispersion of the interaction-renormalized bands in the mBZ. The former is roughly 1/3 of the width of the mBZ while the latter is approximately aligned along a moiré crystallographic axis, but the energy cost for varying away from these values (the 'stiffness' of the IKS order) is quite small. In a sense one can view this softness as being partly responsible for the robustness of the IKS state, since it allows the ordered state to respond to variations in external parameters such as substrate potential, twist angle, or interaction strength by adjusting q while keeping its other properties essentially unchanged.
The above results are obtained in Sec. III by focusing initially on fillings ν = ±2 where the IKS order is especially simple to describe. We broaden the scope of our analysis in Sec. IV to also consider ν = ±1 and ν = ±3, which differ primarily in the spin structure and the stability of IKS order against competing phases. Interestingly, the time-reversal symmetric IKS order provides a way to obtain insulating states with zero Chern number at the odd integer fillings. Such states are difficult to obtain within the QHFM formalism, and were observed experimentally in Refs. [3,15]. In Sec. V we show that the ν = −3 state provides a 'basis spiral' that serves as a building block for IKS states at other fillings, in a manner that we make quantitatively precise.
We derive a Landau-Ginzburg theory of IKS order in Sec. VI, and use this to link the circular (elliptical) nature of the spiral order in valley coherence to the absence 4 FIG. 1. a) Real-space picture of a IKS state with q = G1/3. The colorplot on the superlattice scale shows the charge density, with darks spots corresponding to AA regions. Black arrows represent the complex IVC order parameter ∼ τxσx +i τyσx . For each of the three inequivalent AA regions, the expectation value of c † A cB + c † B cA on the microscopic graphene bonds is shown. Blue (red) dots correspond to positive (negative) expectation values, and the center of the AA region is marked with a black cross. The different inequivalent AA regions have different approximate √ 3 × √ 3 Kekulé-like patterns on the graphene scale. b) Red (black) line shows BM band structure along a cut in the mBZ for valley K (K ). c) The presence of non-local tunneling, strain and substrate potential breaks various symmetries and affects the dispersion. d) If these single-particle perturbations are weak, the interacting model has an approximate U (4)C=1 × U (4)C=−1 symmetry. Dashed/dotted lines indicate the channel of inter-valley coherence that generically occurs for U (1)V -breaking phases. τx,y denotes any valley-off-diagonal components.
FIG. 2. Summary of phases found in self-consistent Hartree-Fock at all integer fillings ν for different heterostrains and substrate potentials ∆, with non-local tunneling included. Color plot diagnoses IVC order (dark blue indicates unbroken U (1)V symmetry). Dashed red lines indicate approximate phase boundaries, and hatched areas denote absence of a charge gap. Properties of each phase are tabulated in Table I (presence) of subsidiary charge density modulations. We also consider the response of the IKS state to quenched disorder (Sec. VII) and thermal fluctuations (Sec. VIII).
A key point is that although disorder on the microscopic graphene scale can couple to the Kekulé distortion as a random field, this is suppressed in powers of the ratio of the graphene and moiré lattice constants, scaling as θ 4 for small twist angles θ. Consequently, the dominant disorder fluctuations are those that occur on the moiré scale. This justifies our assignment of a definite Kekulé order to each AA-stacking region that defines a superlattice 'site'. Thermal effects are more delicate, owing to a rich set of symmetries broken by the IKS state: namely, valley-U (1), superlattice translation, and three-fold rotation. Since the superlattice translation is broken solely by valley-U (1) charged operators, long range order in both is absent at any temperature T > 0, and is replaced by algebraic correlations, which in turn become exponential decay above a Berezinski-Kosterlitz-Thouless transition temperature T BKT at which vortices in the valley order proliferate. In contrast, the rotational symmetry breaking persists as true long-range T > 0 order, so that (ignoring explicit symmetry breaking from strain) the finitetemperature IKS state has a nematic order up to a finite Ising transition at T N . Depending on the ratio of T N to T BKT , we can have a variety of different thermal melting scenarios based on which order is lost first, though the precise details are subtle and may depend on the ability of the superlattice to pin q to a commensurate value (which is likely weak). The relevant scales for T N and T BKT are similar and are set by the IKS stiffness, and are ∼ 7 K, which is comparable with the experimental temperature scales at which gapped insulators are observed.
We emphasize (Sec. IX) that our results closely match current experiments: most notably, through the absence of spin polarization in the ν = ±2 IKS, the relatively greater robustness (as measured by the charge gap) of insulating states on the electron-doped side (ν > 0), and the ability of the IKS to 'reset' the Chern number to zero at odd integer fillings. The IKS order is a nematic at fi-nite temperature. Thus, we expect the orientational symmetry breaking effects of strain to be heavily enhanced in the IKS state, making it a natural proximate order to the reported nematic superconductors near ν = −2, −3. A subset of experiments find correlated insulators at all integer fillings except at neutrality where they see evidence of a gapless state, and ν = ±1 where weak insulating peaks have been observed [4]. This is readily reconciled with our results -since (see Fig. 2) we find IKS order for all integer ν except ν = 0 -if we assume that the relevant experimental samples are subject to a small amount of heterostrain. (This is a relatively mild assumption given the weak strain needed to stabilize the IKS state and the ν = 0 nematic semimetal.) Direct verification of IKS order is a more challenging goal. Owing to the unusual nature of IKS states, the translational symmetry breaking is invisible to valley-diagonal observables, and a Landau-Ginzburg analysis reveals that the circular spiral order does not trigger a parasitic charge-density wave order. Nevertheless, since intervalley coherence necessarily triggers spatial order on the graphene scale, the associated Kekulé distortion should be apparent, e.g. in the locally-AA regions of the superlattice, but will be modulated at the moiré scale. This order can in principle be detected via scanning probes, though the sensitivity required may be difficult to achieve in the very near term. We close with a discussion of future directions motivated by this work, in Sec. X. We provide details of numerical simulations and additional analysis in five technical appendices.

II. MODEL AND NUMERICAL TECHNIQUES
In this section, we discuss the interacting Bistritzer-MacDonald (BM) model in the presence of strain, substrate potential and non-local tunneling (NLT), and describe our HF calculations. Further details can be found in the Appendices.
We begin with the standard single-particle BM model [5] describing the band structure of two graphene layers l = 1, 2 stacked with a relative twist θ near the magic angle. For concreteness, we orient the coordinate system such that the untwisted monolayer Dirac points lie at k = ±K Dx . The interlayer coupling, which is modulated by the moiré pattern, is parameterized by sublattice-dependent hopping constants w AA = 82.5 meV and w AB = 110 meV. The presence of different coupling constants arises from corrugation effects [58,59] that increase the interlayer spacing in AA stacking regions compared to AB (Bernal) regions.
Throughout this work, we fix θ = 1.08 • , unless stated otherwise. Including both spins and valleys (with corresponding Pauli matrices s µ and τ µ respectively), the BM Hamiltonian has 8 central bands near the neutrality point with narrow bandwidth ∼ 5 meV and large gaps ∼ 30 meV to the remote bands (Fig. 1b). The central band wavefunctions are concentrated at the AA-stacked regions (Fig. 1a), which form the moiré lattice. The model possesses (spinless) TRST = τ xK (whereK denotes complex conjugation) and enjoys emergent D 6 and valley-charge conservation (U (1) V ) symmetries, and an approximate particle-hole symmetry (PHS) [10,11,60]. A related antiunitary symmetryT = τ yK can be defined, which is a signature of the Kramers intervalley coherent (KIVC) phase [20] to be reviewed below. The relevant spatial symmetries of the single-valley BM Hamiltonian areĈ 2zT ,Ĉ 3z , andM , where the latter corresponds to an in-plane two-fold rotation around the x-axis which interchanges the two layers. Spin-orbit coupling is neglected, resulting in a total U (2) K × U (2) K flavor symmetry.

A. Chern basis and effect of a substrate potential
The central bands bear a remarkable resemblance to zero Landau levels in opposite fields (an analogy which is sharpened in the chiral limit w AA = 0 [61]). For a given spin/valley, we can take advantage of the weak dispersion to rotate the pair of central BM bands into a C = ±1 Chern basis by diagonalizing the sublattice operator σ z [20]. Each band carries substantial sublattice polarization (tending to ±1 in the chiral limit), and hence we use σ to also refer to this basis. The Chern number of each of the degenerate bands is tied to its valley according to C = σ z τ z [10,20,62]. [Note that the 'Chern basis' defined by a definite value of C does not coincide with the eigenbasis of the single-particle dispersion in the absence of a substrate potential.] Alignment of TBG with the hBN substrate directly couples to the Chern basis via a sublattice mass ∼ σ z with strength ∆ 10 − 20 meV [63][64][65] (we ignore the additional moiré potential coming from the mismatch between the graphene and hBN lattice constants, although it has recently been argued to be important for explaining some of the experimental features [66,67]), and vio-latesĈ 2z andM . The breaking ofĈ 2zT gaps the Dirac points (Fig. 1c), resulting in the formation of Chern bands. Polarization into a subset of these Chern bands (akin to quantum Hall ferromagnetism) is believed to explain the observation of the anomalous Hall (AH) effect at ν = +3 in aligned samples [14,27]. The substratereconstructed central bands are also used as a starting point for constructing more exotic correlated states [68][69][70][71][72][73].

B. Strain effects
Uniaxial strain of strength = 0.1−0.7% is observed in many TBG samples using STM/STS [74][75][76]. At charge neutrality, this small strain is believed to be an important driving force behind the weakening of symmetry-broken insulators found in numerics at zero strain in favour of semimetallic phases [49,56,77]. In the context of van der Waals homobilayers, it is useful to distinguish between 6 Phase |ν| spin pol. valley pol. U (1)VT = τxKT = τyKTa iT a i e iq·a i τz /2 |C| homostrain, where strain is applied identically to both layers, and heterostrain, where the layers are strained independently. Since homostrain, to first order, does not account for the experimentally observed distortion of the moiré lattice, and has a substantially smaller impact on the electronic structure [78], we focus on heterostrain [56,79], which is also believed to be experimentally relevant. The moiré lattice vectors a 1,2 are deformed depending on the value of the strain ratio and strain angle ϕ with respect to the x-axis. The orthogonal direction is also stretched/compressed due to the Poisson ratio 0.16 [80]. To first order in and θ, the twist angle is unaffected. The BM model is modified by taking into account the deformed superlattice basis vectors, as well as adding an effective layer-dependent vector potential A l (similar to the orbital effect of an in-plane magnetic field [81]). Strain preservesĈ 2z but breaksĈ 3z andM . Hence the Dirac points remain intact, but are unpinned from the K M and K M points and migrate towards the mBZ center [79]. The Dirac points also separate in energy leading to Fermi pockets at charge neutrality (CN), and the overall bandwidth of the central bands increases dramatically (see Fig. 1c).
C. Non-local tunneling and breaking of particle-hole symmetry The standard BM Hamiltonian obeys PHS very wellthe only violations come from small twists in the Dirac cone kinetic terms which are suppressed in θ [11,60]. However, many experiments show pronounced electronhole asymmetry [2,3,14,27,28,82,83], with stronger superconductors on the hole side and more robust insulators on the electron side. We model this PHSbreaking by augmenting the BM model with a non-local interlayer tunneling term [59,[84][85][86]. Consistent with density functional theory calculations, the effect of this term is to make the conduction bands more dispersive than the valence bands [59,84] (Fig. 1c). We use the form of NLT motivated in Ref. [84] and choose values λ 2 = 2λ 1 = 0.18 eVÅ, λ 3 = 0 (see App. A for definitions and a discussion on combining NLT and strain).

D. Hartree-Fock procedure
We perform self-consistent HF calculations on the single-particle Hamiltonian with dual-gate screened Coulomb interactions V (q) = e 2 2 0 r q tanh qd, where the screening length d = 25 nm and relative permittivity r = 10. We neglect terms which scatter electrons between the valleys, as they are suppressed for small θ. Along with electron-phonon scattering, such 'intervalley-Hund's couplings' would weakly break the U (2) K × U (2) K symmetry [20]. In order to avoid double-counting interaction effects, we subtract off a density matrix corresponding to decoupled graphene layers at charge neutrality [20,87]. The results shown the main text were obtained by considering only the central bands as active, with the remote valence/conduction bands frozen to be filled/empty. However we have checked that increasing the number of active bands does not lead to a qualitative change in the results, and mainly leads to a decrease in the bandgap and a slight shift in the phase boundaries. For more details on the effects of adding more active bands, we refer to App. D.
In our numerical simulations, we consider completely general moiré translation symmetry breaking Slater determinants with single-particle density matrices of the general form ĉ † kτ asĉk τ a s = P kτ as;k τ a s , satisfying Tr P = (ν + 4)N 1 N 2 , where N = N 1 N 2 is the number of moiré unit cells, and a, a are BM band indices. We use periodic boundary conditions in both directions which leads to discrete values for the allowed momenta, viz. k, k = n1 For most calculations we enforce co-linearity of the spins and accelerate convergence using the optimal damping algorithm [88,89]. All the states found in our HF simulations (discussed in detail in the following sections) contain at most a single wavevector modulation, meaning that P kτ as;k τ a s = 0 only for k = k, k+q, k−q. As detailed below, whenever translational symmetry breaking occurs we find that the finite-q component of the projector P is entirely off-diagonal in valley space.

A. Numerical Hartree-Fock results
Fig. 3a presents the HF phase diagrams at ν = ±2 in the presence of NLT as a function of both strain and substrate potential. The color scale diagnoses the magnitude of U (1) V -breaking. Without symmetry-breaking perturbations, the lowest energy state is theT -symmetric KIVC state [20] (see also Ref. [92]). It consists of a filled intervalley coherent (IVC) band in each Chern sector, and can be succinctly described by ordering of τ x,y σ y where τ x,y denotes any off-diagonal component in valley space (Fig. 1d). The absence of coherence between opposite Chern sectors sidesteps the energy penalty induced by vortices in the order parameter that would be topologically required for other IVC candidates [64]. At U (2) K × U (2) K level, there is a manifold of degenerate states with different spin polarizations (with maximum 2µ B spin moment per moiré cell), but intervalley-Hund's perturbations will lift this degeneracy.
At a finite substrate potential strength, the optimal state becomes a sublattice-polarized U (1) V -symmetric ferromagnet, which can either be the QAH state ∼ σ z τ z or the valley Hall state (VH) ∼ σ z . These are exactly degenerate at HF level since the VH state is obtained by applyingT on one spin component of the QAH state.
Along the strain axis, we find a first-order transition to a novel phase, which we dub the incommensurate Kekulé spiral (IKS) state, at an experimentally relevant strain ratio of ∼ 0.1 − 0.2% (Fig. 3b). The main characteristic of the IKS state is the breaking of moiré translation symmetry at a single wavevector q. The translation-breaking occurs entirely in the intervalley channel, and is clearly identified in HF by the non-vanishing of the following density matrix elements in the sublattice-polarized basis where spin labels have been omitted. Importantly, the IVC occurs at a single q, leading to a circular intervalley coherent spiral of definite handedness, as there is no symmetry relating the spiral we find in HF to the analogous spiral at −q (Fig. 3c). The IKS state also preserves TRST , and has zero total spin and valley polarization (Table I). Since the spins within each valley are also unpolarized, inclusion of intervalley-Hund's coupling does not lead to qualitative changes. The IKS order persists for fairly large substrate potential strengths. This is expected since the intervalley coherence is flexible enough to polarize onto one sublattice, as evidenced from the fairly constant magnitude of IVC throughout the phase. On the other hand, the KIVC is progressively weakened under increasing sublattice potential, and gives way to U (1) V -preserving ferromagnets since its mechanism relies on inter-sublattice coherence [20]. The strong PHS-breaking effect of NLT manifests in the shifted phase boundaries between ν = −2 and +2. Furthermore, the zero-substrate band gaps (of order 10 meV) of both the KIVC and IKS phases are larger on the electron side than the hole side by 10 − 30%, which is consistent with the experimental trend of more robust insulators at positive fillings.
The numerical phase diagram in Fig. 3a was constructed by restricting to a strain angle ϕ = 0 and periodtripling order along G 1 = |G 1 |x. However, when we relax this requirement, we find that the IKS phase actually consists of a family of spirals which differ only in their ordering wavevector q and are close in energy. Fig. 3c plots the IKS energy relative to best translationsymmetric state, as a function of q 1 (q 2 is fixed to 0, so translation symmetry is maintained along a 2 ). The ideal wavevector q 0 is slightly greater than 1/3 of the mBZ, and evolves weakly with strain magnitude. Hence, the spiral ordering generically occurs at an incommensurate q 0 (see also Fig. 6b). Fig. 3d, where we fully relax the constraints on the wavevector q in our HF calculations, reveals that the dispersion about q 0 is very soft in both directions. Note that we have checked that even without enforcing a particular q, the HF still only converges to a single-q state. We find the energy density of the IKS state to have a term of the Particle-hole breaking is introduced via non-local tunnelling (NLT). System size is 12 × 12, and only central bands are active. b) Phase transition between KIVC and IKS states along the strain axis of a) at ν = −2. In the IKS phase, combined UV (1) rotation by 2π/3 and real-space translation by a1 is a symmetry. c) IKS energy relative to the lowest translation-symmetric state for different strains and spiral wavevectors q along G1. NLT is not included. d) Relative energies of IKS states for different strain angles ϕ (blue axis) and enforced symmetry-breaking wavevectors q (see Eq. 2). Red stars denote q's corresponding to period-tripling along the moiré axes, white circles denote minimum energy wavevectors q0. Non-IKS states that converged to higher energies were discarded. Strain is 0.2%, and NLT is not included. Data points have been transformed to fit on a hexagonal mBZ.
form ρs 2 (q − q 0 ) 2 , from which we estimate the wavevector stiffness to be ρ s ∼ 0.4 meV, without strong spatial anisotropy. In Fig. 3d we show that as the strain angle ϕ rotates, q 0 also changes, but appears to have roughly constant magnitude and predominantly lies near a moiré crystallographic axis. Figs. 3c,d are computed without NLT; including NLT does not affect the qualitative features of these plots.
Before concluding the discussion of our numerical HF results, we want to point out the following subtlety. In the absence ofĈ 3z symmetry (which is broken by strain), the Γ-point of the single-valley BM model is no longer a high symmetry point. From this one might conclude that the choice of Γ in one of the two valleys becomes arbitrary (the Γ-point in the other valley is still fixed by eitherĈ 2z orT ). Making a different choice for Γ does not go without consequences for the IKS state, as this changes the wavevector q at which the inter-valley coherence occurs. For commensurate twist angles, however, there is a preferred Γ-point in the mBZ even in the absence ofĈ 3znamely, it is the point that should fold on top of the Γpoint of the mono-layer graphene BZ (which is fixed bŷ C 2z orT ). From this it is clear that the wavevector q is well-defined for commensurate twist angles, and that the corresponding superlattice translation symmetry is unambiguously broken. In our numerical simulations at incommensurate twist angles, we have always used the same choice for the mBZ Γ-point as in the commensurate twist angle case, such that the location of the mBZ Γ-point varies continuously as a function of θ. However, for incommensurate twist angles a different choice for Γ is possible in principle, and thus the wavevector q of the IKS state becomes 'gauge dependent'. This is consistent with the fact that for incommensurate twist angles, there is strictly speaking no superlattice translation symmetry.

B. Generalized Bloch and Lieb-Schulz-Mattis theorems for IKS states
A defining property of the IKS state, which has circular IVC spiral order, is that it is invariant under the combination of a translation along superlattice vector a i and a valley-U (1) rotation which shifts the IVC angle by a i · q. Let us therefore define modified translation opera-torsT ai ≡T ai e iai·q τz/2 . Because the IKS state preserveŝ T ai , a generalized Bloch theorem applies which states 9 that the single-particle wavefunctions should satisfy ψk(r + a i ) = e ik·ai e −iai·q τz/2 ψk(r). (3) Here,k is a new 'momentum' label restricted to the first mBZ, which differs from the conventional crystal momentum. In particular,k labels real, physical momentã k + τ z q/2 in the two valleys τ z = ±. From Eq. (3), it follows that we can write the single-particle wavefunctions as where uk(r) is the periodic part satisfying uk(r + a i ) = uk(r). As a result, we can define a Hartree-Fock band structure in the mBZ for general IKS states, even if the order wavevector q is incommensurate with the moiré lattice. We note that a similar observation has previously been made for incommensurate circular spin spiral states [99,100].
Another, but closely related, consequence of theT ai symmetry is that the IKS state can only have a non-zero energy gap (ignoring the Goldstone modes) at integer fillings -unless it breaks additional symmetries or develops non-trivial topological order. To see why this is the case, first add a small perturbation of the form to the Hamiltonian. This perturbation preservesT ai , but explicitly breaks the valley-U (1) symmetry. As a result, the Goldstone mode of the IKS state acquires a small gap ∆ G ∝ |h|. Next, we invoke a generalized Lieb-Schulz-Mattis (LSM) theorem which states that the IKS state with gapped Goldstone modes can have a unique ground state on the cylinder geometry which is separated by a non-zero energy gap from all other states in the spectrum only if the charge per unit cell is integer. To show that such a generalized LSM theorem indeed holds, one can simply use the standard adiabatic flux-insertion argument put forward by Oshikawa [101,102]. The redefinition of the translation symmetry operator by multiplying it with e iai·q τz/2 does not change this argument, as the additional factor commutes with the electric-charge U (1) symmetry 6 . In general, one expects that the gapless states which occur at non-integer fillings (excluding 6 One subtlety, however, is that in the simplest version of the adiabatic flux-insertion setup, one wants to have the stripes run along the axis of the cylinder, which means that the system has to be compactified in the direction parallel to q. Forcing states with incommensurate q on the cylinder in this way will cause them to adjust q to take on the closest commensurate value compatible with the periodic boundary conditions. However, because we can make the cylinder circumference arbitrarily big, and thus the shift in q arbitrarily small, this will not change whether or not the system has a gap. topological order and additional symmetry breaking) will have a vanishing charge gap, meaning that it is possible to create well-separated particle-hole pairs with arbitrarily small energy.

C. Structure and energetics of the IKS state
Microscopically, theT -invariant IVC order of the IKS state induces a Kekulé-like pattern on the graphene scale ( Fig. 1), with orientation determined by the local IVC angle θ IVC . Since the symmetry-breaking occurs predominantly within the central bands (see Appendix D, specifically Fig. D.3), the Kekulé or √ 3 × √ 3 pattern, which triples the graphene unit cell, is strongest in the AA regions where the flat band wavefunctions are spatially localized. However, the finite-q character of the IKS state means that the microscopic Kekulé-like patterns differ between different AA regions, as dictated by the combined moiré lattice translation and valley-U (1) rotation symmetryT ai noted above (Fig. 3b). For q = G 1 /3, the system forms stripes along the a 2 direction where the graphene-scale Kekulé pattern is the same. Since the translation-breaking order is purely IVC, with no −q or higher harmonic components, there is no additional charge reconstruction at the moiré scale (see Section VI).
Further insight into the properties of the IKS state can be gained by analyzing its momentum-resolved singleparticle density matrix in more detail. Figure 4b plots the strength of the IVC in momentum space, showing that it is close to the maximum value 1/ √ 2 throughout most of the mBZ. The exceptions are at two lobes in the mBZ, where the electron populations N τ (k) in the two valleys Fig. 4a are close to 0 or 2. The total occupation at each k is 2 consistent with an insulating state. The strong momentum-dependence of the IKS state sets it apart from previously studied mean-field phases [20]. From our numerics we find that the same coherence structure is repeated for both spin species. Therefore, henceforth we consider spin to simply be a spectator degree of freedom, an assumption which is further validated by the 'basis spiral' analysis in Section IV.
The locations of IVC-depletion provide strong clues as to the mechanism underlying IKS formation. In Figure 4c, we calculate for each valley the HF spectrum of the lower band of the self-consistent symmetry-preserving semimetal (SM). This captures the major momentumdependent effects that strain and interactions have on the band structure. The dispersions of the two valleys are related by TRS. All Dirac points lie above E F at ν = −2. Near Γ M , there is a region of very high energy (red) that coincides with one of the Dirac points. There is also a region of low energy (dark blue) lying in some other region of the mBZ. Because of TRS, the low/high energy lobes (indicated by hatched/dotted regions) in the two valleys are related by k → −k.
We now sketch an intuitive picture for how these dispersion features influence the parameters of the IKS or- . System size is 48 × 16, strain is 0.2%, substrate is ∆ = 0 meV, and NLT is not included. c) Dispersions of the lower band for the metastable symmetry-preserving self-consistent SM. Hatched (dotted) lobes, situated near low (high) energy regions, will be predominantly filled (empty) in the IKS state. Black lines indicate the Fermi surface, and black stars mark Dirac point locations. d,e) Schematic construction of the IKS state-a relative momentum boost of the valleys by q allows the lobes to overlap each other. Regions not within the lobes participate strongly in IVC. f) HF spectrum for the same density matrix used in c) except that exchange has been neglected. g) ν = +2 HF phase diagram in the strain-twist angle plane, and the strain-relative permittivity plane. System size is 12 × 12.
der. Figure 4d,e demonstrates that coupling the two valleys at a finite q can pairwise align a high energy lobe with a corresponding low energy lobe in the other valley. In these momentum regions, the system will choose to polarize into the energetically favorable valley (Fig. 4a). Elsewhere, substantial valley hybridization is induced. In this way, the IKS state is able to maximize IVC while respecting the prominent characteristics of the band dispersion. Eachk is equally populated, allowing for an insulating state. Note that attempting to induce IVC at q = 0 instead runs into issues-a large portion (∼ 4×lobe area) of the mBZ would be unable to participate in the IVC since the lobes have small overlap. Furthermore, the total electron occupations would vary as a function ofk, meaning the state cannot be insulating.
This perspective naturally explains the strong kdependence of IVC and the slow variation of the IKS energy with q. The somewhat diffuse features of Fig. 4c mean that for nearby q, the locations/shapes of the lobes only change slightly, leading to a small and roughly isotropic wavevector stiffness. A simple estimate for the ideal wavevector q 0 can be made by connecting the minimum energy momentum in valley K with the maximum energy peak in valley K. The predicted q 0 is broadly consistent with HF results of the IKS state for a range of strain angles ϕ (for details see App. D, in particular Fig. D.7). We emphasize that this scenario opens a gap at E F but, unlike most of the translation-invariant insulators, does not rely on gapping out the Dirac points, which remain high in energy above E F . Instead the k-dependent IVC hybridizes the two valleys at finite q and pulls the occupied band below the rest of the states (Fig. 5a).
We can construct a simple ansatz for the IKS projector in the absence of substrate alignment, which matches the HF numerics extremely well. With σ denoting the sublattice-polarized basis, we define two mutually commuting SU (2) Lie algebras γ = (σ x , τ z σ y , τ z σ z ) and η = (τ x σ x , τ y σ x , τ z ), in terms of which C = γ z [21]. We partially fix the gauge by requiring thatĈ 2z acts as τ x σ x andT as τ xK -the remaining gauge freedom acts 11 as e iϕ k τzσz . The spin-singlet IKS state at ν = −2 can then be parameterized by the projector wherek labels the eigenvalues ofT ai as discussed in the previous subsection, the gauge-variant nk is entirely in the x − y plane, and an identity matrix in spin space is implicit. Across most of the mBZ, the vector mk lies in-plane with a constant angle that can be changed by a global U (1) V -rotation. At the lobes, mk orients towards the poles, reflecting the valley polarization in these momentum regions (recall that η z = τ z ) (Fig. 5b). Since γ x , γ y in Eq. (6) anticommute with the Chern number γ z , there is both inter-Chern and intra-Chern IVC in the IKS state with equal magnitude. This implies that the IVC significantly entangles bands with opposite Chern number, in contrast to the usual U (4) × U (4) ferromagnets found in previous mean-field studies [20]. Importantly, this distinguishes the type of IVC order found here from the uniform T IVC state of Ref. [20] (which also preserves TRS at even integer fillings). In terms of symmetries, the only differences between the T IVC and IKS states is that the former is translationally invariant and cannot be made spin singlet by opposite spin rotations in the two valleys. Both states can be madeĈ 2zT -symmetric by a suitable global valley-U (1) rotation. The IKS projector at ν = +2 can be found by particle-hole conjugation. The construction outlined above, involving the nesting of features of a parent symmetry-preserving band structure, is suggestive of a weak-coupling instability. However the U (1) V -breaking coherence occurs nearly everywhere in the mBZ, instead of just the lobe boundaries. Furthermore, the Fermi surfaces of the SM or the noninteracting BM model generically bear little relation to the momentum structure in the IKS state. Indeed, both strain and interactions play a vital role in renormalizing the central bands and setting the stage for symmetrybreaking phases-the non-interacting BM bands have a total bandwidth ∼ 5 meV, which broaden to ∼ 15 meV in the presence of strain (breakingĈ 3z and shifting the Dirac points up/down), and finally 50 meV with the inclusion of interactions. Strengthening the Coulomb interaction (by reducing r ) favours the strong-coupling states (Fig. 4g). Strain thus effectively tunes the system from strong coupling, where Chern-diagonal ferromagnetic states dominate, to intermediate coupling where other phases (such as the IKS order) emerge that violate the U (4) × U (4) hierarchy. In the language of Ref. [92], strain does not significantly impact the quality of the 'flat-metric condition' (see App. D, in particular Fig. D.17), which is used to prove that the uniform U (4) × U (4) ferromagnet is an exact ground state of the pure interaction Hamiltonian at all integer ν. Instead, strain substantially increases the dispersion, thereby undermining the validity of the perturbative analysis about the ferromagnetic states, and allowing for alternative states such as the IKS to come in. On the other hand, twist angle, which weakly influences the non-interacting central band dispersion, does not have a significant impact on the phase diagram (Fig. 4g).
This strong-to weak-coupling crossover can also be viewed through the lens of direct versus exchange energy [103]. The intra-Chern states at small strain are stabilized by exchange, in analogy with quantum Hall ferromagnetism. At larger strains, including just the Hartree piece of the interaction already recreates the key features of the band renormalization, as shown in Fig. 4f. As verified in Section IV, the IKS state is more competitive farther away from charge neutrality, in harmony with the larger Hartree peak (dip) expected for increasing hole (electron) doping [28]. All particles will feel this increased Hartree potential, while exchange effects are only applicable between electrons of the same flavor. We caution though that this direct-exchange dichotomy is not so clear-cut in practice-separation of the IKS energy into its components reveals that both Hartree and Fock contributions change significantly with comparable magnitude as a function of q. Also, exchange does significantly perturb the band structure of the self-consistent SM, and its inclusion is necessary to obtain reasonable predictions for q 0 . This implies that a proper treatment of both terms is required to adequately capture the physics of TBG for realistic parameters (further details in App. D, in particular Fig. D.8 and D.11).
We now turn to other nonzero fillings, focusing on the same departures from the idealized BM model, fixing NLT as above but exploring the phase structure in the strain-substrate plane. As we show, IKS order appears to be a ubiquitous feature for relatively modest strain, and largely independently of substrate strength.
For ν = ±3, a significantly smaller strain ( ∼ 0.02 − 0.05%) is able to destabilize the QAH state found upon adding interactions to the unperturbed BM Hamiltonian (Fig. 6a). Depending on details of the system parameters, it is possible to nucleate an intervening spin-valley polarized nematic semimetal (see App. D, and in particular Fig. D.2). At zero substrate potential and strain, an alternative period-doubling stripe which preservesĈ 2zT and U (1) V has also been proposed [49]. With the parameters used in the main text, we find a direct transition from the QAH to an IKS state under increasing strain. The IKS state at ν = ±3 has the same symmetry properties as the one at ν = ±2, except that it can carry a spin polarization at U (2) K × U (2) K level (of maximum 1µ B per moiré cell). For (anti-)ferromagnetic Hund's coupling, the spin moments carried by the two valleys will (anti-)align. The dispersion curve in Fig. 6b is also similar (compare Fig. 3c), but it has a deeper minimum, and hence a larger IVC stiffness, due to the increased renormalization of the bands. Given the tiny strains required for the IKS to beat the QAH at ν = ±3, a natural question that arises is whether the IKS can in principle be the ground state at zero strain. We emphasize that there is no fundamental reason that prohibits this scenario from occurring; for strong enough perturbations about the U (4) × U (4) limit, the strong-coupling insulators can be superseded by states outside the QHFM paradigm. While our Hartree-Fock calculations suggest this is not the case for our choice of parameters, the fact that the IKS can be obtained self-consistently without strain (see App. D, in particular Fig. D.14 which shows that energy cost is less than 1 meV per unit cell) is a strong indication that the IKS remains a highly competitive state. Altering details of the model Hamiltonian or choosing different parameters (e.g. increasing the relative permittivity) may well tilt the balance in favor of IKS.
At ν = ±1, the weak substrate potential region of the U (1) V -symmetric QAH phase has non-zero IVC order [15,91,92] (Fig. 6c). The transition from the QAH and QAH+IVC states to the IKS state now occurs at a larger strain ( ∼ 0.15 − 0.2%). The effect of intervalley-Hund's terms is similar to that at ν = ±3, i.e. a net spin polarization of (0) 1 for (anti-)ferromagnetic coupling. In contrast to the other integer fillings, at ν = ±1, the IKS state never develops a charge gap. For completeness, in Fig. 6d we also present the phase diagram at ν = 0, which shows no indications of translation symmetry breaking. The KIVC gives way to a symmetric SM at finite strain [56], and a VH insulator at finite substrate.
All numerically obtained IKS states preserve spinless TRS T , implying that the Chern number C vanishes. This fact is remarkable for the odd fillings, since con-ventional spin-valley polarized ferromagnets can only accommodate phases with odd C. However, recent experiments have shown the existence of even-C gapped phases extending down to zero magnetic field at odd fillings [28]. One possible route to achieving this is by folding the mBZ in half and forming period-2 charge order, as theoretically argued by some authors [28]. Each Chern band (a finite sublattice splitting was considered) splits into a |C| = 1 and 0 mini-band, and a variety of different C states can be obtained by selectively polarizing these. Our work proposes a fundamentally distinct scenario, relying instead on IVC to produce the requisite C = 0 bands, and on moiré translation breaking to minimize the energy. The IKS order is agnostic to the presence of substrate alignment, and is a natural robust insulating candidate for experiments where strain is often an external confounding factor. Furthermore, as explained in detail in the next section, translation-breaking phases with nonzero Chern number can be obtained by 'stacking' a phase with IKS order with other translation-symmetric phases to achieve the requisite band filling. Characterizations of the moiré charge order or strain in the sample of Ref. [28] would help determine which theoretical scenario is operative there.

V. RELATIONSHIP BETWEEN IKS STATES AT DIFFERENT FILLINGS
Since the IKS states have similar properties at all nonzero integer fillings, we expect them to be closely related. To make the connection explicit, we consider the |ν| = 3 IKS state as a 'basis spiral'. We start with the ν = −3 IKS state with spin-polarization enforced for simplicity. In order to construct a ν = −2 IKS, we take two copies of the ν = −3 basis spiral in order to obtain a spinunpolarized IKS. The same construction is possible at positive filling by particle-hole conjugation. The notion of the |ν| = 2 stripe as two copies of the |ν| = 3 stripe is consistent with the relative scale of the U (1) V -breaking order parameter in Fig. 3a being √ 2 times that in Fig. 6a. For the ν = −1 IKS state, we note that the translational symmetry breaking is entirely in one spin sector, whereas the other spin sector has the same symmetries as the VH state at ν = −2. This motivates the following construction: We start with the spin-polarized VH state at ν = −2 and add to it the ν = −3 IKS state in the opposite spin sector. Consistent with this fact, the U (1) V -breaking order parameter has the same magnitude in the IKS phases at |ν| = 1 and 3.
We show in Fig. 7 that the trial states for the ν = −2 IKS state based on this construction have energies that are very close to the self-consistent HF solution at those fillings. We also find that HF simulations using these trial states as initial inputs converge very quickly to the selfconsistent IKS ground state at that filling, demonstrating that the trial states have the correct correlations. For ν = −1 the trial state energies are not as close to those of the self-consistent HF IKS state due to the fact the IKS state is not insulating at this filling (see App. D).

VI. LANDAU-GINZBURG ANALYSIS
In this section, we explore the interplay between Kekulé spiral states and charge order, by generalizing the Landau-Ginzburg analysis of Ref. [104]. This has two primary motivations. First, charge order can be detected in a larger variety of experimental probes: for instance highresolution scanning single-electron transistors (SETs) can detect charge modulation, as can scanning tunneling microscopy (STM) measurements (which can also directly access the moiré-scale-modulated Kekulé distortion characteristic of the IKS state.) Second, charge order is more readily pinned by external potentials, and so an IKS state with a subsidiary charge order is likely to respond more strongly to quenched disorder and also to experience stronger commensuration effects.
Our Landau-Ginzburg construction involves the following order parameters 7 where N is the number of moiré unit cells, andψ † k,τ,σ creates an electron in one of the active Chern bands with Chern number C = τ z σ z . We have also partially fixed the gauge by requiring thatĈ 2zT acts as σ x K andT as τ x K. The notation I refers to either 'IVC' or 'Isospin'. Under the symmetries of strained TBG, the order parameters transform as where I q = (I x q , I y q ), R(φ) is a 2 × 2 rotation matrix, and we have used that I −q = I * q and ρ −q = ρ * q . Note that bothĈ 2z andT act differently on the isospin vector I than on conventional spin-1/2 degrees of freedom. An interesting consequence of these unusual symmetry actions is that the bilinear i ij I i q I j * q with ij = − ji and xy = 1, respects all symmetries of strained TBG, including superlattice translations.
As a disclaimer, we point out that the Landau-Ginzburg free energy we study below is not invariant under threefold rotations. One justification is that strain breaksĈ 3z . However, since the strain is very small, a three-fold rotationally symmetric functional (plus small anisotropies) should actually still be an appropriate starting point. Instead, a better justification for a 'unidirectional' free energy is that we are interested only in uni-directional physics here. In other words, we can think 7 Note that there are other operators we could use instead of I x/y . For example, J x/y = 1 N k ψ † k+q τ x/yψk and K x/y = ± 1 N k ψ † k+q τ y/x σzψ k have transform similarly to I. However, using all these operators in the Landau-Ginzburg free energy would not change the conclusions of our analysis -so for simplicity, and because |I| = √ J 2 + K 2 in the IKS state, we consider only I. of our free energy functional as descending from a parent functional which isĈ 3 symmetric, but which breaks the threefold rotation symmetry spontaneously. The free energy we write down is then obtained by expanding the parent functional around one of the three valleys.
A Landau-Ginzburg free energy consistent with alll the above considerations is given by: where I ≡ I q and ρ ≡ ρ 2q . The terms with coefficients r x , λ 3 and λ 4 are not present in the analysis of Ref. [104].
We will now argue that despite the presence of these additional terms, the physical conclusions of Ref. [104] survive. From now on, we will normalize the order parameters such that U ρ = U I = 1, and require that U x < 1 in order for the free energy to be bounded from below. We will also consider the case with λ 2 = λ 3 = λ 4 = 0, as these terms only quantitatively change the physics we want to discuss here. As a first step, we choose a coordinate system such that the order parameters can be written as whereê i are two orthogonal unit vectors, and 0 ≤ α ≤ π/4. If α = 0 then the IVC order is collinear, whereas if α = π/4 the Kekulé spiral is perfectly circular, which is what we find in HF simulations. For intermediate values of α the Kekulé spiral has a non-zero eccentricity tan α.
Using the above expressions for the order parameters, the free energy can be written as +2λ 1 |I| 2 |ρ| cos 2α cos θ.
Next,we minimize F with respect to θ in the presence of non-zero |I|. From this we find that θ is either 0 or π, with the optimal value being determined by the sign of λ 1 . Note that for θ = 0 or π, the charge order preserves theĈ 2z symmetry (as does the circular Kekulé spiral state). Now, from minimizing F with respect to α, it follows that α = π/4 only if |ρ| = 0. Minimizing F with respect to |ρ| gives us the reverse implication: if F is minimized for |ρ| = 0, then this implies that α = π/4 (and also that r ρ > 0). Combining the above implications, we conclude that there is no charge order if and only if the IKS state is perfectly circular. This interplay between charge order and eccentricity suggests interesting possibilities to experimentally detect the IKS order. For example, let us consider a situation where the TBG sample is strained in a direction orthogonal to the spiral wavevector q. Via the distortion of the graphene bond hoppings, this will introduce some non-zero eccentricity for the IKS state and thus charge order with half the period of the Kekulé spiral. Conversely, if non-zero charge order is induced via an inhomogeneous electrostatic potential with wavevector 2q, then the IKS state will respond by changing the amplitude of the Kekulé pattern differently in inequivalent AA regions.
Minimizing F exactly with respect to α is cumbersome, requiring us to minimize the expression Since we are only interested in the parameter regime close to where the circular spiral state is lowest in energy, let us set 2α = π 2 + δ and expand in powers of δ. Keeping terms to second order in δ, Eq. (17) becomes Minimizing this expression with respect to δ, we find whence we can write the free energy as where the dots stand for terms which do not involve the charge density (note that if one keeps higher orders of δ, then the coefficient of |ρ| 4 will also receive a correction).
It is now clear that there will be a second order phase transition to a phase with non-zero charge order when As discussed above, in the charge ordered phase the Kekulé spiral necessarily becomes an elliptical spiral. It thus follows that by tuning a single parameter it is possible to go from the circular Kekulé spiral phase to an elliptical Kekulé spiral phase via a second order phase transition. As a direction for future work, it would be interesting to identify experimental knobs that can drive such a transition.

VII. QUENCHED DISORDER
For the KIVC state, which occurs at very small strain, quenched disorder is relatively innocuous. The reason is that this state breaks the physical time-reversal symmetry, while physical impurities in graphene are timereversal symmetric, which implies that they cannot couple as random-field disorder to the KIVC order. The IKS state, on the other hand, preserves time-reversal symmetry and is consequently less protected against disorder. For example, because the IKS state has a Kekulé pattern on every AA region, graphene-scale bond disorder couples to the IVC order as a random field.
Motivated by this observation, we now investigate the effect of bond disorder on the graphene scale in more microscopic detail (we focus on bond disorder for concrete-ness -the analysis for graphene-scale potential disorder is similar). In particular, we consider the following disorder Hamiltonian: where R runs over the positions of the A sites, and δ 1 , δ 2 and δ 3 connect each A site to its three neighbouring B sites. In the BM basis, the disorder Hamiltonian becomes where A is the area of the sample, and Here, τ, l, σ respectively denote valley, layer and sublattice, |u(k) τ,n is the periodic part of the BM Bloch states, G are moiré reciprocal lattice vectors, and K is the position of the Dirac point corresponding to the τ = + valley in the graphene Brillouin zone. We have also defined q τ τ = q + (τ − τ )K and f (k) = j=1,2,3 e iδj ·k . Note that at this point we have performed an exact unitary transformation, so n and n run over all the BM bands, and not only over the active bands.
Our HF simulations indicate that all IVC order is car-ried by the active bands, indicating the subspace relevant to studying the order parameter. Let us therefore perform a Schrieffer-Wolff transformation to obtain the effect of H dis on this subspace. We define G ij q (k) ≡ P i G q (k)P j , where P 0 (P 1 ) is the projector onto the active (remote) bands. Assuming that the Fermi energy lies within the active bands, the Schrieffer-Wolff transformed disorder Hamiltonian (up to second order) can be written as follows: where the sum overm runs over the remote bands, and E n,k are the energies of the BM bands. In what follows, we are only interested in the part of H dis which couples as a random field to the IKS state, i.e. the part of H dis which is off-diagonal in the valley indices.
Let us first consider the first-order term in H dis . An important observation is that the flat band wavefunctions vary slowly on the moiré scale. Specifically, the BM eigenvectors satisfy if n is one of the two flat bands. Here, G 1 a basis vector of the moiré reciprocal lattice. This implies that ||G 00 q (k)|| ≈ 0 if |q| 4|G 1 |. From this we conclude that the short-wavelength disorder does not couple directly to the IKS order parameter. Furthermore, long-wavelength disorder, which does couple to the IKS state via the first order term, gets suppressed by a factor θ 2 . Physically, this is because the density of electrons occupying the central bands is very small.
Because the first order term in the Schrieffer-Wolff transformed disorder Hamiltonian is inoperative for short-wavelength bond disorder, let us consider the second order term next. The relevant second order processes which hop an electron from valley τ = + to valley τ = − are shown schematically in Fig. 8. In the figure, the numbers next to each arrow indicate the order in which the virtual hopping processes take place. There are two processes which involve an intermediate electron in the remote conduction bands, and two processes which involve an intermediate hole in the remote valence bands. Note that in Eq. (26), the energy denominators do not have absolute value signs. This means that the energy denominators involving remote conduction band energies will be positive, while energy denominators involving remote valence band energies will be negative. This difference in sign comes from the fact that the virtual electron processes generate terms in the second order Schrieffer-Wolff Hamiltonian of the conventional form ψ † ψ, while the virtual hole processes produce terms of the form ψψ † . The difference in sign thus comes from the fermion anticommutation relations.
From Eq. (26) we can estimate the magnitude of the second order term to be where t dis is the typical energy scale of the bond disorder, and ∆ −1 is the inverse gap between the active bands and the remote bands averaged over the mini-Brillouin zone. The energy gain of pinning to shortwavelength bond disorder is thus reduced by a factor of θ 2 t dis ∆ −1 compared to the energy gain of pinning to long-wavelength disorder. As a result, for physically relevant disorder strengths we expect the Imry-Ma domains associated with short-wavelength bond disorder to be much larger than the size of the AA regions. For this reason, we ignore short-wavelength disorder and focus only on long-wavelength or moiré-scale random-field disorder. We discuss the effect of moiré-scale quenched disorder on the finite-temperature physics in the following section.

VIII. FINITE-TEMPERATURE PHASE TRANSITIONS
The IKS order not only breaks the valley-charge conservation symmetry, but also superlattice translation and three-fold rotation symmetries. It is important to distinguish from the outset the different ways in which the translation and rotation symmetries are broken. In particular, because the IKS states are invariant under a combination of superlattice translation and valley U (1) rotation, only local operators with a non-zero valley charge can detect the translational symmetry breaking. A corollary of this observation is that the translational symmetry-breaking order is replaced by quasi-long range order at non-zero temperature, and is completely lost once vortices of the IVC order proliferate. The rotational symmetry breaking, on the other hand, can be detected by operators which have zero valley charge. This means that the nematic order (ignoring the explicit symmetrybreaking effects of strain), can persist as true long-range order at non-zero temperatures.
Let us now consider the different possible ways for the IKS order to disappear at finite temperature. We first note that the destruction of IKS order is likely to be driven by the proliferation of fluctuating defects in the order parameter, and not by the unbinding of excitons. The reason is that the binding energy of the condensed inter-valley coherent excitons is expected to be of the order of the Coulomb scale (∼ 20 meV), whereas the IKS stiffness extracted from our numerical Hartree-Fock calculations is significantly smaller (∼ 0.4 meV).
Taking the above considerations into account, there are three different possible scenarios for the IKS order to disappear. In the first scenario, a Berezinskii-Kosterlitz-Thouless (BKT) transition occurs at a temperature T BKT , at which the IKS angle disorders via the unbinding of vortex-antivortex pairs. The nematic order, however, persists above T BKT and disappears at a higher temperature T N 8 . In the second scenario, the algebraic IKS order and the nematic order disappear simultaneously via a single phase transition, which most likely is first order. In the third scenario, the nematic order disappears inside the region of the phase diagram with algebraic IKS order, such that T N < T BKT . In App. E, we show that all these different scenarios can indeed be realized for the special case of commensurate IVC spiral states. In the commensurate case, we find that the finitetemperature physics of IKS states is closely related to that of frustrated or generalized XY models, which also harbor both quasi-long range order and discrete symmetry breaking [105][106][107][108][109][110][111]. Exactly which of the three scenarios is realized depends on the ratio ρ s /σ DW , where ρ s is the IVC stiffness, and σ DW is the domain wall tension between two different ground states related byĈ 3 . The physical picture which arises in our study of the commensurate models also suggests that the third scenario with T N < T BKT is excluded in the incommensurate case relevant for TBG. But regardless of which of the three scenarios occurs experimentally, we expect all transition temperatures to be of the order T BKT πρ s /2 ∼ 7 K.
In our discussion of the finite-temperature phase diagram so far we have ignored moiré-scale quenched random-field disorder.
Because TBG is a twodimensional material, reintroducing quenched disorder will, strictly speaking, destroy all ordered phases (breaking both discrete and continuous symmetries). However, if the moiré-scale disorder is sufficiently weak the finite-temperature phase transitions we discussed for the disorder-free case should still leave detectable imprints on measurable quantities. In particular, for the case of discrete symmetry breaking, it is known that the size of the Imry-Ma domains is exponentially large in (σ DW /h) 2 , where σ DW again is the domain-wall tension in the ordered state and h is the RMS random-field disorder strength [112,113]. For weak disorder, the Imry-Ma domains will thus be much larger than the relatively small TBG devices currently being used in experiment. In this case, nematicity will survive the presence of quenched disorder as 'vestigial' order of the IVC spiral state, similarly to what was discussed in Ref. [114] for 3D models of charge density wave order in the context of the cuprate superconductors.

IX. EXPERIMENTS
We have demonstrated that IKS order is ubiquitous at nonzero integer fillings in the presence of small amounts of heterostrain, largely independent of substrate or fine details of twist angle. We now argue that this fact, combined with the absence of IKS order at charge neutrality in favour of a gapless nematic semimetal, provides a unified explanation of several recent experiments. Evidently, this picture requires us to posit that small heterostrain is inveitably present in typical experimental samples; this is however not an unreasonable assumption, particularly given the modest heterostrain needed to stabilize IKS order (which is comparable to experimentally observed strains [74][75][76]). Before proceeding, we note that in Ref. [115] we extend the Hartree-Fock analysis and investigate the effects of heterostrain upon doping away from integer fillings and including IKS order. By considering the chemical potential variations and the Fermi surface degeneracies, we obtain results in the strained regime that are consistent with the experimentally measured compressibility traces ('cascade' transitions) and Landau fans [1-3, 14, 24, 27, 28, 82, 116-124].
At even integer fillings, the IKS state can be distinguished from the KIVC state by only probing the spin physics. In particular, in a small magnetic field the KIVC state at ν = ±2 has a local spin moment of 2µ B per moiŕe unit cell 9 , whereas the IKS state at ν = ±2 has a vanishing local moment. An immediate prediction that follows from this observation is that in samples with negligible strain, which would have a strong KIVC gap at neutrality (assuming there is also negligible hBN alignment), the insulators at ν = ±2 should have a non-zero local spin moment, whereas strained samples with semimetallic behaviour at neutrality should have no (spin or orbital) magnetic moment at ν = ±2. By applying a small strain to an initially unstrained sample, one should therefore observe a strong first order transition associated with an abrupt disappearance of the local moment as one enters the IKS phase. If the strain in experiment can both be slowly increased and decreased, hysteretic behaviour should be observed for the local spin moment around the KIVC-IKS transition. We note that Ref. [3] indeed found evidence for spin-unpolarized insulators at ν = ±2 in a sample which is semi-metallic at neutrality, which is consistent with the IKS/nematic SM scenario under the assumption that their samples are heterostrained at the ∼ 0.1 − 0.2% level.
At ν = ±3, the IKS insulators also have a smaller local magnetic moment than the QAH insulators which occur in the absence of strain. As both states are spin polarized in a small field, this difference is now due to the large orbital moment of the QAH insulators, which is absent in the IKS insulator. However, the easiest way to distinguish the IKS state from the QAH state is via the transverse or Hall resistance R xy . This quantity is zero in the IKS state as dictated by the spinless time-reversal symmetry, but takes on a quantized nonzero value in the QAH state. In Refs. [3,15], insulating states were observed at ν = +3 which show Landau fans in magnetotransport measurements that are consistent with a zero Chern number. Given the semi-metallic behaviour at charge neutrality, we thus expect the samples of Ref. [3,15] to be strained and therefore the insulators at ν = +3 to have IKS order.
As discussed in Sec. VIII and Appendix E, the nematic order of the IKS state survives at finite temperature. We therefore predict that all insulators at ν = ±2, ±3 should show strong interaction-induced nematicity, much stronger than what one would naively expect from the small strain present in the sample. This prediction actually fits perfectly with the experimental observations of Ref. [83], where nematicity was observed in the superconducting dome between ν = −3 and ν = −2. Indeed, unless the insulators at integer fillings are separated from the superconducting dome by a strong first order transition, one would generally expect the insulators and the superconductors to either both be isotropic or both be nematic. In Ref. [115], we show that IKS order persists for a finite range of doping around the integers and survives to temperatures 10 much greater than the experimental T c , suggesting that the superconducting dome could indeed inherit physics from the IKS. Furthermore, the ideal ordering wavevector q (which is naturally soft as illustrated in Fig. 3d) changes continuously as a function of density [115], analogous to the evolution of the nematic axis of the superconductor [83]. Such behavior is harder to rationalize for other rotation symmetry-breaking parent states such as stripes. So even without making any assumptions about the nature of the superconducting state, we can interpret the observations of Ref. [83] as indirect evidence for the IKS state.
While the above evidence is reasonable, it remains to a degree circumstantial. A more definitive diagnostic for IKS order is possible in principle, by detecting the Kekulé pattern at the AA regions directly using STM/STS. Kekulé order in mono-layer graphene induced by mobile adatoms (or substrate vacancies) [128,129] has been measured in Ref. [130], whereas Kekulé order induced by large (isotropic) strain [131,132] has been observed in Ref. [133]. However, only a very small fraction of the total number of electrons, i.e. those occupying the central bands, participate in the IKS order we find at nonzero integer fillings in TBG. As a result, the signal coming from the graphene-scale √ 3 × √ 3 Kekulé pattern in the IKS state will be significantly smaller than that seen in the above-mentioned monolayer experiments, and hence could lie below the present-day experimental resolution. (However, the STM studies of Ref. [134,135] were able to detect a Kekulé distortion in monolayer graphene in a high magnetic field at densities that are comparable to those of TBG. The theory here is based on the approximate SU (4) symmetry of the zeroth Landau level, which is close in spirit to the U (4) × U (4) limit of TBG, but the anisotropies and associated mechanisms are different [136].) In Sec. VI, we discussed how charge order could be induced in the IKS state by changing the eccentricity of the Kekulé spiral. This charge order might be easier to detect experimentally than the Kekulé pattern, for example by using high resolution scanning single-electron transistors made of a carbon nanotube [137,138].

X. DISCUSSION
By combining the results obtained in this work at nonzero integer fillings with the findings of Ref. [56] at neutrality, we conclude that by adding a small amount of uniaxial heterostrain to the BM Hamiltonian one obtains from self-consistent Hartree-Fock a global picture of the TBG phase diagram that is consistent with most or even all experimental observations at integer fillings. In particular, at neutrality moderately strained TBG is semimetallic, at ν = ±1 it is metallic but with a significantly lower carrier density than the non-interacting BM model due to strong IKS order, and at ν = ±2, ±3 it becomes an IKS insulator. On top of this, we argued that several other properties of the IKS states (e.g. spin polarization, Chern number, nematicity) fit nicely with many of the experimental observations. The estimated temperatures at which the insulating IKS states should appear (∼ 7 K) also agree remarkably well with experiment. Furthermore as we explain in Ref. [115], the doped descendants of the finite-strain integer orders examined here exhibit Landau fans and compressibility signatures that are consistent with those measured experimentally.
Our results open up several interesting directions for future work. For example, a generalized Pomeranchuk effect has been observed in TBG [139,140], which causes high-entropy insulators to win over metallic states with increasing temperature. There is at present no theory to fit the IKS state into this Pomeranchuk scenario. For the uniform QHFM states, several groups have recently calculated the collective mode spectrum [141][142][143][144]. In this work, we have not attempted to do this for the IKS states, although it would be helpful to obtain a more complete picture of the low-energy physics at the different integer fillings. Arguably the most interesting direction for future work is to investigate the relation between the IKS states and superconductivity. For example, the KIVC state has recently been argued to be a natural parent state for superconductivity in Refs. [21,145]. A future direction would be to adapt this mechanism to the IKS state, or investigate whether local KIVC correlations are strong enough to allow the same mechanism to be operative. It would also be interesting to consider if a similar approach to coupling superconductivity via topological terms in Ref. [48] could be generalized to the 'dominant lobe' physics of the IKS. Refs. [117,146] have subjected TBG to different amounts of Coulomb screening, and found that the strengths of insulating and superconducting orders seem to anti-correlate. Any satisfactory theory of superconductivity should be able to explain this anticorrelation.

ACKNOWLEDGMENTS
We thank Ashvin Vishwanath and Yahui Zhang for interesting discussions on the subtle nature of the translation symmetry breaking order in the IKS state. NB is supported by a senior postdoctoral research fellowship of the Flanders Research Foundation (FWO). We acknowledge support from the European Research Council under the European Union Horizon 2020 Research and Innovation Programme, Grant Agreement No. 804213-TMCS and from EPSRC Grant EP/S020527/1. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data. Fig. D.1a presents phase diagrams where the number of active bands has been increased to six per spin/valley (note that the system size has been decreased to 9 × 9). NLT has not been included, so the phases at negative fillings are similar. Compared with the central-bands-only results, the phase boundaries move around slightly. The IKS at ν = ±1 is not fully-formed since the extra bands renormalize the central bands further.   D.2 presents a phase diagram at ν = +3 at a different set of parameters that stabilizes an intervening nematic phase between the QAH and IKS at zero substrate. The nematic preserves the combinationĈ 2zT , thereby leaving the Dirac points intact. This gapless state is also fully spin and valley polarized. In fact for the figures in the main text, the nematic is lower in energy than the QAH above a certain value of strain that is, however, larger than the critical strain that leads to the IKS. In particular, the predicted ideal IKS wavevector q 0 (which connects the low-energy minimum of valley K to the high-energy peak of valley K) would be around the same for all cases. dA(∇θ) 2 , where dA is an area element and θ is the IVC angle, we extract the stiffnesses ρ i and the offsets of the minimum of the dispersion q 0i along the two directions for different values of strain. We do so by fitting a parabola to the dispersions and extracting ρ i and q 0i according to the formula Fig. D.11 shows the energy difference between the translationally invariant state and the IKS at a given wavevector broken down into the different energy contributions. The IKS pays an exchange penalty compared to the translationally invariant state, however it has favourable direct and kinetic energy. Fig. D.12 shows the energy of the IKS trial state at ν = −1 compared to both the optimal HF IKS state at that filling and the translationally invariant state. The trial state construction at ν = −1 works less well than at ν = −2. Fig. D.13 shows the phase diagram in strain-w AA space. The translationally invariant ferromagnetic states are favoured as one tends to the chiral limit w AA = 0.
In Fig. D.14 we show the energy of the IKS HF state compared to the translationally invariant state. Even at zero strain, one can find the IKS solution, albeit it is a metastable state that is higher in energy than the translationally invariant ground state by about 0.7meV.
In Fig. D.15, we show the substrate-strain phase diagrams for a different interaction subtraction scheme (see Section C 1). Here we use the 'average' scheme P 0 = 1 2 I (in the central band subspace), which matches the scheme used in Ref. [92] in the limit of zero strain. Hence this choice retains the approximate particle-hole symmetry (we neglect NLT). The phase diagrams are similar to those in Fig. 2. Note that the gaps of the QAH state at ν = 3 are significantly smaller than those of the graphene scheme used in the main text. For our choice of twist angle and interlayer tunneling, the bare BM bands are very narrow. Hence with NLT, the state at zero strain and substrate can sometimes differ from the ones shown in Fig. D.15. However within the resolution of the phase diagram, the states at any finite strain/substrate revert to the expected phases.
In Fig. D.16, we show the substrate-strain phase diagrams for the CN subtraction scheme proposed in Ref. [77], where P 0 is constructed by filling the lower BM band. Hence this choice retains the approximate particle-hole symmetry (in the absence of NLT). One quantitative difference with the graphene scheme is that the strain scales for transitions to the IKS are significantly lower (almost two orders of magnitude). The properties of the phases themselves, and their relative positions, are largely unchanged. Note that the gaps of the QAH state at ν = ±3 are significantly smaller than those of the graphene scheme used in the main text. Also at ν = −3, the state at zero strain and substrate is no longer fully valley polarized, and possesses some IVC.
Given the significantly smaller strain scales of the CN scheme, we briefly comment on why this may be happening and why we believe this scheme is not an appropriate choice (at least for finite strains). We note that the physical properties of the subtraction projectors P 0 for the graphene and 'average' schemes are largely independent of strain. The only dependence arises from the changing Hilbert subspace of the central bands. Owing to the large gap to the remote bands at zero strain, this effect is weak. On the other hand, the CN projector P 0 closely tracks the evolution of the BM bands as a function of strain. For example, P 0 changes dramatically depending on the positions of the Dirac points. This likely explains the contrasting behavior to the graphene and 'average' schemes. Further evidence against the CN scheme is the heightened fragility of the ν = +3 QAH phase against tiny amounts of strain in Fig. D.16, even in the presence of large substrate coupling. This is inconsistent with experiments which have observed a robust zero-field QAH phase at this filling with substrate alignment. We anticipate that fluctuations beyond HF will not significantly change this property of the CN scheme, for reasons that we now sketch. Since the QAH is a translationally-invariant ferromagnet (similar to QHFM states), relaxation of the single Slater-determinant constraint is unlikely to substantially improve its energetics. On the other hand, the IKS has a soft ordering wavevector q and its projector has significant variations across the mBZ, suggesting it may benefit comparatively more from beyond-HF corrections. However we caution that this is a qualitative argument; it would be interesting to verify, e.g. via DMRG, whether the CN scheme indeed persistently produces phase structures that are challenging to reconcile with experimental inputs, such as the absence of the QAH phase for reasonable values of strain and substrate potential.
In Fig. D.17, we plot the form factor λ for momentum transfer −G 1 (for G = 0 the result is trivial due to Bloch function orthonormality). This is relevant for examining the flat-metric condition λ τ,a,b (k; 0, G) = ξ(G)δ a,b used in part of the analysis in Refs. [92,93,143]. For the our parameters, the flat-metric condition is already moderately violated at zero strain for this shell of reciprocal lattice momenta. Strain shuffles the form factors around, but does not significantly change the scale of mBZ variation.   for different strain angles. The predicted q0 connects the minimum energy momentum from the lower band in valley K of the SM, to the maximum energy momentum from the lower band in valley K. ∆q is defined as the Euclidean distance between the predicted and the actual q0, normalized to the length of G1. Strain is 0.2%, filling is ν = −2, NLT is not included, and system size is 24 × 24 (12 × 12) for the SM (IKS). Blue crosses: Exchange is included, and band structure is from a self-consistent SM. Orange triangles: Exchange is not included, and band structure is from a self-consistent SM. Green circles: Exchange is not included, but the state considered is calculated with exchange included (i.e. same state as blue crosses).   D.11. IKS dispersion at ν = +2, for strain = 0.25%, strain angle ϕ = 0 • , system size 12 × 12, and no NLT. We show the energy difference between the IKS state at a given wavevector and the translationally invariant state broken down into the different energy contributions: Total energy (Etot), kinetic energy (E kin ), direct energy (ED) and exchange energy (EE).   D.14. Energy of the q1 = |G1|/3 IKS state (red) compared to the lowest energy translationally invariant state (blue). In order to find the IKS HF state when it is not the ground state, we use an annealing procedure. We decrease the strain value starting from 0.3% in steps of 0.02% and use the IKS HF state from the previous strain run as an input state. The energy of the input state (green) is close in energy to the lowest energy IKS HF state at that strain value (yellow). System size is 12 × 12, ∆ = 0meV and NLT is included. D.15. Summary of phases found in self-consistent Hartree-Fock at all non-negative integer fillings ν for different heterostrains and substrate potentials ∆, without non-local tunneling. Compared to the main text (Fig. 2) which used a graphene subtraction scheme, the calculations here used the 'average' scheme (see discussion in App. D). The axes as are the same as where K = πJ/T , and h a is an integer-valued 'height' field living on the dual lattice. In the above expression, the index µ runs over the values x and y, ∆ µ is a lattice derivative, and we have introduced the notation s a,µ = s a,a+eµ . The sum over the Ising spins can now be done exactly in a straightforward way, and one finds that it enforces the constraint which implies that the even and odd values for the height field decouple. Focusing on the even sector, we can define a new integer valued fieldh as h = 2h. After a change of variable, we obtain withK = K/4. Because the BKT transition occurs whenK ∼ 2, we arrive at the conclusion that the fluctuations of the Ising spins reduce the KT temperature by a factor of four. Contrary to what one might have expected, disordering the Ising variables does not completely disorder the XY spins. This is because the XY spins, even with fluctuating ferro-and anti-ferromagnetic interactions, still want to be collinear. So after disordering the Ising spins, the correlation function of cos(θ) will decay exponentially, whereas the correlation function of cos(2θ) still decays algebraically.

b. Sketch of the complete phase diagram
Having explored the regime where the Ising spins are completely disordered, let us now attempt to sketch the complete phase diagram of the coupled Ising-XY model. We will fix the XY coupling J, which sets the overall energy scale, and ask what happens when we tune J I and the temperature T .
The case where J I = 0 was solved exactly above, and it was found that there is a 'half-KT' transition at a temperature which is four times smaller than the critical temperature for the conventional KT transition. We now consider what happens if we turn on a small Ising coupling J I . As shown in Fig. E.2 (a), this introduces a small domain wall energy of 8J I per unit length for a π-domain wall of the XY spins. The reason that the domain wall energy for the XY spins is not of order J, is that we can lower its energy by flipping all the Ising spins along the domain wall path (see figure). A domain wall for the Ising spins, on the other hand, costs an energy of 4(J + J I ) per unit length because it introduces some frustration for the XY spins and is thus much more costly than the XY domain wall in the regime where J I J. As a result, for small J I we expect that the π-domain walls of the XY spins will proliferate first at a small temperature T ∼ J I . This triggers an Ising phase transition where the cos(θ) variable becomes disordered, but cos(2θ) retains its algebraic order. In more physical terms, this Ising transition is a deconfinement transition of half vortices, which are connected by strings of π-domain walls for the XY spins with a corresponding string tension of 8J I . At higher temperatures, there will be a half-KT transition characterized by the unbinding of half-vortices and half-antivortices.
If J I J, then the Ising spins remain frozen up to very high temperatures. So with increasing temperature, there will first be a conventional KT transition where vortex-antivortex pairs unbind. Although the XY spins are disordered above the KT temperature, the system still breaks the C 4 symmetry. This will be detectable via operators like cos(θ (i,j) − θ (i+1,j) ) and cos(θ (i,j) − θ (i,j+1) ), which are invariant under global rotations of the XY spins, but do detect the breaking of lattice rotation symmetry. At a higher temperature T ∼ J I , there is an Ising transition which restores the C 4 symmetry.
The With C4-breaking staggered magnetic field for the Ising spins. 'D' is the completely disordered phase, 'Dn' is the phase with disordered XY spins but with non-zero nematic order, A is the conventional algebraic phase with quasi-long range order, andÃ is the algebraic phase with deconfined half-vortices and thus exponentially decaying correlations of cos nθ with odd n. was studied in previous works [147,148]. Interestingly, it was found that all these transition lines merge in a single Ising transition line as shown in Fig. E.3(a) (in Refs. [149][150][151], it was argued that the multi-critical points where the Ising and KT lines meet correspond to a CFT with central charge c = 3/2 and emergent supersymmetry). Physically, the single Ising transition line at J I ∼ J is a simultaneous deconfinement and unbinding transition of the halfvortices. Along this line, the temperature is between the critical temperatures of the half-KT and the conventional KT transitions, so from the moment the half vortices are no longer confined in pairs to form full vortices, they also simultaneously unbind from their anti-vortex partners [147,148]. Finally, we should comment on the effects of a small field which explicitly breaks the Ising or C 4 symmetry (which would correspond to strain in TBG). In the regime where J I J, the main effect of a non-zero staggered field ±h i (s i,x − s i,y ) for the Ising spins is to increase the string tension between the half-vortices from 8J I to 8J I + 2h. As a result, the Ising transition between the conventional algebraic phase and the algebraic phase with deconfined half-vortices will shift to higher temperatures. In the regime where J I J, the sharp Ising transition will become a crossover. Slightly above the KT transition, the system should still exhibit much stronger nematicity than expected from the small symmetry-breaking field (small strain). The resulting phase diagram with non-zero staggered field is shown in Fig. E.3(b).

Triangular lattice model
The Hamiltonian of the triangular lattice model consists of a sum of two terms: where the sum in the first term is over nearest neighbour sites. Similar to the square lattice model, the triangular lattice model has an XY spin on every site. The quantity ∆θ ij is defined for nearest neighbour XY spins using the bond orientations shown in Fig. E.4(a). In particular, ∆θ ij is equal to either θ i − θ j or to θ j − θ i , depending on the orientation of the bond. The convention is that if an arrow on a particular bond (see  With these definitions, it is straightforward to see that the ground states of H t consist of uni-directional XY spiral states running along one of the three lattice directions. The spirals have a period of three, as follows from the factor 2π/3 in Eq. (E8). The Hamiltonian is invariant under three-fold rotations, so the XY spirals spontaneously break this symmetry. Similar to the square lattice model, the C 3 breaking happens via the edge degrees of freedom s ij . Once the edge spins pick a direction, they induce a uni-directional spiral order for the XY spins via the term 2πs ij /3 inside the cosine interaction.
We want to point out that H t is not invariant under θ i → −θ i , even if one also simultaneously takes s ij → −s ij . This is because H s is not invariant under s ij → −s ij , as can be seen from the set of zero energy plaquettes in Fig. E.4(b). For our purposes, this does not form a problem because TBG has no symmetry which only changes the sign of the IVC angle. Despite what is suggested at first sight by the use of the bond orientations in Fig. E.4(a), H t does have a C 2z symmetry, which acts as θ r → −θ −r and s r → s −r . This is consistent with how C 2z acts on the IVC angle in TBG. As pointed out in the main text, the IVC spiral states preserve the C 2z symmetry.
Despite the different lattice geometry, and the different period of the XY spirals, the physics of the triangular lattice model is very similar to that of the square lattice model discussed above. For ∆ s J, domain walls between different ground states related by a translation or threefold rotation have an energy on the order of ∆ s per unit length. This is illustrated in Fig. E.5. As a result, at a temperature T 3 ∼ ∆ s , the nematic order is lost and correlation functions of cos nθ decay exponentially if n is not a multiple of three. The operators cos 3nθ, on the other hand, retain their algebraic order. The XY spins become completely disordered above a temperature T KT /3 ∼ πJ/18, which is smaller than the conventional KT temperature by a factor of 9. For ∆ s J, the nematic order survives to very high temperatures T 3 ∼ ∆ s . The XY spins disorder at a lower temperature T KT ∼ πJ/2 via a conventional KT transition.
The remaining question is how the regions of the phase diagram with ∆ s J and ∆ s J are connected at intermediate ∆ s ∼ J. An interesting possibility suggested by the square lattice model is that there is a single second order transition in the universality class of the 3−state Potts model which separates the algebraic phase from the completely disordered phase. One reason to consider this possibility is that the physical picture underlying the existence of such a direct second order transition in the square lattice model can be applied to the triangular lattice model as well. In this picture, the direct second order transition would be a simultaneous deconfinement and unbinding transition of fractional vortices which are 1/3 of a conventional XY vortex.
We want to point out that similar physics was discussed in a very recent paper which used a generalized XY model to study hexatic-nematic liquid crystals [152]. In particular, the authors of that work found numerical evidence for continuous transitions in the 3−state Potts universality class, which in our triangular lattice model would correspond to the nematic transition at T 3 ∼ ∆ s for ∆ s J, and the fractional vortex deconfinement transition at T 3 ∼ ∆ s for ∆ s J. In the generalized XY model of Ref. [152], however, all phase transitions seem to merge in a single multi-critical point and there is no evidence for a direct second order transition between the algebraic phase and the completely disordered phase (more results are needed to confirm this). We leave a detailed numerical Monte Carlo study of the triangular lattice model phase diagram, and its connection to the phase diagram of hexatic-nematic liquid crystals, for future work.