# Competing Generalized Wigner Crystal States in Moiré Heterostructures

Shu Fay Ung<EMAIL_ADDRESS>Department of Chemistry, Columbia University,
New York, NY, United States  Joonho Lee<EMAIL_ADDRESS>Department of
Chemistry, Columbia University, New York, NY, United States  Department of
Chemistry and Chemical Biology, Harvard University, MA, United States  David
R. Reichman<EMAIL_ADDRESS>Department of Chemistry, Columbia
University, New York, NY, United States

###### Abstract

We present a comprehensive study of Wigner crystals across a broad range of
filling factors for system sizes up to 162 holes employing Hartree-Fock theory
and explicitly correlated wave function approaches. While we find broad
agreement with the behavior observed in experiments and classical Monte Carlo
simulations, we highlight the fact that the Hartree-Fock energy landscape
appears to be remarkably complex, with many competing states, both ordered and
disordered, separated by energies of a fraction of $\sim$ 1 meV/hole. We
explicitly demonstrate that these states are truly metastable by performing a
stability analysis at the Hartree-Fock level. Correlated wave function methods
reveal correlation energies that are very small, but are large enough to tip
the balance of state ordering found at the Hartree-Fock level.

#### Introduction.

The interest in a novel class of physical systems—so-called moiré systems—has
recently gained tremendous momentum in the condensed-matter physics and
materials science communities [1]. These systems are created from 2D layers
stacked with a small misalignment such that a characteristic long wavelength
moiré pattern forms. The associated moiré period parametrizes the kinetic and
interaction energy scales in the stacked system, providing a control knob to
alter the electronic Hamiltonian. In particular, the interaction strength
between charges can be tuned by changing the twist angle or the layer
constituents, rendering such systems a versatile platform for studying strong-
correlation physics. Experiments on moiré superlattices have revealed a rich
phase diagram of strongly-correlated electronic states, including Mott
insulators, generalized Wigner crystals (GWCs), stripe phases, and liquid
crystals by varying the interaction strength together with the electron
filling factor [2, 3, 4, 5, 6, 7].

Despite the large moiré unit cell, the low-energy physics can be described by
an emergent single-particle Hamiltonian (i.e. the continuum model) sharing the
periodicity of the moiré superlattice. We focus on transition metal
dichalcogenide (TMD) _heterobilayers_ , where carriers near the Fermi level
are localized in one layer only and experience a periodic moiré potential due
to the spatial variation in the valence-band maximum of the other layer [8].
Since the valence-band extrema located at the $K$ and $K^{\prime}$ momentum
points are decoupled and exhibit a large spin-splitting caused by spin-orbit
interactions, spin indices are locked with valley degrees of freedom, giving
rise to an effective two-fold spin degeneracy in the continuum model [9].

The accurate description of strongly-correlated states, however, requires the
treatment of many-body interactions. A natural approach is to incorporate
long-range Coulomb interactions into the continuum model description (i.e. the
_interacting_ continuum model) [10, 11, 12]. This framework resembles that of
the two-dimensional electron gas (2DEG) with a pinning moiré potential, which
may also harbor exotic physics arising from Coulomb frustration [13] and
glassy dynamics [14], in addition to displaying well-known 2DEG phenomena such
as standard Wigner crystallization [15, 16, 17].

Previous theoretical work on GWCs have employed Hartree-Fock (HF) [12, 18, 19,
20, 21], density functional theory [22], classical and quantum Monte Carlo
[23, 24, 25, 22], and exact diagonalization [26, 27] to probe the ground state
phase diagram. While classical Monte Carlo is reasonable deep within the
crystalline regime, a comprehensive picture requires a quantum treatment
beyond mean-field theory to capture strong correlation effects. Such studies
have largely focused on the magnetic behaviour at filling factors
$\nu=1/3,2/3,$ and 1.

In this Letter, we present a comprehensive quantum mechanical study of GWCs in
hole-doped TMD heterobilayers at ten distinct filling factors up to $\nu=1$.
We first characterize the quantum energy landscape of Hartree-Fock solutions
for systems up to 162 holes, demonstrating that multiple metastable states (as
defined by a positive-definite spectrum of the Hartree-Fock Hessian matrix)
generically exist. Such states are very close to each other in energy and may
be disordered or ordered. We then demonstrate that correlated wave function
approaches may be used to accurately estimate the correlation energies of
these states. Such methods have reasonable scaling with system size, and thus
may be used to study large moiré systems. We show that correlation energies
are very small, but can be sufficient to alter the energetic ordering of the
states revealed by Hartree-Fock theory. Hartree atomic units are used unless
stated otherwise.

#### Model and methods.

The low-energy physics of TMD heterobilayers is described by a continuum model
constructed from the topmost valence band [9]. Carriers near the Fermi level
are localized in one layer only, with the effect of the other layer appearing
through local variations in the band gap that correspond to moiré length-scale
variations in the local atomic registry [8]. This manifests in the carriers
experiencing a moiré potential $\Delta(\vec{r})$ that is periodic with respect
to the superlattice, giving rise to the single-particle moiré Hamiltonian [9]

$\displaystyle\hat{h}_{i}$
$\displaystyle=-\frac{\nabla_{i}^{2}}{2m^{*}}+\Delta(\vec{r}_{i}),$ (1)

where $m^{*}$ is the carrier effective mass. To the lowest order in a harmonic
expansion, we have

$\Delta(\vec{r})=2V_{m}\sum_{j=1}^{3}\cos{\left(\vec{g}_{j}\cdot\vec{r}+\phi\right)},$
(2)

where $V_{m}$ is the potential strength, $\phi$ characterizes the shape of
$\Delta(\vec{r})$ in the moiré unit cell, and $\vec{g}_{j}$ are the moiré
reciprocal lattice vectors in the first shell which satisfy $C_{3}$ symmetry.
The emergence of $\Delta(\vec{r})$ leads to zone-folding of the monolayer
Brillouin zones, yielding smaller moiré Brillouin zones (mBZ) that host flat
bands. $V_{m}$ and $\phi$ are intrinsic material properties that can be
obtained from fitting to _ab initio_ band structures [9, 10].

Taking the carriers to be holes, we set $m^{*}=-0.45m_{e}$ ($m_{e}$ is the
bare electron mass), $V_{m}=-15$ meV, and $\phi=45^{\circ}$ for WSe2/WS2. At
zero twist, the 4% lattice mismatch between the WSe2 and WS2 monolayers
produces a moiré superlattice with a lattice constant of $L_{m}=8.3$ nm [10].
Spin-valley locking [9] implies an effective two-fold degeneracy in the model,
thus we identify the filling factor $\nu=2$ (2 carriers per moiré unit cell)
as “full-filling”. We only examine $\nu\leq 1$ since interactions renormalize
bands more strongly with higher charge densities, leading to remote band
mixing that demands caution with respect to the single-band continuum model
for $\nu>1$ [27].

The interacting continuum model is obtained by adding Coulomb interactions

$\hat{H}=\sum_{i=1}^{N}\hat{h}_{i}+\frac{1}{2}\sum_{ij=1}^{N}\frac{1}{\epsilon|\vec{r}_{i}-\vec{r}_{j}|}.$
(3)

Here, we used a conventional $1/\epsilon$ screened Coulomb potential with
$\epsilon$ the effective dielectric constant, although the effects of gate-
screening [18, 27, 25] and a Rytova-Keldysh-type screening [28, 29] can also
be considered. See Appendix B in the Supplemental Material for a discussion.

Our study employed the generalized Hartree-Fock (GHF) approximation, where
spin-orbitals $\\{\chi_{i}(\vec{r})\\}$ comprising the Slater determinant
$\mathinner{|{\Psi_{\text{HF}}}\rangle}$ possess both spin-up and spin-down
character, i.e. they retain a 2-component spinor structure

$\chi_{i}(\vec{r})=\begin{bmatrix}\chi_{i}^{\uparrow}(\vec{r})\\\\[5.0pt]
\chi_{i}^{\downarrow}(\vec{r})\end{bmatrix},$ (4)

where $\chi_{i}^{\sigma}(\vec{r})$ are spatial orbitals in the spin-$\sigma$
sector. By variationally minimizing the Hartree-Fock energy
$E_{\text{HF}}=\mathinner{\langle{\Psi_{\text{HF}}}|}\hat{H}\mathinner{|{\Psi_{\text{HF}}}\rangle}$
with respect to the orbitals $\\{\chi_{i}\\}$, we obtain a set of integro-
differential equations that can be self-consistently solved for the optimal
$\\{\chi_{i}\\}$. We represented the $\chi_{i}^{\sigma}$ in a plane-wave basis
and performed supercell calculations at the $\Gamma$-point of the mBZ. By
imposing periodic boundary conditions only over the supercell, we can
naturally obtain solutions that are expected to break the superlattice
translational symmetry at fractional filling factors [3, 5, 7, 24, 25, 27]. At
each filling factor, we found at least one GHF solution starting from
solutions of the continuum model with matrix elements in the mixed-spin
sectors of the initial density matrix $\mathbf{P}$ (i.e.
$\mathbf{P}^{\uparrow\downarrow}$ and $\mathbf{P}^{\downarrow\uparrow}$) drawn
from a standard Gaussian distribution $\mathcal{N}(0,1)$, normalized by
$10^{-3}\times\max|\mathbf{P}|$. This “noisy” continuum model guess served to
unbias solutions away from collinearity. In addition, we selectively located
other GHF minima via physically motivated initial guesses. We tested the
stability of each GHF solution by computing the Hartree-Fock Hessian matrix,
as discussed in Appendix K of the Supplemental Material. A stable solution
corresponds to one with a positive-definite spectrum; states with one or more
negative eigenvalues correspond to saddle points of index 1 or higher on the
Hartree-Fock landscape of solutions. We considered systems with up to 162
holes and, where possible, performed finite size analyses as detailed in the
Supplemental Material.

To go beyond mean-field theory and incorporate electron correlation effects,
we leveraged several state-of-the-art approaches that have not found common
usage in the study of moiré systems. We studied systems smaller than the
largest ones for which GHF solutions are found, but which are too large to
study by exact diagonalization. We do this as a compromise between choosing
system sizes that are: (i) too small such that it is difficult to parse the
distinction from true correlation effects; and (ii) too large to have
reference solutions to compare with. To find (near) exact energy values from
our distinct GHF solutions, we employed the heat-bath configuration
interaction (HCI) approach [30, 31] and compared these values to those
obtained via the coupled-cluster singles-doubles (CCSD) and the more accurate
but expensive coupled-cluster singles-doubles plus perturbative triples
(CCSD(T)) approaches [32, 33, 34]. It should be noted that CCSD and perhaps
CCSD(T) scale in a manner that would enable their usage on the largest systems
we investigate in this work, making them attractive approaches for future
investigations of correlation effects in moiré systems in general [35].

