# Topological heavy fermions in magnetic field

Keshav Singh1,2 , Aaron Chew3 , Jonah Herzog-Arbeitman3 , B. Andrei
Bernevig3,4,5 and Oskar Vafek1,2
1National High Magnetic Field Laboratory, Tallahassee, Florida, 32310
2Department of Physics, Florida State University, Tallahassee, Florida, 32306
3Department of Physics, Princeton University, Princeton, NJ 08544
4Donostia International Physics Center, P. Manuel de Lardizabal 4, 20018
Donostia-San Sebastian, Spain
5IKERBASQUE, Basque Foundation for Science, Bilbao, Spain

###### Abstract

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.

††preprint: APS/THF in B

##  I Introduction

Since the discovery of the remarkable phase diagram of the magic angle twisted
bilayer graphene (MATBG) [1, 2], substantial effort [3, 4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 30, 36, 37, 38, 39, 40, 41, 42, 43, 44, 20, 45, 37,
46, 47, 48, 49, 50, 51] 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}\pm 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 $\sim$13nm in MATBG has revealed a sequence of
broken symmetry Chern insulators yielding a plethora of $\mathbf{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 $\mathbf{B}$ field. Theoretical studies have previously
focused on non-interacting [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 $\mathbf{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
$\mathbf{B}=0$ as shown in Figs.(1),(8),(9) at $\nu=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
$\frac{\Phi}{\Phi_{0}}$, with $\Phi$ being the flux through the unit cell and
$\Phi_{0}$ being the flux quantum $\frac{hc}{e}$, but we focus our analysis on
the $\frac{1}{q}$ flux sequence and low $\mathbf{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
$\mathbf{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.

Figure 1: Interacting heavy fermion Hofstadter spectra (b) contrasted with
zero field strong coupling spectra (a) for flat band ($M$=0) limit at
$w_{0}/w_{1}=0.7$ at CNP. We fix $m_{max}$ = $\lceil\frac{q-3}{2}\rceil$. The
$\mathbf{B}=0$ energies at $\Gamma$, labelled using
$\mathcal{E}_{\pm}=\pm\frac{U_{1}}{4}+\sqrt{\frac{U^{2}_{1}}{16}+\gamma^{2}}$,
are recovered in $\mathbf{B}\rightarrow 0$ limit of our theory. The value of
parameters used are $J=18.27meV$, $U_{1}=51.72meV$, $\gamma=-39.11meV$,
$v^{\prime}_{\ast}=1.624eV.\mathring{A}$, $v_{\ast}=-4.483eV.\mathring{A}$ and
$\lambda=0.3792L_{m}$. We set $k_{1,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]

$\displaystyle\hat{H}_{0}=\sum_{\mid\mathbf{k}\mid<\Lambda_{c}}\sum_{\tau
s}\sum_{aa^{\prime}=1}^{4}H^{c,\tau}_{aa^{\prime}}(\mathbf{k})\tilde{c}^{\dagger}_{\mathbf{k}a\tau
s}\tilde{c}_{\mathbf{k}a^{\prime}\tau s}+$
$\displaystyle\\!\\!\sum_{\mid\mathbf{k}\mid<\Lambda_{c}}\\!\\!\sum_{\tau
s}\sum_{a=1}^{4}\sum_{b=1}^{2}\left(e^{-\frac{1}{2}\mathbf{k}^{2}\lambda^{2}}H^{cf,v}(\mathbf{k})_{ab}\tilde{c}^{\dagger}_{\mathbf{k}a\tau
s}\tilde{f}_{\mathbf{k}b\tau s}\\!\\!+\\!\\!\text{h.c.}\right).$

Here $\Lambda_{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 $\mathbf{B}=0$ Hamiltonian,
$\lambda\approx 0.38L_{m}$ acts as the damping factor for $c$-$f$
hybridization and $L_{m}$ is the moiré period(mBZ). $\tau=+1(-1)$ represents
graphene valley $\mathbf{K}(\mathbf{K}^{\prime})$ and $s$ spin
$\uparrow,\downarrow$. The $c$-$c$ and $c$-$f$ couplings are

$\displaystyle
H^{c,\tau}_{aa^{\prime}}(\mathbf{k})=\left(\begin{array}[]{cc}0_{2\times
2}&v_{\ast}(\tau k_{x}\sigma_{0}+ik_{y}\sigma_{z})\\\ v_{\ast}(\tau
k_{x}\sigma_{0}-ik_{y}\sigma_{z})&M\sigma_{x}\end{array}\right)$ (4)
$\displaystyle
H^{cf,\tau}_{ab}(\mathbf{k})=\left(\begin{array}[]{c}\gamma\sigma_{0}+v^{\prime}_{\ast}(\tau
k_{x}\sigma_{x}+k_{y}\sigma_{y})\\\ 0_{2\times 2}\end{array}\right),$ (7)

where the parameters $v_{\ast}$, $v^{\prime}_{\ast}$, $\gamma$ 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
$\mathbf{B}$ solution, we first focus on the charge neutrality point (CNP)
wherein for valley polarised (VP) states, the mean-field interaction reads

$\displaystyle V=-\\!\\!\sum_{\mid\mathbf{k}\mid<\Lambda_{c}\tau
s}\\!\\!\tau\left(\sum_{a=3}^{4}\frac{J}{2}\tilde{c}^{\dagger}_{\mathbf{k}a\tau
s}\tilde{c}_{\mathbf{k}a\tau
s}+\sum_{b=1}^{2}\frac{U_{1}}{2}\tilde{f}^{\dagger}_{\mathbf{k}b\tau
s}\tilde{f}_{\mathbf{k}b\tau s}\right).$ (8)

The narrow band filling factors $\nu=\pm 1$ and $\pm 2$ are discussed in the
section III.3.

Figure 2: The evolution of the interacting flat band Hofstadter spectra is
shown as the $c-f$ coupling is increased from 0 to 1 via the scaling factor
$\beta$ defined via replacement $H^{\tau}_{cf}\rightarrow\beta H^{\tau}_{cf}$.
Valley $\mathbf{K}^{\prime}$ is shown and $w_{0}/w_{1}=0.7$ at CNP. For
illustration we have fixed $m_{max}$ = 5. (a) At $\beta=0$, the $c$-fermion
Landau levels and $2q$ $f$-fermions at $U_{1}/2$ are decoupled. (b) As the
$c-f$ coupling is tuned to $35$% of its full value, $2m_{max}+3$ of the $f$
modes start coupling with the $c$ Landau levels and the remaining
$2q-(2m_{max}+3)$ remain at $U_{1}/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 $\mathbf{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
$\frac{J}{2}$ to $\frac{U_{1}}{2}$, is well gapped from the remote magnetic
subbands emerging out of $\mathbf{B}\rightarrow 0$ energies $\mathcal{E}_{+}$
and $-\mathcal{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$.

The real space basis for $c$’s is composed of four $\mathbf{k}\cdot\mathbf{p}$
Bloch states $\tilde{\Psi}_{\mathbf{k}a\tau}(\mathbf{r})$ at $\Gamma$ point in
mBZ. The localized $f$ electron Wannier functions are denoted by
$W_{\mathbf{R},b\tau}(\mathbf{r})$; they behave as $p_{x}\pm ip_{y}$ orbitals
sitting on the AA stacking sites spanned by moiré triangular lattice vector
$\mathbf{R}$, and
$W_{\mathbf{\mathbf{R}},b\tau}(\mathbf{r})=W_{\mathbf{0},b\tau}(\mathbf{r}-\mathbf{R})$.
The $c$-$c$ and $c$-$f$ couplings in Eqs.(4)-(7) are related to the real space
basis as

$H^{c,\tau}_{aa^{\prime}}(\mathbf{k})=\int
d^{2}\mathbf{r}\tilde{\Psi}^{\ast}_{\mathbf{k}a\tau}(\mathbf{r})H^{\tau}_{BM}(\mathbf{p})\tilde{\Psi}^{\ast}_{\mathbf{k}a^{\prime}\tau}(\mathbf{r}),$
(9)
$\frac{e^{-\frac{1}{2}\mathbf{k}^{2}\lambda^{2}}{2}}{\sqrt{N}}H^{cf,\tau}_{ab}(\mathbf{k})=\int
d^{2}\mathbf{r}e^{i\mathbf{k}\cdot\mathbf{r}}\tilde{\Psi}^{\ast}_{\Gamma
a\tau}(\mathbf{r})H^{\tau}_{BM}(\mathbf{p})W_{\mathbf{0},b\tau}(\mathbf{r}),\\\
$ (10)

where $N$ is number of moiré unit cells and the well known BM Hamiltonian [85]
(also see SM-A), $H^{\tau}_{BM}(\mathbf{p})$, is invariant under translation
by integer multiples of moiré lattice vectors
$\mathbf{L}_{1}=L_{m}(\frac{\sqrt{3}}{2},\frac{1}{2})$ and
$\mathbf{L}_{2}=L_{m}(0,1)$. The corresponding primitive moiré reciprocal
lattice vectors are $\mathbf{g}_{1}=\frac{4\pi}{\sqrt{3}L_{m}}(1,0)$ and
$\mathbf{g}_{2}=\frac{2\pi}{\sqrt{3}L_{m}}(-1,\sqrt{3})$.

In the presence of out-of-plane magnetic field, employed via Landau gauge
vector potential $\mathbf{A}=(0,Bx,0)$,
$H^{\tau}_{BM}(\mathbf{p}-\frac{e}{c}\mathbf{A})$ is still invariant under the
translation by $\mathbf{L}_{2}$, but a translation by $\mathbf{L}_{1}$ needs
to be accompanied by a gauge transformation

$\displaystyle
f(\mathbf{r})\rightarrow\hat{t}_{\mathbf{L}_{1}}f(\mathbf{r})=\exp\left(i\frac{eB}{\hbar
c}L_{1x}y\right)f(\mathbf{r}-\mathbf{L}_{1}).$ (11)

Thus, if $f(\mathbf{r})$ is an eigenstate of
$H^{\tau}_{BM}(\mathbf{p}-\frac{e}{c}\mathbf{A})$, then so is
$\hat{t}_{\mathbf{L}_{1}}f(\mathbf{r})$ with the same eigenvalue. Translations
by $\mathbf{L}_{2}$ are generated by
$\hat{t}_{\mathbf{L}_{2}}f(\mathbf{r})=f(\mathbf{r}-\mathbf{L}_{2})$. Then
$\hat{t}_{\pm\mathbf{L}_{2}}\hat{t}_{\mathbf{L}_{1}}=e^{\mp 2\pi
i\frac{\phi}{\phi_{0}}}\hat{t}_{\mathbf{L}_{1}}\hat{t}_{\mathbf{L}_{2}}$,
where $\phi_{0}=\frac{hc}{e}$ and the flux through the unit cell is
$\phi=BL_{1x}L_{m}$. For $\phi/\phi_{0}=p/q$, with $p$ and $q$ relatively
prime integers, the magnetic translation operators obey
$\left[\hat{t}^{q}_{\mathbf{L}_{2}},\hat{t}_{\mathbf{L}_{1}}\right]=0$.

We construct our $\mathbf{B}\neq 0$ states by using irreps of the magnetic
translation group (MTG) and the hybrid Wannier method [82, 86]. Because they
are topologically trivial and well localized, a complete basis for heavy
fermions can be generated from the two $\mathbf{B}=0$ $f$-fermion Wannier
states as

$\eta_{b\tau
k_{1}k_{2}}(\mathbf{r})=\frac{1}{\sqrt{\mathcal{N}}}\sum_{s,n\in\mathbb{Z}}e^{2\pi
i(sk_{1}+nk_{2})}\hat{t}^{s}_{\mathbf{L}_{1}}\hat{t}^{n}_{\mathbf{L}_{2}}W_{\mathbf{0},b\tau}(\mathbf{r}).$
(12)

Here $k_{1,2}\in[0,1)$ (see SM-C) and the normalization factor
$\mathcal{N}=s_{tot}n_{tot}$, where $s_{tot}$ and $n_{tot}$ denote the total
count of $s$ and $n$ (see SM-C). The $c$-fermion basis can in turn be expanded
as the product of Landau levels (LLs) and the $\mathbf{B}=0$ Bloch states at
the $\Gamma$-point, a result obtained when the $\mathbf{k}\cdot\mathbf{p}$
method is extended to $\mathbf{B}\neq 0$ [87] . In order to use the same
quantum numbers as for $f$s, we also project the $c$’s LL wavefunctions
$\Phi_{m}$ onto the (orthonormal) irreps of MTG

$\chi_{k_{1}k_{2}m}(\mathbf{r})=\frac{1}{\sqrt{\ell
L_{m}}}\frac{1}{\sqrt{\mathcal{N}}}\sum_{s\in\mathbb{Z}}e^{2\pi
isk_{1}}\hat{t}^{s}_{\mathbf{L}_{1}}\Phi_{m}(\mathbf{r},k_{2}\mathbf{g}_{2}).$
(13)

Here $k_{1}\in[0,1)$, but unlike in Eq.(12), $k_{2}\in[0,\frac{p}{q})$ (see
SM-C) and

$\displaystyle\Phi_{j}(\mathbf{r},k_{2}\mathbf{g}_{2})=e^{2\pi
ik_{2}\frac{y}{L_{m}}}\varphi_{j}\left(x-k_{2}\frac{2\pi\ell^{2}}{L_{m}}\right),$
(14)

where the harmonic oscillator (h.o.) wavefunctions
$\varphi_{m}(x)=e^{-x^{2}/2\ell^{2}}H_{m}(x/\ell)/\pi^{\frac{1}{4}}\sqrt{2^{m}m!}$
with Hermite polynomials $H_{m}$, and $\ell^{2}=\hbar c/(eB)$. The $k_{2}$
induced offset in the h.o. wavefunctions is thus $qk_{2}L_{1x}/p$. Note that
$\Phi_{m}$ is an eigenstate of $\hat{t}_{\pm\mathbf{L}_{2}}$, and

$\hat{t}_{\pm\mathbf{L}_{2}}\chi_{k_{1}k_{2}m}(\mathbf{r})=e^{\mp 2\pi
ik_{2}}\chi_{[k_{1}\mp\frac{p}{q}]_{1}k_{2}m}(\mathbf{r}),$ (15)

where $[b]_{a}$ represents $b$ modulo $a$ with $a>0$. Since the
$q\mathbf{L}_{2}$ translations break up the $k_{2}$ domain into units of width
$\frac{1}{q}$, from here on we use the label $k=(k_{1},k_{2})$ with
$k_{1}\in[0,1)$ and we fix $k_{2}\in[0,\frac{1}{q})$. The original $k_{2}$
domains are then accessed using labels $r^{\prime}\in\\{0,\ldots q-1\\}$ and
$r\in\\{0,\ldots,p-1\\}$, denoting the magnetic strip
$[\frac{r^{\prime}}{q},\frac{r^{\prime}+1}{q})$ and
$[\frac{r}{q},\frac{r+1}{q})$ along $\mathbf{g}_{2}$ for $\eta$ and $\chi$
respectively. The $\mathbf{B}\neq 0$ Hamitonian cannot mix states with
different $k$. The relabelled finite field basis is $\eta_{b\tau
kr^{\prime}}(\mathbf{r})$ and $\chi_{krm}(\mathbf{r})$ and their annihilation
operators are $f_{b\tau kr^{\prime}}$ and $c_{a\tau krm}$.

Having assembled the low energy basis at $\mathbf{B}\neq 0$, we now expand the
fields for a given spin as

$\displaystyle\psi_{\tau}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\sum_{k\in[0,1)\otimes[0,\frac{1}{q})}\left(\sum_{b=1}^{2}\sum_{r^{\prime}=0}^{q-1}\eta_{b\tau
kr^{\prime}}(\mathbf{r})f_{b\tau kr^{\prime}}+\right.$ (16)
$\displaystyle\left.\sum_{a=1}^{4}\sum_{m=0}^{m_{a,\tau}}\sum_{r=0}^{p-1}\Psi_{a\tau}(\mathbf{r})\chi_{krm}(\mathbf{r})c_{a\tau
krm}\right),$

where $\Psi_{a\tau}(\mathbf{r})=\sqrt{NA_{uc}}\tilde{\Psi}_{\Gamma
a\tau}(\mathbf{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 $\tau$, denoted by $m_{a,\tau}$, 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, $\Psi_{a\tau}(\mathbf{r})\chi_{krm}(\mathbf{r})$ and
$\eta_{b\tau kr^{\prime}}(\mathbf{r})$ are eigenstates of
$\hat{t}_{\mathbf{L}_{1}}$ and $\hat{t}^{q}_{\mathbf{L}_{2}}$ with eigenvalues
$e^{-2\pi ik_{1}}$ and $e^{-2\pi iqk_{2}}$ respectively. This guarantees
states with different $k$ to be orthogonal. Overlaps of states with the same
$k$ but different $r,r^{\prime}$ 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 $\ell^{2}\mathbf{g}^{2}$
as long as the LL index $m$ is held below an upper cutoff $m_{max}\lesssim
q/2$, where $\mathbf{g}$ is the reciprocal moiré lattice vector.

The single-particle THFM at $\mathbf{B}\neq 0$ is then

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\sum_{\tau}\int
d^{2}\mathbf{r}\psi^{\dagger}_{\tau}(\mathbf{r})H^{\tau}_{BM}\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)\psi_{\tau}(\mathbf{r})$
(17) $\displaystyle=$
$\displaystyle\sum_{\tau}H^{\tau}_{cc}+\sum_{\tau}\left(H^{\tau}_{cf}+\text{h.c.}\right),$
(18)

where the $cc$ couplings are

$\displaystyle H^{\tau}_{cc}$ $\displaystyle=$
$\displaystyle\sum_{k\in[0,1)\otimes[0,\frac{1}{q})}\sum_{a,a^{\prime}=1}^{4}\sum_{m=0}^{m_{a,\tau}}\sum_{m^{\prime}=0}^{m_{a^{\prime},\tau}}\sum_{r,\tilde{r}=0}^{p-1}$
(19)
$\displaystyle\tilde{h}^{\tau}_{[amr],[a^{\prime}m^{\prime}\tilde{r}]}(k)c^{\dagger}_{a\tau
krm}c_{a^{\prime}\tau k\tilde{r}m^{\prime}},$

with matrix elements that are found to be the same as those obtained by the
direct minimal substitution in (4) and expanding in LL basis, as is expected
from $\mathbf{k}\cdot\mathbf{p}$[87]:

$\displaystyle\tilde{h}^{\tau}_{[amr],[a^{\prime}m^{\prime}\tilde{r}]}(k)=\delta_{r\tilde{r}}\left(\begin{array}[]{cc}0_{2\times
2}&h^{\tau,c}_{mm^{\prime}}\\\
\sigma_{x}h^{\tau,c}_{mm^{\prime}}\sigma_{x}&M\delta_{mm^{\prime}}\sigma_{x}\end{array}\right)_{aa^{\prime}},$
(22)

where the Pauli matrix $\sigma_{x}$ acts on the $c$ orbitals and

$\displaystyle
h^{+1,c}_{mm^{\prime}}=i\frac{\sqrt{2}v_{\ast}}{\ell}\left(\begin{array}[]{cc}-\sqrt{m^{\prime}}\delta_{m+1,m^{\prime}}&0\\\
0&\sqrt{m}\delta_{m,m^{\prime}+1}\end{array}\right)$ (25)

with $h^{-1,c}_{mm^{\prime}}=-\sigma_{x}h^{+1,c}_{mm^{\prime}}\sigma_{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 $\mathbf{B}$, we plot their and $f$’s spectrum
with the term in Eq.(8) included but turn off the $c$-$f$ coupling by setting
the prefactor of $H^{\tau}_{cf}$ to zero in Fig.(2a). For each spin
projection, there are two $f$s per moiré unit cell at the energy
$\pm\frac{U_{1}}{2}$ at $\tau=\mp 1$ intersected by the LLs with strong
$\mathbf{B}$ dependence due to $c$’s fast $v_{*}$. The total number of $f$
states is $\mathbf{B}$-independent while each LL contains
$N\frac{\phi}{\phi_{0}}$ states. We expect that $H^{\tau}_{cf}$ hybridizes
these modes within an energy $\sim\gamma$, but in order to understand how such
hybridization leaves behind a band of states whose total number is independent
of $\mathbf{B}$ –because its total Chern number vanishes[88]– and which LLs
may decouple from $f$s, we need to carefully analyze $H^{\tau}_{cf}$. As shown
below, it is not simply minimally substituted (10). Instead

$\displaystyle H^{\tau}_{cf}$ $\displaystyle=$
$\displaystyle\sum_{k\in[0,1)\otimes[0,\frac{1}{q})}\sum_{a=1}^{4}\sum_{b=1}^{2}\sum_{m=0}^{m_{a,\tau}}\sum_{r=0}^{p-1}\sum_{r^{\prime}=0}^{q-1}h^{\tau}_{[amr],[br^{\prime}]}(k)$
(26) $\displaystyle c^{\dagger}_{a\tau krm}f_{b\tau kr^{\prime}},$

where, using the magnetic translation symmetry of
$H^{\tau}_{BM}\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)$ and the free sum
over the $s$, we find

$\displaystyle
h^{\tau}_{[amr],[br^{\prime}]}(k)=\frac{s_{tot}}{\sqrt{\mathcal{N}}}\sum_{n\in\mathbb{Z}}e^{2\pi
i(k_{2}+\frac{r^{\prime}}{q})n}\times$ $\displaystyle\int
d^{2}\mathbf{r}(\hat{t}^{n}_{-\mathbf{L}_{2}}\chi_{krm}(\mathbf{r}))^{\ast}\Psi^{*}_{a\tau}(\mathbf{r})H^{\tau}_{BM}\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)W_{\mathbf{0},b\tau}(\mathbf{r}).$

Because $\mathbf{A}$ acts on a well-localized function centered at
$\mathbf{r}=0$, its contribution can be neglected at low $\mathbf{B}$ (we have
confirmed this numerically in SM-D.3). The dominant term can be obtained using
the inverse Fourier transform of Eq.(10) which gives
$\Psi^{\ast}_{a\tau}(\mathbf{r})H^{\tau}_{BM}\left(\mathbf{p}\right)W_{\mathbf{0},b\tau}(\mathbf{r})$.
Substituting into (II) thus reduces $h^{\tau}_{[amr],[br^{\prime}]}(k)$ to
integrals over a gaussian and shifted 1D h.o. wavefunctions which can be
evaluated using results in Ref.[89] (see also SM-D.2).

##  III Analysis of the $\frac{1}{q}$ flux sequence

It is particularly revealing to analyze the case $p=1$, i.e. the
$\phi/\phi_{0}=1/q$ sequence. Then $r=0$, because, as discussed above, $r$
ranges from $0$ to $p-1$. Making use of (15) in (II) and substituting the
explicit expression for $\chi$ from (13), allows us to perform the sum over
$n$. We find

$\displaystyle
h^{\tau}_{[am0],[br^{\prime}]}(k)=\sqrt{\frac{L_{1x}}{\ell}}\sum_{j\in\mathbb{Z}}e^{-2\pi
i(r^{\prime}+jq)k_{1}}\times$ $\displaystyle\int
d^{2}\mathbf{r}\left(\hat{t}^{r^{\prime}+jq}_{\mathbf{L}_{1}}\Phi_{m}(\mathbf{r},k_{2}\mathbf{g}_{2})\right)^{\ast}H_{ab}^{cf,\tau}\left(\frac{1}{i}\frac{\partial}{\partial\mathbf{r}}\right)\frac{e^{-\frac{\mathbf{r}^{2}}{2\lambda^{2}}}}{2\pi\lambda^{2}}.$

The above expression can be visualized as an overlap between a 2D localized
heavy state with size $\lambda$ 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^{\prime},j$ is this integral significant, note that the h.o. wavefunction
is localized in the $x$ direction about
$\left(r^{\prime}+jq+k_{2}q\right)L_{1x}$, and its width is $\sim
2\sqrt{2m+1}\sqrt{q}$ unit cells. In addition, the combination
$\frac{r^{\prime}}{q}+j+k_{2}$ controls the period of oscillation in the
$y$-direction set by $1/(\frac{r^{\prime}}{q}+j+k_{2})$ times the unit cell
size. The integer $r^{\prime}+jq$ thus determines the unit cell to which the
h.o. is shifted, and, because $k_{2}2\pi\ell^{2}/L_{m}=k_{2}qL_{1x}$, the
value of $k_{2}q\in[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.

Figure 3: Schematic representation of the $cf$ coupling for
$\phi/\phi_{0}=1/q$ in Eq.(III). Each tick represents (the 1D projection of) a
moire unit cell, illustrating the overlap between a 2D localized heavy state
with size $\lambda$ 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^{\prime}$ determines the momentum
absorbed by the LL as well as the unit cell to which the h.o. is shifted and
$k_{2}q$ fine tunes the shift within the unit cell. The index $j$ then
determines $q$-unit-cell periodic revival of the h.o. The black parabolas
mimic a quadratic potential to accompany the h.o. wavefunctions.

For example, if $r^{\prime}=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
$\gamma$ will then be significant. Since the spatial extent of the h.o. state
in the $x$-direction is $\sim 2\sqrt{2m+1}\sqrt{q}$ unit cells –which at low
$\mathbf{B}$ is much longer than the localized heavy state even when we
consider that it oscillates and has $m$-nodes– the 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}\lesssim q/2$. If we keep $j=0$ but increase $r^{\prime}$ 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 $\gamma$ 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^{\prime}$ past $\sim\sqrt{2m+1}\sqrt{q}$ results in an exponential
suppression of the hybridization integral due to the large off-set. Therefore,
for the values of $r^{\prime}$ 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.(III) and
replace

$r^{\prime}+jq\rightarrow\text{sgn}_{+}\left(\frac{q}{2}-r^{\prime}\right)\text{min}[r^{\prime},q-r^{\prime}]\equiv
r^{\prime}_{q},$ (29)

where $\text{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$. It also increases the number of nodes in the 1D h.o. wavefunctions with
faster oscillations in the $x$-direction. This results in stronger averaging
when convolving with the localized state and the suppression of the overall
amplitude of the hybridization. Similar averaging occurs upon increasing
$\mathbf{B}$. The suppression of the overall amplitude is reflected in the
decrease of the singular values of the $(m,r^{\prime})$-matrix
$h^{\tau}_{[1m0],[1r^{\prime}]}(k)$ shown in Fig.(5) after rescaling by
$\gamma$, which enters our analysis below. The derivatives in Eq.(III),
appearing in matrix elements $h^{\tau}_{[1(2)m0],[2(1)r^{\prime}]}$, act on
the localized heavy function to change its spatial symmetry from $s$ to
$p_{x,y}$-like. Moreover an integration by parts and expressing the
derivatives via h.o. raising and lowering operators allows us to relate these
cases to the analysis without derivatives(see SM-D.2.2,D.2.3).

The qualitative trends discussed above are captured by our closed form
expression for $h^{\tau}_{[am0],[br^{\prime}]}(k)$. This non-derivative
hybridization for the $\frac{1}{q}$ sequence can be written as

$h^{\tau}_{[am0],[ar^{\prime}]}(k)=\gamma
I^{0}_{m,r^{\prime}}(k),\;\;\text{for}\;a=1,2.$ (30)

Here

$\displaystyle I^{0}_{m,r^{\prime}}(k)=\sqrt{\frac{L_{1x}}{\ell}}e^{i\pi
r^{\prime}_{q}\left(k_{2}-2k_{1}\right)}e^{i\pi
r^{\prime}_{q}(r^{\prime}_{q}-1)\frac{1}{2q}}$ $\displaystyle
e^{-2\pi^{2}\frac{\lambda^{2}}{L^{2}_{m}}\left(k_{2}+\frac{r^{\prime}_{q}}{q}\right)^{2}}\mathcal{F}_{m}\left(\lambda,(r^{\prime}_{q}+qk_{2})L_{1x}\right),$
(31)

and

$\displaystyle\mathcal{F}_{m}(\lambda,x_{0})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{\frac{1}{4}}\sqrt{2^{m}m!}}\sqrt{\frac{\ell^{2}}{\ell^{2}+\lambda^{2}}}e^{-\frac{x_{0}^{2}}{2(\ell^{2}+\lambda^{2})}}$
(32)
$\displaystyle\times\mathcal{H}_{m}\left(\frac{-2x_{0}\ell}{\ell^{2}+\lambda^{2}},-1+\frac{2\lambda^{2}}{\ell^{2}+\lambda^{2}}\right).$

The two variable Hermite polynomials[89] are given by

$\displaystyle\mathcal{H}_{m}(x,y)=m!\sum_{k=0}^{\lfloor\frac{m}{2}\rfloor}\frac{x^{m-2k}y^{k}}{(m-2k)!k!},$
(33)

where $\lfloor m\rfloor$ denotes the floor function at $m$. Their relation to
the Hermite polynomials used above is $H_{m}(x)=\mathcal{H}_{m}(2x,-1)$. In
the limit $\lambda\rightarrow 0$, the 2D heavy localized state becomes the
Dirac $\delta$-function, and as expected based on the above discussion, we
indeed recover $\mathcal{F}_{m}(\lambda\rightarrow
0,x_{0})=\varphi_{m}(-x_{0})$. We found that keeping the full form of
$\mathcal{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
$\mathcal{F}_{m}$(see also Fig.(12) in SM-E). The exponential suppression
factor multiplying $\mathcal{F}_{m}$ in the Eq.(III) comes from the
$y$-integration. Although the suppression depends on the shifted position of
the h.o. states, its dependence on $r^{\prime}$ is weaker than in
$\mathcal{F}_{m}$ which comes from the $x$-integration, as expected based on
the discussion above.

The off-diagonal hybridization can be expressed as

$\displaystyle h^{\tau}_{[1m0],[2r^{\prime}]}(k)$ $\displaystyle=$
$\displaystyle v^{\prime}_{\ast}\left(\tau
I^{x}_{m,r^{\prime}}(k)-iI^{y}_{m,r^{\prime}}(k)\right)$ (34) $\displaystyle
h^{\tau}_{[2m0],[1r^{\prime}]}(k)$ $\displaystyle=$
$\displaystyle\left(h^{\tau}_{[1m0],[2r^{\prime}]}(k)\right)^{\ast}.$ (35)

Using the exact recursion relations between $\mathcal{H}_{m}$ and its
derivatives, $I^{x}$ and $I^{y}$ can be expressed in terms of $I^{0}$ as

$I^{x}_{m,r^{\prime}}(k)=\frac{i}{\sqrt{2}\ell}\left(\sqrt{m}I^{0}_{m-1,r^{\prime}}(k)-\sqrt{m+1}I^{0}_{m+1,r^{\prime}}(k)\right)$
(36) $\displaystyle I^{y}_{m,r^{\prime}}(k)$ $\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2}\ell}\left(\left(1+\frac{\lambda^{2}}{\ell^{2}}\right)\sqrt{m+1}I^{0}_{m+1,r^{\prime}}(k)\right.$
(37) $\displaystyle+$
$\displaystyle\left.\left(1-\frac{\lambda^{2}}{\ell^{2}}\right)\sqrt{m}I^{0}_{m-1,r^{\prime}}(k)\right).$

We add interactions at CNP via the mean field coefficients $J$ and $U_{1}$.
This is accomplished by promoting the operators $\tilde{c}_{\mathbf{k}a\tau}$
and $\tilde{f}_{\mathbf{k}b\tau}$ in Eq.(8) to $c_{a\tau krm}$ and $f_{b\tau
kr^{\prime}}$ respectively, with a sum over all finite $\mathbf{B}$ quantum
numbers. The numerically determined strong coupling Hofstadter spectra using
the above results for the matrix elements for $\tau=-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.

Figure 4: Comparison of THFM Hofstadter spectrum (including dispersion of the
flat bands) and strong coupling projected BM Hofstadter spectrum at
$w_{0}/w_{1}=0.7$. (a) THFM spectrum at with $M=3.248meV$ for
$m_{max}=\lceil\frac{q-3}{2}\rceil$. (b) Exact spectrum in the flat band limit
using the gauge-invariant basis of magnetic translation group irreps. Details
of the method may be found in SM-J

.

Neglecting $\frac{\lambda^{2}}{\ell^{2}}$ in the prefactors of $I^{0}$’s in
Eq.(37) is an approximation of same order as neglecting $\mathbf{A}$ in
$H^{\tau}_{BM}$ in Eq.(II). We confirm it numerically to be an excellent
approximation in SM-D.3; note that we do not set
$\frac{\lambda^{2}}{\ell^{2}}$ to zero in $I^{0}$ for the reasons mentioned
above. Thus the finite $\mathbf{B}$ $c$-$f$ coupling can be accurately re-
expressed as

$\displaystyle
h^{\tau}_{[am0],[br^{\prime}]}(k)\approx\left(\begin{array}[]{c}\gamma
I^{0}_{m,r^{\prime}}(k)\sigma_{0}+h^{\tau,cf}_{m,r^{\prime}}(k)\\\ 0_{2\times
2}\end{array}\right)_{ab},$ (40)

where $h^{+1,cf}_{m,r^{\prime}}(k)=$

$\displaystyle
i\frac{\sqrt{2}v^{\prime}_{*}}{\ell}\left(\begin{array}[]{cc}0&\sqrt{m}I^{0}_{m-1,r^{\prime}}(k)\\\
-\sqrt{m+1}I^{0}_{m+1,r^{\prime}}(k)&0\end{array}\right),$ (43)

and
$h^{-1,cf}_{m,r^{\prime}}(k)=-\sigma_{x}h^{+1,cf}_{m,r^{\prime}}(k)\sigma_{x}$,
where $\sigma_{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.

###  III.1 Change of basis for the $f$-modes

We will find it useful to first perform a singular value decomposition

$I^{0}_{m,r^{\prime}}(k)=U_{mm^{\prime}}\Sigma_{m^{\prime}\tilde{r}}V_{\tilde{r}r^{\prime}},$
(44)

where $U,V$ are are $m_{a,\tau}\times m_{a,\tau}$ and $q\times q$ unitary
matrices respectively, and the $m_{a,\tau}\times q$ rectangular matrix
$\Sigma$ 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,

$\displaystyle h^{\tau}_{[1m0],[1r^{\prime}]}c^{\dagger}_{1\tau k0m}f_{1\tau
kr^{\prime}}=\gamma\Sigma_{m^{\prime},\tilde{r}}\bar{c}^{\dagger}_{1\tau
k0m^{\prime}}\bar{f}_{1\tau k\tilde{r}},$ (45)

where $\bar{c}^{\dagger}_{1\tau k0m^{\prime}}=c^{\dagger}_{1\tau
k0m}U_{mm^{\prime}}$ and $\bar{f}_{1\tau
k\tilde{r}}=V_{\tilde{r}r^{\prime}}f_{1\tau kr^{\prime}}$ are the new modes.
In the above, the $\bar{f}_{1\tau k\tilde{r}}$ modes with
$\tilde{r}>m_{1,\tau}$ 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 $\gamma\Sigma_{m^{\prime},\bar{r}}=\gamma\Sigma_{m^{\prime}}$, where
$\Sigma_{m}$ denotes the $m^{th}$ singular value.

The columns of matrix $U$ are the eigenvectors of the following matrix

$\displaystyle\Lambda^{0}_{mm^{\prime}}(k)=\sum_{r^{\prime}=0}^{q-1}I^{0}_{mr^{\prime}}(k)I^{0^{\ast}}_{m^{\prime}r^{\prime}}(k).$
(46)

Since $r^{\prime}$ is being summed over, we are allowed to shift its range to
$-\lfloor\frac{q}{2}\rfloor,-\lfloor\frac{q}{2}\rfloor+1,\ldots,q-1-\lfloor\frac{q}{2}\rfloor$
and thus set $j=0$ in Eq.(29). Hence we can set $r^{\prime}_{q}=r^{\prime}$ in
Eq.(III) for $I^{0}_{m,r^{\prime}}$. Since $\mathcal{F}_{m}$ decays
exponentially at $\pm\lfloor\frac{q}{2}\rfloor$, in the low field limit, we
can replace the bounds of this sum by $\pm\infty$. Using Eq.(III) we have

$\displaystyle\Lambda^{0}_{mm^{\prime}}(k)\approx\frac{L_{1x}}{\ell}\sum_{r^{\prime}=-\infty}^{\infty}e^{-4\pi^{2}\frac{\lambda^{2}}{L^{2}_{2}}\left(k_{2}+\frac{r^{\prime}}{q}\right)^{2}}\times$
$\displaystyle\mathcal{F}_{m}\left(\lambda,\left(r^{\prime}+k_{2}q\right)L_{1x}\right)\mathcal{F}_{m^{\prime}}\left(\lambda,\left(r^{\prime}+k_{2}q\right)L_{1x}\right).$
(47)

Using the Dirac comb identity,
$\sum_{r^{\prime}\in\mathbb{Z}}\delta(\rho-r^{\prime})=\sum_{t\in\mathbb{Z}}e^{2\pi
it\rho}$, we can convert this summation to an integral as

$\displaystyle\Lambda^{0}_{mm^{\prime}}(k)$ $\displaystyle\approx$
$\displaystyle\frac{L_{1x}}{\ell}\sum_{t\in\mathbb{Z}}\int_{-\infty}^{\infty}d\rho
e^{i2\pi
t\rho}e^{-4\pi^{2}\frac{\lambda^{2}}{L^{2}_{2}}\left(k_{2}+\frac{\rho}{q}\right)^{2}}$
(48)
$\displaystyle\mathcal{F}_{m}\left(\lambda,\left(\rho+k_{2}q\right)L_{1x}\right)\mathcal{F}_{m^{\prime}}\left(\lambda,\left(\rho+k_{2}q\right)L_{1x}\right)$
$\displaystyle=$ $\displaystyle\frac{1}{\ell}\sum_{t\in\mathbb{Z}}e^{-i2\pi
tk_{2}q}\int_{-\infty}^{\infty}d\rho e^{i2\pi
t\rho/L_{1x}}e^{-\frac{\lambda^{2}}{\ell^{2}}\frac{\rho^{2}}{\ell^{2}}}$
$\displaystyle\mathcal{F}_{m}\left(\lambda,\rho\right)\mathcal{F}_{m^{\prime}}\left(\lambda,\rho\right).$

Because $2\pi t/L_{1x}$ is a large wavevector compared to the length scale at
which the h.o. wavefunctions vary, $\sim\ell$, i.e. $\frac{L_{1x}}{2\pi
t}\ll\ell$ at low field, the overlaps are exponentially suppressed in
$(\frac{2\pi t\ell}{L_{1x}})^{2}$. We can thus set $t=0$ in the above sum.
Moreover since the exponential factor contains $\frac{\lambda^{2}}{\ell^{2}}$,
the off-diagonal terms can be neglected and thus $\Lambda^{0}_{mm^{\prime}}$
is diagonal up to a good approximation. This implies that the matrix $U$ is
close to identity and the low $\mathbf{B}$ singular values can be given as

$\Sigma_{m}=\sqrt{\Lambda^{0}_{mm}}.$ (49)

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 $\Sigma_{m}$ can be computed from
Eq.(48) using two index Hermite polynomials[89](also see SM-D.2)

$\displaystyle\Sigma_{m}=\left(\frac{1}{\sqrt{\xi(\kappa)}}\frac{1}{2^{m}m!}\mathcal{H}_{mm}\left(0,\frac{\kappa^{6}}{\xi(\kappa)};0,\frac{\kappa^{6}}{\xi(\kappa)}\mid\frac{2}{\xi(\kappa)}\right)\right)^{\frac{1}{2}},$
(50)

where
$\kappa^{2}=\frac{\lambda^{2}}{\ell^{2}}=(\frac{\phi}{\phi_{0}})2\pi\lambda^{2}/A_{u.c}$
and

$\displaystyle\xi(\kappa)=(1+\kappa^{2}+\kappa^{4})(1+\kappa^{2}),$ (51)
$\displaystyle\mathcal{H}_{mn}\left(x,y;w,z|\beta\right)=$
$\displaystyle\sum^{min(m,n)}_{k=0}\frac{m!n!\beta^{k}}{(n-k)!(m-k)!k!}\mathcal{H}_{m-k}(x,y)\mathcal{H}_{n-k}(w,z).$

Figure 5: First 25 analytical singular values obtained via Eq.(50) compared to
the numerically obtained singular values of $I^{0}$ defined in the Eq.(III) at
$\lambda=0.3792L_{m}$

In the $\mathbf{B}\rightarrow 0$ limit, upto first order in flux, the singular
values can be approximated as

$\Sigma_{m}\approx 1-\left(m+\frac{1}{2}\right)\kappa^{2}.$ (53)

Intuitively this shows that $\Sigma_{m}$ dependence on $m$ grows linearly with
$m$ even at negligibly small fluxes. The comparison of analytical singular
values given by Eq.(50) with the ones obtained numerically from Eq.(III) is
presented in Fig.(5), limiting to the values given in Eq.(53) as shown in
SM-E.

Figure 6: The non-interacting heavy fermion Hofstadter spectra for flat band
THFM at $w_{0}/w_{1}=0.7$. For illustration, we have fixed $m_{max}=5$ so that
the $\mathbf{B}\rightarrow 0$ energies for remote magnetic subbands, i.e.
$\pm\gamma$, are tractable.

Having the closed form for the singular values, we hereby reformulate the
finite $\mathbf{B}$ Hamiltonian (Eq.(17)) in terms of the modes that are
allowed to hybridize. The analysis for $\tau=+1$ is presented below, for
$\tau=-1$ in SM-F.1. Using the SVD decomposition discussed we can re-write the
matrix elements as

$\displaystyle
h^{1}_{[1m0],[1r^{\prime}]}c^{\dagger}_{11k0m}f_{11kr^{\prime}}=\gamma\Sigma_{m^{\prime},\bar{r}}c^{\dagger}_{11k0m}U_{mm^{\prime}}$
$\displaystyle\times V_{\bar{r},r^{\prime}}f_{11kr^{\prime}}$
$\displaystyle=\sum_{m=0}^{m_{1,+1}}\sum_{\bar{r}=0}^{m_{1,+1}}\gamma\Sigma_{m}\delta_{m\bar{r}}c^{\dagger}_{11k0m}\bar{f}_{11k\bar{r}},$
(54)

where we used the fact that $U$ is an identity matrix and the rectangular
matrix $\Sigma_{m,r^{\prime}}$ is non-zero only along its main diagonal.

$\displaystyle
h^{1}_{[2m0],[1r^{\prime}]}c^{\dagger}_{21k0m}f_{11kr^{\prime}}=-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\times$
$\displaystyle\sum_{m=0}^{m_{2,+1}}\sqrt{m+1}c^{\dagger}_{21k0m}\Sigma_{m+1,\bar{r}}V_{\bar{r},r^{\prime}}f_{11kr^{\prime}}=$
$\displaystyle-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\sum_{m=0}^{m_{2,+1}}\sum_{\bar{r}=0}^{m_{2,+1}+1}\sqrt{m+1}\Sigma_{m+1}\delta_{m+1,\bar{r}}c^{\dagger}_{21k0m}\bar{f}_{11k\bar{r}}.$
(55)

Because $m_{2,+1}+1=m_{1,+1}=m_{max}+1$, the upper bound on the $\bar{r}$
summation is the same in (54) and (55). Similarly

$\displaystyle
h^{1}_{[1m0],[2r^{\prime}]}c^{\dagger}_{11k0m}f_{21kr^{\prime}}=$
$\displaystyle
i\sqrt{2}\sqrt{m}\frac{v^{\prime}_{\ast}}{\ell}\sum_{m=0}^{m_{1,+1}}\sum_{\bar{r}=0}^{m_{1,+1}-1}\Sigma_{\bar{r}}\delta_{m-1,\bar{r}}c^{\dagger}_{11k0m}\bar{f}_{21k\bar{r}},$
(56) $\displaystyle
h^{1}_{[2m0],[2r^{\prime}]}c^{\dagger}_{21k0m}f_{21kr^{\prime}}=$
$\displaystyle\gamma\sum_{m=0}^{m_{2,+1}}\sum_{\bar{r}=0}^{m_{2,+1}}\Sigma_{m}\delta_{m,\bar{r}}c^{\dagger}_{21k0m}\bar{f}_{21k\bar{r}}.$
(57)

Because $m_{1,+1}-1=m_{2,+1}=m_{max}$, the upper bound on $\bar{r}$ summation
is the same in Eq.(56) and Eq.(57). Thus for a given $k$, out of the available
$2q$ $f-$modes, only $2m_{max}+3$ couple; recall that $m_{max}\lesssim q/2$.
The mean-field interactions for $f-$ fermion modes in new $\bar{f}$ basis can
be given as

$V^{f,\tau=+1}=\sum_{k\in[0,1)\otimes[0,\frac{1}{q})}V^{f,\tau=+1}_{coupled}+V^{f,\tau=+1}_{decoupled},$
(58)
$V^{f,\tau=+1}_{coupled}=-\frac{U_{1}}{2}\left(\sum_{m=0}^{m_{max}+1}\bar{f}^{\dagger}_{11km}\bar{f}_{11km}+\sum_{m=0}^{m_{max}}\bar{f}^{\dagger}_{21km}\bar{f}_{21km}\right)$
(59)
$V^{f,\tau=+1}_{decoupled}=-\frac{U_{1}}{2}\sum_{b=1}^{2}\sum_{m^{\prime}=m_{max}+\bar{b}}^{q-1}\bar{f}^{\dagger}_{b1km^{\prime}}\bar{f}_{b1km^{\prime}},$
(60)

where $\bar{1}(\bar{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

$H^{\tau=+1}_{coupled}=\sum_{k}\sum_{\alpha,\alpha^{\prime}=1}^{6}\sum_{m=0}^{m_{\alpha}}\sum_{m^{\prime}=0}^{m_{\alpha^{\prime}}}\Xi_{m\alpha,m^{\prime}\alpha^{\prime}}d^{\dagger}_{m\alpha}(k)d_{m^{\prime}\alpha^{\prime}}(k),$
(61)

where $m_{\alpha=1,\ldots,4}=m_{\alpha,+1}$, $m_{5}=m_{max}+1$ and
$m_{6}=m_{max}$, and

$\displaystyle
d^{\dagger}_{m\alpha}(k)=\left(c^{\dagger}_{11k0m},c^{\dagger}_{21k0m},c^{\dagger}_{31k0m},c^{\dagger}_{41k0m},\bar{f}^{\dagger}_{11km},\bar{f}^{\dagger}_{21km}\right)_{\alpha}.$
(62)

We now define an operator

$\displaystyle\hat{h}^{+1}_{\alpha,\alpha^{\prime}}=\left(\begin{array}[]{cccccc}0&0&-i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}&0&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\hat{a}^{\dagger}\Sigma(\hat{a}^{\dagger}\hat{a})\\\
0&0&0&i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}^{\dagger}&-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\hat{a}\Sigma(\hat{a}^{\dagger}\hat{a})&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})\\\
i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}^{\dagger}&0&-\frac{J}{2}&M&0&0\\\
0&-i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}&M&-\frac{J}{2}&0&0\\\
\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\Sigma(\hat{a}^{\dagger}\hat{a})\hat{a}^{\dagger}&0&0&-\frac{U_{1}}{2}&0\\\
-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\Sigma(\hat{a}^{\dagger}\hat{a})\hat{a}&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&0&0&0&-\frac{U_{1}}{2}\end{array}\right)_{\alpha,\alpha^{\prime}},$
(69)

where $\hat{a}$ is a simple h.o. lowering operator in terms of which matrix
$\Xi_{m\alpha,m^{\prime}\alpha^{\prime}}$ can be expressed as

$\Xi_{m\alpha,m^{\prime}\alpha^{\prime}}=\langle
m|\hat{h}^{+1}_{\alpha,\alpha^{\prime}}|m^{\prime}\rangle.$ (70)

Here $|m\rangle$ is a simple h.o. eigenstate and $\Sigma(m)=\Sigma_{m}$.

Note that a naive minimal substitution into Eq.(II), with $\lambda=0$ would
reproduce Eq.(69), except the singular values would get replaced by $1$, which
corresponds to the $\mathbf{B}\rightarrow 0$ limit of our analysis. The
decoupling of $2-(2m_{max}+3)/q$ modes per moiré unit cell per spin and the
fairly strong $\mathbf{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
$\mathbf{B}$ in the narrow band strong coupling window set by
$\left(\frac{J}{2},\frac{U_{1}}{2}\right)$, we can explain the state counting
analytically by presenting an exact solution to the operator
$\hat{h}_{\alpha,\alpha^{\prime}}$ in the flat band limit $M=0$.

###  III.2 Analytical solution for flat band, U(4) symmetric, THFM in
$\mathbf{B}$

For $M=0$, the exact solutions to the eigenstates of the operator in Eq.(69)
are presented below. The $\mathbf{B}$ field independent $-J/2$ Landau level
energy shown in Fig.(1) comes from the anomalous $c$-mode

$\theta_{1}=\left[0,0,\ket{0},0,0,0\right]^{T}.$ (71)

The rest of the problem can be solved using the following ansätze:

$\displaystyle\theta_{3}$ $\displaystyle=$
$\displaystyle\left[c^{(3)}_{1}\ket{0},0,c^{(3)}_{3}\ket{1},0,c^{(3)}_{5}\ket{0},0\right]^{T},$
(72) $\displaystyle\theta_{5}$ $\displaystyle=$
$\displaystyle\left[c^{(5)}_{1}\ket{1},c^{(5)}_{2}\ket{0},c^{(5)}_{3}\ket{2},0,c^{(5)}_{5}\ket{1},c^{(5)}_{6}\ket{0}\right]^{T},$
$\displaystyle\theta_{6_{m}}$ $\displaystyle=$
$\displaystyle\left[c^{(6_{m})}_{1}\ket{m},c^{(6_{m})}_{2}\ket{m-1},c^{(6_{m})}_{3}\ket{m+1},\right.$
(74)
$\displaystyle\left.c^{(6_{m})}_{4}\ket{m-2},c^{(6_{m})}_{5}\ket{m},c^{(6_{m})}_{6}\ket{m-1}\right]^{T},$

where $m\in\\{2,\ldots,m_{max}+1\\}$. $c^{(\beta)}_{\alpha}$ denotes the
coefficient of corresponding h.o state at index $\alpha$ in the 6-component
spinor in Eq.(62) and $\beta$ labels the ansatz index $\theta_{\beta}$. Using
the above, we can set up the eigen-equation and solve for the corresponding
coefficients.

The ansätze $\theta_{3}$ and $\theta_{5}$ yield the $3\times 3$ and $5\times
5$ Hermitian matrices, whose eigenvectors are $c^{(3)}_{\alpha}$ and
$c^{(5)}_{\alpha}$, respectively:

$\displaystyle
h^{+1}_{3}=\left(\begin{array}[]{ccc}0&-i\frac{\sqrt{2}v_{*}}{\ell}&\gamma\Sigma_{0}\\\
i\frac{\sqrt{2}v_{*}}{\ell}&-\frac{J}{2}&0\\\
\gamma\Sigma_{0}&0&-\frac{U_{1}}{2}\end{array}\right),$ (78)

$h^{+1}_{5}=\left(\begin{array}[]{ccccc}0&0&-i\frac{2v_{*}}{\ell}&\gamma\Sigma_{1}&i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{0}\\\
0&0&0&-i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{1}&\gamma\Sigma_{0}\\\
i\frac{2v_{*}}{\ell}&0&-\frac{J}{2}&0&0\\\
\gamma\Sigma_{1}&i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{1}&0&-\frac{U_{1}}{2}&0\\\
-i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{0}&\gamma\Sigma_{0}&0&0&-\frac{U_{1}}{2}\end{array}\right).$
(79)

Similarly, the ansatz $\theta^{m}_{6}$ yields the following $6\times 6$
Hermitian matrix for each $m$, whose eigenvectors are $c^{(6,m)}_{\alpha}$:

$\displaystyle h^{+1,m}_{6}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}0&0&-i\sqrt{2m+2}\frac{v_{*}}{\ell}&0&\gamma\Sigma_{m}&i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m-1}\\\
0&0&0&i\sqrt{2m-2}\frac{v_{*}}{\ell}&-i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m}&\gamma\Sigma_{m-1}\\\
+i\sqrt{2m+2}\frac{v_{*}}{\ell}&0&-\frac{J}{2}&0&0&0\\\
0&-i\sqrt{2m-2}\frac{v_{*}}{\ell}&0&-\frac{J}{2}&0&0\\\
\gamma\Sigma_{m}&i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m}&0&0&-\frac{U_{1}}{2}&0\\\
-i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m-1}&\gamma\Sigma_{m-1}&0&0&0&-\frac{U_{1}}{2}\end{array}\right).$
(86)

