Topological heavy fermions in magnetic field

The recently introduced topological heavy fermion model (THFM) provides a means for interpreting the low-energy electronic degrees of freedom of the magic angle twisted bilayer graphene as hybridization amidst highly dispersing topological conduction and weakly dispersing localized heavy fermions. In order to understand the Landau quantization of the ensuing electronic spectrum, a generalization of THFM to include the magnetic field B is desired, but currently missing. Here we provide a systematic derivation of the THFM in B and solve the resulting model to obtain the interacting Hofstadter spectra for single particle charged excitations. While naive minimal substitution within THFM fails to correctly account for the total number of magnetic subbands within the narrow band i.e., its total Chern number, our method—based on projecting the light and heavy fermions onto the irreducible representations of the magnetic translation group— reproduces the correct total Chern number. Analytical results presented here offer an intuitive understanding of the nature of the (strongly interacting) Hofstadter bands.

The recently introduced topological heavy fermion model (THFM) provides a means for interpreting the low-energy electronic degrees of freedom of the magic angle twisted bilayer graphene as hybridization amidst highly dispersing topological conduction and weakly dispersing localized heavy fermions. In order to understand the Landau quantization of the ensuing electronic spectrum, a generalization of THFM to include the magnetic field B is desired, but currently missing. Here we provide a systematic derivation of the THFM in B and solve the resulting model to obtain the interacting Hofstadter spectra for single particle charged excitations. While naive minimal substitution within THFM fails to correctly account for the total number of magnetic subbands within the narrow band i.e. its total Chern number, our method -based on projecting the light and heavy fermions onto the irreducible representations of the magnetic translation group-reproduces the correct total Chern number. Analytical results presented here offer an intuitive understanding of the nature of the (strongly interacting) Hofstadter bands.

