††thanks: These authors contributed equally††thanks: These authors contributed
equally††thanks: These authors contributed equally

# Supplementary Information for:
Hofstadter states and reentrant charge order in a semiconductor moiré lattice

Carlos R. Kometter  Geballe Laboratory for Advanced Materials, Stanford, CA
94305, USA  Department of Physics, Stanford University, Stanford, CA 94305,
USA  Jiachen Yu  Geballe Laboratory for Advanced Materials, Stanford, CA
94305, USA  Department of Applied Physics, Stanford University, Stanford, CA
94305, USA  Trithep Devakul  Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA  Aidan P. Reddy  Department of
Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA  Yang Zhang  Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA  Benjamin A. Foutty  Geballe
Laboratory for Advanced Materials, Stanford, CA 94305, USA  Department of
Physics, Stanford University, Stanford, CA 94305, USA  Kenji Watanabe
Research Center for Functional Materials, National Institute for Material
Science, 1-1 Namiki, Tsukuba 305-0044, Japan  Takashi Taniguchi  International
Center for Materials Nanoarchitectonics, National Institute for Material
Science, 1-1 Namiki, Tsukuba 305-0044, Japan  Liang Fu  Department of Physics,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Benjamin E. Feldman<EMAIL_ADDRESS>Geballe Laboratory for Advanced
Materials, Stanford, CA 94305, USA  Department of Physics, Stanford
University, Stanford, CA 94305, USA  Stanford Institute for Materials and
Energy Sciences, SLAC National Accelerator Laboratory, Menlo Park, CA 94025,
USA

This PDF file includes:

Materials and Methods

Supplementary Text

Supplementary Figs. S1 to S13

###### Contents

  1. I Materials and Methods
    1. I.1 Sample fabrication
    2. I.2 Single-electron transistor (SET) fabrication and measurement
    3. I.3 Determination of stacking configuration
    4. I.4 Determination of the twist angle and filling factors
    5. I.5 Gap extraction
  2. II Comparison between transport and compressibility measurements
  3. III Spatial dependence
  4. IV Density function theory (DFT) and Hartree-Fock (HF) calculations
  5. V Phenomenological Stoner model
  6. VI Hubbard model
  7. VII Singlet-triplet quantum dot transition
  8. VIII Single-particle Hofstadter spectrum
  9. IX Gaps of the Hofstadter states

##  I Materials and Methods

###  I.1 Sample fabrication

Figure S1: Device image. Optical microscope image of the MoSe2/WSe2 device.
Platinum contacts to the sample, gold ‘contact’ gates to locally dope the
heterobilayer, and the back gate electrode are all indicated. Scale bar: 10
$\mu$m.

The twisted MoSe2/WSe2 device (Fig. S1) was fabricated using standard dry
transfer techniques Wang  _et al._ (2013). The MoSe2 and WSe2 monolayers, few-
layer graphite, and hexagonal boron nitride (hBN) flakes were exfoliated from
bulk crystals. All pickup steps were performed with poly(bisphenol A
carbonate) (PC) / polydimethylsiloxane (PDMS) stamps. First, pre-patterned Pt
contacts (5 nm thickness) were evaporated using standard electron beam
(e-beam) lithography and metallization techniques onto a 30-nm hBN previously
stacked on few-layer graphite. Annealing was performed at 400$\degree$C for 8
hours before and after pre-patterning to eliminate any polymer residue. A
30-nm hBN flake was then used to pick up monolayers of MoSe2 and WSe2, which
were deposited atop the pre-patterned contacts to complete the stack. Finally,
e-beam lithography, CHF3/O2 etching, and metal deposition were used to pattern
local ‘contact’ gates above the Pt contact area, and to make electrical
contact to the Pt leads and graphite gate with Cu/Au (2 nm/100 nm). The
‘contact’ gates served to locally dope the transition metal dichalcogenide
(TMD) to achieve Ohmic contact Movva  _et al._ (2015).

###  I.2 Single-electron transistor (SET) fabrication and measurement

The single-electron transistor (SET) sensor was fabricated by evaporating
aluminium onto the apex of a pulled quartz rod. The size of the apex, and thus
the lateral dimension of the SET, is estimated to be 50-100 nm based on
scanning electron microscope imaging of tips fabricated using the same
procedure. It was brought approximately 50 nm above the sample surface,
resulting in an overall spatial resolution of about 100-150 nm. The scanning
SET measurements were performed in a Unisoku USM 1300 scanning probe
microscope with a customized microscope head. An AC excitation of
$V_{\textrm{G,AC}}=4-8$ mV at frequency $f_{\textrm{G}}=233.33$ Hz was applied
to the back gate, and an AC excitation $V_{\textrm{2D,AC}}=5$ mV at frequency
$f_{\textrm{2D}}=167.77$ Hz was applied to the sample. We then measured the
inverse compressibility
$\text{d}\mu/\text{d}n\propto\frac{I_{\textrm{G}}/V_{\textrm{G,AC}}}{I_{\textrm{2D}}/V_{\textrm{2D,AC}}}$,
where $I_{\textrm{G}}$ and $I_{\textrm{2D}}$ are demodulated from the SET
current through the SET probe using standard lock-in techniques. A DC offset
voltage $V_{\textrm{2D,DC}}$ is further applied to the sample to maintain the
SET at its maximum sensitivity point within a Coulomb blockade oscillation
fringe chosen to be near the ‘flat-band’ condition where the tip does not gate
the sample. This minimizes tip-induced doping and provides an independent
direct DC measurement of $\mu(n)$. If the sample or contacts are highly
resistive such that the inverse measurement frequency is short compared to the
RC charging time of the sample, the AC measurement of inverse compressibility
can be artificially enhanced. The DC measurement is noisier but less
susceptible to these effects. For the data presented in this manuscript, the
AC and DC measurements agree within the noise, indicating sufficient charging
even in the AC scheme. We therefore obtain $\mu(n)$ by integrating the lower-
noise AC measurement of $\text{d}\mu/\text{d}n$. To reduce contact resistance
and facilitate sample charging, a large negative voltage $V_{\textrm{CG}}$ is
applied to the ‘contact’ gates. All measurements were performed at temperature
$T=1.6$ K unless specified otherwise.

###  I.3 Determination of stacking configuration

Figure S2: Second harmonic generation (SHG) measurements. (A) Zoom-in of the
device with WSe2 (orange) and MoSe2 (green) monolayers outlined. The colored
circles indicate the locations where SHG measurements were taken. (B) SHG
signal from each isolated monolayer and the heterobilayer region as a function
of polarization angle $\phi$. The larger signal from the heterobilayer
indicates that the constituent layers have AA stacking.

Second harmonic generation (SHG) can be used to determine the relative twist
between monolayer WSe2 and monolayer MoSe2 due to their lack of inversion
symmetry Kumar  _et al._ (2013); Li  _et al._ (2013). By tracking the
reflected signal at the second harmonic frequency of the incoming signal as
the linear polarizer rotates, crystallographic orientations of monolayer
portions of the device are identified as shown in Fig. S2. The SHG signal also
clarifies the stacking of the heterobilayer. For a bilayer, the second
harmonic fields add constructively for AA stacking and destructively for AB
stacking. The strong SHG signal over the heterobilayer portion of the device
demonstrates that it has AA stacking configuration (Fig. S2B).

###  I.4 Determination of the twist angle and filling factors

