# Dynamical defects in a two-dimensional Wigner crystal:
self-doping and kinetic magnetism

Kyung-Su Kim  Corresponding author.
<EMAIL_ADDRESS>Department of Physics, Stanford University, Stanford, CA
93405  Ilya Esterlis  Department of Physics, University of Wisconsin-Madison,
Madison, Wisconsin 53706, USA  Chaitanya Murthy  Department of Physics,
Stanford University, Stanford, CA 93405  Steven A. Kivelson  Department of
Physics, Stanford University, Stanford, CA 93405

###### Abstract

We study the quantum dynamics of interstitials and vacancies in a two-
dimensional Wigner crystal (WC) using the semi-classical instanton method that
is asymptotically exact at low density, i.e., in the $r_{s}\to\infty$ limit.
The interstitial dynamics favors local ferromagnetism, while vacancy dynamics
mediates antiferromagnetism. The resulting semi-classical expression for the
interstitial energy vanishes at $r_{s}=r_{\rm mit}\approx 100$, signaling a
possible self-doping instability to a partially melted WC for some range of
$r_{s}$ below $r_{\rm mit}$. We thus propose that there exists a ferromagnetic
“metallic electron crystal” phase of the two-dimensional electron gas at
intermediate densities between a low density insulating WC and a high density
Fermi fluid.

##  I Introduction

Despite its prime importance in the field of condensed matter physics, some
basic aspects remain unsettled concerning the physics of the two-dimensional
electron gas (2DEG) at intermediate densities where various forms of “strongly
correlated electron fluids” can arise. The ideal 2DEG is governed by the
simple Hamiltonian

$\displaystyle
H=\sum_{i}\frac{\textbf{p}_{i}^{\,2}}{2m}+\sum_{i<j}\frac{e^{2}}{4\pi\epsilon}\frac{1}{|\textbf{r}_{i}-\textbf{r}_{j}|},$
(1)

with a single dimensionless parameter, $r_{s}=a_{0}/a_{\rm B}$, characterizing
the ratio of the typical interaction strength to the kinetic energy. Here,
$a_{0}=1/\sqrt{\pi n}$ is the average interparticle distance, $n$ is the
electron density, and $a_{\mathrm{B}}=4\pi\epsilon\hbar^{2}/me^{2}$ is the
effective Bohr radius. The phases of the 2DEG in the weak and strong coupling
limits are well-understood: it forms a paramagnetic Fermi liquid (FL) for
small $r_{s}$ (weak coupling) and a Wigner crystal (WC) for large $r_{s}$
(strong coupling) [1]. The present study addresses the intermediate coupling
regime near the quantum metal-insulator transition (MIT). Landmark numerical
studies suggested that the MIT occurs as a direct transition from a Fermi
liquid to an insulating WC at $r_{s}=r_{\rm melt}^{*}\approx 31$ [2, 3, 4].
However, recent experiments [5, 6, 7, 8] suggest that the actual transition
may be more complex.

Figure 1:  Conjectured $T=0$ phases of a clean 2DEG as a function of
$1/r_{s}\propto\sqrt{n}$: WC (Ferro) = fully polarized ferromagnetic WC; WC
(Antiferro) = WC with some form of antiferromagnetism (or a spin liquid
phase); Metallic electron crystal (MeC, Ferro) = metallic electron crystal
characterized by more than one electron per unit cell and partial (or full)
spin polarization; FL (Para) = paramagnetic Fermi liquid. The phase transition
at $r_{\rm F}^{\rm wc}\approx 175$ [9] is due to the change of dominant
exchange interactions from ferromagnetic to antiferromagnetic and is likely to
be first order. $r_{\rm mit}\approx 100$ indicates the “true” metal-insulator
transition due to interstitial self-doping proposed in this paper, and is
distinct from $r_{\rm melt}$ below which the crystalline order vanishes.
$r_{\rm melt}$ is expected to be smaller than the value for a direct FL–WC
transition from quantum Monte Carlo calculations, $r_{\rm melt}^{*}\approx 31$
[4], due to the existence of the intermediate MeC phase. See Sec. V for a
detailed discussion of the conjectured phase diagram.

Apart from the charge ordering, there is another subtle issue regarding the
magnetism. In the FL regime, the paramagnetic state seems to be most favored
[4]. Deep within the WC phase, the magnetism is determined by various ring-
exchange processes. The exchange coefficients can be calculated using the
semi-classical instanton approximation [10, 11, 12, 13, 14], the validity of
which is well-tested by a numerically exact path integral Monte Carlo
calculation [9]. These calculations show that the WC is a ferromagnet for
large enough $r_{s}>r_{\rm F}^{\rm wc}\approx 175$ [9] and some form of a
(highly frustrated) antiferromagnet [12] below $r_{\rm F}^{\rm wc}$ (Fig. 1).
However, the predicted energy scale for the ring-exchange processes is too
small to account for the insulating ferromagnetic phase observed in the large
$r_{s}$ regime of AlAs quantum wells [5, 7], prompting some of the present
authors to propose a (higher-temperature) kinetic mechanism for such a phase
mediated by interstitial hopping processes [15].

In this paper, using the semi-classical instanton approximation, we carry out
a comprehensive study of the quantum dynamics of an interstitial and a vacancy
defect (Fig. 2), two point defects of a WC with the smallest classical
energies [16, 17] 111Note that the energetics of these point defects are also
studied using the path integral Monte Carlo method in Ref. [53], but the
contributions from defect hopping and exchange processes are ignored.. We
first review the formulation of the standard instanton technique and apply it
to derive effective Hamiltonians describing various exchange and defect
hopping processes illustrated in Fig. 3 (Sec. II). While the interstitial
dynamics induces ferromagnetic correlations in the WC [15], vacancy dynamics
favors some form of antiferromagnetism. This predicted particle-hole asymmetry
of magnetism could be experimentally tested by controlled doping of a WC that
is commensurately locked to a weak periodic substrate potential, as discussed
in Sec. III. In Sec. IV, we estimate the energy of an interstitial and a
vacancy from the semi-classical calculations. Interestingly, the semi-
classical expression for the interstitial energy, when extrapolated to a large
but finite $r_{s}$, vanishes around $r_{s}=r_{\rm mit}\approx 100$, signaling
a possible self-doping instability to partially melt the WC. From this, we
propose the existence of a ferromagnetic metallic electron crystal (MeC) phase
as an intermediate phase of the 2DEG (Sec. V). Figure 1 summarizes the main
findings of our study. We conclude with a remark on the fate of the phase
diagram in the presence of weak quenched disorder in Sec. VI.

##  II The Semi-classical approximation

We first review the standard semi-classical instanton method as applied to the
ideal 2DEG (1) in the large $r_{s}$ limit. The exact partition function of the
(fermionic) 2DEG is

$\displaystyle Z=\int d^{2N}\textbf{r}_{0}\sum_{P\in
S_{N}}\frac{(-1)^{P}}{N!}\sum_{\boldsymbol{\sigma}}\left<P\textbf{r}_{0},P\boldsymbol{\sigma}\right|e^{-\beta
H}\left|\textbf{r}_{0},\boldsymbol{\sigma}\right>,$ (2)
$\displaystyle\left<P\textbf{r}_{0},P\boldsymbol{\sigma}\right|e^{-\beta
H}\left|\textbf{r}_{0},\boldsymbol{\sigma}\right>=\delta_{\boldsymbol{\sigma},P\boldsymbol{\sigma}}\left<P\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>,$ (3)
$\displaystyle\left<\textbf{r}_{0}^{\prime}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>=\int_{\tilde{\textbf{r}}(0)=\tilde{\textbf{r}}_{0}}^{\tilde{\textbf{r}}(\tilde{\beta})=\tilde{\textbf{r}}_{0}^{\prime}}D\tilde{\textbf{r}}(\tau)e^{-\sqrt{r_{s}}S},$
(4) $\displaystyle
S=\int_{0}^{\tilde{\beta}}d\tau\left[\frac{1}{2}\left(\frac{d{\tilde{\textbf{r}}}}{d\tau}\right)^{2}+V(\tilde{\textbf{r}})-V_{0}\right],$
(5) $\displaystyle
V(\textbf{r})\equiv\sum_{i<j}\frac{1}{|\textbf{r}_{i}-\textbf{r}_{j}|},$ (6)

where $\textbf{r}(\tau)\equiv\\{\textbf{r}_{i}(\tau)\\}$ are the positions of
$N$ electrons in imaginary time,
$\textbf{r}_{0}\equiv\\{\textbf{r}_{i}(\tau=0)\\}$ are their initial
positions, $\boldsymbol{\sigma}\equiv\\{\sigma_{i}=\,\uparrow,\downarrow\\}$
are their respective spin indices, and $\beta=1/k_{B}T$ is the inverse
temperature. The sum over $N!$ permutations, $P,$ of the coordinates and the
sign factor $(-1)^{P}$ encode the fermionic exchange statistics. For bosonic
particles, one should merely substitute $(-1)^{P}\rightarrow+1$. The 2DEG
Hamiltonian (1) does not act on the electron spins, hence the
$\delta_{\boldsymbol{\sigma},P\boldsymbol{\sigma}}$ factor in the second line
above. The third and fourth lines are the path integral representation of the
$N$-electron propagator. The action is rescaled to make the $r_{s}$ dependence
manifest by introducing dimensionless coordinates,
$\tilde{\textbf{r}}\equiv{\textbf{r}}/a_{0}$, and dimensionless imaginary
time, $\tau$. Correspondingly, $\tilde{\beta}\equiv\beta E^{*}$ is a
dimensionless inverse temperature, where $E^{*}\equiv e^{2}/(4\pi\epsilon
a_{\rm B}r_{s}^{3/2})$. The path integral measure is also defined as an
integration over the dimensionless coordinate $\tilde{\textbf{r}}(\tau)$. The
minimum potential energy $V_{0}={\rm
min}_{\tilde{\textbf{r}}}V(\tilde{\textbf{r}})$ is subtracted for later
convenience 222When considering exchange processes in a pure WC, $V_{0}$
corresponds to the classical WC energy. For tunneling processes involving a
defect, $V_{0}$ corresponds to the classical energy of the defect.. The
Coulomb interaction (last line) is computed numerically using the standard
Ewald method. As usual, the presence of a uniform neutralizing positively-
charged background is assumed. Henceforth, we will drop tildes from the
rescaled coordinates to simplify the notation:
$\tilde{\textbf{r}}\to\textbf{r}.$ We focus on the zero temperature phase of
the problem, and hence will always take $\beta\to\infty$ in the end.

