Quantum Monte Carlo study of the phase diagram of the two-dimensional uniform electron liquid

Abstract

We present a study of spin-unpolarized and spin-polarized two-dimensional uniform electron liquids using variational and diffusion quantum Monte Carlo (VMC and DMC) methods with Slater-Jastrow-backflow trial wave functions. Ground-state VMC and DMC energies are obtained in the density range $1\leq r_{\text{s}}\leq 40$. Single-particle and many-body finite-size errors are corrected using canonical-ensemble twist-averaged boundary conditions and extrapolation of twist-averaged energies to the thermodynamic limit of infinite system size. System-size-dependent errors in Slater-Jastrow-backflow DMC energies caused by partially converged VMC energy minimization calculations are discussed. We find that, for $1\leq r_{\text{s}}\leq 5$, optimizing the backflow function at each twist lowers the twist-averaged DMC energy at finite system size. However, nonsystematic system-size-dependent effects remain in the DMC energies, which can be partially removed by extrapolation from multiple finite system sizes to infinite system size. We attribute these nonsystematic effects to the close competition between fluid and defected crystal phases at different system sizes at low density. The DMC energies in the thermodynamic limit are used to parameterize a local spin density approximation correlation functional for inhomogeneous electron systems. Our zero-temperature phase diagram shows a single transition from a paramagnetic fluid to a hexagonal Wigner crystal at $r_{\text{s}}=35(1)$, with no region of stability for a ferromagnetic fluid.

The two-dimensional (2D) uniform electron liquid (UEL) is an important model system in condensed matter physics and plasma physics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. The 2D UEL model consists of a system of electrons moving in 2D in a uniform, inert, neutralizing background, interacting via the Coulomb potential. 2D UELs are realized in numerous semiconductor heterostructures, e.g., a 2D UEL can be observed at the interface of metal-oxide-semiconductor structures [11], metal-oxide-semiconductor field-effect transistors [12, 13], quantum wells [14, 13], and MgZnO/ZnO heterostructures [15].

According to the 2D UEL phase diagram reported in Ref. 9, the ground state is a non-spin-polarized (paramagnetic) liquid from high density down to the point of Wigner crystallization. However, in the absence of crystallization, lowering the density would eventually make the paramagnetic fluid phase unstable with respect to a spin-polarized (ferromagnetic) fluid. An experimental study of nonequilibrium transport in low-density 2D electron systems at zero external magnetic field suggests that a fully spin polarized fluid is stable before crystallization [16], which was also reported by some previous theoretical studies [3, 8]. On the other hand, more recent experimental work has not found a region of stability for the ferromagnetic fluid [15]. In this work we update the real-space quantum Monte Carlo (QMC) calculations reported in Ref. 9 to include additional long-range backflow terms in trial wave functions and we look at a broader range of system sizes, resulting in a small modification to the crystallization density. Nevertheless, in agreement with the earlier work, our present QMC calculations yield no density range for the stability of fully or partially spin-polarized 2D UELs.

We have carried out QMC calculations for 2D UELs with spin polarizations of $\zeta=0$, 0.5, and 1 within the density range $1\leq r_{\text{s}}\leq 40$, where $r_{\text{s}}$ is the radius of the circle that contains one electron on average in units of the Bohr radius. We have focused on reducing finite-size (FS) errors [17, 18, 19, 20] by using twist-averaging (TA) and extrapolation techniques [21, 2, 22]. In this work we discuss errors in DMC calculations near the crystallization density introduced during wave function (WF) optimization, whose effects may be exaggerated by FS extrapolation.