The strong coupling spectrum obtained using the above matrices is shown in
Fig.(1), which also includes the spectrum at $\tau=-1$. One way to obtain the
spectrum at $\tau=-1$ is by replacing $J,U\rightarrow-J,-U$ in the above
decoupled matrices(see SM-F.1). Also note that the decoupled $f-$modes and the
anomalous $c-$mode at $\tau=-1$ then form the field independent
$\frac{U_{1}}{2}$ and $\frac{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 $\mathbf{B}\rightarrow 0$
energies for $J=U_{1}=0$ are $0$ for the flat band, and $\pm\gamma$ for the
remote bands, because $\Sigma_{m}\rightarrow 1$. Note that these are the non-
interacting zero-field energies at $\Gamma$ point in mBZ. For any non-zero
$\mathbf{B}$, the total zero mode count can be understood as follows:
$h^{\tau,m}_{6}$ being bipartite contributes two zero modes (even if we do not
set $\Sigma_{m}=1$). Including all $m_{max}$ such matrices gives $2m_{max}$
zero modes. $h^{\tau}_{3}$, $h^{\tau}_{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 $\mathbf{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 $\mathbf{B}$ because the
remote bands emanate out of $\pm\gamma$ at $\mathbf{B}=0$.

For $J\neq 0$ and $U_{1}\neq 0$, at valley $\tau$, the $2q-(2m_{max}+3)$ of
the decoupled $f$’s have the energy $-\tau U_{1}/2$. The $2m_{max}+2$ modes
“cluster” around the $\mathbf{B}\rightarrow 0$ energies of the decoupled
blocks: $-\tau\frac{J}{2}$ and
$-\tau\frac{U_{1}}{4}\pm\sqrt{\frac{U^{2}_{1}}{16}+\gamma^{2}}$ (marked by
$\pm\mathcal{E}_{\mp\tau}$ in Fig.(2)), where all three are singly degenerate
for $h^{\tau}_{3}$, last two are doubly degenerate for $h^{\tau}_{5}$ and all
three are doubly degenerate for $h^{\tau}_{6}$. Including the anomalous $-\tau
J/2$ level (Eq.(71)) and the $2q-(2m_{max}+3)$ decoupled $f$ modes at
$-\tau\frac{U_{1}}{2}$, there are $2q$ modes in total at each $k$ within the
strong coupling spectrum at each valley, corresponding to $2$ states per moiré
unit cell per spin and consistent with the projected strong coupling results
[82, 83].

###  III.3 Integer Fillings $\nu\mathbf{=\pm 1,\pm 2}$

The above analysis can straightforwardly be applied to other integer fillings.
In this section we illustrate it by extending the formalism to nonzero
fillings $\nu=\pm 1,\pm 2$. Since the mean field interactions are spin
dependent for integer fillings, we introduce the spin quantum number label,
$s=\uparrow\downarrow$, to the $c$ and $f$ fermion annihilation operators as
$c_{a\tau krms}$ and $f_{b\tau kr^{\prime}s}$ in order to differentiate
between the spin-up and spin-down sectors.

Figure 7: Illustrative plot at $m_{max}=5$, $w_{0}/w_{1}=0.7$, valley
$\mathbf{K}$ and spin $\downarrow$ in flat band limit, $M=0$, for $\nu=-1$,
showing (a) $m_{max}=5$ magnetic subbands emanating out of
$-W_{3}+\frac{J}{2}$ for Chern $-1$ band and (b) $m_{max}+3=8$ magnetic
subbands emanating out of energy $-W_{3}-\frac{J}{2}$ for the Chern $+1$ band,
which includes the field independent level $-W_{3}-\frac{J}{2}$ of decoupled
$c$ mode. Including the $q-(m_{max}+1)=q-6$ decoupled $f$ modes with energy
$-\frac{1}{2}U_{1}-6U_{2}$ for the Chern $-1$ band and $q-(m_{max}+2)=q-7$
decoupled $f$ modes with energy $-\frac{3}{2}U_{1}-6U_{2}$ for the Chern $+1$
band, we have in total $(q-6)+5=q-1$ and $(q-7)+8=q+1$ magnetic subbands for
Chern $\mp 1$ bands respectively as we should.

For $\nu=-1$, we discuss the sector valley $\mathbf{K}$ spin $\downarrow$ for
the VP state [56], wherein at zero-field, the charge $\pm 1$ excitations
occupy well gapped Chern $\mp 1$ states marked by red in Fig(8a). The other
spin-valley sectors at $\nu=-1$, namely valley $\mathbf{K}$ spin $\uparrow$
and valley $\mathbf{K}^{\prime}$ spin $\uparrow\downarrow$(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 $\mathbf{B}$, the Chern $\pm 1$ bands
gain and lose states respectively [88], with a total state count pinned to $2$
per moiré unit cell. The spectrum at $\nu=+1$ is related by particle-hole
symmetry. The interactions at valley $\mathbf{K}$ and spin $\downarrow$ read
as

$\displaystyle V^{\tau=+1,\downarrow}_{\nu=\pm
1}=\nu\sum_{k}\bigg{(}\sum_{a=1}^{4}\sum_{m=0}^{m_{a,+1}}\sum_{r=0}^{p-1}W_{a}c^{\dagger}_{a1krm\downarrow}c_{a1krm\downarrow}$
$\displaystyle+\sum_{a=3,4}\sum_{m=0}^{m_{a,+1}}\sum_{r=0}^{p-1}(-1)^{a+1}\frac{J}{2}c^{\dagger}_{a1krm\downarrow}c_{a1krm\downarrow}$
$\displaystyle+\sum_{b=1,2}\sum_{r^{\prime}=0}^{q-1}\left(\frac{2+(-1)^{b+1}}{2}U_{1}+6U_{2}\right)f^{\dagger}_{b1kr^{\prime}\downarrow}f_{b1kr^{\prime}\downarrow}\bigg{)}.$
(87)

Here $W_{a\in\\{1\ldots 4\\}}$ and $U_{2}$ are mean field coefficients with
$W_{1}=W_{2}$ and $W_{3}=W_{4}$ [56]. In the $\bar{f}$ basis we can re-write
the interaction for $f$-fermion modes as

$\displaystyle V^{f,\tau=+1,\downarrow}_{\nu=\pm 1}$ $\displaystyle=$
$\displaystyle\sum_{k\in[0,1)\otimes[0,\frac{1}{q})}V^{f,\tau=+1,\downarrow,\nu=\pm
1}_{coupled}$ (88) $\displaystyle+V^{f,\tau=+1,\downarrow,\nu=\pm
1}_{decoupled}$

where

$\displaystyle V^{f,\tau=+1,\downarrow,\nu=\pm 1}_{coupled}=$
$\displaystyle\nu\left(\frac{3}{2}U_{1}+6U_{2}\right)\left(\sum_{m=0}^{m_{max}+1}\bar{f}^{\dagger}_{11km\downarrow}\bar{f}_{11km\downarrow}\right)$
(89)
$\displaystyle+\nu\left(\frac{1}{2}U_{1}+6U_{2}\right)\left(\sum_{m=0}^{m_{max}}\bar{f}^{\dagger}_{21kms}\bar{f}_{21km\downarrow}\right),$
(90) $\displaystyle V^{f,\tau=+1,s,\nu=\pm
2}_{decoupled}=\nu\sum_{b=1}^{2}\sum_{m^{\prime}=m_{max}+\bar{b}}^{q-1}$
$\displaystyle\left(\frac{2+(-1)^{b+1}}{2}U_{1}+6U_{2}\right)\bar{f}^{\dagger}_{b1km^{\prime}\downarrow}\bar{f}_{b1km^{\prime}\downarrow},$
(91)

where $\bar{1}(\bar{2})=2,1$. Note that out of the available $2q$ $f$ modes,
$q-(m_{max}+2)$ are decoupled with energy
$\nu\left(\frac{3}{2}U_{1}+6U_{2}\right)$ and $q-(m_{max}+1)$ are decoupled
with energy $\nu\left(\frac{1}{2}U_{1}+6U_{2}\right)$, 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

$H^{\tau=+1,\downarrow,\nu=\pm
1}_{coupled}=\sum_{k}\sum_{\alpha,\alpha^{\prime}=1}^{6}\sum_{m=0}^{m_{\alpha}}\sum_{m^{\prime}=0}^{m_{\alpha^{\prime}}}\Xi^{\downarrow,\nu=\pm
1}_{m\alpha,m^{\prime}\alpha^{\prime}}d^{\dagger}_{m\alpha\downarrow}(k)d_{m^{\prime}\alpha^{\prime}\downarrow}(k),$
(92)

where $m_{\alpha=1,\ldots,4}=m_{\alpha,+1}$, $m_{5}=m_{max}+1$ and
$m_{6}=m_{max}$, and

$d^{\dagger}_{m\alpha\downarrow}(k)=\left(c^{\dagger}_{11k0m\downarrow},c^{\dagger}_{21k0m\downarrow},c^{\dagger}_{31k0m\downarrow},c^{\dagger}_{41k0m\downarrow},\bar{f}^{\dagger}_{11km\downarrow},\bar{f}^{\dagger}_{21km\downarrow}\right)_{\alpha},$
(93)

with

$\displaystyle\Xi^{\downarrow,\nu=\pm
1}_{m\alpha,m^{\prime}\alpha^{\prime}}=\langle
m|\hat{h}^{+1,\downarrow,\nu=\pm
1}_{\alpha,\alpha^{\prime}}|m^{\prime}\rangle,$ (94)

where the operators $\hat{h}^{+1,\downarrow,\nu=\pm
1}_{\alpha,\alpha^{\prime}}$ are given as

$\displaystyle\hat{h}^{+1,\downarrow,\nu=\pm
1}_{\alpha,\alpha^{\prime}}=\left(\begin{array}[]{cccccc}\nu
W_{1}&0&-i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}&0&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\hat{a}^{\dagger}\Sigma(\hat{a}^{\dagger}\hat{a})\\\
0&\nu
W_{1}&0&i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}^{\dagger}&-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\hat{a}\Sigma(\hat{a}^{\dagger}\hat{a})&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})\\\
i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}^{\dagger}&0&\nu(W_{3}+\frac{J}{2})&M&0&0\\\
0&-i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}&M&\nu(W_{3}-\frac{J}{2})&0&0\\\
\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\Sigma(\hat{a}^{\dagger}\hat{a})\hat{a}^{\dagger}&0&0&\nu(\frac{3}{2}U_{1}+6U_{2})&0\\\
-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\Sigma(\hat{a}^{\dagger}\hat{a})\hat{a}&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&0&0&0&\nu(\frac{1}{2}U_{1}+6U_{2})\end{array}\right)_{\alpha,\alpha^{\prime}}.$
(101)

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.(71) forms the field independent $\nu(W_{3}+\frac{J}{2})$
level. The remaining eigenstates can be obtained using the ansätze given in
Eqs.(72)-(74). Setting up eigen-equation for these ansätze yield us the
following $3\times 3$, $5\times 5$ and $m_{max}$ $6\times 6$ matrices:

$\displaystyle h^{+1,\nu=\pm 1}_{3}=\left(\begin{array}[]{ccc}\nu
W_{1}&-i\frac{\sqrt{2}v_{*}}{\ell}&\gamma\Sigma_{0}\\\
i\frac{\sqrt{2}v_{*}}{\ell}&\nu(W_{3}+\frac{J}{2})&0\\\
\gamma\Sigma_{0}&0&\nu(\frac{3}{2}U_{1}+6U_{2})\end{array}\right),$ (105)
$h^{+1,\nu=\pm 1}_{5}=\left(\begin{array}[]{ccccc}\nu
W_{1}&0&-i\frac{2v_{*}}{\ell}&\gamma\Sigma_{1}&i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{0}\\\
0&\nu
W_{1}&0&-i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{1}&\gamma\Sigma_{0}\\\
i\frac{2v_{*}}{\ell}&0&\nu(W_{3}+\frac{J}{2})&0&0\\\
\gamma\Sigma_{1}&i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{1}&0&\nu(\frac{3U_{1}}{2}+6U_{2})&0\\\
-i\frac{\sqrt{2}v_{*}^{\prime}}{\ell}\Sigma_{0}&\gamma\Sigma_{0}&0&0&\nu(\frac{U_{1}}{2}+6U_{2})\end{array}\right),$
(106) $\displaystyle h^{+1,m,\nu=\pm 1}_{6}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}\nu
W_{1}&0&-i\sqrt{2m+2}\frac{v_{*}}{\ell}&0&\gamma\Sigma_{m}&i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m-1}\\\
0&\nu
W1&0&i\sqrt{2m-2}\frac{v_{*}}{\ell}&-i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m}&\gamma\Sigma_{m-1}\\\
+i\sqrt{2m+2}\frac{v_{*}}{\ell}&0&\nu(W_{3}+\frac{J}{2})&0&0&0\\\
0&-i\sqrt{2m-2}\frac{v_{*}}{\ell}&0&\nu(W_{3}-\frac{J}{2})&0&0\\\
\gamma\Sigma_{m}&i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m}&0&0&\nu(\frac{3U_{1}}{2}+6U_{2})&0\\\
-i\sqrt{2m}\frac{v_{*}^{\prime}}{\ell}\Sigma_{m-1}&\gamma\Sigma_{m-1}&0&0&0&\nu(\frac{U_{1}}{2}+6U_{2})\end{array}\right),$
(113)

