# Revisiting Bloch electrons in magnetic field: Hofstadter physics via hybrid
Wannier states

Xiaoyu Wang  National High Magnetic Field Lab, Tallahassee, FL 32310  Oskar
Vafek  National High Magnetic Field Lab, Tallahassee, FL 32310  Department of
Physics, Florida State University, Tallahassee, Florida 32306, USA

###### Abstract

We revisit the Hofstadter butterfly for a subset of topologically trivial
Bloch bands arising from a continuum free electron Hamiltonian in a periodic
lattice potential. We employ the recently developed procedure – which was
previously used to analyze the case of topologically non-trivial bands [Phys.
Rev. B 106, L121111 (2022)] – to construct the finite field Hilbert space from
the zero-field hybrid Wannier basis states. Such states are Bloch extended
along one direction and exponentially localized along the other. The method is
illustrated for square and triangular lattice potentials and is shown to
reproduce all the main features of the Hofstadter spectrum obtained from a
numerically exact Landau level expansion method. In the regime when magnetic
length is much longer than the spatial extent of the hybrid Wannier state in
the localized direction we recover the well known Harper equation. Because the
method applies to both topologically trivial and non-trivial bands, it
provides an alternative and efficient approach to moiré materials in magnetic
field.

##  I Introduction

The discovery of superconductivity and correlated insulating states in the
twisted bilayer graphene (TBG) [1, 2] has invigorated the study of various 2D
moiré electronic materials [3, 4, 5, 6, 7, 8, 9, 10]. The moiré superlattice
is generated by stacking 2D layered structures either with a small twist or
via microscopic lattice mismatch. Despite the large moiré unit cell containing
as many as $\sim 10,000$ atoms, the low energy physics is dominated by only a
few isolated narrow Bloch bands formed due to the moiré superlattice
potential, motivating theoretical studies that focus on these low energy
degrees of freedom [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22].

The large moiré unit cell has also enabled the study of magnetic field effects
in such systems in the regime with considerable fraction of one full magnetic
flux quantum per moiré unit cell. Electronic interaction effects intermixed
with strong magnetic fields have been studied in various moiré materials,
revealing not only Landau level degeneracies indicative of the symmetry-
breaking phases at zero magnetic field, but also novel field-induced
insulating states that can carry finite Chern numbers [23, 24, 25, 26, 27, 28,
29, 30, 31, 32]. These experimental results in turn motivate further
theoretical studies of Hofstadter physics [33, 34, 35, 36, 37, 38].

Traditionally, Bloch electrons in magnetic field $\mathbf{B}$ have been
studied either by solving the continuum Hamiltonian
$\hat{H}(\mathbf{r},\hat{\mathbf{p}}+\frac{e}{c}\mathbf{A})$, where the
magnetic vector potential satisfies $\mathbf{B}=\nabla\times\mathbf{A}$, or
via Peierl’s phase substitution of hopping amplitudes [39, 40, 41, 42]. The
latter is justified if the subset of Bloch bands of interest at zero field is
amenable to a tight-binding description i.e. there is no topological
obstruction to Wannierization. In the first approach, one calculates the
matrix elements of $\hat{H}$, including the periodic lattice potential, in the
Hilbert space spanned by Landau level (LL) wavefunctions (obtained without the
periodic lattice potential), and diagonalizes the resulting matrix. To achieve
numerical convergence within an energy window $W$ an upper LL index cutoff
$N_{c}\sim\lambda W/\hbar\omega_{c}$ is needed, where $\omega_{c}=eB/m_{e}c$
is the cyclotron frequency, $m_{e}$ is the bare electron mass, and $\lambda$
is a number that increases with the strength of the periodic lattice potential
$V$. For example, Ref. [33] has pointed out that an upper LL cutoff of
$N_{c}\sim 25\phi_{0}/\phi$ is needed to faithfully reproduce the narrow band
Hofstadter spectrum of TBG where $\phi_{0}=hc/e$ is the magnetic flux quantum,
and $\phi$ is the magnetic flux per moiré unit cell. This is a computationally
intensive procedure especially at low $\mathbf{B}$ when $N_{c}$ becomes large.
In the second approach, the magnetic field effects are accounted for via
Peierl’s substitution, i.e. replacing the intersite hopping amplitude $t_{ij}$
with $t_{ij}\exp\left(-i\frac{e}{\hbar
c}\int_{\mathbf{r}_{i}}^{\mathbf{r}_{j}}\mathrm{d}\mathbf{r}\cdot\mathbf{A}(\mathbf{r})\right)$.
The Peierl’s substitution has been used extensively in the literature due to
its simplicity in addressing Hofstadter physics. However, it is unclear how to
generalize the Peierl’s substitution to a subset of Bloch bands where 2D
exponentially localized and symmetric Wannier orbitals cannot be constructed
[43, 44, 45].

In an earlier work of ours [46], to address the Hofstadter physics in TBG, we
proposed a procedure for constructing the narrow band Hilbert space at a
rational magnetic flux ratio $\phi/\phi_{0}=p/q$ by projecting the zero-field
hybrid Wannier basis states (hWS) onto eigenstates of the magnetic translation
group (MTG). Such hWSs are Bloch extended along one direction and
exponentially localized along the other, and can always be constructed without
topological obstruction [47, 48, 49]. The pair of hWSs within a valley of TBG
carry $\pm 1$ Chern numbers, which are manifested in the intra-moiré-unit-cell
shift of the averaged position along the localization direction when the Bloch
wavenumber along the extended direction is changed [50]. We demonstrated that
the wavefunctions generated from this projection procedure have a good overlap
with the exact wavefunctions obtained using the LL approach at low
$\mathbf{B}$, while being a much more efficient numerical procedure than the
LL approach for addressing interaction effects [46].

In this work, we present a detailed discussion of the procedure developed in
Ref. [46] and apply it to revisit the non-interacting Hofstadter spectra for
square and triangular lattice potentials, where the lowest energy Bloch bands
at zero field are topologically trivial. We make quantitative comparisons to
the exact LL approach to demonstrate the projection method’s regime of
validity, as well as derive the Peierl’s substitution and Harper equation
studied extensively in the literature. The paper is organized as follows: In
section II we briefly discuss the LL approach and MTG eigenstates. In section
III we elaborate on the hWS and how to construct complete and orthonormal set
of MTG eigenstates for a subset of Bloch bands at a finite magnetic field. We
derive the Peierl’s substitution as a limiting case when the magnetic length
[defined below Eq. (1)] is much longer than the spatial support of the hWS,
and we rederive the Harper equation directly using the hWS approach. The
summary is provided in section IV.

##  II Landau level approach

We begin with a brief review of the Landau level (LL) approach for addressing
the magnetic field effects on electrons moving in a 2D periodic lattice
potential. For notational convenience, we introduce the two lattice vectors as
$\mathbf{a}_{1}$ and $\mathbf{a}_{2}$, and two reciprocal lattice vectors as
$\mathbf{g}_{1}$ and $\mathbf{g}_{2}$. They satisfy the relation
$\mathbf{a}_{i}\cdot\mathbf{g}_{j}=2\pi\delta_{ij}$. For a generic 2D lattice
$\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ are not required to be orthogonal [see
e.g. Fig. 1(d)]. We further work with the Landau gauge
$\mathbf{A}=Bx\mathbf{e}_{y}$ where $B$ is the magnetic field in the out-of-
plane $\mathbf{e}_{z}$-direction, and $\mathbf{e}_{y}$ ($\mathbf{e}_{x}$) is
along the direction parallel (perpendicular) to $\mathbf{a}_{2}$;
$\mathbf{e}_{i=x,y,z}$ are unit vectors. The LL wavefunctions are given by

$\ket{\psi_{n}(k_{y})}=\frac{1}{\sqrt{N_{2}a_{2}}}e^{ik_{y}y}\hat{T}(-k_{y}\ell^{2}\mathbf{e}_{x})\ket{n},$
(1)

where $N_{2}a_{2}$ is the length of the system along the $\mathbf{e}_{y}$
direction, $\ell\equiv\sqrt{\frac{\hbar c}{eB}}$ is the magnetic length,
$k_{y}\in\mathbb{R}$ is the momentum quantum number such that
$\hat{p}_{y}\ket{\psi_{n}(k_{y})}=\hbar k_{y}\ket{\psi_{n}(k_{y})}$, and
$\ket{n}$ is the $n$-th eigenstate for a 1D harmonic oscillator,

$\bra{\mathbf{r}}\ket{n}=\frac{1}{\pi^{1/4}\sqrt{2^{n}n!}}e^{-\frac{x^{2}}{2\ell^{2}}}H_{n}\left(\frac{x}{\ell}\right),$
(2)

where $H_{n}(x)$ is the Hermite polynomial. The operator
$\hat{T}(\mathbf{r}_{0})=e^{-i\mathbf{r}_{0}\cdot\hat{\mathbf{p}}/\hbar}$
generates a translation by $-\mathbf{r}_{0}$, i.e. for a general function
$f(\mathbf{r})$ we have

$\hat{T}(\mathbf{r}_{0})f(\mathbf{r})=f(\mathbf{r}-\mathbf{r}_{0}).$ (3)

The LL degeneracy argument proceeds as usual: consider a system of area
$N_{1}\mathbf{a}_{1}\times N_{2}\mathbf{a}_{2}$, such that it extends along
$\mathbf{a}_{i=1,2}$ by $N_{i=1,2}\in\mathbb{Z}$ unit cells. For open boundary
conditions along $\mathbf{a}_{1}$, the quantum number $k_{y}$ must satisfy
$(k_{y}\ell^{2})_{\text{max}}-(k_{y}\ell^{2})_{\text{min}}=N_{1}a_{1x}$. For
periodic boundary condition along $\mathbf{a}_{2}$, the separation between
adjacent wavevectors is $\delta k_{y}=\frac{2\pi}{N_{2}a_{2}}$. The total LL
degeneracy is then
$\mathcal{N}=\left[(k_{y})_{\text{max}}-(k_{y})_{\text{min}}\right]/\delta
k_{y}=(N_{1}N_{2})\frac{|\mathbf{a}_{1}\times\mathbf{a}_{2}|}{2\pi\ell^{2}}$,
or equivalently, the LL degeneracy per unit cell is given as

$\frac{\mathcal{N}}{N_{1}N_{2}}=\frac{\phi}{\phi_{0}},$ (4)

where $\phi_{0}=hc/e$ is the magnetic flux quantum, and
$\phi=B|\mathbf{a}_{1}\times\mathbf{a}_{2}|$ is the magnetic flux through a
unit cell.