Slater-Jastrow-backflow (SJB) WFs [23, 24] with plane wave orbitals $\exp(i\mathbf{G}\cdot\mathbf{r})$ were used in our QMC calculations for fluid phases. Details of the SJB WF have been explained in our recent studies [25, 26, 27]; the principle difference from Ref. 9 is the inclusion of a long-range two-body backflow term ${\bm{\pi}}$ that is expanded in a plane-wave basis. This term lowers the variational energy in finite cells, and is intended to make the treatment of two-body correlations in different simulation cells more consistent, to aid extrapolation to the thermodynamic limit. We used the variational and diffusion quantum Monte Carlo (VMC and DMC) methods [28] as implemented in the casino package [29]. To impose fermionic antisymmetry in our DMC calculations we have used the fixed-node approximation [30], where the nodal surface is forced to be the same as that of the trial WF. The variational parameters in the Jastrow factor and backflow functions were optimized within VMC by minimizing either the energy variance [31, 32] or the mean absolute deviation of the local energies from the median [33], then minimizing the energy expectation value [31, 34]. We used a target population of 2560 walkers in our DMC calculations. At $r_{\text{s}}=1$, 2, 5, 10, 15, 20, 25, 30, 35, and 40 we used DMC time steps of $\tau=0.005$, 0.01, 0.05, 1, 1.5, 2, 3, 4, 5, and 6 Ha^-1, respectively. Random errors due to TA are larger than the time-step bias in each case. Our DMC energies obtained at various system sizes using a hexagonal simulation cell are reported in the Supplemental Material [35].

Except where otherwise stated, our Jastrow factors and backflow functions were optimized at zero twist (i.e., pure periodic boundary conditions). For paramagnetic 2D UELs with densities $r_{\text{s}}=1$, 2, and 5, and system size $N=146$, we tested separately optimizing the Jastrow factor and backflow function at each twist. The results are listed in Table 1. We have used 30 random twists in our analysis. Our analysis shows that separately optimizing the Jastrow factor and backflow function at each twist reduces the TA DMC energy by $-0.225(9)\times 10^{-3}/r_{\text{s}}$ Ha/el. for paramagnetic 2D UELs at $N=146$. These results demonstrate the existence of significant, quasirandom FS errors in TA energies in which the same backflow function is used for all twists.

Table 1: TA DMC energies of paramagnetic 146-electron 2D UELs. The DMC calculations were performed using either the same Jastrow factor and backflow function, optimized at zero twist, for all twists (DMC₀) or SJB wave functions separately optimized at each twist (DMC₁).

$r_{\text{s}}$	TA DMC energy (Ha/el.)		Difference (Ha/el.)
	DMC0	DMC1
1	$-0.210731(7)$	$-0.210957(5)$	$-0.225(8)\times 10^{-3}$
2	$-0.258378(9)$	$-0.258501(6)$	$-0.12(1)\times 10^{-3}$
5	$-0.149595(4)$	$-0.149637(2)$	$-0.042(5)\times 10^{-3}$

WF optimization is challenging at large system sizes ($N\gtrsim 150$) and at densities near the solid-liquid phase transition ($r_{\text{s}}\gtrsim 20$) [35]. Our final WFs are obtained by energy minimization, but usually this requires a reasonable first approximation to the WF, obtained using a different method. Although unreweighted variance minimization exhibits superior numerical stability to energy minimization at high and intermediate densities, we found numerical instabilities in variance minimization for systems with large $N$ and low density. Numerical instability in variance minimization can occur when the nodes of the trial WF are altered during minimization. Near the transition density we find that numerical instabilities become more severe at large system size due to the increased likelihood of producing a wave function for the “wrong” phase (e.g., a defected crystal wave function instead of a fluid wave function) and the growth of complexity of the energy landscape in parameter space. Minimizing a more robust measure of the spread of the local energies, such as the mean absolute deviation (MAD) from the median local energy [33], provides a much lower variational energy (Table 2), because the MAD is less sensitive to divergent local energies [29].

