---
author:
- Yu Saito
- Jingyuan Ge
- Louk Rademaker
- Kenji Watanabe
- Takashi Taniguchi
- 'Dmitry A. Abanin'
- 'Andrea F. Young'
bibliography:
- 'references.bib'
title: Hofstadter subband ferromagnetism and symmetry broken Chern insulators in twisted bilayer graphene
---

**In bilayer graphene rotationally faulted to $\bm{\theta\approx 1.1^\circ}$, interlayer tunneling and rotational misalignment conspire to create a pair of low energy flat bands[@bistritzer_moire_2011-1] that have been found to host a variety of insulating, superconducting, and magnetic states at partial filling[@cao_correlated_2018; @cao_unconventional_2018; @yankowitz_tuning_2019; @lu_superconductors_2019; @sharpe_emergent_2019; @serlin_intrinsic_2020]. Most work to date has focused on the zero magnetic field phase diagram, with magnetic field used as a probe of the $\bm{B = 0}$ band structure. Here, we show that twisted bilayer graphene (tBLG) in a magnetic fields of several Tesla hosts a cascade of ferromagnetic Chern insulators with Chern number $\bm{|C| = 1, 2}$ and $\bm{3}$. We argue that the emergence of the Chern insulators is driven by the interplay of the moiré superlattice with the magnetic field, which endow the flat bands with a substructure of topologically nontrivial subbands characteristic of the Hofstadter butterfly[@hofstadter_energy_1976; @thouless_quantized_1982]. The new phases can be accounted for in a Stoner picture[@nomura_quantum_2006] in which exchange interactions favor polarization into one or more spin- and valley-isospin flavors; in contrast to conventional quantum Hall ferromagnets, however, electrons polarize into between one and four copies of a single Hofstadter sub-band with Chern number $\bm{C = -1}$[@bistritzer_moire_2011-1; @moon_energy_2012; @zhang_landau_2019]. In the case of the $\bm{C = \pm 3}$ insulators in particular, magnetic field catalyzes a first order phase transition from the spin- and valley-unpolarized $\bm{B=0}$ state into the ferromagnetic state. Distinct from other moiré heterostructures[@dean_hofstadters_2013; @ponomarenko_cloning_2013; @hunt_massive_2013], tBLG realizes the strong-lattice limit of the Hofstadter problem and hosts Coulomb interactions that are comparable to the full bandwidth $\bm{W}$ and are consequently much stronger than the width of the individual Hofstadter subbands. In our experimental data, the dominance of Coulomb interactions manifests through the appearance of Chern insulating states with spontaneously broken superlattice symmetry at half filling of a $\bm{C=-2}$ subband[@wang_evidence_2015; @spanton_observation_2018]. Our experiments show that that tBLG may be an ideal venue to explore the strong interaction limit within partially filled Hofstadter bands.**

The energy spectrum of a two dimensional electron subjected simultaneously to a magnetic field and periodic potential is described by a fractal structure known as the Hofstadter butterfly[@hofstadter_energy_1976]. For a system with $n$ electronic bands at zero magnetic field, the Hofstadter spectrum at a magnetic flux per unit cell of $\Phi/\Phi_0=p/q$ hosts $n\times q$ ‘subbands’ (here $\Phi_0=h/e$ is the magnetic flux quantum and $p$ and $q$ are integers with greatest common divisor 1). Remarkably, these subbands are each characterized by a nonzero, integer-valued topological index, known as a Chern number $C$, which describes their contribution to the quantized Hall conductivity when an integer number of subbands are filled[@thouless_quantized_1982]. Owing to their high sample quality and the large area of the superlattice, moiré van der Waals heterostructures have provided a versatile materials system for exploring the physics of the Hofstadter subbands[@dean_hofstadters_2013; @ponomarenko_cloning_2013; @hunt_massive_2013], including the observation of correlation driven Chern insulators at fractional filling of a Hofstadter subband that either spontaneously break the superlattice symmetry[@wang_evidence_2015] or fractionalize electrons into anyons[@spanton_observation_2018]. In typical moiré heterostructures of graphene and hexagonal boron nitride (hBN), however, accessing the Hofstadter regime experimentally requires strong magnetic fields. This is tied to the fact that the superlattice potential induced by the moiré pattern is weak, and provides only a small perturbation to the conventional Landau levels. In this limit, Hofstadter subbands form within a single Landau level, and are restricted to an energy bandwidth per Landau level proportional to $W\propto \exp\left(-\frac{\pi}{2}\frac{\Phi_0}{\Phi}\right)$[@dean_fractional_2020]. Achieving separations between subbands that exceed the energy scales related to disorder or temperature thus requires a significant fraction of the magnetic flux quantum to be threaded through each superlattice unit cell.