We approach this problem using a semi-classical instanton approximation, which
is asymptotically exact in the $r_{s}\to\infty$ (strong coupling) limit. In
Sec. II.1, we briefly review the semi-classical derivation of ring-exchange
processes in the WC. In Sec. II.2 and II.3, we consider tunneling processes
involving a single interstitial and vacancy, respectively, and derive the
corresponding effective Hamiltonians describing their dynamics. The
application of the semi-classics to a bosonic system is addressed in Sec.
II.4.

Figure 2: (a) A classical centered interstitial and (b) a vacancy
configuration. Small black arrows are drawn to indicate the positions of the
interstitial (a) and the vacancy [left panel of (b)]. The vacancy
configuration has $D_{2}$ symmetry, and not the full $D_{6}$ symmetry of the
underlying WC; therefore, a vacancy has three possible orientations $\alpha$.
We introduce a pictorial notation for the vacancy for later convenience.
Figure 3: Tunneling processes considered in this paper. (a) WC exchange
processes. (b) Exchange processes involving an interstitial. (c) Interstitial
hopping processes. (d) Exchange processes involving a vacancy. (e) Vacancy
hopping processes. In (b,c), black arrows indicate the positions of
interstitials. In (e), a black (cyan) oval denotes an initial (final) vacancy
configuration corresponding to each vacancy hopping process.
$t_{11},t_{12},t_{22},t_{23}$ exhaust all the nearest-neighbor vacancy hopping
processes; others are related to one of these by symmetry. Panels (a-c) are
adapted from Ref. [15].

###  II.1 Wigner crystal ring-exchange processes

In the $r_{s}\to\infty$ limit, the classical ground state manifold consists of
a triangular lattice WC with $2^{N}$-fold degeneracy in spin states. The
lifting of this degeneracy and the nature of the resulting magnetic order is
determined for $1\ll r_{s}<\infty$ by WC ring-exchange processes. Various
ring-exchange processes correspond to distinct instanton solutions of the
action and can be calculated via the dilute instanton approximation [20, 11,
10, 12, 13, 14, 15], which we briefly review below. (See Refs. [13, 14] for
more details.) The result is an effective spin Hamiltonian expressed as a sum
over all ring-exchange processes:

$\displaystyle
H_{\mathrm{eff}}^{\mathrm{wc}}=-\sum_{a}(-1)^{P_{a}}\,{J_{a}}\,(\mathcal{\hat{P}}_{a}+\mathcal{\hat{P}}^{-1}_{a}),$
(7)

where the semi-classical calculation gives a leading-order large $r_{s}$
asymptotic expression for $J_{a}$. Here, $\mathcal{\hat{P}}_{a}$ is the
permutation operator corresponding to the permutation $P_{a}$, and can be
decomposed as a product of two-particle exchange operators. The two-particle
exchange operators, in turn, can be written in terms of spin operators as
$\mathcal{\hat{P}}_{(i,j)}=2(\vec{S}_{i}\cdot\vec{S}_{j}+\frac{1}{4})$.

To illustrate how this works, recall the familiar problem of the semi-
classical calculation of the tunnel splitting in a symmetric double-well
potential [21, 22, 23]. For large enough $\beta$ such that
$\beta\hbar\omega_{0}\gg 1$, the excited states in each well can be neglected.
(Here, $\omega_{0}$ is the oscillation frequency in either well.) In this
limit, we obtain asymptotic relations

$\displaystyle\left<\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>\sim|\psi(\textbf{r}_{0})|^{2}e^{-\beta
E_{0}}\cosh(\beta\Delta),$ (8)
$\displaystyle\left<-\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>\sim|\psi(\textbf{r}_{0})|^{2}e^{-\beta
E_{0}}\sinh(\beta\Delta),$ (9)

where the minima of the two wells are at $\pm\textbf{r}_{0}$,
$|\psi(\textbf{r}_{0})|^{2}=|\psi(-\textbf{r}_{0})|^{2}$ is the probability
density of the wave function at these positions, and $E_{0}$ and $2\Delta$
are, respectively, the mean energy and the splitting between the even and odd
parity ground states. The right-hand side of each expression is obtained by
inserting the resolution of the identity on the left-hand side.

In viewing this same problem from the path integral perspective in the semi-
classical limit, one first solves for the instanton path—the smallest action
path that begins at the bottom of one well and ends at the bottom of the
other. The net duration (in imaginary time) of this tunneling event is of
order $\omega_{0}^{-1}$. We then sum over multiple such instanton events to
obtain an expression of the same form as above, where the diagonal (off-
diagonal) propagator in Eq. 8 (Eq. 9) contains all the terms with an even
(odd) number of events. Expanding these expressions in power series, one sees
that the typical number of tunneling events is $\sim\beta\Delta$ and the mean
imaginary time interval between them is of order $\hbar/\Delta$. Note that in
the semi-classical limit $\hbar/\Delta\gg\omega_{0}^{-1}$, the instantons are
dilute and hence effectively non-interacting (see Fig. 4). Looked at another
way, for a range of temperature such that $\hbar\omega_{0}\gg T\gg\Delta$,
where multiple instanton events can be neglected, we can compute $\Delta$ as

$\displaystyle\Delta=\beta^{-1}\,\frac{\left<-\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{\text{1-inst}}}{\left<\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{\text{0-inst}}},$ (10)

where the subscripts designate the number of instanton events.

Figure 4: An example of a multi-instanton configuration for the double well
potential shown in the inset. The “size” of an instanton in imaginary time is
$\sim 1/\hbar\omega_{0}$ and the “distance” between them is $\sim 1/\Delta$.

The analysis is somewhat more complicated but structurally similar for the
present problem. Consider the propagator
$\left<P_{a}\textbf{r}_{0}\right|e^{-\beta H}\left|\textbf{r}_{0}\right>$
where $\textbf{r}_{0}$ is an initial WC configuration and $P_{a}$ is the
permutation corresponding to a particular ring-exchange process [see Fig.
3(a)]. In the semi-classical (large $r_{s}$) limit, this propagator is again
expressible as a weighted sum over multi-instanton contributions. For
temperatures such that $\hbar\Omega\gg T\gg J_{a}$, where $J_{a}$ is the
tunnel splitting corresponding to the process $P_{a}$, the propagator is
dominated (up to symmetry) by a single “$a$” instanton contribution associated
with the path $\textbf{r}^{(a)}(\tau)$ with the smallest action subject to the
boundary conditions $\textbf{r}^{(a)}(0)=\textbf{r}_{0}$ and
$\textbf{r}^{(a)}(\tilde{\beta})=P_{a}\textbf{r}_{0}$. Here,
$\hbar\Omega/2\sim r_{s}^{-3/2}$ is the zero-point energy of the WC, while
$J_{a}$ is exponentially small in $\sqrt{r_{s}}$ at large $r_{s}$. The
single-$a$-instanton contribution to the propagator can be expressed as

$\displaystyle\left<P_{a}\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{a,\text{1-inst}}$ $\displaystyle\approx
e^{-\sqrt{r_{s}}S_{a}}\int_{\delta\textbf{r}(0)=\textbf{0}}^{\delta\textbf{r}(\tilde{\beta})=\textbf{0}}D\delta\textbf{r}(\tau)\exp\\!\left[-\frac{\sqrt{r_{s}}}{2}\int_{0}^{\tilde{\beta}}\delta\textbf{r}(\tau)^{T}{\bf{\hat{M}}}^{(a)}(\tau)\delta\textbf{r}(\tau)\right]$
$\displaystyle=e^{-\sqrt{r_{s}}S_{a}}\left(\det\big{[}\sqrt{r_{s}}\
{\bf{\hat{M}}}^{(a)}(\tau)\big{]}\right)^{-1/2},$ (11)
$\displaystyle\hat{M}^{(a)}_{ij}(\tau)$
$\displaystyle\equiv\frac{\delta^{2}S}{\delta r^{(a)}_{i}(\tau)\,\delta
r^{(a)}_{j}(\tau)}=-\delta_{ij}\frac{\partial^{2}}{\partial\tau^{2}}+\partial_{i}\partial_{j}V\big{[}\textbf{r}^{(a)}(\tau)\big{]},$
(12)

where $S_{a}\equiv S\big{[}\textbf{r}^{(a)}(\tau)\big{]}$ with the trajectory
$\textbf{r}^{(a)}(\tau)$ satisfying $\delta
S\big{[}\textbf{r}^{(a)}(\tau)\big{]}=0$, and
$\delta\textbf{r}(\tau)\equiv\textbf{r}(\tau)-\textbf{r}^{(a)}(\tau)$ is the
fluctuation coordinate. Fluctuations are treated within a harmonic
approximation around the semi-classical path. In Eq. 12, the derivative
$\partial_{i}$ is with respect to the normalized coordinates. Note that
${\bf\hat{M}}^{(a)}$ has a zero eigenvalue solution
$\dot{\textbf{r}}^{(a)}(\tau)$ corresponding to the translation in imaginary
time, which has to be treated with care [22, 21, 23, 14, 13]. Separating the
zero mode contribution from the full determinant, one obtains

$\displaystyle\left<P_{a}\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{a,\text{1-inst}}$ (13)
$\displaystyle=\beta\,\frac{e^{2}}{4\pi\epsilon a_{\rm
B}r_{s}^{3/2}}\sqrt{\frac{S_{a}}{2\pi}}\cdot e^{-\sqrt{r_{s}}S_{a}}\left({\rm
det}^{\prime}\big{[}\sqrt{r_{s}}\,{\bf\hat{M}}^{(a)}(\tau)\big{]}\right)^{-\frac{1}{2}},$

where the prime denotes that the zero eigenvalue must be omitted in the
calculation of the determinant. Note that since an instanton is a localized
object with a characteristic size $\Delta\tau=\Omega^{-1}$, one can neglect
the exponentially small correction from its tail provided $\beta\hbar\Omega\gg
1$.