We found that more than a dozen cycles of VMC configuration-generation and energy minimization are often required to converge adequately in the regime of large $N$ and large $r_{\text{s}}$, as shown in Fig. 1. The magnitude of the effect is significant on the scale of the FS error, even if the subsequent use of DMC weakens the dependence on the backflow function. A related observation is that the TA VMC and DMC energies per particle are often higher than suggested by extrapolation from smaller system sizes (see Figs. 2 and 3). In finite UELs of a given density and spin polarization, either Wigner crystal-like or Fermi fluid-like states may be energetically favored at different system sizes $N$, different simulation cell shapes, and different twists. Fluid energies fluctuate as a function of $N$ due to shell-filling effects in reciprocal space, while crystal energies fluctuate as a function of $N$ due to spatial commensurability effects in real space. Furthermore, there may be many near-degenerate low-energy states, e.g., due to different defected Wigner crystals in a particular cell. Even if a fluid trial wave function is used, it may be energetically advantageous for the system to approximate a floating Wigner crystal when the wave function is optimized. Conversely, at particular supercell sizes, it may be advantageous for a Wigner crystal wave function to try to approximate a fluid excited state. For certain simulation cells it may be especially difficult to find the lowest-energy fluid or crystal wave function. We can try to tip the balance in favor of our selected phase by using “magic” numbers of electrons (such that the occupied plane-wave orbitals form a closed-shell configuration in reciprocal space) for fluid phases, and by using square numbers of electrons for crystals, such that defect-free lattices are commensurate with the simulation cell. However, this is not guaranteed to prevent wave-function optimization giving an approximation to an excited state of the wrong phase. Although a Jastrow factor $\exp(J)$ cannot truly introduce nodes, if $J$ becomes large and negative then it may introduce “quasinodes” that affect the DMC energy over any reasonable imaginary time scale.

Figure 1: SJB VMC energy against energy minimization cycle for paramagnetic 2D UELs with density parameters $r_{\text{s}}=30$ (top panel) and 35 (bottom panel). The systems studied had $N=218$, 242, 254, and 302 electrons. The initial WF was optimized by minimizing the mean absolute deviation from the median local energy. The real WF was optimized at zero twist. The energies of the 2D UELs at $r_{\text{s}}=30$ were obtained using hexagonal and square simulation cells.

Figure 1 shows the energy minimization process for the paramagnetic liquid phase with system sizes $N=218$, 242, 254, and 302 and densities of $r_{\text{s}}=30$ and 35. The SJB WFs used for energy minimization were initially optimized by MAD minimization. The problem of slow convergence with optimization cycle can be observed in Fig. 1 and in the Supplemental Material [35]. These figures also show occasional jumps suggestive of wave-function parameters moving from one local minimum-energy configuration to another.

Residual canonical-ensemble TA errors in the Hartree–Fock kinetic and exchange energies around the crystallization density are orders of magnitude smaller than the fluctuations in the DMC energy shown in Fig. 3; hence residual momentum quantization effects in the TA energy cannot explain the nonsystematic system-size dependence. Other sources of nonsystematic FS error, such as Ruderman-Kittel oscillations being forced to be commensurate with the simulation cell, may be present in Fig. 2. Such nonsystematic FS errors are at least partially averaged out by extrapolation to the thermodynamic limit.