###  II.1 Eigenstates of magnetic translation group

The single-electron Hamiltonian
$\hat{H}(\mathbf{r},\hat{\mathbf{p}}+\frac{e}{c}\mathbf{A})$ is invariant
under discrete magnetic translations
$\commutator{\hat{H}}{\hat{t}(\mathbf{a}_{i})}=0$, where
$\hat{t}(\mathbf{a}_{i=1,2})$ are discrete magnetic translation operators
along the $\mathbf{a}_{i=1,2}$ directions. In the Landau gauge, they are given
as:

$\displaystyle\hat{t}({\mathbf{a}_{1}})$
$\displaystyle=e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}(\mathbf{a}_{1}),$
(5) $\displaystyle\hat{t}(\mathbf{a}_{2})$
$\displaystyle=\hat{T}(\mathbf{a}_{2}),$ (6)

where we defined the wavevector associated with magnetic scattering

$\mathbf{q}_{\phi}=\frac{2\pi}{a_{2}}\frac{\phi}{\phi_{0}}\mathbf{e}_{y}=\frac{a_{1x}}{\ell^{2}}\mathbf{e}_{y}.$
(7)

These operators satisfy:

$\hat{t}({\mathbf{a}_{2}})\hat{t}({\mathbf{a}_{1}})=e^{i2\pi\frac{\phi}{\phi_{0}}}\hat{t}(\mathbf{a}_{1})\hat{t}(\mathbf{a}_{2}).$
(8)

Therefore, if ${\phi}/{\phi_{0}}={p}/{q}$ where $p$ and $q$ are coprime
integers,

$\commutator{\hat{t}(\mathbf{a}_{1})}{\hat{t}^{q}(\mathbf{a}_{2})}=0.$ (9)

As a result, eigenstates of
$\hat{H}(\mathbf{r},\hat{\mathbf{p}}+\frac{e}{c}\mathbf{A})$ can be chosen to
be simultaneous eigenstates of $\hat{t}(\mathbf{a}_{1})$ and
$\hat{t}^{q}(\mathbf{a}_{2})$.

It is straightforward to show that a complete and orthonormal set of MTG basis
states can be constructed from the LL wavefunctions as:

$\ket{\Psi_{n,r}(\mathbf{k})}=\frac{1}{\sqrt{N}}\sum_{s=-\infty}^{\infty}e^{i2\pi
k_{1}s}\hat{t}^{s}(\mathbf{a}_{1})\ket{\psi_{n}\left(\frac{2\pi}{a_{2}}(k_{2}+\frac{r}{q})\right)},$
(10)

where ${N}$ is a normalization factor. The wavevector
$\mathbf{k}=k_{1}\mathbf{g}_{1}+k_{2}\mathbf{g}_{2}$ resides in the magnetic
Brillouin zone defined as $k_{1}\in[0,1)$ and $k_{2}\in[0,\frac{1}{q})$, and
the integer $r$ labels the magnetic strip $[\frac{r-1}{q},\frac{r}{q})$ along
the $\mathbf{g}_{2}$ direction. The independent basis states are defined for
$r=0,\dots,p-1$, because (see appendix A):

$\ket{\Psi_{n,r+p}(\mathbf{k})}=e^{i2\pi\left(k_{1}-(k_{2}+\frac{r+p}{q})\frac{a_{1y}}{a_{2}}\right)}\ket{\Psi_{n,r}(\mathbf{k})}.$
(11)

It is straightforward to check that the states in Eq. (10) are MTG
eigenstates, i.e.,

$\displaystyle\hat{t}(\mathbf{a}_{1})\ket{\Psi_{n,r}(\mathbf{k})}$
$\displaystyle=e^{-i2\pi k_{1}}\ket{\Psi_{n,r}(\mathbf{k})},$ (12)
$\displaystyle\hat{t}^{q}(\mathbf{a}_{2})\ket{\Psi_{n,r}(\mathbf{k})}$
$\displaystyle=e^{-i2\pi qk_{2}}\ket{\Psi_{n,r}(\mathbf{k})}.$ (13)

They also satisfy the orthonormality condition:

$\bra{\Psi_{n_{1},r_{1}}(\mathbf{k})}\ket{\Psi_{n_{2},r_{2}}(\mathbf{p})}=\delta_{n_{1},n_{2}}\delta_{r_{1},r_{2}}\delta_{\mathbf{k},\mathbf{p}}.$
(14)

In the absence of periodic lattice potential, the LL degeneracy per unit cell
is given by $p\cdot 1\cdot\frac{1}{q}$ consistent with Eq. (4). Here $p$ comes
from the degeneracy of quantum number $r$, and a fully occupied magnetic
Brillouin zone corresponds to $1\cdot\frac{1}{q}$ of the zero-field Brillouin
zone occupied, i.e. the fraction of one particle per unit cell.

Note that $\hat{t}(\mathbf{a}_{2})$ acts non-trivially on the LL-based MTG
eigenstates, and because

$\hat{t}(\mathbf{a}_{2})\ket{\Psi_{n,r}(\mathbf{k})}=e^{-i2\pi\left(k_{2}+\frac{r}{q}\right)}\ket{\Psi_{n,r}(\mathbf{k}+\frac{p}{q}\mathbf{g}_{1})},$
(15)

$\hat{t}(\mathbf{a}_{2})\ket{\Psi_{n,r}(\mathbf{k})}$ is an MTG eigenstate at
wavevector $\mathbf{k}+\frac{p}{q}\mathbf{g}_{1}$.

Without loss of clarity from now on we use $\ket{\psi_{n}(k_{2}+\frac{r}{q})}$
for convenience, to denote a LL wavefunction with wavenumber
$\frac{2\pi}{a_{2}}\left(k_{2}+\frac{r}{q}\right)$.

###  II.2 Matrix elements of the Hamiltonian

We study the matrix elements of the continuum single electron Hamiltonian
given by:

$\hat{H}(\mathbf{r},\hat{\mathbf{p}}+\frac{e}{c}\mathbf{A})=\frac{(\hat{\mathbf{p}}+\frac{e}{c}\mathbf{A})^{2}}{2m_{e}}+V(\mathbf{r}),$
(16)

where $e>0$ is the electric charge, and
$V(\mathbf{r})=\sum_{\mathbf{g}}V_{\mathbf{g}}e^{i\mathbf{g}\cdot\mathbf{r}}$
is the periodic lattice potential, where
$\mathbf{g}=m\mathbf{g}_{1}+n\mathbf{g}_{2}$ are reciprocal lattice vectors,
$m,n\in\mathbb{Z}$.

The matrix elements of the kinetic energy in the LL-based MTG eigenstate basis
[Eq. (10)] can be straightforwardly calculated by rewriting
$\hat{\pi}_{x}=\frac{\hbar}{\sqrt{2}\ell}(a+a^{\dagger})$ and
$\hat{\pi}_{y}=\frac{\hbar}{i\sqrt{2}\ell}(a^{\dagger}-a)$, where
$\hat{\vec{\pi}}=\hat{\mathbf{p}}+\frac{e}{c}\mathbf{A}$ is the canonical
momentum, and $a$ is the harmonic oscillator lowering operator.

The matrix elements of a general operator of the form
$\hat{O}_{\mathbf{q}}=\mathcal{O}_{\mathbf{q}}e^{i\mathbf{q}\cdot\mathbf{r}}$
can be calculated as follows [46]:

$\begin{split}&(\hat{O}_{\mathbf{q}})_{nr_{1},mr_{2}}(\mathbf{k},\mathbf{p})\equiv\bra{\Psi_{n,r_{1}}(\mathbf{k})}\hat{O}_{\mathbf{q}}\ket{\Psi_{m,r_{2}}(\mathbf{p})}\\\
=&\mathcal{O}_{\mathbf{q}}\delta_{p_{1},[k_{1}-q_{1}]_{1}}\sum_{s=-\infty}^{\infty}\delta_{\tilde{p}_{y}-sq_{\phi},\tilde{k}_{y}-q_{y}}\\\
\times&e^{i2\pi
p_{1}s}e^{-is\tilde{p}_{y}a_{1y}}e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{-iq_{x}\tilde{k}_{y}\ell^{2}}e^{\frac{i}{2}q_{x}q_{y}\ell^{2}}\\\
\times&\bra{n}e^{c_{-}a+c_{+}a^{\dagger}}\ket{m},\end{split}$ (17)

where we have defined $c_{\pm}=i\frac{\ell}{\sqrt{2}}(q_{x}\mp iq_{y})$, and:

$\tilde{\mathbf{k}}=\mathbf{k}+\frac{r_{1}}{q}\mathbf{g}_{2},\
\tilde{\mathbf{p}}=\mathbf{p}+\frac{r_{2}}{q}\mathbf{g}_{2}.$ (18)

$\tilde{k}_{y}$ and $\tilde{p}_{y}$ are defined as
$\tilde{\mathbf{k}}\cdot\mathbf{e}_{y}$ and
$\tilde{\mathbf{p}}\cdot\mathbf{e}_{y}$, respectively. The notation $[b]_{a}$
represents $b$ modulo $a$, with $a>0$. The expression in the last line is
calculated as:

$\begin{split}&\bra{n}e^{c_{-}a+c_{+}a^{\dagger}}\ket{m}\\\
=&\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}\end{split}$ (19)

where

$L_{n}^{k}(x)=\sum_{m=0}^{n}(-x)^{m}\frac{(n+k)!}{(n-m)!(k+m)!m!}$

is the associated Laguerre polynomial. Note that any operator of the form
$e^{i\mathbf{q}\cdot\mathbf{r}}$ is a dense matrix in the LL indices
$\\{m,n\\}$. This poses numerical challenges at low magnetic flux ratios when
the upper LL cutoff is large.

The eigenstates and eigenenergies of the single electron Hamiltonian in a
magnetic field and periodic lattice potential can now be solved by
diagonalizing the matrix Hamiltonian in the LL-based MTG basis. We briefly
discuss the degeneracy of energy levels. Consider an energy eigenstate at a
momentum $\mathbf{k}$ inside the magnetic Brillouin zone, such that:

$\hat{H}\ket{\tilde{\Psi}_{n}(\mathbf{k})}=\varepsilon_{n,\mathbf{k}}\ket{\tilde{\Psi}_{n}(\mathbf{k})}.$
(20)

Making use of the magnetic translation operator $\hat{t}(\mathbf{a}_{2})$,
observe that