#### HF results.

We first present a survey of the behavior of GHF solutions across $\nu$ for
the interacting continuum model. The system sizes we considered are such that
the GWCs are commensurate with the supercell, and at each $\nu$ satisfy
$28\leq N\leq 162$ (see Table 1, Appendix D in the Supplemental Material for a
table of the system sizes considered).

Figure 1: (a) Charge gaps $\Delta_{c}$ relative to the moiré scale interaction
strength $U_{m}=1/\epsilon L_{m}$ at various fractional filling factors $\nu$,
where $L_{m}$ is the moiré superlattice constant. The inset shows experimental
measurements of the 2s exciton resonance energy $E_{\text{2s}}$ used to probe
the charge gap in untwisted WSe2/WS2 heterobilayers [3]. All solutions
depicted were obtained using the continuum model guess as discussed in the
main text, except for the red cross which marks $\Delta_{c}/U$ calculated from
the ordered stripe solution. At the level of HF, we observe an asymmetry in
$\Delta_{c}$ about $\nu=1/2$ which qualitatively agrees with experiment and
previous theoretical work. (b)-(d) The charge density at select filling
factors plotted over the supercell. The indices $n_{1},n_{2}\in\mathbb{Z}^{+}$
represent a point $\vec{r}=n_{1}\vec{L}_{m,1}+n_{2}\vec{L}_{m,2}$ in the
supercell, where $\vec{L}_{m,i}$ are the primitive moiré lattice vectors.
Disordered configurations were generally observed across $\nu$, although
ordered phases were also seen at $\nu=1/7,1/3,2/3,$ and 1. At $\nu=1/2$, the
continuum model guess yielded a disordered phase in contrast to the ordered
phase obtained from the stripe guess (Fig. 2a).

Fig. 1a shows charge gaps $\Delta_{c}$ obtained either from the largest system
in a series of finite-sized calculations, or values extrapolated to the
thermodynamic limit where appropriate (see Appendix H in the Supplemental
Material for details). At each $\nu$, $\Delta_{c}$ of only ordered solutions
were considered if found; the charge gaps otherwise originate from disordered
solutions. Larger values of $\Delta_{c}$ can occur in ordered solutions
(indicated at $\nu=1/2$ in Fig. 1a), but require tuning of the initial guess
to find.

We plot $\Delta_{c}$ normalized by the moiré scale interaction strength
$U_{m}=1/\epsilon L_{m}$. As in previous experiments [3, 7] and theoretical
work [27], we observed an asymmetry in $\Delta_{c}$ about $\nu$ = 1/2 at the
level of Hartree-Fock theory. The charge gaps for $\nu<1/2$ are larger than
for $\nu>1/2$ due to the larger Coulomb-to-kinetic-energy ratio at lower
average charge densities, which is reminiscent of the low-density behaviour in
the 2DEG.

Compared to the exact diagonalization results of Morales-Durán et. al. [27]
obtained with small system sizes ($N\leq 16$ for $\nu\leq 1$), we found larger
charge gaps and more stable GWCs 111The discrepancy with the charge gaps
reported in [27] appears not to arise entirely due to the use of Hartree-Fock,
but due to differences in the model (i.e. whether a projection onto continuum
model orbitals was involved), the system size, the basis set size, and the
manner in which the charge gap was calculated (e.g. via the HOMO-LUMO gap or
the quasiparticle gap $E(N+1)+E(N-1)-2E(N)$) as well. For larger values of
$W/V_{m}$ (i.e. smaller $V_{m}$) we found a closing of the gap towards
metallic behavior compared to the values reported in Fig. 1. Our
$(L_{m},m^{*},V_{m})$ parameter values give $W_{m}/V_{m}=0.16$, where
$W_{m}=1/m^{*}L_{m}^{2}$ is the moiré scale kinetic energy. While they
predicted metallic phases across fractional filling factors at this value of
$W_{m}/V_{m}$, we found insulating states with larger $\Delta_{c}/U_{m}$
values than those computed in their insulating states.

We also observed a manifestation of finite-size effects in the spin-
collinearity of solutions. While we only found collinear polarized Hartree-
Fock solutions for $\nu<1$ in the largest $N$-particle systems studied at each
$\nu$ when initiating from the continuum model guess as described above, the
smallest $N$-particle system found may be non-collinear (see Appendix J in the
Supplemental Material for details on spin-collinearity determination). This
suggests that the determination of the spin properties of low-lying states may
be sensitive to finite-size effects, at least at the Hartree-Fock level. It
should be noted that we have occasionally found our initial guess to be
insufficient in biasing solutions away from ferromagnetic behavior, and that
lower energy solutions may indeed be non-collinear. As the primary goal of
this work is to examine charge orders, we do not believe any of our main
conclusions are affected.

Using initial guesses from the continuum model (as described in the previous
section), the GHF solutions we found consist of localized charge densities
across all $\nu$ which are generally of a disordered nature. Figs. 1b-d only
show the densities for $\nu=2/5,1/2,$ and $2/3$, while Appendix E in the
Supplemental Material provides the remaining real- and Fourier-space figures
as well as further discussion. In the cases shown in Fig. 1, $C_{3}$ symmetry
is generally preserved in Fourier space, implying that the charge densities
also obey the $C_{3}$ symmetry of the moiré superlattice (with $\nu=2/5$
perhaps the exception). The superlattice translational symmetry, on the
contrary, can be broken. The ordered phases we observed occur at
$\nu=1/7,1/3,2/3,$ and 1, which exhibit triangular lattice and honeycomb
charge orders (Fig. 1d) as seen in experiment [5] and previous theoretical
work [37, 7, 27]. Disordered states at other filling factors also bear
resemblance to MC predictions, e.g. the “columnar dimer crystal” at $\nu=2/5$
(Fig. 1b) [25]. Notably, our survey revealed a disordered (labyrinthine)
stripe configuration at $\nu=1/2$, in contrast to the ordered stripe phase
often discussed in the literature 222A similar labyrinthine solution was also
obtained via classical Monte Carlo simulations in the Extended Data of [7].
The fragility of the ordered stripe state has been noted experimentally by Li
et. al. [5], who highlight the sensitivity of this phase to lattice strain and
local inhomogeneities. In the absence of such experimental considerations,
however, we posit that our disordered solution is one of several nearly
degenerate local minima (i.e. _stable_ solutions) that differ slightly in
their real-space configurations. As discussed in more detail below, we can
also find low-lying ordered states at filling fractions where the continuum
model guess yields disordered configurations. This is accomplished via
selection of physically motivated initial guesses which are also stable.

Figure 2: The charge density $\rho(\vec{r})$ at filling factor $\nu=1/2$.
(a)-(c) HF solutions displaying an ordered stripe state, disordered stripe
(labyrinthine) state, and an inhomogeneous honeycomb structure. (d)-(f) CCSD
solutions. Comparing the top and bottom rows, we observed that quantum
fluctuations induced by correlation delocalizes $\rho(\vec{r})$. Note that
$\rho(\vec{r})$ in (a) and (d) are shifted by $L_{m}$ relative to one another,
while (e) and (f) are rotated by $60^{\circ}$ relative to each other.

#### Nearly degenerate solutions.

We now focus on $\nu=1/2$ to probe the existence of nearly degenerate yet
distinct HF solutions by supplying different initial guesses. In an attempt to
recover the ordered stripe state, we explicitly constructed a polarized
(ordered) stripe guess in addition to the continuum model guess employed
previously. The results presented here were obtained with $N=18$ and
$N_{\mathrm{basis}}=211$; larger system and basis set size studies can be
found in Appendix F of the Supplemental Material.

Figs. 2a-c show the charge configurations obtained from each initialization.
Unsurprisingly, the solution from the ordered stripe guess is also an ordered
stripe configuration (Fig. 2a), whereas the continuum model guess yielded a
disordered (labyrinthine) stripe state (Fig. 2b). While several other HF
solutions were also discovered throughout our investigation, we here only
present one additional solution that displays an inhomogeneous density across
the honeycomb lattice (Fig. 2c). These solutions can also be clearly
distinguished in Fourier space (Figs. 3a-c)—notably, the ordered stripe phase
exhibits peaks that obey $C_{2}$ symmetry instead of the $C_{3}$ symmetry of
the moiré superlattice. The energy per particle at the GHF level for these
solutions are $\varepsilon_{\mathrm{ordered}}=-41.233$ meV,
$\varepsilon_{\mathrm{disordered}}=-41.465$ meV, and
$\varepsilon_{\mathrm{honeycomb}}=-41.114$ meV, respectively. Taken together,
these results suggest the existence of a rough energy landscape with many
nearly degenerate local minima surrounding the ordered phase.

Figure 3: Fourier-space plots of the charge densities obtained at $\nu=1/2$.
(a)-(f) Fast Fourier transform (FFT) of the corresponding densities in Fig. 2
with the average subtracted, i.e.
$|\rho(\vec{k})|=|\mathrm{FFT}\left\\{\rho(\vec{r})-\langle\rho\rangle\right\\}|$.
The indices $\kappa_{1},\kappa_{2}\in\mathbb{Z}$ represent a point
$\vec{k}=\kappa_{1}\vec{G}_{1}+\kappa_{2}\vec{G}_{2}$ in reciprocal space,
where $\vec{G}_{i}$ are the primitive reciprocal lattice vectors associated
with the supercell. (g)-(i) FFT peaks extracted from (a)-(c), plotted directly
in reciprocal space $\vec{k}=(k_{1},k_{2})$. Large dots denote where the FFT
signal is maximum.

#### Further study with correlated methods.

To investigate the effect of quantum fluctuations at $\nu=1/2$, we performed
correlated calculations starting from the mean-field solutions found
previously. Again considering the case $N=18$ and $N_{\mathrm{basis}}=211$, we
leveraged only spin-collinear variants of HCI, CCSD, and CCSD(T) since the HF
solutions we investigated were spin-collinear. We will use “vHCI” and “vHCI +
PT2” to distinguish between the variational component of HCI and the full
algorithm (variational plus perturbative components). A brief discussion of
these approaches can be found in Appendix C of the Supplemental Material.