We investigated this issue further. First, we investigated how the initial SJB WF optimization depends on the number of configurations and the twist. The results are summarized in Table 2. They indicate that the optimized energy for a large system at low density does not depend strongly on the number of configurations or the twist wavevector, but does depend on the optimization method. Because the local energy diverges at nodes, the unweighted variance and MAD landscapes depend strongly on the sampled configurations. Nevertheless, it is clear that MAD minimization provides much better results than variance minimization and hence better starting points for energy minimization. Second, we compared the FS behavior with Slater-Jastrow (SJ) and SJB WFs. The fixed-node SJ-DMC energy is independent of the Jastrow factor. Figure 2 shows that the TA SJ-DMC energies behave in a nonsystematic manner as a function of system size, and that this behavior is further exaggerated by the inclusion of backflow. This demonstrates (i) that the nonsystematic behavior is not simply a result of difficulties optimizing wave functions and (ii) that applying analytic FS corrections [17, 18, 19] to results obtained at a single fixed cell size is unreliable because nonsystematic FS effects cannot be removed by this approach. Third, we examined the static structure (SF) factor and momentum-density (MD) of a paramagnetic 2D UEL with $N=254$ at $r_{\text{s}}=30$ obtained by VMC at a single twist. The results are presented in the Supplemental Material [35]. We have not found significant WF-optimization-dependent or system-size dependent anomalies in the SF and MD (e.g., the MD retains its discontinuity at the Fermi wavevector); unfortunately one cannot straightforwardly check for the wrong phase emerging, as the SJB wave function does not truly permit a change of phase. Fourth, we compared the energies of 2D UELs at $r_{\text{s}}=30$ obtained using hexagonal and square simulation cells with system sizes $62\leq N\leq 302$ [35]. The Wigner crystal energies used in our phase diagram (Fig. 4) were calculated using hexagonal simulation cells. Since the Madelung energy of a square lattice is higher, using a square lattice would have made crystal-like wave functions less competitive. Fifth, we performed SJ-VMC, SJ-DMC, and DMC with only the Slater trial wave function (S-DMC) calculations for paramagnetic Fermi fluids at $r_{\text{s}}=30$ in square cells with $N=90$ electrons. We chose three different initial wave functions for the optimization (Starting points 1–3) and optimized the Jastrow factor by energy minimization at either $\Gamma$ (${\bf k}_{\text{s}}={\bf 0}$) or the Baldereschi point [${\bf k}_{\text{s}}=(1/4){\bf b}_{1}+(1/4){\bf b}_{2}$, where ${\bf b}_{1}$ and ${\bf b}_{2}$ are the supercell reciprocal lattice vectors]. The DMC calculations were performed at time steps of 2 and 8 Ha^-1, with populations in inverse proportion to the time step, and the energies were extrapolated linearly to zero time step and infinite population. Hartree-Fock and VMC energies are shown in Table 3, while DMC energies with and without a Jastrow factor are shown in Table 4.

Table 2: Optimized energy of paramagnetic 2D UELs at $r_{\text{s}}=30$ with $N=302$ electrons. Only polynomial and plane-wave two-body Jastrow terms were used. The energies were calculated using different numbers of configurations $N_{\text{conf}}$. Two different optimization methods, unreweighted variance minimization (“varmin”) and MAD minimization (“madmin”), were used. Three twists, with fractional coordinates ${\bf k}_{0}=(0,0)$, ${\bf k}_{1}=(1/3,1/3)$, and ${\bf k}_{2}=(1/4,1/2)$, were used. All the energies are in Ha/el. The energies which are lower than $-0.031$ Ha/el. are highlighted in bold.

$N_{\text{conf}}$	Twist ${\bf k}_{0}$		Twist ${\bf k}_{1}$		Twist ${\bf k}_{2}$
	“varmin”	“madmin”	“varmin”	“madmin”	“varmin”	“madmin”
$1920$	$-0.027955(7)$	$-0.019491(8)$	$-0.025725(2)$	${\bf-0.0315169(4)}$	$-0.025009(2)$	${\bf-0.0314244(5)}$
$3840$	$-0.019497(8)$	$-0.019508(8)$	$-0.028892(8)$	$-0.019487(9)$	$-0.0307900(8)$	${\bf-0.0315429(4)}$
$7680$	$-0.0307833(9)$	${\bf-0.0314710(4)}$	$-0.01950(1)$	$-0.019496(9)$	$-0.029414(3)$	${\bf-0.0315166(4)}$
$9600$	$-0.025907(2)$	${\bf-0.0315087(4)}$	$-0.019502(9)$	$-0.019496(8)$	$-0.028948(2)$	$-0.01948(1)$
$11520$	$-0.019503(8)$	${\bf-0.0315034(4)}$	$-0.030772(3)$	${\bf-0.0315255(4)}$	$-0.0308548(6)$	$-0.01949(1)$
$15360$	$-0.019492(8)$	$-0.019523(9)$	${\bf-0.0313126(5)}$	$-0.01948(2)$	$-0.030773(1)$	$-0.019507(9)$
$19200$	$-0.027983(2)$	$-0.019519(8)$	$-0.030738(2)$	$-0.019494(9)$	$-0.01947(1)$	${\bf-0.0315176(5)}$

