# Narrow bands in magnetic field and strong-coupling Hofstadter spectra

Xiaoyu Wang  National High Magnetic Field Laboratory, Tallahassee, Florida,
32310, USA  Oskar Vafek<EMAIL_ADDRESS>National High Magnetic Field
Laboratory, Tallahassee, Florida, 32310, USA  Department of Physics, Florida
State University, Tallahassee, Florida 32306, USA

###### Abstract

We develop a new, efficient, and general method to determine the Hofstadter
spectrum of isolated narrow bands. The method works for topological as well as
for trivial narrow bands by projecting the zero $\mathbf{B}$-field hybrid
Wannier states – which are localized in one direction and Bloch extended in
another direction – onto a representation of the magnetic translation group in
the Landau gauge. We then apply this method to find, for the first time, the
Hofstadter spectrum for the exact single particle charged excitations in the
strong coupling limit of the magic angle twisted bilayer graphene at the
charge neutrality point and at $|\nu|=2$ down to low magnetic fields when the
flux through the moiré unit cell is only $\sim 1/25$ of the electronic flux
quantum i.e. $\sim 1$T at the first magic angle. The resulting spectra provide
a means to investigate Landau quantization of the quasiparticles even if their
dispersion is interaction induced.

The rise of moiré materials [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16] has brought into focus the challenge to understand the physics of
correlated narrow bands subject to quantizing magnetic field $\mathbf{B}$ [17,
18, 19, 20, 21, 22, 23]. Such narrow bands can be topologically non-trivial
even at $\mathbf{B}=0$, as is the case for the magic angle twisted bilayer
graphene (MATBG)[24, 25, 26]. Moreover, for a moiré period $\sim 13$nm, as in
MATBG, the magnetic flux through the unit cell, $\phi$, can readily become
comparable or even exceed the flux quantum $\phi_{0}=hc/e$ using existing high
field magnets, so that the interplay of strong correlation and Hofstadter
physics can be realized in a laboratory [27, 28, 15, 29].

The traditional way to determine the non-interacting Hofstadter spectrum in
the MATBG is to minimally couple the magnetic vector potential $\mathbf{A}$ to
the continuum Bistritzer-MacDonald (BM) Hamiltonian [30] and then to expand it
in the Landau level (LL) basis [17, 18, 19, 20]. Although this provides a
reliable method, it requires a large upper cutoff on the LL index [19] in
order to converge, particularly at low $\mathbf{B}$, or close to simple
rational values of $\phi/\phi_{0}=p/q$ where the LL basis method becomes
prohibitively computationally expensive. This is because many Landau quantized
remote bands are effectively kept together with the Landau quantized narrow
bands of interest. Equivalently, at low $\mathbf{B}$, the real space shape of
the narrow band wavefunctions –with peaks in the local density of states at
the moiré triangular lattice sites– is mainly determined by the interlayer
tunneling ($w_{0,1}$) induced periodic potential and a superposition of a
large number of LLs is needed in order to recover such real space structure.
If one is then interested in interaction induced phenomena within the
resulting narrow bands a more efficient method is desirable.

The new method introduced here avoids the above mentioned difficulties. We
illustrate it at low $\mathbf{B}$, but the method is readily generalizable to
vicinity of simple fractions $p/q$. Thus, we first solve the $\mathbf{B}=0$
problem using standard (efficient) methods and find the hybrid Wannier states
for the $\mathbf{B}=0$ narrow bands [31, 26, 32, 33]. Such states are
exponentially localized in one direction and Bloch extended in another, say
$y$-direction [32] (see Fig. 1). We stress that even if the band is
topologically non-trivial, there is no obstruction to 1D exponential
localization. The key insight is that at $\mathbf{B}\neq 0$, for the hybrid
Wannier state (WS) centered at and near the origin, the Landau gauge vector
potential ${\bf A}=Bx\hat{{\bf y}}$ can be treated perturbatively, because the
region in real space where $\mathbf{A}$ is large gets suppressed by the
exponential localization of the hybrid WS (see Fig. 1). Moreover, the discrete
translation symmetry along the $y$-direction used in constructing the hybrid
WSs is preserved by such $\mathbf{A}$. Next, we generate the rest of the basis
by projecting the hybrid WSs centered at and near the origin onto a
representation of the magnetic translation group. This gives us two quantum
numbers, $k_{1}\in[0,1)$ and $k_{2}\in[0,1/q)$, associated with magnetic
translations by two non-collinear vectors $\mathbf{L}_{1}$ and
$q\mathbf{L}_{2}$ (Fig. 1). States with different $k_{1}$ and $k_{2}$ are then
guaranteed to be orthogonal. Because in the original ($\mathbf{B}=0$)
Brillouin zone $k_{2}$ belonged to a larger range $[0,1)$, we generate $q$
states for each starting hybrid WS at the same $k_{1}\in[0,1)$ and
$k_{2}\in[0,1/q)$ when $\mathbf{B}\neq 0$. Thus, for each $\mathbf{B}=0$
narrow band (of which there are two per valley and spin in MATBG) and for each
hybrid WS center described by a discrete index $n_{0}$, we have $q$ states.
The resulting states at the same $k_{1}\in[0,1)$ and $k_{2}\in[0,1/q)$ then
typically are not orthogonal, but by adjusting the range of $n_{0}$, the set
of states can be readily made overcomplete and thus span the $\mathbf{B}\neq
0$ narrow band Hilbert space. A simple procedure involving diagonalization of
the overlap matrix and keeping the $2q$ largest overlap eigenvalues (per spin
and valley) is then applied to obtain $2q$ orthogonal states within the MATBG
narrow bands at $\mathbf{B}\neq 0$. For MATBG and at low $\mathbf{B}$ we find
that the largest $2q$ overlap eigenvalues are clearly separated by a gap from
the remaining small eigenvalues, and that the $2q$ orthogonal states have an
almost perfect support by the $\mathbf{B}\neq 0$ narrow bands only (see Fig.
2b and 2c).

If we use this method on a topologically trivial narrow band, then a single
value of $n_{0}=0$ (i.e. only the hybrid WS centered at the origin) is
sufficient and none of the overlap eigenvalues become small even when
$\phi=\phi_{0}$.

Figure 1: Left: illustrative real space probability density of a hybrid
Wannier state $\ket{w_{c}(n_{0},k_{2}\mathbf{g}_{2})}$, with Chern index
$c=+1$, $n_{0}=0$ and $k_{2}=0$, and the Landau gauge magnetic vector
potential $\mathbf{A}=Bx\hat{\mathbf{y}}$. Moiré unit cell primitive vectors
are $\mathbf{L}_{1,2}$. Right: moiré Brillouin zone and reciprocal lattice
vectors $\mathbf{g}_{1,2}$. $K_{t,b}$ denote the Dirac point from the top and
bottom layers of the twisted bilayer graphene.

On the other hand, for the topologically non-trivial narrow bands of MATBG, we
need to keep at least two starting states with $n_{0}=0$ and $n_{0}=\pm 1$
(for either sign) in order to obtain complete orthogonal basis spanning the
$\mathbf{B}\neq 0$ narrow bands. This is a direct consequence of the non-
trivial topology of the $\mathbf{B}=0$ narrow band Hilbert space, spanned by a
band with Chern number $+1$ and a band with Chern number $-1$, one of which is
then deficient by $p$ anomalous sub-bands while the other has an excess of $p$
sub-bands when $\mathbf{B}\neq 0$ 111In Supplemental Material we provide a
derivation of $p$ anomalous magnetic sub-bands at flux $p/q$ in the chiral
limit, and provide numerical evidence both in and away from the chiral limit.
[35, 36]. We confirm this by studying the sublattice polarization of the
resulting states in Fig. 3 and analytical arguments in the chiral limit
presented in SM.

Our new basis can now be readily applied to finding the $\mathbf{B}\neq 0$
single electron or single hole excitation spectra in the strong coupling
problem by using the method introduced in Refs. [37, 38]. Note that even at
$\mathbf{B}\neq 0$, the 2-fold rotation about the out-of-plane axis $C_{2}$,
the particle-hole $P$ [26, 19] and the valley $U(1)$ conservation symmetries
of the BM Hamiltonian are preserved at any $w_{0}/w_{1}$; the time reversal
symmetry $T$ is of course broken by $\mathbf{B}$. Therefore, $C_{2}P$
guarantees that if $\Psi_{\mathbf{K},m,k_{1},k_{2}}(\mathbf{r})$ is an
eigenstate of
$\hat{H}^{\mathbf{K}}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right)$ defined via
Eq. (S70) below with an eigenvalue $E_{\mathbf{K},m,k_{1},k_{2}}$, then
$-i\mu_{y}\sigma_{x}e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}\Psi_{\mathbf{K},m,k_{1},k_{2}}(\mathbf{r})$
is an opposite valley eigenstate of
$\hat{H}^{\mathbf{K}^{\prime}}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right)$
with an eigenvalue $-E_{\mathbf{K},m,k_{1},k_{2}}$. The Pauli matrices
$\sigma$ and $\mu$ act in the sublattice and layer spaces, respectively.
Eliminating the remote magnetic sub-bands using the RG procedure introduced in
Ref. [37] therefore still results in the residual Coulomb interaction
projected onto the $\mathbf{B}\neq 0$ narrow band Hilbert space to be of the
form expressed in Eq. (9). Moreover, ignoring the Zeeman effect, $C_{2}P$
guarantees that the spin valley $U(4)$ symmetry [39, 40, 41, 23] is still
present even at $\mathbf{B}\neq 0$. We can therefore follow the double
commutator method outlined in Refs. [37, 38] in order to find the spectrum of
the single particle or single hole excitations at $\mathbf{B}\neq 0$. The
solutions of the Eq. (11) for two-gate screened Coulomb interaction,
$V_{\mathbf{q}}=\frac{2\pi
e^{2}}{\epsilon|\mathbf{q}|}\tanh\left(\frac{|\mathbf{q}|\xi}{2}\right)$, with
the gate separation $\xi=L_{m}$ are shown in the Fig. 3 for the charge neutral
point (CNP, i.e. $\nu=0$), together with their $\mathbf{B}=0$ density of
states. The results at $|\nu|=2$ for the heavy and light mass sides are shown
in the Fig.4. Below we provide details of the calculations which lead to the
stated results.

Figure 2: (a) Hofstadter spectrum for the non-interacting BM Hamiltonian
$\hat{H}_{BM}^{\mathbf{K}}(p_{x},p_{y}-\frac{eB}{c}x)$ calculated using Landau
level basis at magic angle ${w_{1}}/{v_{F}k_{\theta}}=0.586$ and
${w_{0}}/{w_{1}}=0.7$. The horizontal axes in (b) and (c) are the overlaps
between the $\mathbf{B}=0$ narrow bands hybrid Wannier states projected onto
the $\mathbf{B}\neq 0$ representation of the magnetic translation group
$|V_{a}\rangle$, and the exact magnetic subband states at various energies
obtained using the LL basis
$\sum_{a}|\langle\Psi_{\mathbf{K},m}(k_{1},k_{2})|V_{\mathbf{K},a}(k_{1},k_{2})\rangle|^{2}$;
the vertical axes are the sub-band index $m$. The states in between the dashed
lines belong to the $\mathbf{B}\neq 0$ narrow bands shown as blue in (a),
demonstrating that at low $\mathbf{B}$, $|V_{a}\rangle$ have support almost
exclusively within the $\mathbf{B}\neq 0$ narrow band Hilbert space.
${w_{0}}/{w_{1}}=0.7$ in (b) and (c) is at the chiral limit
${w_{0}}/{w_{1}}=0$.

To obtain the narrow band Hilbert space, we start by considering the BM model
at $\mathbf{B}\neq 0$ in Landau gauge
$\hat{H}^{\mathbf{K}}_{BM}\left(p_{x},p_{y}-\frac{eB}{c}x\right)$ where at the
valley $\mathbf{K}$