$\begin{split}\hat{H}\left(\hat{t}(\mathbf{a}_{2})\ket{\tilde{\Psi}_{n}(\mathbf{k})}\right)&=\hat{t}(\mathbf{a}_{2})\hat{H}\ket{\tilde{\Psi}_{n}(\mathbf{k})}\\\
&=\varepsilon_{n,\mathbf{k}}\left(\hat{t}(\mathbf{a}_{2})\ket{\tilde{\Psi}_{n}(\mathbf{k})}\right).\end{split}$
(21)

The state $\hat{t}(\mathbf{a}_{2})\ket{\tilde{\Psi}_{n}(\mathbf{k})}$ is
therefore also an energy eigenstate at $\varepsilon_{n,\mathbf{k}}$, but at
the wavevector $\mathbf{k}+\frac{p}{q}\mathbf{g}_{1}$. As a result each energy
level is at least $q$-fold degenerate, and therefore the dispersion of the
resulting magnetic subbands is effectively restricted to the wavevector domain
$[0,\frac{1}{q})\times[0,\frac{1}{q})$.

The LL approach is an exact method, limited in practice only by the truncation
of the upper LL index. In the free electron case, retaining $N_{c}$ LLs allows
an accurate representation of the energetics up to
$W\approx\frac{1}{\lambda}N_{c}\frac{\hbar^{2}}{m_{e}\ell^{2}}$, where
$\lambda>1$ is a parameter dependent on the strength of the lattice potential
as discussed in the introduction. In the low field limit, the LLs become
dense, and a larger LL index is therefore necessary to describe the magnetic
subbands emanating from the $\mathbf{B}=0$ Bloch bands up to the energy $W$.
This makes the LL approach both inefficient and not intuitive to study the low
field physics, and an alternative approach that bridges the zero field and
finite field Hilbert space is preferable. As we show in the next section, this
is achieved by projecting the hybrid Wannier basis states (hWSs) – which form
the basis of the $\mathbf{B}=0$ Hilbert space – onto representations of the
MTG.

##  III Hybrid Wannier approach

###  III.1 Hybrid Wannier states at $\mathbf{B}=0$

Figure 1:  Hybrid Wannier states $w_{n_{0},k_{2}}(\mathbf{r})$ of the lowest
energy Bloch band for a square lattice (a,b) and triangular lattice (c,d)
potential. In both cases we used the potential
$V(\mathbf{r})=-V_{0}\sum_{\mathbf{g}}e^{i\mathbf{g}\cdot\mathbf{r}}$. We set
$\hbar^{2}/m_{e}a^{2}=1$, where $m_{e}$ is the bare electron mass, and $a$ the
lattice constant. For square lattice,
$\mathbf{g}=\pm\mathbf{g}_{1},\pm\mathbf{g}_{2}$, whereas for triangular
lattice,
$\mathbf{g}=\pm\mathbf{g}_{1},\pm\mathbf{g}_{2},\pm(\mathbf{g}_{1}+\mathbf{g}_{2})$.
The states are well localized on the lattice sites, and their spatial support
narrows with increasing lattice potential strength.

In the absence of the magnetic field, the energy eigenstates are given by the
Bloch states

$\psi_{n,\mathbf{k}}(\mathbf{r})=\frac{1}{\sqrt{\mathcal{A}}}e^{i\mathbf{k}\cdot\mathbf{r}}u_{n,\mathbf{k}}(\mathbf{r}),$
(22)

where $\mathcal{A}$ is the area of the 2D system, $n$ is the band label,
$\mathbf{k}=k_{1}\mathbf{g}_{1}+k_{2}\mathbf{g}_{2}$ is the crystal momentum
in the first Brillouin zone [for convenience we choose it to be
$k_{1},k_{2}\in[0,1)\times[0,1)$], and $u_{n,\mathbf{k}}(\mathbf{r})$ is
periodic under discrete lattice translations $\hat{T}(\mathbf{a}_{i=1,2})$.
The notation $\psi_{n,\mathbf{k}}(\mathbf{r})$ should not be confused with the
LL wavefunctions discussed in the previous section. One can construct
spatially localized basis states by performing unitary transformations on the
(extended) Bloch energy eigenstates. For a subset of Bloch bands of interest,
we can construct hWSs from the Bloch states even if these bands have non-
trivial topology, provided that the gap to the adjacent bands does not close.
hWSs are exponentially localized in one direction (say $\mathbf{a}_{1}$), and
Bloch extended along the other (say $\mathbf{a}_{2}$). In Fig. 1 we show
several examples of hWSs for the lowest energy Bloch band on a square and
triangular lattice. These states can be constructed [48, 50] by diagonalizing
the periodic version of the position operator
$\exp(i\delta\mathbf{k}\cdot\mathbf{r})$ projected onto the Bloch basis from
the desired energy bands, where
$\delta\mathbf{k}=\frac{1}{N_{1}}\mathbf{g}_{1}$ hybridizes Bloch states at
$\mathbf{k}$ and $\mathbf{k}+\delta\mathbf{k}$. The resulting hWS that is
exponentially localized near a column at $n_{0}\mathbf{a}_{1}$ can be
expressed as:

$\ket{w_{\alpha}(n_{0},k_{2})}=\frac{1}{\sqrt{N_{1}}}\sum_{n\in{\rm
subset}}\sum_{k_{1}}e^{-i2\pi
k_{1}n_{0}}U_{n,\alpha}(\mathbf{k})\ket{\psi_{n,\mathbf{k}}},$ (23)

where the summation over $n$ is over a subset of Bloch bands of interest, and
$U(\mathbf{k})$ is a unitary matrix at every $\mathbf{k}$. Under discrete
translations the hWSs satisfy

$\displaystyle\hat{T}(\mathbf{a}_{1})\ket{w_{\alpha}(n_{0},k_{2})}$
$\displaystyle=\ket{w_{\alpha}(n_{0}+1,k_{2})},$ (24)
$\displaystyle\hat{T}(\mathbf{a}_{2})\ket{w_{\alpha}(n_{0},k_{2})}$
$\displaystyle=e^{-i2\pi k_{2}}\ket{w_{\alpha}(n_{0},k_{2})}.$ (25)

For topologically non-trivial bands, the hWSs contain information about non-
zero Chern number of a Bloch band, which is manifested in the non-trivial
evolution of the averaged position
$\langle\mathbf{r}\cdot\mathbf{g}_{1}\rangle/|\mathbf{g}_{1}|$ within the hWSs
when $k_{2}$ is continuously increased from 0 to 1. One example of such hWSs
is the pair of narrow bands in TBG for a given valley and spin. For more
details we refer interested readers to Refs. [50, 46].

###  III.2 MTG eigenstates from hybrid Wannier states

The set of hWSs for all Bloch bands forms a complete basis even in finite
magnetic field. However it is not useful if we are only interested in the
Hofstadter physics of a subset of Bloch bands.

It is tempting, but wrong, to take the subset of the $\mathbf{B}=0$ Bloch
bands as a basis of the corresponding subset of the $\mathbf{B}\neq 0$ states.
Note that the correct projector at $\mathbf{B}\neq 0$,
$\hat{P}_{B}=\sum_{n}^{\prime}\sum_{\mathbf{k}}\ket{\tilde{\Psi}_{n}(\mathbf{k})}\bra{\tilde{\Psi}_{n}(\mathbf{k})}$
[see Eq. (20)], where $n$ is summed over the subset of magnetic subbands of
interest, is invariant under any integer multiple of magnetic translations,
i.e.,
$\commutator{\hat{P}_{B}}{\hat{t}^{s_{1}}(\mathbf{a}_{1})\hat{t}^{s_{2}}(\mathbf{a}_{2})}=0$,
where $s_{1},s_{2}\in\mathbb{Z}$. However, the $\mathbf{B}=0$ projector,
$\hat{P}=\sum_{n\in\mathrm{subset}}\sum_{\mathbf{k}}\ket{\psi_{n,\mathbf{k}}}\bra{\psi_{n,\mathbf{k}}}$,
is not invariant under $\hat{t}^{s_{1}}(\mathbf{a}_{1})$, because

$\displaystyle\hat{t}^{s_{1}}(\mathbf{a}_{1})\hat{P}\hat{t}^{-s_{1}}(\mathbf{a}_{1})=$
(26)
$\displaystyle\sum_{mm^{\prime}\in\mathrm{fullset}}\sum_{\mathbf{k}}\ket{\psi_{m,\mathbf{k}}}\bra{\psi_{m^{\prime},\mathbf{k}}}\sum_{n\in\mathrm{subset}}U_{mn}(\mathbf{k})U_{nm^{\prime}}^{\dagger}(\mathbf{k}).$

Here we defined

$\hat{t}^{s_{1}}(\mathbf{a}_{1})\ket{\psi_{n,\mathbf{k}}}=\sum_{m\in{\rm
fullset}}U_{mn}(\mathbf{k})\ket{\psi_{m,\mathbf{k}}}.$ (27)

Due to the restriction on $n$, the right-hand-side of Eq. (26) is not equal to
$\hat{P}$, the narrow-band projector at $\mathbf{B}=0$. In other words, the
$y$-dependent phase in $\hat{t}(\mathbf{a}_{1})$ takes the states outside of
the subset of the zero-field bands of interest. This problem is severe even at
low $\mathbf{B}$ for sufficiently large $s_{1}$. Therefore $\hat{P}$ is not a
projector onto the states of interest in finite magnetic fields.

In order to construct the correct finite-field Hilbert space, we first
construct MTG eigenstates from the hWSs

$\ket{W_{\alpha,r}(\mathbf{k})}=\frac{1}{\sqrt{N_{1}}}\sum_{s=-\frac{N_{1}}{2}}^{\frac{N_{1}}{2}}e^{i2\pi
k_{1}s}\hat{t}^{s}({\mathbf{a}_{1}})\ket{w_{\alpha}(0,k_{2}+\frac{r}{q})}.$
(28)

Unlike in the previous subsection, here $k_{2}\in[0,\frac{1}{q})$, and
$r=0,\dots,q-1$. The choice of hWSs at $n_{0}=0$ is motivated by the fact that
the vector potential in the Landau gauge vanishes at the origin. Therefore,
$\ket{w_{\alpha}(0,k_{2}+\frac{r}{q})}$ must have a large overlap onto the
subbset of the $\mathbf{B}\neq 0$ Hilbert space with the similar energy at
small $B$, i.e. onto the magnetic subbands emanating from the $\mathbf{B}=0$
bands of interest. The rest of the basis can be conveniently obtained by using
the magnetic translation operator which moves hWSs along the $\mathbf{a}_{1}$
direction, while also attaching a phase due to the vector potential.