Electrons in the tBLG flat bands are strongly localized and should consequently realize a strong-lattice limit of the Hofstadter butterfly in finite magnetic field. In this limit, topologically nontrivial subbands may be separated by a significant fraction of the total bandwidth even for $\Phi \ll \Phi_0$[@bistritzer_moire_2011; @moon_energy_2012; @zhang_landau_2019]. However, discussions in the literature of the magnetotransport properties of tBLG have focused on relating the observed magnetoresistance features to the enigmatic $B=0$ phase diagram, without regard for the interplay of the magnetic flux with the periodic potential. Here, we show that magnetic fields can indeed drive new symmetry breaking transitions, producing magnetoresistance features at low magnetic fields that are not correlated with the realized $B=0$ ground states. Fig. \[fig:1\]a shows longitudinal resistivity ($\rho_\mathrm{xx}$) data from a $\theta=1.12^\circ$ tBLG device (see Fig. \[fig:device\_contact\]) as a function of magnetic field ($B$). We present data at 4 K to suppress all but the most energetically robust features; Hall resistance ($\rho_\mathrm{xy}$) data and lower temperature $\rho_\mathrm{xx}$ data are available in Figs. \[fig:Full LL supple\] and \[fig:Hall\_landaufan\]. At $B = 0$, our device shows familiar features of flat band tBLG: correlated insulators or resistivity peaks at $\nu=\pm2$, $\pm3$, and $0$ and a robust superconducting state for $\nu\lesssim-2$ (see Figs. \[fig:1\]c and \[fig:SC\]). Here $\nu$ indicates the number of electrons per superlattice unit cell. However, at $B\gtrsim 5$ T (for $\nu<0$) and $B\gtrsim2$ T (for $\nu>0$), this behavior abruptly changes, giving way to a series of Chern insulators, indicated in blue in Fig. \[fig:1\]b. They are characterized by $\rho_\mathrm{xx}$ minima and integer quantized Hall conductivity that increases in magnitude as the band is filled.

Within the Hofstadter picture, gapped states follow linear trajectories according to a Diophantine equation, $\nu = tn_\mathrm{\Phi} + s$ ($s, t \in \mathbb{Z}$) in the space of $\nu$ and $n_\mathrm{\Phi} \equiv \Phi/\Phi_0$[@wannier_result_1978], the magnetic flux per moire unit cell. Here $s$ is the Bloch band filling index, which encodes the number of electrons bound to each lattice unit cell, while $t$ is the total Chern number associated with a given gap, and consequently the quantized Hall conductivity[@streda_quantised_1982; @thouless_quantized_1982]. Starting from $\nu=-4$ and moving into the flat band, at $n_\mathrm{\Phi} = 0.3$ the observed sequence follows $(t,s)=(-1,-3), (-2,-2), (-3,-1)$ and $(-4,0)$. This pattern is consistent with sequentially filling four bands with Chern number $C=-1$, each of which binds a single electron per moiré unit cell.

The observed transition between low and high magnetic field regimes can be understood from the evolution of the tBLG flat bands in a magnetic field. At $B=0$, tBLG hosts two low-energy bands (for a given spin and valley flavor) connected by two Dirac point that remain gapless absent symmetry breaking by substrate potentials[@po_origin_2018; @kang_symmetry_2018; @koshino_maximally_2018]. A key feature of the $B=0$ band structure is the triangular form of the density of states, shown in Fig. \[fig:1\]e. It was recently argued[@zondiner_cascade_2020] that exchange interactions within such a band lead to a peculiar form of Stoner ferromagnetism in which exchange interactions favor collective states that are symmetric between two, three, or four isospin flavors. In our device, the low $B$ quantum oscillations (see Fig. \[fig:low field fan\]) are qualitatively consistent with this mechanism, showing signatures of Dirac fermions with fully broken flavor symmetry for $\nu<-3$, twofold symmetry for $-3<\nu<-2$, and fourfold symmetry between $-2<\nu<0$ (see insets to Fig. \[fig:1\]c and \[fig:low field fan\]) Notably, in our device, and indeed in the majority of published transport data in tBLG devices, symmetry breaking is absent at $\nu=\pm1$, as evidenced by the absence of low-$B$ quantum oscillations near those fillings. This effect is not explained within the model of Reference , but may arise due to the reconstruction of the low energy bands that is expected to occur at partial fillings due to the Coulomb interactions[@guinea_electrostatic_2018].