On the other hand, the diagonal propagator in the zero instanton sector
$\left<\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{\text{0-inst}}$ can be obtained by making a
harmonic approximation of $V$ around $\textbf{r}_{0}$:

$\displaystyle\left<\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{\text{0-inst}}\approx\left(\det\big{[}\sqrt{r_{s}}\,{\bf\hat{M}}^{(0)}(\tau)\big{]}\right)^{-\frac{1}{2}},$
(14) $\displaystyle\
{\hat{M}}^{(0)}(\tau)\equiv-\delta_{ij}\frac{\partial^{2}}{\partial\tau^{2}}+\partial_{i}\partial_{j}V(\textbf{r}_{0}).$
(15)

Normalizing the propagator in the one instanton sector by that in the zero
instanton sector, as in Eq. 10, one obtains

$\displaystyle J_{a}$
$\displaystyle=\beta^{-1}\,\frac{\left<P_{a}\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{a,\text{1-inst}}}{\left<\textbf{r}_{0}\right|e^{-\beta
H}\left|\textbf{r}_{0}\right>|_{\text{0-inst}}}$
$\displaystyle=\frac{e^{2}}{4\pi\epsilon a_{\rm
B}}\cdot\frac{{A_{a}}}{{r_{s}^{5/4}}}\sqrt{\frac{S_{a}}{2\pi}}\
e^{-\sqrt{r_{s}}S_{a}}>0,$ (16) $\displaystyle A_{a}$
$\displaystyle=\left[\frac{{\rm
det}^{\prime}\big{(}-\partial_{\tau}^{2}+V^{\prime\prime}\big{[}\textbf{r}^{(a)}(\tau)\big{]}\big{)}}{\det\big{(}-\partial_{\tau}^{2}+V^{\prime\prime}(\textbf{r}_{0})\big{)}}\right]^{-\frac{1}{2}},$
(17)

where $A_{a}$ is called a “fluctuation determinant,” calculated in the
normalized coordinates with $r_{s}=1$, and the $\beta\to\infty$ limit is
implicitly taken in the end. In the second line, the extra factor of
$r_{s}^{1/4}$ comes from the normalization of the determinant

$\displaystyle\left(\frac{{\rm
det}^{\prime}\big{[}\sqrt{r_{s}}\,{\bf\hat{M}}^{(a)}(\tau)\big{]}}{\det\big{[}\sqrt{r_{s}}\,{\bf\hat{M}}^{(0)}(\tau)\big{]}}\right)^{-\frac{1}{2}}=r_{s}^{1/4}\left(\frac{{\rm
det}^{\prime}\big{[}{\bf\hat{M}}^{(a)}(\tau)\big{]}}{\det\big{[}{\bf\hat{M}}^{(0)}(\tau)\big{]}}\right)^{-\frac{1}{2}}.$

Hence, $A_{a}$ (17) (and also $S_{a}$) are dimensionless numbers with no
$r_{s}$ dependence. In Eq. 17, $V^{\prime\prime}$ denotes the hessian matrix
of $V$. We refer readers to Appendix A for the details of the numerical
calculation of $S_{a}$ and $A_{a}$. For the ring-exchange processes
illustrated in Fig. 3(a), we quote the results for $S_{a}$ and $A_{a}$ from
Ref. [14]: $S_{2}=1.64$, $A_{2}=1.30$; $S_{3}=1.53$, $A_{3}=1.10$;
$S_{4}=1.66$, $A_{4}=1.24$; $S_{5}=1.91$, $A_{5}=1.57$; $S_{6}=1.78$,
$A_{6}=1.45$. Our calculations, and those of Ref. [13], agree with these
values.

The remaining issue concerns the sign factor $(-1)^{P_{a}}$ that enters
$H_{\rm eff}^{\rm wc}$ in Eq. 7, which is due to the anti-symmetry of the
many-body electronic wave function (see chapter V of Ref. [11] for an
explanation). As recognized by Thouless [20], this implies that a ring-
exchange process involving an even (odd) number of electrons mediates an
antiferromagnetic (ferromagnetic) interaction.

Table 1: Dimensionless actions $S_{a}$ and fluctuation determinants $A_{a}$
for tunneling processes illustrated in Fig. 3(b-e) calculated in this paper.
The parentheses in the first and fourth columns denote $(N_{\rm move},M)$,
where $N_{\rm move}$ is the number of electrons that are allowed to adjust
their positions during minimization and $M$ is the number of time slices for
the discretized tunneling paths (i.e., there are $M-1$ intermediate
configurations). Processes for a centered interstitial (vacancy) are
calculated in a hexagonal supercell with $12\times 12+1$ ($10\times 10-1$)
electrons starting and ending at fully relaxed defect configurations.

###  II.2 Processes involving a single interstitial

Tunneling processes involving a single centered interstitial (CI) defect [Fig.
2(a)] were first considered in Ref. [15]. We refine the results obtained there
in two ways: (1) by improving the estimate of the classical action (which is
done by solving the classical equations of motion for a finite sized system
with periodic boundary conditions) using a hexagonal, instead of a
rectangular, supercell with $12\times 12+1$ electrons; and (2) by explicitly
calculating the fluctuation determinants $A_{a}$ rather than simply making
dimensional estimates. Fig. 3(b-c) show the tunneling processes considered in
this paper with the corresponding $S_{a}$ and $A_{a}$ listed in Table 1. The
hopping matrix elements $t_{a}>0$ are again expressed in terms of $S_{a}$ and
$A_{a}$ as in Eq. 16. Note that four hopping processes have smaller actions
than those of exchange processes and hence are more important when $r_{s}\gg
1$. The effective Hamiltonian in the presence of a dilute concentration of
interstitials is [15]

$\displaystyle H_{\mathrm{eff}}^{\rm i,hop}$
$\displaystyle=-t_{2}\sum_{\begin{subarray}{c}(n,j,n^{\prime})\\\
\in(t_{2}\>\mathrm{path})\end{subarray}}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}c^{\dagger}_{n,\sigma^{\prime}}f^{\dagger}_{j,\sigma}f_{j,\sigma^{\prime}}c_{n^{\prime},\sigma}$
$\displaystyle\ -t_{2}^{\prime}\sum_{\begin{subarray}{c}(n,j,n^{\prime})\\\
\in(t_{2}^{\prime}\>\mathrm{path})\end{subarray}}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}c^{\dagger}_{n,\sigma^{\prime}}f^{\dagger}_{j,\sigma}f_{j,\sigma^{\prime}}c_{n^{\prime},\sigma}$
$\displaystyle\
-t_{2}^{\prime\prime}\sum_{\begin{subarray}{c}(n,j,n^{\prime})\\\
\in(t_{2}^{\prime\prime}\>\mathrm{path})\end{subarray}}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}c^{\dagger}_{n,\sigma^{\prime}}f^{\dagger}_{j,\sigma}f_{j,\sigma^{\prime}}c_{n^{\prime},\sigma}$
$\displaystyle\
-t_{1}\sum_{\left<n,n^{\prime}\right>}\sum_{\sigma=\uparrow,\downarrow}c^{\dagger}_{n,\sigma}c_{n^{\prime},\sigma}\,+\,\left[U=\infty\right]+\cdots,$
(18)

where $f_{j\sigma}^{\dagger}$ ($c_{n,\sigma}^{\dagger}$) is the creation
operator of electrons that live on the WC sites $j$ (triangular plaquette
centers $n$) and the $U=\infty$ condition precludes any double occupancy. The
omitted terms correspond to hopping processes other than those shown in Fig.
3(c) and direct (elastic) interactions between interstitials [24]. Modified
exchange interactions near an interstitial are understood to be absorbed in
$H_{\rm eff}^{\rm wc}$ (7).

###  II.3 Processes involving a single vacancy

The classical vacancy defect has $D_{2}$ symmetry instead of the full $D_{6}$
symmetry of the underlying triangular lattice [17] [see Fig. 2(b)]. Therefore,
associated with each location of a vacancy, there are 3 inequivalent
orientations related by $C_{6}$ rotations. We will denote these by an index
$\alpha=1,2,3$; $\alpha=2$ and $3$ are related to $\alpha=1$ by
$\mathcal{C}_{6}$ and $\mathcal{C}_{6}^{2}$ respectively.

We considered tunneling processes involving a single vacancy defect as
illustrated in Fig. 3(d-e), with their corresponding values of $S_{a}$ and
$A_{a}$ listed in Table 1. The calculation is done in a hexagonal supercell
containing $10\times 10-1$ electrons. Again, matrix elements
$\Delta,t_{a},J_{a,{\rm V}}>0$ are given by Eq. 16 in terms of $S_{a}$ and
$A_{a}$. Note that, as in the interstitial case, the tunnel barriers
(determined by $S_{a}$) for hopping processes are smaller than those for
exchange processes. In particular, the tunnel barrier for $t_{11}$ is
anomalously smaller than that for any other processes, implying that the
dynamics of the vacancy is predominantly uni-directional in the
$r_{s}\to\infty$ limit 333Such a peculiar defect dynamics has also been
predicted to occur in solid helium [54]. 444Such a “restricted mobility” would
in turn lead to an exponentially slow (in $\sqrt{r_{s}}$) thermalization in
the presence of a dilute vacancy concentration..

The resulting effective Hamiltonian describing the dynamics of vacancies can
be written straightforwardly as follows. First, corresponding to each
orientation $\alpha$ of a vacancy, we introduce a (hard-core) bosonic operator
$b_{i,\alpha}^{\dagger}$ that suitably relaxes the positions of the WC
electrons near the vacancy site $i$ to the associated configuration that
minimizes the (classical) Coulomb energy. Then, with the definitions
$\textbf{b}_{i}^{\dagger}\equiv[b_{i,1}^{\dagger},b_{i,2}^{\dagger},b_{i,3}^{\dagger}]$
and

$\displaystyle\mathfrak{D}\equiv\begin{bmatrix}0&\Delta&\Delta\\\
\Delta&0&\Delta\\\ \Delta&\Delta&0\end{bmatrix},\
\Upsilon\equiv\begin{bmatrix}t_{11}&t_{12}&t_{12}\\\ t_{12}&t_{22}&t_{23}\\\
t_{12}&t_{23}&t_{22}\end{bmatrix},\ \mathcal{C}_{6}=\begin{bmatrix}0&1&0\\\
0&0&1\\\ 1&0&0\end{bmatrix},$ (19)