It is straightfoward to check that $\ket{W_{\alpha,r}(\mathbf{k})}$ are indeed
MTG eigenstates, i.e.,

$\displaystyle\hat{t}({\mathbf{a}_{1}})\ket{W_{\alpha,r}(\mathbf{k})}$
$\displaystyle=e^{-i2\pi k_{1}}\ket{W_{\alpha,r}(\mathbf{k})},$ (29)
$\displaystyle\hat{t}^{q}({\mathbf{a}_{2}})\ket{W_{\alpha,r}(\mathbf{k})}$
$\displaystyle=e^{-i2\pi qk_{2}}\ket{W_{\alpha,r}(\mathbf{k})}.$ (30)

Figure 2:  Typical absolute value of the matrix elements of the overlap matrix
[Eq. (32)], calculated at flux $\phi/\phi_{0}=1/10$ for (a) square and (b)
triangular lattice potentials respectively with $V_{0}=3$. We set the energy
scale $\hbar^{2}/m_{e}a^{2}=1$. (c) Eigenvalues of the overlap matrix. The
unlabeled axes in (a,b) correspond to the basis index $r,r^{\prime}$, and in
(c) the ordered eigenvalue index. If the MTG basis states generated from hWSs
constitute a complete and orthonormal basis, the overlap matrix is an identity
matrix and all eigenvalues are equal to 1. While this is true for square
lattice potential, there are small deviations from orthonormality for the case
of triangular lattice potential. As a result a further orthonormalization
procedure outlined in the text is needed.

To address the completeness and orthonormality of the wavefunctions defined in
Eq. (28), we define an overlap matrix

$\Lambda_{\alpha r_{1},\beta
r_{2}}(\mathbf{k})\equiv\bra{W_{\alpha,r_{1}}(\mathbf{k})}\ket{W_{\beta,r_{2}}(\mathbf{k})}.$
(31)

Note that states with different $\mathbf{k}$ are automatically orthogonal due
to different eigenvalues under $\hat{t}(\mathbf{a}_{i=1,2})$. If these
wavefunctions represent a complete and orthonormal set, the overlap matrix
should be an identity matrix. The overlap matrix can be calculated as follows:

$\begin{split}\Lambda_{\alpha r_{1},\beta r_{2}}(\mathbf{k})=&\sum_{s}e^{i2\pi
k_{1}s}\bra{w_{\alpha}(0,k_{2}+\frac{r_{1}}{q})}\hat{t}^{s}({\mathbf{a}_{1}})\ket{w_{\beta}(0,k_{2}+\frac{r_{2}}{q})}\\\
=&\sum_{s}e^{i2\pi
k_{1}s}e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}\bra{w_{\alpha}(0,\tilde{k}_{2})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{\beta}(s,\tilde{p}_{2})},\end{split}$
(32)

where on second line we defined $\tilde{k}_{2}=k_{2}+\frac{r_{1}}{q}$ and
$\tilde{p}_{2}=k_{2}+\frac{r_{2}}{q}$ ($\tilde{k}_{2},\tilde{p}_{2}\in[0,1)$),
and used the operator identity in Eq. (57).

If the hWSs represent topologically trivial Bloch bands (e.g., Fig. 1), the
overlap of two hWSs is exponentially suppressed unless they have the same
localization center, i.e. unless $s=0$ in Eq. (32). Therefore
$\bra{W_{\alpha,r_{1}}(\mathbf{k})}\ket{W_{\beta,r_{2}}(\mathbf{p})}\approx\delta_{\alpha,\beta}\delta_{\mathbf{k},\mathbf{p}}\delta_{r_{1},r_{2}}$
with exponential accuracy. This is illustrated in the Fig. 2 for the square
and triangular lattice potentials. Even if $\Lambda(\mathbf{k})$ is not an
identity matrix, a complete and orthonormal basis set can be generated by
eigen-decomposition: $U^{\dagger}\Lambda U=D$, and redefining a new set of
basis states as:

$\ket{V_{a}(\mathbf{k})}=\sum_{\alpha,r}\ket{W_{\alpha,r}(\mathbf{k})}\left(U\frac{1}{\sqrt{D}}\right)_{\alpha
r,a}.$ (33)

If hWSs represent topological bands with finite Chern numbers, hWSs with
different localization centers have a finite spatial overlap, and the overlap
matrix in Eq. (32) strongly deviates from an identity matrix. Importantly, as
discussed in Ref. [46], for hWSs with Chern numbers $\pm 1$, the number of
independent MTG eigenstates should increase/decrease according to the Streda
formula. This manifests in the rank deficiency of the overlap matrix. A
procedure that makes use of a different choice of the localization center
$n_{0}$ (still close to $0$ where the $\mathbf{A}$ is small) was proposed in
Ref. [46], where it was shown to resolve this issue.

Figure 3:  The amount of spillover $\eta$ of MTG basis states into the Hilbert
space of remote Hofstadter bands for square and triangular lattice potential.
We set the energy scale $\hbar^{2}/m_{e}a^{2}=1$. The spillover decreases with
increasing lattice potential strength, and increases with increasing magnetic
flux.

The next question is how well do these basis states describe finite field
Hilbert space of interest. To quantify the amount of spillover into remote
magnetic subbands, we can expand Eq. (33) in the exact wavefunctions obtained
based on the LL approach,

$\begin{split}\ket{V_{a}(\mathbf{k})}&=\sum_{n\in\mathrm{active}}M_{n,a}(\mathbf{k})\ket{\tilde{\Psi}_{n}(\mathbf{k})}\\\
&+\sum_{n^{\prime}\in\mathrm{remote}}M_{n^{\prime},a}(\mathbf{k})\ket{\tilde{\Psi}_{n^{\prime}}(\mathbf{k})},\end{split}$
(34)

where $\ket{\tilde{\Psi}_{n}(\mathbf{k})}$ are the exact eigenstates obtained
from the LL approach, and $M_{n,a}(\mathbf{k})$ are the expansion
coefficients. Define:

$\eta_{a}\equiv\frac{1}{N_{1}}\frac{1}{N_{2}/q}\sum_{\mathbf{k}}\sum_{n^{\prime}\in\mathrm{remote}}|M_{n^{\prime},a}(\mathbf{k})|^{2}.$
(35)

Then $\eta_{a}\in[0,1]$ characterizes the amount of spillover into remote
magnetic subbands. In Fig. 3 we plot $\eta=\eta_{a=1}$ versus $\phi/\phi_{0}$
for the lowest energy Bloch band for square and triangular lattice potentials.
The amount of spillover decreases with increasing strength of the lattice
potential, and remains small even for reasonably large magnetic flux.

###  III.3 Peierl’s factor and MTG basis states

To gain a better insight into the MTG basis states obtained from zero field
hWSs, let us consider an example of an isolated and topologically trivial
band, where 2D Wannierization can be achieved. The hWSs are related to 2D
Wannier orbitals via Fourier transform along the $\mathbf{a}_{2}$ axis, i.e.,

$\bra{\mathbf{r}}\ket{w(n_{0},k_{2})}=\frac{1}{\sqrt{N_{2}}}\sum_{\mathbf{R}}\phi_{\mathbf{R}}(\mathbf{r})\delta_{\mathbf{R}\cdot\mathbf{g}_{1},2\pi
n_{0}}e^{ik_{2}\mathbf{g}_{2}\cdot\mathbf{R}}.$ (36)

$\phi_{\mathbf{R}}(\mathbf{r})$ denotes a 2D localized Wannier orbital at site
$\mathbf{R}=m\mathbf{a}_{1}+n\mathbf{a}_{2}$, with the normalization condition
$\int\mathrm{d}^{2}\mathbf{r}\phi_{\mathbf{R}}^{*}(\mathbf{r})\phi_{\mathbf{R}^{\prime}}(\mathbf{r})=\delta_{\mathbf{R},\mathbf{R}^{\prime}}$.
The MTG basis states can then be written as

$\begin{split}\bra{\mathbf{r}}\ket{W_{r}(\mathbf{k})}=\frac{1}{\sqrt{N_{1}N_{2}}}\sum_{\mathbf{R}}e^{i\tilde{\mathbf{k}}\cdot\mathbf{R}}\hat{t}^{m}({\mathbf{a}_{1}})\hat{t}^{n}({\mathbf{a}_{2}})\phi_{\mathbf{0}}(\mathbf{r}),\end{split}$
(37)

where we have defined
$\tilde{\mathbf{k}}=k_{1}\mathbf{g}_{1}+(k_{2}+\frac{r}{q})\mathbf{g}_{2}$.
Compared to the zero-field Bloch states, the MTG states are obtained by acting
on the 2D Wannier orbital at site $\mathbf{R}=\mathbf{0}$ with non-commuting
magnetic translation operators instead of the usual lattice translation
operators.

We make a connection to the Peierl’s factor used in the literature [40], where
the phase $\exp\left(-i\frac{e}{\hbar
c}\int_{\mathbf{R}}^{\mathbf{r}}\mathrm{d}\mathbf{r}^{\prime}\cdot\mathbf{A}(\mathbf{r}^{\prime})\right)$
is attached to the 2D ($\mathbf{B}=0$) Wannier orbital
$\phi_{\mathbf{R}}(\mathbf{r})$ with the integration along the straight line
from $\mathbf{R}$ to $\mathbf{r}$. Expanding out the magnetic translation
operators in Eq. (37):