where $m\in\\{2\ldots m_{max}+1\\}$. The magnetic subbands contributed by the
coupled modes for Chern $+1$ band emanate out of $\mathbf{B}\rightarrow 0$
energy $\nu(W_{3}+\frac{J}{2})$ of the above decoupled matrices, which is
$m_{max}$ fold degenerate for matrix Eq.(113) and singly degenerate for
matrices Eq.(105) and Eq.(106). Including the field independent level
$\nu(W_{3}+\frac{J}{2})$ of the decoupled $c$ mode, we have in total
$m_{max}+3$ magnetic subbands emanating out of $\nu(W_{3}+\frac{J}{2})$ for
Chern $+1$ band. Now recall that we have $q-(m_{max}+2)$ decoupled $f$ modes
with energy $\nu(\frac{3U_{1}}{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 $\mathbf{B}\rightarrow 0$ energy
$\nu(W_{3}-\frac{J}{2})$ of matrix Eq.(113). Including the $q-(m_{max}+1)$
decoupled $f$ modes with energy $\nu(\frac{U_{1}}{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 $\pm 1$ band offers in total $q\pm
1$ magnetic subbands or $1\pm\frac{1}{q}=1\pm\frac{\phi}{\phi_{0}}$ states per
moiré unit cell as expected. The full spectrum in flat band limit is shown in
Fig.(8). Note that the $\mathbf{B}\rightarrow 0$ energy of the decoupled
matrices correspond to the zero field THFM energies at $\Gamma$ in mBZ.

The mean field interaction at filling $\nu=\pm 2$ for VP state at $\tau=+1$
for spin sector $s$ reads[56]

$\displaystyle V^{\tau=+1,s}_{\nu=\pm
2}=\nu\sum_{k}\bigg{(}\sum_{a=1}^{4}\sum_{m=0}^{m_{a,+1}}\sum_{r=0}^{p-1}W_{a}c^{\dagger}_{a1krms}c_{a1krms}$
$\displaystyle+\sigma_{s}\sum_{a=3,4}\sum_{m=0}^{m_{a,+1}}\sum_{r=0}^{p-1}\frac{J}{4}c^{\dagger}_{a1krms}c_{a1krms}$
$\displaystyle+\sum_{b=1,2}\sum_{r^{\prime}=0}^{q-1}\left(\frac{4+\sigma_{s}}{4}U_{1}+6U_{2}\right)f^{\dagger}_{b1kr^{\prime}s}f_{b1kr^{\prime}s}\bigg{)},$
(114)

where $\sigma_{s}=\pm 1$ for $s=\uparrow\downarrow$ respectively. In the
$\bar{f}$ basis we can re-write the interaction for $f$-fermion modes as

$\displaystyle V^{f,\tau=+1,s}_{\nu=\pm 2}$ $\displaystyle=$
$\displaystyle\sum_{k\in[0,1)\otimes[0,\frac{1}{q})}V^{f,\tau=+1,s,\nu=\pm
2}_{coupled}$ (115) $\displaystyle+V^{f,\tau=+1,s,\nu=\pm 2}_{decoupled}$

where

$\displaystyle V^{f,\tau=+1,s,\nu=\pm
2}_{coupled}=\nu\left(\frac{4+\sigma_{s}}{4}U_{1}+6U_{2}\right)$
$\displaystyle\times\left(\sum_{m=0}^{m_{max}+1}\bar{f}^{\dagger}_{11kms}\bar{f}_{11kms}+\sum_{m=0}^{m_{max}}\bar{f}^{\dagger}_{21kms}\bar{f}_{21kms}\right),$
(116) $\displaystyle V^{f,\tau=+1,s,\nu=\pm
2}_{decoupled}=\nu\left(\frac{4+\sigma_{s}}{4}U_{1}+6U_{2}\right)$
$\displaystyle\times\sum_{b=1}^{2}\sum_{m^{\prime}=m_{max}+\bar{b}}^{q-1}\bar{f}^{\dagger}_{b1km^{\prime}s}\bar{f}_{b1km^{\prime}s},$
(117)

where $\bar{1}(\bar{2})=2,1$. Note that we still have $2-(2m_{max}+3)/q$
decoupled $f$ states per moiré unit cell per spin. The coupled modes can then
be described by

$H^{\tau=+1,s,\nu=\pm
2}_{coupled}=\sum_{k}\sum_{\alpha,\alpha^{\prime}=1}^{6}\sum_{m=0}^{m_{\alpha}}\sum_{m^{\prime}=0}^{m_{\alpha^{\prime}}}\Xi^{s,\nu=\pm
2}_{m\alpha,m^{\prime}\alpha^{\prime}}d^{\dagger}_{m\alpha
s}(k)d_{m^{\prime}\alpha^{\prime}s}(k),$ (118)

where $m_{\alpha=1,\ldots,4}=m_{\alpha,+1}$, $m_{5}=m_{max}+1$ and
$m_{6}=m_{max}$, and

$d^{\dagger}_{m\alpha
s}(k)=\left(c^{\dagger}_{11k0ms},c^{\dagger}_{21k0ms},c^{\dagger}_{31k0ms},c^{\dagger}_{41k0ms},\bar{f}^{\dagger}_{11kms},\bar{f}^{\dagger}_{21kms}\right)_{\alpha}.$
(119)

with

$\displaystyle\Xi^{s,\nu=\pm 2}_{m\alpha,m^{\prime}\alpha^{\prime}}=\langle
m|\hat{h}^{+1,s,\nu=\pm 2}_{\alpha,\alpha^{\prime}}|m^{\prime}\rangle.$ (120)

where the operators $\hat{h}^{+1,s,\nu=\pm 2}_{\alpha,\alpha^{\prime}}$ are
given as

$\displaystyle\hat{h}^{+1,s,\nu=\pm
2}_{\alpha,\alpha^{\prime}}=\left(\begin{array}[]{cccccc}\nu
W_{1}&0&-i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}&0&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\hat{a}^{\dagger}\Sigma(\hat{a}^{\dagger}\hat{a})\\\
0&\nu
W_{1}&0&i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}^{\dagger}&-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\hat{a}\Sigma(\hat{a}^{\dagger}\hat{a})&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})\\\
i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}^{\dagger}&0&\nu(W_{3}+\sigma_{s}\frac{J}{4})&M&0&0\\\
0&-i\sqrt{2}\frac{v_{\ast}}{\ell}\hat{a}&M&\nu(W_{3}+\sigma_{s}\frac{J}{4})&0&0\\\
\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\Sigma(\hat{a}^{\dagger}\hat{a})\hat{a}^{\dagger}&0&0&\nu(\frac{4+\sigma_{s}}{4}U_{1}+6U_{2})&0\\\
-i\sqrt{2}\frac{v^{\prime}_{\ast}}{\ell}\Sigma(\hat{a}^{\dagger}\hat{a})\hat{a}&\gamma\Sigma(\hat{a}^{\dagger}\hat{a})&0&0&0&\nu(\frac{4+\sigma_{s}}{4}U_{1}+6U_{2})\end{array}\right)_{\alpha,\alpha^{\prime}}.$
(127)

Figure 8: Interacting heavy fermion Hofstadter spectra for (b) Valley
$\mathbf{K}$ spin ${\color[rgb]{1,0,0}\downarrow}$, (c)Valley $\mathbf{K}$
spin ${\color[rgb]{0,0,1}\uparrow}$ and (d)Valley $\mathbf{K}^{\prime}$ spin
${\color[rgb]{0,1,0}\uparrow}$ contrasted with (a) zero field spectrum at
filling $\nu=-1$ at $w_{0}/w_{1}=0.7$, with parameters $W_{1}=44.05meV$,
$W_{3}=49.33meV$, $U_{2}=2.656meV$ in the flat band limit $M=0$ with
$m_{max}=\lceil\frac{q-3}{2}\rceil$. The color on panel-(a) labels the spin
and valley sector, while color labelling for panels-(b,c,d) denotes the
f-character.

For $M=0$, the exact eigenstates for the above operators in Eq-(127) can be
solved using the same ansätze given in Eqs.(72)-(74)(see SM-H). The spectrum
is shown in Fig.(9). The decoupled $c$ state given by Eq.(71), forms the field
independent $\pm(2W_{3}+\frac{J}{2})$ and $\pm(2W_{3}-\frac{J}{2})$ levels at
fillings $\nu=\pm 2$ and $\tau=+1$ for spin $\uparrow\downarrow$ sectors
respectively. However the level $\pm(2W_{3}+\frac{J}{2})$ is singly degenerate
as it is not the part of spectrum at $\tau=-1$(see SM-F.3) and the total state
count is pinned to $2$ per moiré unit cell per valley per spin for the narrow
bands(see SM-H).

Figure 9: Interacting heavy fermion Hofstadter spectra for (b) spin sector up
($\uparrow$) and (c) spin sector down ($\downarrow$) contrasted with the (a)
zero-field spectrum at filling $\nu=-2$ at $w_{0}/w_{1}=0.7$, valley
$\mathbf{K}$ in the flat band limit $M=0$ with
$m_{max}=\lceil\frac{q-3}{2}\rceil$. The color on panel-(a) labels the spin
sector, while color labelling for panels-(b,c) denotes the f-character.

##  IV Discussion