For homobilayers and heterobilayers with small lattice mismatch, the moiré
unit cell area diverges as the twist angle approaches zero. Twist angle
disorder can therefore generate local variability in the unit cell area. Thus,
while SHG measurements distinguish between AA and AB stacking, we apply a more
precise method to estimate the local twist angle in SET measurements. We
calibrate the gate-sample capacitance and therefore the conversion between
back gate voltage $V_{\textrm{G}}$ and carrier density $n$ based on the slope
of the Hofstadter states in a magnetic field $B$. We then identify one and two
holes per moiré unit cell ($\nu=-1,-2$) as the densities with the largest
peaks in inverse compressibility (the relative gate voltages they occur at in
different locations with distinct twist angles further confirms this filling
factor assignment). The $B=0$ intercept obtained from extrapolation of the
Hofstadter states also matches these assignments and confirms that the largest
peaks are at integer fillings. From the density at $\nu=-2$, we calculate the
twist angle according to
$n=\frac{4}{\sqrt{3}}\frac{\theta^{2}+\delta^{2}}{a_{\text{Mo}}^{2}}$, where
$\delta=1-\frac{a_{\text{W}}}{a_{\text{Mo}}}$ and
$a_{\text{W(Mo)}}=3.282\text{ \AA}\ (3.290\text{ \AA})$ is the lattice
constant for WSe2 (MoSe2).

###  I.5 Gap extraction

Figure S3: Inverse electronic compressibility $\text{d}\mu/\text{d}n$
background and determination of gaps at fractional fillings. (A) Local
$\text{d}\mu/\text{d}n$ as a function of moiré filling factor $\nu$ (red).
Purple and black dashed lines correspond to Methods 1 and 2, respectively, of
accounting for the negative compressibility background. The purple dashed line
is a linear fit to the overall negative background outside the peaks (Method
1), while the black dashed lines linearly connect the local minima surrounding
the incompressible peaks of the charge-ordered states (Method 2). (B) Chemical
potential $\mu$ as a function of $\nu$ from the raw $\text{d}\mu/\text{d}n$
(red) and after subtracting the purple linear fit from (A) (purple). (C-E) The
gap size $\Delta_{\nu}$ of each charge-ordered state as a function of magnetic
field $B$. These correspond to the step size in chemical potential at each
respective $\nu$ from the raw $\text{d}\mu/\text{d}n$ (C), that obtained using
Method 1 to account for the negative background (D), and that determined from
Method 2 (E).

The thermodynamic gap $\Delta_{\nu}$ of a given incompressible state at
filling $\nu$ is given by the corresponding step in chemical potential. We
extract the energy gaps by integrating $\text{d}\mu/\text{d}n$ with respect to
the density $n$,

$\displaystyle\Delta_{\nu}\equiv\Delta\mu$
$\displaystyle=\mu(n_{+})-\mu(n_{-})$
$\displaystyle=\int^{n_{+}}_{n_{-}}\left[\frac{\textrm{d}\mu}{\textrm{d}n}-\left(\frac{\textrm{d}\mu}{\textrm{d}n}\right)_{B}\right]\text{d}n,$

where $n_{+(-)}$ is the density at the local maximum (minimum) chemical
potential above (below) the state of interest, and
$\left(\frac{\textrm{d}\mu}{\textrm{d}n}\right)_{B}$ is the negative
compressibility background to be removed.

As noted in the main text, there is a strong background of negative
compressibility over large ranges of filling factor and magnetic field. Simply
integrating the raw measured $\text{d}\mu/\text{d}n$ would therefore result in
gaps that are suppressed relative to their true value. To account for the
background and obtain more accurate gap sizes, we implement a method (Method
1) analogous to that followed in Ref. Eisenstein  _et al._ (1994).
Specifically, a linear fit to the $\text{d}\mu/\text{d}n$ background
surrounding the incompressible charge ordered states at fractional filling is
subtracted from the raw $\text{d}\mu/\text{d}n$ (purple dashed line, Fig.
S3A). This method only assumes a slowly varying negative background and is
insensitive to details in the vicinity of the incompressible peaks. The
chemical potential before and after accounting for the negative background is
shown in red and purple, respectively, in Fig. S3B. The corresponding gaps
extracted from the raw and background-corrected chemical potential are plotted
as a function of magnetic field in Fig. S3C and Fig. S3D, respectively.

Finally, for completeness we also show a different method (Method 2) of
accounting for the background. Rather than linearly fitting to the background
far from the peaks, we identify the local minima in $\text{d}\mu/\text{d}n$
adjacent to each incompressible peak (black dashed lines, Fig. S3A).
Integrating the area under the peak but above the line connecting these local
minima yields the energy gaps shown in Fig. S3E. The gaps extracted from the
raw $\text{d}\mu/\text{d}n$ provide an effective lower bound, and those
extracted using Method 2 serve as an upper bound. We believe that Method 1 is
the most physically justified, and present the corresponding gaps in Fig. 3 of
the main text. We emphasize that the method used to account for the negative
background compressibility changes the gap sizes quantitatively but not
qualitatively. Specifically, for all three gap extraction methods, the gaps of
charged ordered states at fractional fillings increase with $B$ at magnetic
fields below the phase transition at 6 T and decrease with $B$ above this
critical field.

##  II Comparison between transport and compressibility measurements

Transport measurements across the device at $B=0$ (Fig. S4A) show prominent
resistance peaks that can be matched to gapped states visible in SET
measurements conducted in an intervening region of the sample (Fig. S4B).
However, the local electronic compressibility measurements are significantly
sharper and reveal additional charge ordered states at fractional fillings.
The mobility edge in transport also occurs further from the band edge than in
SET measurements.

Figure S4: Comparison of transport and $\text{d}\mu/\text{d}n$. (A) Two-probe
resistance $R$ as a function of back gate voltage $V_{\textrm{G}}$ at
temperature $T=4$ K. Inset: Circuit diagram showing the contacts used for the
measurement. (B) Local $\text{d}\mu/\text{d}n$ measurement at a location with
twist angle $\theta=1.35\degree$ at temperature $T=1.6$ K. Inset: The
approximate tip location of this measurement is indicated by the circle.

##  III Spatial dependence

Figure S5: Spatial dependence. (A-D) Measurements of local
$\text{d}\mu/\text{d}n$ as a function of $\nu$ and $B$ at four independent
locations which are each separated by at least 150 nm from each other and from
the data presented in Fig. 1B of the main text. The corresponding local twist
angles are $1.35\degree$ (A), $1.36\degree$ (B), $1.24\degree$ (C), and
$1.44\degree$ (D). A qualitatively similar phase diagram of competing
Hofstadter and charged-ordered states is evident in each location,
demonstrating the generality of the underlying physics. In (D), a constant
background compressibility of $8\times 10^{-12}$ meV cm2 was subtracted from
the raw AC signal, which exhibited a spurious enhancement due to poorer sample
charging. The subtracted background was benchmarked based on a simultaneous DC
measurement, and the uncertainty of this background is $<2\times 10^{-12}$ meV
cm2.

The general pattern of Hofstadter and charge-ordered phases are not unique to
a particular twist angle or location in our device. Measurements of inverse
compressibility as a function of carrier density and magnetic field at five
different locations of the sample demonstrate the robustness of these phases
and their phenomenology (Fig. 1b and Fig. S5). All of these locations show
magnetic field driven phase transitions indicated by an abrupt decrease of the
$\nu=-2$ gap and weakened or absent charge-ordered states between $-2<\nu<-1$
at high field, with near perfect agreement in two independent locations that
have almost identical twist angles (Fig. S5A,B). Other features, such as
regions of negative compressibility and the overarching pattern of where
Hofstadter and charge-ordered states respectively dominate, are linked to
filling factor and occur at similar locations within the phase diagram,
showing that they are caused by intrinsic local physics in the heterobilayer
and do not simply occur at a specific gate voltage or field.