$\bra{\mathbf{r}}\ket{W_{r}(\mathbf{k})}=\frac{1}{\sqrt{N_{1}N_{2}}}\sum_{\mathbf{R}}e^{i\tilde{\mathbf{k}}\cdot\mathbf{R}}e^{i\frac{m(m-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{-im\mathbf{q}_{\phi}\cdot\mathbf{r}}\phi_{\mathbf{R}}(\mathbf{r}).$
(38)

If the Wannier orbitals are exponentially localized near the lattice sites,
and the magnetic length is much longer than their spatial extent, then we can
take $|\mathbf{r}-\mathbf{R}|/\ell\ll 1$ and approximate

$e^{-im\mathbf{q}_{\phi}\cdot\mathbf{r}}\phi_{\mathbf{R}}(\mathbf{r})\approx
e^{-im\mathbf{q}_{\phi}\cdot\mathbf{R}}\left[e^{-i\frac{e}{\hbar
c}\int_{\mathbf{R}}^{\mathbf{r}}\mathrm{d}\mathbf{r}^{\prime}\cdot\mathbf{A}(\mathbf{r}^{\prime})}\phi_{\mathbf{R}}(\mathbf{r})\right],$
(39)

where we used

$\begin{split}&i\frac{e}{\hbar
c}\int_{\mathbf{R}}^{\mathbf{r}}\mathrm{d}\mathbf{r}^{\prime}\cdot\mathbf{A}(\mathbf{r}^{\prime})\\\
=&\frac{i}{\ell^{2}}R_{x}(y-R_{y})+\frac{i}{2\ell^{2}}(x-R_{x})(y-R_{y})\\\
=&im\frac{a_{1x}}{\ell^{2}}(y-R_{y})+\frac{i}{2\ell^{2}}(x-R_{x})(y-R_{y})\\\
=&im\mathbf{q}_{\phi}\cdot(\mathbf{r}-\mathbf{R})+O(|\mathbf{r}-\mathbf{R}|^{2}/\ell^{2}).\end{split}$
(40)

This explicitly demonstrates the link between the MTG basis states and the
Wannier orbitals with a Peierl’s factor.

###  III.4 Matrix elements of the Hamiltonian

Here we work out the matrix elements of the single-electron Hamiltonian in the
MTG basis states defined in Eq. (28). Matrix elements with respect to the
orthonormal basis $\ket{V_{a}(\mathbf{k})}$ can be obtained via the basis
transformation in Eq. (33). Specifically:

$\begin{split}&\bra{W_{\alpha,r_{1}}(\mathbf{k})}\hat{H}\ket{W_{\beta,r_{2}}(\mathbf{k})}=\sum_{s}e^{i2\pi
k_{1}s}\\\
\times&e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}\bra{w_{\alpha}(0,\tilde{k}_{2})}\hat{H}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{\beta}(s,\tilde{p}_{2})},\end{split}$
(41)

where we defined $\tilde{k}_{2}=k_{2}+\frac{r_{1}}{q}$ and
$\tilde{p}_{2}=k_{2}+\frac{r_{2}}{q}$ ($\tilde{k}_{2},\tilde{p}_{2}\in[0,1)$).
This expression is obtained by expanding the MTG basis states using Eq. (28),
and making use of the operator identity (see appendix A):

$\hat{t}^{s}(\mathbf{a}_{1})=e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}^{s}(\mathbf{a}_{1}).$
(42)

We hereby split the full single-electron Hamiltonian into the zero-field
component $\hat{H}_{0}$ and the finite field component
$\hat{H}_{\mathbf{B}}=\hat{H}-\hat{H}_{0}$ respectively, and calculate each
term separately.

####  III.4.1 $\hat{H}_{0}$

For the zero field Hamiltonian, its matrix elements can be straightforwardly
studied as follows:

$\begin{split}I_{1}=&\sum_{s,s_{0},\gamma}e^{i2\pi
k_{1}s}e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}\\\
&\times\bra{w_{\alpha}(0,\tilde{k}_{2})}\hat{H}_{0}\ket{w_{\gamma}(s_{0},\tilde{k}_{2})}\\\
&\times\bra{w_{\gamma}(s_{0},\tilde{k}_{2})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{\beta}(s,\tilde{p}_{2})},\end{split}$
(43)

where we have inserted the narrow band projector
$\sum_{s_{0}\gamma}\ket{w_{\gamma}(s_{0},\tilde{k}_{2})}\bra{w_{\gamma}(s_{0},\tilde{k}_{2})}$.
The second line can be calculated by writing down hWSs in the Bloch band basis
using Eq. (23):

$\begin{split}&\bra{w_{\alpha}(0,\tilde{k}_{2})}\hat{H}_{0}\ket{w_{\gamma}(s_{0},\tilde{k}_{2})}\\\
=&\frac{1}{N_{1}}\sum_{\tilde{k}_{1},n}e^{-i2\pi\tilde{k}_{1}s_{0}}\left[U^{\dagger}_{\alpha,n}(\tilde{\mathbf{k}})\varepsilon_{n,\tilde{\mathbf{k}}}U_{n,\gamma}(\tilde{\mathbf{k}})\right],\end{split}$
(44)

where in this section we have redefined
$\tilde{\mathbf{k}}=\tilde{k}_{1}\mathbf{g}_{1}+\tilde{k}_{2}\mathbf{g}_{2}$
($\tilde{k}_{1},\tilde{k}_{2}\in[0,1)$). It leads to 1D hopping between hWSs
at different localization centers with the same wavevector $\tilde{k}_{2}$.
The last line in Eq. (43) can be calculated as:

$\begin{split}&\bra{w_{\gamma}(s_{0},\tilde{k}_{2})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{\beta}(s,\tilde{p}_{2})}\\\
=&\frac{1}{N_{1}}\sum_{\tilde{k}_{1},\tilde{p}_{1}}e^{i2\pi\tilde{k}_{1}s_{0}}e^{-i2\pi\tilde{p}_{1}s}\\\
&\times\sum_{n,m}U^{\dagger}_{\gamma,n}(\tilde{\mathbf{k}})\bra{\psi_{n}(\tilde{\mathbf{k}})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{\psi_{m}(\tilde{\mathbf{p}})}U_{m,\beta}(\tilde{\mathbf{p}}).\end{split}$
(45)

It represents scattering between Bloch eigenstates by wavevector
$s\mathbf{q}_{\phi}$. Therefore, the matrix elements of the zero-field
Hamiltonian in the MTG basis represent a hop along the $\mathbf{a}_{1}$
direction, followed by a magnetic scattering by wavevector
$s\mathbf{q}_{\phi}$. Note that if $\mathbf{a}_{1}\cdot\mathbf{a}_{2}\neq 0$,
i.e., non-orthogonal unit cell vectors such as triangular or honeycomb
lattice, magnetic scattering hybridizes Bloch states with wavevectors
$\tilde{\mathbf{k}}=\tilde{k}_{1}\mathbf{g}_{1}+\tilde{k}_{2}\mathbf{g}_{2}$
and
$\tilde{\mathbf{p}}=\tilde{p}_{1}\mathbf{g}_{1}+\tilde{p}_{2}\mathbf{g}_{2}$
such that:

$\tilde{p}_{1}=[\tilde{k}_{1}+\frac{sp}{q}\frac{a_{1y}}{a_{2}}]_{1},\
r_{2}=[r_{1}+sp]_{q}.$ (46)

Figure 4: Hofstadter spectra obtained using exact LL approach (blue) and the
hWS approach (red). Upper panels are square lattice and lower panels are
triangular lattice. A lattice potential strength $V_{0}=3$ is used. We set the
energy scale $\hbar^{2}/m_{e}a^{2}=1$. For LL approach we studied the $p/q$
sequence for $q=48$ and $p=1,\dots 48$. For hWS approach we studied the
sequence $q=96$ and $p=1,\dots 96$. The numbers in figures (a,b,d,e) label the
Chern number of the corresponding energy gaps.

####  III.4.2 $\hat{H}_{\mathbf{B}}$

The finite field term in the Landau gauge is given by
$\hat{H}_{\mathbf{B}}(\mathbf{r})=\frac{\hbar
x\hat{p}_{y}}{m_{e}\ell^{2}}+\frac{\hbar^{2}x^{2}}{2m_{e}\ell^{4}}$. It
contains polynomials of the coordinates, and grows when moving away from the
axes origin where $x=0$. Below we work out the matrix elements for the term
$\frac{\hbar x\hat{p}_{y}}{m_{e}\ell^{2}}$, and the second term can be
calculated in an analogous fashion. Specifically,

$\begin{split}I_{2}=&\sum_{s}e^{i2\pi
k_{1}s}e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}\\\
\times&\frac{\hbar^{2}}{m_{e}\ell^{2}}\bra{w_{\alpha}(0,\tilde{k}_{2})}\frac{x\hat{p}_{y}}{\hbar}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{\beta}(s,\tilde{p}_{2})}.\end{split}$
(47)

We first note that since the hWSs are exponentially suppressed away from their
localization centers, for $s\gg 0$, the term is exponentially small regardless
of the operator. Using the estimate $\langle x\rangle\sim a_{1x}$ and
$\langle\hat{p}_{y}\rangle\sim 2\pi\hbar/a_{2}$,
$I_{2}\sim\frac{\hbar^{2}}{m_{e}\ell^{2}}$, i.e., cyclotron frequency. At low
fields such that $p/q\ll 1$, the finite field term is smaller than the zero
field term by ratio $p/q$. Therefore, unless there are further band flattening
effects for the zero-field term (e.g., TBG at magic angle), the finite field
term is negligible.

The second line of Eq. (47) can be calculated by expanding in the Bloch basis:

$\begin{split}&\frac{\hbar}{m_{e}\ell^{2}N_{1}}\sum_{\tilde{k}_{1}\tilde{p}_{1}}e^{-i2\pi\tilde{p}_{1}s}e^{i\frac{m(m-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{r}}\\\
\times&\sum_{n,m}U^{\dagger}_{\alpha,n}(\tilde{\mathbf{k}})U_{m,\beta}(\tilde{\mathbf{p}})\bra{\psi_{n}(\tilde{\mathbf{k}})}\frac{x\hat{p}_{y}}{\hbar}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{\psi_{m}(\tilde{\mathbf{p}})},\end{split}$
(48)

where second line can be explicitly calculated in the plane wave basis:

$\begin{split}&\bra{\psi_{n}(\tilde{\mathbf{k}})}\frac{x\hat{p}_{y}}{\hbar}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{\psi_{m}(\tilde{\mathbf{p}})}\\\
=&\sum_{\mathbf{g},\mathbf{g}^{\prime}}u_{n\mathbf{g}}^{*}(\tilde{\mathbf{k}})u_{m\mathbf{g}^{\prime}}(\tilde{\mathbf{p}})\frac{1}{\mathcal{A}}\int\mathrm{d}^{2}\mathbf{r}\\\
&e^{-i(\mathbf{g}+\tilde{\mathbf{k}})\cdot\mathbf{r}}\left[x\left(g_{y}+\tilde{k}_{y}\right)\right]e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}e^{i(\mathbf{g}^{\prime}+\tilde{\mathbf{p}})\cdot\mathbf{r}}.\end{split}$
(49)

Here we have expanded the periodic part of the Bloch states in a Fourier
series:
$u_{n\tilde{\mathbf{k}}}(\mathbf{r})=\sum_{\mathbf{g}}e^{i\mathbf{g}\cdot\mathbf{r}}u_{n\mathbf{g}}(\tilde{\mathbf{k}})$,
and $\mathcal{A}=N_{1}N_{2}a_{1x}a_{2y}$ is the area of the system. It is
important that the real space integral $\int\mathrm{d}^{2}\mathbf{r}[...]$
must be placed in a box to avoid the revivals of the hWSs along the
$\mathbf{e}_{x}$ direction at the boundary of a torus. In practice we choose
the integration domain as
$x\in[-\frac{N_{1}a_{1x}}{2},-\frac{N_{1}a_{1x}}{2})$ (see appendix B for
details).

###  III.5 Comparing Hofstadter spectra

In Fig. 4 we show the comparison of the Hofstadter spectra for the lowest
energy Bloch band with square [Figs. 4(a,b,c)] and triangular [Figs. 4(d,e,f)]
lattice potentials respectively, calculated using the exact LL approach and
the approximate hWS approach. For the parameters used in the calculation, the
hWSs are well localized within a unit cell (see Fig. 1), and there is a good
quantitative agreement of the Hofstadter spectra up to $\phi/\phi_{0}\approx
0.2$. At higher magnetic flux, the MTG states generated from hWSs have a
larger spillover onto remote magnetic subbands, and lead to an overall upward
shift of the Hofstadter spectra when compared to LL calculations. Despite the
quantitative differences in the Hofstadter spectra, the prominent magnetic
subbands and the Chern numbers associated with gaps are correctly captured via
the hWS method all the way to $\phi/\phi_{0}=1$, as indicated by the labels in
Figs. 4(a,b,d,e).

###  III.6 Harper equation

Eq. (41) gives a general procedure for calculating the matrix elements of the
single-electron Hamiltonian for a subset of Bloch bands isolated from the
rest, by making use of the hWSs at zero field and projecting onto
representations of the MTG at finite field. Here we elaborate on the case of a
square lattice potential and an isolated band with trivial topology. As
mentioned, in this case the hWSs are just the 1D Fourier transforms of the
Bloch eigenstates with a smooth gauge in the Brillouin zone and the relation
to 2D Wannier orbitals is described in Eq. (36). Below we show that the Harper
equation of Hofstadter [42] is recovered as long as the magnetic length $\ell$
is much longer than the localization length $\xi$ of the Wannier orbitals. We
express the left hand side of Eq. (45) using 2D exponentially localized
Wannier orbitals,

$\begin{split}&\bra{w(s_{0},\tilde{k}_{2})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w(s,\tilde{p}_{2})}\\\
=&\frac{1}{N_{2}}\sum_{n_{1},n_{2}}\left[\int\mathrm{d}^{2}\mathbf{r}\phi_{\mathbf{R}_{1}}^{*}(\mathbf{r})\phi_{\mathbf{R}_{2}}(\mathbf{r})e^{-is\mathbf{q}_{\phi}\cdot(\mathbf{r}-\mathbf{R}_{2})}\right]e^{-is\mathbf{q}_{\phi}\cdot\mathbf{R}_{2}}e^{-i2\pi\tilde{k}_{2}n_{1}}e^{i2\pi\tilde{p}_{2}n_{2}},\end{split}$
(50)

where we defined $\mathbf{R}_{1}=s_{0}\mathbf{a}_{1}+n_{1}\mathbf{a}_{2}$ and
$\mathbf{R}_{2}=s\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}$. Since Wannier orbitals
are exponentially localized, $s_{0}$ and $s_{1}$ have to be close to each
other. Making the assumption that the 1D hopping is short ranged, Eq. (44)
forces $s_{0}$ to be near 0, therefore $s$ is also close to 0. As a result,
the factor $e^{-is\mathbf{q}_{\phi}\cdot(\mathbf{r}-\mathbf{R}_{2})}$ is
slowly varying over the lengthscale $\sim\xi$, the spatial extent of the
exponentially localized Wannier orbital, and can be Taylor expanded in the
real space integration in the bracket. The integral can therefore be expanded
in power series of $\xi/\ell$ as,

$\begin{split}&\int\mathrm{d}^{2}\mathbf{r}\phi_{\mathbf{R}_{1}}^{*}(\mathbf{r})\phi_{\mathbf{R}_{2}}(\mathbf{r})e^{-is\mathbf{q}_{\phi}\cdot(\mathbf{r}-\mathbf{R}_{2})}\\\
\approx&\int\mathrm{d}^{2}\mathbf{r}\phi_{\mathbf{R}_{1}}^{*}(\mathbf{r})\phi_{\mathbf{R}_{2}}(\mathbf{r})\left[1-is\mathbf{q}_{\phi}\cdot(\mathbf{r}-\mathbf{R}_{2})+\ldots\right]\\\
=&\delta_{s_{0},s}\delta_{n_{1},n_{2}}+O(a_{1x}\xi/\ell^{2}).\end{split}$ (51)

As a result, Eq. (50) becomes:

$\begin{split}&\bra{w(s_{0},\tilde{k}_{2})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w(s,\tilde{p}_{2})}\\\
=&\delta_{s_{0},s}\delta_{[r_{1}+sp]_{q},r_{2}}+O(a_{1x}\xi/\ell^{2}),\end{split}$
(52)

where in the last line we have substituted $\tilde{k}_{2}=k_{2}+r_{1}/q$ and
$\tilde{p}_{2}=k_{2}+r_{2}/q$. Keeping leading order term in the above
expression, and combining with Eq. (44) we get the following simplification
for Eq. (43):

$\begin{split}I_{1}\approx&\frac{1}{N_{1}}\sum_{s}e^{i2\pi
k_{1}s}\sum_{\tilde{k}_{1}}e^{-i2\pi\tilde{k}_{1}s}\varepsilon_{\tilde{\mathbf{k}}}\delta_{[r_{1}+sp]_{q},r_{2}}.\end{split}$
(53)

For well localized Wannier states the dispersion can be approximated by
nearest-neighbor hoppings, i.e.,
$\varepsilon_{\tilde{\mathbf{k}}}=2E_{0}[\cos(2\pi\tilde{k}_{1})+\cos(2\pi\tilde{k}_{2})]$,
and as a result:

$I_{1}\approx
E_{0}\left[2\cos(2\pi\tilde{k}_{2})\delta_{r_{1},r_{2}}+\sum_{s=\pm 1}e^{i2\pi
k_{1}s}\delta_{[r_{1}+sp]_{q},r_{2}}\right].$ (54)

We compare the above expression to the Harper equation obtained in Ref. [42]:

$g(m+1)+g(m-1)+2\cos(\frac{2\pi p}{q}m-2\pi k_{1})g(m)=\frac{E}{E_{0}}g(m),$
(55)

where we have replaced the notation $\nu$ used in Ref. [42] with $2\pi k_{1}$,
and $\alpha$ with $p/q$. Performing Fourier series expansion of the
wavefunction
$g(m)=\sum_{k_{2},r_{2}}e^{i2\pi(k_{2}+\frac{r_{2}}{q})m}g_{k_{2},r_{2}}$
produces the eigen-equation:

$\begin{split}&Eg_{k_{2},r_{1}}\\\
=&\sum_{r_{2}}E_{0}\left[2\cos(2\pi\tilde{k}_{2})\delta_{r_{1},r_{2}}+\sum_{s=\pm
1}e^{i2\pi
k_{1}s}\delta_{[r_{1}+sp]_{q},r_{2}}\right]g_{k_{2},r_{2}}.\end{split}$ (56)

The second line is precisely Eq. (54), explicitly recovering the Harper
equation as a limiting case of the finite magnetic field problem when the
magnetic length is much longer than the spatial extent of the Wannier orbital.

##  IV Summary

Using the magnetic translation group projection method recently developed in
the Ref. [46] based on the $\mathbf{B}=0$ hybrid Wannier states (hWS), we re-
examined the Hofstadter physics of the $\mathbf{B}\neq 0$ magnetic subbands
corresponding to a subset of energetically isolated, and topologiclaly
trivial, Bloch bands. Employing the continuum electron moving in 2D square and
triangular lattice potentials as examples, we demonstrated that the method
works well up to moderate strengths of magnetic flux per unit cell.
Importantly, it naturally bridges the zero field and finite field Hilbert
spaces. We also recovered the Harper equation [42] from the MTG states
generated from hWSs and showed that its regime of validity is determined by
the ratio of magnetic length to the size of exponentially localized Wannier
orbitals.

Although here we only applied this projection method to the non-interacting
Hamiltonians, the hWS procedure was shown to be useful in studying interacting
Hofstadter problem [46]. This is primarily due to two reasons: (1) At low
magnetic fields, the procedure can faithfully construct the correct Hilbert
space without having to increase the internal dimensions unlike the
conventional LL-based approach, where the upper LL index cutoff increases with
decreasing field to achieve numerical convergence. (2) The plane wave nature
of the hWSs makes the matrix representation of operators of the form
$\mathcal{O}_{\mathbf{q}}e^{i\mathbf{q}\cdot\mathbf{r}}$ sparse. By comparison
they are dense matrices in the LL wavefunction basis. The method can therefore
efficiently address the effects of Coulomb interactions in a finite magnetic
field in the presence of a periodic potential.

##  V Acknowledgements