Figure 4: Energy per particle $\varepsilon$ relative to the HF ordered stripe
solution, denoted with $\varepsilon_{\mathrm{ref}}$. We use vHCI to denote
only the variational component of HCI, while vHCI + PT2 denotes both the
variational and perturbative components, i.e. the full HCI algorithm. The
energy differences are $<1$ meV across all methods. Whereas the energetic
ordering found by HF and all correlated approaches are inconsistent, all
correlated approaches yield orderings and energy values in remarkable
agreement with each other.

Our results show that the complexity of the energy landscape persists even
after introducing electron correlation. Fig. 4 illustrates the small energy
scales associated with the solutions. For reference, the relevant energy
scales defining our system are $U_{m}\sim V_{m}\sim\mathcal{O}(10)$ meV and
$W_{m}\sim\mathcal{O}(1)$ meV. Both the correlation energy and energy
differences between the three solutions are $\lesssim 1$ meV. Similar energies
obtained across vHCI, vHCI+PT2, CCSD, and CCSD(T) moreover provide confidence
that these values are nearly exact. The three correlated solutions remain
distinct, which suggests that the metastability associated with the HF
solutions is retained. Each correlated approach qualitatively predicts charge
densities that agree with each other and HF—panels (d)-(f) of Figs. 2 and 3
demonstrate the delocalization of charge density induced by quantum
fluctuations which still preserve the essential Fourier-space structure from
HF.

Nevertheless, the corrections produced by charge correlation is non-trivial.
The minuscule correlation energy is sufficient to alter the energetic ordering
of states from HF, despite indicating that a single Slater determinant is
adequate in providing a quantitative description of each state. For example,
while HF predicts the lowest-energy state of the three we have isolated to be
the disordered stripe configuration, all correlated methods predict it to be
the honeycomb solution instead.

#### Conclusions.

In this Letter we have numerically studied GWCs in hole-doped TMD
heterobilayers at various fractional filling factors using HF and correlated
methods. As a concrete example, we focused on the untwisted WSe2/WS2
heterobilayer embedded in a hBN dielectric environment, which is an
experimentally popular system for probing GWCs. Our results show qualitative
agreement with experiment in the asymmetry of charge gaps about $\nu=1/2$ and
in the charge orderings observed, particularly at $\nu=1/3,2/3,$ and 1. The
multiple and potentially disordered solutions we discovered suggested a rugged
energy manifold of nearly degenerate states, thus we focused on $\nu=1/2$ to
probe this. Indeed, the predicted energy differences between solutions were
$<1$ meV at the level of HF and beyond. Utilizing HCI, CCSD, and CCSD(T), we
showed that while quantum fluctuations delocalize the charge densities
produced by mean field theory, the densities still bear strong similarities
across the methods. These observations reveal the complexity of the moiré
energy landscape—defined by several nearly degenerate metastable states—in
both mean-field and correlated theories. Future work will be devoted to a
systematic and complete characterization of the moiré energy landscape as well
as addressing its implications for experiments.

The Supplemental Material for this work is available online. We acknowledge
funding support from NSF CHE–2245592 and thank Y. Zhang for help concerning
the implementation of the model and A. Millis for useful discussions.
Computations in this paper were primarily executed on the Ginsburg High
Performance Computing Cluster at Columbia University. Some calculations were
run on the FASRC Cannon cluster supported by the FAS Division of Science
Research Computing Group at Harvard University. This work also used Bridges-2
Regular Memory (RM) nodes at the Pittsburgh Supercomputing Center and the
Expanse cluster at the San Diego Supercomputer Center through allocation
CHE220021 from the Advanced Cyberinfrastructure Coordination Ecosystem:
Services & Support (ACCESS) program, which is supported by National Science
Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296.

## References

  * Andrei  _et al._ [2021] E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, Nature Reviews Materials 6, 201 (2021).
  * 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).
  * 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).
  * 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).
  * Li  _et al._ [2021] H. Li, S. Li, E. C. Regan, D. Wang, W. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, M. F. Crommie, and F. Wang, Nature 597, 650 (2021).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. Iqbal 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).
  * 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).
  * Zhang  _et al._ [2017] C. Zhang, C.-P. Chuu, X. Ren, M.-Y. Li, L.-J. Li, C. Jin, M.-Y. Chou, and C.-K. Shih, Science Advances 3, e1601459 (2017).
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. Macdonald, Physical Review Letters 121 (2018), 10.1103/PhysRevLett.121.026402, arXiv: 1804.03151.
  * Zhang  _et al._ [2020a] Y. Zhang, N. F. Yuan, and L. Fu, Physical Review B 102, 201115 (2020a), arXiv: 1910.14061.
  * Zhang  _et al._ [2020b] Y. Zhang, H. Isobe, and L. Fu,  (2020b), arXiv: 2005.04238.
  * Hu and MacDonald [2021] N. C. Hu and A. H. MacDonald, Physical Review B 104, 214403 (2021).
  * Jamei  _et al._ [2005] R. Jamei, S. Kivelson, and B. Spivak, Physical Review Letters 94, 056805 (2005).
  * Schmalian and Wolynes [2000] J. Schmalian and P. G. Wolynes, Physical Review Letters 85, 836 (2000).
  * Tanatar and Ceperley [1989] B. Tanatar and D. M. Ceperley, Physical Review B 39, 5005 (1989).
  * Drummond and Needs [2009] N. D. Drummond and R. J. Needs, Physical Review Letters 102, 1 (2009).
  * Loos and Gill [2016] P.-F. Loos and P. M. W. Gill, WIREs Comput Mol Sci 6, 410 (2016).
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. Das Sarma, Physical Review B 102, 201104 (2020).
  * Pan and Das Sarma [2021] H. Pan and S. Das Sarma, Physical Review Letters 127, 096802 (2021).
  * Pan and Das Sarma [2022] H. Pan and S. Das Sarma, Physical Review B 105, 041109 (2022).
  * Kaushal  _et al._ [2022] N. Kaushal, N. Morales-Durán, A. H. MacDonald, and E. Dagotto, Communications Physics 5, 1 (2022).
  * Yang  _et al._ [2023] Y. Yang, M. Morales, and S. Zhang,  (2023), arXiv:2306.14954 [cond-mat].
  * Padhi  _et al._ [2021] B. Padhi, R. Chitra, and P. W. Phillips, Physical Review B 103, 125146 (2021).
  * Zhang  _et al._ [2021] Y. Zhang, T. Liu, and L. Fu, Physical Review B 103, 155142 (2021).
  * Matty and Kim [2022] M. Matty and E.-A. Kim, Nature Communications 13, 7098 (2022).
  * Morales-Durán  _et al._ [2021] N. Morales-Durán, P. Potasz, and A. H. MacDonald, Physical Review B 103, L241110 (2021), arXiv:2011.13558 [cond-mat].
  * Morales-Durán  _et al._ [2022] N. Morales-Durán, P. Potasz, and A. H. MacDonald,  (2022), arXiv:2210.15168 [cond-mat].
  * Rytova [1967] N. S. Rytova,  (1967), arXiv:1806.00976 [cond-mat].
  * Keldysh [1979] L. V. Keldysh, JETP Lett. 29 (1979).
  * Holmes  _et al._ [2016] A. A. Holmes, N. M. Tubman, and C. J. Umrigar, Journal of Chemical Theory and Computation 12, 3674 (2016).
  * Sharma  _et al._ [2017] S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi, and C. J. Umrigar, Journal of Chemical Theory and Computation 13, 1595 (2017).
  * Coester and Kümmel [1960] F. Coester and H. Kümmel, Nuclear Physics 17, 477 (1960).
  * Čížek [1966] J. Čížek, The Journal of Chemical Physics 45, 4256 (1966).
  * Čižek and Paldus [1971] J. Čižek and J. Paldus, International Journal of Quantum Chemistry 5, 359 (1971).
  * Faulstich  _et al._ [2022] F. M. Faulstich, K. D. Stubbs, Q. Zhu, T. Soejima, R. Dilip, H. Zhai, R. Kim, M. P. Zaletel, G. K.-L. Chan, and L. Lin,  (2022), arXiv:2211.09243 [cond-mat, physics:physics].
  * Note [1] The discrepancy with the charge gaps reported in [27] appears not to arise entirely due to the use of Hartree-Fock, but due to differences in the model (i.e. whether a projection onto continuum model orbitals was involved), the system size, the basis set size, and the manner in which the charge gap was calculated (e.g. via the HOMO-LUMO gap or the quasiparticle gap $E(N+1)+E(N-1)-2E(N)$) as well. For larger values of $W/V_{m}$ (i.e. smaller $V_{m}$) we found a closing of the gap towards metallic behavior compared to the values reported in Fig. 1.
  * Liu  _et al._ [2021] E. Liu, T. Taniguchi, K. Watanabe, N. M. Gabor, Y.-T. Cui, and C. H. Lui, Physical Review Letters 127, 037402 (2021).
  * Note [2] A similar labyrinthine solution was also obtained via classical Monte Carlo simulations in the Extended Data of [7].
  * Bernevig  _et al._ [2021] B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Physical Review B 103, 205413 (2021).
  * Szabo and Ostlund [1996] A. Szabo and N. L. Ostlund, “Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory,” (1996).
  * Shavitt and Bartlett [2009] I. Shavitt and R. J. Bartlett,  _Many-body methods in chemistry and physics: MBPT and coupled-cluster theory_ , Cambridge molecular science (Cambridge University Press, Cambridge ; New York, 2009).
  * Goings  _et al._ [2015] J. J. Goings, F. Ding, M. J. Frisch, and X. Li, The Journal of Chemical Physics 142, 154109 (2015).
  * Epifanovsky  _et al._ [2021] E. Epifanovsky _et al._ , The Journal of Chemical Physics 155, 084801 (2021).
  * Pulay [1980] P. Pulay, Chemical Physics Letters 73, 393 (1980).
  * Pulay [1982] P. Pulay, Journal of Computational Chemistry 3, 556 (1982).
  * Van Voorhis and Head-Gordon [2002] T. Van Voorhis and M. Head-Gordon, Molecular Physics 100, 1713 (2002).
  * Sun [2015] Q. Sun, Journal of Computational Chemistry 36, 1664 (2015).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Computational Molecular Science 8, e1340 (2018).
  * Sun  _et al._ [2020] Q. Sun _et al._ , The Journal of Chemical Physics 153, 024109 (2020).
  * Smith  _et al._ [2017] J. E. T. Smith, B. Mussard, A. A. Holmes, and S. Sharma, Journal of Chemical Theory and Computation 13, 5468 (2017).
  * Note [3] In the sense that swapping the occupied and vacant sites in the $\nu=1/7$ configuration gives the $\nu=6/7$ configuration.
  * Xing  _et al._ [2021] X. Xing, X. Li, and L. Lin,  (2021), arXiv:2108.00206 [physics].
  * Bonsall and Maradudin [1977] L. Bonsall and A. A. Maradudin, Physical Review B 15, 1959 (1977).
  * Yang  _et al._ [2020] Y. Yang, V. Gorelov, C. Pierleoni, D. M. Ceperley, and M. Holzmann, Physical Review B 101, 085115 (2020).
  * Hunt  _et al._ [2018] R. J. Hunt, M. Szyniszewski, G. I. Prayogo, R. Maezono, and N. D. Drummond, Physical Review B 98, 075122 (2018).
  * Small  _et al._ [2015] D. W. Small, E. J. Sundstrom, and M. Head-Gordon, The Journal of Chemical Physics 142, 094112 (2015).
  * Sharada  _et al._ [2015] S. M. Sharada, D. Stück, E. J. Sundstrom, A. T. Bell, and M. Head-Gordon, Molecular Physics 113, 1802 (2015).