the effective Hamiltonian in the presence of a dilute concentration of
vacancies is

$\displaystyle H_{\rm eff}^{\rm
v,hop}=-\sum_{i}\Bigg{[}\textbf{b}^{\dagger}_{i}\mathfrak{D}\textbf{b}_{i}+\sum_{\delta=\pm\textbf{e}_{1}}f^{\dagger}_{i,\sigma}f_{i+\delta,\sigma}\textbf{b}^{\dagger}_{i+\delta}\Upsilon\textbf{b}_{i}$
$\displaystyle+\sum_{\delta=\pm\textbf{e}_{2}}f^{\dagger}_{i,\sigma}f_{i+\delta,\sigma}\textbf{b}^{\dagger}_{i+\delta}\
\mathcal{C}_{6}^{-1}\Upsilon\mathcal{C}_{6}\ \textbf{b}_{i}$
$\displaystyle+\sum_{\delta=\pm\textbf{e}_{3}}f^{\dagger}_{i,\sigma}f_{i+\delta,\sigma}\textbf{b}^{\dagger}_{i+\delta}\
\mathcal{C}_{6}^{-2}\Upsilon\mathcal{C}_{6}^{2}\
\textbf{b}_{i}\Bigg{]}+\left[U=\infty\right]+\cdots,$ (20)

where $\textbf{e}_{1}=[1,0]$, $\textbf{e}_{2}=[1/2,\sqrt{3}/2]$,
$\textbf{e}_{3}=[-1/2,\sqrt{3}/2]$, and $f_{i,\sigma}^{\dagger}$ is again the
creation operator of an electron living at the WC site $i$. The first term
describes on-site orientation-mixing processes corresponding to $\Delta$; the
second term describes vacancy hopping processes in the $\pm\textbf{e}_{1}$
directions; and the third (fourth) term describes vacancy hopping processes in
the $\pm\textbf{e}_{2}$ ($\pm\textbf{e}_{3}$) directions, which can be related
to the second term by $\mathcal{C}_{6}$ ($\mathcal{C}_{6}^{2}$) rotation [see
Fig. 3(e)]. The omitted terms correspond to hopping processes other than those
shown in Fig. 3(e) and direct (elastic) interactions between vacancies [24].
The modified exchange interactions are again absorbed in $H_{\rm eff}^{\rm
wc}$ (7).

###  II.4 Two-Dimensional Bose Gas

For Coulomb-interacting bosonic particles, one merely needs to substitute
$(-1)^{P_{a}}\to+1$ in $H_{\rm eff}^{\rm wc}$ (7) without changing the forms
of $H_{\rm eff}^{\rm i,hop}$ and $H_{\rm eff}^{\rm v,hop}$ (II.2 and II.3).
The consequence is that all ring-exchange processes and interstitial and
vacancy hopping processes mediate ferromagnetism. This is a special case of a
more general result that the ground state of an interacting multi-component
bosonic system is a fully polarized ferromagnet [27, 28].

##  III Kinetic magnetism

Here, we discuss the magnetic correlations induced by defect hopping processes
555The potential importance of defect hopping processes in the magnetism of a
quantum crystal was pointed out earlier in Refs. [20, 55]..

As shown in Ref. [15], the interstitial dynamics due to $H_{\rm eff}^{\rm
i,hop}$ (II.2) favors a fully polarized ferromagnet. In the presence of small
antiferromagnetic WC exchange interactions, a single interstitial can only
delocalize in a finite region, forming a large ferromagnetic polaron of size
$\sim a_{0}^{2}\sqrt{t/J}$, where $t$ and $J$ are appropriate sums of
interstitial hopping matrix elements and WC antiferromagnetic exchange
coefficients, respectively. At $r_{s}\approx 45$, we estimate that a single
interstitial polarizes about $\sim 40$ WC spins.

On the other hand, it is known that the dynamics of a single hole in the
$U=\infty$ Hubbard model on a non-bipartite lattice leads to some form of
antiferromagnetism [30, 31, 32, 33]; therefore, a single vacancy mediates
antiferromagnetic correlations around it. In the presence of competing
exchange interactions of the underlying WC, a vacancy forms a finite-sized
antiferromagnetic polaron.

Such a striking particle-hole asymmetry of magnetism may be directly tested by
controlled doping of a WC in the presence of a smoothly varying weak external
periodic potential. Consider a weak commensurate potential that has minima
$-W<0$ at the triangular WC sites. When the density is tuned away (but not too
far away) from the commensurate value, the defect-doped commensurately locked
WC can be shown to have a smaller energy than the pure incommensurate WC as
follows: The defect-doped commensurate WC has an energy per electron $\Delta
E_{\rm comm}/N\approx-W+E_{\rm def}\cdot|\delta|+O(W|\delta|,\delta^{2})$ as
compared to the pure incommensurate WC, where $\delta$ is the ratio of defect
electrons to the total number of electrons and $E_{\rm def}=E_{\rm CI}$
($E_{\rm V}$) is the energy of an interstitial (vacancy) defect in the absence
of the external potential. Therefore, for a range of doping $-W/E_{\rm
V}<\delta<W/E_{\rm CI}$, the system will form a commensurately locked defect-
doped WC phase. Such a phase, in turn, is characterized by ferromagnetism
(antiferromagnetism) for interstitial doping $\delta>0$ (vacancy doping
$\delta<0$). Such a proposal may be experimentally tested in certain Moiré
systems that support a commensurately locked WC phase [34, 35, 36, 37, 38]
666Such a particle-hole asymmetry was also pointed out in a related, but
distinct, system in the presence of a Zeeman field [56]..

We close this section with a remark on the (fermionic) 2DEG with an
approximate $SU(N)$ symmetry. In such a case, the analysis in Sec. II can be
straightforwardly generalized by introducing more than two components of the
spin indices $\sigma_{i}$ in the effective Hamiltonians (7, II.2, II.3). As in
the $SU(2)$ case [15], the interstitial dynamics can be proven to lead to a
fully polarized $SU(N)$ ferromagnet, whereas the vacancy dynamics will induce
local antiferromagnetism as suggested from the generalized “counter-Nagaoka”
theorem [40]. Such considerations would apply to a multi-valley 2DEG with an
isotropic effective mass tensor in the limit of small band filling. For
example, a hexagonal system with band minima located at $K$ and $K^{\prime}$
points of the Brillouin zone retains an emergent $SU(4)$ system in the limit
of dilute filling of the bands when various band structure effects (e.g.,
trigonal warping) can be ignored 777The enlarged $SU(4)$ symmetry of such a
system is only an approximate symmetry in the limit of dilute filling of the
bands. This approximate $SU(4)$ symmetry is broken by small band structure
effects, which qualitatively change the low-energy spectrum. For example, the
full $SU(4)$ polarization induced by interstitial hopping processes will
result in $3$ magnon modes with $\omega\sim k^{2}$ low-energy spectrum. Taking
into account small band structure effects, the symmetry is reduced to $SU(2)$,
inducing a gap (proportional to the degree of symmetry breaking) to $2$ magnon
modes. This in turn leads to a difference in the thermodynamic behavior at
temperatures above and below the gap scale. Similar symmetry breaking effects
in the low-energy spectrum would also be present in the antiferromagnetic
phases. .

##  IV Energy of Wigner crystal defects

Figure 5: The energy of an interstitial ($E_{\rm CI}$) and a vacancy ($E_{\rm
V}$) as a function of $r_{s}$ in the fully spin polarized sector from the
semi-classical expressions (21, 25). $E_{\rm CI}(r_{s})$ [$E_{\rm V}(r_{s})$]
crosses zero around $r_{\rm mit}\approx 100$ [$r_{s}\approx 28$]. $E_{\rm V}$
is an overestimate because the vacancy kinetic energy would be further reduced
in an antiferromagnetic background by $O(t_{\rm V})$. As such, the value of
$r_{s}$ at which $E_{\rm V}(r_{s})$ crosses zero is expected to be larger than
the estimate $r_{s}\approx 28$ obtained from Eq. 25.

In this section, we estimate the energy of an interstitial and a vacancy from
the results of the semi-classical calculations.

If the small antiferromagnetic exchange interactions of the underlying WC are
neglected, the dilute concentration of interstitials in $H_{\rm eff}^{\rm
i,hop}$ leads to fully polarized ferromagnetism [15]. In such a case, $H_{\rm
eff}^{\rm i,hop}$ (II.2) reduces to an interstitial hopping problem on a
honeycomb lattice. The minimum possible kinetic energy for a single
interstitial is achieved when it is fully delocalized with quasi-momentum
$\textbf{k}=0$. Including the classical Coulomb energy and the zero-point
vibrational energy, the minimum interstitial energy (in units of
$e^{2}/4\pi\epsilon a_{\rm B}$) is

$\displaystyle E_{\rm CI}=\frac{C_{1,{\rm CI}}}{r_{s}}+\frac{C_{3/2,{\rm
CI}}}{r_{s}^{3/2}}-\left(3t_{1}+3t_{2}+6t_{2}^{\prime}+6t_{2}^{\prime\prime}\right)+\cdots,$
(21)

where the hopping matrix elements implicitly depend on $r_{s}$. Here, we have
neglected terms corresponding to higher order perturbative corrections (i.e,
higher powers of $r_{s}^{-1/2}$) from phonon anharmonicity as well as non-
perturbative corrections to the semi-classical instanton approximation. We
calculated $C_{1,{\rm CI}}$ and $C_{3/2,{\rm CI}}$ for supercells up to size
$28\times 28+1$; extrapolation to an infinite supercell size gives $C_{1,{\rm
CI}}=0.0769$ and $C_{3/2,{\rm CI}}=-0.295$ 888 Our results for $C_{1,{\rm
CI}}$ and $C_{3/2,{\rm CI}}$ agree with Ref. [17]. Note the factor $2$
difference between our values and those of Ref. [17] due to the difference in
the unit of energy..