$\displaystyle\hat{H}^{\mathbf{K}}_{BM}(p_{x},p_{y})$ $\displaystyle=$
$\displaystyle\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\left(\mathbf{p}+\hbar\mathbf{q}_{1}\right)\end{array}\right).$
(3)

The Hamiltonian in valley $\mathbf{K}^{\prime}$ can be obtained by first
applying time reversal to $\hat{H}^{\mathbf{K}}_{BM}(p_{x},p_{y})$ followed by
the minimal substitution $p_{y}\rightarrow p_{y}-\frac{eB}{c}x$. The Pauli
matrices $\sigma$ act in the sublattice space 222We ignore the small rotation
of $\sigma$ matrices which was shown to lead to negligible effects on the
narrow band Hilbert space.. The interlayer hopping functions are
$T(\mathbf{r})=\sum_{j=1}^{3}T_{j}e^{-i\mathbf{q}_{j}\cdot\mathbf{r}}$ where
$\mathbf{q}_{1}=k_{\theta}(0,-1)$,
$\mathbf{q}_{2,3}=k_{\theta}\left(\pm\frac{\sqrt{3}}{2},\frac{1}{2}\right)$,
$k_{\theta}=\frac{8\pi}{3a_{0}}\sin\frac{\theta}{2}=4\pi/(3L_{m})$,
$a_{0}\approx 0.246$nm, $L_{m}$ is the period of the moiré lattice, and
$T_{j+1}=w_{0}1_{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)$,
where $1_{n}$ is an $n\times n$ unit matrix. At $\mathbf{B}=0$,
$\hat{H}^{\mathbf{K}}_{BM}$ is invariant under discrete translations by any
integer multiple of
$\mathbf{L}_{1}=L_{m}\left(\frac{\sqrt{3}}{2},\frac{1}{2}\right)$ and
$\mathbf{L}_{2}=L_{m}\left(0,1\right)$. At $\mathbf{B}\neq 0$ and in the
chosen gauge $\hat{H}_{BM}^{\mathbf{K}}$ 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\psi(\mathbf{r})\rightarrow\hat{t}_{\mathbf{L}_{1}}\psi(\mathbf{r})=e^{i\frac{eB}{\hbar
c}L_{1x}y}\psi(\mathbf{r}-\mathbf{L}_{1}).$ (4)

Thus, if $\psi(\mathbf{r})$ is an eigenstate then so is $e^{i\frac{eB}{\hbar
c}L_{1x}y}\psi(\mathbf{r}-\mathbf{L}_{1})$. Translations by $\mathbf{L}_{2}$
are generated by
$\hat{t}_{\mathbf{L}_{2}}\psi(\mathbf{r})=\psi(\mathbf{r}-\mathbf{L}_{2})$.
Then $\hat{t}_{\mathbf{L}_{2}}\hat{t}_{\mathbf{L}_{1}}=e^{-2\pi
i\phi/\phi_{0}}\hat{t}_{\mathbf{L}_{1}}\hat{t}_{\mathbf{L}_{2}}$, where
$\phi_{0}=\frac{hc}{e}$ and $\phi=BL_{1x}L_{m}$. If $\phi/\phi_{0}=p/q$, with
$p$ and $q$ relatively prime integers, then
$\left[\hat{t}^{q}_{\mathbf{L}_{2}},\hat{t}_{\mathbf{L}_{1}}\right]=0.$

The $\mathbf{B}=0$ hybrid WSs, $|w_{\pm}(n,k\mathbf{g}_{2})\rangle$, can be
chosen to be eigenstates of the periodic position operator
$\hat{O}=\hat{P}e^{-i\frac{1}{N_{1}}\mathbf{g}_{1}\cdot\mathbf{r}}\hat{P}$,
projected using $\hat{P}$ onto the $\mathbf{B}=0$ narrow band Hilbert space
studied (for details see Ref.[32]); here $N_{1}$ is a large integer. The
eigenvalues $e^{-2\pi i\frac{1}{N_{1}}\left(n+\langle
x_{\pm}\rangle_{k}/|\mathbf{L}_{1}|\right)}$ give the Wilson loops [31, 26,
32, 33] for the Chern $+1$ and Chern $-1$ hybrid WSs. These states are
localized along $\mathbf{L}_{1}$ and Bloch extended along $\mathbf{L}_{2}$, as
illustrated in the Fig. 1 As shown in Ref. [32], they satisfy,

$\displaystyle\hat{t}_{\mathbf{L}_{1}}|w_{\pm}(n,k_{2}\mathbf{g}_{2})\rangle$
$\displaystyle=$ $\displaystyle e^{i\frac{eB}{\hbar
c}L_{1x}y}|w_{\pm}(n+1,k_{2}\mathbf{g}_{2})\rangle$ (5)
$\displaystyle\hat{t}_{\mathbf{L}_{2}}|w_{\pm}(n,k_{2}\mathbf{g}_{2})\rangle$
$\displaystyle=$ $\displaystyle e^{-2\pi
ik_{2}}|w_{\pm}(n,k_{2}\mathbf{g}_{2})\rangle.$ (6)

We construct our basis for the narrow band at $\mathbf{B}\neq 0$ by projecting
$|w_{\pm}(n_{0},k_{2}\mathbf{g}_{2})\rangle$ onto representation of the
magnetic translation group (MTG). We include in our set a range of $n_{0}$’s
near $0$ as

$|W_{\pm}(k_{1},k_{2};n_{0})\rangle=\frac{1}{\sqrt{N}}\sum_{s=-\infty}^{\infty}e^{2\pi
isk_{1}}\hat{t}^{s}_{\mathbf{L}_{1}}|w_{\pm}(n_{0},k_{2}\mathbf{g}_{2})\rangle,$
(7)

with normalization factor $N$ and for $k_{1}\in[0,1)$ and temporarily letting
$k_{2}\in[0,1)$. The results in Figs.2b, 2c and 3 include $n_{0}=0$ and $1$.
Note that $|W_{\pm}(k_{1},k_{2};n_{0})\rangle$ are simultaneous 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. Thus the
$q\mathbf{L}_{2}$ translations break up the $k_{2}$ domain into $q$ pieces of
equal width $1/q$. Therefore, we let
$\left|W_{\pm}\left(k_{1},k_{2}+{l}/{q};n_{0}\right)\right\rangle$,
permanently fix $k_{2}\in[0,1/q)$, and let $l=0,1,\ldots q-1$. For different
values of $k_{1}$ and $k_{2}$ in their respective domains
$\left|W_{\pm}\left(k_{1},k_{2}+{l}/{q};n_{0}\right)\right\rangle$’s are
orthogonal because they have different eigenvalues under
$\hat{t}_{\mathbf{L}_{1}}$ and $\hat{t}^{q}_{\mathbf{L}_{2}}$. For the same
$k_{1}$ and $k_{2}$, but different $l$ (and different $n_{0}$) the states
$\left|W_{\pm}\left(k_{1},k_{2}+{l}/{q};n_{0}\right)\right\rangle$’s are in
general not orthogonal. To orthogonalize them we diagonalize the overlap
matrix $M_{ab}=\langle W_{a}|W_{b}\rangle=\left(U^{\dagger}DU\right)_{ab}$
where $D$ is diagonal. In the above we combined $l$, the Chern number index
$c=\pm$, and $n_{0}$ into a single index $a$ for each $k_{1}$ and $k_{2}$,
whose dependence we temporarily suppress. Then we let

$|V_{a}\rangle=\sum_{b}|W_{b}\rangle U^{\dagger}_{ba}\frac{1}{\sqrt{D_{a}}}$
(8)

where $b$ runs over all the indices but $a$ runs only over the $2q$ largest
eigenvalues $D_{a}$. As demonstrated in Fig. 2b and Fig. 2c, at low
$\mathbf{B}$, the $2q$ orthogonal states $|V_{a}(k_{1},k_{2})\rangle$ at each
$k_{1}$ and $k_{2}$ now form the basis spanning almost exclusively only the
$\mathbf{B}\neq 0$ narrow bands. At larger $\mathbf{B}$, we find a spillover
into the remote bands; for the results presented the spillover is negligible.

Figure 3: Landau level spectrum at magic angle in the strong coupling limit
for (a) $w_{0}/w_{1}=0$ and (b) $w_{0}/w_{1}=0.7$. The gray lines denote the
respective density of states $\mathcal{N}(E)$ at $\mathbf{B}=0$, and the color
of each magnetic sub-band denotes the average value of its sublattice
polarization. Blue (red) denotes purely A (B) sublattice polarization.

Next, we turn to the excitations in the strong coupling limit. In this case
the Hamiltonian consists of only the interaction
$V(\mathbf{r}-\mathbf{r}^{\prime})$ projected onto the $\mathbf{B}\neq 0$
narrow band basis. As described earlier, the $C_{2}P$ symmetry guarantees that
the dominant term takes the form

$\hat{H}_{int}=\frac{1}{2}\int d\mathbf{r}\int
d\mathbf{r}^{\prime}V(\mathbf{r}-\mathbf{r}^{\prime})\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}^{\prime}),$
(9)

where $\delta\rho(\mathbf{r})=\rho(\mathbf{r})-\bar{\rho}(\mathbf{r})$.
Restoring the indices on our $\mathbf{B}\neq 0$ narrow band basis functions
$\langle\mathbf{r}|V_{a}\rangle$, the projected density operator is
$\rho(\mathbf{r})=$

$\sum_{k_{1},k_{2},a}\sum_{k^{\prime}_{1},k^{\prime}_{2},a^{\prime}}V^{\dagger}_{\mathbf{K},a}(k_{1},k_{2};\mathbf{r})V_{\mathbf{K},a^{\prime}}(k^{\prime}_{1},k^{\prime}_{2};\mathbf{r})\mathfrak{d}^{\dagger}_{a,k_{1},k_{2}}\mathfrak{d}_{a^{\prime},k^{\prime}_{1},k^{\prime}_{2}}.$
(10)

We arranged the fermion creation operators with $2q$ discrete quantum numbers
$a$ and momentum $k_{1}$,$k_{2}$ into 4-component “spinor”,
$\mathfrak{d}^{\dagger}_{a,k_{1},k_{2}}=\left(d^{\dagger}_{\uparrow,\mathbf{K};a,k_{1},k_{2}},d^{\dagger}_{\downarrow,\mathbf{K};a,k_{1},k_{2}},d^{\dagger}_{\uparrow,\mathbf{K}^{\prime};C_{2}P[a,k_{1},k_{2}]},d^{\dagger}_{\downarrow,\mathbf{K}^{\prime};C_{2}P[a,k_{1},k_{2}]}\right)$.
The $U(4)$ manifold can be generated from a valley polarized state, which is
an eigenstate of $\rho(\mathbf{r})$ with the eigenvalue equal to
$\bar{\rho}(\mathbf{r})$ at CNP, where it takes the form, say,
$|\Phi_{\nu=0}\rangle=\prod_{a,k_{1},k_{2},\sigma=\uparrow,\downarrow}d^{\dagger}_{\sigma,\mathbf{K};a,k_{1},k_{2}}|0\rangle$.
Excitations can be created using an operator $X$ (see Ref. [37]) and their
strong coupling eigenenergies can be read off from the equation

$\displaystyle EX|\Phi_{\nu}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int
d^{2}\mathbf{r}d^{2}\mathbf{r}^{\prime}V(\mathbf{r}-\mathbf{r}^{\prime})\left[\rho(\mathbf{r}),\left[\rho(\mathbf{r}^{\prime}),X\right]\right]|\Phi_{\nu}\rangle$
(11) $\displaystyle+$ $\displaystyle\int
d^{2}\mathbf{r}d^{2}\mathbf{r}^{\prime}V(\mathbf{r}-\mathbf{r}^{\prime})\left[\rho(\mathbf{r}),X\right]\delta\bar{\rho}(\mathbf{r}^{\prime})|\Phi_{\nu}\rangle,$