## Supplemental Material

##  Appendix A Model details

The interacting continuum model is described by the Hamiltonian

$\displaystyle\hat{H}=\sum_{i=1}^{N}\left\\{-\frac{\nabla_{i}^{2}}{2m^{*}}+\Delta(\vec{r}_{i})\right\\}+\frac{1}{2}\sum_{ij=1}^{N}\frac{1}{\epsilon|\vec{r}_{i}-\vec{r}_{j}|},$
(5)

where $m^{*}$ is the carrier effective mass, $\Delta$ is the moiré potential,
and $\epsilon$ is the effective dielectric constant. In our calculations, we
instead work with a scaled Hamiltonian $\hat{H}^{\prime}$ defined by

$\displaystyle\hat{H}^{\prime}$
$\displaystyle=\sum_{i=1}^{N}\left\\{-\frac{\nabla_{i}^{\prime
2}}{2}+\Delta^{\prime}(\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}_{i})\right\\}+\frac{1}{2}\sum_{ij=1}^{N}\frac{1}{|\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}_{i}-\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}_{j}|},$ (6) $\displaystyle\hat{H}^{\prime}$
$\displaystyle=\alpha\hat{H},\qquad\Delta^{\prime}=\alpha\Delta,\qquad\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}=\beta\vec{r},$ (7) $\displaystyle\alpha$
$\displaystyle=\frac{\epsilon^{2}}{m^{*}},\qquad\beta=\frac{m^{*}}{\epsilon},$
(8)

which effectively maps the interacting continuum model to a 2-dimensional
electron gas (2DEG) in the scaled moiré potential $\Delta^{\prime}$ [11]. This
form allows the interacting continuum model to be treated with density
functional theory (DFT) employing exchange-correlation functionals derived
from the 2DEG. Note that if $\mathinner{|{\Psi^{\prime}}\rangle}$ is an
eigenstate of $\hat{H}^{\prime}$ with energy $E^{\prime}$, then

$\displaystyle\hat{H}^{\prime}\mathinner{|{\Psi^{\prime}}\rangle}=\alpha\hat{H}\mathinner{|{\Psi^{\prime}}\rangle}$
$\displaystyle=E^{\prime}\mathinner{|{\Psi^{\prime}}\rangle},$
$\displaystyle\implies\hat{H}\mathinner{|{\Psi^{\prime}}\rangle}$
$\displaystyle=\frac{E^{\prime}}{\alpha}\mathinner{|{\Psi^{\prime}}\rangle},$
(9)

i.e. $\mathinner{|{\Psi^{\prime}}\rangle}$ is also an eigenstate of the
original Hamiltonian $\hat{H}$ with energy $E=E^{\prime}/\alpha$. Thus, we
divided energies obtained from the scaled Hamiltonian by $\alpha$ to recover
energies associated with the physical moiré system.

##  Appendix B The screened Coulomb potential

Figure 5: The setup for deriving the screened Coulomb potential.

In experiments on TMD moiré systems that we consider, the setup often involves
encapsulating the moiré system between dielectric materials and metallic gates
(Fig. 5). The encapsulating material and gates alter the electrostatic
environment for a charge in the TMD layer, resulting in a screened Coulomb
interaction between in-plane charges. Here, we derive a screened Coulomb
potential that is a generalization of the double gate-screened [27, 35, 39]
and Rytova-Keldysh [28, 29] potentials considered in moiré and TMD materials,
respectively.

Regarding the setup in Fig. 5, our task is to find the electric potential
$\varphi_{2}(\vec{r};\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})$ at a point
$\vec{r}=(x,y,z)$ due to a charge $q$ located at $\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}=(0,0,z^{\prime})$, where
$z,z^{\prime}\in[-l/2,l/2]$. Let $\varphi_{i}$ and $\epsilon_{i}$ denote the
potential and dielectric constant in layer $i$, respectively. Specializing to
the case where $\epsilon_{1}=\epsilon_{3}$, the electrostatic equations to
solve are