We have put forward a generalization of THFM in finite B. Although the
formalism applies to any rational value of $\frac{\Phi}{\Phi_{0}}$, the
physical nature of hybridization amidst the heavy $f$ and topological $c$
fermions is particularly revealing for the $\frac{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 $\frac{1}{q}$ of the primitive reciprocal lattice
momentum $\mathbf{g}_{2}$ causes a shift by a moiré unit cell. Because the
momentum is defined only modulo $\mathbf{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 $c$-$f$ hybridization along the $\frac{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 $\mathbf{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 $\pm\frac{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 $\frac{U_{1}}{2}$ around $q\approx
80$, because the theory starts “losing” LLs with increasing field, reminiscent
of the results obtained in [82]. Although there are $2-\frac{2m_{max}+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\neq 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 other integer fillings,
candidate groundstates [90, 91], and Hofstadter-scale fluxes where reentrant
many-body and topological effects are at play [78, 92, 93, 49, 94, 95, 96, 97,
98, 99], is also left for the future work.

###### Acknowledgements.

The authors thank Xiaoyu Wang and Zhi-da Song for valuable conversations, and
to Dumitru Călugăru for computational advice. J. H.-A. is supported by a Hertz
Fellowship and by ONR Grant No. N00014- 20-1-2303. B.A.B is supported by the
DOE Grant No. DE-SC0016239 and by the EPiQS Initiative, Grant GBMF11070. A.C.
was supported by Grant No. GBMF8685 towards the Princeton theory program and
by the Gordon and Betty Moore Foundation through the EPiQS Initiative, Grant
GBMF11070. Further sabbatical support for A.C., J.H.A and B.A.B was provided
by the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No. 101020833), the Schmidt
Fund for Innovative Research, Simons Investigator Grant No. 404513. O.V. was
supported by NSF Grant No. DMR-1916958 and is partially funded by the Gordon
and Betty Moore Foundation’s EPiQS Initiative Grant GBMF11070, National High
Magnetic Field Laboratory through NSF Grant No. DMR-1157490 and the State of
Florida.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018b).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Intrinsic quantized anomalous hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * Xu and Balents [2018] C. Xu and L. Balents, Topological superconductivity in twisted multilayer graphene, Phys. Rev. Lett. 121, 087001 (2018).
  * Po  _et al._ [2018a] H. C. Po, H. Watanabe, and A. Vishwanath, Fragile topology and wannier obstructions, Phys. Rev. Lett. 121, 126402 (2018a).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions, Phys. Rev. B 98, 085435 (2018).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Phys. Rev. X 8, 031087 (2018).
  * Po  _et al._ [2018b] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of mott insulating behavior and superconductivity in twisted bilayer graphene, Phys. Rev. X 8, 031089 (2018b).
  * Isobe  _et al._ [2018] H. Isobe, N. F. Q. Yuan, and L. Fu, Unconventional superconductivity and density waves in twisted bilayer graphene, Phys. Rev. X 8, 041041 (2018).
  * Wu  _et al._ [2018] F. Wu, A. H. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene, Phys. Rev. Lett. 121, 257001 (2018).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong coupling phases of partially filled twisted bilayer graphene narrow bands, Phys. Rev. Lett. 122, 246401 (2019).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of magic angles in twisted bilayer graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Lian  _et al._ [2019] B. Lian, Z. Wang, and B. A. Bernevig, Twisted bilayer graphene: A phonon-driven superconductor, Phys. Rev. Lett. 122, 257002 (2019).
  * Wu and Das Sarma [2020] F. Wu and S. Das Sarma, Collective excitations of quantum anomalous hall ferromagnets in twisted bilayer graphene, Phys. Rev. Lett. 124, 046403 (2020).
  * Zhang  _et al._ [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Correlated insulating phases of twisted bilayer graphene at commensurate filling fractions: A hartree-fock study, Phys. Rev. B 102, 035136 (2020).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. 124, 097601 (2020).
  * Cea and Guinea [2020] T. Cea and F. Guinea, Band structure and insulating states driven by coulomb interaction in twisted bilayer graphene, Phys. Rev. B 102, 045107 (2020).
  * Kang and Vafek [2020] J. Kang and O. Vafek, Non-abelian dirac node braiding and near-degeneracy of correlated phases at odd integer filling in magic-angle twisted bilayer graphene, Phys. Rev. B 102, 035161 (2020).
  * Abouelkomsan  _et al._ [2020] A. Abouelkomsan, Z. Liu, and E. J. Bergholtz, Particle-hole duality, emergent fermi liquids, and fractional chern insulators in moiré flatbands, Phys. Rev. Lett. 124, 106803 (2020).
  * Bultinck  _et al._ [2020a] N. Bultinck, S. Chatterjee, and M. P. Zaletel, Mechanism for anomalous hall ferromagnetism in twisted bilayer graphene, Phys. Rev. Lett. 124, 166601 (2020a).
  * Soejima  _et al._ [2020] T. Soejima, D. E. Parker, N. Bultinck, J. Hauschild, and M. P. Zaletel, Efficient simulation of moiré materials using the density matrix renormalization group, Phys. Rev. B 102, 205111 (2020).
  * Ledwith  _et al._ [2020] P. J. Ledwith, G. Tarnopolsky, E. Khalaf, and A. Vishwanath, Fractional chern insulator states in twisted bilayer graphene: An analytical approach, Phys. Rev. Res. 2, 023237 (2020).
  * Repellin  _et al._ [2020] C. Repellin, Z. Dong, Y.-H. Zhang, and T. Senthil, Ferromagnetism in narrow bands of moiré superlattices, Phys. Rev. Lett. 124, 187601 (2020).
  * Bultinck  _et al._ [2020b] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Phys. Rev. X 10, 031034 (2020b).
  * Minhal Gardezi  _et al._ [2020] S. Minhal Gardezi, H. Pirie, W. Dorrell, and J. E. Hoffman, Acoustic twisted bilayer graphene, arXiv e-prints , arXiv:2010.10037 (2020), arXiv:2010.10037 [cond-mat.mes-hall] .
  * Vafek and Kang [2020a] O. Vafek and J. Kang, Renormalization group study of hidden symmetry in twisted bilayer graphene with coulomb interactions, Phys. Rev. Lett. 125, 257602 (2020a).
  * Xie and MacDonald [2021] M. Xie and A. H. MacDonald, Weak-field hall resistivity and spin-valley flavor symmetry breaking in magic-angle twisted bilayer graphene, Phys. Rev. Lett. 127, 196401 (2021).
  * Wang  _et al._ [2021] J. Wang, J. Cano, A. J. Millis, Z. Liu, and B. Yang, Exact landau level description of geometry and interaction in a flatband, Phys. Rev. Lett. 127, 246403 (2021).
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Twisted bilayer graphene. iv. exact insulator ground states and phase diagram, Phys. Rev. B 103, 205414 (2021).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Charged skyrmions and topological origin of superconductivity in magic-angle graphene, Science advances 7, eabf5299 (2021).
  * Chen  _et al._ [2021] B.-B. Chen, Y. D. Liao, Z. Chen, O. Vafek, J. Kang, W. Li, and Z. Y. Meng, Realization of topological Mott insulator in a twisted bilayer graphene lattice model, Nature Communications 12, 5480 (2021), arXiv:2011.07602 [cond-mat.str-el] .
  * Bernevig  _et al._ [2021a] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. i. matrix elements, approximations, perturbation theory, and a $k\cdot{}p$ two-band model, Phys. Rev. B 103, 205411 (2021a).
  * Song  _et al._ [2021a] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Twisted bilayer graphene. ii. stable symmetry anomaly, Phys. Rev. B 103, 205412 (2021a).
  * Bernevig  _et al._ [2021b] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Twisted bilayer graphene. iii. interacting hamiltonian and exact symmetries, Phys. Rev. B 103, 205413 (2021b).
  * Bernevig  _et al._ [2021c] B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Twisted bilayer graphene. v. exact analytic many-body excitations in coulomb hamiltonians: Charge gap, goldstone modes, and absence of cooper pairing, Phys. Rev. B 103, 205415 (2021c).
  * Xie  _et al._ [2021a] F. Xie, A. Cowsik, Z.-D. Song, B. Lian, B. A. Bernevig, and N. Regnault, Twisted bilayer graphene. vi. an exact diagonalization study at nonzero integer filling, Phys. Rev. B 103, 205416 (2021a).
  * Sheffer and Stern [2021] Y. Sheffer and A. Stern, Chiral magic-angle twisted bilayer graphene in a magnetic field: Landau level correspondence, exact wave functions, and fractional chern insulators, Phys. Rev. B 104, L121405 (2021).
  * Vafek and Kang [2021] O. Vafek and J. Kang, Lattice model for the coulomb interacting chiral limit of magic-angle twisted bilayer graphene: Symmetries, obstructions, and excitations, Phys. Rev. B 104, 075143 (2021).
  * Da Liao  _et al._ [2021] Y. Da Liao, J. Kang, C. N. Breiø, X. Y. Xu, H.-Q. Wu, B. M. Andersen, R. M. Fernandes, and Z. Y. Meng, Correlation-induced insulating topological phases at charge neutrality in twisted bilayer graphene, Phys. Rev. X 11, 011014 (2021).
  * Liu  _et al._ [2021] S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality, Phys. Rev. Res. 3, 013033 (2021).
  * Kang  _et al._ [2021] J. Kang, B. A. Bernevig, and O. Vafek, Cascades between light and heavy fermions in the normal state of magic-angle twisted bilayer graphene, Phys. Rev. Lett. 127, 266402 (2021).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Kekulé spiral order at all nonzero integer fillings in twisted bilayer graphene, Phys. Rev. X 11, 041063 (2021).
  * Zhang  _et al._ [2022] X. Zhang, G. Pan, B.-B. Chen, H. Li, K. Sun, and Z. Y. Meng, Quantum Monte Carlo sign bounds, topological Mott insulator and thermodynamic transitions in twisted bilayer graphene model, arXiv e-prints , arXiv:2210.11733 (2022), arXiv:2210.11733 [cond-mat.str-el] .
  * Yu  _et al._ [2022a] J. Yu, M. Xie, F. Wu, and S. Das Sarma, Euler Obstructed Cooper Pairing in Twisted Bilayer Graphene: Nematic Nodal Superconductivity and Bounded Superfluid Weight, arXiv e-prints , arXiv:2202.02353 (2022a), arXiv:2202.02353 [cond-mat.supr-con] .
  * Faulstich  _et al._ [2022] F. M. Faulstich, K. D. Stubbs, Q. Zhu, T. Soejima, R. Dilip, H. Zhai, R. Kim, M. P. Zaletel, G. Kin-Lic Chan, and L. Lin, Interacting models for twisted bilayer graphene: a quantum chemistry approach, arXiv e-prints , arXiv:2211.09243 (2022), arXiv:2211.09243 [cond-mat.str-el] .
  * Călugăru  _et al._ [2022] D. Călugăru, N. Regnault, M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, A. Yazdani, O. Vafek, and B. A. Bernevig, Spectroscopy of twisted bilayer graphene correlated insulators, Phys. Rev. Lett. 129, 117602 (2022).
  * Hofmann  _et al._ [2022] J. S. Hofmann, E. Khalaf, A. Vishwanath, E. Berg, and J. Y. Lee, Fermionic monte carlo study of a realistic model of twisted bilayer graphene, Phys. Rev. X 12, 011061 (2022).
  * Kim  _et al._ [2022] S.-W. Kim, S. Jeon, M. Jip Park, and Y. Kim, Replica Higher-Order Topology of Hofstadter Butterflies in Twisted Bilayer Graphene, arXiv e-prints , arXiv:2204.08087 (2022), arXiv:2204.08087 [cond-mat.mes-hall] .
  * Xie  _et al._ [2023] F. Xie, J. Kang, B. A. Bernevig, O. Vafek, and N. Regnault, Phase diagram of twisted bilayer graphene at filling factor $\nu$ =$\pm$3, Phys. Rev. B 107, 075156 (2023), arXiv:2209.14322 [cond-mat.str-el] .
  * Lau and Coleman [2023] L. L. H. Lau and P. Coleman, Topological Mixed Valence Model for Twisted Bilayer Graphene, arXiv e-prints , arXiv:2303.02670 (2023), arXiv:2303.02670 [cond-mat.str-el] .
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Failure of nielsen-ninomiya theorem and fragile topology in two-dimensional systems with space-time inversion symmetry: application to twisted bilayer graphene at magic angle, Physical Review X 9, 021013 (2019).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Faithful tight-binding models and fragile topology of magic-angle bilayer graphene, Physical Review B 99, 195455 (2019).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, All magic angles in twisted bilayer graphene are topological, Physical review letters 123, 036401 (2019).
  * Song  _et al._ [2021b] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Twisted bilayer graphene. ii. stable symmetry anomaly, Physical Review B 103, 205412 (2021b).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Magic-angle twisted bilayer graphene as a topological heavy fermion problem, Physical review letters 129, 047601 (2022).
  * Călugăru  _et al._ [2023] D. Călugăru, M. Borovkov, L. L. Lau, P. Coleman, Z.-D. Song, and B. A. Bernevig, Tbg as topological heavy fermion: Ii. analytical approximations of the model parameters, arXiv preprint arXiv:2303.03429 (2023).
  * Hu  _et al._ [2023] H. Hu, B. A. Bernevig, and A. M. Tsvelik, Kondo Lattice Model of Magic-Angle Twisted-Bilayer Graphene: Hund’s Rule, Local-Moment Fluctuations, and Low-Energy Effective Theory, arXiv e-prints , arXiv:2301.04669 (2023), arXiv:2301.04669 [cond-mat.str-el] .
  * Xie  _et al._ [2019] Y. Xie, B. Lian, B. Jäck, X. Liu, C.-L. Chiu, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene, Nature (London) 572, 101 (2019), arXiv:1906.09274 [cond-mat.mes-hall] .
  * Rozen  _et al._ [2021] A. Rozen, J. M. Park, U. Zondiner, Y. Cao, D. Rodan-Legrain, T. Taniguchi, K. Watanabe, Y. Oreg, A. Stern, E. Berg, _et al._ , Entropic evidence for a pomeranchuk effect in magic-angle graphene, Nature 592, 214 (2021).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, _et al._ , Electronic correlations in twisted bilayer graphene near the magic angle, Nature physics 15, 1174 (2019).
  * Kerelsky  _et al._ [2018] A. Kerelsky, L. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, _et al._ , Magic angle spectroscopy, arXiv preprint arXiv:1812.08776 (2018).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, _et al._ , Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, _et al._ , Cascade of phase transitions and dirac revivals in magic-angle graphene, Nature 582, 203 (2020).
  * Vafek and Kang [2020b] O. Vafek and J. Kang, Renormalization group study of hidden symmetry in twisted bilayer graphene with coulomb interactions, Physical Review Letters 125, 257602 (2020b).
  * Bernevig  _et al._ [2021d] B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Twisted bilayer graphene. v. exact analytic many-body excitations in coulomb hamiltonians: Charge gap, goldstone modes, and absence of cooper pairing, Physical Review B 103, 205415 (2021d).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Strongly correlated chern insulators in magic-angle twisted bilayer graphene, Nature 588, 610 (2020).
  * Wong  _et al._ [2020] D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Cascade of electronic transitions in magic-angle twisted bilayer graphene, Nature (London) 582, 198 (2020), arXiv:1912.06145 [cond-mat.mes-hall] .
  * Lu  _et al._ [2021] X. Lu, B. Lian, G. Chaudhary, B. A. Piot, G. Romagnoli, K. Watanabe, T. Taniguchi, M. Poggio, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, Multiple flat bands and topological Hofstadter butterfly in twisted bilayer graphene close to the second magic angle, Proceedings of the National Academy of Science 118, e2100006118 (2021), arXiv:2006.13963 [cond-mat.mes-hall] .
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, _et al._ , Unconventional sequence of correlated chern insulators in magic-angle twisted bilayer graphene, Nature Physics 17, 1210 (2021).
  * Das  _et al._ [2021] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Symmetry-broken Chern insulators and Rashba-like Landau-level crossings in magic-angle bilayer graphene, Nature Physics 17, 710 (2021), arXiv:2007.13390 [cond-mat.str-el] .
  * Yu  _et al._ [2022b] J. Yu, B. A. Foutty, Z. Han, M. E. Barber, Y. Schattner, K. Watanabe, T. Taniguchi, P. Phillips, Z.-X. Shen, S. A. Kivelson, and B. E. Feldman, Correlated Hofstadter spectrum and flavour phase diagram in magic-angle twisted bilayer graphene, Nature Physics 18, 825 (2022b), arXiv:2108.00009 [cond-mat.str-el] .
  * Saito  _et al._ [2021] Y. Saito, J. Ge, L. Rademaker, K. Watanabe, T. Taniguchi, D. A. Abanin, and A. F. Young, Hofstadter subband ferromagnetism and symmetry-broken chern insulators in twisted bilayer graphene, Nature Physics 17, 478 (2021).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Chern insulators, van hove singularities and topological flat bands in magic-angle twisted bilayer graphene, Nature materials 20, 488 (2021).
  * Xie  _et al._ [2021b] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, _et al._ , Fractional chern insulators in magic-angle twisted bilayer graphene, Nature 600, 439 (2021b).
  * Yu  _et al._ [2022] J. Yu, B. A. Foutty, Z. Han, M. E. Barber, Y. Schattner, K. Watanabe, T. Taniguchi, P. Phillips, Z.-X. Shen, S. A. Kivelson, _et al._ , Correlated hofstadter spectrum and flavour phase diagram in magic-angle twisted bilayer graphene, Nature Physics 18, 825 (2022).
  * Das  _et al._ [2022] I. Das, C. Shen, A. Jaoui, J. Herzog-Arbeitman, A. Chew, C.-W. Cho, K. Watanabe, T. Taniguchi, B. A. Piot, B. A. Bernevig, and D. K. Efetov, Observation of reentrant correlated insulators and interaction-driven fermi-surface reconstructions at one magnetic flux quantum per moiré unit cell in magic-angle twisted bilayer graphene, Phys. Rev. Lett. 128, 217701 (2022).
  * Hejazi  _et al._ [2019] K. Hejazi, C. Liu, and L. Balents, Landau levels in twisted bilayer graphene and semiclassical orbits, Physical Review B 100, 035115 (2019).
  * Bistritzer and MacDonald [2011a] R. Bistritzer and A. H. MacDonald, Moiré butterflies in twisted bilayer graphene, Phys. Rev. B 84, 035440 (2011a).
  * Lian  _et al._ [2018] B. Lian, F. Xie, and B. A. Bernevig, The Landau Level of Fragile Topology, arXiv e-prints , arXiv:1811.11786 (2018), arXiv:1811.11786 [cond-mat.mes-hall] .
  * Wang and Vafek [2022] X. Wang and O. Vafek, Narrow bands in magnetic field and strong-coupling hofstadter spectra, Physical Review B 106, L121111 (2022).
  * Herzog-Arbeitman  _et al._ [2022a] J. Herzog-Arbeitman, A. Chew, and B. A. Bernevig, Magnetic bloch theorem and reentrant flat bands in twisted bilayer graphene at $2\pi$ flux, Phys. Rev. B 106, 085140 (2022a).
  * Herzog-Arbeitman  _et al._ [2022b] J. Herzog-Arbeitman, A. Chew, D. K. Efetov, and B. A. Bernevig, Reentrant correlated insulators in twisted bilayer graphene at 25 t ($2\pi$ flux), Phys. Rev. Lett. 129, 076401 (2022b).
  * Bistritzer and MacDonald [2011b] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011b).
  * Wang and Vafek [2023] X. Wang and O. Vafek, Revisiting bloch electrons in magnetic field: Hofstadter physics via hybrid wannier states, arXiv preprint arXiv:2303.16347 (2023).
  * Luttinger and Kohn [1955] J. M. Luttinger and W. Kohn, Motion of electrons and holes in perturbed periodic fields, Physical Review 97, 869 (1955).
  * Středa [1982] P. Středa, Theory of quantised hall conductivity in two dimensions, Journal of Physics C: Solid State Physics 15 (1982).
  * Babusci  _et al._ [2012] D. Babusci, G. Dattoli, and M. Quattromini, On integrals involving hermite polynomials, Applied Mathematics Letters 25, 1157 (2012).
  * Nuckolls  _et al._ [2023] K. P. Nuckolls, R. L. Lee, M. Oh, D. Wong, T. Soejima, J. P. Hong, D. Călugăru, J. Herzog-Arbeitman, B. A. Bernevig, K. Watanabe, T. Taniguchi, N. Regnault, M. P. Zaletel, and A. Yazdani, Quantum textures of the many-body wavefunctions in magic-angle graphene, arXiv e-prints , arXiv:2303.00024 (2023), arXiv:2303.00024 [cond-mat.mes-hall] .
  * Wang  _et al._ [2022] X. Wang, J. Finney, A. L. Sharpe, L. K. Rodenbach, C. L. Hsueh, K. Watanabe, T. Taniguchi, M. A. Kastner, O. Vafek, and D. Goldhaber-Gordon, Unusual magnetotransport in twisted bilayer graphene from strain-induced open Fermi surfaces, arXiv e-prints , arXiv:2209.08204 (2022), arXiv:2209.08204 [cond-mat.mes-hall] .
  * Herzog-Arbeitman  _et al._ [2020] J. Herzog-Arbeitman, Z.-D. Song, N. Regnault, and B. A. Bernevig, Hofstadter Topology: Noncrystalline Topological Materials at High Flux, Phys. Rev. Lett. 125, 236804 (2020), arXiv:2006.13938 [cond-mat.mes-hall] .
  * Shaffer  _et al._ [2021] D. Shaffer, J. Wang, and L. H. Santos, Theory of Hofstadter superconductors, Phys. Rev. B 104, 184501 (2021), arXiv:2108.04831 [cond-mat.supr-con] .
  * Johnstone  _et al._ [2022] D. Johnstone, M. J. Colbrook, A. E. B. Nielsen, P. Öhberg, and C. W. Duncan, Bulk localized transport states in infinite and finite quasicrystals via magnetic aperiodicity, Phys. Rev. B 106, 045149 (2022), arXiv:2107.05635 [cond-mat.str-el] .
  * Chen and Yao [2022] Y.-F. Chen and D.-X. Yao, Fragile topological phase on the triangular kagome lattice and its bulk-boundary correspondence, arXiv e-prints , arXiv:2202.11736 (2022), arXiv:2202.11736 [cond-mat.mes-hall] .
  * Herzog-Arbeitman  _et al._ [2022] J. Herzog-Arbeitman, Z.-D. Song, L. Elcoro, and B. A. Bernevig, Hofstadter Topology with Real Space Invariants and Reentrant Projective Symmetries, arXiv e-prints , arXiv:2209.10559 (2022), arXiv:2209.10559 [cond-mat.mes-hall] .
  * Fang and Cano [2022] Y. Fang and J. Cano, Symmetry indicators in commensurate magnetic flux, arXiv e-prints , arXiv:2209.10557 (2022), arXiv:2209.10557 [cond-mat.mes-hall] .
  * Guan  _et al._ [2022] Y. Guan, O. V. Yazyev, and A. Kruchkov, Reentrant magic-angle phenomena in twisted bilayer graphene in integer magnetic fluxes, Phys. Rev. B 106, L121115 (2022).
  * Chen  _et al._ [2023] Z. Y. Chen, Z. Zhang, S. A. Yang, and Y. X. Zhao, Classification of time-reversal-invariant crystals with gauge structures, Nature Communications 14, 743 (2023), arXiv:2302.00187 [cond-mat.mes-hall] .
  * Zak [1964a] J. Zak, Magnetic translation group, Phys. Rev. 134, A1602 (1964a).
  * Zak [1964b] J. Zak, Magnetic translation group. ii. irreducible representations, Phys. Rev. 134, A1607 (1964b).

###### Contents

  1. I Introduction
  2. II Hamiltonian and the basis states
  3. III Analysis of the $\frac{1}{q}$ flux sequence
    1. III.1 Change of basis for the $f$-modes
    2. III.2 Analytical solution for flat band, U(4) symmetric, THFM in $\mathbf{B}$
    3. III.3 Integer Fillings $\nu\mathbf{=\pm 1,\pm 2}$
  4. IV Discussion
  5. A Bistritzer-MacDonald model
  6. B Magnetic Translation Operators
    1. B.1 Landau gauge magnetic translation group identities
  7. C Finite $\mathbf{B}$ basis for the $c$ and $f$ fermions
  8. D Computation of Matrix elements
    1. D.1 c-c Coupling
    2. D.2 c-f Coupling
      1. D.2.1 Non-derivative or the zeroth-order coupling $\gamma$
      2. D.2.2 $k_{x}$ coupling
      3. D.2.3 $k_{y}$ coupling
      4. D.2.4 Minimal coupling
      5. D.2.5 Closed form expression for $c$-$f$ coupling at finite field
    3. D.3 Low field Analysis of $c$-$f$ coupling
  9. E Singular Values for $I^{0}$
  10. F Discussion for $\tau=-1$
    1. F.1 CNP
    2. F.2 $\nu=\pm 1$
    3. F.3 $\nu=\pm 2$
  11. G $\nu=\pm 1$, Valley $\mathbf{K}$, Spin $\uparrow$
  12. H $\nu=\pm 2$, Valley $\mathbf{K}$, Spin $\uparrow\downarrow$
  13. I Review of Magnetic Bloch’s Theorem at $2\pi$ Flux
  14. J Rational Flux
    1. J.1 Groundstate Charge Density
    2. J.2 Many-body Charge Excitations

##  Appendix A Bistritzer-MacDonald model

The Bistritzer-Macdonald (BM) model of magic-angle twisted bilayer
graphene[85] in the valley $\mathbf{K}$$(\tau=+1$) is

$\displaystyle
H_{BM}^{\mathbf{K}}(\mathbf{p})=\left(\begin{array}[]{cc}v_{F}\sigma\cdot\mathbf{p}&T(\mathbf{r})e^{i\mathbf{q}_{1}\cdot\mathbf{r}}\\\
e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}T^{\dagger}(\mathbf{r})&v_{F}\sigma\cdot(\mathbf{p}+\hbar\mathbf{q}_{1})\end{array}\right),$
(130)

where $\sigma$ acts in sublattice space and $\mathbf{p}=0$ denotes the moiré
Dirac point $K_{M}$. The interlayer hopping functions are

$\displaystyle
T(\mathbf{r})=\sum_{j=1}^{3}T_{j}e^{-i\mathbf{q}_{1}\cdot\mathbf{r}},$ (131)

where $\mathbf{q}_{1}=k_{\theta}(0,-1)$,
$\mathbf{q}_{2,3}=k_{\theta}(\pm\frac{\sqrt{3}}{2},\frac{1}{2})$,
$k_{\theta}=\frac{8\pi}{3a_{0}}\sin\frac{\theta}{2}=\frac{4\pi}{3L_{m}}$,
$a_{0}\approx 0.246$nm, $L_{m}=\frac{a_{0}}{2\sin{\theta/2}}$ is moiré period.
$\theta=1.05^{\circ}$ in this work.

$\displaystyle
T_{j+1}=w_{0}I_{2}+w_{1}\left(\cos\left(\frac{2\pi}{3}j\right)\sigma_{x}+\sin\left(\frac{2\pi}{3}j\right)\sigma_{y}\right)$
(132)

where $I_{n}$ is the $n\times 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
$\mathbf{L}_{1}=L_{m}(\frac{\sqrt{3}}{2},\frac{1}{2})$ and
$\mathbf{L}_{2}=L_{m}(0,1)$. Thus if $f(\mathbf{r})$ is an eigenstate, then so
is

$\displaystyle\hat{T}_{\mathbf{L}_{1,2}}f(\mathbf{r})=f(\mathbf{r}-\mathbf{L}_{1,2})$
(133)

where $\hat{T}_{\mathbf{L}_{1,2}}$ are the usual discrete translation
operators.

##  Appendix B Magnetic Translation Operators

At $\mathbf{B}\neq 0$ in the Landau gauge $\mathbf{A}=(0,Bx)$, the BM model in
flux is found via minimally substitution:
$H^{\mathbf{K}}_{BM}\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)=H^{\mathbf{K}}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right)$.
The finite $\mathbf{B}$ BM model at $\mathbf{K}^{\prime}$$(\tau=-1$) can be
obtained similarly by applying time reversal to Eq.(130) followed by the
minimal substitution. However at $\mathbf{B}\neq 0$, $H^{\tau}_{BM}$ is still
invariant under the translation by $\mathbf{L}_{2}$, but a translation by
$\mathbf{L}_{1}$ needs to be accompanied by a gauge transformation as

$\displaystyle e^{i\frac{eB}{\hbar
c}L_{1x}y}H^{\tau}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x+\frac{eB}{c}L_{1x}\right)e^{-i\frac{eB}{\hbar
c}L_{1x}y}=H^{\tau}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right).$ (134)