where we extended the result to include $\nu=\pm 2$ fillings [38, 43]; the
valley polarized states $|\Phi_{\nu}\rangle$ are eigenstates of
$\delta\rho(\mathbf{r})$ with an eigenvalue $\delta\bar{\rho}(\mathbf{r})$
[44]. The eigenenergies of the strong coupling single particle or single hole
excitations can now be determined from diagonalizing a $2q\times 2q$ matrix
for each $k_{1}$ and $k_{2}$. Their degeneracy is determined by considering
the action of $X$ on $|\Phi_{\nu}\rangle$.

The resulting spectra at CNP are shown in the right panel of Fig. 3a for the
chiral limit $w_{0}/w_{1}=0$ and the right panel of Fig. 3b for
$w_{0}/w_{1}=0.7$; the spectra at $\nu=2$ are shown in Fig. 4. We clearly see
that despite being at strong coupling the excitations’ spectra are Landau
quantized in $\mathbf{B}\neq 0$. In the chiral limit (Fig. 3a), the degeneracy
of the low lying excitations limits to $4$ at low $\mathbf{B}$ due to spin and
sublattice degrees of freedom, the latter taking on values $\pm 1$ as marked
by the blue and red colors. Because they originate from $\mathbf{B}=0$ Chern
bands with opposite total Chern numbers, the $B$ sublattice sector has $q-1$
sub-bands while the $A$ sublattice sector has $q+1$ sub-bands for the $1/q$
sequence shown. Note that at small $\mathbf{B}$ there is a small splitting
between the low lying opposite sublattice polarized strong coupling sub-bands
due to broken $C_{2}T$ symmetry and that this splitting increases with
increasing $\mathbf{B}$. A similar conclusion has been reached in a recent
theoretical work [23], which reported energy splitting of the charge-$\pm 1$
excitations at full flux $\phi/\phi_{0}=1$. Also note the opposite evolution
of the sub-bands emanating from the $\mathbf{B}=0$ van Hove singularities.
Many of the features are reproduced at $w_{0}/w_{1}=0.7$, except the smaller
mean value of the sublattice polarization (as marked by the color scheme), and
larger splitting between the low lying magnetic sub-bands. Interestingly, the
sizable splitting between the light fermion LLs seen for $w_{0}/w_{1}=0.7$ in
Figs. 3b and 4d even at small $\phi/\phi_{0}$ would give rise to prominent LL
filling factors $|\nu_{LL}|=0,2$ at CNP, and $\nu_{LL}=0,1$ on the light mass
side of $\nu=2$, as observed in Ref. [15] without invoking moiré translational
symmetry breaking.

Figure 4: Landau level spectrum of charge $\pm 1$ excitations at $\nu=2$ at
magic angle in the strong coupling limit. Heavy hole (a,c) and light electron
(b,d) excitations for $w_{0}/w_{1}=0$ and $w_{0}/w_{1}=0.7$ respectively [38,
43]. The color scale for sublattice polarization is the same as in Fig. 3.

Published STM spectroscopy data [45] at $\mathbf{B}\neq 0$ show only results
from regions of various devices with the values of heterostrain $0.1\%-0.4\%$.
It is known that even such small values of strain dramatically increase the
non-interacting narrow band width [46, 47], making the kinetic energy
comparable or larger than the Coulomb interaction scale $e^{2}/\epsilon
L_{m}$, and stabilizing energetically proximate nematic state [48, 32, 47].
Therefore, the available STM data [45] at $\mathbf{B}\neq 0$ may not be in the
limit dominated by the Coulomb interaction complicating the direct comparison
with the strong coupling result presented here. The spectroscopic measurements
on magic angle devices at $\mathbf{B}\neq 0$ with negligible strain would
therefore be highly desirable.

The method presented here is general, and can be used to find the Hofstadter
spectrum at larger $\phi/\phi_{0}$ by starting with simple fractions
$1/\bar{q}$, where $\bar{q}$ is a small integer such as $1$,$2$ or $3$ and
where the LL based calculation is manageable, building the hybrid WSs for the
$2\bar{q}$ Hofstadter bands, and then projecting onto the representation of
the magnetic translation group for $\phi/\phi_{0}$ away from $1/\bar{q}$. Such
generalizations, as well as strain effects will be presented in future work.

###### Acknowledgements.

We thank B. Andrei Bernevig, Jonah Herzog-Arbeitman, Jian Kang for helpful
discussions. X.W. acknowledges financial support from National MagLab through
Dirac fellowship, which is funded by the National Science Foundation (Grant
No. DMR-1644779) and the state of Florida. O.V. was supported by NSF Grant No.
DMR-1916958.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Nature 572, 95 (2019).
  * 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, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019).
  * Jiang  _et al._ [2019] Y. Jiang, X. Lai, K. Watanabe, T. Taniguchi, K. Haule, J. Mao, and E. Y. Andrei, Nature 573, 91 (2019).
  * 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, Science 363, 1059 (2019).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Nature Physics 15, 1174 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * 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, Nature 572, 101 (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, E. Berg, P. Jarillo-Herrero, and S. Ilani, Nature 582, 203 (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, Nature 582, 198 (2020).
  * Serlin  _et al._ [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Science 367, 900 (2020).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Nature 583, 375 (2020).
  * Liu  _et al._ [2021a] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. I. A. Li, Science 371, 1261 (2021a).
  * 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, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Nature Physics 17, 1210 (2021).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Nature Materials 20, 488 (2021).
  * Bistritzer and MacDonald [2011a] R. Bistritzer and A. H. MacDonald, Phys. Rev. B 84, 035440 (2011a).
  * Moon and Koshino [2014] P. Moon and M. Koshino, Phys. Rev. B 90, 155406 (2014).
  * Hejazi  _et al._ [2019] K. Hejazi, C. Liu, and L. Balents, Phys. Rev. B 100, 035115 (2019).
  * Zhang  _et al._ [2019] Y.-H. Zhang, H. C. Po, and T. Senthil, Phys. Rev. B 100, 125104 (2019).
  * Lian  _et al._ [2020] B. Lian, F. Xie, and B. A. Bernevig, Phys. Rev. B 102, 041402 (2020).
  * Herzog-Arbeitman  _et al._ [2020] J. Herzog-Arbeitman, Z.-D. Song, N. Regnault, and B. A. Bernevig, Phys. Rev. Lett. 125, 236804 (2020).
  * Herzog-Arbeitman  _et al._ [2021] J. Herzog-Arbeitman, A. Chew, D. K. Efetov, and B. A. Bernevig, (2021), arXiv:2111.11434 .
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018).
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Phys. Rev. X 9, 021013 (2019).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Phys. Rev. Lett. 123, 036401 (2019).
  * Dean  _et al._ [2013] C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim, Nature 497, 598 (2013).
  * Saito  _et al._ [2021] Y. Saito, J. Ge, L. Rademaker, K. Watanabe, T. Taniguchi, D. A. Abanin, and A. F. Young, Nature Physics 17, 478 (2021).
  * Finney  _et al._ [2021] J. Finney, A. L. Sharpe, E. J. Fox, C. L. Hsueh, D. E. Parker, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, C. R. Dean, A. Vishwanath, M. Kastner, and D. Goldhaber-Gordon, (2021), arXiv:2105.01870 .
  * Bistritzer and MacDonald [2011b] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011b).
  * Yu  _et al._ [2011] R. Yu, X. L. Qi, A. Bernevig, Z. Fang, and X. Dai, Phys. Rev. B 84, 075119 (2011).
  * Kang and Vafek [2020] J. Kang and O. Vafek, Phys. Rev. B 102, 035161 (2020).
  * Soejima  _et al._ [2020] T. Soejima, D. E. Parker, N. Bultinck, J. Hauschild, and M. P. Zaletel, Phys. Rev. B 102, 205111 (2020).
  * Note [1] In Supplemental Material we provide a derivation of $p$ anomalous magnetic sub-bands at flux $p/q$ in the chiral limit, and provide numerical evidence both in and away from the chiral limit.
  * Popov and Milekhin [2021] F. K. Popov and A. Milekhin, Phys. Rev. B 103, 155150 (2021).
  * Sheffer and Stern [2021] Y. Sheffer and A. Stern, Phys. Rev. B 104, L121405 (2021).
  * Vafek and Kang [2020] O. Vafek and J. Kang, Phys. Rev. Lett. 125, 257602 (2020).
  * Bernevig  _et al._ [2021a] B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Phys. Rev. B 103, 205415 (2021a).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Phys. Rev. Lett. 122, 246401 (2019).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Phys. Rev. X 10, 031034 (2020).
  * Bernevig  _et al._ [2021b] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205413 (2021b).
  * Note [2] We ignore the small rotation of $\sigma$ matrices which was shown to lead to negligible effects on the narrow band Hilbert space.
  * Kang  _et al._ [2021] J. Kang, B. A. Bernevig, and O. Vafek, Phys. Rev. Lett. 127, 266402 (2021).
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Phys. Rev. B 103, 205414 (2021).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 588, 610 (2020).
  * Bi  _et al._ [2019] Z. Bi, N. F. Q. Yuan, and L. Fu, Phys. Rev. B 100, 035448 (2019).
  * Parker  _et al._ [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Phys. Rev. Lett. 127, 027601 (2021).
  * Liu  _et al._ [2021b] S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, Phys. Rev. Research 3, 013033 (2021b).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Bernevig  _et al._ [2021c] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205411 (2021c).
  * Vafek and Melikyan [2006] O. Vafek and A. Melikyan, Phys. Rev. Lett. 96, 167005 (2006).

Supplemental Materials for “Narrow bands in magnetic field and strong-coupling
Hofstadter spectra”

##  Appendix A Index theorem and analytical results in the chiral limit
$w_{0}/w_{1}=0$.

The non-interacting results in this section appeared previously in Refs.[35,
36]; the strong coupling results relevant for the main text are new. We
include the re-derivation of the former for convenience.

Let us adopt the symmetric gauge $\mathbf{A}=\frac{1}{2}B\left(-y,x,0\right)$.

$\displaystyle\hat{H}^{\mathbf{K}}_{BM}\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}v_{F}\sigma\cdot\left(\mathbf{p}-\frac{e}{c}\mathbf{A}\right)&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\left(\mathbf{p}+\mathbf{q}_{1}-\frac{e}{c}\mathbf{A}\right)\end{array}\right);$
(S3) $\displaystyle
v_{F}\sigma\cdot\left(\mathbf{p}+\mathbf{q}_{1}-\frac{e}{c}\mathbf{A}\right)$
$\displaystyle=$ $\displaystyle\hbar
v_{F}\left(\begin{array}[]{cc}0&\frac{1}{i}\frac{\partial}{\partial
x}-\frac{\partial}{\partial y}+ik_{\theta}+i\frac{eB}{2\hbar c}(x-iy)\\\
\frac{1}{i}\frac{\partial}{\partial x}+\frac{\partial}{\partial
y}-ik_{\theta}-i\frac{eB}{2\hbar c}(x+iy)&0\end{array}\right)$ (S6)
$\displaystyle=$ $\displaystyle\hbar
v_{F}\left(\begin{array}[]{cc}0&\frac{2}{i}\frac{\partial}{\partial
z}+ik_{\theta}+i\frac{eB}{2\hbar c}\bar{z}\\\
\frac{2}{i}\frac{\partial}{\partial\bar{z}}-ik_{\theta}-i\frac{eB}{2\hbar
c}z&0\end{array}\right)$ (S9)