$\displaystyle\nabla_{\vec{r}}^{2}\ \varphi_{1}(\vec{r};\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$ $\displaystyle=0,$ (10)
$\displaystyle\nabla_{\vec{r}}^{2}\ \varphi_{2}(\vec{r};\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$ $\displaystyle=-\frac{4\pi
q}{\epsilon_{2}}\delta(\vec{r}-\vec{r}\mkern 2.0mu\vphantom{r}^{\prime}),$
(11) $\displaystyle\nabla_{\vec{r}}^{2}\ \varphi_{3}(\vec{r};\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$ $\displaystyle=0,$ (12)

with boundary conditions (BCs; dropping the $\vec{r}\mkern
2.0mu\vphantom{r}^{\prime}$ parameter for notational clarity)

$\displaystyle\varphi_{1}\left(\vec{\rho},\frac{d}{2}\right)$
$\displaystyle=0,$ (13)
$\displaystyle\varphi_{1}\left(\vec{\rho},\frac{l}{2}\right)$
$\displaystyle=\varphi_{2}\left(\vec{\rho},\frac{l}{2}\right),$ (14)
$\displaystyle\varphi_{2}\left(\vec{\rho},-\frac{l}{2}\right)$
$\displaystyle=\varphi_{3}\left(\vec{\rho},-\frac{l}{2}\right),$ (15)
$\displaystyle\varphi_{3}\left(\vec{\rho},-\frac{d}{2}\right)$
$\displaystyle=0,$ (16)
$\displaystyle\epsilon_{1}\frac{\partial\varphi_{1}}{\partial
z}\left(\vec{\rho},\frac{l}{2}\right)$
$\displaystyle=\epsilon_{2}\frac{\partial\varphi_{2}}{\partial
z}\left(\vec{\rho},\frac{l}{2}\right),$ (17)
$\displaystyle\epsilon_{2}\frac{\partial\varphi_{2}}{\partial
z}\left(\vec{\rho},-\frac{l}{2}\right)$
$\displaystyle=\epsilon_{1}\frac{\partial\varphi_{3}}{\partial
z}\left(\vec{\rho},-\frac{l}{2}\right).$ (18)

$\vec{\rho}\in\mathbb{R}^{2}$ is a vector in the $xy$-plane. Furthermore, we
assume $d\ll l$ so we can treat the layers as possessing infinite extent in
the $xy$-plane. Fourier transforming (10)-(12) in $\vec{\rho}$, we obtain a
system of second-order partial differential equations (PDEs)

$\displaystyle\frac{\partial^{2}}{\partial
z^{2}}\varphi_{1}(\vec{k},z)-k^{2}\varphi_{1}(\vec{k},z)$ $\displaystyle=0,$
(19) $\displaystyle\frac{\partial^{2}}{\partial
z^{2}}\varphi_{2}(\vec{k},z)-k^{2}\varphi_{2}(\vec{k},z)$
$\displaystyle=-\frac{4\pi q}{\epsilon_{2}}\delta(z-z^{\prime}),$ (20)
$\displaystyle\frac{\partial^{2}}{\partial
z^{2}}\varphi_{3}(\vec{k},z)-k^{2}\varphi_{3}(\vec{k},z)$ $\displaystyle=0.$
(21)

Let us define the operator $\mathcal{L}=\frac{\partial^{2}}{\partial
z^{2}}-k^{2}$. Since $\varphi_{2}(\vec{k},z)$ is the solution to an
_inhomogeneous_ PDE with _inhomogeneous_ BCs, we can write it as

$\displaystyle\varphi_{2}(\vec{k},z)$
$\displaystyle=\tilde{\varphi}_{2}(\vec{k},z)+G(\vec{k},z),$ (22)

where $\tilde{\varphi}_{2}(\vec{k},z)$ is the solution to the _homogeneous_
PDE $\mathcal{L}\varphi(\vec{k},z)=0$ with _inhomogeneous_ BCs, while
$G(\vec{k},z)$ is the solution to the _inhomogeneous_ PDE
$\mathcal{L}\varphi(\vec{k},z)=-(4\pi q/\epsilon_{2})\delta(z-z^{\prime})$
with _homogeneous_ BCs. The screened potential we obtain from solving the
system of PDEs is given by

$\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$
$\displaystyle=\Xi\Bigg{\\{}\left[e^{k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}+e^{-k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)-\sinh{[k(d-l)]}\right]$
$\displaystyle\qquad+\left[e^{k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}+e^{-k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)+\sinh{[k(d-l)]}\right]\Bigg{\\}},$
(23) $\displaystyle\Xi$ $\displaystyle=-\frac{4\pi
q}{\epsilon_{1}k\sinh{(kl)}}\left\\{\left(\frac{\epsilon_{2}}{\epsilon_{1}}+1\right)^{2}\sinh{(kd)}-\left(\frac{\epsilon_{2}}{\epsilon_{1}}-1\right)^{2}\sinh{[k(d-2l)]}-2\left[\left(\frac{\epsilon_{2}}{\epsilon_{1}}\right)^{2}-1\right]\sinh{(kl)}\right\\}^{-1},$
(24) $\displaystyle G(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})$
$\displaystyle=-\frac{2\pi
q}{\epsilon_{2}k\sinh{(kl)}}\left\\{\cosh{[k(z+z^{\prime})]}-\cosh{[k\left(|z-z^{\prime}|-l\right)]}\right\\}.$
(25)

###  B.1 Double gate-screened potential

In the limit $\epsilon_{1}=\epsilon_{2}$, the BCs match those of the double
gate-screened potential and so we recover this solution. Setting
$\epsilon_{1}=\epsilon_{2}$ in (23)-(25),

$\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$
$\displaystyle=\Xi\Bigg{\\{}\left[e^{k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}+e^{-k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}\right]\left[1-\cosh{[k(d-l)]}-\sinh{[k(d-l)]}\right]$
$\displaystyle\qquad+\left[e^{k(z-l/2)}\sinh{\left[k\left(\frac{l}{2}-z^{\prime}\right)\right]}+e^{-k(z+l/2)}\sinh{\left[k\left(\frac{l}{2}+z^{\prime}\right)\right]}\right]\left[1-\cosh{[k(d-l)]}+\sinh{[k(d-l)]}\right]\Bigg{\\}},$
(26) $\displaystyle\Xi$ $\displaystyle=-\frac{4\pi q}{\epsilon
k\sinh{(kl)}}\left\\{4\sinh{(kd)}\right\\}^{-1},$ (27) $\displaystyle
G(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})$
$\displaystyle=-\frac{2\pi q}{\epsilon
k\sinh{(kl)}}\left\\{\cosh{[k(z+z^{\prime})]}-\cosh{[k\left(|z-z^{\prime}|-l\right)]}\right\\}.$
(28)

Since Coulomb interactions in the GWC states occur over the length scale of
the moiré lattice constant $L_{m}\gg l$, we are interested only in the
asymptotic behaviour of $\varphi_{2}(\vec{r};\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$ for distances $\rho\gg l$. Using the fact that
$k\propto 1/\rho\ll 1/l$ in this case, it is sufficient for us to know the
Fourier component $\varphi_{2}(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$ for $kl\ll 1$. It then follows that
$k|z|,k|z^{\prime}|\ll 1$ and $k|z-z^{\prime}|\approx 0$ since $-l/2\leq
z,z^{\prime}\leq l/2$. Taylor expanding $\tilde{\varphi}_{2}$ and $G$ to first
order in $kl$ and making the approximation $z\approx z^{\prime}$ (note that we
still have $k|z|\sim\mathcal{O}(kl)$), we obtain

$\displaystyle\tilde{\varphi}_{2}(\vec{k},z)$
$\displaystyle\approx\Xi\left\\{2kl(1-\cosh{kd})+(kl)^{2}\sinh{kd}-4(kz)^{2}\sinh{kd}\right\\}$
$\displaystyle=-\frac{\pi q}{\epsilon
k\sinh{kd}}\left\\{2(1-\cosh{kd})+kl\sinh{kd}-4\frac{(kz)^{2}}{kl}\sinh{kd}\right\\},$
(29) $\displaystyle G(\vec{k},z)$ $\displaystyle\approx-\frac{2\pi q}{\epsilon
k}\left[2\frac{(kz)^{2}}{kl}-\frac{1}{2}kl\right],$ (30)
$\displaystyle\implies\varphi_{2}(\vec{k},z)$ $\displaystyle\approx-\frac{\pi
q}{\epsilon
k\sinh{kd}}\left\\{2(1-\cosh{kd})+\cancel{kl\sinh{kd}}-\cancel{4\frac{(kz)^{2}}{kl}\sinh{kd}}+\cancel{4\frac{(kz)^{2}}{kl}\sinh{kd}}-\cancel{kl\sinh{kd}}\right\\}$
$\displaystyle=\frac{2\pi q}{\epsilon
k}\cdot\tanh{\left(\frac{kd}{2}\right)}.$ (31)

###  B.2 Rytova-Keldysh potential

In the limit $d\to\infty$, the BCs are similar to the setup in [28] and we
recover Eq. (5) therein. We first shift our coordinates $z\to z-l/2$ and
$z^{\prime}\to z^{\prime}-l/2$ in (23)-(25) to match the coordinate system of
[28]:

$\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$
$\displaystyle=\Xi\Bigg{\\{}\left[e^{kz}\sinh{(kz^{\prime})}-e^{-k(z-l)}\sinh{\left[k(z^{\prime}-l)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)-\sinh{[k(d-l)]}\right]$
$\displaystyle\qquad+\left[e^{-kz}\sinh{(kz^{\prime})}-e^{k(z-l)}\sinh{\left[k(z^{\prime}-l)\right]}\right]\left[\frac{\epsilon_{2}}{\epsilon_{1}}\left(1-\cosh{[k(d-l)]}\right)+\sinh{[k(d-l)]}\right]\Bigg{\\}},$
(32) $\displaystyle\Xi$ $\displaystyle=-\frac{4\pi
q}{\epsilon_{1}k\sinh{(kl)}}\left\\{\left(\frac{\epsilon_{2}}{\epsilon_{1}}+1\right)^{2}\sinh{(kd)}-\left(\frac{\epsilon_{2}}{\epsilon_{1}}-1\right)^{2}\sinh{[k(d-2l)]}-2\left[\left(\frac{\epsilon_{2}}{\epsilon_{1}}\right)^{2}-1\right]\sinh{(kl)}\right\\}^{-1},$
(33) $\displaystyle G(\vec{k},z;\vec{r}\mkern 2.0mu\vphantom{r}^{\prime})$
$\displaystyle=-\frac{2\pi
q}{\epsilon_{2}k\sinh{(kl)}}\left\\{\cosh{[k(z+z^{\prime}-l)]}-\cosh{[k\left(|z-z^{\prime}|-l\right)]}\right\\}.$
(34)

Applying the hyperbolic relations

$\displaystyle\cosh{[k(d-l)]}$
$\displaystyle=\cosh{(kd)}\cosh{(kl)}-\sinh{(kd)}\sinh{(kl)},$ (35)
$\displaystyle\sinh{[k(d-l)]}$
$\displaystyle=\sinh{(kd)}\cosh{(kl)}-\cosh{(kd)}\sinh{(kl)},$ (36)
$\displaystyle\sinh{[k(d-2l)]}$
$\displaystyle=\sinh{(kd)}\cosh{(2kl)}-\cosh{(kd)}\sinh{(2kl)},$ (37)

and rewriting in terms of

$\displaystyle\delta$
$\displaystyle=\frac{\epsilon_{2}/\epsilon_{1}-1}{\epsilon_{2}/\epsilon_{1}+1},$
(38)

we have upon taking $d\to\infty$:

$\displaystyle\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$ $\displaystyle\approx\frac{2\pi
q}{\epsilon_{2}k\sinh{(kl)}}\cdot\frac{1+\delta}{e^{2kl}-\delta^{2}}\cdot
e^{kl}\Big{\\{}e^{kl}\cosh{[k(z+z^{\prime}-l)]}-\cosh{[k(z-z^{\prime})]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\delta[\cosh{k(z-z^{\prime})}]-e^{-kl}\cosh{[k(z+z^{\prime}-l)]}\Big{\\}},$
(39) $\displaystyle\implies\varphi_{2}(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$
$\displaystyle=\tilde{\varphi}_{2}(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})+G(\vec{k},z;\vec{r}\mkern
2.0mu\vphantom{r}^{\prime})$ $\displaystyle\approx\frac{2\pi
q}{\epsilon_{2}k}\left\\{e^{-k|z-z^{\prime}|}+\frac{2\delta}{e^{2kl}-\delta^{2}}\left[\delta\cosh{[k(z-z^{\prime})]}+e^{kl}\cosh[k(z+z^{\prime}-l)]\right]\right\\}.$
(40)

##  Appendix C Correlated methods

Any antisymmetric $N$-particle state $\mathinner{|{\Psi}\rangle}$ can be
represented in the basis consisting of a reference Slater determinant
$\mathinner{|{\Phi}\rangle}$ and all of its excitations, i.e.

$\displaystyle\mathinner{|{\Psi}\rangle}$
$\displaystyle=\left(c_{0}+\sum_{ia}c_{i}^{a}\hat{a}_{a}^{\dagger}\hat{a}_{i}+\sum_{\begin{subarray}{c}i<j\\\
a<b\end{subarray}}c_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_{i}\hat{a}_{j}+\cdots\right)\mathinner{|{\Phi}\rangle},$
(41)

where $c_{i}^{a},c_{ij}^{ab},\dots$ are coefficients and
$\hat{a}_{i}^{\dagger}(\hat{a}_{i})$ are creation (annihilation) operators
defined over the orbital basis set $\\{\phi_{i}\\}$ from which
$\mathinner{|{\Phi}\rangle}$ is constructed [40]. This representation becomes
exact as the orbital basis approaches completeness. In principle, the
coefficients can be determined through exact diagonalization, but the problem
quickly becomes intractable as the Hilbert space scales exponentially with
system size. Correlated wave function methods thus seek to approximate
$\eqref{ci_wfn}$ starting from some known reference
$\mathinner{|{\Phi}\rangle}$—often obtained from Hartree-Fock theory (hence by
_correlation_ , we refer to contributions beyond Hartree-Fock). The heat-bath
configuration interaction (HCI) [30, 31] and coupled-cluster (CC) theories
[32, 33, 34] are two such methods that we employed in our study.

###  C.1 Heat-bath configuration interaction

Given a Hamiltonian $\hat{H}$, HCI finds an approximation to the ground state
$\mathinner{|{\Psi}\rangle}$ consisting of a selected set of the most
important determinants only. This is achieved in two stages: first, a
variational wave function is constructed from a set of chosen determinants;
second, perturbative corrections to the variational energy determine
additional determinants to be added. The selection of determinants in each
stage is tuned by two parameters, $\epsilon_{1}$ and $\epsilon_{2}$, which
control the tradeoff between speed and accuracy. We use vHCI and vHCI + PT2 to
denote the variational stage and its correction with perturbation theory (i.e.
the full HCI algorithm), respectively.

vHCI consists of iteratively expanding the selected space (which usually only
includes the Hartree-Fock determinant initially) while minimizing the
selected-space energy. At each iteration, new determinants
$\mathinner{|{\Phi^{\prime}_{i}}\rangle}$ connected to the current space by
nonzero Hamiltonian matrix elements
$H_{ij}=\mathinner{\langle{\Phi^{\prime}_{i}}|}\hat{H}\mathinner{|{\Phi_{j}}\rangle}$
are added according to the importance measure

$f(\mathinner{|{\Phi^{\prime}_{i}}\rangle})=\max_{j}(|H_{ij}c_{j}|),$ (42)

where $\\{c_{j}\\}$ are determinant coefficients. The variational stage is
summarized as follows [30]:

  1. 1.

Diagonalize the Hamiltonian in the selected space. Let $\\{c_{j}\\}$ denote
the determinant coefficients of the lowest eigenvector.

  2. 2.

Generate all exterior determinants
$\\{\mathinner{|{\Phi^{\prime}_{i}}\rangle}\\}$ that satisfy
$|H_{ij}c_{j}|>\epsilon_{1}$ for at least one determinant
$\mathinner{|{\Phi_{j}}\rangle}$ in the current space. Add these determinants
to the selected space.

  3. 3.

Repeat Steps 1-2 until the number of new determinants generated is $<1\%$ of
the number of determinants in the selected space.

The generation of $\\{\mathinner{|{\Phi^{\prime}_{i}}\rangle}\\}$ in Step 2 is
achieved via _deterministic heat-bath sampling_. Starting from a reference
determinant $\mathinner{|{\Phi_{j}}\rangle}$, only determinants
$\mathinner{|{\Phi^{\prime}_{i}}\rangle}$ connected to
$\mathinner{|{\Phi_{j}}\rangle}$ through matrix elements
$\mathinner{\langle{\Phi^{\prime}_{i}}|}\hat{H}\mathinner{|{\Phi_{j}}\rangle}$
that satisfy
$|\mathinner{\langle{\Phi^{\prime}_{i}}|}\hat{H}\mathinner{|{\Phi_{j}}\rangle}|>\epsilon$
are sampled. The threshold $\epsilon$ here is chosen to be
$\epsilon=\epsilon_{1}/|c_{j}|$. Duplicates are removed as the sampling
algorithm is applied to each determinant in the current space.

Once the variational wave function
$\mathinner{|{\Psi_{0}}\rangle}=\sum_{i}c_{i}\mathinner{|{\Phi_{i}}\rangle}$
is obtained, the partitioning of the Hamiltonian used in the perturbative
stage is defined as

$\displaystyle\hat{H}=\hat{H}_{0}+\hat{V},$
$\displaystyle\hat{H}_{0}=\sum_{ij}H_{ij}\mathinner{|{\Phi_{i}}\rangle}\mathinner{\langle{\Phi_{j}}|},\quad\hat{H}_{0}\mathinner{|{\Psi_{0}}\rangle}=E_{0}\mathinner{|{\Psi_{0}}\rangle}.$
(43)

The second-order correction to the energy is thus

$\displaystyle\Delta E_{2}$
$\displaystyle=\sum_{a}\frac{\left(\sum_{i}H_{ai}c_{i}\right)^{2}}{E_{0}-E_{a}},$
(44)

where the sum is over new determinants
$\\{\mathinner{|{\Phi^{\prime}_{a}}\rangle}\\}$ that give the matrix elements
$H_{ai}=\mathinner{\langle{\Phi^{\prime}_{a}}|}\hat{H}\mathinner{|{\Phi_{i}}\rangle}$
and $E_{a}=H_{aa}$. Since (44) is expensive to evaluate, vHCI + PT2
approximates $\Delta E_{2}$ using only terms in the sum that satisfy
$|H_{ai}c_{i}|>\epsilon_{2}$, i.e.

$\displaystyle\Delta
E_{2}\approx\sum_{a}\frac{\left(\sum_{i,|H_{ai}c_{i}|>\epsilon_{2}}H_{ai}c_{i}\right)^{2}}{E_{0}-E_{a}}.$
(45)

Similar to vHCI, the set $\\{\mathinner{|{\Phi^{\prime}_{a}}\rangle}\\}$ is
generated by applying the heat-bath sampling algorithm to each determinant in
$\mathinner{|{\Psi_{0}}\rangle}$ using a threshold of
$\epsilon=\epsilon_{2}/|c_{i}|$.

###  C.2 Coupled-cluster theory

Coupled-cluster (CC) theory [41] assumes an ansatz for the $N$-particle ground
state of the form

$\mathinner{|{\Psi_{\text{CC}}}\rangle}=e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle},$
(46)

where $\mathinner{|{\Psi_{0}}\rangle}$ is a reference Slater determinant
(often taken to be the Hartree-Fock solution). The _cluster operator_
$\hat{T}$ is given by

$\displaystyle\hat{T}$
$\displaystyle=\hat{T}_{1}+\hat{T}_{2}+\cdots+\hat{T}_{N},$ (47)
$\displaystyle\hat{T}_{1}$
$\displaystyle=\sum_{ia}t_{i}^{a}\hat{a}_{a}^{\dagger}\hat{a}_{i},$ (48)
$\displaystyle\hat{T}_{2}$ $\displaystyle=\sum_{\begin{subarray}{c}i<j\\\
a<b\end{subarray}}t_{ij}^{ab}\hat{a}_{a}^{\dagger}\hat{a}_{b}^{\dagger}\hat{a}_{i}\hat{a}_{j},$
(49) $\displaystyle\qquad\vdots$

where $\hat{T}_{i}$ are cluster operators associated with $i$-particle
excitations and $t_{i}^{a},t_{ij}^{ab},\dots$ are _cluster amplitudes_ to be
determined. These amplitudes together with the CC energy can be computed
through a projection of the Schrödinger equation onto
$\mathinner{|{\Psi_{0}}\rangle}$ and its excitations:

$\displaystyle\mathinner{\langle{\Psi_{0}}|}e^{-\hat{T}}\left(\hat{H}-E\right)e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle}$
$\displaystyle=0,$ (50)
$\displaystyle\mathinner{\langle{\Psi_{0}}|}\hat{a}_{i}^{\dagger}\hat{a}_{a}e^{-\hat{T}}\left(\hat{H}-E\right)e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle}$
$\displaystyle=0,$ (51)
$\displaystyle\mathinner{\langle{\Psi_{0}}|}\hat{a}_{j}^{\dagger}\hat{a}_{i}^{\dagger}\hat{a}_{b}\hat{a}_{a}e^{-\hat{T}}\left(\hat{H}-E\right)e^{\hat{T}}\mathinner{|{\Psi_{0}}\rangle}$
$\displaystyle=0,$ (52) $\displaystyle\ \ \vdots$ (53)

The resulting set of coupled polynomial equations (i.e. the _coupled-cluster
equations_) can be solved (often iteratively) for the energy and amplitudes.

Note that $\mathinner{|{\Psi_{\text{CC}}}\rangle}$ with $\hat{T}$ untruncated
is simply a non-linear reparametrization of (41), i.e. it is equivalent to
exact diagonalization. The power of the CC approach, however, lies in the
truncated form of the ansatz. The exponential form allows higher order
excitations of $\mathinner{|{\Psi_{0}}\rangle}$ to be included in the ansatz
despite truncating $\hat{T}$. One of the most widely used variants is the
truncation $\hat{T}=\hat{T}_{1}+\hat{T}_{2}$, known as coupled-cluster
singles-doubles (CCSD). The accuracy of CCSD can be naïvely improved by
including $\hat{T}_{3}$, but this approach is often too costly. A more
feasible way to improve CCSD solutions is to include contributions from
$\hat{T}_{3}$ perturbatively, leading to the coupled-cluster singles-doubles
plus perturbative triples (CCSD(T)) method which is often called the “gold
standard” of quantum chemistry.

##  Appendix D Computational details

###  D.1 GHF

The spinor structure of GHF orbitals breaks $\hat{S}^{2}$ and $\hat{S}_{z}$
symmetry, hence $\mathinner{|{\Psi_{\text{GHF}}}\rangle}$ is not constrained
to be an eigenstate of $\hat{S}^{2}$ and $\hat{S}_{z}$. Sacrificing symmetries
of the exact solution provides the additional flexibility necessary to obtain
the true variational minimum, which can be advantageous in studying
magnetically frustrated systems that require a non-spin-collinear description
[42]. As previous work [27] on TMD heterobilayers has suggested GWC ground
states with non-collinear spins at some fractional filling factors, GHF is an
appropriate mean-field method for our problem.

We performed all GHF calculations using a development version of Q-Chem [43].
Solutions were optimized through a hybrid scheme of first using the Direct
Inversion in the Iterative Subspace (DIIS) method [44, 45] for
$\mathcal{O}(1)$ iterations before switching to the Geometric Direct
Minimization (GDM) algorithm [46]. This approach combines DIIS’s ability to
efficiently converge to the _global minimum_ during the early iterations
together with GDM’s ability to traverse challenging energy landscape features
and robustly arrive at a _local minimum_. We executed supercell calculations
composed of $\sqrt{N_{\text{cell}}}\times\sqrt{N_{\text{cell}}}$ unit cells at
the $\Gamma$-point of the mBZ, which is analogous to performing unit cell
calculations at $N_{k}=N_{\text{cell}}$ $\vec{k}$-points adequately sampled
from the mBZ. Periodic boundary conditions are only imposed across the
supercell, which allow us to naturally obtain solutions where discrete
translational symmetry on the moiré superlattice is broken. To break $C_{3}$
symmetry, however, requires explicit encoding of the broken symmetry in the
initial guess.

For our survey of GHF solutions across filling factors in Figs. 1, 6, 8 and 7,
we started with a “noisy” spin-unpolarized initial guess constructed from
continuum model orbitals, i.e. eigenstates of (1) (henceforth referred to as
the “continuum model guess”). By “noisy”, we mean that matrix elements of the
spin $\uparrow\downarrow$ and $\downarrow\uparrow$ blocks of the density
matrix $\mathbf{P}$ are randomly sampled from the standard Gaussian
distribution $\mathcal{N}(0,1)$, normalized by
$10^{-3}\times\max|\mathbf{P}|$. Table 1 lists the filling factors and system
sizes we considered in this survey. At $\nu=1/2$, we additionally considered a
spin-polarized stripe guess constructed from orthogonal Bloch orbitals of the
form

$\psi_{\vec{k}}(\vec{r})\propto\sum_{i=1}^{N}e^{i\vec{k}\cdot\vec{R}_{i}}\phi_{i}(\vec{r}),$
(54)

where $N$ is the particle number, $\phi_{i}(\vec{r})$ are Gaussians centered
at stripe lattice sites $\vec{R}_{i}$, and $\vec{k}$ are crystal momenta that:

  1. (i)

are restricted to the first Brillouin zone defined by the reciprocal space of
the stripe lattice;

  2. (ii)

satisfy the periodic boundary conditions imposed on the supercell.

To evaluate the basis set convergence of solutions, we performed a series of
calculations with increasing basis set size $N_{\text{basis}}$ starting from a
small $N_{\text{basis}}^{\text{min}}$. The converged solution at each
$N_{\text{basis}}$ was input as the initial guess for the next calculation
using a larger $N_{\text{basis}}$, where the $N_{\text{basis}}^{\text{min}}$
solutions were first ensured to be stable (see Appendix K for details on
stability analysis).

Filling factor $\nu$ | Sytem size $N$
---|---
1/7 | 7, 28
1/4 | 4, 16, 36, 49, 64, 81
1/3 | 3, 12, 27, 48, 75, 108
2/5 | 10, 40
1/2 | 8, 18, 32, 50, 72, 98, 162
3/5 | 15, 60
2/3 | 6, 24, 54, 96, 150
3/4 | 12, 27, 48, 75, 108, 147
6/7 | 42, 168
1 | 9, 16, 36, 64, 81, 144
Table 1: The filling factors and system sizes studied with GHF.

###  D.2 Correlated methods

We performed correlated calculations using only spin-collinear variants of
HCI, CCSD, and CCSD (T). CCSD and CCSD(T) calculations were carried out using
PySCF [47, 48, 49], whereas HCI calculations with $\sim\mathcal{O}(10^{6})$
determinants and parameters $\epsilon_{1},\epsilon_{2}\approx 10^{-7}$ were
executed using a locally modified version of Dice [30, 31, 50].

Recall that our HF calculations utilized a complex plane wave basis
$\left\\{e^{i\vec{G}_{p}\cdot\vec{r}}\right\\}$ which yields 2-electron
integral tensors $\innerproduct{pq}{rs}$ possessing a 4-fold symmetry, i.e.

$\displaystyle\innerproduct{pq}{rs}$
$\displaystyle=\innerproduct{rs}{pq}=\innerproduct{sr}{qp}=\innerproduct{qp}{sr}.$
(55)

Quantum chemistry codes including PySCF, however, often assume the basis sets
to be real and hence $\innerproduct{pq}{rs}$ to be 8-fold symmetric, i.e.

$\displaystyle\innerproduct{pq}{rs}$
$\displaystyle=\innerproduct{rq}{ps}=\innerproduct{ps}{rq}=\innerproduct{rs}{pq}$
$\displaystyle=\innerproduct{qp}{sr}$
$\displaystyle=\innerproduct{sp}{qr}=\innerproduct{qr}{sp}=\innerproduct{sr}{qp}.$
(56)

Thus, to interface with such codes, we applied a unitary transformation to
rotate the plane wave basis into a real-valued basis composed of sines and
cosines (up to a normalization factor):

$\displaystyle\begin{bmatrix}e^{i\vec{G}_{p}\cdot\vec{r}}\\\
e^{-i\vec{G}_{p}\cdot\vec{r}}\end{bmatrix}$
$\displaystyle\xlongrightarrow{\quad\mathbf{U}\quad}\sqrt{2}\begin{bmatrix}\cos{(\vec{G}_{p}\cdot\vec{r})}\\\
\sin{(\vec{G}_{p}\cdot\vec{r})}\end{bmatrix},$ (57) $\displaystyle\mathbf{U}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{bmatrix}1&-i\\\ 1&i\end{bmatrix}.$
(58)

The resulting 2-electron integral tensor has 8-fold symmetry recovered.

##  Appendix E GHF charge densities across filling factors

Figure 6: GHF charge density $\rho(\vec{r})$ at filling factors $\nu\leq 1$
obtained from the unpolarized continuum model guess. The densities are plotted
over the supercell and the grid lines show the moiré unit cell boundaries. The
indices $n_{1},n_{2}\in\mathbb{Z}^{+}$ represent a point
$\vec{r}=n_{1}\vec{L}_{m,1}+n_{2}\vec{L}_{m,2}$ in the supercell, where
$\vec{L}_{m,i}$ are the primitive moiré lattice vectors. Figure 7: Fast
Fourier transform (FFT) of the charge densities in Fig. 6 with the average
density subtracted, i.e.
$|\rho(\vec{k})|=|\mathrm{FFT}\left\\{\rho(\vec{r})-\langle\rho\rangle\right\\}|$.
The indices $\kappa_{1},\kappa_{2}\in\mathbb{Z}$ represent a point
$\vec{g}=\kappa_{1}\vec{G}_{1}+\kappa_{2}\vec{G}_{2}$ in reciprocal space,
where $\vec{G}_{i}$ are the primitive reciprocal lattice vectors associated
with the supercell. Figure 8: Peaks extracted from Fig. 8. We use the term
“primary” to describe peaks with the largest magnitudes, while “secondary”
describes any feature of secondary importance. $\vec{G}_{m}$ denotes the
primitive moiré reciprocal vector.

In this section, we provide further discussion on the charge densities obtain
from the unpolarized continuum model guess. Figs. 6 and 7 show the largest
basis set charge densities across filling factors and their Fast Fourier
transforms (FFTs). Fig. 8 depicts Fourier-space peaks extracted from Fig. 7
(see Appendix G for FFT details). As indicated by the $C_{3}$ symmetry seen in
reciprocal space (Fig. 7), all solutions we found possess configurations that
obey the $C_{3}$ symmetry of the moiré superlattice. We did, however, find
states with expanded unit cells that break the superlattice translational
symmetry. This situation is possible because:

  1. (i)

the continuum model guess preserves $C_{3}$ symmetry;

  2. (ii)

the supercell constructed from
$\sqrt{N_{\text{cell}}}\times\sqrt{N_{\text{cell}}}$ unit cells preserves
$C_{3}$ symmetry; but

  3. (iii)

the periodic boundary conditions imposed only over the supercell allow the
superlattice translational symmetry to be broken.

The charge densities at $\nu=1/7,\ 1/3$, and $2/3$ exhibit the clearest GWC
charge-orders with expanded unit cells. At $\nu=1/3$ and $2/3$, FFT peaks at
$\vec{k}$ where $|\vec{k}|/|\vec{G}_{m}|=1/\sqrt{3}$ are observed, reflecting
a unit cell defined by the vectors

$\displaystyle\vec{L}_{1}$
$\displaystyle=2\vec{L}_{m,1}-\vec{L}_{m,2},\quad\vec{L}_{2}=\vec{L}_{m,1}+\vec{L}_{m,2}.$
(59)

Such a cell is constructed from occupied (vacant) moiré sites for $\nu=1/3$
($2/3$). At $\nu=1/7$, FFT peaks occur at $\vec{k}$ where
$|\vec{k}|/|\vec{G}_{m}|=1/\sqrt{7}$, indicating a unit cell defined by the
vectors

$\displaystyle\vec{L}_{1}$
$\displaystyle=3\vec{L}_{m,1}-\vec{L}_{m,2},\quad\vec{L}_{2}=\vec{L}_{m,1}+2\vec{L}_{m,2}.$
(60)

In contrast to the case with $\nu=1/3$ and $2/3$, the $\nu=1/7$ and $6/7$
charge densities we obtained are not dual 333In the sense that swapping the
occupied and vacant sites in the $\nu=1/7$ configuration gives the $\nu=6/7$
configuration to each other in Fig. 6 and can also be distinguished in Fourier
space (Figs. 7-8). This is likely because our solution at $\nu=6/7$ is one of
several local minima with similar configurations that are close in energy.

At $\nu=1/4,\ 2/5,\ 1/2,\ 3/5$, and $3/4$, $\rho(\vec{r})$ exhibits charge
leakage and delocalization between moiré superlattice sites due to quantum
fluctuations. While the primary peaks originating from the moiré superlattice
is the most pronounced feature in the FFT, extraction of peaks associated with
secondary features was less straightforward for these values of $\nu$—the
plotted values are our best estimates. At $\nu=1/4$, FFT peaks at $\vec{k}$
where $|\vec{k}|/|\vec{G}_{m}|=1/2$ suggest the presence of an expanded unit
cell defined by the vectors

$\displaystyle\vec{L}_{1}$
$\displaystyle=2\vec{L}_{m,1},\quad\vec{L}_{2}=2\vec{L}_{m,2}.$ (61)

At $\nu=1/2$, the similarity in the FFT peak locations with those of $\nu=2/3$
implies an underlying honeycomb structure. Finally, although there is greater
uncertainty in the secondary peak locations for $\nu=3/4,\ 2/5,$ and $3/5$, a
reasonable conclusion is that they also harbor similar periodic features as
the $\nu=2/3$ configuration.

##  Appendix F Ordered stripe phase at $\nu=1/2$

Figure 9: The ordered stripe phase at $\nu=1/2$ for various system sizes $N$.
Figure 10: Convergence of the energy difference between the ordered and
disordered (labyrinthine) stripe phases with system and basis set size. The
quantity $2|\vec{G}_{\mathrm{max}}/\vec{G}_{m}|$ is a measure of the basis set
size, where $|\vec{G}_{\mathrm{max}}|$ is the maximum momentum included in our
plane wave basis set. (For reference: $N_{\mathrm{basis}}=211$ for $N=18$
corresponds to $2|\vec{G}_{\mathrm{max}}/\vec{G}_{m}|\approx 2.5$.) As $N$
increases, the ordered stripe configuration becomes lower in energy but the
difference
$|\varepsilon_{\mathrm{ordered}}-\varepsilon_{\mathrm{disordered}}|$
decreases.

Fig. 9 shows the ordered stripe phase at various system sizes $N$ obtained by
explicitly constructing an initial guess with the desired order. We were able
to attain solutions breaking $C_{3}$ symmetry by breaking the symmetry in our
initial guess. Fig. 10 illustrates the convergence of the energy difference
per particle between the ordered and disordered stripe phases with system and
basis set size. As $N$ increases, the ordered stripe becomes lower in energy
but the two phases become more nearly degenerate, which is suggestive of
glassy behaviour. We verified that all solutions were stable in the charge
sector, although there were instabilities in the spin sector associated with
transitions from a polarized to an unpolarized state at
$2|\vec{G}_{\mathrm{max}}/\vec{G}_{m}|\approx 2.5$. Nevertheless, the energy
differences per hole between the stable and unstable solutions were $<1$ meV.

##  Appendix G Fast Fourier transform (FFT) of the charge density

Prior to computing the FFT of the charge densities $\rho(\vec{r})$ in Fig. 6,
we subtracted the average density $\langle\rho\rangle$ to remove the
$\vec{k}=0$ Fourier component. Fig. 7 shows the magnitude of the Fourier
components
$|\rho(\vec{k})|=|\mathrm{FFT}\left\\{\rho(\vec{r})-\langle\rho\rangle\right\\}|$.
At $\nu=1/7,\ 1/3,\ 2/3,\ 6/7,$ and $1$, the FFT exhibits clean peaks that can
be easily extracted. At $\nu=1/4,\ 2/5,\ 1/2,\ 3/5,$ and $3/4$, however, in
addition to the clear peaks located at the moiré reciprocal lattice sites, the
FFT displays secondary disperse features in reciprocal space associated with
charge delocalization in real space. We estimated the locations of possible
secondary peaks suppressed by the disperse features but nonetheless can still
be gauged by eye. An example is shown for $\nu=1/4$ in Fig. 11(a), where we
first removed the primary peaks at the moiré reciprocal lattice sites before
applying a Gaussian blur filter to smooth out $|\rho(\vec{k})|$ near the
secondary peaks. The locations (marked by the red dots) where
$|\rho(\vec{k})|$ is above some threshold value are then extracted. In cases
where there are clusters of similar $|\rho(\vec{k})|$ above the threshold
value, as for $\nu=3/4$ in Fig. 11(b), we identified the cluster center by
$k$-means clustering.

(a) $\nu=1/4$

(b) $\nu=3/4$

Figure 11: Illustration of our FFT secondary peak extraction by: (i) removing
the primary peaks at the moiré reciprocal lattice sites; (ii) applying a
Gaussian blur filter to smooth out $|\rho(\vec{k})|$ near the secondary peaks;
and (iii) extracting the locations where $|\rho(\vec{k})|$ is above some
threshold value directly or via $k$-means clustering.

##  Appendix H Estimating the charge gaps

The charge gap $\Delta_{c}$ is defined as

$\displaystyle\Delta_{c}$ $\displaystyle=\lim_{N\to\infty}\Delta_{c}(N),$ (62)
$\displaystyle\Delta_{c}(N)$ $\displaystyle=E(N+1)+E(N-1)-2E(N)$ (63)
$\displaystyle=-\left[E(N)-E(N+1)\right]+\left[E(N-1)-E(N)\right]$
$\displaystyle=-\text{EA}+\text{IP},$ (64)

where $\text{EA}=E(N)-E(N+1)$ is the electron affinity and
$\text{IP}=E(N-1)-E(N)$ is the ionization potential. Assuming that the
orbitals in the $N$ and $N\pm 1$ ground states are identical within Hartree-
Fock theory, Koopman’s theorem [40] allows us to identify

EA
$\displaystyle=-\varepsilon_{\text{LUMO}},\quad\text{IP}=-\varepsilon_{\text{HOMO}},$
(65)

where $\varepsilon_{\text{LUMO}}$ ($\varepsilon_{\text{HOMO}}$) is the lowest
unoccupied (highest occupied) molecular orbital. This gives

$\displaystyle\Delta_{c}(N)\approx\varepsilon_{\text{LUMO}}-\varepsilon_{\text{HOMO}}.$
(66)

In the thermodynamic limit, the approximation in (66) should become an exact
equality since the addition or removal of a single particle will have minimal
effect on the ground state orbitals. We estimated $\Delta_{c}$ using (66) in
our Madelung-corrected charge gap (70), which—compared to (63)—has the
advantage of:

  1. (i)

being cheaper to compute since it only requires one self-consistent loop.

  2. (ii)

ensuring $\Delta_{c}\geq 0$.

At $\nu=1/4,\ 1/3,\ 1/2,\ 2/3,\ 3/4$, and $1$ where there are sufficient data
points, we assumed the extrapolation formula for the Madelung-corrected charge
gaps (see Appendix I) to be:

$\displaystyle\tilde{\Delta}_{c}(N)=\frac{\gamma}{N^{3/2}}+\Delta_{c}(\infty),$
(67)

where $\gamma$ is a parameter to be fitted. At $\nu=1/7,\ 2/5,\ 3/5,$ and
$6/7$ where there are only two data points, we estimated $\Delta_{c}(\infty)$
using $\tilde{\Delta}_{c}(N)$ computed with the largest system size. A large
enough $N_{\text{basis}}^{\text{max}}$ was used such that all
$\tilde{\Delta}_{c}(N)$ are converged to $<0.5$ meV in Fig. 1.

##  Appendix I Finite-size errors

Finite-size errors inevitably arise in the study of bulk systems through
practical calculations performed on a finite number $N$ of particles. While
periodic boundary conditions help to reduce finite-size effects, non-
negligible errors still remain—especially for singular long-ranged potentials
such as the Coulomb interaction [52]. We therefore employed corrections
involving the triangular lattice Madelung energy, given by [53]

$\varepsilon_{M}=\frac{-3.921034}{2\epsilon\sqrt{\pi}}\cdot
r_{s}^{-1}N^{-1/2},$ (68)

where $\epsilon$ and $r_{s}$ are the effective dielectric constant and the
Wigner-Seitz radius, respectively. Such correction schemes were introduced to
remove spurious “self-interactions” between particles and their periodic
images that are artifacts of the periodic boundary conditions. We also
extrapolated our results to the thermodynamic limit where appropriate.

For a $N$-particle system, the Madelung-corrected total energy per particle
$\varepsilon(N)=E(N)/N$ and charge gap $\Delta_{c}(N)=E(N+1)+E(N-1)-2E(N)$ are
computed as [54, 52]

$\displaystyle\tilde{\varepsilon}(N)$
$\displaystyle=\varepsilon(N)+\varepsilon_{M},$ (69)
$\displaystyle\tilde{\Delta}_{c}(N)$
$\displaystyle=\Delta_{c}(N)+2|\varepsilon_{M}|,$ (70)

which improves the scaling in the finite-size errors from
$\mathcal{O}(N^{-1/2})$ to $\mathcal{O}(N^{-1})$ for $\varepsilon$ [52] and to
$\mathcal{O}(N^{-3/2})$ for $\Delta_{c}$ [55]. Our extrapolations of
$\tilde{\varepsilon}$ and $\tilde{\Delta}_{c}$ to $N\to\infty$ thus assume
$\mathcal{O}(N^{-1})$ and $\mathcal{O}(N^{-3/2})$ scalings, respectively. As
illustrative examples, Fig. 12 shows the convergence of $\tilde{\varepsilon}$
and $\tilde{\Delta}_{c}$ with $N$ at $\nu=1/2,3/4,$ and $1$. We found that the
convergence is $\nu$-dependent and observed that the variation in
$\tilde{\varepsilon}$ as a function of $N$ is generally smaller than in
$\tilde{\Delta}_{c}$. The values do not converge monotonically because they
may be associated with disordered configurations that differ slightly from one
another.

Figure 12: Convergence of the Madelung-corrected total energy per particle
$\tilde{\varepsilon}$ and charge gap $\tilde{\Delta}_{c}$ with system size
$N$. The values obtained using the largest basis set with the unpolarized
continuum model guess are shown.

##  Appendix J Spin-collinearity test

We evaluated the spin collinearity of our GHF solutions following the approach
of [56]. A wave function $\mathinner{|{\Psi}\rangle}$ is considered to be
collinear if it is an eigenstate of an axial spin operator
$\hat{S}_{c}=\vec{c}\cdot\hat{\vec{S}}$ associated with some Cartesian
coordinate, where $\vec{c}\in\mathbb{R}^{3}$ is a unit vector and
$\hat{\vec{S}}=(\hat{S}_{x},\hat{S}_{y},\hat{S}_{z})$ is the usual many-body
spin vector operator. In other words, $\mathinner{|{\Psi}\rangle}$ is
collinear if it satisfies the relation

$\displaystyle\hat{S}_{c}\mathinner{|{\Psi}\rangle}=s_{c}\mathinner{|{\Psi}\rangle},$
(71) $\displaystyle s_{c}\in\\{-N/2,-N/2+1,\dots,N/2-1,N/2\\}.$

The collinearity test for a state consists of computing the quantities $\mu$
and $\varepsilon_{0}$, defined as

$\displaystyle\mu$
$\displaystyle=\langle(\hat{S}_{c})^{2}\rangle-\langle\hat{S}_{c}\rangle^{2},$
(72) $\displaystyle\varepsilon_{0}$
$\displaystyle=\norm{\langle\hat{\vec{S}}\rangle}=\norm{\left(\langle\hat{S}_{x}\rangle,\langle\hat{S}_{y}\rangle,\langle\hat{S}_{z}\rangle\right)},$
(73)

where $\norm{\cdot}$ denotes the Euclidean norm. The conditions for
collinearity are thus:

  1. (i)

$\mu=0$.

  2. (ii)

$\varepsilon_{0}\in\\{-N/2,-N/2+1,\dots,N/2-1,N/2\\}.$

Condition (i) is obvious since it can easily be proved that $\mu=0$ if and
only if $\mathinner{|{\Psi}\rangle}$ is an eigenstate of $\hat{S}_{c}$.
Condition (ii) is less straightforward and we refer the reader to [56] for its
justification. In practice, we deem $\mu\lesssim 10^{-4}$ as describing a
collinear state. Fig. 13 shows $\mu$ and $\epsilon_{0}$ computed for the GHF
solutions initialized with the continuum model guess in our initial survey
across $\nu$.

Figure 13: The quantities $\mu$ and $\varepsilon_{0}$ of GHF solutions
initialized from the unpolarized continuum model guess. While the smallest
$N$-particle system we studied at each $\nu$ may not be collinear, the largest
$N$-particle system is always collinear. This suggests that the emergence of
non-collinear solutions could be an artifact of finite-size effects. In the
order of increasing $\nu$, the smallest (largest) $N$ we considered were 7
(28), 4 (81), 3 (108), 10 (40), 8 (162), 15 (60), 6 (150), 12 (147), 42 (168),
and 9 (144).

##  Appendix K Stability analysis

The self-consistent procedure used in GHF only guarantees solutions to be
stationary points. To determine whether they are local minima or saddle
points, we inspect the lowest eigenvalues of the electronic Hessian following
the approach of [57]. Diagonalization is executed via a finite-difference
implementation of the Davidson method, where for a vector $\vec{v}$ in the
subspace, the Hessian-vector product is computed as

$\mathbf{H}\vec{v}\approx\left[\frac{\nabla E(\mathbf{C}_{+})-\nabla
E(\mathbf{C}_{-})}{2\xi}\right]^{\mathsf{T}}.$ (74)

$\nabla E$ is the analytic energy gradient, while $\mathbf{C}_{\pm}$ are the
GHF solution coefficient matrices perturbed in opposite directions determined
by the finite step size $\xi$. If the solution is deemed unstable, we can
displace the solution along the lowest eigenvector and use this as an initial
guess in a new GHF calculation.