Beyond this general agreement, we note qualitative differences at larger and
smaller twist angles (Fig. S5C,D). At smaller twist angles, we resolve fewer
Hofstadter states emanating from $s=-3$, and the $t=1$, $s=-2$ state is
absent. In addition, the $\nu=-3$ gap is stronger, appears for a larger range
in $B$, and persists all the way down to $B=0$ (Fig. S5C). In contrast, at
larger twist angles, we observe an extra $t=-1$, $s=-1$ state (Fig. S5D) and a
$\nu=-2$ gap that closes fully at intermediate fields. Thus, the twist angle
tunes the competition between Hofstadter and charge-ordered states. We ascribe
this to differences in the electronic bandwidth, which is highly sensitive to
the twist angle. Smaller twist angles produce flatter bands, favoring charge-
ordered phases, while the more dispersive bands at larger twist angles favor
Hofstadter states.

Figure S6: Magnetic field dependence of charge-ordered states at other
locations. (A,B) Gaps at fractional fillings extracted from Fig. 1B in the
main text. (C,D) Gaps at fractional fillings extracted from Fig. S5B. Gray
regions indicate fields at which we cannot extract quantitatively accurate
gaps because the lock-in amplifier overloaded during measurement.

In Fig. S6, we show the magnetic field dependence of the thermodynamic gaps of
charge-ordered states in two different locations (those shown in Fig. 1B of
the main text and in Fig. S5B) that have similar twist angle to Fig. 3C but
are sufficiently far away to represent independent measurements. Contradictory
to theoretical predictions for a single band Hubbard model in a triangular
lattice Pan  _et al._ (2020); Morales-Durán  _et al._ (2022), the $\nu=-3/2$
gap is generally larger than those at $-4/3$ and $-5/3$ in our device. This
contrasts with several previous measurements in other materials systems where
states at multiples of $\nu=1/3$ are strongest Tang  _et al._ (2020); Li  _et
al._ (2021a); Jin  _et al._ (2021); Li  _et al._ (2021b). However, we note
that a larger gap at $\nu=-3/2$ has also been previously reported in
measurements of aligned WS2/WSe2 heterobilayers between $1<|\nu|<2$ Xu  _et
al._ (2020); Huang  _et al._ (2021).

Details of the magnetic field evolution of charge-ordered states at fractional
fillings depend on exact location. The increase in the charge gap at
$\nu=-3/2$ as a function of magnetic field is reproducible in multiple
locations. However, it is harder to identify generic behavior for other
fillings. In particular, the gaps at $-4/3$ and $-5/3$ in Fig. S6B increase up
to $B=1.6$ T beyond which the gap dependence is more ambiguous. In Fig. S6D,
there is not a large change in gap size over the magnetic fields where we can
extract reliable data. We note that while spin physics has been predicted for
ground states on a triangular moiré lattice Wu  _et al._ (2018a); Pan  _et
al._ (2020); Morales-Durán  _et al._ (2022), the very large $U/t\approx
O(10^{2})$ implies that the energy scale associated with the spin physics,
$J=t^{2}/U\approx 1$ mK, is much smaller than the experimental temperature,
and is therefore unlikely to be relevant to the gap scaling with magnetic
field. An applied magnetic field squeezes the wave function of an electron in
a harmonic potential well, which decreases $|t|$ by reducing the inter-well
wavefunction hybridization. This effect could in principle lead to an increase
in the gaps of the charge-ordered states with field. However, for realistic
parameters of this system, the cyclotron frequency $\omega_{c}=eB/(2\pi
m^{*})\ll\omega_{0}$ is very small ($\omega_{0}$ is the natural frequency of
the harmonic potential), indicating that this effect is also negligible.

##  IV Density function theory (DFT) and Hartree-Fock (HF) calculations

We study the TMD heterobilayer WSe2/MoSe2 with a small twist angle starting
from AA stacking, where every metal (W) or chalcogen (Se) atom on the top
layer is aligned with the same type of atom on the bottom layer. Within a
local region of a twisted bilayer, the atom configuration is identical to that
of an untwisted bilayer (the lattice mismatch between WSe2 and MoSe2 is less
than 0.2%), where one layer is laterally shifted relative to the other layer
by a corresponding displacement vector ${\bm{d}}_{0}$. For this reason, the
moiré band structures of twisted TMD bilayers can be constructed from a family
of untwisted bilayers at various ${\bm{d}}_{0}$, all having $1\times 1$ unit
cell. Our analysis thus starts from untwisted bilayers.

In particular,
${\bm{d}}_{0}=0,\left(-{\bm{a}}_{1}+{\bm{a}}_{2}\right)/3,\left({\bm{a}}_{1}+{\bm{a}}_{2}\right)/3$,
where ${\bm{a}}_{1,2}$ is the primitive lattice vector for untwisted bilayers,
correspond to three high-symmetry stacking configurations of untwisted TMD
bilayers, which we refer to as MM, XM, MX. In MM (MX) stacking, the M atom on
the top layer is locally aligned with the M (X) atom on the bottom layer,
likewise for XM. The bilayer structure in these stacking configurations is
invariant under three-fold rotation around the $z$ axis.

Density functional calculations are performed using generalized gradient
approximation with SCAN+rVV10 van der Waals density functional Peng  _et al._
(2016), as implemented in the Vienna ab initio Simulation Package Kresse and
Furthmüller (1996). Pseudopotentials are used to describe the electron-ion
interactions. We first construct the zero-twisted angle WSe2/MoSe2 bilayer at
MM, MX, and XM lateral configurations with vacuum spacing larger than 20 Å to
avoid artificial interaction between the periodic images along the $z$
direction. The structure relaxation is performed with force on each atom less
than 0.001 eV/Å. We use $12\times 12\times 1$ for structure relaxation and
self-consistent calculation. The more accurate SCAN+rVV10 van der Waals
density functional gives the relaxed layer distances as 6.46 Å, 6.45 Å and
6.92 Å for MX, XM and MM stacking structures, respectively. By calculating the
work function from electrostatic energy of converged charge density, we obtain
the band structure of MM, MX, and XM-stacked bilayers, with reference energy
$E=0$ chosen to be the absolute vacuum level, shown in Fig S7.

Figure S7: DFT band structure for high symmetry stacking structures.

The low-energy moiré band structure for the K valley valence bands is
described by the continuum model Hamiltonian Wu  _et al._ (2018b),

$\mathcal{H}_{\bm{g}^{\prime},\bm{g}}(\bm{k})\equiv-\frac{\hbar^{2}(\bm{k}+\bm{g})^{2}}{2m^{*}}\delta_{\bm{g}^{\prime},\bm{g}}+V\sum_{j=1}^{3}(e^{i\phi}\delta_{\bm{g}^{\prime},\bm{g}+\bm{g}_{j}}+e^{-i\phi}\delta_{\bm{g}^{\prime},\bm{g}-\bm{g}_{j}}).$
(1)

We use $a_{M}=13.5$ nm, $m^{*}=0.5m_{e}$, and $(V,\phi)=(12.3$
meV,$-125.1^{\circ})$ as obtained from the DFT bands at the three stacking
regions. In the limit $\frac{\hbar^{2}}{2m^{*}a_{M}^{2}}/V\rightarrow\infty$,
the model becomes an electron gas and, in the opposite limit, it approaches an
array of isolated quantum dots. The long moiré wavelength of our sample
$a_{M}=13.5$ nm puts it closer to the quantum dot array limit with
$\frac{\hbar^{2}}{2m^{*}a_{M}^{2}}/V=0.034$.