An out of plane magnetic field significantly reconstructs the low energy bands. Figure \[fig:1\]f shows the evolution of electronic structure as a function of magnetic field, while Fig. \[fig:1\]g shows the calculated density of states (DOS) at $n_\mathrm{\Phi}=2/5$, corresponding to $B=12$ T for our device. At low magnetic fields, the magnetic field reconstructs the two flat bands into three groups of subbands (for each spin and valley) having net Chern numbers $-1$, $+2$, and $-1$ respectively. This structure is generic to all tBLGs, with the $C=2$ band originating in the zero energy Landau level in the two constituent graphene monolayers and the $C=-1$ bands ensuring that the combined bands maintain net Chern number of zero. Even in moderate fields the DOS is characterized by three effective subbands that contrast to the two bands at $B=0$. The cascade of Chern insulators is well explained by the sequential filling of four spin and valley-projected copies of the lower $C = -1$ subbands, as illustrated in Fig \[fig:1\]d and as emerges from a simplified analysis of a Stoner model that accounts for the magnetic-field induced change in bandwidth (see Supplementary information and Fig. \[Fig:SupplTheory2\]). This picture is mirrored for electron doping, and is also consistent with the Landau fan originating from the CNP, which shows the strongest features at $(t,s)=(\pm4,0),(\pm2,0)$, and $(0,0)$, consistent with successive fillings of the central $C=2$.