where $z=x+iy$, $\bar{z}=x-iy$, $2\frac{\partial}{\partial
z}=\frac{\partial}{\partial x}-i\frac{\partial}{\partial y}$ and
$2\frac{\partial}{\partial\bar{z}}=\frac{\partial}{\partial
x}+i\frac{\partial}{\partial y}$. In the chiral limit,

$\displaystyle T(\mathbf{r})$ $\displaystyle=$ $\displaystyle
w_{1}\left(\begin{array}[]{cc}0&e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{-i\frac{2\pi}{3}}e^{-i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{i\frac{2\pi}{3}}e^{-i\mathbf{q}_{3}\cdot\mathbf{r}}\\\
e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{i\frac{2\pi}{3}}e^{-i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{-i\frac{2\pi}{3}}e^{-i\mathbf{q}_{3}\cdot\mathbf{r}}&0\end{array}\right).$
(S12)

In the above, $\hat{H}^{\mathbf{K}}_{BM}$ acts on
$\left(A_{top},B_{top},A_{bot},B_{bot}\right)$. Consider the unitary
transformation after which we have $\hat{H}^{\mathbf{K}}_{BM}$ act on
$\left(A_{top},A_{bot},B_{top},B_{bot}\right)$ as

$\displaystyle\hat{H}^{\mathbf{K}}_{BM}$ $\displaystyle\rightarrow$
$\displaystyle\left(\begin{array}[]{cc}0&\mathcal{D}^{\dagger}\\\
\mathcal{D}&0\end{array}\right),\;\;\mathcal{D}=\left(\begin{array}[]{cc}-i\hbar
v_{F}\left(2\frac{\partial}{\partial\bar{z}}+\frac{z}{2\ell^{2}}\right)&w_{1}U(\mathbf{r})e^{i\mathbf{q}_{1}\cdot\mathbf{r}}\\\
w_{1}U(-\mathbf{r})e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}&-i\hbar
v_{F}\left(2\frac{\partial}{\partial\bar{z}}+k_{\theta}+\frac{z}{2\ell^{2}}\right)\end{array}\right).$
(S17)

Here $\ell^{2}=\hbar c/(eB)$ and
$U(\mathbf{r})=e^{-i\mathbf{q}_{1}\cdot\mathbf{r}}+e^{i\frac{2\pi}{3}}e^{-i\mathbf{q}_{2}\cdot\mathbf{r}}+e^{-i\frac{2\pi}{3}}e^{-i\mathbf{q}_{3}\cdot\mathbf{r}}$.
Therefore, any state of the form

$\displaystyle
e^{-\frac{1}{4\ell^{2}}\bar{z}z}f(z)\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})\\\
0\end{array}\right)$ (S20)

and

$\displaystyle
e^{-\frac{1}{4\ell^{2}}\bar{z}z}f(z)\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})\\\
0\end{array}\right)$ (S23)

is a normalizable zero energy solution for an analytic $f(z)$, because
$\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})$ and
$\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})$ are the exact zero
energy states at $\mathbf{B}=0$; such a $\mathbf{B}=0$ state can always be
found even away from the magic angle [49]. Note that these states live
entirely on the $A$-sublattice and that unlike in the $\mathbf{B}=0$ (see
Ref.[49]), there is no normalizable solution on the $B$-sublattice.

Now, $f(z)\in\left(1,z,z^{2},\ldots,z^{N}\right)$, where $N+1=N_{\phi}$ is the
degeneracy of the Landau level. To show that the two states in Eq.(S20) are
linearly independent we need to show that the equation

$\displaystyle\sum_{n=0}^{N}c_{n}z^{n}e^{-\frac{1}{4\ell^{2}}\bar{z}z}\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})\\\
0\end{array}\right)+\sum_{n=0}^{N}c^{\prime}_{n}z^{n}e^{-\frac{1}{4\ell^{2}}\bar{z}z}\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})\\\
0\end{array}\right)=0$ (S28)

has a solution for all $\mathbf{r}$ only if $c_{n}=c^{\prime}_{n}=0$ for all
$n$. To do so, we note that even at $\mathbf{B}\neq 0$, we have
$C^{\prime}_{2}T=\mu_{x}K(x\rightarrow-x)$ symmetry and the unitary $P$
symmetry, whose combination also changes the sign of
$\hat{H}^{\mathbf{K}}_{BM}$:

$\displaystyle
PC^{\prime}_{2}T:\mu_{z}\hat{H}^{\mathbf{K}^{*}}_{BM}(x,-y)\mu_{z}=-\hat{H}^{\mathbf{K}}_{BM}(x,y).$
(S29)

This means that the zero modes can be chosen to be eigenstates of
$PC^{\prime}_{2}T$. To find out the parity of the zero modes at
$\mathbf{K}_{m}$ and $\mathbf{K}^{\prime}_{m}$, we note that at $\mathbf{B}=0$
the $\mathbf{k}$-points $\mathbf{K}_{m}$ and $\mathbf{K}^{\prime}_{m}$ are
related by the particle hole symmetry
$P=i\mu_{y}e^{i\mathbf{q}_{1}\cdot\mathbf{r}}\hat{\mathcal{I}}$ where
$\hat{\mathcal{I}}\psi(\mathbf{r})=\psi(-\mathbf{r})$. Therefore, we can
choose
$\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})=\hat{P}\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})$.
We see that $P$ anti-commutes with $PC^{\prime}_{2}T$ and does not change the
sublattice. Therefore, the parity of
$\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})$ and
$\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})=\hat{P}\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})$
must be opposite under $PC^{\prime}_{2}T$.

We can also see this explicitly from a perturbative “tripod model” solution of
Ref.[50]. Up to a normalization, the approximate zero mode of interest at
$\mathbf{K}_{m}$ is

$\displaystyle\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})=\sum_{\mathbf{g}}\left(\begin{array}[]{c}A^{top}_{n,\mathbf{g}}(-\mathbf{q}_{1})\\\
A^{bot}_{n,\mathbf{g}}(-\mathbf{q}_{1})\end{array}\right)e^{i\mathbf{g}\cdot\mathbf{r}}\approx\left(\begin{array}[]{c}1+e^{-i(\mathbf{g}_{1}+\mathbf{g}_{2})\cdot\mathbf{r}}+e^{-i\mathbf{g}_{2}\cdot\mathbf{r}}\\\
-i\frac{\hbar v_{F}k_{\theta}}{w_{1}}\end{array}\right)$ (S34)

and at $\mathbf{K}^{\prime}_{m}$ it is

$\displaystyle\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})=\sum_{\mathbf{g}}\left(\begin{array}[]{c}A^{top}_{n,\mathbf{g}}(0)\\\
A^{bot}_{n,\mathbf{g}}(0)\end{array}\right)e^{i\mathbf{g}\cdot\mathbf{r}}\approx\left(\begin{array}[]{c}i\frac{\hbar
v_{F}k_{\theta}}{w_{1}}\\\
1+e^{i(\mathbf{g}_{1}+\mathbf{g}_{2})\cdot\mathbf{r}}+e^{i\mathbf{g}_{2}\cdot\mathbf{r}}\end{array}\right)$
(S39)

where
$\mathbf{g}_{1}+\mathbf{g}_{2}=k_{\theta}\left(\frac{\sqrt{3}}{2},\frac{3}{2}\right)$
and $\mathbf{g}_{2}=k_{\theta}(-\frac{\sqrt{3}}{2},\frac{3}{2})$;
$k_{\theta}=\frac{4\pi}{3L_{m}}$. We see that

$\displaystyle
PC^{\prime}_{2}T:\;\;\mu_{z}\Psi^{*chiral}_{\mathbf{K}_{m}}(x,-y)$
$\displaystyle=$ $\displaystyle\Psi^{chiral}_{\mathbf{K}_{m}}(x,y)$ (S40)
$\displaystyle\;\;\mu_{z}\Psi^{*chiral}_{\mathbf{K}^{\prime}_{m}}(x,-y)$
$\displaystyle=$ $\displaystyle-\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(x,y),$
(S41)

in other words, they have opposite parity under $PC^{\prime}_{2}T$.

Because $PC^{\prime}_{2}T$ is a linear operator, going back to the equation
defining the linear independence, we can apply $PC^{\prime}_{2}T$ to both
sides of the Eq.(S28) to find

$\displaystyle
PC^{\prime}_{2}T:\;\;\sum_{n=0}^{N}c_{n}z^{n}e^{-\frac{1}{4\ell^{2}}\bar{z}z}\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})\\\
0\end{array}\right)-\sum_{n=0}^{N}c^{\prime}_{n}z^{n}e^{-\frac{1}{4\ell^{2}}\bar{z}z}\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})\\\
0\end{array}\right)=0.$ (S46)

Adding and subtracting Eqs.(S28) and (S46) we find

$\displaystyle\sum_{n=0}^{N}c_{n}z^{n}e^{-\frac{1}{4\ell^{2}}\bar{z}z}\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}_{m}}(\mathbf{r})\\\
0\end{array}\right)=0,$ (S49)
$\displaystyle\sum_{n=0}^{N}c^{\prime}_{n}z^{n}e^{-\frac{1}{4\ell^{2}}\bar{z}z}\left(\begin{array}[]{c}\Psi^{chiral}_{\mathbf{K}^{\prime}_{m}}(\mathbf{r})\\\
0\end{array}\right)=0.$ (S52)

Since the spinors and the gaussian factors are non-zero, the above hold only
if

$\displaystyle\sum_{n=0}^{N}c_{n}z^{n}=0,$ (S53)
$\displaystyle\sum_{n=0}^{N}c^{\prime}_{n}z^{n}=0.$ (S54)

But polynomials of different degrees are linearly independent, as can be seen
by taking successive derivatives and showing that the only way these equations
are satisfied for all $z$ is if each coefficient vanishes identically.

This proves that we have two Landau levels worth of zero modes in the chiral
limit at $\mathbf{B}\neq 0$ at a general twist angle. This number is
equivalent to having two exact zero modes for $k_{1}\in[0,1)$ and
$k_{2}\in[0,\frac{p}{q})$.

Because in the chiral limit

$\\{\hat{H}^{\mathbf{K}}_{BM},1_{2}\sigma_{z}\\}=0$ (S55)

and because the sublattice polarization eigevalues are $\pm 1$, by the index
theorem (see e.g. Refs.[36, 51]) we have

$\mbox{Tr}\left[\mathcal{P}1_{2}\sigma_{z}\right]=n_{+}-n_{-},$ (S56)

where $\mathcal{P}$ is the projector onto the narrow band Hilbert space, and
where $n_{+}$ is the number of zero energy modes with $1_{2}\sigma_{z}$
eigenvalue $+1$ and $n_{-}$ is the number of zero energy modes with
$1_{2}\sigma_{z}$ eigenvalue $-1$. But we found all of the zero energy modes
and they are sublattice $A$ polarized. Therefore, for a given $k_{1}\in[0,1)$
and $k_{2}\in[0,1/q)$ we have
$\mbox{Tr}\left[\mathcal{P}1_{2}\sigma_{z}\right]=2p$ for any finite
$\mathbf{B}$ and independent of the twist angle (as long as BM model applies).
Note that at $\mathbf{B}=0$,
$\mbox{Tr}\left[\mathcal{P}1_{2}\sigma_{z}\right]=0$ because for every
$A$-sublattice polarized zero mode there is a $B$-sublattice polarized zero
mode[49]. Therefore, $\mbox{Tr}\left[\mathcal{P}1_{2}\sigma_{z}\right]$ is
discontinuous at $\mathbf{B}=0$.