The lowest several non-interacting moiré bands in our sample originate from
two-dimensional harmonic oscillator-like orbitals on a triangular lattice of
$XM$ stacking sites. As holes are doped into moiré bands, it is possible that
their mutual Coulomb repulsion could cause them to localize to other regions
of the moiré unit cell, a scenario known as charge transfer Zhang  _et al._
(2020). The small value of $\frac{\hbar^{2}}{2m^{*}a_{M}^{2}}/V$ in our system
suggests that holes will remain localized to triangular lattice sites between
filling factors $-2<\nu<-1$. To test this, we calculate the self-consistent,
unrestricted HF ground state at spin-projected filling factors
$\nu_{\uparrow}=-1,\nu_{\downarrow}=0$ and at
$\nu_{\uparrow}=\nu_{\downarrow}=-1$. In the HF approximation, the continuum
model Hamiltonian becomes Hu and MacDonald (2021)

$\displaystyle\begin{split}\mathcal{H}_{s^{\prime}\bm{g}^{\prime},s\bm{g}}(\bm{k})&\equiv\delta_{s^{\prime},s}[-\frac{\hbar^{2}(\bm{k}+\bm{g})^{2}}{2m^{*}}\delta_{\bm{g}^{\prime},\bm{g}}+V\sum_{j=1}^{3}(e^{i\phi}\delta_{\bm{g}^{\prime},\bm{g}+\bm{g}_{j}}+e^{-i\phi}\delta_{\bm{g}^{\prime},\bm{g}-\bm{g}_{j}})]\\\
&+\frac{1}{A}\sum_{\bm{k}^{\prime},\bm{g}^{\prime\prime}}[\delta_{s^{\prime},s}\sum_{s^{\prime\prime}}V_{s,s^{\prime\prime}}(\bm{g}^{\prime}-\bm{g})\rho_{s^{\prime\prime},\bm{g}^{\prime}+\bm{g}^{\prime\prime};s^{\prime\prime},\bm{g}+\bm{g}^{\prime\prime}}(\bm{k}^{\prime})-V_{s^{\prime},s}(\bm{g}^{\prime\prime}+\bm{k}^{\prime}-\bm{k})\rho_{s^{\prime},\bm{g}^{\prime}+\bm{g}^{\prime\prime};s,\bm{g}+\bm{g}^{\prime\prime}}(\bm{k}^{\prime})].\end{split}$
(2)

The first line of this expression contains the single-particle terms and the
second line contains the interaction terms. The first interaction term is the
local Hartree potential an the second is the Fock non-local exchange
potential. Here $s^{\prime}$ and $s$ are spin indices,
$\rho_{s^{\prime}\bm{g}^{\prime};s\bm{g}}\equiv\bra{\Psi_{HF}}c^{\dagger}_{s^{\prime},\bm{k}+\bm{g^{\prime}}}c_{s,\bm{k}+\bm{g}}\ket{\Psi_{HF}}$
is the single-particle density matrix in the basis of simultaneous $\bm{p}$,
$S_{z}$ eigenstates, and $\ket{\Psi_{HF}}$ is the single-slater-determinant
ground state. We perform the calculation with an $9\times 9$ Monkhorst-Pack
Brillouin zone mesh and the lowest $13$ shells of reciprocal lattice vectors.
We do not explore states with broken translation symmetry which are unlikely
at integer filling factors. We fix the $S_{z}$ eigenvalue of our state at the
outset of the calculation so that
$\rho_{s^{\prime},\bm{g}^{\prime};s,\bm{g}}\propto\delta_{s^{\prime},s}$. We
use the single-gate-screened Coulomb interaction in the image charge
approximation,
$V_{s^{\prime},s}(\bm{q})=\frac{e^{2}}{2\epsilon\epsilon_{0}}\frac{(1-e^{-2dq})}{q}$
where $d=30$ nm is the sample-gate distance and $q\equiv|\bm{q}|$.
Additionally, we choose $\epsilon=5$ to account for the hexagonal boron
nitride dielectric environment of the superlattice.

In Fig. S8A we plot the non-interacting continuum model band structure. The
top (flat) band originates from an $s$ orbital and the next pair of
(dispersive) bands originate from two $p$ orbitals at each triangular lattice
site. We note that the bandwidth is affected by screening from the underlying
filled bands. Hartree terms tend to weaken the moiré potential and therefore
lead to more delocalized wavefunctions and stronger inter-site orbital
hybridization, causing the bandwidth to be increased and filling-dependent.

In Fig. S8B, we plot the hole number density of the HF ground state at spin-
projected filling factors $\nu_{\uparrow}=-1,\nu_{\downarrow}=0$, and in Fig.
S8C, the same quantity at $\nu_{\uparrow}=\nu_{\downarrow}=1$. In both cases,
the charge density within the moiré unit cell is concentrated to the
triangular lattice sites of $V_{0}(\bm{r})$ maxima. The charge density within
the triangular lattice site at $\nu_{\uparrow}=\nu_{\downarrow}=-1$ is spread
out compared to when $\nu_{\uparrow}=-1,\nu_{\downarrow}=0$ because of the
repulsion between the two holes on each site. Although we do not enforce
$\psi_{n\uparrow{\bm{k}}}=\psi_{n\downarrow{\bm{k}}}$ (in other words, our HF
calculation is unrestricted), our self-consistent solution in this case
satisfies this condition. These results serve as additional evidence that our
system is appropriately modeled as a triangular lattice and that the charge
transfer scenario does not occur at least for filling factors $\nu\geq-2$.

Figure S8: Hartree-Fock (HF) study of charge configuration between filling
factors $\nu=-1$ and $\nu=-2$. (A) Non-interacting continuum model band
structure. (B) Hole number density of the HF ground state at
$\nu_{\uparrow}=-1,\nu_{\downarrow}=0$ and (C) at
$\nu_{\uparrow}=\nu_{\downarrow}=-1$. The hole number density in both cases is
localized to the $XM$ triangular lattice sites within the moiré unit cell.

Finally, we comment on the possible role of $\Gamma$ valley bands. At the $MX$
and $XM$ regions, the DFT band structure (Fig. S7) suggests that the $\Gamma$
pocket is competitive in energy with the $K$ pocket. While a continuum model
describing both $K$ and $\Gamma$ valleys show that the first band is indeed
from $K$ valley, the first $\Gamma$ valley band is only $\sim 10$ meV away. We
note that the competition of $\Gamma$ and $K$ valley bands in DFT is very
delicate, and depends sensitively on details such as choice of van der Waal
density functional and relaxed lattice structure. In our analyses, we have
assumed that the $K$ valley bands are energetically favored and the $\Gamma$
band can be safely ignored. This is consistent with prior reports that have
shown the valence band extrema reside at the $K$ valley for AA stacking
configuration Wilson  _et al._ (2017); Seyler  _et al._ (2019); Tran  _et al._
(2019); Barré  _et al._ (2022). We cannot conclusively rule out the
possibility that the $\Gamma$ bands are energetically favored in our
experiment, but we remark that the underlying physical picture would remain
largely unchanged. The valley structure can be probed directly in future
experiments via magnetic circular dichroism or through the application of an
in-plane magnetic field.

For completeness, we additionally remark that an alternative (although less
likely) scenario for the origin of the dispersive and flat bands we observe
could arise from competition between the $\Gamma$ and $K$ valleys. In this
scenario, the $\Gamma$ valley bands, which have larger effective mass and feel
a stronger moiré potential, play the role of the flat bands, while the $K$
valley bands play the role of the dispersive bands. However, this scenario
requires fine tuning (the first $\Gamma$ and $K$ bands must be very close in
energy) and fails to explain some observations. First, the presence of
observable LL gaps in the dispersive band is at odds with the flatness of the
first $K$ band (without significant interaction broadening). Indeed, in other
moiré TMD systems with smaller moiré period, charge order is generally favored
even though the smaller moiré unit cell corresponds to larger bandwidth
relative to interaction energy. Additionally, the first $K$ band, associated
with an $s$ orbital, should display more prominent hole-like LL gaps (see Fig.
S12), which is at odds with the observation of only electron-like LL gaps in
the first dispersive band.