Figure 2: TA DMC energies for paramagnetic ($\zeta=0$) 2D UELs with density parameter $r_{\text{s}}=30$ obtained using SJ and SJB wave functions. The differences between the SJB and SJ energies for $N=218$, 242, 254, and 302 are $-0.041(2)$, $-0.042(4)$, $-0.068(1)$, and $-0.037(2)$ mHa/el., respectively.

Table 3: Hartree-Fock (HF) and SJ-VMC energies for 90-electron Fermi fluids in a square simulation cell, with the wave function being optimized from three different random starting points at the $\Gamma$ point and at the Baldereschi point of the square simulation cell.

${\bf k}_{\text{s}}$	Start	Energy (Ha/el.)
		HF	TA HF	SJ-VMC	TA SJ-VMC
$\Gamma$	1	$-0.0196518$	$-0.0196034$	$-0.0317373(2)$	$-0.0317344(2)$
$\Gamma$	2	$-0.0196518$	$-0.0196034$	$-0.0251398(8)$	$-0.025177(2)$
$\Gamma$	3	$-0.0196518$	$-0.0196034$	$-0.0271506(7)$	$-0.02716(1)$
Bald.	1	$-0.0194515$	$-0.0196034$	$-0.0317356(2)$	$-0.0317350(2)$
Bald.	2	$-0.0194515$	$-0.0196034$	$-0.02495(1)$	$-0.02500(1)$
Bald.	3	$-0.0194515$	$-0.0196034$	$-0.0277466(6)$	$-0.027726(7)$

Table 4: As Table 3, but using S-DMC and SJ-DMC calculations.

${\bf k}_{\text{s}}$	Start	Energy (Ha/el.)
		S-DMC	TA S-DMC	SJ-DMC	TA SJ-DMC
$\Gamma$	1	$-0.03151(2)$	$-0.03151(2)$	$-0.031919(1)$	$-0.031903(1)$
$\Gamma$	2	$-0.03151(2)$	$-0.03151(2)$	$-0.02676(7)$	$-0.02683(4)$
$\Gamma$	3	$-0.03151(2)$	$-0.03151(2)$	$-0.02848(2)$	$-0.02849(2)$
Bald.	1	$-0.03151(2)$	$-0.03159(2)$	$-0.031899(2)$	$-0.031898(1)$
Bald.	2	$-0.03151(2)$	$-0.03159(2)$	$-0.02670(2)$
Bald.	3	$-0.03151(2)$	$-0.03159(2)$	$-0.02901(2)$	$-0.02901(2)$

As expected, the choice of twist for optimization has only a small effect on the fluid energy at this low density. Indeed, twist averaging only has a small effect on the fluid energy. Starting points 2 and 3 produce wave functions that give a significantly lower VMC energy than the Hartree-Fock energy, i.e., they introduce correlations that lower the energy. But the energy obtained is not as low as the energy obtained in starting point 1. Hence there are issues with local minima.

The S-DMC energy obtained with the HF wave function is lower than the SJ-DMC energy obtained using the wave functions from starting points 2 and 3, despite the fact that the HF energy is higher than the SJ-VMC energies. Hence these Jastrow factor introduce features (“quasinodes”) that are difficult for DMC to remove (i.e., would require much smaller time steps). Since the HF wave function clearly corresponds to a fluid phase, it is reasonable to conclude that the SJ wave functions from starting points 2 and 3 do not correspond to a fluid phase. The S-DMC energy is in fairly good agreement with the SJ-DMC energy obtained with the good Jastrow factor obtained in starting point 1. The remaining small difference is very likely because of the need to extrapolate the S-DMC energy to zero time step a bit more carefully. Hence the Jastrow factor obtained in starting point 1 is very likely to describe a fluid phase.

Overall the picture is consistent with the idea that optimization of the wave function can get trapped in local minima that are not fluid-like, and DMC is unable to repair the situation when this happens. Some system sizes are more prone to problems with wave function optimization than others, due to the relative availability of crystal-like local minima. A disturbing feature is that one can optimize a wave function and lower the VMC energy, and yet the DMC energy is worse than not having a Jastrow factor at all. Having a more complex trial wave function is not a guarantee that the calculation is improved: it can make it easier to fall into unwanted local minima. All these conclusions apply more strongly in the case of SJB wave functions, which directly affect the DMC energy, even in the limit of zero time step.