Figure S1: (a) Trace of sublattice polarization for the magic angle in chiral
limit $w_{0}/w_{1}=0$ (red) and for $w_{0}/w_{1}=0.7$ (blue) as a function of
magnetic flux $\phi$ through the moire unit cell in units of the flux quantum
$\phi_{0}=hc/e$, where $\mathcal{P}$ is the projector onto the narrow bands in
magnetic field. Due to the $C_{2}T$ symmetry at $\mathbf{B}=0$,
Tr$[\mathcal{P}1_{2}\sigma_{z}]$ must vanish at $\phi/\phi_{0}=0$; note the
discontinuous jump of the trace of the projected sublattice polarization at
non-zero $\phi/\phi_{0}$. (b) Spectrum of charge $\pm 1$ excitations at CNP in
the chiral limit for the $\phi/\phi_{0}=1/q$ sequence, for a single momentum
$\mathbf{k}=(0,0)$ in the magnetic Brillouin zone. There are $q+1$ A
sublattice polarized (blue) states and $q-1$ B sublattice polarized (red)
states.

We demonstrate the above result numerically in the Fig. S1(a). As seen,
although the sublattice polarization is no longer perfect at
$\frac{w_{0}}{w_{1}}\neq 0$, the discontinuity at $\mathbf{B}=0$ persists.

In the strong coupling chiral limit the effective Hamiltonian for the single
particle excitations commutes with $1_{2}\sigma_{z}$. Therefore, all of the
strong coupling magnetic sub-bands can be chosen to be eigenstates of
$1_{2}\sigma_{z}$ with eigenvalues $+1$ or $-1$. Since we have just proved
that the narrow band Hilbert space onto which we projected the interaction has
an extra $2p$ sublattice A polarized states at each $k_{1}\in[0,1)$ and
$k_{2}\in[0,1/q)$, we must have $q+p$ sublattice $A$ polarized magnetic sub-
bands and $q-p$ sublattice $B$ polarized sub-bands as demonstrated in Fig.
S1(b), where for clarity we show the spectrum for a single value of $k_{1}$
and $k_{2}$ along the $\phi/\phi_{0}=1/q$ sequence.

Throughout the main text and the supplementary, we have defined the “magic
angle” as the condition [49]:

$\frac{w_{1}}{v_{F}k_{\theta}}\equiv 0.586,\ w_{1}=96.056\text{meV},\
\frac{\hbar v_{F}}{a}=2135.4\text{meV}.$ (S57)

Here $a\approx 2.46\AA$ is the graphene lattice constant. These parameter
choices place the magic angle at $\theta=1.05^{\circ}$.

##  Appendix B Hybrid Wannier approach in finite magnetic field and numerical
procedure

###  B.1 Hybrid Wannier wavefunction representation of $\mathbf{B}=0$ narrow
band

In twisted bilayer graphene, hybrid Wannier wavefunctions have been
constructed as basis states for describing the narrow band physics of the
Bistritzer-MacDonald (BM) Hamiltonian at $\mathbf{B}=0$. Unlike exponentially
localized and symmetric Wannier states in both directions, there are no
topological obstructions to constructing hybrid Wannier states, which are
exponentially localized along one direction and Bloch extended along the
other. Detailed discussions of hybrid Wannier states have been given in Refs.
[26, 32]. Here we merely outline the general procedure for constructing hybrid
Wannier states. We begin by solving for the two narrow band energy eigenstates
of the BM Hamiltonian per valley and spin (Eq. (S70)), and construct them to
also be eigenstates of the $C_{2}T$ operator with eigenvalue $+1$; this fixes
their phase up to a sign. Next, the hybrid Wannier states
$\ket{w_{c}(n_{1},k_{2}\mathbf{g}_{2})}$ are constructed as eigenstates of the
projected (periodic) position operator
$\hat{P}e^{-i\frac{\mathbf{g}_{1}}{N_{1}}\cdot\mathbf{r}}\hat{P}$, where
$N_{1}$ is a large integer that discretizes the momentum space along the
$\mathbf{g}_{1}=\frac{4\pi}{\sqrt{3}L_{m}}(1,0)$ direction, and $\hat{P}$ is
the projector onto the $\mathbf{B}=0$ narrow bands. The hybrid Wannier states
are labeled by the Chern number $c=\pm 1$, momentum $k_{2}\in[0,1)$ along the
$\mathbf{g}_{2}\equiv\frac{4\pi}{\sqrt{3}L_{m}}(-\frac{1}{2},\frac{\sqrt{3}}{2})$
direction, and index $n_{1}\in\mathbb{Z}$ denoting the exponential
localization of the hybrid Wannier states near the real space strip
$\mathbf{r}=(n_{1}L_{1x},y)$. They can be represented as a one-dimensional
Fourier transform of the Chern Bloch states
$\ket{\Psi_{c}(k_{1}\mathbf{g}_{1},k_{2}\mathbf{g}_{2})}$, smooth in $k_{1}$,
as [32]:

$\ket{w_{c}(n_{1},k_{2}\mathbf{g}_{2})}=\frac{1}{\sqrt{N_{1}}}\sum_{k_{1}=0}^{1-\frac{1}{N_{1}}}e^{-i2\pi
k_{1}n_{1}}\ket{\Psi_{c}(k_{1}\mathbf{g}_{1},k_{2}\mathbf{g}_{2})},$ (S58)

where $k_{1},k_{2}\in[0,1)$. The hybrid Wannier states in the opposite valley
are related by time reversal symmetry.

###  B.2 Magnetic translation group eigenstates generated by hybrid Wannier
states

Figure S2: Eigenvalues of the overlap matrix
$\Lambda_{W}(c_{1}l_{1}n_{1},c_{2}l_{2}n_{2};k_{1},k_{2})$ for a fixed
$k_{1}=k_{2}=0$. Parameter choices are $q=13$ and $\theta=1.05^{\circ}$. (a,b)
correspond to the chiral limit $w_{0}/w_{1}=0$, and (c,d) correspond to
$w_{0}/w_{1}=0.7$.

In finite magnetic field, we choose the Landau gauge
$\mathbf{A}=|\mathbf{B}|x\hat{\mathbf{y}}$, and generate eigenstates of the
magnetic translation group (MTG) via:

$\ket{W_{c}(k_{1},k_{2};n_{0})}=\frac{1}{\sqrt{N}}\sum_{s_{1}=-\infty}^{\infty}e^{i2\pi
k_{1}s_{1}}\hat{t}_{\mathbf{L}_{1}}^{s_{1}}\ket{w_{c}(n_{0},k_{2}\mathbf{g}_{2})},$
(S59)

where $\hat{t}_{\mathbf{L}_{1,2}}$ are generators of magnetic translations by
the primitive vectors of the moire unit cell
$\mathbf{L}_{1}=L_{m}(\frac{\sqrt{3}}{2},\frac{1}{2})$ and
$\mathbf{L}_{2}=L_{m}(0,1)$. They are given by:

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

where $\hat{T}_{\mathbf{L}_{1,2}}$ are usual discrete translation operators,
defined via their action on a function $\psi(\mathbf{r})$ as
$\hat{T}_{\mathbf{L}_{1,2}}\psi(\mathbf{r})=\psi(\mathbf{r}-\mathbf{L}_{1,2})$,
and we have defined a magnetic translation wavevector:

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

For rational fluxes $\phi/\phi_{0}=p/q$, it is straightforward to show that
$\commutator{\hat{t}_{\mathbf{L}_{1}}}{\hat{t}_{\mathbf{L}_{2}}^{q}}=0$, and:

$\hat{t}_{\mathbf{L}_{1}}\ket{W_{c}(k_{1},k_{2};n_{0})}=e^{-i2\pi
k_{1}}\ket{W_{c}(k_{1},k_{2};n_{0})},\
\hat{t}_{\mathbf{L}_{2}}^{q}\ket{W_{c}(k_{1},k_{2};n_{0})}=e^{-i2\pi
qk_{2}}\ket{W_{c}(k_{1},k_{2};n_{0})}.$ (S62)

Note however that in general
$\commutator{\hat{t}_{\mathbf{L}_{1}}}{\hat{t}_{\mathbf{L}_{2}}}\neq 0$. We
can relabel the MTG eigenstates as:

$\ket{W_{c}(k_{1},k_{2}+\frac{l}{q};n_{0})},\ k_{1}\in[0,1),\
k_{2}\in[0,1/q),\ l=0,\dots q-1.$ (S63)

Here the magnetic Brillouin zone is defined by the magnetic strip
$[0,1)\times[0,1/q)$. MTG eigenstates labeled by different $k_{1}$,$k_{2}$
quantum numbers in the magnetic Brillouin zone are orthogonal. For a given
index $n_{0}$, there are $2q$ states generated from hybrid Wannier states, and
are labeled by the Chern number $c$ and the index $l$ of magnetic strips along
the $\mathbf{g}_{2}$ direction.

Due to the nontrivial band topology encoded in the hybrid Wannier
wavefunctions, the MTG eigenstates defined in Eq. (S59) are not guaranteed to
be orthonormal. Consider the overlap matrix at a given momentum
$k_{1}$,$k_{2}$:

$\begin{split}&\Lambda_{W}(c_{1}l_{1}n_{1},c_{2}l_{2}n_{2};k_{1},k_{2})\\\
\equiv&\bra{W_{c_{1}}(k_{1},k_{2}+\frac{l_{1}}{q};n_{1})}\ket{W_{c_{2}}(k_{1},k_{2}+\frac{l_{2}}{q};n_{2})}\\\
=&\sum_{s_{1}}e^{i2\pi
k_{1}s_{1}}\bra{w_{c_{1}}(n_{1},(k_{2}+\frac{l_{1}}{q})\mathbf{g}_{2})}\hat{t}_{\mathbf{L}_{1}}^{s_{1}}\ket{w_{c_{2}}(n_{2},(k_{2}+\frac{l_{2}}{q})\mathbf{g}_{2})}\\\
=&\sum_{s_{1}}e^{i2\pi
k_{1}s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}\bra{w_{c_{1}}(n_{1},(k_{2}+\frac{l_{1}}{q})\mathbf{g}_{2})}e^{is_{1}\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{c_{2}}(n_{2}+s_{1},(k_{2}+\frac{l_{2}}{q})\mathbf{g}_{2})}\\\
=&\frac{1}{N_{1}}\sum_{s_{1},\bar{k}_{1},\bar{p}_{1}}e^{i2\pi
k_{1}s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-i2\pi\bar{k}_{1}n_{1}}e^{i2\pi\bar{p}_{1}(n_{2}+s_{1})}\bra{\Psi_{c_{1}}(\bar{k}_{1}\mathbf{g}_{1},(k_{2}+\frac{l_{1}}{q})\mathbf{g}_{2})}e^{is_{1}\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{\Psi_{c_{2}}(\bar{p}_{1}\mathbf{g}_{1},(k_{2}+\frac{l_{2}}{q})\mathbf{g}_{2})}\\\
\equiv&\frac{1}{N_{1}}\sum_{s_{1},\bar{k}_{1},\bar{p}_{1}}e^{i2\pi
k_{1}s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-i2\pi\bar{k}_{1}n_{1}}e^{i2\pi\bar{p}_{1}(n_{2}+s_{1})}M(c_{1}\bar{\mathbf{k}},c_{2}\bar{\mathbf{p}};s_{1}),\end{split}$
(S64)

where on the last line for notational simplicity we have defined
$\bar{\mathbf{k}}=\bar{k}_{1}\mathbf{g}_{1}+(k_{2}+\frac{l_{1}}{q})\mathbf{g}_{2}$,
$\bar{\mathbf{p}}=\bar{p}_{1}\mathbf{g}_{1}+(k_{2}+\frac{l_{2}}{q})\mathbf{g}_{2}$,
and $M(c_{1}\bar{\mathbf{k}},c_{2}\bar{\mathbf{p}};s_{1})$ as the matrix
elements of $e^{is_{1}\mathbf{q}_{\phi}\cdot\mathbf{r}}$ in the Chern Bloch
basis. The magnetic wavevector $\mathbf{q}_{\phi}$ hybridizes Chern Bloch
states at different wavevectors satisfying the following two Diophantine
equations:

$\bar{k}_{1}=[\bar{p}_{1}+s_{1}\frac{p}{2q}]_{1},\ l_{1}=[l_{2}+s_{1}p]_{q},$
(S65)

where $[\dots]_{n}$ denotes the modulus with respect to integer $n$.

Figure S3: Measures of the goodness of hybrid Wannier approach on representing
the exact $\mathbf{B}\neq 0$ narrow band wavefunctions for a few values of
$w_{0}/w_{1}$ , twist angle $\theta$ at $\mathbf{k}=(0,0)$ of the magnetic
Brillouin zone. Top panel, the SVD eigenvalues of the overlap matrix
$\Lambda_{LL}(a,n;k_{1},k_{2})\equiv\bra{V_{a}(k_{1},k_{2})}\ket{\Psi_{n}(k_{1},k_{2})}$,
where $\ket{\Psi_{n}(k_{1},k_{2})}$ ($n=1\dots 2q$) are exact wavefunctions
calculated from the Landau-level wavefunction based approach. Bottom panel,
the weight of $\ket{V_{a}(k_{1},k_{2})}$ inside the narrow bands, defined as
$\text{Tr}[\Lambda_{LL}\Lambda_{LL}^{\dagger}]/2q$. Weight $=1$ indicates that
$\ket{V_{a}(k_{1},k_{2})}$ is related to the exact narrow band wavefunctions
via a unitary rotation. The exact $\mathbf{B}\neq 0$ narrow band eigenstates
using the Landau level approach are calculated with upper cutoff on the Landau
level index $n_{LL}=25q$.

We first discuss the completeness of the basis states. Note that if we were to
fix the indices $n_{1}=n_{2}=0$, the $2q$ MTG eigenstates are not linearly
independent. This is illustrated in Fig. S2(a,c), where we show the
eigenvalues of $\Lambda_{W}$ at momentum $k_{1}=k_{2}=0$. At
$\phi/\phi_{0}=p/q$, there are $2q-p$ linearly independent states. An
analytical proof of the state-deficiency is given in Sec. I of the SM for the
chiral limit $w_{0}/w_{1}=0$ using an index theorem. On the other hand, a
complete basis set can be generated by enlarging the $n_{1,2}$ range and
choosing $n_{1},n_{2}\in\\{0,1\\}$ for the trial MTG eigenstates. As
illustrated in Fig. S2(b,d), this procedure generates an overcomplete basis
set (i.e., number of independent states greater than $2q$). We therefore
choose $2q$ states with largest overlap eigenvalues to represent the narrow
band Hilbert space in a finite magnetic field. Note however, that only at low
magnetic fields are the $2q$ largest eigenvalues well separated from the
remaining $2q$ states by a well defined spectral gap, i.e. the spillover from
remote bands is small at low magnetic fields. However, at larger $\mathbf{B}$
the overlap spectral gap closes (e.g. Fig. S2(d) with
$\frac{\phi}{\phi_{0}}=\frac{2}{13}$), due to significant spillover from the
remote bands, making the $\mathbf{B}=0$ hybrid Wannier construction less
reliable. This is also seen in Fig. S4 where we make a quantitative comparison
of the Hofstadter spectra for the non-interacting BM Hamiltonian calculated
using the hybrid Wannier approach and more conventional Landau-level based
approach.

The desired orthonormalized basis set $\\{\ket{V_{a}(k_{1},k_{2})},\ a=1,\dots
2q\\}$ is thus obtained via:

$\ket{V_{a}(k_{1},k_{2})}\equiv\sum_{c,l,n_{0}}\ket{W_{c}(k_{1},k_{2}+\frac{l}{q};n_{0})}U(cln_{0},a;k_{1},k_{2})\frac{1}{\sqrt{D_{a}(k_{1},k_{2})}},$
(S66)

where $D_{a}$ are the largest $2q$ eigenvalues of $\Lambda_{W}$ for any given
$k_{1}$,$k_{2}$, and $U$ is a rectangular matrix satisfying:

$\sum_{c_{1}l_{1}n_{1},c_{2}l_{2}n_{2}}U^{*}(c_{1}l_{1}n_{1},a;k_{1},k_{2})\Lambda_{W}(c_{1}l_{1}n_{1},c_{2}l_{2}n_{2};k_{1},k_{2})U(c_{2}l_{2}n_{2},b;k_{1},k_{2})=\delta_{a,b}D_{a}(k_{1},k_{2}).$
(S67)

One can straightforwardly check the orthonormality condition
$\bra{V_{a}(k_{1},k_{2})}\ket{V_{b}(k_{1}^{\prime},k_{2}^{\prime})}=\delta_{k_{1},k_{1}^{\prime}}\delta_{k_{2},k_{2}^{\prime}}\delta_{a,b}$.

To quantify how well $\ket{V_{a}(k_{1},k_{2})}$ describes the $\mathbf{B}\neq
0$ narrow band Hilbert space, we follow Ref. [19] and obtain the narrow band
eigenstates $\ket{\Psi_{n}(k_{1},k_{2})}$ by expanding the BM Hamiltonian in
the Landau level basis of monolayer graphene. We define the overlap matrix
between states generated via these two procedures :

$\Lambda_{LL}(a,n;k_{1},k_{2})\equiv\bra{V_{a}(k_{1},k_{2})}\ket{\Psi_{n}(k_{1},k_{2})}.$
(S68)

If the hybrid Wannier approach generates exact eigenstates, then for each
$k_{1}$ and $k_{2}$ the matrix $\Lambda_{LL}(a,n;k_{1},k_{2})$ is unitary of
size $2q\times 2q$. In Fig. S3 we show both the SVD eigenvalues of
$\Lambda_{LL}$ and the projected weight
($\text{Tr}[\Lambda_{LL}\Lambda_{LL}^{\dagger}]/2q$) for a few twist angles
and ratios of $w_{0}/w_{1}$. Observe that as the magnetic field decreases, the
spillover of $\ket{V_{a}}$ into remote bands also decreases, and $\ket{V_{a}}$
extrapolates to the exact narrow band wavefunctions in the
$\mathbf{B}\rightarrow 0$ limit.

###  B.3 Matrix elements of the non-interacting BM Hamiltonian in the MTG
eigenstates

The matrix elements of the BM Hamiltonian in the orthonormalized MTG
eigenstate basis is given by:

$\begin{split}&\bra{V_{a}(k_{1},k_{2})}\hat{H}^{\mathbf{K}}_{BM}(\mathbf{p}-\frac{e}{c}\mathbf{A})\ket{V_{b}(k_{1},k_{2})}\\\
=&\frac{1}{\sqrt{D_{a}(k_{1},k_{2})}}U^{\dagger}(a,c_{1}l_{1}n_{1};k_{1},k_{2})H_{BM}(c_{1}l_{1}n_{1},c_{2}l_{2}n_{2};k_{1},k_{2})U(c_{2}l_{2}n_{2},b;k_{1},k_{2})\frac{1}{\sqrt{D_{b}(k_{1},k_{2})}},\end{split}$
(S69)

where repeated indices are summed over, and:

$\begin{split}&H_{BM}(c_{1}l_{1}n_{1},c_{2}l_{2}n_{2};k_{1},k_{2})\\\
\equiv&\bra{W_{c_{1}}(k_{1},k_{2}+\frac{l_{1}}{q};n_{1})}\hat{H}^{\mathbf{K}}_{BM}(\mathbf{p}-\frac{e}{c}\mathbf{A})\ket{W_{c_{2}}(k_{1},k_{2}+\frac{l_{2}}{q};n_{2})}\\\
=&\sum_{s_{1}}e^{i2\pi
k_{1}s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}\bra{w_{c_{1}}(n_{1},(k_{2}+\frac{l_{1}}{q})\mathbf{g}_{2})}\hat{H}^{\mathbf{K}}_{BM}(\mathbf{p}-\frac{e}{c}\mathbf{A})e^{is_{1}\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{c_{2}}(n_{2}+s_{1},(k_{2}+\frac{l_{2}}{q})\mathbf{g}_{2})}.\end{split}$
(S70)

The matrix elements of the BM Hamiltonian can be split into two terms

$\hat{H}^{\mathbf{K}}_{BM}(\mathbf{p}-\frac{e}{c}\mathbf{A})=\hat{H}^{\mathbf{K}}_{BM}(\mathbf{p})-\frac{\hbar
v_{F}}{\ell}(1_{2}\sigma_{y})\frac{x}{\ell},$ (S71)

where $\ell\equiv\sqrt{\frac{\hbar c}{eB}}$ is the magnetic length.

We stress that due to exponential localization of the hybrid Wannier states
along the $\mathbf{L}_{1}$ direction, the matrix elements of the BM
Hamiltonian in Eq. (S70) are non-neglibile only if $n_{1}$ and $n_{2}+s_{1}$
are close to each other. Since both $n_{1},n_{2}\in\\{0,1\\}$, this constrains
the summation over $s_{1}$ to a few moiré lattice constants. In practice we
are able to achieve numerical convergence for $s_{1}\in[-4,4]$. Therefore, $x$
in the matrix element is $O(L_{m})$ making the second term $O(\hbar
v_{F}L_{m}/\ell^{2})$. The first term in Eq. (S71) is nominally $O(\hbar
v_{F}k_{\theta})$ except very near the magic angle where there is an
additional suppression of the bandwidth by a factor of $\eta\sim 1/40$ (in the
chiral limit magic angle $\eta$ vanishes). The second term is therefore
nominally smaller than the first term at weak magnetic fields by a factor of
order $O(L^{2}_{m}/(4\ell^{2}))$, except near the magic angle where this
factor has an extra enhancement by $1/\eta$.

The matrix elements of the $\mathbf{B}=0$ term are calculated as follows:

$\begin{split}&H_{BM}^{(1)}(c_{1}l_{1}n_{1},c_{2}l_{2}n_{2};k_{1},k_{2})\\\
=&\frac{1}{N_{1}}\sum_{s_{1}\bar{k}_{1}\bar{p}_{1}}e^{i2\pi
k_{1}s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-i2\pi\bar{k}_{1}n_{1}}e^{i2\pi\bar{p}_{1}(n_{2}+s_{1})}\sum_{c_{3}}\varepsilon_{c_{1},c_{3}}(\bar{\mathbf{k}})M(c_{1}\bar{\mathbf{k}},c_{2}\bar{\mathbf{p}};s_{1}),\end{split}$
(S72)

where $\varepsilon_{c_{1},c_{3}}(\bar{\mathbf{k}})$ is the matrix elements of
the zero field BM Hamiltonian in the Chern Bloch basis, and $M$ is defined in
Eq. (S64).

The matrix elements for the vector potential term is calculated as:

$\begin{split}&H_{BM}^{(2)}(c_{1}l_{1}n_{1},c_{2}l_{2}n_{2};k_{1},k_{2})\\\
=&-\frac{\hbar
v_{F}}{\ell}\frac{1}{N_{1}}\sum_{s_{1}\bar{k}_{1}\bar{p}_{1}}e^{i2\pi
k_{1}s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-i2\pi\bar{k}_{1}n_{1}}e^{i2\pi\bar{p}_{1}(n_{2}+s_{1})}N(c_{1}\bar{\mathbf{k}},c_{2}\bar{\mathbf{p}};s_{1}),\end{split}$
(S73)

where we have defined:

$N(c_{1}\bar{\mathbf{k}},c_{2}\bar{\mathbf{p}};s_{1})\equiv\frac{1}{N_{1}L_{1,x}}\frac{1}{qN_{2}L_{2}}\int_{-N_{1}L_{1,x}/2}^{N_{1}L_{1,x}/2}\mathrm{d}x\int_{0}^{qN_{2}L_{2}}\mathrm{d}y\Psi_{c_{1}\bar{\mathbf{k}}}^{\dagger}(\mathbf{r})\left[\sigma_{y}\mu_{0}\frac{x}{\ell}e^{is_{1}\mathbf{q}_{\phi}\cdot\mathbf{r}}\right]\Psi_{c_{2}\bar{\mathbf{p}}}(\mathbf{r}).$
(S74)

Here we write it explicitly as a real space integral. We chose the momentum
space mesh to be $k_{1}=\frac{m_{1}}{N_{1}}$ where $m_{1}=0,\dots,N_{1}-1$,
and $k_{2}=\frac{m_{2}}{N_{2}q}$ where $m_{2}=0,\dots,N_{2}-1$. The
$x$-integration range in Eq. (S74) is necessary to avoid revivals of the
hybrid Wannier states since
$\ket{w_{c}(n+N_{1},k_{2}\mathbf{g}_{2})}=\ket{w_{c}(n,k_{2}\mathbf{g}_{2})}$.
Eq. (S74) is calculated by expressing the Chern Bloch states in the plane wave
basis, and perform real-space integration accordingly.

Figure S4: Comparison of the Hofstadter spectra of the non-interacting BM
Hamiltonian calculated via the hybrid Wannier approach (blue) and conventional
Landau level approach (gray). (a) is away from magic angle at
$\theta=1.38^{\circ}$, and (b) is at the magic angle.

In Fig. S4 we show the Hofstadter spectra calculated using the hybrid Wannier
approach in comparison with the exact spectra calculated using the Landau
level approach discussed in Ref. [19]. The spectra based on hybrid Wannier
approach are calculated for $\phi/\phi_{0}=p/83$, where $p=1\dots 20$. For the
Landau level approach, we only compute the $1/q$ sequence for $q=2\dots 18$.
To achieve numerical convergence we choose the upper cutoff for Landau levels
at $n_{LL}=25q$ (total number of Landau levels kept is $25q$ including the
zeroth Landau level). Observe that the Hofstadter spectra calculated using the
hybrid Wannier approach faithfully reproduce the exact results, with better
agreement as the magnetic field is reduced. At larger magnetic fields, the
spillover effect from the remote bands is significant, making the hybrid
Wannier approach unreliable. This occurs at $\frac{\phi}{\phi_{0}}\approx
0.157$ for parameters used in Fig. S4(b), which coincides roughly where the
spectral gap of the overlap matrix closes, shown in Fig. S2(d). Therefore, the
spectral gap of the overlap matrix provides a measure of how good the
$\mathbf{B}=0$ hybrid Wannier approach is in describing the narrow band
Hilbert space at $\mathbf{B}\neq 0$.

We note that due to the exponential localization of the hybrid Wannier states
along the $\mathbf{L}_{1}$-direction, in practice when calculating the matrix
elements in the hybrid Wannier approach, we restrict the summation over
$s_{1}$ from $-4$ to $4$. We checked that numerical convergence has been
achieved for $\Lambda_{W}$, matrix elements of BM Hamiltonian, as well as
matrix elements of the strong coupling Hamiltonian. Furthermore we choose a
momentum mesh such that $N_{1}=qN_{2}=2qn$ where we adjust integer $n$ to
achieve numerical convergence with the momentum mesh.

###  B.4 Matrix elements of the projected density operator

In momentum space, the Coulomb interaction takes the form:

$\hat{H}_{int}=\frac{1}{2}\sum_{\mathbf{q}}V(\mathbf{q})\delta\rho_{\mathbf{q}}\delta\rho_{-\mathbf{q}}$
(S75)

where the projected electron density operator in valley $\mathbf{K}$ is given
by:

$\hat{\rho}_{\mathbf{q}}=\sum_{a,k_{1},k_{2};b,p_{1},p_{2}}\bra{V_{a}(k_{1},k_{2})}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{V_{b}(p_{1},p_{2})}\mathfrak{d}^{\dagger}_{a,k_{1},k_{2}}\mathfrak{d}_{b,p_{1},p_{2}}$
(S76)

The density operator in valley $\mathbf{K}^{\prime}$ is related via $C_{2}P$
symmetry.

The background charge term $\bar{\rho}_{\mathbf{q}}$ can be calculated from
valley $\mathbf{K}$ (guaranteed by $C_{2}P$ symmetry), and is given as:

$\bar{\rho}_{\mathbf{q}}=2\sum_{m,n\in\mathbb{Z}}\delta_{\mathbf{q},m\mathbf{g}_{1}+n\mathbf{g}_{2}}\sum_{a,k_{1},k_{2}}\bra{V_{a}(k_{1},k_{2})}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{V_{a}(k_{1},k_{2})}$
(S77)

where $2$ comes from spin degeneracy (neglecting Zeeman splitting). Eqs. (S76)
and (S77) are calculated in a similar fashion to the matrix elements of the BM
Hamiltonian discussed previously. We choose $m,n\in\\{0,\pm 1,\pm 2,\pm 3\\}$
for numerical convergence.

##  Appendix C Solving for the dispersion of charged excitations in the strong
coupling limit using Landau-level based approach

###  C.1 MTG eigenstates generated from Landau level states

In the main text we have focused our attention on the hybrid Wannier approach
and solving for the $\mathbf{B}\neq 0$ dispersion of charged excitations. Here
we also present a calculation based on the Landau-level approach [17, 18, 19,
20], which is much more computationally demanding at low $\mathbf{B}$. Later
we provide a consistency check between Landau-level approach and the hybrid
Wannier approach.

We begin with a brief discussion of the Landau level eigenstates of the Dirac
Hamiltonian of monolayer graphene. For simplicity we consider the following
Dirac Hamiltonian in a magnetic field:

$\hat{H}^{\mathbf{K}}_{l}(\mathbf{p}-\frac{e\mathbf{A}}{c})=v_{F}\left[\sigma_{x}(p_{x}-K_{l,x})+\sigma_{y}(p_{y}-K_{l,y}-x/\ell^{2})\right].$
(S78)

Here $l=1,2$ is the layer index, and $\mathbf{K}_{l}=(K_{l,x},K_{l,y})$ is the
position of the Dirac cone in the reciprocal space. The eigenstates of the
Dirac Hamiltonian are solved by going to the harmonic oscillator basis:
$x=\frac{\ell}{\sqrt{2}}(a+a^{\dagger})$, and
$p_{x}=\frac{1}{i\sqrt{2}\ell}(a-a^{\dagger})$. The particle-hole symmetric
Landau level eigenstates are given as:

$\bra{\mathbf{r}}\ket{\psi^{(l)}_{n\gamma}(k_{2})}=e^{iK_{l,x}x}e^{i\frac{2\pi}{L_{m}}k_{2}y}\frac{1}{\sqrt{2}}\begin{pmatrix}\phi_{n}(x-\tilde{k}_{2,l}\ell^{2})\\\
-i\gamma\phi_{n-1}(x-\tilde{k}_{2,l}\ell^{2})\end{pmatrix},$ (S79)

where $\epsilon_{n\gamma}=\frac{v_{F}}{\ell}\gamma\sqrt{2n}$ is the energy of
the Dirac Hamiltonian, labeled by $n=1,2,\dots$, and $\gamma=\pm 1$
corresponds to positive and negative energy solutions. In addition, there is
an anomalous zero energy state given by

$\bra{\mathbf{r}}\ket{\psi^{(l)}_{0}(k_{2})}=e^{iK_{l,x}x}e^{i\frac{2\pi}{L_{m}}k_{2}y}\begin{pmatrix}\phi_{n}(x-\tilde{k}_{2,l}\ell^{2})\\\
0\end{pmatrix}$ (S80)

which lives on the A sublattice. $\phi_{n}(x)$ is the eigenfunction of
$a^{\dagger}a$, and is given to be:

$\phi_{n}(x)=\frac{1}{\pi^{1/4}}\frac{1}{\sqrt{2^{n}n!}}e^{-x^{2}/2\ell^{2}}H_{n}(x/\ell),$
(S81)

where $H_{n}(x)$ is the Hermite polynomial. The shift in the position for a
given momentum $k_{2}\mathbf{g}_{2}$ is given by:

$\tilde{k}_{2,l}\ell^{2}=\left(\frac{2\pi}{L_{m}}k_{2}-K_{l,y}\right)\ell^{2}.$
(S82)

In the Landau level basis, the eigenstates of the MTG are generated as:

$\begin{split}\ket{\Psi_{n\gamma}^{(l)}(k_{1},k_{2})}&=\sum_{s_{1}=-\infty}^{\infty}e^{i2\pi
k_{1}s_{1}}\hat{t}_{\mathbf{L}_{1}}^{s_{1}}\ket{\psi^{(l)}_{n\gamma}(k_{2})}\\\
&=\sum_{s_{1}}e^{i2\pi(k_{1}-\frac{k_{2}}{2})s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-is_{1}K_{l,x}L_{1,x}}\ket{\psi^{(l)}_{n\gamma}(k_{2}+s_{1}\frac{\phi}{\phi_{0}})}.\end{split}$
(S83)

It is straightforward to check that:

$\displaystyle\ket{\Psi_{n\gamma}^{(l)}(k_{1}+1,k_{2})}$
$\displaystyle=\ket{\Psi_{n\gamma}^{(l)}(k_{1},k_{2})},$ (S84)
$\displaystyle\ket{\Psi_{n\gamma}^{(l)}(k_{1},k_{2}+\frac{\phi}{\phi_{0}})}$
$\displaystyle=e^{iK_{l,x}L_{1,x}}e^{-i2\pi(k_{1}-\frac{k_{2}}{2})}\ket{\Psi_{n\gamma}^{(l)}(k_{1},k_{2})},$
(S85)
$\displaystyle\bra{\Psi_{n\gamma}^{(l)}(k_{1},k_{2})}\ket{\Psi_{n^{\prime}\gamma^{\prime}}^{(l^{\prime})}(p_{1},p_{2})}$
$\displaystyle=\delta_{l,l^{\prime}}\delta_{k_{1},k_{2}}\delta_{p_{1},p_{2}}\delta_{n,n^{\prime}}\delta_{\gamma,\gamma^{\prime}}.$
(S86)

Therefore, the MTG eigenstates defined in
$(k_{1},k_{2})\in[0,1)\times[0,\frac{\phi}{\phi_{0}})$ form a complete and
orthornomal basis set in a finite magnetic field.

The exact eigenstates of the narrow bands for non-interacting BM Hamiltonian
are solved by computing the matrix elements in these MTG eigenstates.

###  C.2 Strong coupling Hofstadter spectra using Landau level approach

We project the strong coupling Hamiltonian onto the $\mathbf{B}\neq 0$ narrow
band eigenstates, which we label as $\ket{\Psi_{a,\mathbf{k}}}$ where
$a=1,\dots 2q$, and $\mathbf{k}=k_{1}\mathbf{g}_{1}+k_{2}\mathbf{g}_{2}$. The
Hamiltonian can be written as:

$\displaystyle\hat{H}_{int}$ $\displaystyle=\frac{1}{2}\sum
V_{\mathbf{q}}\delta\rho_{\mathbf{q}}\delta\rho_{-\mathbf{q}},$ (S87)
$\displaystyle\delta\rho_{\mathbf{q}}$
$\displaystyle=\sum_{\eta=\mathbf{K},\mathbf{K}^{\prime}}\sum_{s=\uparrow,\downarrow}\sum_{a\mathbf{k},b\mathbf{p}}\bra{\Psi_{\eta
s,a,\mathbf{k}}}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi_{\eta
s,b,\mathbf{p}}}\mathfrak{d}^{\dagger}_{\eta s,a,\mathbf{k}}\mathfrak{d}_{\eta
s,b,\mathbf{p}}-\bar{\rho}_{\mathbf{q}},$ (S88)
$\displaystyle\bar{\rho}_{\mathbf{q}}$
$\displaystyle=\frac{1}{2}\sum_{m,n\in\mathbb{Z}}\delta_{\mathbf{q},m\mathbf{g}_{1}+n\mathbf{g}_{2}}\sum_{\eta,s}\sum_{a,\mathbf{k}}\bra{\Psi_{\eta
s,a,\mathbf{k}}}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi_{\eta
s,a,\mathbf{k}}}.$ (S89)