##  V Phenomenological Stoner model

The set of parameters used in the phenomenological model to generate Fig. 2D
in the main text are as follows: density of states of the dispersive band
$D=1/9t$, effective interaction $U_{fd}=1/3D=3t$, band separation $\Delta=1/D$
(see Fig. 2C), and $g=0.22W/$Tesla. Here, $t$ is the hopping matrix element of
the $p$-orbitals in the tight-binding Hamiltonian and is an adjustable
parameter. Parameters are chosen such that the model quantitatively matches
the experimental phase diagram in Fig. 1B, but we note that a realistic value
of $t\approx 0.2$ meV yields a $g$-factor of 8$\mu_{B}$, which is of similar
magnitude to values from other measurements in the same system Seyler  _et
al._ (2019).

In this model, the width of the reentrant region in which both flat and
dispersive bands are partially filled is constant and given by
$\delta\nu=(2DU_{fd}-1)/(DU_{fd}-1)$. In the experimental phase diagram, the
width of this region is also roughly constant $\delta\nu\approx 1/2$. This
implies that the effective interaction parameter must be $U_{fd}\approx
1/(3D)$. Note, however, that this does not imply that the on-site Hubbard
interaction is $1/(3D)$, as the parameter $U_{fd}$ is generically not equal
(or even proportional) to the Hubbard interaction. We substantiate this
statement through our analysis of an integrable microscopic model in the next
section. The model also predicts that for smaller angle (larger $D$), the
crossover from filling the flat band to filling the dispersive band should
occur at smaller $B$ for given $(\nu_{f},\nu_{d})$. This is also consistent
with our experimental observations in Fig. S5.

Compressibility
$\frac{\textrm{d}\mu}{\textrm{d}\nu}=\frac{\textrm{d}^{2}E}{{\textrm{d}\nu}^{2}}$
can be calculated from the phenomenological model, as is shown in Fig. S9.
Within the region where reentrant charge-ordered states appear, a constant
negative $d\mu/d\nu=\frac{-DU_{fd}^{2}}{1-2DU_{fd}}<0$ is obtained.
Experimentally, the negative compressibility within this region relative to
that associated with the adjacent dispersive bands is evident at base
temperature in Fig. 4A as well as in Fig. S5(A,B). The charge-ordered states
within this reentrant region consist of states originating from both flat and
dispersive bands. The degree to which these distinct states hybridize is an
interesting open question to be explored. If the bands do not hybridize, the
charge-ordered states would involve simultaneous localization of carriers from
both flat and dispersive bands. Alternatively, it is possible that the bands
hybridize and produce composite low-energy states with large effective mass,
i.e. heavy fermions. Finally, we comment on the different behavior we observe
while filling the same band at high field between $-2<\nu<-1$ and at low field
between $-3<\nu<-2$ (Fig. 1B). This difference may be caused by screening from
the underlying filled bands. As noted above, this broadens the dispersion due
to Hartree screening. It may also enhance effective Coulomb screening and
therefore weaken electronic interactions.

Figure S9: Phenomenological Stoner model. $\text{d}\mu/\text{d}n$ calculated
from the phenomenological Stoner model.

##  VI Hubbard model

We analyze a Hubbard model description for the experimental system in the hole
density range $-2\leq\nu\leq-1$ in a magnetic field. In this region, we assume
the the spin-up $s$-orbital flat band is fully filled due to the spin Zeeman
effect, and the relevant degrees of freedom are the spin-down $s$-orbital band
and dispersive $p$-orbital bands, whose competition is tuned by the magnetic
field. The presence of an insulating state at $\nu=-2$ suggests that only a
single $p$ band (the one derived from the spin-up $p$ state with orbital
angular momentum quantum number $\ell_{z}=+1$) is being filled.

We model this system as a two-orbital tight binding model on a triangular
lattice, where the two orbitals correspond to the flat and dispersive bands,
respectively. To obtain a tractable model with interactions, we further take
the limit where the flat band is perfectly flat and the interaction is only
on-site via a Hubbard term. The resulting Hamiltonian is

$\mathcal{H}=-t\sum_{\langle
i,j\rangle}d^{\dagger}_{i}d_{j}+(\Delta+3|t|-gB)\sum_{i}n^{d}_{i}+U\sum_{i}n^{d}_{i}n^{f}_{i}$
(3)

where $d^{\dagger}_{i},d_{i}$ are the fermionic creation and annihilation
operators for the dispersive band on site $i$,
$n^{d}_{i}=d^{\dagger}_{i}d_{i}$, and the sum $\langle i,j\rangle$ is over
nearest neighbor sites on the triangular lattice. Finally, $n^{f}_{i}=0,1$ is
the density of the upper flat band fermions on site $i$ which, due to the
absence of a kinetic term, is a good quantum number. The hopping parameter $t$
represents the kinetic energy of the dispersive band, and is negative for the
$p$ orbital band in question. The magnetic field is modeled as a Zeeman term
$gB$. We have chosen the constant potential offset $(\Delta+3|t|-gB)$ so that
the band bottom is at $\Delta-gB$ in order to make a direct comparison with
the phenomenological Stoner model.

The on-site repulsion $U$ is estimated to be much larger than $t$. Taking a
harmonic approximation for the moiré potential with a moiré period of
$\approx$ 14 nm and effective mass $m^{*}=0.5m_{e}$, leads to an effective
harmonic oscillator lengthscale $\lambda\approx 2$ nm, and corresponding
Hubbard energy scale
$U\approx\frac{e^{2}}{4\pi\epsilon\epsilon_{0}\lambda}\approx(700/\epsilon)$
meV, where $\epsilon$ is the relative dielectric constant. Meanwhile, the non-
interacting bandwidth of the $p$ band from the continuum model is estimated to
be only $W=9|t|\approx 2$ meV. Although the bandwidth will be renormalized due
to Coulomb interaction with the filled $s$ band, it is unlikely to change the
fact that the local Hubbard interaction is much larger than $t$. Thus, we take
the limit $U/t\gg 1$, in which there remains only a single energy scale $|t|$
for total density $n^{f}+n^{d}\leq 1$, where $n^{f,d}=|\nu_{f,d}|$. In this
limit, double occupancy is forbidden. As there is no double occupancy, the
interaction energy between the heavy and light bands is purely kinetic: the
sites on which $n^{f}_{i}=1$ act as ‘defects’ for which the dispersive
fermions must avoid, thus paying the cost in kinetic energy.

The exact ground state(s) of the above model at a fixed $n^{f}$ and $n^{d}$
correspond to states with a particular configuration of $n^{f}_{i}$ which
minimizes the kinetic energy of the dispersive band. In reality, however, this
effect will be overshadowed by long range Coulomb interactions, which will
play a significant role in determining the final charge configuration leading,
among other things, to the formation of incompressible generalized Wigner
crystal states at commensurate fillings. For our current analysis, we take a
middle ground and simply average over many configurations of $n^{f}_{i}$ with
equal probability.