In other words if $\tilde{f}(\mathbf{r})$ is an eigenstate of
$H^{\tau}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right)$, then so is

$\displaystyle\tilde{f}(\mathbf{r})\rightarrow\hat{t}_{\mathbf{L}_{1}}\tilde{f}(\mathbf{r})=e^{i\frac{eB}{\hbar
c}L_{1x}y}\tilde{f}(\mathbf{r}-\mathbf{L}_{1})$ (135)

Thus $\hat{t}_{\mathbf{L}_{1}}=e^{i\frac{eB}{\hbar
c}L_{1x}y}\hat{T}_{\mathbf{L}_{1}}$ is the generator of magnetic translation
by $\mathbf{L}_{1}$. Note that $\hat{t}_{L_{1}}$ can alternatively be
presented as

$\displaystyle\hat{t}_{\mathbf{L}_{1}}=e^{i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}_{\mathbf{L}_{1}}$
(136)

where the magnetic translation wavevector $\mathbf{q}_{\phi}$ is defined as

$\displaystyle\mathbf{q}_{\phi}=\frac{\phi}{\phi_{0}}\left(\frac{1}{2}\mathbf{g}_{1}+\mathbf{g}_{2}\right)$
(137)

where $\mathbf{g}_{1,2}$ are corresponding reciprocal lattice vectors to
$\mathbf{L}_{1,2}$ and are given as:
$\mathbf{g}_{1}=\frac{4\pi}{\sqrt{3}L_{m}}(1,0)$ and
$\mathbf{g}_{2}=\frac{4\pi}{\sqrt{3}L_{m}}(-\frac{1}{2},\frac{\sqrt{3}}{2})$.
Hence

$\hat{t}_{\mathbf{L}_{1}}=e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}},$ (138)

where $L_{2}=\mid\mathbf{L}_{2}\mid=L_{m}$. Magnetic translation by
$\mathbf{L}_{2}$ are generated by
$\hat{t}_{\mathbf{L}_{2}}\tilde{f}(\mathbf{r})=\tilde{f}(\mathbf{r}-\mathbf{L}_{2})$.
We thus have

$\displaystyle\hat{t}_{\mathbf{L}_{2}}\hat{t}_{\mathbf{L}_{1}}\tilde{f}(\mathbf{r})$
$\displaystyle=$ $\displaystyle e^{i\frac{eB}{\hbar
c}L_{1x}(y-L_{2})}\tilde{f}(\mathbf{r}-\mathbf{L}_{1}-\mathbf{L}_{2})$ (139)
$\displaystyle\hat{t}_{\mathbf{L}_{1}}\hat{t}_{\mathbf{L}_{2}}\tilde{f}(\mathbf{r})$
$\displaystyle=$ $\displaystyle e^{i\frac{eB}{\hbar
c}L_{1x}y}\tilde{f}(\mathbf{r}-\mathbf{L}_{1}-\mathbf{L}_{2})$ (140)
$\displaystyle\Rightarrow\hat{t}_{\mathbf{L}_{2}}\hat{t}_{\mathbf{L}_{1}}$
$\displaystyle=$ $\displaystyle e^{-2\pi
i\frac{\phi}{\phi_{0}}}\hat{t}_{\mathbf{L}_{1}}\hat{t}_{\mathbf{L}_{2}}$ (141)

where $\phi_{0}=\frac{hc}{e}$ and the flux through the unit cell is
$\phi=BL_{1x}L_{2}$. If $\phi/\phi_{0}=p/q$, with $p$ and $q$ relatively prime
integers,

$\displaystyle\left[\hat{t}^{q}_{\mathbf{L}_{2}},\hat{t}_{\mathbf{L}_{1}}\right]=0$
(142)

and

$\displaystyle\left[\hat{t}_{\mathbf{L}_{1,2}},H^{\tau}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right)\right]=0.$
(143)

The simultaneous eigenstates of the magnetic translation operators
$\hat{t}_{\mathbf{L}_{1}}$ and $\hat{t}^{q}_{\mathbf{L}_{2}}$ can thus be used
to produce a complete and orthonormal set of basis states for solving the
$\mathbf{B}\neq 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 for
$\hat{t}^{s}_{\mathbf{L}_{1}}$ which relates it to
$\hat{T}^{s}_{\mathbf{L}_{1}}$.

###  B.1 Landau gauge magnetic translation group identities

Using Eq.(138), we have

$\displaystyle\hat{t}^{s}_{\mathbf{L}_{1}}$ $\displaystyle=$
$\displaystyle\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)^{s}=\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\ldots\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)$ (144)
$\displaystyle=$ $\displaystyle\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\ldots\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{2y-L_{1y}}{L_{2}}}\hat{T}^{2}_{\mathbf{L}_{1}}\right)$
(145) $\displaystyle=$ $\displaystyle\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\ldots\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{3y-2L_{1y}-L_{1y}}{L_{2}}}\hat{T}^{3}_{\mathbf{L}_{1}}\right)$
(146) $\displaystyle=$ $\displaystyle\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\ldots\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\right)\left(e^{2\pi
i\frac{\phi}{\phi_{0}}\frac{4y-3L_{1y}-2L_{1y}-L_{1y}}{L_{2}}}\hat{T}^{4}_{\mathbf{L}_{1}}\right).$
(147)

Now, $1+2+3+\ldots+(s-1)=\frac{1}{2}s(s-1)$. Therefore,

$\displaystyle\hat{t}^{s}_{\mathbf{L}_{1}}$ $\displaystyle=$ $\displaystyle
e^{-\pi is(s-1)\frac{p}{q}\frac{L_{1y}}{L_{2}}}e^{2\pi
is\frac{p}{q}\frac{y}{L_{2}}}\hat{T}^{s}_{\mathbf{L}_{1}}$ (148)

where $\phi/\phi_{0}=p/q$. Although the above formula was derived for positive
$s$, it also holds for negative $s$. This can be seen by noting that

$\displaystyle\left(e^{-\pi
is(s-1)\frac{\phi}{\phi_{0}}\frac{L_{1y}}{L_{2}}}e^{i\frac{eB}{\hbar
c}sL_{1x}y}\hat{T}^{s}_{\mathbf{L}_{1}}\right)\left(e^{\pi
is(-s-1)\frac{\phi}{\phi_{0}}\frac{L_{1y}}{L_{2}}}e^{-i\frac{eB}{\hbar
c}sL_{1x}y}\hat{T}^{-s}_{\mathbf{L}_{1}}\right)=1.$ (149)

So, $\hat{t}^{-1}_{\mathbf{L}_{1}}$ denotes the inverse of
$\hat{t}_{\mathbf{L}_{1}}$, and
$\hat{t}^{-s}_{\mathbf{L}_{1}}=\left(\hat{t}^{-1}_{\mathbf{L}_{1}}\right)^{s}$.

##  Appendix C Finite $\mathbf{B}$ basis for the $c$ and $f$ fermions

The $\mathbf{B}=0$ basis for $f$ fermions is constituted by two AA-stacking
localised Wannier states per valley per spin [56]. We denote them as
$W_{\mathbf{0},b\tau}(\mathbf{r})$, for $b\in\\{1,2\\}$ and $\tau=\pm 1$.
Using the recently introduced method [86], we can construct the finite
$\mathbf{B}$ basis by first building hybrid Wannier states out of the zero
field Wannier states followed by projection onto a representation of magnetic
translation group(MTG). For given spin, the hybrid Wannier states read

$\displaystyle w_{b\tau}(\mathbf{r},k_{2}\mathbf{g}_{2})$ $\displaystyle=$
$\displaystyle\sum_{n\in\mathbb{Z}}e^{ik_{2}\mathbf{g}_{2}\cdot(n\mathbf{L}_{2})}\hat{t}^{n}_{\mathbf{L}_{2}}W_{\mathbf{0},b\tau}(\mathbf{r})=\sum_{n\in\mathbb{Z}}e^{2\pi
ik_{2}n}W_{0,b\tau}(\mathbf{r}-n\mathbf{L}_{2}).$ (150)

where $k_{2}\in[0,1)$. These hybrid Wannier states although localized along
the $x$ direction are Bloch-like extended in the $y$ direction and thus
eigenstates of $\hat{t}_{\mathbf{L}_{2}}$

$\hat{t}_{\mathbf{L}_{2}}w_{b\tau}(\mathbf{r},k_{2}\mathbf{g}_{2})=w_{b\tau}(\mathbf{r}-\mathbf{L}_{2},k_{2}\mathbf{g}_{2})=e^{-2\pi
ik_{2}}w_{b\tau}(\mathbf{r},k_{2}\mathbf{g}_{2}).$ (151)

We next construct normalized eigenstates of MTG out of the hybrid Wannier
states as

$\displaystyle\eta_{b\tau k_{1}k_{2}}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\mathcal{N}}}\sum_{s\in\mathbb{Z}}e^{2\pi
isk_{1}}\left(\hat{t}^{s}_{\mathbf{L}_{1}}w_{b\tau}(\mathbf{r},k_{2}\mathbf{g}_{2})\right),$
(152)