Due to the presence of small antiferromagnetic exchange interactions in the
underlying WC, a single interstitial will only delocalize in a finite region,
forming a ferromagnetic polaron. This effectively increases the energy of an
interstitial by a small amount as measured relative to the energy of the WC
999Equation 21 is the energy of an interstitial in the fully polarized sector
as compared to the fully polarized pure WC. As the WC is expected to be
antiferromagnetic at $r_{s}<r_{\rm F}^{\rm wc}\approx 175$, the true
interstitial energy should be calculated by comparing the energy of a single
interstitial in some background antiferromagnetic WC to that of the
antiferromagnetic WC without an interstitial..

Estimating the energy of a vacancy is much more involved due to its kinetic
frustration and the associated antiferromagnetism. We proceed by evaluating
the vacancy energy in the fully spin polarized sector, which will give an
overestimate compared to the true minimum energy in some antiferromagnetic
background within the semi-classical approximation. We first construct a many-
body basis in the single vacancy sector with fully polarized WC spins as

$\displaystyle\left|i,\alpha\right>$ $\displaystyle\equiv
b^{\dagger}_{i,\alpha}f_{i,\uparrow}\,f^{\dagger}_{1,\uparrow}\cdots
f^{\dagger}_{N,\uparrow}\left|\emptyset\right>$
$\displaystyle=(-1)^{i}b^{\dagger}_{i,\alpha}f^{\dagger}_{1,\uparrow}\cdots
f^{\dagger}_{i-1,\uparrow}f^{\dagger}_{i+1,\uparrow}\cdots
f^{\dagger}_{N,\uparrow}\left|\emptyset\right>.$ (22)

Then, we obtain the following three-band tight binding model for a single
vacancy in the fully polarized sector:

$\displaystyle H_{\rm eff}^{\rm v,hop}|_{\rm
FM,1V}=\sum_{i}\Big{[}-\left|i\right>\mathfrak{D}\left<i\right|+\sum_{\delta=\pm\textbf{e}_{1}}\left|{i+\delta}\right>\Upsilon\left<i\right|$
$\displaystyle+\sum_{\delta=\pm\textbf{e}_{2}}\left|i+\delta\right>\,\mathcal{C}_{6}^{-1}\Upsilon\mathcal{C}_{6}\,\left<i\right|+\sum_{\delta=\pm\textbf{e}_{3}}\left|i+\delta\right>\,\mathcal{C}_{6}^{-2}\Upsilon\mathcal{C}_{6}^{2}\,\left<i\right|\Big{]},$
(23)

where we introduced a vectorial notation
$\left|i\right>\equiv[\left|i,1\right>,\left|i,2\right>,\left|i,3\right>]$ and
$\left<i\right|\equiv[\left<i,1\right|,\left<i,2\right|,\left<i,3\right|]^{\rm
T}$. Note the change of signs for the hopping matrix elements in the vacancy
basis (IV) as compared to the electron basis (II.3). This sign-changing is
accountable for the kinetic frustration of a vacancy in the ferromagnetic
background and the resulting antiferromagnetism as discussed in the previous
section [33, 40]. As usual, this Hamiltonian can be block-diagonalized in
quasi-momentum space:

$\displaystyle\left<\textbf{k}\right|H_{\rm eff}^{\rm v,hop}|_{\rm
FM,1V}\left|\textbf{k}\right>=-\mathfrak{D}+2\cos(\textbf{k}\cdot\textbf{e}_{1})\Upsilon$
$\displaystyle\quad\ +2\cos(\textbf{k}\cdot\textbf{e}_{2})\
\mathcal{C}_{6}^{-1}\Upsilon\mathcal{C}_{6}+2\cos(\textbf{k}\cdot\textbf{e}_{3})\
\mathcal{C}_{6}^{-2}\Upsilon\mathcal{C}_{6}^{2}.$ (24)

From this, we numerically find the minimum possible vacancy kinetic energy,
$-t_{\rm V}$, using the values in Table 1. The semi-classical expression for
the vacancy energy in the fully spin polarized sector (in units of
$e^{2}/4\pi\epsilon a_{\rm B}$) is then

$\displaystyle E_{\rm V}=\frac{C_{1,{\rm V}}}{r_{s}}+\frac{C_{3/2,{\rm
V}}}{r_{s}^{3/2}}-t_{\rm V}+\cdots,$ (25)

where, again, the neglected terms correspond to higher order perturbative and
non-perturbative corrections. The true minimum kinetic energy of a single
vacancy in some antiferromagnetic background should be smaller than this
estimate by $O(t_{\rm V}).$ We calculated $C_{1,{\rm V}}$ and $C_{3/2,{\rm
V}}$ for supercells up to $24\times 24-1$; extrapolation to an infinite
supercell size gives $C_{1,{\rm V}}=0.1094$ and $C_{3/2,{\rm V}}=-0.368$.
Whereas our result for $C_{1,{\rm V}}$ agrees with Ref. [17] up to the fourth
digit, our value for $C_{3/2,{\rm V}}$ is slightly different from theirs
101010Our calculation is performed in a hexagonal supercell as opposed to that
of Ref. [17] in a rectangular supercell..

$E_{\rm CI}$ and $E_{\rm V}$ are plotted as a function of $r_{s}$ in Fig. 5.
The semi-classical expression for the interstitial (vacancy) energy, when
extrapolated to a finite $r_{s}$, vanishes around $r_{\rm mit}\approx 100$
($r_{s}\approx 28$). This indicates a possible instability of the WC to
interstitial self-doping for $r_{s}<r_{\rm mit}$.

##  V Intermediate Phases of the 2DEG

Our predicted value of $r_{\rm mit}\approx 100$ is larger than $r_{\rm
melt}^{*}\approx 31$, where $r_{\rm melt}^{*}$ is the value below
which—according to existing variational calculations—the energy of the
paramagnetic Fermi liquid becomes smaller than the energy of a WC (with a
particular assumed antiferromagnetic order) [4]. This suggests that there is a
range of densities, $r_{\rm melt}<r_{s}<r_{\rm mit}$, for which a metallic
electron crystal (MeC) phase with more than one electron per crystalline unit
cell is stable 111111The possibility of a quantum crystal in which the
particle density per crystalline unit cell is different from an integer value
was originally conceived of in Ref. [57] in the context of solid helium..
Here, $r_{\rm melt}$ is the value below which the crystalline order vanishes
(see Fig. 1). To the best of our knowledge, the proposed MeC phase with more
than one electron per crystalline unit cell has not been studied using the
variational quantum Monte Carlo method. A related, but distinct, MeC phase
with less than one electron per crystalline unit cell has been studied in
Refs. [46] and [4]; however, these studies are in disagreement with each
other.

Since an interstitial forms a large ferromagnetic polaron in a WC when the
dominant WC exchange interaction is antiferromagnetic [15], the MeC phase
occurring for $r_{\rm melt}<r_{s}<r_{\rm mit}$ (Fig. 1) is likely accompanied
by partial or full spin polarization. The partially spin polarized MeC phase
near $r_{\rm mit}$ is expected to be characterized by an anomalously large
quasi-particle (interstitial) effective mass due to its accompanying
ferromagnetic cloud. Such massive ferromagnetic polarons have a tendency to
agglomerate, leading to phase separation [47, 48, 49] and rendering the
transition at $r_{\rm mit}$ to be first-order 121212This argument ignores a
subtlety regarding the formation of mesoscale density-modulated
“microemulsion” phases in a (possibly narrow) regime about any putative first-
order transition in the presence of Coulomb interactions [58, 59, 60, 61]..
Note that for single component (spinless) electrons, polaron formation is not
an issue so the self-doping transition may be continuous.

There could well be multiple phases within the MeC region of the phase
diagram. For instance, assuming that the density of interstitials for $r_{s}$
only slightly below $r_{\rm mit}$ is sufficiently small, the resulting state
will be less than fully spin polarized. As $r_{s}$ is further decreased, the
density of interstitials increases with an accompanying increase in spin
polarization, which could eventually lead to full spin polarization within the
MeC phase.

Finally, the transition to a fully melted Fermi fluid occurs when the energy
of the MeC phase crosses that of the liquid phase. Existing variational
quantum Monte Carlo estimates of the critical $r_{s}$ involve comparing the
energy of the liquid to that of the insulating WC. If, as we have suggested,
the MeC has lower energy than the insulating WC in an intermediate $r_{s}$
regime, it would presumably be stable against quantum melting at somewhat
higher densities (smaller $r_{s}$). Thus, this carries with it the likely
implication that $r_{\rm melt}<r_{\rm melt}^{*}\approx 31$. We also note that
the application of an in-plane magnetic field would further stabilize the MeC
phase to even smaller $r_{s}<r_{\rm melt}$, by effectively increasing the
energy of the Fermi fluid (without much affecting the energy of the already
polarized MeC).

##  VI Effects of weak disorder

Before concluding, we remark on the effect of small quenched disorder on the
phase diagram (Fig. 1). Firstly, even weak disorder is expected to destroy any
long-range crystalline order; hence all the electronic crystalline states we
have discussed are defined only in an approximate sense as short-range ordered
states. Also, the ferromagnetic MeC phase is characterized by the reduced
density of mobile electrons and their increased effective mass; hence, even
weak disorder is likely to result in strong localization and destroy the
metallic character of the phase. The resulting disorder-induced intermediate
insulating ferromagnetic phase may be an explanation for the recently observed
insulating ferromagnetism in AlAs quantum wells [5, 7]. Note that such a
proposal predicts the vanishing of ferromagnetism for $r_{s}>r_{\rm mit}$,
giving way to an insulating antiferromagnetic phase with much smaller exchange
energy scale 131313This explanation for insulating ferromagnetism is distinct
from, although related to, a previous proposal [15] based on WC-FL puddle
formation..

## Acknowledgement

We thank Boris Spivak for initial insights which led to this investigation and
Akshat Pandey for collaboration on a previous work. K-S.K acknowledges the
hospitality of the Massachusetts Institute of Technology, where this work was
completed. I.E. acknowledges Eugene Demler, Hongkun Park, Jiho Sung, Pavel
Volkov, Jue Wang, and Yubo Yang for helpful discussions on related work.
K-S.K. and SAK were supported in part by the Department of Energy, Office of
Basic Energy Sciences, Division of Materials Sciences and Engineering, under
contract DE-AC02-76SF00515 at Stanford. I.E. was supported by AFOSR Grant No.
FA9550-21-1-0216 and the University of Wisconsin–Madison. C.M. was supported
in part by the Gordon and Betty Moore Foundation’s EPiQS Initiative through
GBMF8686. Parts of the computing for this project were performed on the
Sherlock computing cluster at Stanford University.