We numerically simulate this Hamiltonian on a $L\times L$ triangular lattice
torus. To do this, we first fix particle numbers $N^{f},N^{d}=1\dots L^{2}$,
with $N^{f}+N^{d}\leq L^{2}$. We then solve for the ground state energy of
$\mathcal{H}$ with $(\Delta+3|t|-gB)=0$ for $N_{\textrm{samp}}$ random
configurations of $n^{f}_{i}\in\\{0,1\\}$, such that
$\sum_{i}n^{f}_{i}=N^{f}$. The ground state energy for each sample is obtained
by diagonalizing the quadratic Hamiltonian for $d$ fermions and filling
$N^{d}$ lowest energy bands. To decrease finite size effects, we also average
over fluxes threaded through the torus. The final result is the configuration-
averaged energy per unit cell $E_{0}(n^{d},n^{f})$, where
$n^{(f,d)}=N^{(f,d)}/L^{2}$. Note that there is a single energy scale set by
$W=9|t|$ in the large $U/t$ limit. For a given total density
$n^{tot}=n^{d}+n^{f}$ and magnetic field $B$, the optimal density distribution
can then be obtained by minimizing
$E(n^{tot},B)=\min_{n^{d}\in[0,n^{tot}]}\left\\{E_{0}(n^{d},n^{tot}-n^{d})+(\Delta+3|t|-gB)n^{d}\right\\}$.
Choosing $\Delta=W$ and $g=0.22W/$Tesla as in the Stoner model results in the
phase diagram shown in Fig. 2F of the main text.

The final result is a relatively sharp transition at a critical field
$B_{c}\approx 6$ T between filling of the flat and dispersive bands. The
critical field varies little with filling, except for a small upturn near
$\nu=-1-n^{tot}=-2$. Near this upturn, compressibility becomes very negative,
as shown in Fig S10A, in excellent qualitative agreement with the experimental
compressibility measurements in this range.

Finally, we can obtain the effective interaction $U_{fd}$ as a function of
$n^{f}$ and $n^{d}$ in our effective model. The energy function $E_{0}$ can be
decomposed as
$E_{0}(n^{d},n^{f})=E_{0}(n^{d},0)+E_{0}(0,n^{f})+U_{\mathrm{eff}}(n^{d},n^{f})n^{d}n^{f}$.
The numerically obtained $U_{\mathrm{eff}}(n^{d},n^{f})/W$ as a function of
the dispersive and flat band density is shown in Fig S10B. The value of
$U_{\mathrm{eff}}/W\approx 0.3$ in a large part of the diagram, which is
consistent with the $U_{fd}=W/3$ suggested by the Stoner model in combination
with the experimental width $\delta\nu=1/2$ of the reentrant region (despite
this model not being directly applicable for $\nu<-2$).

Figure S10: Filling-dependent interaction strength and compressibility in the
Hubbard model. (A) Effective interaction strength $U_{eff}$ calculated as a
function of flat and dispersive band filling $n^{f,d}$. (B) Numerically
evaluated compressibility within the Hubbard model.

##  VII Singlet-triplet quantum dot transition

By expanding the moiré potential near its maxima to leading (quadratic) order,
we approximate the effective Hamiltonian for a hole isolated to a single moiré
potential maximum, as that of a circular harmonic oscillator. The effective
Hamiltonian for two Coulomb-interacting holes isolated to the same potential
well in a magnetic field is then

$\displaystyle\begin{split}H_{QD}&=\hbar\Omega_{+}(a^{\dagger}_{{1+}}a_{1+}+a^{\dagger}_{{2+}}a_{2+}+1)+\hbar\Omega_{-}(a^{\dagger}_{{1-}}a_{1-}+a^{\dagger}_{{2-}}a_{2-}+1)+\frac{e^{2}}{4\pi\epsilon\epsilon_{0}|\bm{r}_{1}-\bm{r}_{2}|}-g\mu_{B}BS_{z}.\end{split}$
(4)

Here $a_{\alpha\pm}\equiv(a_{\alpha x}\mp ia_{\alpha y})/\sqrt{2}$ lowers the
energy of particle $\alpha$ by $\Omega_{\pm}$ and changes its $L_{z}$
eigenvalue by $\pm\hbar$.
$\Omega_{\pm}=\sqrt{\omega^{2}+(\frac{\omega_{c}}{2})^{2}}\pm\frac{\omega_{c}}{2}$,
$\omega_{c}\equiv eB/m^{*}$ is the effective cyclotron frequency, and
$\omega\equiv\sqrt{\frac{16\pi^{2}V\cos{\phi}}{a_{M}^{2}m^{*}}}$ is the
natural harmonic oscillator frequency determined the moiré potential. $g$ is
the effective spin $g$-factor. From the continuum model parameters ($V=12.3$
meV, $m^{*}=0.5m_{e}$, $\phi=-125.1^{\circ}$, $a_{M}=13.5$ nm) we find
$\hbar\omega=40.3$ meV and $l=1.95$ nm.

The singlet ground state has spin and orbital angular momentum quantum numbers
$(S,L)=(0,0)$ and the triplet $(S,L)=(1,1)$. The magnetic field couples to the
two particle state through linear and quadratic terms in $\hbar\omega_{c}$
appearing in $\Omega_{+-}$ and through Zeeman coupling. In the relevant
parameter regime, $\omega_{c}\ll\omega_{0}$. For this reason, the spin singlet
state is decoupled from the magnetic field [with leading corrections of
$\mathcal{O}(\omega_{c}/\omega)^{2}$] and the triplet state is coupled
linearly to the magnetic field through its orbital and spin angular momentum.
In the absence of the Coulomb interaction, the singlet-triplet gap is
$\Delta_{st}(B)\equiv E_{0}(S=1)-E_{0}(S=0)=\hbar\Omega_{-}-g\mu_{B}B$. The
singlet-triplet ground state transition occurs at $B_{c}$ such that
$\Delta_{st}(B_{c})=0$. As the strength of the Coulomb interaction increases,
$\Delta_{st}(B)$ monotonically decreases for all $B$ and thus so too does
$B_{c}$.

To determine $B_{c}$ in the presence of a finite Coulomb interaction, we solve
this two-body problem using exact diagonalization following the approach of
Zeng and MacDonald (2022), to which we refer the reader for a technical
discussion of the calculation. We plot the critical magnetic field for the
singlet-triplet transition $B_{c}(\epsilon^{-1})$ as a function of Coulomb
interaction strength for several values of $g$. In the absence of the Coulomb
interaction, $B_{c}\gg 6$ T where $6$ T is the approximate value of the
observed phase transition in the filling range $-2>\nu>-1$. In the presence of
a finite Coulomb interaction, however, $B_{c}$ monotonically decreases,
approaching values comparable to $6$ T. We note that this singlet-triplet
transition has been studied theoretically and observed experimentally in GaAs
quantum dots Wagner  _et al._ (1992); Ashoori  _et al._ (1993).

Finally, we comment on our observation that $B_{c}$ is smaller in regions of
smaller twist angle. This is consistent with the predicted scaling of $B_{c}$
with twist angle in our quantum dot model. Since, as noted above,
$\omega_{c}\ll\omega_{0}$, $\Delta_{st}(B)$ is a linear function of $B$. In
the absence of Coulomb interactions, $\Delta_{st}(B=0)=\hbar\omega\propto
a_{M}^{-1}\propto\theta$. The interaction strength is characterized by the
dimensionless ratio of energy scales
$\frac{e^{2}/(4\pi\epsilon\epsilon_{0}l)}{\hbar\omega}\propto
a_{M}^{1/2}\propto\theta^{-1/2}$. As a function of interaction strength,
$\Delta_{st}$ decreases monotonically. Since the non-interacting gap decreases
with decreasing $\theta$ and the coupling constant, which decreases the gap,
increases with decreasing $\theta$, $\Delta_{st}(B=0)$ and therefore also
$B_{c}$ monotonically decrease with decreasing $\theta$.