X.W. acknowledges financial support from the Gordon and Betty Moore
Foundation’s EPiQS Initiative Grant GBMF11070. 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, 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).
  * Cao  _et al._ [2020] Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 583, 215 (2020).
  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Nature 579, 353 (2020).
  * Liu  _et al._ [2020] X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. Haei Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Nature 583, 221 (2020).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Nature Materials 19, 861 (2020).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang, Nature 579, 359 (2020).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021).
  * Hao  _et al._ [2021] Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).
  * Kometter  _et al._ [2022] C. R. Kometter, J. Yu, T. Devakul, A. P. Reddy, Y. Zhang, B. A. Foutty, K. Watanabe, T. Taniguchi, L. Fu, and B. E. Feldman, Hofstadter states and reentrant charge order in a semiconductor moiré lattice (2022).
  * Yuan and Fu [2018] N. F. Q. Yuan and L. Fu, Phys. Rev. B 98, 045103 (2018).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * Kang and Vafek [2018] J. Kang and O. Vafek, Phys. Rev. X 8, 031088 (2018).
  * Po  _et al._ [2018a] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018a).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Phys. Rev. Lett. 122, 246401 (2019).
  * Bernevig  _et al._ [2021a] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205411 (2021a).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Phys. Rev. B 103, 205412 (2021).
  * Bernevig  _et al._ [2021b] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205413 (2021b).
  * Lian  _et al._ [2021a] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Phys. Rev. B 103, 205414 (2021a).
  * Bernevig  _et al._ [2021c] B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Phys. Rev. B 103, 205415 (2021c).
  * Balents  _et al._ [2020] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Nature Physics 16, 725 (2020).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Phys. Rev. Lett. 129, 047601 (2022).
  * 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).
  * 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).
  * 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).
  * 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).
  * 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).
  * 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).
  * 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).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Nature Materials 20, 488 (2021).
  * Finney  _et al._ [2022] 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. A. Kastner, and D. Goldhaber-Gordon, Proceedings of the National Academy of Sciences 119, e2118482119 (2022).
  * 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, and B. E. Feldman, Nature Physics 18, 825 (2022).
  * 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._ [2021b] B. Lian, F. Xie, and B. A. Bernevig, Phys. Rev. B 103, L161405 (2021b).
  * Herzog-Arbeitman  _et al._ [2022a] J. Herzog-Arbeitman, A. Chew, D. K. Efetov, and B. A. Bernevig, Phys. Rev. Lett. 129, 076401 (2022a).
  * Parker  _et al._ [2021] D. Parker, P. Ledwith, E. Khalaf, T. Soejima, J. Hauschild, Y. Xie, A. Pierce, M. P. Zaletel, A. Yacoby, and A. Vishwanath, arXiv:2112.13837 (2021).
  * Herzog-Arbeitman  _et al._ [2022b] J. Herzog-Arbeitman, A. Chew, and B. A. Bernevig, Phys. Rev. B 106, 085140 (2022b).
  * Peierls [1933] R. Peierls, Zeitschrift für Physik 80, 763 (1933).
  * Luttinger [1951] J. M. Luttinger, Phys. Rev. 84, 814 (1951).
  * Harper [1955] P. G. Harper, Proceedings of the Physical Society. Section A 68, 874 (1955).
  * Hofstadter [1976] D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * Soluyanov and Vanderbilt [2011] A. A. Soluyanov and D. Vanderbilt, Phys. Rev. B 83, 035108 (2011).
  * Po  _et al._ [2018b] H. C. Po, H. Watanabe, and A. Vishwanath, Phys. Rev. Lett. 121, 126402 (2018b).
  * Wang and Vafek [2020] X. Wang and O. Vafek, Phys. Rev. B 102, 075142 (2020).
  * Wang and Vafek [2022] X. Wang and O. Vafek, Phys. Rev. B 106, L121111 (2022).
  * Sgiarovello  _et al._ [2001] C. Sgiarovello, M. Peressi, and R. Resta, Phys. Rev. B 64, 115202 (2001).
  * Yu  _et al._ [2011] R. Yu, X. L. Qi, A. Bernevig, Z. Fang, and X. Dai, Phys. Rev. B 84, 075119 (2011).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Kang and Vafek [2020] J. Kang and O. Vafek, Phys. Rev. B 102, 035161 (2020).

##  Appendix A Landau gauge magnetic translation group identities

Here we derive an identity for $\hat{t}^{s}(\mathbf{a}_{1})$, which relates to
the zero field discrete translation operator $\hat{T}(\mathbf{a}_{1})$. Note
that:

$\begin{split}\hat{t}^{s}(\mathbf{a}_{1})&=\left[e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}(\mathbf{a}_{1})\right]^{s-2}\left[e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}(\mathbf{a}_{1})\right]\left[e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}(\mathbf{a}_{1})\right]\\\
&=\left[e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}(\mathbf{a}_{1})\right]^{s-2}e^{i\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{-i2\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}^{2}(\mathbf{a}_{1})\\\
&=\left[e^{-i\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}(\mathbf{a}_{1})\right]^{s-3}e^{i(1+2)\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{-i3\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}^{3}(\mathbf{a}_{1})\\\
&=e^{i[1+2+\dots+(s-1)]\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}^{s}(\mathbf{a}_{1})\\\
&=e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\hat{T}^{s}(\mathbf{a}_{1}).\end{split}$
(57)

We also show the MTG basis states, $\ket{\Psi_{n,r}(\mathbf{k})}$ expressed
using LL wavefunctions, form a complete set for $r=0,\dots p-1$. It is
sufficient to show that:

$\begin{split}\ket{\Psi_{n,r+p}(\mathbf{k})}&=e^{i2\pi\left(k_{1}-(k_{2}+\frac{r}{q})\frac{a_{1y}}{a_{2}}\right)}\ket{\Psi_{n,r}(\mathbf{k})},\end{split}$
(58)

namely that states at $r+p$ and $r$ differ by a complex phase factor. Note
that by definition in Eq. (10),

$\ket{\Psi_{n,r+p}(\mathbf{k})}=\frac{1}{\sqrt{N}}\sum_{s=-\infty}^{\infty}e^{i2\pi
k_{1}s}\hat{t}^{s}(\mathbf{a}_{1})\ket{\psi_{n}\left(\frac{2\pi}{a_{2}}(k_{2}+\frac{r}{q}+\frac{p}{q})\right)},$
(59)

and the LL wavefunction is expressed in Eq. (1) as:

$\ket{\psi_{n}\left(\tilde{k}_{y}+\frac{2\pi}{a_{2}}\frac{p}{q})\right)}=e^{i\tilde{k}_{y}y}e^{i\frac{2\pi}{a_{2}}\frac{p}{q}y}\hat{T}(-\frac{2\pi}{a_{2}}\frac{p}{q}\ell^{2}\mathbf{e}_{x})\hat{T}(-\tilde{k}_{y}\ell^{2}\mathbf{e}_{x})\ket{n}.$
(60)

Here we have defined
$\tilde{k}_{y}\equiv\frac{2\pi}{a_{2}}\left(k_{2}+\frac{r}{q}\right)$ for
notational convenience. Note that

$\frac{2\pi}{a_{2}}\frac{p}{q}\ell^{2}\mathbf{e}_{x}=\frac{1}{a_{2}}\cdot\frac{a_{1x}a_{2}}{\ell^{2}}\cdot\ell^{2}\mathbf{e}_{x}=a_{1x}\mathbf{e}_{x},$
(61)

and that

$\hat{T}(-a_{1x}\mathbf{e}_{x})\ket{n}=\hat{T}(-\mathbf{a}_{1})\ket{n}.$ (62)

Making use of the identity:

$\hat{t}(-\mathbf{a}_{1})\equiv\hat{t}^{\dagger}(\mathbf{a}_{1})=\hat{T}^{\dagger}(\mathbf{a}_{1})e^{i\mathbf{q}_{\phi}\cdot\mathbf{r}}=e^{i\frac{2\pi}{a_{2}}\frac{p}{q}a_{1y}}e^{i\frac{2\pi}{a_{2}}\frac{p}{q}y}\hat{T}(-\mathbf{a}_{1}),$
(63)

we can rewrite Eq. (60) as:

$\ket{\psi_{n}\left(\tilde{k}_{y}+\frac{2\pi}{a_{2}}\frac{p}{q})\right)}=e^{-i\frac{2\pi}{a_{2}}\frac{p}{q}a_{1y}}e^{i\tilde{k}_{y}y}\hat{t}(-\mathbf{a}_{1})\hat{T}(-\tilde{k}_{y}\ell^{2}\mathbf{e}_{x})\ket{n}=e^{-i\left(\tilde{k}_{y}+\frac{2\pi}{a_{2}}\frac{p}{q}\right)a_{1y}}\hat{t}(-\mathbf{a}_{1})\ket{\psi_{n}(\tilde{k}_{y})}.$
(64)

Here for the last equality we used
$e^{i\tilde{k}_{y}y}\hat{t}(-\mathbf{a}_{1})=e^{-i\tilde{k}_{y}a_{1y}}\hat{t}(-\mathbf{a}_{1})e^{i\tilde{k}_{y}y}$.
Substituting this back into Eq. (59),

$\begin{split}\ket{\Psi_{n,r+p}(\mathbf{k})}&=e^{-i\frac{2\pi}{a_{2}}\left(k_{2}+\frac{r+p}{q}\right)a_{1y}}\frac{1}{\sqrt{N}}\sum_{s=-\infty}^{\infty}e^{i2\pi
k_{1}s}\hat{t}^{s}(\mathbf{a}_{1})\hat{t}(-\mathbf{a}_{1})\ket{\psi_{n}\left(\frac{2\pi}{a_{2}}(k_{2}+\frac{r}{q})\right)}\\\
&=e^{i2\pi\left(k_{1}-(k_{2}+\frac{r+p}{q})\frac{a_{1y}}{a_{2}}\right)}\frac{1}{\sqrt{N}}\sum_{s=-\infty}^{\infty}e^{i2\pi
k_{1}(s-1)}\hat{t}^{s-1}(\mathbf{a}_{1})\ket{\psi_{n}\left(\frac{2\pi}{a_{2}}(k_{2}+\frac{r}{q})\right)}\\\
&=e^{i2\pi\left(k_{1}-(k_{2}+\frac{r+p}{q})\frac{a_{1y}}{a_{2}}\right)}\ket{\Psi_{n,r}(\mathbf{k})}.\end{split}$
(65)

##  Appendix B Overlap matrix

Here we show how to calculate the overlap matrix between MTG basis states
generated from zero field hybrid Wannier states, Eq. (32). Specifically, we
define the overlap matrix: $\Lambda_{\alpha r_{1},\beta r_{2}}(k_{1},k_{2})$
where $k_{1}\in[0,1)$, $k_{2}\in[0,1/q)$, and $r_{1,2}=0,\dots q-1$. It can be
calculated by switching to the Bloch basis,

$\begin{split}\Lambda_{\alpha r_{1},\beta
r_{2}}(k_{1},k_{2})&=\sum_{s}e^{i2\pi
k_{1}s}e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}\bra{w_{\alpha}(0,k_{2}+\frac{r_{1}}{q})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{w_{\beta}(s,k_{2}+\frac{r_{2}}{q})}\\\
&=\sum_{s}e^{i2\pi
k_{1}s}e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}\frac{1}{N_{1}}\sum_{\bar{k}_{1}\bar{p}_{1}}e^{-i2\pi\bar{p}_{1}s}\bra{\psi_{\alpha}(\bar{k}_{1},k_{2}+\frac{r_{1}}{q})}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}\ket{\psi_{\beta}(\bar{p}_{1},k_{2}+\frac{r_{2}}{q})}.\end{split}$
(66)