##  Appendix A Numerical calculations of $S$ and $A$

In this section, we review a numerical method for calculating $S_{a}$ (5) and
$A_{a}$ (17), closely following Ref. [14]. Although we applied the semi-
classical instanton calculation to the 2DEG specifically, the method outlined
here applies to any system with a general potential $V(\textbf{r})$ with
degenerate minima in the semi-classical limit. We first calculate the
instanton action $S_{a}$ (5) by discretizing a tunneling path:

$\displaystyle S_{a}$
$\displaystyle=\int_{{\textbf{r}}_{0}}^{{\textbf{r}_{0}^{\prime}}}d\textbf{r}\,\sqrt{2\Delta
V(\textbf{r})}$ (26) $\displaystyle\approx\sum_{k=1}^{N_{\rm
time}}\frac{1}{2}|\textbf{r}_{k}-\textbf{r}_{k-1}|\cdot\big{[}\sqrt{2\Delta
V(\textbf{r}_{k})}+\sqrt{2\Delta V(\textbf{r}_{k-1})}\big{]},$

where we defined $\Delta V(\textbf{r})\equiv V(\textbf{r})-V(\textbf{r}_{0})$
and used the semi-classical equation of motion $\ddot{\textbf{r}}=\nabla
V(\textbf{r})$. $\textbf{r}_{0}$ and $\textbf{r}_{0}^{\prime}$ are initial and
final minimum configurations of $V$, respectively,
$\textbf{r}_{k}\equiv\textbf{r}(\tau_{k})$ is the collective coordinate of
particles at time $\tau_{k}$, where
$0\equiv\tau_{0}<\tau_{1}<\tau_{2}<\cdots<\tau_{M}\equiv\tilde{\beta}$, and
$\textbf{r}_{M}\equiv\textbf{r}_{0}^{\prime}$. In order to make the distances
$|\textbf{r}_{k}-\textbf{r}_{k-1}|$ approximately equal, each $\textbf{r}_{k}$
is taken to be constrained in the hyperplane defined by
$(\textbf{r}_{k}-\textbf{r}_{0})\cdot(\textbf{r}_{0}^{\prime}-\textbf{r}_{0})=\frac{k}{M}|\textbf{r}_{0}^{\prime}-\textbf{r}_{0}|$.
Numerical minimization of the discretized action (26) is performed with a
standard optimization package [52]. We will henceforth denote by
$\textbf{r}_{k}$ ($k=0,1,..,M$) the optimized tunneling path for the
$a$-instanton process.

The fluctuation determinant $A_{a}$ captures the Gaussian fluctuations around
the semi-classical path

$\displaystyle
A_{a}=\frac{F^{\prime}\big{[}\textbf{r}^{(a)}(\tau)\big{]}}{F\big{[}\textbf{r}_{0}\big{]}}=\left[\frac{{\rm
det}^{\prime}\big{(}-\partial_{\tau}^{2}+V^{\prime\prime}\big{[}\textbf{r}^{(a)}(\tau)\big{]}\big{)}}{\det\big{(}-\partial_{\tau}^{2}+V^{\prime\prime}(\textbf{r}_{0})\big{)}}\right]^{-\frac{1}{2}},$
(27) $\displaystyle
F\big{[}\textbf{r}(\tau)\big{]}\equiv\int_{\delta\textbf{r}(0)=0}^{\delta\textbf{r}(\tilde{\beta})=0}D\delta\textbf{r}(\tau)\exp\left[-\frac{1}{2}\int_{0}^{\tilde{\beta}}\left(\delta\dot{\textbf{r}}(\tau)^{2}+\delta\textbf{r}(\tau)^{\rm
T}V^{\prime\prime}\big{[}\textbf{r}^{(a)}(\tau)\big{]}\delta\textbf{r}(\tau)\right)\right]=\left<\textbf{0}\right|\mathcal{T}\exp(-\tilde{\beta}\,h[\textbf{r}^{(a)}(\tau)])\left|\textbf{0}\right>,$
(28) $\displaystyle
h[\textbf{r}(\tau)]\equiv-\frac{1}{2}\nabla^{2}+\frac{1}{2}\delta\textbf{r}(\tau)^{\rm
T}V^{\prime\prime}\big{[}\textbf{r}(\tau)\big{]}\delta\textbf{r}(\tau),$ (29)

where $\mathcal{T}\exp(\cdots)$ denotes the imaginary-time-ordered
exponential,
$\delta\textbf{r}(\tau)\equiv\textbf{r}(\tau)-\textbf{r}^{(a)}(\tau)$ is the
fluctuation coordinate, and the primed determinant in the first line is again
computed with the zero mode omitted. $\tilde{\beta}\to\infty$ is implicitly
taken in the end in calculating $A_{a}.$ As discussed below, the calculation
of $A_{a}$ can be done numerically by first computing
${F\big{[}\textbf{r}^{(a)}(\tau)\big{]}}/{F\big{[}\textbf{r}_{0}\big{]}}$ that
includes the zero mode contribution, and then multiplying by the square root
of the smallest eigenvalue (which is exponentially small in $\tilde{\beta}$)
of the operator
$-\partial_{\tau}^{2}+V^{\prime\prime}\big{[}\textbf{r}^{(a)}(\tau)\big{]}$.

$F\big{[}\textbf{r}^{(a)}(\tau)\big{]}$ can be calculated by discretizing the
path integral expression (28). First, we further define the time slices
intermediate to those defined above

$\displaystyle
0<\tau_{1/2}<\tau_{1}<\tau_{1+1/2}<\cdots<\tau_{M-1/2}<\tau_{M}\equiv\tilde{\beta},$
(30)

where each interval,
$\Delta\tau_{k}\equiv\tau_{k+\frac{1}{2}}-\tau_{k-\frac{1}{2}}$
($k=1,\cdots,M-1$), is calculated by inverting the semi-classical equation of
motion

$\displaystyle\Delta\tau_{k}$
$\displaystyle\equiv\int_{\textbf{r}_{k-\frac{1}{2}}}^{\textbf{r}_{k+\frac{1}{2}}}\frac{d\textbf{r}}{\sqrt{2\Delta
V[\textbf{r}^{(a)}(\tau)]}}$ $\displaystyle\approx\frac{1}{\sqrt{2\Delta
V(\textbf{r}_{k})}}\cdot\frac{1}{2}\left(|\textbf{r}_{k+1}-\textbf{r}_{k}|+|\textbf{r}_{k}-\textbf{r}_{k-1}|\right),$
(31)

and analogously for the end intervals,
$\Delta\tau_{0}\equiv\tau_{\frac{1}{2}}\approx\frac{1}{\sqrt{2\Delta
V(\textbf{r}_{0})}}\cdot\frac{1}{2}|\textbf{r}_{1}-\textbf{r}_{0}|$ and
$\Delta\tau_{M}\equiv\tau_{M}-\tau_{M-\frac{1}{2}}\approx\frac{1}{\sqrt{2\Delta
V(\textbf{r}_{M})}}\cdot\frac{1}{2}|\textbf{r}_{M}-\textbf{r}_{M-1}|$. (Note
that the end intervals formally diverge, $\Delta\tau_{0,M}\to\infty$, as
$\tilde{\beta}\to\infty$.) Then, the propagator at each interval can be
approximated by that of the quantum harmonic oscillator (Mehler kernel) of
$h[\textbf{r}^{(a)}(\tau)]\approx h_{k}\equiv h[\textbf{r}_{k}]$

$\displaystyle\left<\delta\textbf{r}_{k+\frac{1}{2}}\right|e^{-\Delta\tau_{k}h_{k}}\left|\delta\textbf{r}_{k-\frac{1}{2}}\right>=\prod_{n=1}^{2N}\left(\sqrt{B^{(a)}_{n,k}}\exp[-S^{(a)}_{n,k}]\right),$
(32) $\displaystyle
S^{(a)}_{n,k}=\frac{A^{(a)}_{n,k}}{2}\left[\left<\textbf{v}_{n,k}\middle|\delta\textbf{r}_{k-\frac{1}{2}}\right>^{2}+\left<\textbf{v}_{n,k}\middle|\delta\textbf{r}_{k+\frac{1}{2}}\right>^{2}\right]-B^{(a)}_{n,k}\left<\textbf{v}_{n,k}\middle|\delta\textbf{r}_{k-\frac{1}{2}}\right>\left<\textbf{v}_{n,k}\middle|\delta\textbf{r}_{k+\frac{1}{2}}\right>,$
(33)

where

$\displaystyle
V^{\prime\prime}(\textbf{r}_{k})\textbf{v}_{n,k}\equiv(\omega_{n,k})^{2}\textbf{v}_{n,k},\
\ \ (n=1,\cdots,2N),$ (34) $\displaystyle
A^{(a)}_{n,k}\equiv\frac{\omega_{n,k}}{\tanh(\omega_{n,k}\Delta\tau_{k})},\
B^{(a)}_{n,k}\equiv\frac{\omega_{n,k}}{\sinh(\omega_{n,k}\Delta\tau_{k})}.$
(35)

Eq. 34 defines normal mode frequencies and eigenmodes at each time slice $k$.
Note that at intermediate times $k\neq 0,M$, $\omega_{n,k}$ is in general
complex. At the end intervals $k=0,M$, one substitutes
$\delta\textbf{r}_{-\frac{1}{2}}\to\delta r_{0}=0$ and
$\delta\textbf{r}_{M+\frac{1}{2}}\to\delta r_{M}=0$ in the above expressions.
(Note that as $\tilde{\beta}\to\infty$, the propagators at the end intervals
approach zero exponentially. However, as we will see below, such contributions
cancel when calculating $A_{a}$ as we are calculating the ratio between two
$F$s.)