Figure S11: Critical magnetic field for singlet-triplet transition of the two-
hole quantum dot. Value of the applied magnetic field at which the two-hole
ground state within a moiré potential well undergoes a singlet-triplet
transition as a function of Coulomb interaction strength $\epsilon^{-1}$
determined by exact diagonalization of two Coulomb-interacting holes in a
circular harmonic potential.

##  VIII Single-particle Hofstadter spectrum

Figure S12: Single-particle Hofstadter spectrum. (A-B) Single-particle and
spinless Hofstadter spectrum calculated from the continuum model for the flat
$s$ (A) and dispersive $p$ (B) bands.

In Fig. S12, we show the single-particle Hofstadter spectrum for the first
three bands of a single spin species, computed from the continuum model. We
work with the continuum model,

$H(\mathbf{k})=-\frac{\mathbf{k}^{2}}{2m}+V(\mathbf{r})$ (5)

where
$V(\mathbf{r})=2v\sum_{i=1,3,5}\cos(\mathbf{g}_{i}\cdot\mathbf{r}+\phi)$. A
magnetic field is introduced via minimal coupling by substituting
$\mathbf{k}\rightarrow\bm{\pi}\equiv\mathbf{k}-\mathbf{A}$, where
$\mathbf{A}=(B/2)(x\hat{\mathbf{y}}-y\hat{\mathbf{x}})$ is the vector
potential in the symmetric gauge. The gauge-invariant operators $\pi_{x}$ and
$\pi_{y}$ satisfy the commutation relation $[\pi_{x},\pi_{y}]=iB$. From these,
we define Landau level raising operator
$a^{\dagger}=(\pi_{x}-i\pi_{y})/\sqrt{2B}$. We also define the operators
$\mathbf{K}=\mathbf{k}+\mathbf{A}$, which are gauge-dependent and satisfy
$[\pi_{i},K_{i}]=0$ and $[K_{x},K_{y}]=-iB$. We work in the eigenbasis of
$n=a^{\dagger}a$ and $K_{x}$, spanned by the states $|n,K_{x}\rangle$ where
$n=0,1,2,\dots$ and $K_{x}\in\mathbb{R}$.

The kinetic energy is diagonal in this basis is independent of $K_{x}$

$-\frac{\bm{\pi}^{2}}{2m}\ket{n,K_{x}}=-\frac{B}{m}\left(n+\frac{1}{2}\right)\ket{n,K_{x}}$
(6)

The potential term can be decomposed into terms of the form
$e^{i\mathbf{Q}\cdot\mathbf{r}}$, which act in this basis as

$e^{i\mathbf{Q}\cdot\mathbf{r}}\ket{n,K_{x}}=\sum_{n^{\prime}}D_{n^{\prime}n}(\frac{Q_{x}+iQ_{y}}{\sqrt{2B}})e^{-i\frac{Q_{y}}{B}\left(K_{x}+\frac{Q_{x}}{2}\right)}\ket{n^{\prime},K_{x}+Q_{x}}$
(7)

where

$D_{n^{\prime}n}(z)=\begin{cases}z^{n^{\prime}-n}e^{-\frac{|z|^{2}}{2}}\sqrt{\frac{n^{\prime}!}{n!}}L^{(n^{\prime}-n)}_{n}(|z|^{2})&n^{\prime}\geq
n\\\
(-z^{*})^{n-n^{\prime}}e^{-\frac{|z|^{2}}{2}}\sqrt{\frac{n!}{n^{\prime}!}}L^{(n-n^{\prime})}_{n^{\prime}}(|z|^{2})&n^{\prime}<n\end{cases}$
(8)

The reciprocal lattice vectors can be expressed as

$\mathbf{g}_{1}=(0,2\Delta_{y});\;\;\mathbf{g}_{3}=(-\Delta_{x},-\Delta_{y});\;\;\mathbf{g}_{5}=(\Delta_{x},-\Delta_{y})$
(9)

where $\Delta_{x}=\sqrt{3}g/2,\Delta_{y}=g/2$, and $g=4\pi/(\sqrt{3}a_{M})$.
In terms of $\Delta_{x,y}$, it becomes clear that the potential term only
mixes $K_{x}$ with $K_{x}\pm\Delta_{x}$, and is periodic under
$K_{x}\rightarrow K_{x}+p\Delta_{x}$ if $B=p\Delta_{y}\Delta_{x}/(2\pi q)$,
for integers $p,q$. We therefore introduce a new basis

$|n,j;x_{0},k_{0}\rangle=\mathcal{N}\sum_{\ell\in\mathbb{Z}}e^{\frac{2\pi
ik_{0}}{p}(x_{0}+j+\ell p)}|n,K_{x}=\Delta_{x}(x_{0}+j+\ell p)\rangle$ (10)

up to a normalization factor, where $x_{0}\in[0,1]$ and $k_{0}\in[0,1]$, and
$j=1,\dots,p$ with $j\equiv j+p$. Each term in the potential
$e^{i\mathbf{Q}\cdot\mathbf{r}}$ can be characterized by two integers
$q_{x}=Q_{x}/\Delta_{x}$ and $q_{y}=Q_{y}/\Delta_{y}$, which act in this basis
as

$e^{i\mathbf{Q}\cdot\mathbf{r}}|n,j;x_{0},k_{0}\rangle=\sum_{n^{\prime}}D_{n^{\prime}n}(\frac{Q_{x}+iQ_{y}}{\sqrt{2B}})e^{-\frac{2\pi
iq}{p}q_{y}(x_{0}+j+\frac{q_{x}}{2})-2\pi
ik_{0}\frac{q_{x}}{p}}\ket{n,j+q_{x};x_{0},k_{0}}$ (11)

In practice, we truncate the Landau level basis to some large but finite
number $n=0,\dots,N_{\mathrm{max}}$, which must be chosen to not affect the
low energy eigenvalues of interest. For each fixed $x_{0}$, $k_{0}$, the
resulting Hamiltonian is $pN_{\mathrm{max}}\times pN_{\mathrm{max}}$
dimensional and can be exactly diagonalized to obtain a set of eigenvalues
$E^{x_{0}k_{0}}_{i}$. The density of states is obtained by averaging over
$x_{0}$ and $k_{0}$ with a factor of $1/(2q)$,

$\mathrm{DOS}(E)=\frac{1}{2qA_{uc}}\left\langle\sum_{i}\delta(E-E^{x_{0}k_{0}}_{i})\right\rangle_{x_{0}k_{0}}$
(12)

where $A_{uc}=\sqrt{3}a_{M}^{2}/2$ is the moiré unit cell area.

The results are shown in Fig. S12 for the first three moiré bands (of a single
spin species) with a moiré period of $a_{M}=13.5$ nm, effective mass
$m=0.5m_{e}$, potential strength $v=12.3$ meV and phase factor
$\phi=-125.1^{\circ}$. For a given $p,q$, we keep $N_{\mathrm{max}}=\lfloor
100q/p\rfloor$ states which is sufficiently large to avoid cutoff effects in
the first three bands, and average over randomly chosen $x_{0},k_{0}\in[0,1]$
until a sufficiently smooth DOS is achieved.

The Hofstadter spectra reveals many features which are consistent with
experimental observations at low density. The first band is derived from the
$s$ orbital of the moiré site and is exceptionally flat. The energy gaps
between Landau levels (Hofstadter subbands) within the first band is on the
order of $0.01$ meV, much too small to be detected experimentally. Indeed, in
the density range $-2\leq\nu\leq 0$ for $B<6$ T where the flat band is being
filled, no Hofstadter gaps are observed.