The most striking features of the Chern insulator cascade are the states at $(t,s)=(\pm 3, \pm 1)$(Fig. \[fig:2\] and \[fig:p\_landaufan\], which are not correlated with a set of low-magnetic field quantum oscillations (Fig. \[fig:low field fan\]). Indeed, low temperature measurements show that these states emerge without warning in the midst of well formed four-component quantum oscillations originating from the charge neutrality point, consistent with previous reports of orbital magnetic states at finite $B$ [@lu_superconductors_2019; @stepanov_interplay_2019]. The sudden phase transition (Fig. \[fig:2\]a and b) and sudden appearance of an activation gap (Fig. \[fig:gap\_temp\]) can be understood within the context of Stoner ferromagnetism that accounts for the Chern character of the $C=-1$ Hofstadter subband. Within this picture, symmetry breaking occurs when the density of states $\mathrm{DOS}>1/U$, where $U$ is the exchange interaction strength. As the magnetic field increases, both the bandwidth of and total number of states within the $C = -1$ subband decrease; numerically, however, this leads to an increase in average DOS (see Fig. \[Fig:SupplTheoryDOS\] and supplementary information), favoring symmetry breaking above some critical $B$. The nature of this first order transition is illustrated in Fig. 2c, which shows a line cut connecting the $(-3,-1)$ and $(-8,0)$ states. At $(-3,-1)$, three $C=-1$ subbands are completely filled and one empty; while at $(-8,0)$ all four flavors are equally occupied, filling four copies of a $C = -2$ band that includes the entire $C=-1$ band except for the highest energy state associated with a single, strain split Landau level in each flavor[@zhang_landau_2019].

In one view, the Chern insulator states can be seen as competing ground states in twisted bilayer graphene at integer band filling and $B=0$[@xie_nature_2020]. It is known that small perturbations such as an aligned hBN substrate can favor Chern insulators at $B=0$, for instance at $\nu=3$[@serlin_intrinsic_2020]. Even absent this splitting, however, a retrospective review shows that nearly all published data containing magnetotransport data show a feature likely associated with the transition we report here to one or more of the Chern insulator states we observe here[@cao_unconventional_2018; @yankowitz_tuning_2019; @lu_superconductors_2019; @stepanov_interplay_2019; @uri_mapping_2019]. These transitions can be generically understood as reflecting the finite orbital moment of the Chern insulator states, which lowers their energy relative to competing $C=0$ insulators or semimetallic ground states in a finite magnetic field. However, we prefer to analyze the observed ferromagnetic Chern insulator states as generalizations of quantum Hall ferromagnets in which the Landau levels are replaced by $C=-1$ subbands of the tBLG Hofstadter spectrum. We note that just as in graphene Landau levels[@dean_fractional_2020], each Chern ferromagnet corresponds to polarization within the combined space of spins and valleys with the precise isospin ordering set by the interplay of higher order effects including the Zeeman effect and the anisotropy of the Coulomb interactions themselves. Interestingly, these anisotropies may differ within the Chern bands as compared to both Landau levels and the $B=0$ flat bands, possibly favoring different broken symmetries at integer fillings.

This approach moreover has the advantage of providing a quantitative description of the electronic substructure of the C=-1 band at high magnetic fields. Figure \[fig:3\]a and b display measured $\rho_\mathrm{xx}$ and $1/\rho_\mathrm{xy}$, respectively, from 6 T to 12 T in a density region corresponding to occupying the third and fourth copies of the lower $C=-1$ bands. The low disorder in our graphite gated sample[@zibrov_tunable_2017] allows us to resolve a number of smaller energy gaps characterized within the $C=-1$ band. First, we resolve the next level of subband structure within the $C=-1$ subband, which is composed of well separated clusters of subbands with net $C=+1$ and $C=-2$ (see Fig. \[fig:3\]c-d). These features are repeated in experimental data, consistent with the scenario of sequential filling of the overall $C=-1$ band described above. In addition, we also observe robust states showing $(t, s) = (-3, -1/2)$ and $(-2, -3/2)$ at half filling of the $C=-2$ band (highlighted in red in Fig. \[fig:3\]c). These states show quantized $1/\rho_\mathrm{xy}=t e^2/h$ as expected from the measured slope $t$, and an activation gap of around 1 K (Fig. \[fig:3\]e, \[fig:sbci\] and \[fig:sbci\_gap\]). Similar states are observed at positive fillings as well (Fig. \[fig:sbci\]).

The quantum number $s$ encodes the total number of electrons bound to a given lattice site at fixed magnetic flux, and is an integer for all states that can be described within a noninteracting Hofstadter butterfly. The simplest mechanism that allows this number to be fractional is for the superlattice symmetry to break, and we associated the observation of states with half integer $s$ to the formation of symmetry broken Chern insulators, in which electronic interactions within a $C=-2$ Hofstadter subband drive spontaneous doubling of the unit cell. While such states were previously observed in hBN-aligned graphene[@wang_evidence_2015; @spanton_observation_2018], they typically appear only at magnetic fields in excess of 16 T, more than twice as high as in tBLG. The relative strength of the SBCI states can be related to the strong lattice limit of the Hofstadter band structure realized in tBLG. Specifically, within the Hofstadter subbands of a weak-superlattice system the Coulomb interaction scale is never stronger than $e^2/\ell_B$ (here $\ell_B=\sqrt{\hbar/eB}$ is the magnetic length), which vanishes at low $B$. In contrast, in tBLG and other moire flat band systems the magnetic-field independent moire wavelength sets the scale for interactions. Interactions can drive the formation of correlated states at partial filling of Hofsatdter subbands in the low field limit, just as they drive correlated states at $B=0$. Our experiment suggests that other exotic ground states at fractional fillings of Chern bands may be accessible in moire flat band systems, potentially even at $B=0$ in the presence of appropriate time-reversal symmetry breaking order. *Note:* During the preparation of this work we became aware of two additional reports reporting some of the same observations[@nuckolls_strongly_2020; @wu_chern_2020].

Methods {#methods .unnumbered}
=======

**Device fabrication**\
The tBLG used in this study were fabricated using a “cut-and-stack” technique[@saito_independent_2020]. Prior to stacking, we first cut graphene into two pieces using AFM to prevent the unintentional strain in tearing graphene. We used a poly(bisphenol A carbonate) (PC)/polydimethylsiloxane (PDMS) stamp mounted on a glass slide to pick up a 40-nm-thick hBN flake (typically 30$-$50 nm) at 90$-$110$^\circ$C, and carefully pick up the 1st half of a pre-cut graphene piece, rotate and pick up again the 2nd half of graphene piece in series at 25 $^\circ$C using the hBN flake. Here we rotated graphene pieces manually by a twist angle of about 1.2$^\circ-$1.3$^\circ$. Finally, the 3-layer stack (hBN-tBLG) is transferred onto another stack for the bottom gate part (hBN-graphite or graphene gate), which is prepared in advance by the same dry transfer process and cleaned by the typical solvent wash using chloroform, acetone, methanol and IPA followed by vacuum annealing (400$^\circ$C for 8 hours) to remove the residue of PC film on the hBN surface. We did neither squeeze the bubbles nor perform any heat annealing after the stack is completed to prevent the relaxation of tBLG. Electrical connections to the tBLG were made by CHF$_3$/O$_3$ etching and deposition of the Cr/Pd/Au (2/15/180 nm) metal edge-contacts[@wang_one-dimensional_2013]. The device used in this study is identical to the Device 5 in the previous study[@saito_independent_2020].

**Transport measurements**\
All transport measurements in this study were carried out in a dilution refrigerator (Bluefors LD400) with a nominal base temperature of 10 mK, which is equipped with a 14 T superconducting magnet and heavy RF and audio frequency filtering with a cutoff frequency of $\sim$ 10 kHz. The temperature dependent measurements were done by controlling the temperature using a heater mounted on a mixing chamber plate. Standard low frequency lock-in techniques with Stanford Research SR860 amplifiers were used to measure the resistivity $\rho_\mathrm{xx}$ and $\rho_\mathrm{xy}$ with an excitation current of 1$-$3 nA at a frequency of 17.777 Hz.

**Twist angle determination**\
The twist angle $\theta$ is determined from the values of charge carrier density at which the insulating states at $n_{\nu = \pm 4}$ are observed, following $n_{\nu = + 4} = 8 \theta^2/\sqrt{3}a^2$ , where $a$= 0.246 nm is the lattice constant of graphene. The values of $n_{\nu = \pm 4}$ are determined from the sequence of quantum oscillations in a magnetic field that project to $n_{\nu = \pm 4}$ or $n_{\nu = \pm 2}$ for devices.

Acknowledgments {#acknowledgments .unnumbered}
===============

The authors acknowledge discussions with E. Berg, I. Protopopov, T. Senthil, and M. Zaletel. The experimental work was primarily supported by the ARO under W911NF-17-1-0323. Y.S. acknowledges the support of the Elings Prize Fellowship from the California NanoSystems Institute at University of California, Santa Barbara. K.W. and T.T. acknowledge support from the Elemental Strategy Initiative conducted by the MEXT, Japan and the CREST (JPMJCR15F3), JST. A.F.Y. acknowledges the support of the David and Lucille Packard Foundation and the Alfred P. Sloan foundation. L.R. was supported by the Swiss National Science Foundation via an Ambizione grant. D.A. acknowledges the support of the Swiss National Science Foundation.

Author Contributions {#author-contributions .unnumbered}
====================

Y.S. and J.G. fabricated tBLG devices. Y.S. performed the measurements and analyzed the data. L.R. and D.A. performed the theoretical calculations. Y.S. and A.F.Y wrote the paper with input from L.R. and D.A. T.T. and K.W. grew the hBN crystals.

Competing interests {#competing-interests .unnumbered}
===================

The authors declare no competing financial interests.

**Supplementary Information for\
Hofstadter subband ferromagnetism and symmetry broken Chern insulators in twisted bilayer graphene**\
Yu Saito, Jingyuan Ge, Louk Rademaker, Kenji Watanabe, Takashi Taniguchi, Dmitry A. Abanin and Andrea F. Young\
Correspondence to: andrea@physics.ucsb.edu\

![image](device_contact.pdf){width="\columnwidth"}

Low field quantum oscillation analysis
--------------------------------------

The low charge- and twist angle disorder of our device allows us to resolve quantum oscillations down to very low magnetic fields. Figure \[fig:low field fan\] shows quantum oscillation data for $-4<\nu<0$. At a magnetic of 0.6 T, quantum oscillations are not visible except originating from the charge neutrality point. Notably, the strongest features are at Landau level filling $\nu_{LL}=\pm4, \pm12, \pm20$ (Fig. \[fig:low field fan\]c), precisely as expected from the noninteracting band structure of twisted bilayer graphene[@bistritzer_moire_2011]. This sequence can be understood as twice that of monolayer graphene, and indeed arises from the two gapless Dirac cones of the constituent graphene monolayers.

This is in marked contrast to published quantum oscillations in other twisted bilayers which uniformly show a sequence of $\pm4$, $\pm8$, $\pm12$, etc. However, we see the emergence of the more typical sequence at higher magnetic fields, for example in the data taken at $B = 2$ T shown in Fig. \[fig:low field fan\]c. By this magnetic field, we also see quantum oscillations near $\nu=-2$ and $\nu=-3$, with sequences of $-2, -4, ,6,...$ and $-1,-2,-3,-4....$, respectively. The transition in behavior may arise from a single particle effect such as strain[@zhang_landau_2019], which may well be ubiquitous in graphene devices at some level. However, conventional quantum Hall ferromagnetism may also be at play, as evidenced by the emergence of additional quantum oscillation minima at $\nu_\mathrm{LL}=1,2,3$ and $6$.

The data can be elegantly explained under an ansatz consistent with that taken in Ref. [@zondiner_cascade_2020], with the modification that no three-fold symmetry state exists. In this picture, $-4<\nu<-3$ hosts a one component Fermi surface, $-3<\nu<-2$ has a two component Fermi surface, and $-2<\nu<0$ has a four component Fermi surface. A key conclusion from our data is that the failure to observe the expected sequences near $\nu=-2$ (where two copies of gapless Dirac, $\nu_\mathrm{LL}=-2,-6,-10$ would be expected) or $\nu=-3$ (where a single Dirac point is expected, $\nu_\mathrm{LL}=-1,3,-5$) are consistent with the breakdown of symmetry observed at neutrality at the magnetic fields where oscillations are observed in the other fans. Cleaner samples, allowing measurements of the $\nu-2$ and $\nu=-3$ Landau fans, can be expected to similarly reveal the underlying Dirac physics.

Continuum model
---------------

We calculated the zero-field dispersion and density of states using the continuum model of twisted bilayer graphene[@bistritzer_moire_2011-1; @koshino_maximally_2018]. Given the large amount of literature written on this model, we give only a succinct summary here.

For a twist angle $\theta$ and the graphene lattice constant $a = 0.246$ nm, the distance between the ${\bf K}_\ell$ points of the two layers $\ell = 1,2$ in momentum space is $k_\theta = \frac{4 \pi}{3a} 2 \sin \theta/2$. The corresponding Moiré reciprocal lattice vectors are $${\bf G}^M_1 = k_\theta ( - \sqrt{3}, -3) /2, \; \; {\bf G}^M_2 = k_\theta (\sqrt{3},0).$$ Consequently, the Moiré lattice constant is $a_M = \frac{a}{2 \sin \theta/2}$. For the intralayer Hamiltonian in the ${\bf K}$-valley we take $$H_\ell ({\bf k}) = -v \left[ R(\pm \theta/2) ({\bf k} - {\bf K}_\ell) \right] \cdot {\bf \sigma}
    \label{Eq:IntraLayer}$$ where $R(\pm \theta/2)$ is a rotation matrix by angle $\theta/2$ for layers $\ell =1,2$, respectively, and $\sigma=(\sigma_x,\sigma_y)$ are the Pauli matrices acting in the subspace of two graphene sublattices. We set the Dirac velocity to be $v = 2.51 a$ eV, where $a$ is the graphene lattice constant. We have taken $\hbar=1$.

The interlayer Moiré coupling $H_M$ allows for scattering between neighboring Moiré Brillouin zones. It is common to write this Hamiltonian in real space. For the ${\bf K}$-valley, it is given by: $$\begin{aligned}
    H_M ({\bf r}) &=& w \begin{pmatrix}
        \alpha & 1 \\ 1 & \alpha 
    \end{pmatrix}     + w \begin{pmatrix}
        \alpha & e^{-i \phi} \\ e^{i \phi} & \alpha 
    \end{pmatrix}e^{i  {\bf G}^M_1 \cdot {\bf r}} 
    + w \begin{pmatrix}
        \alpha &  e^{i \phi}\\ e^{-i \phi} & \alpha
    \end{pmatrix} e^{i  ({\bf G}^M_1 + {\bf G}^M_2) \cdot {\bf r}}\end{aligned}$$ where $\phi = 2\pi/3$. Note that in momentum-space, a term of the form $e^{i  {\bf G}^M_1 \cdot {\bf r}} $ allows scattering from ${\bf k}$ to ${\bf k} + {\bf G}^M_1$. We set $w=110$ meV, and $\alpha = 0.8$ to model the lattice relaxation. We include Moiré mini-Brillouin zones around the ${\bf K}$-point of the monolayer graphene up to a cut-off distance of $\sim 5 |{\bf G}^M|$.

Calculation of Hofstadter butterfly
-----------------------------------

To calculate the effect of a finite magnetic field, we performed the Peierls substitution ${\bf k} \rightarrow {\bf k} - e {\bf A}$. The intralayer Hamiltonian Eq.  is diagonalized in a standard Landau level basis $| L \alpha n y \rangle$, where $L$ is the layer, $\alpha$ the sublattice, $n$ the Landau level index and $y$ the guiding center coordinate. We then write the interlayer coupling $H_M$ in this basis following Refs. [@bistritzer_moire_2011; @zhang_landau_2019; @hejazi_landau_2019; @lian_landau_2020]. The spectrum is calculated for flux per Moiré unit cell $\Phi = \Phi_0 \frac{p}{q}$ with $p,q$ coprime and $q\leq 30$. Note that the spectrum is the same for both spin and valley species, so we only compute it here for one valley, neglecting the Zeeman spin splitting.

The resulting spectrum, known as the Hofstadter butterfly[@hofstadter_energy_1976], turned out to match the experiments best for twist angle $\theta = 1.8^\circ$, see Fig. \[Fig:SupplTheory\]. In the literature (see Refs. ), there is a large variety of shapes of butterflies for the same twist angle. Further experimental study of the exact spectrum might elucidate which model parameters are best suited to describe tBLG. In order to fit the data, we have rescaled the energy in units of the bandwidth. The full Hofstadter butterfly is shown in Fig. \[Fig:SupplTheory\].

![ **The Hofstadter butterfly calculated for $\theta = 1.8^\circ$ and further parameters outlined in the text.** Each state pictured is fourfold (spin and valley) degenerate.[]{data-label="Fig:SupplTheory"}](FigS_FullButterfly.pdf){height="6cm"}

The dominant features of the Hofstadter butterfly are the band that evolves from the zeroth Landau level at small magnetic field, with Chern number +2 per spin/valley, flanked by two $C=-1$ bands to give a net zero Chern number for the fully filled band. The bandwidth of the $C=-1$ band at small fluxes can be approximated by the simple formula $W_{-1}(\Phi) = W_0 - a v_F^* \sqrt{\Phi} - b \Phi / m^*$, where $a,b$ are scaling constants, $v_F^*$ is the effective Dirac velocity at charge neutrality and $m^*$ is the effective mass at the flat band edges. Because of this scaling, the net density of states in the $C=-1$ band is increasing with flux, as can be seen in Fig. \[Fig:SupplTheoryDOS\]. Consequently, for a fixed interaction strength one expects a first order Stoner transition into a spin/valley polarized state as a function of $B$, if one does not already exist at $B=0$.

![ The average density of states of the $C=-1$ band is given by the number of states $n(\Phi) = 1 - \Phi$ divided by its bandwidth $W(\Phi)$. Since the bandwidth decreases faster than the number of states, the net density of states is an increasing function of flux. The continuum approximation here is given by taking the bandwidth $W_{-1}(\Phi) = W_0 - a v_F^* \sqrt{\Phi} - b \Phi / m^*$. \[Fig:SupplTheoryDOS\]](FigS12_DOS.pdf){height="4cm"}

### Disorder

In the absence of disorder, Landau levels are perfectly degenerate and are therefore always susceptible to quantum Hall ferromagnetism. Disorder, however, will broaden the Landau levels. Any realistic assessment of possible symmetry breaking in Landau levels requires a study of the role of disorder.

In tBLG, the dominant form of disorder is twist angle inhomogeneity [@wilson_disorder_2020]. This mostly affects the band edges, as the bandwidth of the flat bands is highly sensitive to twist angle. Here, we model the disorder phenomenologically by add a self-energy that is increasing with energy away from charge neutrality, $\mathrm{Im} \Sigma (\epsilon) = - \gamma_B | \epsilon|$, where $\gamma_B$ is a dimensionless parameter proportional to the strength of disorder. In the final calculations we set $\gamma = 0.001$. Note that this shape of disorder self-energy also applies to single-layer graphene in the presence of short-range potential impurity scattering.

Using this disorder potential, we can convert the Hofstadter spectrum into a density of states $\rho (\Phi, \epsilon)$, which is used in Fig. 1 of the main text.

### Hartree-Fock description of interaction effects {#Sec:HF}

To account for the state with spontaneous polarization, we performed a simplified Hartree-Fock calculation for each flux.As free parameters we introduce separate chemical potentials $\mu_{\sigma \xi}$ for each spin $\sigma = \uparrow, \downarrow$ and valley $\xi = \pm1$. The filling per spin/valley is $$\nu_{\sigma \xi} = \int_0^{\mu_{\sigma \xi}} d\epsilon \; \rho(\epsilon),$$ where $\rho(\epsilon)$ is the density of states obtained from the Hofstadter butterfly described above. The total kinetic energy is given by $$K = \sum_{\sigma \xi} \int_0^{\mu_{\sigma \xi}} d\epsilon \; \epsilon \; \rho(\epsilon).
    \label{HFKinetic}$$ The total interaction energy is a Hartree-Fock approximation of a generic repulsion between the spin and valleys[@zondiner_cascade_2020], $$V = U \sum_{(\sigma \xi) \neq (\sigma' \xi') } \nu_{\sigma \xi} \nu _{\sigma' \xi'}
    \label{HFInteraction}$$ where $U$ denotes the strength of the interactions. We minimized the total energy $K+V$ by varying the chemical potentials $\mu_{\sigma \xi}$, with the constraint that the total filling $\sum_{\sigma \xi} \nu_{\sigma \xi} = \nu$ is constant. We repeated this calculation for a range of densities $\nu \in [-4,4]$ and fluxes $\Phi \in [0,1/2]$.

In Fig. \[Fig:SupplTheory2\] we show the resulting filling per spin/valley as a function of total filling, for three values of the flux: $\Phi = 2/5$, $\Phi = 1/5$ and $\Phi = 1/30$, and $U=0.41W$. In the absence of spontaneous polarization, each spin/valley would be equally occupied. This always occurs when the total filling is $\nu = \pm 4 \Phi$, when all spin/valley degrees of freedom have fully occupied or unoccupied the $C=-1$ Hofstadter miniband.

For small fluxes (Fig. \[Fig:SupplTheory2\]b and c), the spin/valley polarization follows the same pattern as proposed by Ref. [@zondiner_cascade_2020]. Starting at $\nu=-4$, all four spin/valley flavors increase their filling in parallel, when a first order transition occurs that fixes one of the spin/valley flavors to a filling $\nu = -\Phi$, and at the same time ’resets’ the other three fillings to $\nu = -1$. This process of parallel filling followed by a reset continues until the $C=-1$ band is completely filled.

In contrast, at the higher flux $\Phi=2/5$ the filling of the $C=-1$ band occurs sequentially: first one spin/valley flavor completely fills the $C=-1$ band, followed by the second flavor, and so forth. This is a familiar pattern of ferromagnetic states, which occurs, e.g. in the interaction-split zeroth Landau level in graphene.

Note that this sequential filling at higher fluxes is consistent with the experimental results. For example, in between the (-3,-1) and (-4,0) state, one finds a (-2,-1) state. This state can be understood by having three spin/valley degrees of freedom completely filling the $C=-1$ band, whereas one spin-valley flaver partially fills the $C=-1$ band up to the $C=+1$ gap (see Fig. \[fig:3\])

![ \[Fig:SupplTheory2\] Filling per spin/valley as a function of total filling, obtained using Hartree-Fock calculations with $U=0.41W$ for three typical values of the magnetic flux $\Phi = 2/5$, $1/5$ and $1/30$. At low fluxes, the ](FigS11_Phi2_5.pdf "fig:"){width="30.00000%"} ![ \[Fig:SupplTheory2\] Filling per spin/valley as a function of total filling, obtained using Hartree-Fock calculations with $U=0.41W$ for three typical values of the magnetic flux $\Phi = 2/5$, $1/5$ and $1/30$. At low fluxes, the ](FigS11_Phi1_5.pdf "fig:"){width="30.00000%"} ![ \[Fig:SupplTheory2\] Filling per spin/valley as a function of total filling, obtained using Hartree-Fock calculations with $U=0.41W$ for three typical values of the magnetic flux $\Phi = 2/5$, $1/5$ and $1/30$. At low fluxes, the ](FigS11_Phi1_30.pdf "fig:"){width="30.00000%"}