Here $1_{4}$ is the identity operator in the Hilbert space spanned by layer
and sublattice degrees of freedom. For completeness we added the valley
($\eta$) and spin ($s$) indices as subscripts to the narrow band eigenstates.
The background charge density term $\bar{\rho}_{\mathbf{q}}$ can be calculated
entirely in valley $\mathbf{K}$, due to the $C_{2}P$ symmetry relating the
eigenstate wavefunctions in the two valleys. To obtain the energy spectrum of
excitations at integer fillings we apply the double commutator method
discussed in the main text.

The key numerical procedure is to calculate the overlap matrix
$\bra{\Psi_{a,\mathbf{k}}}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi_{b,\mathbf{p}}}$,
which we discuss below. From here onward we also drop the valley and spin
indices for notational convenience. We proceed by constructing the matrix in
the MTG eigenstate basis:

$\bra{\Psi^{(l_{1})}_{n_{1}\gamma_{1},\mathbf{k}}}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi^{(l_{2})}_{n_{2}\gamma_{2},\mathbf{p}}},$
(S90)

and then project onto the narrow bands via:

$\bra{\Psi_{a,\mathbf{k}}}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi_{b,\mathbf{p}}}=U^{\dagger}_{a,n_{1}\gamma_{1}l_{1}}(\mathbf{k})\bra{\Psi^{(l_{1})}_{n_{1}\gamma_{1},\mathbf{k}}}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi^{(l_{2})}_{n_{2}\gamma_{2},\mathbf{p}}}U_{n_{2}\gamma_{2}l_{2},b}(\mathbf{p}),$
(S91)