The second and third bands are derived from the $p$ orbitals with orbital
angular momentum $\ell_{z}=\pm 1$. The magnetic field couples to the orbital
angular momentum leading to a linear splitting of the two bands with $B$. The
lowest energy state after the $s$ orbital bands have been filled is the state
with orbital angular momentum $\ell_{z}=+1$ and spin aligned with $B$. In the
experimental phase diagram, we expect this state to be filled in the region
$-3\leq\nu\leq-2$ for $B<6$ T, and $-2\leq\nu\leq-1$ for $B>6$ T due to the
interaction induced band reordering. In both these regions, only electron-like
Hofstadter gaps with positive slope are observed, i.e. $(t,s)$ with positive
$t$. This is consistent with the Hofstadter spectrum in Fig. S12 which shows
that only electron-like Landau level gaps in the second band are appreciable.
This can be understood as arising to an asymmetry in the dispersion of the
band: the band features a flat band top but dispersive band bottom. Thus, only
Landau levels emanating from the band bottom, i.e. electron-like Landau
levels, are observable.

For filling beyond this, we expect the bands to be strongly renormalized by
interactions with filled states and no longer resemble that of the single
particle band structure. The order in which the bands fill will also
sensitively depend on interactions and spin $g$-factor. The closing and
reopening of gaps at fixed filling $\nu=-3$ and $\\-2.5$ are likely also
related to changes in the underlying band occupation. A comprehensive
theoretical study of the interaction-renormalized band structure at higher
filling, and the corresponding Hofstadter spectrum is beyond the scope of this
work.

##  IX Gaps of the Hofstadter states

Figure S13: Hofstadter gaps for $-4<\nu<-2$. (A,C) Zoom-ins from Figs. 1B and
4A in the main text, which were measured at two independent locations. (B,D)
Gaps extracted from (A) and (C), respectively. (A) was taken at temperature
$T=1.6$ K while (C) at $T=330$ mK and with a smaller tip-sample separation.

We extract the Hofstadter gaps from measurements in two different locations
(Fig. S13A,C) and plot them as Wannier diagrams in Fig. S13B,D. Both plots
show similar quantitative and qualitative trends with magnetic field. The
Hofstadter gaps first increase with field after emerging but then decrease
rapidly, both near the region of reentrant charge order and also when
approaching $\Phi/\Phi_{0}=1/2$. Additionally, the Wannier plot further
demonstrates the asymmetry in the second moiré band dispersion. In particular,
above the reentrant charge order for $-3<\nu<-2$, we observe larger electron-
like than hole-like Hofstadter gaps in agreement with the different effective
masses for top and bottom band edges (this is also evident from the size of
the Hofstadter gaps for $-2<\nu<-1$ at high fields).

## References

  * Wang  _et al._ (2013) L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. A. Muller, J. Guo, P. Kim, J. Hone, K. L. Shepard, and C. R. Dean, Science 342, 614 (2013), publisher: American Association for the Advancement of Science Section: Report.
  * Movva  _et al._ (2015) H. C. P. Movva, A. Rai, S. Kang, K. Kim, B. Fallahazad, T. Taniguchi, K. Watanabe, E. Tutuc, and S. K. Banerjee, ACS Nano 9, 10402 (2015).
  * Kumar  _et al._ (2013) N. Kumar, S. Najmaei, Q. Cui, F. Ceballos, P. M. Ajayan, J. Lou, and H. Zhao, Physical Review B 87, 161403 (2013).
  * Li  _et al._ (2013) Y. Li, Y. Rao, K. F. Mak, Y. You, S. Wang, C. R. Dean, and T. F. Heinz, Nano Letters 13, 3329 (2013).
  * Eisenstein  _et al._ (1994) J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Physical Review B 50, 1760 (1994).
  * Pan  _et al._ (2020) H. Pan, F. Wu, and S. Das Sarma, Physical Review B 102, 201104 (2020).
  * Morales-Durán  _et al._ (2022) N. Morales-Durán, P. Potasz, and A. H. MacDonald, “Magnetism and Quantum Melting in Moir$\backslash$’e-Material Wigner Crystals,”  (2022), arXiv:2210.15168 [cond-mat] .
  * 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).
  * Li  _et al._ (2021a) T. Li, J. Zhu, Y. Tang, K. Watanabe, T. Taniguchi, V. Elser, J. Shan, and K. F. Mak, Nature Nanotechnology 16, 1068 (2021a).
  * Jin  _et al._ (2021) C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, L. Fu, J. Shan, and K. F. Mak, Nature Materials 20, 940 (2021), number: 7 Publisher: Nature Publishing Group.
  * Li  _et al._ (2021b) T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, J. Shan, and K. F. Mak, Nature 597, 350 (2021b).
  * Xu  _et al._ (2020) Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan, Nature 587, 214 (2020).
  * Huang  _et al._ (2021) X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, S.-F. Shi, and Y.-T. Cui, Nature Physics 17, 715 (2021).
  * Wu  _et al._ (2018a) F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Physical Review Letters 121, 026402 (2018a).
  * Peng  _et al._ (2016) H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Physical Review X 6, 041005 (2016).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical review B 54, 11169 (1996).
  * Wu  _et al._ (2018b) F. Wu, T. Lovorn, E. Tutuc, and A. MacDonald, Physical Review Letters 121, 026402 (2018b).
  * Zhang  _et al._ (2020) Y. Zhang, N. F. Q. Yuan, and L. Fu, Physical Review B 102, 201115 (2020).
  * Hu and MacDonald (2021) N. C. Hu and A. H. MacDonald, Physical Review B 104, 214403 (2021).
  * Wilson  _et al._ (2017) N. R. Wilson, P. V. Nguyen, K. Seyler, P. Rivera, A. J. Marsden, Z. P. L. Laker, G. C. Constantinescu, V. Kandyba, A. Barinov, N. D. M. Hine, X. Xu, and D. H. Cobden, Science Advances 3, e1601832 (2017), publisher: American Association for the Advancement of Science.
  * Seyler  _et al._ (2019) K. L. Seyler, P. Rivera, H. Yu, N. P. Wilson, E. L. Ray, D. G. Mandrus, J. Yan, W. Yao, and X. Xu, Nature 567, 66 (2019), number: 7746 Publisher: Nature Publishing Group.
  * Tran  _et al._ (2019) K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, J. Embley, A. Zepeda, M. Campbell, T. Autry, T. Taniguchi, K. Watanabe, N. Lu, S. K. Banerjee, K. L. Silverman, S. Kim, E. Tutuc, L. Yang, A. H. MacDonald, and X. Li, Nature 567, 71 (2019), number: 7746 Publisher: Nature Publishing Group.
  * Barré  _et al._ (2022) E. Barré, O. Karni, E. Liu, A. L. O’Beirne, X. Chen, H. B. Ribeiro, L. Yu, B. Kim, K. Watanabe, T. Taniguchi, K. Barmak, C. H. Lui, S. Refaely-Abramson, F. H. da Jornada, and T. F. Heinz, Science 376, 406 (2022), publisher: American Association for the Advancement of Science.
  * Zeng and MacDonald (2022) Y. Zeng and A. H. MacDonald, Physical Review B 106, 035115 (2022).
  * Wagner  _et al._ (1992) M. Wagner, U. Merkt, and A. Chaplik, Physical Review B 45, 1951 (1992).
  * Ashoori  _et al._ (1993) R. Ashoori, H. Stormer, J. Weiner, L. Pfeiffer, K. Baldwin, and K. West, Physical review letters 71, 613 (1993).