$F\big{[}\textbf{r}^{(a)}(\tau)\big{]}$ can finally be computed by integrating
over the intermediate fluctuation coordinates
$\delta\textbf{r}_{k-\frac{1}{2}}$

$\displaystyle F\big{[}\textbf{r}^{(a)}(\tau)\big{]}$
$\displaystyle=\int\left(\prod_{k=1}^{M}d^{2N}\delta\textbf{r}_{k-\frac{1}{2}}\right)\left<\textbf{0}\middle|e^{-\Delta\tau_{M}h_{M}}\middle|\delta\textbf{r}_{M-\frac{1}{2}}\right>\left<\delta\textbf{r}_{M-\frac{1}{2}}\middle|e^{-\Delta\tau_{M-1}h_{M-1}}\middle|\delta\textbf{r}_{M-\frac{3}{2}}\right>\cdots\left<\delta\textbf{r}_{\frac{1}{2}}\middle|e^{-\Delta\tau_{0}h_{0}}\middle|\textbf{0}\right>$
$\displaystyle=\left(\prod_{k=0}^{M}\prod_{n=1}^{2N}\sqrt{B_{n,k}^{(a)}}\right)\det(\mathcal{M}^{(a)})^{-\frac{1}{2}},$
(36) $\displaystyle\mathcal{M}^{(a)}$
$\displaystyle\equiv\sum_{k=1}^{M}e_{k,k}\otimes(\mathcal{A}_{k-1}^{(a)}+\mathcal{A}_{k}^{(a)})-\sum_{k=1}^{M-1}(e_{k,k+1}+e_{k+1,k})\otimes\mathcal{B}_{k}^{(a)}$
$\displaystyle=\begin{bmatrix}\mathcal{A}_{0}^{(a)}+\mathcal{A}_{1}^{(a)}&-\mathcal{B}_{1}^{(a)}&\textbf{0}&\cdots&\textbf{0}\\\
-\mathcal{B}_{1}^{(a)}&\mathcal{A}_{1}^{(a)}+\mathcal{A}_{2}^{(a)}&-\mathcal{B}_{2}^{(a)}&\cdots&\textbf{0}\\\
\textbf{0}&-\mathcal{B}_{2}^{(a)}&\mathcal{A}_{2}^{(a)}+\mathcal{A}_{3}^{(a)}&\ddots&\vdots\\\
\vdots&\vdots&\ddots&\ddots&-\mathcal{B}_{M-1}^{(a)}\\\
\textbf{0}&\textbf{0}&\cdots&-\mathcal{B}_{M-1}^{(a)}&\mathcal{A}_{M-1}^{(a)}+\mathcal{A}_{M}^{(a)}\end{bmatrix},$
(37) $\displaystyle\mathcal{A}_{k}^{(a)}$
$\displaystyle\equiv\sum_{n=1}^{2N}A_{n,k}^{(a)}\textbf{v}_{n,k}\textbf{v}_{n,k}^{\rm
T},\
\mathcal{B}_{k}^{(a)}\equiv\sum_{n=1}^{2N}B_{n,k}^{(a)}\textbf{v}_{n,k}\textbf{v}_{n,k}^{\rm
T}.$ (38)

Here $\mathcal{M}^{(a)}$ is a real symmetric block tridiagonal matrix,
$e_{i,j}$ is the $M\times M$ matrix with $1$ at the $(i,j)$-th entry with all
other entries $0$, $\otimes$ is the Kronecker product of two matrices and
$\mathcal{A}_{k}^{(a)}$ and $\mathcal{B}_{k}^{(a)}$ are $2N\times 2N$
matrices. [Note that in the present WC problem, one needs to project out two
zero eigenmodes $\textbf{v}_{n,k}$ (for each $k$) corresponding to uniform
translations in the $x$ and $y$ directions; hence $\mathcal{A}_{k}^{(a)}$ and
$\mathcal{B}_{k}^{(a)}$ become $(2N-2)\times(2N-2)$ matrices.]

In calculating $F[\textbf{r}_{0}]$—which essentially is the propagator of a
quantum harmonic oscillator—with the same procedure, one merely substitutes
$h_{k}\to h_{0}$ in every equation Eq. (32–38)

$\displaystyle
F[\textbf{r}_{0}]=\left(\prod_{k=0}^{M}\prod_{n=1}^{2N}\sqrt{B_{n,k}^{(0)}}\right)\det(\mathcal{M}^{(0)})^{-\frac{1}{2}}$
(39) $\displaystyle
A^{(0)}_{n,k}\equiv\frac{\omega_{n,0}}{\tanh(\omega_{n,0}\Delta\tau_{k})},\
\,B^{(0)}_{n,k}\equiv\frac{\omega_{n,0}}{\sinh(\omega_{n,0}\Delta\tau_{k})},$
(40)
$\displaystyle\mathcal{A}_{k}^{(0)}\equiv\sum_{n=1}^{2N}A_{n,k}^{(0)}\textbf{v}_{n,0}\textbf{v}_{n,0}^{\rm
T},\
\,\mathcal{B}_{k}^{(0)}\equiv\sum_{n=1}^{2N}B_{n,k}^{(0)}\textbf{v}_{n,0}\textbf{v}_{n,0}^{\rm
T},$ (41)
$\displaystyle\mathcal{M}^{(0)}\equiv\sum_{k=1}^{M}e_{k,k}\otimes(\mathcal{A}_{k-1}^{(0)}+\mathcal{A}_{k}^{(0)})$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
-\sum_{k=1}^{M-1}(e_{k,k+1}+e_{k+1,k})\otimes\mathcal{B}_{k}^{(0)}.$ (42)

Therefore,

$\displaystyle\frac{F\big{[}\textbf{r}^{(a)}(\tau)\big{]}}{F\big{[}\textbf{r}_{0}\big{]}}=\left(\prod_{k=1}^{M-1}\prod_{n=1}^{2N}{\frac{B_{n,k}^{(a)}}{B_{n,k}^{(0)}}}\right)^{\frac{1}{2}}\left[\frac{\det(\mathcal{M}^{(a)})}{\det(\mathcal{M}^{(0)})}\right]^{-\frac{1}{2}}.$
(43)

Here the product over $k$ runs only from $1$ to $M-1$ because the end interval
contributions ($k=0,M$) of $B_{n,k}^{(a)}$ and $B_{n,k}^{(0)}$ are identical
although they formally approach $0$ as $\tilde{\beta}\to\infty$ [since
$\Delta\tau_{0,M}\to\infty$]. Similarly, one takes
$A_{n,0}^{(a)}=A_{n,0}^{(0)}=A_{n,M}^{(a)}=A_{n,M}^{(0)}=\omega_{n,0}$ in
calculating $\det\mathcal{M}$, as $\Delta\tau\to\infty$ [since
$\tanh(\omega_{n,0}\Delta\tau)\to 1$].

Finally, one needs to divide Eq. 43 by the (formally diverging) zero mode
contribution to $F[\textbf{r}^{(a)}(\tau)]={\rm
det}\big{(}-\partial_{\tau}^{2}+V^{\prime\prime}\big{[}\textbf{r}^{(a)}(\tau)\big{]}\big{)}^{-1/2}$.
For this, we numerically find the smallest $\lambda$ such that

$\displaystyle\frac{1}{F_{\lambda}[\textbf{r}^{(a)}(\tau)]}\equiv\det(-\partial_{\tau}^{2}+V^{\prime\prime}\big{[}\textbf{r}^{(a)}(\tau)\big{]}-\lambda)^{\frac{1}{2}}=0,$
(44)

where the left hand side is calculated similarly as in Eqs. (32–38) with the
substitution $V^{\prime\prime}(\textbf{r}_{k})\to
V^{\prime\prime}(\textbf{r}_{k})-\lambda.$ The fluctuation determinant is then
obtained as

$\displaystyle
A_{a}=\sqrt{\lambda}\cdot\frac{F\big{[}\textbf{r}^{(a)}(\tau)\big{]}}{F\big{[}\textbf{r}_{0}\big{]}}.$
(45)