Figure 3: Extrapolation of TA DMC fluid energies to the thermodynamic limit of infinite system size for paramagnetic ($\zeta=0$) 2D UELs with density parameter $r_{\text{s}}=30$. The red data points are assumed to be trapped in nonglobal minima during optimization and do not correspond to the fluid ground state.

We calculated the TA DMC energy for all the system sizes shown in Fig. 1. The energy at the thermodynamic limit is obtained by a linear extrapolation of the TA DMC energy as a function of $N^{-5/4}$ [10]. Figure 3 shows the extrapolation of TA DMC energies of two sets of data points labeled as “nonglobal minima” and “global minima.” All the TA DMC energies at the infinite system size limit presented in the rest of this paper are obtained using the extrapolation of TA DMC of systems that we believe correspond to the fluid ground state.

Table 5 lists our TA DMC energies at the thermodynamic limit for 2D UELs with spin polarizations $\zeta=0$, 0.5, and 1, and density parameters $1\leq r_{\text{s}}\leq 40$. We compare our results with the previous works of Drummond and Needs [9, 10]. The energies extrapolated to infinite system size are slightly lower than those obtained by Drummond and Needs. There are two major differences between our work and the previous ones. First, we used the two-body plane-wave backflow ${\bm{\pi}}$ term [26, 25], giving additional variational freedom in simulation cells of finite size. Second, we used larger system sizes and hexagonal simulation cells. These differences do not matter in principle, but in practice they affect the extrapolation to infinite system size.

Table 5: DMC energy in Ha/el. in the infinite system-size limit for 2D UELs with spin polarizations $\zeta=0$, 0.5, and 1 and density parameters $1\leq r_{\text{s}}\leq 40$. The results of Drummond and Needs are taken from Refs. 9 and 10.

$r_{\text{s}}$	Present work			Drummond and Needs
	$\zeta=0$	$\zeta=0.5$	$\zeta=1$	$\zeta=0$	$\zeta=1$
1	$-0.21017(6)$	…	$0.12626(4)$	$-0.2104(6)$	…
2	$-0.25846(7)$	…	$-0.19453(1)$	…	…
5	$-0.14998(2)$	…	$-0.143714(2)$	$-0.14963(3)$	…
10	$-0.085568(5)$	…	$-0.084585(6)$	$-0.085399(6)$	…
15	$-0.06138(2)$	…	$-0.059731(1)$	…	…
20	$-0.046388(1)$	$-0.0462813(1)$	$-0.046236(1)$	$-0.046305(4)$	$-0.046213(3)$
25	$-0.037807(3)$	$-0.0377575(1)$	$-0.037754(1)$	$-0.037774(2)$	$-0.037740(2)$
30	$-0.0319383(2)$	$-0.0319153(1)$	$-0.031919(2)$	$-0.031926(1)$	$-0.031913(1)$
35	$-0.0276718(4)$	…	$-0.027668(7)$	$-0.027665(1)$	$-0.027657(1)$
40	$-0.0244226(3)$	…	$-0.0244289(1)$	$-0.024416(1)$	$-0.024416(1)$

We used our DMC energies extrapolated to infinite system size to calculate the phase diagram of the 2D UEL (Fig. 4). The ferromagnetic and antiferromagnetic crystal energies used in our phase diagram calculation are taken from Ref. 9. We find that the paramagnetic Fermi fluid transitions to a hexagonal Wigner crystal at $r_{\text{s}}=35(1)$. The previous work by Drummond and Needs [9] predicted a value of $r_{\text{s}}=31(1)$ for this phase transition. Similar to the previous phase diagram predicted by Drummond and Needs [9], we find no region of stability for the polarized fluid phase with spin polarizations of $\zeta=0.5$ or 1.