I. INTRODUCTION
Since the discovery of the remarkable phase diagram of the magic angle twisted bilayer graphene (MATBG) [1, 2], substantial effort  has been devoted to understanding its rich physics.
The presence of topological flat bands within this system [52][53][54][55] provides a novel platform to study the interplay between strong electron correlations and band topology. The recently introduced topological heavy fermion model (THFM) for MATBG [56][57][58] bridges the contrary signatures of localized [59,60] and delocalized physics [61,62] reported via STM and transport measurements [63,64]. Within THFM the low energy electrons are viewed as a result of the hybridization between heavy p x ± ip y -like Wannier states, localized at the AA stacking sites, and topological conduction fermions, denoted by f and c respectively in analogy to heavy fermion materials [56]. Among its other features, THFM allows for an intuitive explanation of the charged excitation spectra at integer fillings hitherto obtained via strong coupling expansion of projected models [65,66].
The large moiré period of ∼13nm in MATBG has revealed a sequence of broken symmetry Chern insulators yielding a plethora of B induced phases at lower fluxes [63,[67][68][69][70][71][72][73][74][75][76][77] and has showcased, for the first time, reentrant correlated Hofstadter states at magnetic fields as low as 31T [78]. Thus it becomes important to better understand the interplay of correlations and band topology in the presence of a perpendicular B field. Theoretical studies have previously focused on noninteracting [79][80][81] and strong coupling [82][83][84] regimes. Although exact, each employed considerable numerical analysis, restricting a deeper physical understanding of the mechanism for Landau quantization.
In this manuscript, we show how one can understand the Landau quantization of the strong coupling spectra in terms of hybridization amidst Landau levels (LLs) of c fermions and hybrid Wannier states of f fermions. We find that only a particular number of f fermion momentum channels are allowed to hybridize to c fermion LLs, with coupling strength which decreases with increasing B and increasing LL index m. Moreover, through our analysis we can clearly understand the reason why a naive minimal coupling is unable to recover the correct total Chern number of the flat band. In the flat band limit of THFM, our framework allows for an exact solution including the dominant interactions and analytically explains the total Chern number. Although going away from the flat band limit requires numerics, given the simple structure of the our Hamiltonian, we are still able to compute the spectrum to unprecedentedly small fluxes and find it is captured well by the analytical solution in the flat band limit which can be taken all the way to B = 0 as shown in Figs.(1), (8), (9) at ν = 0, −1, −2 respectively. Moreover it provides means to intuitively understand the underlying energetics of the problem as we tune the coupling between the c and f modes from 0% to 100% as shown in Fig.(2).
The formulas we derive are general for any rational value of Φ Φ0 , with Φ being the flux through the unit cell and Φ 0 being the flux quantum hc e , but we focus our analysis on the 1 q flux sequence and low B where the results become particularly transparent. Our analysis as well unveils the physical nature of the anomalous low energy mode which is seen to be almost B-independent, also observed in previous numerics [82], as the anomalous zero-LL of a massless Dirac particle, a key ingredient of the topological heavy fermion picture of MATBG. Although this work deals directly with THFM, our methods apply more generally. 16 + γ 2 , are recovered in B → 0 limit of our theory. The value of parameters used are J = 18.27meV , U1 = 51.72meV , γ = −39.11meV , v * = 1.624eV.Å, v * = −4.483eV.Å and λ = 0.3792Lm. We set k1,2 = 0, although inconsequential as magnetic subbands are Landau levels in this regime and thus do not disperse.

II. HAMILTONIAN AND THE BASIS STATES
The THFM in momentum space is given by [56] H 0 = |k|<Λc τ s (1) Here Λ c is the momentum cutoff for c fermions while f s, whose bandwidth is negligibly small, reside in the entire moire Brillouin zone. The tilde indicates the fermions describing the B = 0 Hamiltonian, λ ≈ 0.38L m acts as the damping factor for c-f hybridization and L m is the moiré period(mBZ). τ = +1(−1) represents graphene valley K(K ) and s spin ↑, ↓. The c-c and c-f couplings As the c − f coupling is tuned to 35% of its full value, 2mmax + 3 of the f modes start coupling with the c Landau levels and the remaining 2q − (2mmax + 3) remain at U1/2. As a result of the level repulsion at the crossings, the Landau levels get pushed down in energy, and start gaining f -character as the B field increases. The level repulsion keeps pushing the levels apart as c − f coupling grows further. (c) At full strength of the c − f coupling, the strong coupling window for narrow bands from J 2 to U 1 2 , is well gapped from the remote magnetic subbands emerging out of B → 0 energies E+ and −E−. The f nature of magnetic subbands within the narrow bands grows with field except for the J/2 level, which is a pure anomalous c-Landau level that does not hybridize with the f . are H c,τ aa (k) = where the parameters v * , v * , γ and M are given in [56]. The interactions effects in THFM are captured via the parameters J and U 1 [56,57]. In order to illustrate our finite B solution, we first focus on the charge neutrality point (CNP) wherein for valley polarised (VP) states, the mean-field interaction reads The narrow band filling factors ν = ±1 and ±2 are discussed in the section III C. The real space basis for c's is composed of four k · p Bloch statesΨ kaτ (r) at Γ point in mBZ. The localized f electron Wannier functions are denoted by W R,bτ (r); they behave as p x ± ip y orbitals sitting on the AA stacking sites spanned by moiré triangular lattice vector R, and W R,bτ (r) = W 0,bτ (r−R). The c-c and c-f couplings in Eqs.
Having assembled the low energy basis at B = 0, we now expand the fields for a given spin as where Ψ aτ (r) = √ N A ucΨΓaτ (r), A uc denotes the moiré unit cell area. Anticipating the appearance of anomalous Dirac LLs for the topological c fermions, we allow for the a dependence of the upper cutoff on the LL index at each valley τ , denoted by m a,τ , with m 1,+1 = m 2,−1 = m max + 1, m 2,+1 = m 1,−1 = m max , m 3,+1 = m 4,−1 = m max + 2 and m 4,+1 = m 3,−1 = m max − 1.
By construction, Ψ aτ (r)χ krm (r) and η bτ kr (r) are eigenstates oft L1 andt q L2 with eigenvalues e −2πik1 and e −2πiqk2 respectively. This guarantees states with different k to be orthogonal. Overlaps of states with the same k but different r, r is discussed in detail in the SM-C, though we note in passing that the orthogonality amidst c-fermion states relies on the fact that their overlaps are exponentially small in 2 g 2 as long as the LL index m is held below an upper cutoff m max q/2, where g is the reciprocal moiré lattice vector.
The single-particle THFM at B = 0 is then where the cc couplings are with matrix elements that are found to be the same as those obtained by the direct minimal substitution in (2) and expanding in LL basis, as is expected from k · p [87]: where the Pauli matrix σ x acts on the c orbitals and with h −1,c mm = −σ x h +1,c mm σ x . For M = 0, we recover the LLs of two massless Dirac particles, with two zero LLs at each valley. Because these modes play a key role in explaining the interacting spectra of THFM in B, we plot their and f 's spectrum with the term in Eq.(4) included but turn off the c-f coupling by setting the prefactor of H τ cf to zero in Fig.(2a). For each spin projection, there are two f s per moiré unit cell at the energy ± U1 2 at τ = ∓1 intersected by the LLs with strong B dependence due to c's fast v * . The total number of f states is B-independent while each LL contains N φ φ0 states. We expect that H τ cf hybridizes these modes within an energy ∼ γ, but in order to understand how such hybridization leaves behind a band of states whose total number is independent of B -because its total Chern number vanishes [88]-and which LLs may decouple from f s, we need to carefully analyze H τ cf . As shown below, it is not simply minimally substituted (6). Instead where, using the magnetic translation symmetry of H τ BM p − e c A and the free sum over the s, we find Because A acts on a well-localized function centered at r = 0, its contribution can be neglected at low B (we have confirmed this numerically in SM-D 3  [89] (see also SM-D 2).

III. ANALYSIS OF THE 1 q FLUX SEQUENCE
It is particularly revealing to analyze the case p = 1, i.e. the φ/φ 0 = 1/q sequence. Then r = 0, because, as FIG. 3. Schematic representation of the cf coupling for φ/φ0 = 1/q in Eq. (20). Each tick represents (the 1D projection of) a moire unit cell, illustrating the overlap between a 2D localized heavy state with size λ siting at the origin (red) and a Landau level (LL) i.e. a 1D harmonic oscillator (h.o.) shifted in the x-direction, its wavefunctions sketched by blue, with a plane wave phase variation in the y-direction (not shown) that depends on the shift. The r determines the momentum absorbed by the LL as well as the unit cell to which the h.o. is shifted and k2q fine tunes the shift within the unit cell. The index j then determines q-unitcell periodic revival of the h.o. The black parabolas mimic a quadratic potential to accompany the h.o. wavefunctions. discussed above, r ranges from 0 to p − 1. Making use of (11) in (19) and substituting the explicit expression for χ from (9), allows us to perform the sum over n. We find The above expression can be visualized as an overlap between a 2D localized heavy state with size λ sitting at the origin and a 1D h.o. shifted in the x-direction with a plane wave phase variation in the y-direction that depends on the shift (see Fig.(3)).
To understand under what choice of m, r , j is this integral significant, note that the h.o. wavefunction is localized in the x direction about (r + jq + k 2 q) L 1x , and its width is ∼ 2 √ 2m + 1 √ q unit cells. In addition, the combination r q + j + k 2 controls the period of oscillation in the ydirection set by 1/( r q +j+k 2 ) times the unit cell size. The integer r + jq thus determines the unit cell to which the h.o. is shifted, and, because k 2 2π 2 /L m = k 2 qL 1x , the value of k 2 q ∈ [0, 1) fine-tunes the shift within the unit cell. The index j then determines q-unit-cell periodic revival of the h.o. also illustrated. For example, if r = j = 0 then the h.o. is centered at the unit cell containing the localized heavy state and the period of oscillations in the y-direction is long compared to the unit cell, encompassing at least q unit cells. The hybridization with the localized heavy state proportional to γ will then be significant. Since the spatial extent of the h.o. state in the x-direction is ∼ 2 √ 2m + 1 √ q unit cells -which at low B is much longer than the localized heavy state even when we consider that it oscillates and has m-nodesthe result of the integration will be approximately given by the value of the h.o. wavefunction at −k 2 qL 1x , up to an overall phase, unless m is close to m max q/2. If we keep j = 0 but increase r to 1 then the h.o. is centered at the unit cell adjacent to the one containing the localized heavy state and the period of oscillations in the y-direction is still long, between q/2 and q unit cells. The hybridization with the localized heavy state proportional to γ will still be significant and the result of the integration will still be approximately given by the value of the h.o. wavefunction but now at −(1 + k 2 q)L 1x , up to an overall phase. We thus see that for fixed j = 0 and any value of k 2 q, increasing r past ∼ √ 2m + 1 √ q results in an exponential suppression of the hybridization integral due to the large off-set. Therefore, for the values of r past q/2, it is the j = −1 revival copy of the h.o. states which gives the dominant contribution. All other values of j can be neglected, allowing us to remove the sum over j in Eq. (20) and replace where sgn + (x) is the usual sign function except at 0 where it evaluates to 1. Increasing the value of m increases the size of the wavefunction which is another way to understand the effect of the bound on m. Here and The two variable Hermite polynomials [89] are given by where m denotes the floor function at m. Their relation to the Hermite polynomials used above is H m (x) = H m (2x, −1). In the limit λ → 0, the 2D heavy localized state becomes the Dirac δ-function, and as expected based on the above discussion, we indeed recover F m (λ → 0, x 0 ) = ϕ m (−x 0 ). We found that keeping the full form of F m is needed in order to achieve accurate results even for the low B range, therefore we do not take this limit when handling F m (see also Fig.(12) in SM-E). The exponential suppression factor multiplying F m in the Eq.(23) comes from the yintegration. Although the suppression depends on the shifted position of the h.o. states, its dependence on r is weaker than in F m which comes from the x-integration, as expected based on the discussion above. The off-diagonal hybridization can be expressed as Using the exact recursion relations between H m and its derivatives, I x and I y can be expressed in terms of I 0 as We add interactions at CNP via the mean field coefficients J and U 1 . This is accomplished by promoting the operatorsc kaτ andf kbτ in Eq.(4) to c aτ krm and f bτ kr respectively, with a sum over all finite B quantum numbers. The numerically determined strong coupling Hofstadter spectra using the above results for the matrix elements for τ = −1 is shown in Fig.(4a). To study the effects of including kinetic energy within the THFM, Fig(4b)shows a comparison with the Hofstadter excitation spectrum of the Coulomb interacting BM model with the band kinetic energy projected away. The gauge-invariant formalism introduced by Ref. [83] is used to compute the exact strong coupling excitation spectrum to fluxes reaching q = 60 using the sparse matrix properties of the basis (general and efficient formulae may be found in SM-J). The kinetic energy, neglected in the strong coupling calculations, decreases the charge gap but preserves the low-lying anomalous Landau level. Neglecting λ 2 2 in the prefactors of I 0 's in Eq. (29) is an approximation of same order as neglecting A in H τ BM in Eq. (19). We confirm it numerically to be an excellent approximation in SM-D 3; note that we do not set λ 2 2 to zero in I 0 for the reasons mentioned above. Thus the finite B c-f coupling can be accurately re-expressed as and h −1,cf m,r (k) = −σ x h +1,cf m,r (k)σ x , where σ x acts in orbital space of c and f fermions. The above form for c-f coupling sets up a platform to understand the problem in terms of coupled and decoupled modes of f , and eventually paves the way for an analytical solution.
A. Change of basis for the f -modes We will find it useful to first perform a singular value decomposition where U, V are are m a,τ ×m a,τ and q ×q unitary matrices respectively, and the m a,τ × q rectangular matrix Σ contains the singular values along the main diagonal and zeros elsewhere; summation convention on repeated indices is used throughout this section unless explicitly stated. Through U and V , we can gain insight into the form of f and c modes allowed to hybridize. For example, wherec † 1τ k0m = c † 1τ k0m U mm andf 1τ kr = Vr r f 1τ kr are the new modes. In the above, thef 1τ kr modes withr > m 1,τ decouple from c's (recall that m 1,+1 = m max + 1 and m 1,−1 = m max ). The remaining f modes couple to c modes with the strength γΣ m ,r = γΣ m , where Σ m denotes the m th singular value.
The columns of matrix U are the eigenvectors of the following matrix Since r is being summed over, we are allowed to shift its range to − q 2 , − q 2 + 1, . . . , q − 1 − q 2 and thus set j = 0 in Eq. (21). Hence we can set r q = r in Eq. (23) for I 0 m,r . Since F m decays exponentially at ± q 2 , in the low field limit, we can replace the bounds of this sum by ±∞. Using Eq.(23) we have Using the Dirac comb identity, r ∈Z δ(ρ − r ) = t∈Z e 2πitρ , we can convert this summation to an integral as Because 2πt/L 1x is a large wavevector compared to the length scale at which the h.o. wavefunctions vary, ∼ , i.e. L1x 2πt at low field, the overlaps are exponentially suppressed in ( 2πt L1x ) 2 . We can thus set t = 0 in the above sum. Moreover since the exponential factor contains λ 2 2 , the off-diagonal terms can be neglected and thus Λ 0 mm is diagonal up to a good approximation. This implies that the matrix U is close to identity and the low B singular values can be given as Note that the information of magnetic quantum number k gets dissolved automatically due to the low field limit where the magnetic sub-bands are not k dispersive. A compact expression for Σ m can be computed from Eq.(36) using two index Hermite polynomials [89](also see SM-D 2) where κ 2 = λ 2 2 = ( φ φ0 )2πλ 2 /A u.c and In the B → 0 limit, upto first order in flux, the singular values can be approximated as Intuitively this shows that Σ m dependence on m grows linearly with m even at negligibly small fluxes. The comparison of analytical singular values given by Eq. (38) with the ones obtained numerically from Eq.(23) is  Having the closed form for the singular values, we hereby reformulate the finite B Hamiltonian (Eq.(13)) in terms of the modes that are allowed to hybridize. The analysis for τ = +1 is presented below, for τ = −1 in SM-F 1. Using the SVD decomposition discussed we can re-write the matrix elements as where we used the fact that U is an identity matrix and the rectangular matrix Σ m,r is non-zero only along its main diagonal.
Because m 1,+1 −1 = m 2,+1 = m max , the upper bound on r summation is the same in Eq. (44) and Eq. (45). Thus for a given k, out of the available 2q f −modes, only 2m max + 3 couple; recall that m max q/2. The meanfield interactions for f − fermion modes in newf basis can be given as where1(2) = 2(1). Note that there are 2q − (2m max + 3) decoupled f modes for each k. Physically this corresponds to 2 − (2m max + 3)/q states per unit cell.
The coupled modes can then be described as where m α=1,...,4 = m α,+1 , m 5 = m max + 1 and m 6 = m max , and We now define an operator whereâ is a simple h.o. lowering operator in terms of which matrix Ξ mα,m α can be expressed as Here |m is a simple h.o. eigenstate and Σ(m) = Σ m . Note that a naive minimal substitution into Eq.(1), with λ = 0 would reproduce Eq.(51), except the singular values would get replaced by 1, which corresponds to the B → 0 limit of our analysis. The decoupling of 2 − (2m max + 3)/q modes per moiré unit cell per spin and the fairly strong B and m dependence of singular values thus gets completely overlooked by the naive minimal substitution, which therefore fails to recover the correct state count for narrow bands to be 2 states per moiré unit cell per spin per valley. Although we get the correct total Chern number 0, i.e. total 2 states per moiré unit cell per spin independent of B in the narrow band strong coupling window set by J 2 , U1 2 , we can explain the state counting analytically by presenting an exact solution to the operatorĥ α,α in the flat band limit M = 0.
The rest of the problem can be solved using the following ansätze: 3 |1 , 0, c 3 |2 , 0, c where m ∈ {2, . . . , m max +1}. c The ansätze θ 3 and θ 5 yield the 3 × 3 and 5 × 5 Hermitian matrices, whose eigenvectors are c (3) α and c (5) α , respectively: (58) Similarly, the ansatz θ m 6 yields the following 6 × 6 Hermitian matrix for each m, whose eigenvectors are c The strong coupling spectrum obtained using the above matrices is shown in Fig.(1), which also includes the spectrum at τ = −1. One way to obtain the spectrum at τ = −1 is by replacing J, U → −J, −U in the above decoupled matrices(see SM-F 1). Also note that the decoupled f −modes and the anomalous c−mode at τ = −1 then form the field independent U1 2 and J 2 levels respectively.
The state counting within the spectrum can now be easily understood from the aforementioned analysis. As shown in Fig.(6), the B → 0 energies for J = U 1 = 0 are 0 for the flat band, and ±γ for the remote bands, because Σ m → 1. Note that these are the non-interacting zerofield energies at Γ point in mBZ. For any non-zero B, the total zero mode count can be understood as follows: h τ,m 6 being bipartite contributes two zero modes (even if we do not set Σ m = 1). Including all m max such matrices gives 2m max zero modes. h τ 3 , h τ 5 and the anomalous c level contribute one zero mode each, giving us a total of 2m max + 3 zero modes. Including the 2q − (2m max + 3) of the decoupled f modes, at each k we recover sum total of 2q zero modes in the non-interacting case for each valley, giving the total of 2 states per moiré unit cell per spin independent of B, i.e. the total Chern number 0. Note that the magnetic subbands within the narrow band window must remain separated by a gap from the remote subbands for a non-zero range of B because the remote bands emanate out of ±γ at B = 0.
The above analysis can straightforwardly be applied to other integer fillings. In this section we illustrate it by extending the formalism to nonzero fillings ν = ±1, ±2.
Since the mean field interactions are spin dependent for integer fillings, we introduce the spin quantum number label, s =↑↓, to the c and f fermion annihilation operators as c aτ krms and f bτ kr s in order to differentiate between the spin-up and spin-down sectors. of decoupled c mode. Including the q − (mmax + 1) = q − 6 decoupled f modes with energy − 1 2 U1 − 6U2 for the Chern −1 band and q − (mmax + 2) = q − 7 decoupled f modes with energy − 3 2 U1 − 6U2 for the Chern +1 band, we have in total (q − 6) + 5 = q − 1 and (q − 7) + 8 = q + 1 magnetic subbands for Chern ∓1 bands respectively as we should.
For ν = −1, we discuss the sector valley K spin ↓ for the VP state [56], wherein at zero-field, the charge ±1 excitations occupy well gapped Chern ∓1 states marked by red in Fig(8a). The other spin-valley sectors at ν = −1, namely valley K spin ↑ and valley K spin ↑↓(degenerate) are discussed in SM-F 2,G. The red band in the energy window (−30, −50) meV has Chern number −1 [56], while the dispersive red band in energy window (−55, −100)meV has Chern number +1 [56]. Our formalism straightforwardly shows that in presence of magnetic field B, the Chern ±1 bands gain and lose states respectively [88], with a total state count pinned to 2 per moiré unit cell. The spectrum at ν = +1 is related by particle-hole symmetry. The interactions at valley K and spin ↓ read as Here W a∈{1...4} and U 2 are mean field coefficients with W 1 = W 2 and W 3 = W 4 [56]. In thef basis we can re-write the interaction for f -fermion modes as where1(2) = 2, 1. Note that out of the available 2q f modes, q − (m max + 2) are decoupled with energy ν 3 2 U 1 + 6U 2 and q − (m max + 1) are decoupled with energy ν 1 2 U 1 + 6U 2 , i.e. a total of 2 − (2m max + 3)/q decoupled f modes per moiré unit cell. The coupled modes can then be described by where m α=1,...,4 = m α,+1 , m 5 = m max + 1 and m 6 = m max , and where the operatorsĥ +1,↓,ν=±1 α,α are given aŝ The eigenstates of the above operator are exactly solvable for the flat band limit, i.e. M = 0 and the spectra is shown in Fig.(8b). The decoupled c fermion given in Eq.(53) forms the field independent ν(W 3 + J 2 ) level. The remaining eigenstates can be obtained using the ansätze given in Eqs.(54)- (56). Setting up eigenequation for these ansätze yield us the following 3 × 3, 5 × 5 and m max 6 × 6 matrices: h +1,m,ν=±1 where m ∈ {2 . . . m max + 1}. The magnetic subbands contributed by the coupled modes for Chern +1 band emanate out of B → 0 energy ν(W 3 + J 2 ) of the above decoupled matrices, which is m max fold degenerate for matrix Eq.(71) and singly degenerate for matrices Eq.(69) and Eq.(70). Including the field independent level ν(W 3 + J 2 ) of the decoupled c mode, we have in total m max + 3 magnetic subbands emanating out of ν(W 3 + J 2 ) for Chern +1 band. Now recall that we have q − (m max + 2) decoupled f modes with energy ν( 3U1 2 + 6U 2 ) for the Chern +1 band, which gives in total q + 1 magnetic subbands for the Chern +1 band. Similarly the coupled modes contribute m max magnetic subbands for Chern −1 band which emanate out of the m max fold degenerate B → 0 energy ν(W 3 − J 2 ) of matrix Eq. (71). Including the q −(m max +1) decoupled f modes with energy ν( U1 2 + 6U 2 ) for the Chern −1 band, we have in total q − 1 magnetic subbands for the Chern −1 band. An illustration for m max = 5 is shown in Fig.(7). Our method thus captures the fact that Landau quantisation of Chern ±1 band offers in total q ± 1 magnetic subbands or 1 ± 1 q = 1 ± φ φ0 states per moiré unit cell as expected. The full spectrum in flat band limit is shown in Fig.(8). Note that the B → 0 energy of the decoupled matrices correspond to the zero field THFM energies at Γ in mBZ.

IV. DISCUSSION
We have put forward a generalization of THFM in finite B. Although the formalism applies to any rational value of Φ Φ0 , the physical nature of hybridization amidst the heavy f and topological c fermions is particularly revealing for the 1 q sequence. In the Landau gauge, the momentum imparted by the heavy hybrid Wannier states onto the c's LLs is converted into the shift of the position of the LL wavefunctions, each unit of 1 q of the primitive reciprocal lattice momentum g 2 causes a shift by a moiré unit cell. Because the momentum is defined only modulo g 2 , the LL wavefunction shift is also defined only modulo q unit cells, causing "revivals" illustrated in Fig.3. Together with a change of the oscillation rate of the wavefunction in the perpendicular direction under the shift, this process controls the strength of the cf hybridization along the 1 q sequence. Its consequence on the spectrum, including interactions at CNP, are illustrated in the Fig.(2) as the c-f hybridization is turned on from zero to its full strength. The finite B analytical solution for the flat U (4) symmetric THFM provides an intuitive picture of the mechanism for Landau quantization of the strong coupling spectra of MATBG at integer fillings in terms of the decoupled f modes and coupled c-f modes, all the way to zero magnetic field. It also provides a deeper understanding of the nature of the ± J 2 level at CNP, observed in numerics before [82], as the anomalous zero-LL of a massless Dirac particle, a key ingredient of the topological heavy fermion picture of MATBG. For q dependent m max , we see the a peeling away of levels close to U1 2 around q ≈ 80, because the theory starts "losing" LLs with increasing field, reminiscent of the results obtained in [82]. Although there are 2 − 2mmax+3 q decoupled f −modes per unit cell per spin at CNP, the total number of states in the narrow band strong coupling window still remains pinned to 2 per unit cell per spin, as expected for a total Chern number 0. Even though the full M = 0 problem requires numerical analysis, we are able to probe till fluxes at least as low as 1/145, which was not possible through the the framework of strong coupling expansion. We moreover argue that the overall physical picture should stay unchanged, as M anyways is the smallest energy scale in the problem. Throughout the text we neglected the spin Zeeman effect as it leads to a much smaller energy splitting than the orbital effect, the former is only a few Kelvin at the highest fields considered here while the latter is several meV, so at least an order of magnitude larger. The effect of renormalization of mean field parameters in magnetic field is yet to be incorporated in our framework. A full analysis for . The color on panel-(a) labels the spin sector, while color labelling for panels-(b,c) denotes the f-character.
other integer fillings, candidate groundstates [90,91], and Hofstadter-scale fluxes where reentrant many-body and topological effects are at play [49,78,[92][93][94][95][96][97][98][99], is also left for the future work.   The Bistritzer-Macdonald (BM) model of magic-angle twisted bilayer graphene [85] in the valley K(τ = +1) is where σ acts in sublattice space and p = 0 denotes the moiré Dirac point K M . The interlayer hopping functions are where where I n is the n × n unit matrix. We will present results for the case w 0 /w 1 = 0.7. The BM model is invariant under translation by integer multiples of moiré lattice vectors L 1 = L m ( whereT L1,2 are the usual discrete translation operators.

Appendix B: Magnetic Translation Operators
At B = 0 in the Landau gauge A = (0, Bx), the BM model in flux is found via minimally substitution: The finite B BM model at K (τ = −1) can be obtained similarly by applying time reversal to Eq.(A1) followed by the minimal substitution. However at B = 0, H τ BM is still invariant under the translation by L 2 , but a translation by L 1 needs to be accompanied by a gauge transformation as In other words iff (r) is an eigenstate of H τ BM p x , p y − eB c x , then so is L1 is the generator of magnetic translation by L 1 . Note thatt L1 can alternatively be presented ast where the magnetic translation wavevector q φ is defined as where g 1,2 are corresponding reciprocal lattice vectors to L 1,2 and are given as: where L 2 =| L 2 |= L m . Magnetic translation by L 2 are generated byt L2f (r) =f (r − L 2 ). We thus havê where φ 0 = hc e and the flux through the unit cell is φ = BL 1x L 2 . If φ/φ 0 = p/q, with p and q relatively prime integers, The simultaneous eigenstates of the magnetic translation operatorst L1 andt q L2 can thus be used to produce a complete and orthonormal set of basis states for solving the B = 0 BM model. As we will see below, depending on whether they are conduction or heavy fermions there will be p states or q states per momentum point. It will be helpful to derive an identity fort s L1 which relates it toT s L1 .
where k 2 ∈ [0, 1). These hybrid Wannier states although localized along the x direction are Bloch-like extended in the y direction and thus eigenstates oft L2 t L2 w bτ (r, k 2 g 2 ) = w bτ (r − L 2 , k 2 g 2 ) = e −2πik2 w bτ (r, k 2 g 2 ). (C2) We next construct normalized eigenstates of MTG out of the hybrid Wannier states as where k 1 ∈ [0, 1) and normalization N = n tot s tot . Here s tot and n tot denote the total count of s and n. Note that although these summations should be unbounded, we require these cutoffs for defining the normalization of our states and intermediate steps. They eventually cancel in the final formulas are are not physical. The normalisation is justified in Eqs.(C7)-(C9). The domains of k 1,2 can be understood via noting that η bτ k1k2 is periodic under k 1,2 → k 1,2 + 1. Under magnetic translationst L1 andt q L2 , we then havê Thus states with k 2 that differ by 1/q have the same eigenvalue undert q L2 . This is because the qL 2 translations break up the k 2 domains into units of width 1/q. To make it apparent from the quantum numbers, we relabel the states as Thus η bτ kr (r) with different k quantum numbers are guaranteed to be orthogonal. Moreover we have q states for given k, corresponding to r ∈ {0 . . . q − 1} for each b ∈ {1, 2}, i.e. a total of 2q heavy fermion states for given k. For same k and different r the overlaps for given valley arê Now since the B = 0 Wannier states are well localised and orthonormal, we havê and thus,ˆd Henceforth we have a complete orthonormal basis for the heavy fermions in finite B, constituted by 2q states for given k, i.e. 2 states per moiré unit cell per valley per spin, which clearly implies that these are total Chern 0 states. Now let us discuss the basis for c fermions. Remember that at B = 0, the basis for c fermions is constituted by four k · p Bloch states at the Γ point in the moiré Brillouin Zone(mBZ), per valley per spin. We denote them byΨ Γaτ for a ∈ {1 . . . 4} and τ = ±1.
We promote the B = 0 c fermion basis to finite B using a result obtained when the k · p method is extended to finite B [87], which prescribes the ansatz for finite field k · p states to be Landau level(LL) coefficients on top of the zero field k · p states. Henceforth, we define the finite B basis for k · p Bloch states as an expansion over product of LLs and the B = 0 k · p Bloch states at the Γ point in mBZ. However in order to use the same quantum numbers label as finite B f fermion basis we first project the Landau gauge LL onto the representation of MTG as where k 1 ∈ [0, 1) and unlike before, k 2 ∈ [0, p/q) as justified in Eqs.(C13)-(C16). The functions ϕ m are mth 1D harmonic oscillator states that determine the Landau gauge Landau levels where H m is the Hermite polynomial and 2 = c/(eB). The normalization of the LL MTG states is explained in Eqs.(C23)-(C28). The domains for the quantum numbers k 1,2 can be understood by noting that χ k1k2m is periodic under k 1 → k 1 + 1 and up to a phase under k 2 → k 2 + p/q as shown below now using the fact, 2π 2 L2 p q = L 1x we have The domain for k 2 can also be understood in the context of the fact that χ k1k2m are essentially LLs and thus behave as Chern number +1 states in B = 0. In order to illustrate it, let us discuss the total number of available states for χ k1k2m . We start by noting that under magnetic translationst L1 andt q If we have a system size N 1 L 1 and qN 2 L 2 , then So, the total number of states is N 1 Here A uc =ẑ · (L 1 × L 2 ) denotes the area of moiré unit cell. Thus qA uc is the area of magnetic unit cell and total area A tot = N 1 N 2 qA uc . This is a well known result for the degeneracy of the Landau level. Let us now discuss the orthonormalization for the LL MTG eigenstates. The overlaps amidst the MTG LL are given as: where in Eq.(C24)-(C25), we have used the fact that y-integral implies k 2 −k 2 = (s −s) p q . Now since k 2 , k 2 ∈ [0, p q ), the integral evaluates to n tot L 2 δ k2,k 2 δ s,s . Moreover each harmonic oscillator function in the x-integral can be shifted by the same amount since the arguments coincide. Therefore the x-integral evaluates to δ mm . Eventually in Eq.(C26)-(C27), the s summation implies k 1 − k 1 = integer. But since k 1 , k 1 ∈ [0, 1), we must have k 1 = k 1 and the summation evaluates to s tot δ k1,k 1 . Therefore χ's are orthogonal and normalized. To have the same quantum number label as that of B = 0 f fermion states, we relabel the MTG LLs as Thus for given k, we have p states, corresponding to r ∈ {0 . . . p − 1}, i.e. p q states per moiré unit cell per valley per spin as expected. Finally we have the finite B c fermion basis states as Ψ aτ χ krm , where with N being the total number of moiré unit cells. The √ N A uc factor is required due to normalization of χ. Note that since Ψ aτ is a Bloch state at Γ,T where i ∈ {1, 2}. Thus Ψ aτ χ krm are MTG eigenstates with same eigenvalue as that of χ krm . Now let us discuss the orthonormalization for these states. Since in this work we focus on 1 q sequence, for the following discussion on orthonormality we set index r in χ krm to zero and work with χ k0m . The Bloch periodicity for Ψ Γaτ (r) allows us to expand it as where G = n 1 g 1 + n 2 g 2 , n 1,2 ∈ Z, denote moiré reciprocal lattice vectors and can be found in [57]. Now using Eq.(C30), the overlap matrix can be given as where 1/q is the flux per unit cell per flux quantum Φ 0 and m max denotes the upper cutoff on the indices m and m . The integral can be computed as: The delta function constraint tells us where g 2 = g y L2 2π ∈ Z. Note that getting δ k 1(2) ,k 1 (2 ) from the overlap is a direct consequence of Eqs.(C17)-(C18) for 1 q sequence. The above integral can be evaluated as where p x is momentum operator along x direction later changed to harmonic oscillator(h.o.) basis: x = √ 2 (a + a † ), and p x = i We also know the identity where L m−n n (x) is the associated Laguerre polynomial,

c-c Coupling
The c − c coupling is given as where M +1 µ = v F I 2 ⊗ σ, M −1 µ = −v F I 2 ⊗σ with σ = (σ x , σ y ) andσ = (σ x , −σ y ). T τ denotes the remaining factors in BM Hamiltonian. So, where the matrix ε τ is given in Eq.(D15) is k independent because Ψ aτ is defined at Γ. So, Since the factors ε τ a a Ψ * aτ (r)Ψ a τ (r) and Ψ * aτ (r)M τ µ Ψ a τ (r) are periodic with respect to primitive moiré lattice vectors L 1 and L 2 , we can perform a Fourier expansion as where The matrices τ a,a and M τ,µ aa can be obtained from the zero field cc coupling matrix where the Pauli matrices σ act in the orbital space of c fermions. Since t L1,2 , p µ − e c A µ (r) = 0, we have

27
The action of the two components of p µ − e c A µ (r) in Landau gauge can be calculated as Substituting back to χ, and using the orthogonality of χ's and Eq.(D15), we finally have, Note that a straightforward canonical substitution in Eq.(D14), k x + ik y → √ 2Bâ, would yield us the same finite field c − c coupling, whereâ denotes the LL lowering operator for the corresponding c fermion LLs.

c-f Coupling
The c − f coupling at finite B is given as Note that the B-field term eventually turns out negligible because it acts on a well localized function and so at small B it is exponentially suppressed in the region where the vector potential is appreciable.
In order to proceed we use the fact that [56] d 2 re −ik ·r Ψ * aτ (r)H τ BM (p x , p y ) W 0bτ (r) ≈ A uc e −k 2 λ 2 /2 H cf,τ ab (k ), which is √ N A uc bigger than in [56] because of choice of normalisation of χ discussed in previous section. Using we have Moreover since H τ BM (p) is linear in p, we can approximatê where M τ ab can be read of from Eq.(D37) to be Thus we have Therefore, L 2 e 2πis p q y L 2 e 2πi(k2+ r q ) L 2 e −πis(s−1) p q L 1y L 2 e 2πi(k2+ r q +s p where the sum over n above led to the Diophantine equation Plugging in the normalization factors results in Note that the y integrals are elementary. At small B and not too large m, the x integrals we be performed by taking the advantage of the mismatch in length scales i.e. λ and Taylor expand the harmonic oscillator functions about the origin. However at larger B we may need to obtain the integrals with higher accuracy. To this end we will find it useful to employ two-variable Hermite polynomial H n (x, y) result from [89]. It will be useful to first discuss the integrals required for casting Eq.(D56) into a closed form.
We start by defining the two-variable Hermite polynomial as Let y = x/ and y 0 = x 0 / . Then, as expected because in the limit λ → 0, the harmonic oscillator wavefunction does not vary significantly within the spatial extent of the heavy fermion gaussian, which is set by λ. Now, let The y integrals in Eq.(D56) can be performed usinĝ Using the integrals provided in Eq.(D66) and Eq.(D68), we now proceed to evaluate Eq.(D56). L 2 e πis(s−1) p q L 1y L 2 L 2 e πis(s−1) p q L 1y L 2 e −2π 2 (k2+ r q +s p q ) L 2 e πis(s−1) p q L 1y L 2 L 2 e πis(s−1) p q L 1y L 2 The minimal coupling term, eM u A µ /c, can be expressed in terms of the non-derivative coupling I 0 [mr],[r ] (k 1 , k 2 ) and k y coupling I y [mr],[r ] (k 1 , k 2 ) as shown in this section. Let us call this term I A , given as L 2 e πis(s−1) p q L 1y L 2 We can re-express the factor xϕ(x − x 0 ) as Moreover, using recursion relation for Hermite polynomials Substituting Eq.(D81) into Eq.(D78), we have: Using the fact 2 = qL1xL2 2πp , note that Thus we have 3. Low field Analysis of c-f coupling For analyzing the low B structure for c-f hybridization matrix, given in Eq.(D86), we will find it useful to first discuss recursion relations for F m (λ, x 0 ) defined in Eq.(D66). We start by considering the generating function for H m (x, y) (defined in Eq.(D57)) Using the generating function, we can derive the following recursion relations L 2 e πis(s−1) p q L 1y L 2 L 2 e πis(s−1) p q L 1y . Now using the fact 2πp 2 qL1xL2 = 1 and the recursion relation given in Eq.(D93), we have where I 0 [mr],[r ] is defined in Eq.(D70). Using the above derived form for I y [mr], [r ] in Eq(D97), the minimal coupling term, given in Eq.(D84), can be rewritten as Note that this is an O( λ 2 2 ) term. At low B, i.e λ 2 2 , this term is negligible as anticipated. This is confirmed numerically as shown in Fig.(11a). Note that we however still need to keep the full λ 2 2 dependence within F m for accuracy in low B range as shown discussed in Sec.(E)- Fig.(12). Keeping this in mind, let us now dwell into a low B limit wherein we can drop the O( λ 2 2 ) term in the recursion relation Eq.(D93), i.e.
Moreover using Eq.(D97), I y [mr], [r ] in this limit reads    In this section, we study the singular values of I 0 , whose strength controls the coupling between the heavy and conduction fermions in finite B. The singular values of I 0 [mr],[r ] for the 1 q sequence can be approximated via the integral F m (λ, ρ) defined in D66 can re-expressed as: Naively one would expect that dropping O( λ 2 2 ) in F m (λ, x 0 ), i.e. replacing F m (λ, x 0 ) by ϕ(−x 0 ), to be good enough approximation at low field. However as discussed in main text, the fairly strong m dependence of singular values gets compromised under such an approximation, as shown in Fig.(12a).
Let ρ/ = x and κ = λ/ . Then Now, using the result from [89], we havê and the two variable Hermite polynomial H n (x, y) is defined in Eq.(D57), we have where ξ(κ) is given as Note that where A uc is moire unit cell area. Since Σ 2 m is an analytic function in κ 2 , we have an expression for the singular values continuously down to zero flux. Also note that as φ i.e. approaching 1 at zero flux with a negative slope slope of m + 1 2 2πλ 2 Auc as shown in Fig.(13c). Note that even at extremely low flux, the m dependence of singular values grow linearly.
It would prove feasible to discuss the asymptote of Σ 2 m for large values of m. Let us start expanding the four variable Hermite polynomial in Eq.(E6).
where for brevity we denote κ 6 ξ(κ) and 2 ξ(κ) by y and β respectively. Then where we divided the sum over k into k = 0, k = 1 . . . Using Eq.(E16), we have Now considerĀ where in Eq.(E19), we use the Stirling's approximation for large m, i.e. ln(m!) = mln(m) − m. Similarly consider where in Eq.(E23) we again have used the Stirling's approximation. Now substituting back Eq.(E20) and Eq.(E25) in Eq.(E17) we have the following asymptote: • If m is even (E26) • If m is odd (E27) Fig.(13b) shows that the asymptotic expression for singular values, given in Eq.(E26) and Eq.(E26) , exactly matches the analytical singular values, given in Eq.(E8), for all m = 0 (since the asymptotic formulation does not hold for m = 0). In the Fig.(13b), m = 0 singular value on the asymptotic branch has been calculated using the analytical value, i.e. Σ 0 = 1 √ ξ(κ) . Appendix F: Discussion for τ = −1 In this section we discuss finite B THFM at τ = −1. Note that summation convention on repeated indices is implied unless explicitly stated. Using the SVD decomposition discussed we can re-write the matrix elements as where we used the fact that U is an identity matrix and the rectangular matrix Σ m,r is non-zero only along its main diagonal. Similarly Because m 2,−1 − 1 = m 1,−1 = m max , the upper bound on ther summation is the same in (F1) and (F2). Similarly 14. Interacting heavy fermion Hofstadter spectra for (b) CNP at K contrasted with (a) corresponding zero-field spectra at w0/w1 = 0.7 in the flat band limit Because m 1,−1 + 1 = m 2,−1 = m max + 1, the upper bound onr summation is the same in Eq.(F3) and Eq.(F4). Thus for a given k, out of the available 2q f −modes, only 2m max + 3 couple.

CNP
The interactions inf basis at CNP [56] can be given as Note that there are 2q − (2m max + 3) decoupled f modes for each k. Physically this corresponds to 2 − (2m max + 3)/q states per moiré unit cell. The coupled modes can then be described as where m α=1,...,4 = m α,−1 , m 5 = m max and m 6 = m max + 1, and The following operator can be defined for τ = −1 whereâ is a simple h.o. lowering operator in terms of which matrixΞ mα,m α can be expressed as Here |m is a simple h.o. eigenstate and Σ(m) = Σ m . We now discuss the the exact solutions to the eigenstates of the operator in Eq.(F10) in flat band limit, M = 0. The B field independent J/2 Landau level comes from the anomalous c-mode The rest of the problem can be solved using the following ansätze: 2 |0 , 0, c 2 |1 , 0, c 4 |2 , c 5 |0 , c where m ∈ {2, . . . , m max + 1}. c The ansätzeθ 3 andθ 5 yield the 3×3 and 5×5 Hermitian matrices, whose eigenvectors are c α and c (5) α , respectively: Similarly, the ansatzθ m 6 yields the following 6 × 6 Hermitian matrix for each m, whose eigenvectors are c The magnetic subbands within the narrow bands emanate out of the B → 0 energy eigenvalue of the above decoupled matrices, J 2 , which is m max fold degenerate for matrix in Eq.(F18) and singly degenerate for matrices in Eq.(F16) and Eq.(F17). Including the decoupled c mode, we have m max + 3 magnetic modes emanating out of this m max + 3 fold degenerate B → 0 energy eigenvalue. Now recall that we have 2q − (2m max + 3) decoupled f modes with energy U1 2 . Thus in total we have 2q magnetic modes within the narrow bands, which corresponds to 2 states per moiré unit cell per spin. The spectrum for flat band limit has been shown in Fig.(14b). Note that the B → 0 energies recovered by the decoupled matrices are the corresponding zero field energies of THFM at Γ in mBZ.

ν = ±1
The interactions at ν = ±1 for valley K spin ↑↓ [56] is given as W a∈{1...4} and U 2 are mean field coefficients with W 1 = W 2 and W 3 = W 4 [56]. The interaction for f modes in thef basis can then be given as Yet again there are 2q − (2m max + 3) decoupled f modes for each k. Physically this corresponds to 2 − (2m max + 3)/q states per moiré unit cell for each spin projection. The coupled modes can then be described by where m α=1,...,4 = m α,−1 , m 5 = m max and m 6 = m max + 1, and  can be solved as eigenvectors of the following 3 × 3, 5 × 5 and m max 6 × 6 Hermitian matrices respectively: The magnetic subbands within the narrow bands emanate out of the B → 0 energy eigenvalue of the above decoupled matrices, ν(W 3 − J 2 ), which is m max fold degenerate for matrix in Eq.(F29) and singly degenerate for matrices in Eq.(F27) and Eq.(F28). Including the decoupled c mode, we have m max + 3 magnetic modes emanating out of this m max + 3 fold degenerate B → 0 energy eigenvalue. Now recall that we have 2q − (2m max + 3) decoupled f modes with energy ν( U1 2 + 6U 2 ). Thus in total we have 2q magnetic modes within the narrow bands, which corresponds to 2 states per moiré unit cell for each spin projection. The spectrum for flat band limit has been shown in Fig.(15a). Note that the B → 0 energies recovered by the decoupled matrices are the corresponding zero field energies of THFM at Γ in mBZ.

ν = ±2
The interactions at ν = ±2 for valley K spin ↑↓ [56] is given as The interaction for f modes in thef basis can then be given as Yet again there are 2q −(2m max +3) decoupled f modes for each k. Physically this corresponds to 2−(2m max +3)/q states per moiré unit cell for each spin projection. The coupled modes can then be described by where m α=1,...,4 = m α,−1 , m 5 = m max and m 6 = m max + 1, and where the operatorsĥ −1,s,ν=±2 α,α are given aŝ whereâ is a simple h.o. lowering operator with |m being a simple h.o. eigenstate and Σ(m) = Σ m . The exact eigenstates for the above operator are exactly solvable in flat band limit, M = 0. The field independent ν(W 3 − J 4 ) level is formed by the decoupled anomalous c level given in Eq.(F12). The rest of eigenstates can be solved using the ansätze given in Eq.(F13)-Eq.(F15). The corresponding coefficients c (3) α , c (5) α and c (6,m) α can be solved as eigenvectors of the following 3 × 3, 5 × 5 and m max 6 × 6 Hermitian matrices respectively: h −1,ν=±2 The magnetic subbands within the narrow bands emanate out of the B → 0 energy eigenvalue of the above decoupled matrices, ν(W 3 − J 4 ), which is m max fold degenerate for matrix in Eq.(F40) and singly degenerate for matrices in Eq.(F38) and Eq.(F39). Including the decoupled c mode, we have m max + 3 magnetic modes emanating out of this m max + 3 fold degenerate B → 0 energy eigenvalue. Now recall that we have 2q − (2m max + 3) decoupled f modes with energy ν( 3U1 4 + 6U 2 ). Thus in total we have 2q magnetic modes within the narrow bands, which corresponds to 2 states per moiré unit cell for each spin projection. The spectrum for flat band limit has been shown in Fig.(16b). Note that the B → 0 energies recovered by the decoupled matrices are the corresponding zero field energy of THFM at Γ in mBZ.
Appendix G: ν = ±1, Valley K, Spin ↑ Here we discuss the finite B THFM for valley K spin ↑. The interactions read as [56] V τ =+1,s=↑ The interaction for f modes in thef basis can then be given as where m α=1,...,4 = m α,+1 , m 5 = m max + 1 and m 6 = m max , and with Ξ ↑,ν=±1 mα,m α = m|ĥ +1,↑,ν=±1 where the operatorsĥ +1,↑,ν=±1 α,α are given aŝ whereâ is a simple h.o. lowering operator with |m being a simple h.o. eigenstate and Σ(m) = Σ m . The exact eigenstates for the above operator are exactly solvable in flat band limit, M = 0. The B field independent level ν(W 3 + J/2) is formed by the anomalous c-mode The rest of the problem can be solved using the following ansätze: 3 |1 , 0, c 2 |0 , c 3 |2 , 0, c where m ∈ {2 . . . m max + 1}. c can thus be solved as eigenvectors of the following 3 × 3, 5 × 5 and m max 6 × 6 Hermitian matrices respectively: h +1,ν=±1 h +1,m,ν=±1 The magnetic subbands within the narrow bands emanate out of the B → 0 energy eigenvalue of the above decoupled matrices, ν(W 3 + J 2 ), which is m max fold degenerate for matrix in Eq.(G15) and singly degenerate for matrices in Eq.(G13) and Eq.(G14). Including the decoupled c mode, we have m max + 3 magnetic modes emanating out of this m max + 3 fold degenerate B → 0 energy eigenvalue. Now recall that we have 2q − (2m max + 3) decoupled f modes with energy ν( 3U1 2 + 6U 2 ). Thus in total we have 2q magnetic modes within the narrow bands, which corresponds to 2 states per moiré unit cell. The spectrum for flat band limit has been shown in Fig.(15b). Note that the B → 0 energies recovered by the decoupled matrices are the corresponding zero field energies of THFM at Γ in mBZ.
. Upon canonical substitution, the Hamiltonian in flux reads An example band structure at flux φ = 2π/3 calculated with our technique is depicted in Fig.(18).

Groundstate Charge Density
A key feature of the HF model of TBG is the concentration of charge at the moiré unit cell center in the noninteracting model. This is one indication that the Wannierization must contain orbitals at the 1a position. To check the validity of the HF approximation in small flux, we will compute the charge density of the groundstate in flux φ = 2π q , for this section restricting to p = 1 for convenience. We will compute the quantity where |GS = j k,n,ηj ,sj γ † k,n,ηj ,sj |0 is the Slater determinant state filling all the 2q Hofstadter flat bands and n(r) = l,α,η,s c † l,α,η,s (r)c l,α,η,s (r) is the local charge operator summed over layer l, sublattice α, valley η, and spin s. The vectors G q = 2πZg 1 + 2πZg 2 /q are the reciprocal lattice vectors of the MBZ (note that G q · a 2 can be fractional). The momentum space density operator ρ is defined as ρ q+Gq =ˆd 2 re −i(q+Gq)·r l,α,η,s c † l,α,η,s (r)c l,α,η,s (r).
We first find the matrix elements of this operator in terms of Landau levels. Define the creation operator ψ † k,m,α,l,η,s as the operator that creates the state |k, m, α, l, η, s . Following the calculations from Ref. [83] gives the following FIG. 19. Real space cuts along the L1 lattice vector. All fluxes have many-body charge densities that are peaked at the 1a location, but the 0 flux is the most peaked while the 2π is the most spread. As flux is increased, the wavefunction weight is increasingly spread, but remains strongly localized at 1a. This validates the use of the HF Hamiltonian in flux.
where M, N sum over all the Hofstadter bands and the arrow indicates a discretization of the the magnetic BZ. Acting on the groundstate, we have switching back to the continuum normalization in the second line. As q ranges over the MBZ, the continuum delta function δ(q 1 )δ(q 2 ) appears with a prefactor of (2π) 2 /q, so that the following integral over the MBZ is normalized to 1:ˆ[ 0,2π)×[0,2π/q) d 2 q (2π) 2 /q (2π) 2 q δ(q 1 )δ(q 2 ) = 1.
We now prove that only GS| ρ Gq |GS is nonzero, for G q a reciprocal lattice vector of the zero-flux (non-magnetic) BZ: that is, GS| ρ Gq |GS is nonzero only if G q = 2πZg 1 + 2πZg 2 . This follows from the periodicity of the form factors: which is derived from the form factor definition: 1 q δ(k 1 − k 1 )δ(k 2 − k 2 )e iξ Gq (k+g1/2πq) H Gq mn = k + g 1 /2πq, m|e −2πiGq·r |k + g 1 /2πq, n (J51) = 1 q δ(k 1 − k 1 )δ(k 2 − k 2 )e −iq(π(Gq)1(Gq)2+ ij (Gq)ikj −2π(Gq)2/q) [e i ij 2π(Gq)iZj ] mn (J52) = e 2πiGq·L2 k, m|e −2πiGq·r |k , n , or k + g 1 /2πq, m|e −2πiGq·r |k + g 1 /2πq, n = e 2πiGq·L2 k, m|e −2πiGq·r |k , n j=0 e ijGq·a2 = q if G q · L 2 ∈ Z, and otherwise vanishes: that is, only if G q is a reciprocal lattice vector. This observation allows significant speedup in numerical calculations, as G q vectors that are not reciprocal lattice vectors do not need to be calculated in the expression for the many-body charge density. The integrals in Eq.(J47) will appear throughout the many-body calculations, so we find it convenient to define By proving that n G vanishes if G is not a reciprocal lattice vector, we have shown that the charge distribution at rational flux is periodic over the physical unit cell, not just over the (extended) magnetic unit cell. A figure illustrating the expectation value over a slice in real space along the L 1 direction is depicted in Fig.(19). We see that even in the presence of flux, the charge density remains localized at the 1a position, though the spread does increase slightly with flux.

Many-body Charge Excitations
It was realized in [82] that the double commutator calculation could be extended to TBG in flux. In this appendix, we detail a variation of this calculation using our gauge invariant formalism [83,84], which makes the derivation of the charge excitations very simple. First we recall some basic results from the strong coupling theory of TBG: the double-commutator technique of [27,36] will be used extensively in this section, along with results from [35,83]. Starting with a filled-band groundstate ansatz [12,25,30] at filling ν (ν = 0 is the charge neutrality point) given by Working at p = 1, the equation breaks up the momentum q by restricting q into a vector in the MBZ, then adding a multiple of g 2 /q, and a reciprocal lattice vector G, so that the total momentum q + 2πj q g 2 + G. At charge neutrality, ν = 0, |Ψ ν is a groundstate of H int because H int is positive semi-definite and ρ q+ 2π q jb2+G |Ψ ν = ν 2 (2π) 2 q δ(q 1 )δ(q 2 )δ j,0 qn G |Ψ ν (J63) and thus H int |Ψ ν=0 = 0. This convention justifies the choice of the projected interaction Hamiltonian. More generally, we observe that where (2π) 2 δ 2 (0) is the total number of unit cells. Thus the many-body energy density of |Ψ ν is given by the quantity in parentheses in Eq.(J64). A double-commutator calculation [27,36] gives the charge ±1 excitations in terms of the following 2q × 2q effective Hamiltonians (the φ = 2π/q dependence is suppressed in M (k, q) and n G ): where q sums over the N states in the full Brillouin zone, and the chemical potential µ = (E min,+1 − E min,−1 )/2 is chosen so that the minima of the charge ±1 spectra at µ = 0, denoted E min,+± , are the same. Physically, this choice of µ corresponds to setting zero energy in the middle of the excitation gap. We refer to the first term in R ± (k) as the Fock term. Since V (q) > 0 and M † M is a positive-semi-definite matrix, the Fock term is positive semi-definite. It also does not depend on the filling ν. Note that the phase e iξq(k) in the form factors M (k, q) cancels in the M † M product, and so does not need to be evaluated. To further speed up the numerical calculations, n G can be precomputed since it does not depend on k. Most importantly, the computation of M (k, q) can be carried out with sparse matrix multiplication of e i ij qiZj on U (k). We find that even for very large q, the computations of R ± (k) are feasible taking the Landau level cutoff to be 30q and truncating the G sum after three shells. A plot of the many-body charge excitations calculated with strong coupling is shown in Fig.(20); this is compared to the spectrum obtained via heavy fermions to excellent agreement in the main text.