where $k_{1}\in[0,1)$ and normalization $\mathcal{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.(156)-(158). The domains of $k_{1,2}$ can be
understood via noting that $\eta_{b\tau k_{1}k_{2}}$ is periodic under
$k_{1,2}\rightarrow k_{1,2}+1$. Under magnetic translations
$\hat{t}_{\mathbf{L}_{1}}$ and $\hat{t}^{q}_{\mathbf{L}_{2}}$, we then have

$\displaystyle\hat{t}_{\mathbf{L}_{1}}\eta_{b\tau k_{1}k_{2}}$
$\displaystyle=$ $\displaystyle e^{-2\pi ik_{1}}\eta_{b\tau k_{1}k_{2}}$ (153)
$\displaystyle\hat{t}^{q}_{\mathbf{L}_{2}}\eta_{b\tau k_{1}k_{2}}$
$\displaystyle=$ $\displaystyle e^{-2\pi iqk_{2}}\eta_{b\tau k_{1}k_{2}}.$
(154)

Thus states with $k_{2}$ that differ by $1/q$ have the same eigenvalue under
$\hat{t}^{q}_{\mathbf{L}_{2}}$. 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

$\displaystyle\eta_{b\tau kr^{\prime}}(\mathbf{r})=\eta_{b\tau
k_{1}k_{2}+\frac{r^{\prime}}{q}}(\mathbf{r})\;\;\text{,}\;\;k=(k_{1},k_{2})\in[0,1)\otimes[0,1/q)\;\;\text{,}\;\;r^{\prime}=0,1\ldots
q-1.$ (155)

Thus $\eta_{b\tau kr^{\prime}}(\mathbf{r})$ with different $k$ quantum numbers
are guaranteed to be orthogonal. Moreover we have $q$ states for given $k$,
corresponding to $r^{\prime}\in\\{0\ldots q-1\\}$ for each $b\in\\{1,2\\}$,
i.e. a total of $2q$ heavy fermion states for given $k$. For same $k$ and
different $r^{\prime}$ the overlaps for given valley are

$\displaystyle\int d^{2}\mathbf{r}\eta^{*}_{b\tau
kr^{\prime}_{1}}(\mathbf{r})\eta_{b^{\prime}\tau
kr^{\prime}_{2}}(\mathbf{r})=$
$\displaystyle\frac{1}{\mathcal{N}}\sum_{s_{1}s_{2}\in\mathbb{Z}}\sum_{n_{1}n_{2}\in\mathbb{Z}}e^{2\pi
i(s_{2}-s_{1})k_{1}}e^{-2\pi i(k_{2}+\frac{r^{\prime}_{1}}{q})n_{1}}e^{2\pi
i(k_{2}+\frac{r^{\prime}_{2}}{q})n_{2}}e^{\pi
is_{1}(s_{1}-1)\frac{p}{q}\frac{L_{1y}}{L_{2}}}e^{-\pi
is_{2}(s_{2}-1)\frac{p}{q}\frac{L_{1y}}{L_{2}}}$ $\displaystyle\times\int
d^{2}\mathbf{r}e^{2\pi
i(s_{2}-s_{1})\frac{p}{q}\frac{y}{L_{2}}}W^{*}_{\mathbf{0}b\tau}(\mathbf{r}-s_{1}\mathbf{L}_{1}-n_{1}\mathbf{L}_{2})W_{\mathbf{0}b^{\prime}\tau}(\mathbf{r}-s_{2}\mathbf{L}_{1}-n_{2}\mathbf{L}_{2}).$
(156)

Now since the $\mathbf{B}=0$ Wannier states are well localised and
orthonormal, we have

$\displaystyle\int d^{2}\mathbf{r}e^{2\pi
i(s_{2}-s_{1})\frac{p}{q}\frac{y}{L_{2}}}W^{*}_{\mathbf{0}b\tau}(\mathbf{r}-s_{1}\mathbf{L}_{1}-n_{1}\mathbf{L}_{2})W_{\mathbf{0}b^{\prime}\tau}(\mathbf{r}-s_{2}\mathbf{L}_{1}-n_{2}\mathbf{L}_{2})\approx\delta_{bb^{\prime}}\delta_{s_{1}s_{2}}\delta_{n_{1}n_{2}}$
(157)

and thus,

$\displaystyle\int d^{2}\mathbf{r}\eta^{*}_{b\tau
kr^{\prime}_{1}}(\mathbf{r})\eta_{b^{\prime}\tau
kr^{\prime}_{2}}(\mathbf{r})\approx\delta_{bb^{\prime}}\frac{1}{\mathcal{N}}\sum_{s_{1}\in\mathbb{Z}}\sum_{n_{1}\in\mathbb{Z}}e^{2\pi
i\frac{r^{\prime}_{2}-r^{\prime}_{1}}{q}n_{1}}=\delta_{r^{\prime}_{1}r^{\prime}_{2}}\delta_{bb^{\prime}}\frac{s_{tot}n_{tot}}{\mathcal{N}}=\delta_{r^{\prime}_{1}r^{\prime}_{2}}\delta_{bb^{\prime}}.$
(158)

Henceforth we have a complete orthonormal basis for the heavy fermions in
finite $\mathbf{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
$\mathbf{B}=0$, the basis for $c$ fermions is constituted by four
$\mathbf{k}\cdot\mathbf{p}$ Bloch states at the $\Gamma$ point in the moiré
Brillouin Zone(mBZ), per valley per spin. We denote them by
$\tilde{\Psi}_{\Gamma a\tau}$ for $a\in\\{1\ldots 4\\}$ and $\tau=\pm 1$.

We promote the $\mathbf{B}=0$ $c$ fermion basis to finite $\mathbf{B}$ using a
result obtained when the $\mathbf{k}\cdot\mathbf{p}$ method is extended to
finite $\mathbf{B}$[87], which prescribes the ansatz for finite field
$\mathbf{k}\cdot\mathbf{p}$ states to be Landau level(LL) coefficients on top
of the zero field $\mathbf{k}\cdot\mathbf{p}$ states. Henceforth, we define
the finite $\mathbf{B}$ basis for $\mathbf{k}\cdot\mathbf{p}$ Bloch states as
an expansion over product of LLs and the $\mathbf{B}=0$
$\mathbf{k}\cdot\mathbf{p}$ Bloch states at the $\Gamma$ point in mBZ. However
in order to use the same quantum numbers label as finite $\mathbf{B}$ $f$
fermion basis we first project the Landau gauge LL onto the representation of
MTG as

$\displaystyle\chi_{k_{1}k_{2}m}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\ell
L_{2}}}\frac{1}{\sqrt{\mathcal{N}}}\sum_{s\in\mathbb{Z}}e^{2\pi
isk_{1}}\left(\hat{t}^{s}_{\mathbf{L}_{1}}e^{2\pi
ik_{2}\frac{y}{L_{2}}}\varphi_{m}\left(x-k_{2}\frac{2\pi\ell^{2}}{L_{2}}\right)\right)$
(159)

where $k_{1}\in[0,1)$ and unlike before, $k_{2}\in[0,p/q)$ as justified in
Eqs.(162)-(165). The functions $\varphi_{m}$ are $m$th 1D harmonic oscillator
states that determine the Landau gauge Landau levels

$\displaystyle\varphi_{m}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi^{\frac{1}{4}}}\frac{1}{\sqrt{2^{m}m!}}e^{-x^{2}/2\ell^{2}}H_{m}(x/\ell)$
(160)

where $H_{m}$ is the Hermite polynomial and $\ell^{2}=\hbar c/(eB)$. The
normalization of the LL MTG states is explained in Eqs.(172)-(177). The
domains for the quantum numbers $k_{1,2}$ can be understood by noting that
$\chi_{k_{1}k_{2}m}$ is periodic under $k_{1}\rightarrow k_{1}+1$ and up to a
phase under $k_{2}\rightarrow k_{2}+p/q$ as shown below

$\displaystyle\chi_{k_{1}k_{2}+\frac{p}{q}m}=\frac{1}{\sqrt{\ell
L_{2}}}\frac{1}{\sqrt{\mathcal{N}}}\sum_{s\in\mathbb{Z}}e^{2\pi
isk_{1}}\left(\hat{t}^{s}_{\mathbf{L}_{1}}e^{2\pi
i(k_{2}+\frac{p}{q})\frac{y}{L_{2}}}\varphi_{m}\left(x-(k_{2}+\frac{p}{q})\frac{2\pi\ell^{2}}{L_{2}}\right)\right),$
(161)

now using the fact, $\frac{2\pi\ell^{2}}{L_{2}}\frac{p}{q}=L_{1x}$ we have

$\displaystyle\chi_{k_{1}k_{2}+\frac{p}{q}m}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{\ell
L_{2}}}\frac{1}{\sqrt{\mathcal{N}}}\sum_{s\in\mathbb{Z}}e^{2\pi
isk_{1}}\left(\hat{t}^{s}_{\mathbf{L}_{1}}e^{2\pi
i(k_{2}+\frac{p}{q})\frac{y}{L_{2}}}\varphi_{m}\left(x-L_{1x}-k_{2}\frac{2\pi\ell^{2}}{L_{2}}\right)\right)$
(162) $\displaystyle=$ $\displaystyle e^{2\pi
ik_{2}\frac{L_{1y}}{L_{2}}}\frac{1}{\sqrt{\ell
L_{2}}}\frac{1}{\sqrt{\mathcal{N}}}\sum_{s\in\mathbb{Z}}e^{2\pi
isk_{1}}\left(\hat{t}^{s}_{\mathbf{L}_{1}}e^{2\pi
i\frac{p}{q}\frac{y}{L_{2}}}\hat{T}_{\mathbf{L}_{1}}\left(e^{2\pi
ik_{2}\frac{y}{L_{2}}}\varphi_{m}\left(x-k_{2}\frac{2\pi\ell^{2}}{L_{2}}\right)\right)\right)$
(163) $\displaystyle=$ $\displaystyle e^{2\pi
ik_{2}\frac{L_{1y}}{L_{2}}}\frac{1}{\sqrt{\ell
L_{2}}}\frac{1}{\sqrt{\mathcal{N}}}\sum_{s\in\mathbb{Z}}e^{2\pi
isk_{1}}\left(\hat{t}^{s+1}_{\mathbf{L}_{1}}e^{2\pi
ik_{2}\frac{y}{L_{2}}}\varphi_{m}\left(x-k_{2}\frac{2\pi\ell^{2}}{L_{2}}\right)\right)$
(164) $\displaystyle=$ $\displaystyle e^{-2\pi
i\left(k_{1}-k_{2}\frac{L_{1y}}{L_{2}}\right)}\chi_{k_{1}k_{2}m}$ (165)

The domain for $k_{2}$ can also be understood in the context of the fact that
$\chi_{k_{1}k_{2}m}$ are essentially LLs and thus behave as Chern number $+1$
states in $\mathbf{B}\neq 0$. In order to illustrate it, let us discuss the
total number of available states for $\chi_{k_{1}k_{2}m}$. We start by noting
that under magnetic translations $\hat{t}_{\mathbf{L}_{1}}$ and
$\hat{t}^{q}_{\mathbf{L}_{2}}$

$\displaystyle\hat{t}_{\mathbf{L}_{1}}\chi_{k_{1}k_{2}m}$ $\displaystyle=$
$\displaystyle e^{-2\pi ik_{1}}\chi_{k_{1}k_{2}m},$ (166)
$\displaystyle\hat{t}^{q}_{\mathbf{L}_{2}}\chi_{k_{1}k_{2}m}$ $\displaystyle=$
$\displaystyle e^{-2\pi iqk_{2}}\chi_{k_{1}k_{2}m}.$ (167)

If we have a system size $N_{1}\mathbf{L}_{1}$ and $qN_{2}\mathbf{L}_{2}$,
then

$\displaystyle\hat{t}^{N_{1}}_{\mathbf{L}_{1}}\chi_{k_{1}k_{2}m}=\chi_{k_{1}k_{2}m}$
(168) $\displaystyle\Rightarrow$ $\displaystyle e^{-2\pi
ik_{1}N_{1}}=1;\;\Rightarrow\;\;k_{1}=0,\frac{1}{N_{1}},\frac{2}{N_{1}},\ldots,1-\frac{1}{N_{1}}$
(169)

Similarly,

$\displaystyle\hat{t}^{qN_{2}}_{\mathbf{L}_{2}}\chi_{k_{1}k_{2}m}=\chi_{k_{1}k_{2}m}$
(170) $\displaystyle\Rightarrow$ $\displaystyle e^{-2\pi
ik_{2}qN_{2}}=1;\;\Rightarrow\;\;k_{2}=0,\frac{1}{qN_{2}},\frac{2}{qN_{2}},\ldots,\frac{p}{q}-\frac{1}{qN_{2}}$
(171)