Figure 4: DMC energy extrapolated to infinite system size as a function of $r_{\text{s}}$ for 2D UELs with spin polarizations $\zeta=0$, 0.5, and 1. Our results (“Prs-wrk”) are compared with those of Drummond and Needs [9], which include paramagnetic ($\zeta=0$) and ferromagnetic fluids ($\zeta=1$), and ferromagnetic and antiferromagnetic hexagonal Wigner crystals. The Madelung energy of a hexagonal lattice, $-v_{\text{M0}}/r_{\text{s}}$, where $v_{\text{M0}}=-0.50751467391482663$ is the Madelung constant at $r_{\text{s}}=1$, has been subtracted from all the DMC energies.

The correlation energy $E_{\text{c}}$ is approximately given by the difference between the TA DMC energy at the thermodynamic limit and the Hartree-Fock energy. Our correlation energy results are shown in Table 6.

Table 6: SJB-DMC correlation energy in the limit of infinite system size for 2D UELs with spin polarizations $\zeta=0$ and 1.

$r_{\text{s}}$	Correlation energy (Ha/el.)
	$\zeta=0$	$\zeta=1$
1	$-0.10996(6)$	$-0.02491(4)$
2	$-0.08335(7)$	$-0.02012(1)$
5	$-0.04993(2)$	$-0.013948(2)$
10	$-0.03054(5)$	$-0.009702(6)$
15	$-0.02359(2)$	$-0.007587(1)$
20	$-0.017627(1)$	$-0.006295(1)$
25	$-0.014598(3)$	$-0.005401(1)$
30	$-0.012487(2)$	$-0.004736(2)$
35	$-0.0109311(4)$	$-0.004232(7)$
40	$-0.0097298(3)$	$-0.0038332(1)$

We fit our TA DMC correlation energies to the Padé function [3]

$$E_{\text{c}}(r_{\text{s}},\zeta)=a_{\zeta}\times\frac{1+b_{\zeta}r_{\text{s}}^{1/2}}{1+b_{\zeta}r_{\text{s}}^{1/2}+c_{\zeta}r_{\text{s}}+d_{\zeta}r_{\text{s}}^{3/2}},$$ (1)

which has the correct asymptotic behavior at high and low densities [3], i.e., $E_{\text{c}}(r_{\text{s}})=a-acr_{\text{s}}+O(r_{\text{s}}^{3/2})$ at small $r_{\text{s}}$, while $E_{\text{c}}(r_{\text{s}}\rightarrow\infty)=(ab/d)r_{\text{s}}^{-1}+a(1/d-bc/d^{2})r_{\text{s}}^{-3/2}+O(r_{\text{s}}^{-2})$ at large $r_{\text{s}}$. The fitting parameters $a$, $b$, $c$, and $d$ are listed in Table 7.

Table 7: Fitting parameters for the SJB-DMC correlation energy [Eq. (1)] of paramagnetic ($\zeta=0$) and ferromagnetic ($\zeta=1$) 2D UELs.

Parameter	$\zeta=0$	$\zeta=1$
$a$ (Ha/el.)	$-0.1842$	$-0.0426$
$b$	$3.2320$	$28.8732$
$c$	$1.6027$	$16.5990$
$d$	$1.2541$	$4.6858$

In summary, we have used SJB WFs to perform VMC and DMC calculations for 2D UELs with spin polarizations $\zeta=0$, 0.5, and 1 within the density range $1\leq r_{\text{s}}\leq 40$. We have corrected single-particle and many-body FS errors by canonical-ensemble TA and extrapolation of the TA energies to infinite system size. We find that separately optimizing the Jastrow factor and backflow function at each twist improves the TA DMC energy for the density range $1\leq r_{\text{s}}\leq 5$ and possibly beyond. Optimization of the trial wave function is challenging near the transition density because the complexity of the energy landscape makes it difficult to find the ground state corresponding to the Fermi fluid. We predict that the paramagnetic 2D UEL transforms to a Wigner crystal at $r_{\text{s}}=35(1)$. We have not found any region of stability for 2D UELs with spin polarizations $\zeta=0.5$ or 1.