## References

  * Wigner [1934] E. Wigner, On the interaction of electrons in metals, Phys. Rev. 46, 1002 (1934).
  * Tanatar and Ceperley [1989] B. Tanatar and D. M. Ceperley, Ground state of the two-dimensional electron gas, Physical Review B 39, 5005 (1989).
  * Attaccalite  _et al._ [2002] C. Attaccalite, S. Moroni, P. Gori-Giorgi, and G. B. Bachelet, Correlation energy and spin polarization in the 2D electron gas, Physical Review Letters 88, 256601 (2002).
  * Drummond and Needs [2009] N. Drummond and R. Needs, Phase diagram of the low-density two-dimensional homogeneous electron gas, Physical Review Letters 102, 126402 (2009).
  * Hossain  _et al._ [2020] M. S. Hossain, M. Ma, K. V. Rosales, Y. Chung, L. Pfeiffer, K. West, K. Baldwin, and M. Shayegan, Observation of spontaneous ferromagnetism in a two-dimensional electron system, Proceedings of the National Academy of Sciences 117, 32244 (2020).
  * Hossain  _et al._ [2021] M. S. Hossain, M. Ma, K. Villegas-Rosales, Y. Chung, L. Pfeiffer, K. West, K. Baldwin, and M. Shayegan, Spontaneous valley polarization of itinerant electrons, Physical Review Letters 127, 116601 (2021).
  * Kim and Kivelson [2021] K.-S. Kim and S. A. Kivelson, Discovery of an insulating ferromagnetic phase of electrons in two dimensions, Proceedings of the National Academy of Sciences 118 (2021).
  * Falson  _et al._ [2022] J. Falson, I. Sodemann, B. Skinner, D. Tabrea, Y. Kozuka, A. Tsukazaki, M. Kawasaki, K. von Klitzing, and J. H. Smet, Competing correlated states around the zero-field Wigner crystallization transition of electrons in two dimensions, Nature Materials 21, 311 (2022).
  * Bernu  _et al._ [2001] B. Bernu, L. Cândido, and D. M. Ceperley, Exchange frequencies in the 2D Wigner crystal, Phys. Rev. Lett. 86, 870 (2001).
  * Roger [1984] M. Roger, Multiple exchange in ${}^{3}\mathrm{He}$ and in the Wigner solid, Phys. Rev. B 30, 6432 (1984).
  * Roger  _et al._ [1983] M. Roger, J. Hetherington, and J. Delrieu, Magnetism in solid ${}^{3}\mathrm{He}$ , Reviews of Modern Physics 55, 1 (1983).
  * Chakravarty  _et al._ [1999] S. Chakravarty, S. Kivelson, C. Nayak, and K. Voelker, Wigner glass, spin liquids and the metal-insulator transition, Philosophical Magazine B 79, 859 (1999).
  * Katano and Hirashima [2000] M. Katano and D. S. Hirashima, Multiple-spin exchange in a two-dimensional Wigner crystal, Phys. Rev. B 62, 2573 (2000).
  * Voelker and Chakravarty [2001] K. Voelker and S. Chakravarty, Multiparticle ring exchange in the Wigner glass and its possible relevance to strongly interacting two-dimensional electron systems in the presence of disorder, Physical Review B 64, 235125 (2001).
  * Kim  _et al._ [2022] K.-S. Kim, C. Murthy, A. Pandey, and S. A. Kivelson, Interstitial-induced ferromagnetism in a two-dimensional Wigner crystal, Physical Review Letters 129, 227202 (2022).
  * Fisher  _et al._ [1979] D. S. Fisher, B. I. Halperin, and R. Morf, Defects in the two-dimensional electron solid and implications for melting, Phys. Rev. B 20, 4692 (1979).
  * Cockayne and Elser [1991] E. Cockayne and V. Elser, Energetics of point defects in the two-dimensional Wigner crystal, Phys. Rev. B 43, 623 (1991).
  * Note [1] Note that the energetics of these point defects are also studied using the path integral Monte Carlo method in Ref. [53], but the contributions from defect hopping and exchange processes are ignored.
  * Note [2] When considering exchange processes in a pure WC, $V_{0}$ corresponds to the classical WC energy. For tunneling processes involving a defect, $V_{0}$ corresponds to the classical energy of the defect.
  * Thouless [1965] D. Thouless, Exchange in solid 3He and the Heisenberg Hamiltonian, Proceedings of the Physical Society (1958-1967) 86, 893 (1965).
  * Coleman [1988] S. Coleman, Aspects of symmetry: selected Erice lectures, Cambridge University Press (1988).
  * Altland and Simons [2010] A. Altland and B. D. Simons, Condensed matter field theory, Cambridge University Press (2010).
  * Zinn-Justin [2021] J. Zinn-Justin, Quantum field theory and critical phenomena, Oxford University Press (2021).
  * [24] The interaction between defects is small in proportion to a polynomial power of $1/d$, where $d$ is the distance between them, and can be neglected if the defect concentration is dilute enough.
  * Note [3] Such a peculiar defect dynamics has also been predicted to occur in solid helium [54].
  * Note [4] Such a “restricted mobility” would in turn lead to an exponentially slow (in $\sqrt{r_{s}}$) thermalization in the presence of a dilute vacancy concentration.
  * Eisenberg and Lieb [2002] E. Eisenberg and E. H. Lieb, Polarization of interacting bosons with spin, Physical review letters 89, 220403 (2002).
  * Yang and Li [2003] K. Yang and Y.-Q. Li, Rigorous proof of pseudospin ferromagnetism in two-component bosonic systems with component-independent interactions, International Journal of Modern Physics B 17, 1027 (2003).
  * Note [5] The potential importance of defect hopping processes in the magnetism of a quantum crystal was pointed out earlier in Refs. [20, 55].
  * Haerter and Shastry [2005] J. O. Haerter and B. S. Shastry, Kinetic antiferromagnetism in the triangular lattice, Physical Review Letters 95, 087202 (2005).
  * Sposetti  _et al._ [2014] C. N. Sposetti, B. Bravo, A. E. Trumper, C. J. Gazza, and L. O. Manuel, Classical antiferromagnetism in kinetically frustrated electronic models, Physical Review Letters 112, 187204 (2014).
  * Zhu  _et al._ [2022] Z. Zhu, D. Sheng, and A. Vishwanath, Doped Mott insulators in the triangular-lattice Hubbard model, Physical Review B 105, 205110 (2022).
  * Kim [2023] K.-S. Kim, Exact hole-induced resonating-valence-bond ground state in certain $U=\infty$ Hubbard models, Physical Review B 107, L140401 (2023).
  * 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, Simulation of Hubbard model physics in WSe2/WS2 Moiré superlattices, Nature 579, 353 (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, Correlated insulating states at fractional fillings of Moiré superlattices, Nature 587, 214 (2020).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang, Mott and generalized Wigner crystal states in WSe2/WS2 Moiré superlattices, Nature 579, 359 (2020).
  * Morales-Durán  _et al._ [2021] N. Morales-Durán, A. H. MacDonald, and P. Potasz, Metal-insulator transition in transition metal dichalcogenide heterobilayer Moiré superlattices, Phys. Rev. B 103, L241110 (2021).
  * Yang  _et al._ [2023] Y. Yang, M. Morales, and S. Zhang, Metal-insulator transition in transition metal dichalcogenide heterobilayer: accurate treatment of interaction (2023), arXiv:2306.14954 [cond-mat.str-el] .
  * Note [6] Such a particle-hole asymmetry was also pointed out in a related, but distinct, system in the presence of a Zeeman field [56].
  * Kim and Katsura [2023] K.-S. Kim and H. Katsura, Exact hole-induced $SU(N)$ flavor-singlets in certain $U=\infty$ $SU(N)$ Hubbard models, arXiv preprint arXiv:2307.09500 (2023).
  * Note [7] The enlarged $SU(4)$ symmetry of such a system is only an approximate symmetry in the limit of dilute filling of the bands. This approximate $SU(4)$ symmetry is broken by small band structure effects, which qualitatively change the low-energy spectrum. For example, the full $SU(4)$ polarization induced by interstitial hopping processes will result in $3$ magnon modes with $\omega\sim k^{2}$ low-energy spectrum. Taking into account small band structure effects, the symmetry is reduced to $SU(2)$, inducing a gap (proportional to the degree of symmetry breaking) to $2$ magnon modes. This in turn leads to a difference in the thermodynamic behavior at temperatures above and below the gap scale. Similar symmetry breaking effects in the low-energy spectrum would also be present in the antiferromagnetic phases.
  * Note [8] Our results for $C_{1,{\rm CI}}$ and $C_{3/2,{\rm CI}}$ agree with Ref. [17]. Note the factor $2$ difference between our values and those of Ref. [17] due to the difference in the unit of energy.
  * Note [9] Equation 21 is the energy of an interstitial in the fully polarized sector as compared to the fully polarized pure WC. As the WC is expected to be antiferromagnetic at $r_{s}<r_{\rm F}^{\rm wc}\approx 175$, the true interstitial energy should be calculated by comparing the energy of a single interstitial in some background antiferromagnetic WC to that of the antiferromagnetic WC without an interstitial.
  * Note [10] Our calculation is performed in a hexagonal supercell as opposed to that of Ref. [17] in a rectangular supercell.
  * Note [11] The possibility of a quantum crystal in which the particle density per crystalline unit cell is different from an integer value was originally conceived of in Ref. [57] in the context of solid helium.
  * Falakshahi and Waintal [2005] H. Falakshahi and X. Waintal, Hybrid phase at the quantum melting of the Wigner crystal, Physical Review Letters 94, 046801 (2005).
  * Arovas  _et al._ [2022] D. P. Arovas, E. Berg, S. A. Kivelson, and S. Raghu, The Hubbard model, Annual Review of Condensed Matter Physics 13, 239 (2022).
  * Emery  _et al._ [1990] V. Emery, S. Kivelson, and H. Lin, Phase separation in the $t-J$ model, Physical Review Letters 64, 475 (1990).
  * Eisenberg  _et al._ [2002] E. Eisenberg, R. Berkovits, D. A. Huse, and B. Altshuler, Breakdown of the Nagaoka phase in the two-dimensional $t-{J}$ model, Physical Review B 65, 134437 (2002).
  * Note [12] This argument ignores a subtlety regarding the formation of mesoscale density-modulated “microemulsion” phases in a (possibly narrow) regime about any putative first-order transition in the presence of Coulomb interactions [58, 59, 60, 61].
  * Note [13] This explanation for insulating ferromagnetism is distinct from, although related to, a previous proposal [15] based on WC-FL puddle formation.
  * Mogensen and Riseth [2018] P. K. Mogensen and A. N. Riseth, Optim: A mathematical optimization package for Julia, Journal of Open Source Software 3, 615 (2018).
  * Cândido  _et al._ [2001] L. Cândido, P. Phillips, and D. Ceperley, Single and paired point defects in a 2d Wigner crystal, Physical review letters 86, 492 (2001).
  * Andreev [1976] A. F. Andreev, Diffusion in quantum crystals, Soviet Physics Uspekhi 19, 137 (1976).
  * Spivak and Zhou [2000] B. Spivak and F. Zhou, Ferromagnetic correlations in quasi-one-dimensional conducting channels, Physical Review B 61, 16730 (2000).
  * Davydova  _et al._ [2023] M. Davydova, Y. Zhang, and L. Fu, Itinerant spin polaron and metallic ferromagnetism in semiconductor Moiré superlattices, Physical Review B 107, 224420 (2023).
  * Andreev and Lifshitz [1969] A. F. Andreev and I. M. Lifshitz, Quantum Theory of Defects in Crystals, Soviet Journal of Experimental and Theoretical Physics 29, 1107 (1969).
  * Spivak [2003] B. Spivak, Phase separation in the two-dimensional electron liquid in MOSFET’s, Physical Review B 67, 125205 (2003).
  * Spivak and Kivelson [2004] B. Spivak and S. A. Kivelson, Phases intermediate between a two-dimensional electron liquid and Wigner crystal, Physical Review B 70, 155114 (2004).
  * Jamei  _et al._ [2005] R. Jamei, S. Kivelson, and B. Spivak, Universal aspects of Coulomb-frustrated phase separation, Physical Review Letters 94, 056805 (2005).
  * Spivak and Kivelson [2006] B. Spivak and S. A. Kivelson, Transport in two dimensional electronic micro-emulsions, Annals of Physics 321, 2071 (2006).