where $U$ is a rectangular part of the unitary matrix that diagonalizes the BM
Hamiltonian. Repeated indices are summed over.

The matrix element defined in Eq. (S90) is calculated as follows:

$\begin{split}&\bra{\Psi^{(l_{1})}_{n_{1}\gamma_{1},\mathbf{k}}}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\Psi^{(l_{2})}_{n_{2}\gamma_{2},\mathbf{p}}}\\\
=&\sum_{s_{1},s_{1}^{\prime}}e^{-i2\pi k_{1}s_{1}}e^{i2\pi
p_{1}s_{1}^{\prime}}e^{i2\pi
q_{1}s_{1}}\bra{\psi^{(l_{1})}_{n_{1}\gamma_{1}}(k_{2})}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\hat{t}_{\mathbf{L}_{1}}^{s_{1}^{\prime}-s_{1}}\ket{\psi^{(l_{2})}_{n_{2}\gamma_{2}}(p_{2})}\\\
=&\sum_{s_{1}}\delta_{[p_{1}+q_{1}]_{1},k_{1}}e^{i2\pi
p_{1}s_{1}}\bra{\psi^{(l_{1})}_{n_{1}\gamma_{1}}(k_{2})}1_{4}e^{i\mathbf{q}\cdot\mathbf{r}}\hat{t}_{\mathbf{L}_{1}}^{s_{1}}\ket{\psi^{(l_{2})}_{n_{2}\gamma_{2}}(p_{2})}\\\
=&\delta_{l_{1},l_{2}}\sum_{s_{1}}\delta_{[p_{1}+q_{1}]_{1},k_{1}}e^{i2\pi(p_{1}-\frac{p_{2}}{2})s_{1}}e^{-i\frac{s_{1}(s_{1}-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{L}_{1}}e^{-is_{1}K_{l,x}L_{1,x}}\bra{\psi^{(l_{1})}_{n_{1}\gamma_{1}}(k_{2})}1_{2}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\psi^{(l_{1})}_{n_{2}\gamma_{2}}(p_{2}+s_{1}\frac{\phi}{\phi_{0}})}.\end{split}$
(S92)

Here we only calculate the expression when both Landau level indices $n_{1}$
and $n_{2}$ are non-zero. The case where either is zero can be calculated
straightforwardly. The expectation value of
$1_{2}e^{i\mathbf{q}\cdot\mathbf{r}}$ operator in the Landau level basis of
layer $l$ is given as follows:

$\begin{split}&\bra{\psi^{(l)}_{n_{1}\gamma_{1}}(k_{2})}1_{2}e^{i\mathbf{q}\cdot\mathbf{r}}\ket{\psi^{(l)}_{n_{2}\gamma_{2}}(p_{2}+s_{1}\phi/\phi_{0})}\\\
=&\frac{1}{2}\delta_{k_{2},p_{2}+s_{1}\frac{\phi}{\phi_{0}}+q_{2}}\int_{-\infty}^{\infty}\mathrm{d}x\begin{pmatrix}\phi_{n_{1}}(x-\tilde{k}_{2,l}\ell^{2}),&i\gamma_{1}\phi_{n_{1}-1}(x-\tilde{k}_{2,l}\ell^{2})\end{pmatrix}1_{2}e^{iq_{x}x}\begin{pmatrix}\phi_{n_{2}}(x-(\tilde{k}_{2,l}-q_{y})\ell^{2})\\\
-i\gamma_{2}\phi_{n_{2}-1}(x-(\tilde{k}_{2,l}-q_{y})\ell^{2})\end{pmatrix}\\\
=&\frac{1}{2}\delta_{k_{2},p_{2}+s_{1}\frac{\phi}{\phi_{0}}+q_{2}}\left[\int_{-\infty}^{\infty}\mathrm{d}x\phi_{n_{1}}(x-\tilde{k}_{2,l}\ell^{2})e^{iq_{x}x}\phi_{n_{2}}(x-(\tilde{k}_{2,l}-q_{y})\ell^{2})+\right.\\\
&\left.\gamma_{1}\gamma_{2}\int_{-\infty}^{\infty}\mathrm{d}x\phi_{n_{1}-1}(x-\tilde{k}_{2,l}\ell^{2})e^{iq_{x}x}\phi_{n_{2}-1}(x-(\tilde{k}_{2,l}-q_{y})\ell^{2})\right].\end{split}$
(S93)

Note that:

$\begin{split}&\int_{-\infty}^{\infty}\mathrm{d}x\phi_{n_{1}}(x-\tilde{k}_{2,l}\ell^{2})e^{iq_{x}x}\phi_{n_{2}}(x-(\tilde{k}_{2,l}-q_{y})\ell^{2})\\\
=&e^{iq_{x}\tilde{k}_{2,l}\ell^{2}}\int_{-\infty}^{\infty}\mathrm{d}x{\phi_{n_{1}}(x)}e^{iq_{x}x}{\phi_{n_{2}}(x+q_{y}\ell^{2})}\\\
=&e^{iq_{x}\tilde{k}_{2,l}\ell^{2}}\int_{-\infty}^{\infty}\mathrm{d}x\phi_{n_{1}}(x)e^{iq_{x}x}e^{ip_{x}q_{y}\ell^{2}}\phi_{n_{2}}(x)\\\
=&e^{iq_{x}\tilde{k}_{2,l}\ell^{2}}e^{-\frac{i}{2}q_{x}q_{y}\ell^{2}}\int_{-\infty}^{\infty}\mathrm{d}x\phi_{n_{1}}(x)e^{iq_{x}x+iq_{y}\ell^{2}p_{x}}\phi_{n_{2}}(x)\\\
=&e^{iq_{x}\tilde{k}_{2,l}\ell^{2}}e^{-\frac{i}{2}q_{x}q_{y}\ell^{2}}\int_{-\infty}^{\infty}\mathrm{d}x\phi_{n_{1}}(x)e^{c_{-}a+c_{+}a^{\dagger}}\phi_{n_{2}}(x),\end{split}$
(S94)

where we have used $e^{X}e^{Y}=e^{X+Y+\frac{1}{2}\commutator{X}{Y}+\dots}$,
and:

$c_{\pm}=i\frac{\ell}{\sqrt{2}}(q_{x}\pm iq_{y}).$ (S95)

The expectation value in the harmonic oscillator basis is calculated as
follows:

$\int_{-\infty}^{\infty}\mathrm{d}x{\phi_{n}(x)}e^{c_{-}a+c_{+}a^{\dagger}}{\phi_{m}(x)}=\begin{cases}e^{\frac{1}{2}c_{+}c_{-}}\sqrt{\frac{m!}{n!}}(c_{+})^{n-m}L_{m}^{n-m}(-c_{+}c_{-})&\text{
for }n\geq m,\\\
e^{\frac{1}{2}c_{+}c_{-}}\sqrt{\frac{n!}{m!}}(c_{-})^{m-n}L_{n}^{m-n}(-c_{+}c_{-})&\text{
for }n<m,\end{cases}$ (S96)

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

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

Figure S5: $|\bar{\rho}_{\mathbf{g}_{1}}|$ normalized by its $\mathbf{B}=0$
value, calculated using hybrid Wannier approach (gray) and Landau level
approach with various upper Landau level cutoffs.

###  C.3 Comparing Landau level approach and hybrid Wannier approach

As has been pointed out by Ref. [19], to achieve numerical convergence for the
Hofstadter spectra of the non-interacting BM Hamiltonian, the number of Landau
levels to be kept is roughly $25q$, when $\phi/\phi_{0}=1/q$. This makes the
calculation of strong coupling Hofstadter spectra prohibitively costly at low
magnetic fields. On the other hand, as has already been demonstrated earlier
in Fig. S4, in the hybrid Wannier approach, only two sets of Chern states are
needed, making it much more computationally viable.

In Fig. S5 we show the magnetic field dependence of
$|\bar{\rho}_{\mathbf{g}_{1}}|$ defined in Eq. (S89), calculated using the
Landau-level approach for a sequence of upper Landau level cutoffs $n_{LL}$,
as well as the hybrid Wannier approach. This shows that to achieve numerical
convergence in the Landau level approach a significant number of Landau levels
needs to be kept. Therefore, it is preferable to use the hybrid Wannier
approach at low magnetic fields.

Figure S6: Comparing the strong coupling Hofstadter spectra at magic angle
with $w_{0}/w_{1}=0$ (a-c) and $w_{0}/w_{1}=0.7$ (d-f). In every figure, the
left panel is the hybrid Wannier approach and right panel is Landau-level
based approach. The three columns are charge $\pm 1$ excitations at charge
neutral point, charge $+1$ at $\nu=2$, and charge $-1$ at $\nu=2$
respectively.

In Fig. S6 we compare the strong coupling Hofstadter spectra for these two
approaches at a single momentum point $(k_{1},k_{2})=(0,0)$. The hybrid
Wannier approach shows results for ${\phi}/{\phi_{0}}=1/4,\dots,1/25$. Due to
numerical stability and convergence issues at low fields, we only show the
Landau level approach calculation for ${\phi}/{\phi_{0}}=1/4,\dots,1/8$ with
upper Landau level cutoff of $25q$. Due to remote band spillover effects at
these values of $\phi/\phi_{0}$ for the hybrid Wannier approach (see for
instance Fig. S3), full quantitative comparison should not be expected.
Nevertheless, the two approaches display the same qualitative Landau
quantization of the strong coupling energy dispersions, including Landau level
degeneracies, the sublattice polarization, the opposite energetic evolution of
the sublattice polarized bands with magnetic field close to the van Hove
singularities of the $\mathbf{B}=0$ dispersions, thereby providing an
additional confirmation of the validity of the hybrid Wannier method
introduced in this work.