On the second line we introduced the Bloch states
$\ket{\psi_{\alpha}(\mathbf{k})}=\sum_{n}U_{n,\alpha}(\mathbf{k})\ket{\psi_{n}(\mathbf{k})}$,
with the distinction that the subscript $\alpha$ labels Fourier transforms of
the hWS, whereas $n$ labels Bloch eigenstates. Second line is calculated in
the plain wave basis, by writing down
$\bra{\mathbf{r}}\ket{\psi_{\alpha}(\mathbf{k})}=\sum_{\mathbf{g}}e^{i(\mathbf{k}+\mathbf{g})\cdot\mathbf{r}}u_{\mathbf{g}}(\alpha\mathbf{k})$.
We arrive at the following expression:

$\begin{split}\Lambda_{\alpha r_{1},\beta
r_{2}}(k_{1},k_{2})&=\sum_{s}e^{i2\pi
k_{1}s}e^{i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}\frac{1}{N_{1}}\sum_{\bar{k}_{1}\bar{p}_{1}}\sum_{\mathbf{g},\mathbf{g}^{\prime}}e^{-i2\pi\bar{p}_{1}s}u^{*}_{\mathbf{g}}(\alpha,\bar{k}_{1},k_{2}+\frac{r_{1}}{q})u_{\mathbf{g}^{\prime}}(\beta,\bar{p}_{1},k_{2}+\frac{r_{2}}{q})\\\
&\times\frac{1}{\mathcal{A}}\int\mathrm{d}^{2}\mathbf{r}e^{-i(\bar{\mathbf{k}}+\mathbf{g})\cdot\mathbf{r}}e^{-is\mathbf{q}_{\phi}\cdot\mathbf{r}}e^{i(\bar{\mathbf{p}}+\mathbf{g}^{\prime})\cdot\mathbf{r}}.\end{split}$
(67)

The second line gives the constraint
$\bar{\mathbf{k}}+\mathbf{g}+s\mathbf{q}_{\phi}=\bar{\mathbf{p}}+\mathbf{g}^{\prime}$,
which is equivalent to:

$\bar{k}_{1}+l_{1}+\frac{sp}{q}\frac{a_{1y}}{a_{2}}=\bar{p}_{1}+l_{1}^{\prime},\
\frac{r_{1}}{q}+l_{2}+\frac{sp}{q}=\frac{r_{2}}{q}+l_{2}^{\prime}.$ (68)

##  Appendix C Real space integration of $x$ and $x^{2}$

The single-electron Hamiltonian in a finite magnetic field contains terms such
as $x\hat{p}_{y}$ and $x^{2}$, we hereby give analytical expressions for their
real space integrations. In evaluating the matrix elements for the
Hamiltonian, boundary terms do not matter, as the exponential localization of
the hWSs guarantees that the main contributing terms to $\sum_{s}$ in Eq. (41)
are restricted to the vicinity of $s=0$. We therefore choose the real space
domain to contain $N_{1}$ and $N_{2}$ unit cells in the $\mathbf{a}_{1}$ and
$\mathbf{a}_{2}$ directions, and place the origin of the coordinate system to
be at the center of the real space domain. Specifically the integral is
performed in the following manner:

$\frac{1}{\mathcal{A}}\int\mathrm{d}^{2}\mathbf{r}[...]=\frac{1}{N_{1}N_{2}}\frac{1}{a_{1x}a_{2}}\int_{-\frac{N_{1}}{2}a_{1x}}^{\frac{N_{1}}{2}a_{1x}}\mathrm{d}x\int_{\frac{a_{1y}}{a_{1x}}x-\frac{N_{2}}{2}a_{2}}^{\frac{a_{1y}}{a_{1x}}x+\frac{N_{2}}{2}a_{2}}\mathrm{d}y[...].$
(69)

As a result,

$\frac{1}{\mathcal{A}}\int\mathrm{d}^{2}\mathbf{r}\left(\frac{x}{a_{1x}}\right)e^{i\mathbf{q}\cdot\mathbf{r}}=\delta_{q_{y},0}(-i)\frac{\tilde{q}_{x}\cos(\frac{N_{1}}{2}\tilde{q}_{x})-\frac{2}{N_{1}}\sin(\frac{N_{1}}{2}\tilde{q}_{x})}{\tilde{q}_{x}^{2}}.$
(70)

and

$\frac{1}{\mathcal{A}}\int\mathrm{d}^{2}\mathbf{r}\left(\frac{x}{a_{1x}}\right)^{2}e^{i\mathbf{q}\cdot\mathbf{r}}=\delta_{q_{y},0}\frac{4\tilde{q}_{x}\cos(\frac{N_{1}}{2}\tilde{q}_{x})+(-\frac{8}{N_{1}}+N_{1}\tilde{q}_{x})\sin(\frac{N_{1}}{2}\tilde{q}_{x})}{2\tilde{q}_{x}^{3}}.$
(71)

where we have defined $\tilde{q}_{x}\equiv q_{x}a_{1x}$.

##  Appendix D Overlap between MTG basis states generated from LL approach and
hWS approach

To demonstrate the accuracy of the hWS approach we need to calculate its
overlap with the exact wavefunctions which can be obtained using LL approach,
as discussed in Eq. (34). Here we calculate the following expression:

$\Lambda_{\Psi,W}(nr_{1},\alpha
r_{2};\mathbf{k})\equiv\bra{\Psi_{n,r_{1}}(\mathbf{k})}\ket{W_{\alpha,r_{2}}(\mathbf{k})},$
(72)

where $n$ runs over LL indices, $\alpha$ runs over hWS indices,
$r_{1}=0\dots,p-1$, and $r_{2}=0,\dots q-1$. Using the definitions Eqs. (10)
and (28), we have:

$\begin{split}\Lambda_{\Psi,W}(nr_{1},\alpha
r_{2};\mathbf{k})&=\sum_{s}e^{-i2\pi
k_{1}s}\bra{\psi_{n}(k_{2}+\frac{r_{1}}{q})}\hat{t}_{\mathbf{a}_{1}}^{-s}\ket{w_{\alpha}(0,k_{2}+\frac{r_{2}}{q})}\\\
&=\sum_{s}e^{-i2\pi
k_{1}s}e^{-i\frac{s(s-1)}{2}\mathbf{q}_{\phi}\cdot\mathbf{a}_{1}}e^{i2\pi(k_{2}+\frac{r_{1}}{q})\frac{a_{1y}}{a_{2}}s}\frac{1}{\sqrt{N_{1}}}\sum_{\bar{k}_{1}}\bra{\psi_{n}(k_{2}+\frac{r_{1}}{q}-\frac{sp}{q})}\ket{\psi_{\alpha}(\bar{k}_{1},k_{2}+\frac{r_{2}}{q})}.\end{split}$
(73)

We evaluate the overlap between a LL wavefunction with a Bloch wavefunction as
follows:

$\begin{split}&\bra{\psi_{n}(k_{2}+\frac{r_{1}}{q}-\frac{sp}{q})}\ket{\psi_{\alpha}(\bar{k}_{1},k_{2}+\frac{r_{2}}{q})}\\\
=&\frac{1}{N_{2}a_{2}}\frac{1}{\sqrt{N_{1}a_{1x}\ell}}\sum_{\mathbf{g}}u_{\alpha\mathbf{g}}\int\mathrm{d}^{2}\mathbf{r}\phi_{n}\left(\frac{x}{\ell}+\frac{2\pi}{a_{2}}(k_{2}+\frac{r_{1}-sp}{q})\ell\right)e^{-i\frac{2\pi}{a_{2}}(k_{2}+\frac{r_{1}-sp}{q})y}e^{i(\bar{\mathbf{k}}+\mathbf{g})\cdot\mathbf{r}}\\\
=&\frac{1}{\sqrt{N_{1}a_{1x}\ell}}\sum_{\mathbf{g}}u_{\alpha\mathbf{g}}\delta_{\frac{r_{1}-sp}{q},\frac{r_{2}}{q}+l_{2}}\int\mathrm{d}x\phi_{n}\left(\frac{x}{\ell}+(\bar{k}_{y}+g_{y})\ell\right)e^{i(\bar{k}_{x}+g_{x})x}\\\
=&\frac{1}{\sqrt{N_{1}a_{1x}\ell}}\sum_{\mathbf{g}}u_{\alpha\mathbf{g}}\delta_{\frac{r_{1}-sp}{q},\frac{r_{2}}{q}+l_{2}}e^{-i(\bar{k}_{x}+g_{x})(\bar{k}_{y}+g_{y})\ell^{2}}\int\mathrm{d}x\phi_{n}\left(\frac{x}{\ell}\right)e^{i(\bar{k}_{x}+g_{x})x}\\\
=&\frac{1}{\sqrt{N_{1}}}\sqrt{\frac{\ell}{a_{1x}}}\sum_{\mathbf{g}}u_{\alpha\mathbf{g}}\delta_{\frac{r_{1}-sp}{q},\frac{r_{2}}{q}+l_{2}}e^{-i(\bar{k}_{x}+g_{x})(\bar{k}_{y}+g_{y})\ell^{2}}\sqrt{2\pi}(i)^{n}\phi_{n}\left((\bar{k}_{x}+g_{x})\ell\right)\end{split}$
(74)

Here we have defined
$\bar{\mathbf{k}}=\bar{k}_{1}\mathbf{g}_{1}+(k_{2}+\frac{r_{2}}{q})\mathbf{g}_{2}\equiv\bar{k}_{x}\mathbf{e}_{x}+\bar{k}_{y}\mathbf{e}_{y}$.
On the second line we used the following definitions for the Bloch and LL
wavefunctions:

$\displaystyle\bra{\mathbf{r}}\ket{\psi_{\alpha}(\mathbf{k})}$
$\displaystyle=\frac{1}{\sqrt{N_{2}a_{2}}}\frac{1}{\sqrt{N_{1}a_{1x}}}\sum_{\mathbf{g}}u_{\alpha\mathbf{g}}(\mathbf{k})e^{i(\mathbf{k}+\mathbf{g})\cdot\mathbf{r}},$
(75) $\displaystyle\bra{\mathbf{r}}\ket{\psi_{n}(k_{2})}$
$\displaystyle=\frac{1}{\sqrt{N_{2}a_{2}}}\frac{1}{\sqrt{\ell}}\phi_{n}\left(\frac{x}{\ell}+\frac{2\pi}{a_{2}}k_{2}\ell\right)e^{i\frac{2\pi}{a_{2}}k_{2}y},$
(76)

where:

$\phi_{n}\left(\frac{x}{\ell}\right)=\frac{1}{\pi^{1/4}\sqrt{2^{n}n!}}e^{-\frac{x^{2}}{2\ell^{2}}}H_{n}\left(\frac{x}{\ell}\right).$
(78)