So, the total number of states is
$N_{1}\frac{p}{q}qN_{2}=\frac{\Phi}{\Phi_{0}}N_{1}qN_{2}=\frac{\Phi}{\Phi_{0}}N_{1}qN_{2}\frac{L_{1x}L_{2}}{A_{uc}}=\frac{\Phi}{\Phi_{0}}\frac{A_{tot}}{A_{uc}}=\frac{B}{\Phi_{0}}A_{tot}=\frac{\Phi_{tot}}{\Phi_{0}}$.
Here $A_{uc}=\hat{z}\cdot(\mathbf{L}_{1}\times\mathbf{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:

$\displaystyle\int
d^{2}\mathbf{r}\chi^{*}_{k_{1}k_{2}m}(\mathbf{r})\chi_{k_{1}^{\prime}k_{2}^{\prime}m^{\prime}}(\mathbf{r})=$
(172) $\displaystyle\frac{1}{\ell
L_{2}\mathcal{N}}\sum_{s\in\mathbb{Z}}\sum_{s^{\prime}\in\mathbb{Z}}e^{-2\pi
is\left(k_{1}-\left(k_{2}-\frac{p}{q}\right)\frac{L_{1y}}{L_{2}}\right)}e^{i\pi
s(s+1)\frac{p}{q}\frac{L_{1y}}{L_{2}}}e^{2\pi
is^{\prime}\left(k^{\prime}_{1}-\left(k^{\prime}_{2}-\frac{p}{q}\right)\frac{L_{1y}}{L_{2}}\right)}e^{-i\pi
s^{\prime}(s^{\prime}+1)\frac{p}{q}\frac{L_{1y}}{L_{2}}}$ $\displaystyle\int
dye^{-iy\left(s\frac{p}{q}+k_{2}\right)\frac{2\pi}{L_{2}}}e^{iy\left(s^{\prime}\frac{p}{q}+k^{\prime}_{2}\right)\frac{2\pi}{L_{2}}}\int
dx\varphi_{m}\left(x-\left(s\frac{p}{q}+k_{2}\right)\frac{2\pi}{L_{2}}\ell^{2}\right)\varphi_{m^{\prime}}\left(x-\left(s^{\prime}\frac{p}{q}+k^{\prime}_{2}\right)\frac{2\pi}{L_{2}}\ell^{2}\right)$
(173) $\displaystyle=\frac{1}{\ell
L_{2}\mathcal{N}}\sum_{s\in\mathbb{Z}}\sum_{s^{\prime}\in\mathbb{Z}}e^{-2\pi
is\left(k_{1}-\left(k_{2}-\frac{p}{q}\right)\frac{L_{1y}}{L_{2}}\right)}e^{i\pi
s(s+1)\frac{p}{q}\frac{L_{1y}}{L_{2}}}e^{2\pi
is^{\prime}\left(k^{\prime}_{1}-\left(k^{\prime}_{2}-\frac{p}{q}\right)\frac{L_{1y}}{L_{2}}\right)}e^{-i\pi
s^{\prime}(s^{\prime}+1)\frac{p}{q}\frac{L_{1y}}{L_{2}}}$ $\displaystyle
n_{tot}L_{2}\delta_{s,s^{\prime}}\delta_{k_{2},k_{2}^{\prime}}\ell\delta_{m,m^{\prime}}$
(174) $\displaystyle=\frac{1}{L_{2}\mathcal{N}}\sum_{s\in\mathbb{Z}}e^{2\pi
is(k_{1}-k_{1}^{\prime})}n_{tot}L_{2}\delta_{k_{2},k_{2}^{\prime}}\delta_{m,m^{\prime}}$
(175)
$\displaystyle=\frac{1}{L_{2}\mathcal{N}}s_{tot}\delta_{k_{1},k_{1}^{\prime}}n_{tot}L_{2}\delta_{k_{2},k_{2}^{\prime}}\delta_{m,m^{\prime}}$
(176)
$\displaystyle=\delta_{k_{1},k_{1}^{\prime}}\delta_{k_{2},k_{2}^{\prime}}\delta_{m,m^{\prime}}.$
(177)

where in Eq.(173)-(174), we have used the fact that $y$-integral implies
$k_{2}-k^{\prime}_{2}=(s^{\prime}-s)\frac{p}{q}$. Now since
$k_{2},k^{\prime}_{2}\in[0,\frac{p}{q})$, the integral evaluates to
$n_{tot}L_{2}\delta_{k_{2},k_{2}^{\prime}}\delta_{s,s^{\prime}}$. 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 $\ell\delta_{mm^{\prime}}$. Eventually in Eq.(175)-(176), the $s$ summation
implies $k_{1}-k^{\prime}_{1}=\text{integer}$. But since
$k_{1},k^{\prime}_{1}\in[0,1)$, we must have $k_{1}=k^{\prime}_{1}$ and the
summation evaluates to $s_{tot}\delta_{k_{1},k^{\prime}_{1}}$. Therefore
$\chi$’s are orthogonal and normalized. To have the same quantum number label
as that of $\mathbf{B}\neq 0$ $f$ fermion states, we relabel the MTG LLs as

$\displaystyle\chi_{krm}(\mathbf{r})=\chi_{k_{1}k_{2}+\frac{r}{q}m}(\mathbf{r})\;\;\text{,}\;\;k=(k_{1},k_{2})\in[0,1)\otimes[0,1/q)\;\;\text{,}\;\;r=0,1\ldots
p-1.$ (178)

Thus for given $k$, we have $p$ states, corresponding to $r\in\\{0\ldots
p-1\\}$, i.e. $\frac{p}{q}$ states per moiré unit cell per valley per spin as
expected. Finally we have the finite $\mathbf{B}$ c fermion basis states as
$\Psi_{a\tau}\chi_{krm}$, where

$\Psi_{a\tau}=\sqrt{NA_{uc}}\tilde{\Psi}_{\Gamma a\tau}$ (179)

with $N$ being the total number of moiré unit cells. The $\sqrt{NA_{uc}}$
factor is required due to normalization of $\chi$. Note that since
$\Psi_{a\tau}$ is a Bloch state at $\Gamma$,

$\displaystyle\hat{T}_{\mathbf{L}_{i}}\Psi_{a\tau}(\mathbf{r})=\Psi_{a\tau}(\mathbf{r})$
(180)
$\displaystyle\implies\hat{t}_{\mathbf{L}_{i}}\left(\Psi_{a\tau}(\mathbf{r})\chi_{krm}\right)(\mathbf{r})=\Psi_{a\tau}\hat{t}_{\mathbf{L}_{i}}\chi_{krm}(\mathbf{r}),$
(181)

where $i\in\\{1,2\\}$. Thus $\Psi_{a\tau}\chi_{krm}$ are MTG eigenstates with
same eigenvalue as that of $\chi_{krm}$. Now let us discuss the
orthonormalization for these states. Since in this work we focus on
$\frac{1}{q}$ sequence, for the following discussion on orthonormality we set
index $r$ in $\chi_{krm}$ to zero and work with $\chi_{k0m}$. The Bloch
periodicity for $\Psi_{\Gamma a\tau}(\mathbf{r})$ allows us to expand it as

$\displaystyle\Psi_{\Gamma
a\tau}(\mathbf{r})=\frac{1}{\sqrt{NA_{uc}}}\sum_{\mathbf{G}}e^{i\mathbf{G}\cdot\mathbf{r}}C_{a\tau\mathbf{G}},$
(182)

where $\mathbf{G}=n_{1}\mathbf{g}_{1}+n_{2}\mathbf{g}_{2}$,
$n_{1,2}\in\mathbb{Z}$, denote moiré reciprocal lattice vectors and

$\displaystyle C_{a\tau\mathbf{G}}=\frac{1}{\sqrt{NA_{uc}}}\int
d^{2}\mathbf{r}e^{-i\mathbf{G}\cdot\mathbf{r}}\tilde{\Psi}_{\Gamma
a\tau}(\mathbf{r})$ (183)

can be found in [57]. Now using Eq.(179), the overlap matrix can be given as

$\displaystyle O^{(q,m_{max})}_{[m^{\prime},a^{\prime}][m,a]}=\int
d^{2}\mathbf{r}\left(\Psi_{a^{\prime}\tau}\chi_{k^{\prime}0m^{\prime}}\right)^{\ast}\Psi_{a\tau}\chi_{k0m}=\sum_{\mathbf{G}_{1},\mathbf{G}_{2}}C^{\ast}_{a^{\prime}\tau\mathbf{G}_{1}}C_{a\tau\mathbf{G}_{2}}\int
d^{2}\mathbf{r}e^{i(\mathbf{G}_{2}-\mathbf{G}_{1})\cdot\mathbf{r}}\chi^{\ast}_{k^{\prime}0m^{\prime}}\chi_{k0m},$
(184)

where $1/q$ is the flux per unit cell per flux quantum $\Phi_{0}$ and
$m_{max}$ denotes the upper cutoff on the indices $m$ and $m^{\prime}$. The
integral can be computed as:

$\displaystyle\int
d^{2}\mathbf{r}e^{i\mathbf{g}\cdot\mathbf{r}}\chi^{*}_{k^{\prime}0m^{\prime}}(\mathbf{r})\chi_{k0m}(\mathbf{r})$
$\displaystyle=$ $\displaystyle\frac{1}{\ell
L_{2}\mathcal{N}}\sum_{s,s^{\prime}\in\mathbb{Z}}e^{2\pi
i(k_{1}s-k^{\prime}_{1}s^{\prime})}\int d^{2}\mathbf{r}(e^{2\pi
ik_{2}^{\prime}\frac{y}{L_{2}}}\varphi_{m^{\prime}}(x-k^{\prime}_{2}\frac{2\pi\ell^{2}}{L_{2}}))^{*}t^{-s^{\prime}}_{L_{1}}e^{i\mathbf{g}\cdot\mathbf{r}}t^{s}_{L_{1}}(e^{2\pi
ik_{2}\frac{y}{L_{2}}}\varphi_{m}(x-k_{2}\frac{2\pi\ell^{2}}{L_{2}}))$
$\displaystyle=$
$\displaystyle\delta_{k_{1},k^{\prime}_{1}}\frac{s_{tot}}{\ell
L_{2}\mathcal{N}}\sum_{s\in\mathbb{Z}}e^{2\pi isk_{1}}\int
d^{2}\mathbf{r}(e^{2\pi
ik_{2}^{\prime}\frac{y}{L_{2}}}\varphi_{m^{\prime}}(x-k^{\prime}_{2}\frac{2\pi\ell^{2}}{L_{2}}))^{*}e^{i\mathbf{g}\cdot\mathbf{r}}t^{s}_{L_{1}}(e^{2\pi
ik_{2}\frac{y}{L_{2}}}\varphi_{m}(x-k_{2}\frac{2\pi\ell^{2}}{L_{2}}))$
$\displaystyle=$
$\displaystyle\delta_{k_{1},k^{\prime}_{1}}\sum_{s\in\mathbb{Z}}\delta_{k_{2}+\frac{s}{q},k_{2}^{\prime}-\frac{L_{2}g_{y}}{2\pi}}e^{2\pi
is(k_{1}-\frac{k_{2}}{2})}e^{-i\pi s(s-1)\frac{1}{2q}}\frac{1}{\ell}\int
dx\varphi_{m^{\prime}}(x-k^{\prime}_{2}\frac{2\pi\ell^{2}}{L_{2}})e^{ig_{x}x}\varphi_{m}(x-(k_{2}+\frac{s}{q})\frac{2\pi\ell^{2}}{L_{2}})$

The delta function constraint tells us

$\displaystyle\delta_{k_{2}+\frac{s}{q},k_{2}^{\prime}-\frac{L_{2}g_{y}}{2\pi}}=\delta_{k_{2},k^{\prime}_{2}}\delta_{s,-qg_{2}},$
(186)

where $g_{2}=g_{y}\frac{L_{2}}{2\pi}\in\mathbb{Z}$. Note that getting
$\delta_{k_{1(2)},k_{1^{\prime}(2^{\prime})}}$ from the overlap is a direct
consequence of Eqs.(166)-(167) for $\frac{1}{q}$ sequence. The above integral
can be evaluated as

Figure 10: Figures showing the diagonal values of Overlap matrix
$O^{(q,m_{max})}_{[m^{\prime},a^{\prime}][m,a]}$ and frobenius norm of
$O^{(q,m_{max})}-\mathbb{I}_{m_{max}}$ $\forall$ $m_{max}=0\ldots q-1$ for
(a,b) $q=30$, (c,d) $q=50$, (e,f) $q=60$ and (g,h) $q=70$ respectively. For
higher values of $q$, the value of $m_{max}$ for which diagonal values start
deviating from 1 have been (i) plotted using $\bullet$.

$\displaystyle=$ $\displaystyle\frac{1}{\ell}\int
dx\varphi_{m^{\prime}}(x-k_{2}\frac{2\pi\ell^{2}}{L_{2}})e^{ig_{x}x}\varphi_{m}(x-(k_{2}-g_{2})\frac{2\pi\ell^{2}}{L_{2}})$
$\displaystyle=$
$\displaystyle\frac{1}{\ell}e^{ig_{x}k_{2}\frac{2\pi\ell^{2}}{L_{2}}}\int
dx\varphi_{m^{\prime}}(x)e^{ig_{x}x}\varphi_{m}(x+g_{2}\frac{2\pi\ell^{2}}{L_{2}})$
$\displaystyle=$
$\displaystyle\frac{1}{\ell}e^{ig_{x}k_{2}\frac{2\pi\ell^{2}}{L_{2}}}\int
dx\varphi_{m^{\prime}}(x)e^{ig_{x}x}e^{ig_{y}\ell^{2}p_{x}}\varphi_{m}(x)$
$\displaystyle=$
$\displaystyle\frac{1}{\ell}e^{ig_{x}k_{2}\frac{2\pi\ell^{2}}{L_{2}}}e^{-\frac{i}{2}g_{x}g_{y}\ell^{2}}\int
dx\varphi_{m^{\prime}}(x)e^{c_{+}a^{\dagger}+c_{-}a}\varphi_{m}(x)$
$\displaystyle=$ $\displaystyle
e^{ig_{x}k_{2}\frac{2\pi\ell^{2}}{L_{2}}}e^{-\frac{i}{2}g_{x}g_{y}\ell^{2}}\bra{m^{\prime}}e^{c_{+}a^{\dagger}+c_{-}a}\ket{m},$

where $p_{x}$ is momentum operator along $x$ direction later changed to
harmonic oscillator(h.o.) basis: $x=\frac{\ell}{\sqrt{2}}(a+a^{\dagger})$, and
$p_{x}=\frac{i}{\sqrt{2\ell}}(a^{\dagger}-a)$. The h.o. kets are defined as
$\langle m\mid x\rangle=\frac{1}{\sqrt{\ell}}\varphi_{m}(x)$, such that
$\langle m\mid n\rangle=\delta_{m,n}$. Eventually we use the BCH formula,
$e^{X}e^{Y}=e^{X+Y+\frac{1}{2}[X,Y]+\ldots}$ , and

$\displaystyle c_{\pm}=\frac{i\ell}{\sqrt{2}}(g_{x}\pm ig_{y}).$ (188)

We also know the identity

$\displaystyle\bra{m^{\prime}}e^{c_{+}a^{\dagger}+c_{-}a}\ket{m}=\begin{cases}e^{\frac{1}{2}c_{+}c_{-}}\sqrt{\frac{m!}{n!}}(c_{+})^{n-m}L^{n-m}_{m}(-c_{+}c_{-})&\text{for}\leavevmode\nobreak\
\leavevmode\nobreak\ n\geq m,\\\
e^{\frac{1}{2}c_{+}c_{-}}\sqrt{\frac{n!}{m!}}(c_{-})^{m-n}L^{m-n}_{n}(-c_{+}c_{-})&\text{for}\leavevmode\nobreak\
\leavevmode\nobreak\ n<m,\end{cases}$ (189)

where $L^{m-n}_{n}(x)$ is the associated Laguerre polynomial,

$\displaystyle
L^{m}_{N}(x)=\sum_{k=0}^{N}\frac{(N+m)!}{(N-k)!(m+k)!k!}(-x)^{k}.$ (190)

As shown in Eq.(189), the overlaps are exponentially small in $(c_{+}c_{-})$,
i.e. $(g\ell)^{2}$, and thus negligible. However, beyond an upper-bound on
$m$, the orthonormality begins to fail. Fig.(10-a)-(10-h) show that the
diagonal values of overlap matrix in Eq.(184) start deviating from $1$ beyond
some upper-bound on $m$ exactly when the Frobenius norm of $O-\mathbb{I}$,
with $\mathbb{I}$ being identity matrix of dimension same as that of $O$,
starts deviating from $0$. Thus the deviation of diagonal values from $1$ is a
good indicative of the breaking up of orthonormality. In this work we choose
$m_{max}=\lceil\frac{q-3}{2}\rceil$, which is well below the upper-bound as
shown in Fig.(10-i).

##  Appendix D Computation of Matrix elements

###  D.1 c-c Coupling

The $c-c$ coupling is given as

$\displaystyle H^{\tau}_{cc}$ $\displaystyle=$
$\displaystyle\sum_{k\in[0,1)\otimes[0,\frac{1}{q})}\sum_{a,a^{\prime}=1}^{4}\sum_{m=0}^{m_{a,\tau}}\sum_{m^{\prime}=0}^{m_{a^{\prime},\tau}}\sum_{r,\tilde{r}=0}^{p-1}\tilde{h}^{\tau}_{[amr],[a^{\prime}m^{\prime}\tilde{r}]}(k)c^{\dagger}_{a\tau
krm}c_{a^{\prime}\tau k\tilde{r}m^{\prime}},$ (191)

where

$\displaystyle\tilde{h}^{\tau}_{[amr],[a^{\prime}m^{\prime}r^{\prime}]}(k_{1},k_{2})$
$\displaystyle=$ $\displaystyle\int
d^{2}\mathbf{r}\Psi^{*}_{a\tau}(\mathbf{r})\chi^{*}_{krm}(\mathbf{r})H^{\tau}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right)\Psi_{a^{\prime}\tau}(\mathbf{r})\chi_{kr^{\prime}m^{\prime}}(\mathbf{r}).$
(192)

Since the BM Hamiltonian is linear in $\mathbf{p}-\frac{e}{c}\mathbf{A}$, we
can write it as

$\displaystyle H^{\tau}_{BM}\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)$
$\displaystyle=$
$\displaystyle\mathcal{M}^{\tau}_{\mu}\left(p_{\mu}-\frac{e}{c}A_{\mu}(\mathbf{r})\right)+\mathcal{T}^{\tau}(\mathbf{r})$
(193)

where $\mathcal{M}^{+1}_{\mu}=v_{F}\mathbb{I}_{2}\otimes\mathbf{\sigma}$,
$\mathcal{M}^{-1}_{\mu}=-v_{F}\mathbb{I}_{2}\otimes\mathbf{\bar{\sigma}}$ with
$\mathbf{\sigma}=(\sigma_{x},\sigma_{y})$ and
$\bar{\sigma}=(\sigma_{x},-\sigma_{y})$. $\mathcal{T}^{\tau}$ denotes the
remaining factors in BM Hamiltonian. So,

$\displaystyle
H^{\tau}_{BM}\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)\Psi_{a^{\prime}\tau}(\mathbf{r})\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})=\mathcal{M}^{\tau}_{\mu}\left(p_{\mu}-\frac{e}{c}A_{\mu}(\mathbf{r})\right)\Psi_{a^{\prime}}(\mathbf{r})\chi_{k_{1},k_{2}+\frac{r^{\prime}}{q},m^{\prime}}(\mathbf{r})+\mathcal{T}^{\tau}(\mathbf{r})\Psi_{a^{\prime}\tau}(\mathbf{r})\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})$
(194)
$\displaystyle=\left(\mathcal{M}^{\tau}_{\mu}p_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})+\mathcal{T}^{\tau}(\mathbf{r})\Psi_{a^{\prime}\tau}(\mathbf{r})\right)\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})+\mathcal{M}^{\tau}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})\left(\left(p_{\mu}-\frac{e}{c}A_{\mu}(\mathbf{r})\right)\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})\right)$
(195)
$\displaystyle=\varepsilon^{\tau}_{a^{\prime\prime}a^{\prime}}\Psi_{a^{\prime\prime}\tau}(\mathbf{r})\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})+\mathcal{M}^{\tau}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})\left(p_{\mu}-\frac{e}{c}A_{\mu}(\mathbf{r})\right)\chi_{krm^{\prime}}(\mathbf{r}),$
(196)

where the matrix $\varepsilon^{\tau}$ is given in Eq.(213) is $\mathbf{k}$
independent because $\Psi_{a\tau}$ is defined at $\Gamma$. So,

$\displaystyle\tilde{h}^{\tau}_{[amr],[a^{\prime}m^{\prime}r^{\prime}]}(k_{1},k_{2})=\int
d^{2}\mathbf{r}\Psi^{*}_{a\tau}(\mathbf{r})\varepsilon^{\tau}_{a^{\prime\prime}a^{\prime}}\Psi_{a^{\prime\prime}\tau}(\mathbf{r})\chi^{*}_{krm}(\mathbf{r})\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})$
(197) $\displaystyle+\int
d^{2}\mathbf{r}\Psi^{*}_{a\tau}(\mathbf{r})\mathcal{M}^{\tau}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})\chi^{*}_{krm}(\mathbf{r})\left(p_{\mu}-\frac{e}{c}A_{\mu}(\mathbf{r})\right)\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})$
(198)

Since the factors
$\varepsilon^{\tau}_{a^{\prime\prime}a^{\prime}}\Psi^{*}_{a\tau}(\mathbf{r})\Psi_{a^{\prime\prime}\tau}(\mathbf{r})$
and
$\Psi^{*}_{a\tau}(\mathbf{r})\mathcal{M}^{\tau}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})$
are periodic with respect to primitive moiré lattice vectors $\mathbf{L}_{1}$
and $\mathbf{L}_{2}$, we can perform a Fourier expansion as

$\displaystyle\varepsilon^{\tau}_{a^{\prime\prime}a^{\prime}}\Psi^{*}_{a\tau}(\mathbf{r})\Psi_{a^{\prime\prime}\tau}(\mathbf{r})$
$\displaystyle=$
$\displaystyle\varepsilon^{\tau}_{a^{\prime\prime}a^{\prime}}\sum_{\mathbf{g}}e^{i\mathbf{g}\cdot\mathbf{r}}\left(\frac{1}{A_{tot}}\int
d^{2}\mathbf{r}^{\prime}e^{-i\mathbf{g}\cdot\mathbf{r}^{\prime}}\Psi^{*}_{a\tau}(\mathbf{r}^{\prime})\Psi_{a^{\prime\prime}\tau}(\mathbf{r}^{\prime})\right)$
(199)
$\displaystyle\Psi^{*}_{a\tau}(\mathbf{r})\mathcal{M}^{\tau}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{g}}e^{i\mathbf{g}\cdot\mathbf{r}}\left(\frac{1}{A_{tot}}\int
d^{2}\mathbf{r}^{\prime}e^{-i\mathbf{g}\cdot\mathbf{r}^{\prime}}\Psi^{*}_{a\tau}(\mathbf{r}^{\prime})\mathcal{M}^{\tau}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r}^{\prime})\right)$
(200)

As shown in previous section, at small $B$ the overlaps of the MTG LL with
$e^{i\mathbf{g}\cdot\mathbf{r}}$ are suppressed by a factor of
$\exp\left(-\mathbf{g}^{2}\ell^{2}/4\right)$. So at small enough $B$ (large
enough $\ell$) the overlaps with non-zero $\mathbf{g}$ components are
negligible. Keeping only the $\mathbf{g}=0$ Fourier component then gives

$\displaystyle\tilde{h^{\tau}}_{[amr],[a^{\prime}m^{\prime}r^{\prime}]}(k_{1},k_{2})\approx\left(\frac{1}{A_{tot}}\int
d^{2}\mathbf{r}\Psi^{*}_{a\tau}(\mathbf{r})\varepsilon^{\tau}_{a^{\prime\prime}a^{\prime}}\Psi_{a^{\prime\prime}\tau}(\mathbf{r})\right)\left(\int
d^{2}\mathbf{r}\chi^{*}_{km}(\mathbf{r})\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})\right)+$
$\displaystyle\left(\frac{1}{A_{tot}}\int
d^{2}\mathbf{r}\Psi^{*}_{a\tau}(\mathbf{r})\mathcal{M^{\tau}}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})\right)\left(\int
d^{2}\mathbf{r}\chi^{*}_{krm}(\mathbf{r})\left(p_{\mu}-\frac{e}{c}A_{\mu}(\mathbf{r})\right)\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})\right)+\mathcal{O}\left(\exp\left(-\mathbf{g}^{2}\ell^{2}/4\right)\right)$
(201)
$\displaystyle=\varepsilon^{\tau}_{aa^{\prime}}\delta_{rr^{\prime}}\delta_{mm^{\prime}}+\mathcal{M}^{\tau,\mu}_{aa^{\prime}}\left(\int
d^{2}\mathbf{r}\chi^{*}_{krm}(\mathbf{r})\left(p_{\mu}-\frac{e}{c}A_{\mu}(\mathbf{r})\right)\chi_{kr^{\prime}m^{\prime}}(\mathbf{r})\right)+\mathcal{O}\left(\exp\left(-\mathbf{g}^{2}\ell^{2}/4\right)\right),$
(202)

where

$\displaystyle\mathcal{M}^{\tau,\mu}_{aa^{\prime}}=\frac{1}{A_{tot}}\int
d^{2}\mathbf{r}\Psi^{*}_{a\tau}(\mathbf{r})\mathcal{M}^{\tau}_{\mu}\Psi_{a^{\prime}\tau}(\mathbf{r})$
(203)

The matrices $\epsilon^{\tau}_{a,a^{\prime}}$ and
$\mathcal{M}^{\tau,\mu}_{aa^{\prime}}$ can be obtained from the zero field
$cc$ coupling matrix

$\displaystyle
H^{c,\tau}_{aa^{\prime}}(\mathbf{k})=\left(\begin{array}[]{cc}0_{2\times
2}&v_{\ast}(\tau k_{x}\sigma_{0}+ik_{y}\sigma_{z})\\\ v_{\ast}(\tau
k_{x}\sigma_{0}-ik_{y}\sigma_{z})&M\sigma_{x}\end{array}\right)$ (206)

as

$\displaystyle M^{\tau,x}_{aa^{\prime}}=\frac{\tau
v_{\ast}}{\hbar}\left(\begin{array}[]{cc}0&\sigma_{0}\\\
\sigma_{0}&0\end{array}\right)_{aa^{\prime}}\;\;\text{,}\;\;M^{\tau,y}_{aa^{\prime}}=\frac{iv_{\ast}}{\hbar}\left(\begin{array}[]{cc}0&\sigma_{z}\\\
-\sigma_{z}&0\end{array}\right)_{aa^{\prime}}\;\;\text{,}\;\;\varepsilon^{\tau}_{aa^{\prime}}=M\left(\begin{array}[]{cc}0&0\\\
0&\sigma_{x}\end{array}\right)_{aa^{\prime}}.$ (213)

where the Pauli matrices $\sigma$ act in the orbital space of $c$ fermions.
Since