{'fraction_non_alphanumeric': 0.02585853408508457, 'fraction_numerical': 0.011814454126089185, 'mean_word_length': 4.50289098857707, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Diffusion quantum Monte Carlo study of three-dimensional Wigner crystals

We report diffusion quantum Monte Carlo calculations of three-dimensional Wigner crystals in the density range r_s=100-150. We have tested different types of orbital for use in the approximate wave functions but none improve upon the simple Gaussian form. The Gaussian exponents are optimized by directly minimizing the diffusion quantum Monte Carlo energy. We have carefully investigated and sought to minimize the potential biases in our Monte Carlo results. We conclude that the uniform electron gas undergoes a transition from a ferromagnetic fluid to a body-centered-cubic Wigner crystal at r_s=106+/-1. The diffusion quantum Monte Carlo results are compared with those from Hartree-Fock and Hartree theory in order to understand the role played by exchange and correlation in Wigner crystals. We also study"floating"Wigner crystals and give results for their pair-correlation functions.


I. INTRODUCTION
Electrons in a uniform potential are expected to crystallize at low densities to minimize their interaction energy. 1,2 In this paper we investigate Wigner crystals in three dimensions using quantum Monte Carlo (QMC) methods within the variational (VMC) and diffusion (DMC) approaches. The DMC method is currently the most accurate available for calculating the zero-temperature ground state energy of extended quantum mechanical systems.
There has been some debate about the density at which the transition from a Fermi fluid to a Wigner crystal should occur in three dimensions. In their pioneering DMC study of the phases of the electron gas, Ceperley and Alder 3 obtained a transition density of r s = 100 ± 20, 4 but a more recent DMC study gave r s = 65 ± 10. 5 Moreover, highly accurate DMC energies for the low-density fluid have recently become available 6 , which may further modify predictions of the transition density. The primary goal of this work is to provide highly accurate DMC energies for three-dimensional Wigner crystals and to use them in conjunction with the fluid data of Zong et al. 6 to predict an accurate value for the transition density.
To achieve sufficient accuracy we have carefully studied the possible sources of error in our calculations, including finite size effects, timestep errors and population control errors. We have used three procedures for optimizing the trial wave functions: minimization of the variance of the energy within VMC 7,8 ; minimization of the VMC energy; and minimization of the DMC energy.
We compare our results with recently-published fully-self-consistent unrestricted Hartree-Fock (HF) calculations 9 , and with the results of a simple version of Hartree theory, which allows us to understand the effects of exchange and correlation in Wigner crystals.
Finally, we discuss "floating" Wigner crystals in which the homogeneous and isotropic nature of the ground state is restored, relating their properties to those of a "fixed" crystal. DMC results for the pair correlation functions (PCFs) of floating Wigner crystals are compared with those of the fluid phases.

II. QUANTUM MONTE CARLO METHODS
We have performed both VMC and DMC calculations using the CASINO package. 10 In VMC, expectation values are calculated using an approximate many-body wave function, the integrals being performed by a Monte Carlo method. The approximate wave function normally contains a number of variable parameters, whose values are obtained by an optimization procedure.
In the DMC method 3,11 the imaginary time Schrödinger equation is used to evolve an ensemble of electronic configurations towards the ground state. The fermionic symmetry is maintained by the fixed-node approximation 12 , in which the nodal surface of the wave function is constrained to equal that of an approximate wave function. We will refer to the approximate wave functions used in VMC and DMC as trial wave functions. The fixed-node DMC energy provides a variational upper bound on the ground state energy with an error that is second order in the error in the nodal surface. 13,14 The trial wave function introduces importance sampling and controls both the statistical efficiency and the final accuracy that can be obtained. In DMC the final accuracy depends on the nodal surface of the trial wave function via the fixed-node approximation, while in VMC the final accuracy depends on the entire trial wave function, so that VMC energies are more sensitive to the quality of the trial wave function than DMC energies. Apart from the fixed-node approximation, DMC results may be subject to bias from the use of the short-time approximation (finite timestep errors), population control errors, and effects from the finite size of the simulation cell. We have made strenuous attempts to reduce these errors: see Sec. IV. The statistical errors in our QMC data are estimated using the blocking method 15 to eliminate the effects of serial correlation.

III. TRIAL WAVE FUNCTIONS
A. The Slater-Jastrow form We have used trial wave functions of the standard Slater-Jastrow form, where D ↑ and D ↓ are Slater determinants of up and down spin single-particle orbitals and R ↑ and R ↓ denote the coordinates of the up and down spin electrons, respectively, and e J is the Jastrow correlation factor.

B. Jastrow factors
We have used Jastrow factors of the form where with F ij = √ 2A if electrons i and j have the same spin and F ij = √ A if the electrons have opposite spins. This term satisfies the electron-electron cusp conditions. 11,16 The constant L 0 is set to 0.3 of the Wigner-Seitz radius of the simulation cell and A is a free parameter. The u 0 term is set to zero for r ij greater than the Wigner-Seitz radius, resulting in a small discontinuity in the Jastrow factor of less than 2 × 10 −5 in magnitude. To investigate possible bias from this discontinuity we compared the values of the two standard estimators of the kinetic energy involving the gradient and Laplacian of the trial wave function from a very long VMC run at r s = 100. The estimators agreed to within the statistical error, which was smaller than in our final DMC runs.
The second term in the Jastrow factor is given by where T l is the lth Chebyshev polynomial, L ′ is the Wigner-Seitz radius of the simulation cell and B ′ and the α l are parameters to be determined.

C. Orbitals for the Slater determinants
The Slater determinants for the crystalline phases were formed from localized non-orthogonal single-particle orbitals centered on the lattice sites of a body-centered cubic (BCC) crystal. A BCC crystal is expected to be favored in the low-density limit because it has the lowest Madelung energy. Throughout, we use φ(r) to denote a spatial orbital centered on the origin. Periodic orbitals for use in a simulation of a finite system subject to periodic boundary conditions are constructed for each lattice point in the simulation cell by summing over all the replicas of φ centered on that lattice point. Clearly, if all the individual orbitals are periodic then their Slater determinant is too. We use a Jastrow factor containing only homogeneous terms, see Sec. III B, with the differences between electron coordinates being evaluated under the minimum image convention. Hence the overall wave function is periodic.
In their VMC and DMC studies of Wigner crystals, Ceperley 17 and Ceperley and Alder 3 used Gaussian orbitals: They determined C by variational methods, with the Jastrow factor being optimized simultaneously. Ortiz et al. 5 used exponentials of two-parameter Padé functions: and determined the values of C and k and the parameters in their Jastrow factor by minimizing the variance of the energy within VMC. We have also investigated two new types of orbital for Wigner crystals. A straightforward generalization of Eq. (5) is to use a linear combination of Gaussian orbitals: We have also considered orbitals based on an expansion in the eigenstates of a simple harmonic oscillator. An orbital is constructed by multiplying the simple Gaussian function by a polynomial. For a BCC lattice with identical orbitals on each site the polynomial should have the full symmetry of the lattice, i.e., φ(r) = exp(−Cr 2 ) 1 + αr 2 + βr 4 + γ(x 2 y 2 + x 2 z 2 + y 2 z 2 ) + O(r 6 ) , where r = (x, y, z). This orbital has considerable flexibility at small r.

D. Optimization of the trial wave functions
Parameters in the trial wave functions may be optimized in a variety of ways. In principle the DMC energy depends only on the nodal surface of the wave function, which is determined by the form of the orbital. It is therefore best to minimize the DMC energy directly with respect to the parameters in the orbitals, but this is a costly and laborious procedure which we have carried out only for the Gaussian parameter C of the simple Gaussian orbital. In principle the DMC energy does not depend on the Jastrow factor, so it cannot be optimized in this fashion.
We first studied the simple Gaussian and Padé forms of Eqs. (5) and (6), using energy variance minimization to optimize the C parameter and the parameters in Jastrow factor, and minimization of the VMC energy with respect to variations in k, but we found the optimal value of k to be very close to zero. The Padé orbital seems to offer little advantage over a simple Gaussian orbital at the densities studied (100 ≤ r s ≤ 150).
We optimized the expansion coefficients and Gaussian exponents of the linear combination of Gaussian orbitals form of Eq. (7), together with a Jastrow factor, using VMC energy variance minimization. However, again, it was found that in practice this orbital offers no advantage over a single Gaussian function at the densities studied.
We attempted to optimize the α parameter in the harmonic oscillator form of Eq. (8) for a Wigner crystal at r s = 100, but the resulting wave function gave a DMC energy within the error bars of the one obtained by setting α = 0.
We have therefore adopted the simple Gaussian orbital of Eq. (5) for our main calculations. We have adopted the following procedure to optimize the trial wave functions. The A parameter in the Jastrow factor was optimized by minimization of the VMC energy, the parameters in the S 1 part of the Jastrow factor by energy variance minimization, and the C parameter in the Gaussian orbitals by minimization of the DMC energy. These minimizations were performed in turn until the changes in the parameters were negligible. We found that the DMC-optimized exponents obey C ≈ 0.11r −3/2 s , which gives rather smaller values than those used by Ceperley 3,17 (C ≈ 0.2r −3/2 s ) and considerably smaller than those predicted by HF theory or the simple Hartree model (C ≈ 0.5r −3/2 s ), see Sec. VII.

IV. ACCURACY OF THE DMC RESULTS
Our DMC algorithm is essentially that of Umrigar et al. 18 Here we explore the sources of error in our DMC calculations and justify our choices of the parameters for the final production runs, which are summarized in Sec. IV D. Unless otherwise stated, we use simple Gaussian orbitals throughout this section.

A. Finite size effects
We used periodic boundary conditions and the Ewald interaction energy to reduce the finite size effects. We tested the convergence of the Ewald sums and found that truncation errors were less than 10 −3 of the statistical error bars on the final DMC runs.
The energy per electron at a given density depends on the number of electrons in the simulation cell. We wish to obtain the energy per electron in the limit that the number of electrons per simulation cell goes to infinity. Two approaches have been used previously when dealing with finite size effects in QMC simulations of Wigner crystals. Ortiz et al. 5 used large simulation cells and found the finite size errors to be less than their statistical error bars of 8.5 × 10 −6 a.u. per electron for numbers of electrons in excess of 500. Ceperley 3,17 , on the other hand, used smaller system sizes in conjunction with an extrapolation scheme.
Because we wish to work to very high accuracy, we use quite large simulations cells and the extrapolation formula derived by Ceperley 17 : where the constants E ∞ and E N are the total energies per electron of the infinite system and the system with N electrons, and b is independent of both r s and N .
Starting from a two-parameter Padé orbital and corresponding Jastrow factor optimized in a 64-electron simulation cell, we attempted to further optimize the wave function in a 216-electron unit cell. This attempt did not lead to a lowering of the VMC energy, suggesting that the 64-electron simulation cell is adequate for optimization purposes; this size of cell was used in all of our subsequent optimization runs.

B. Population control biasing
The use of a finite population of configurations results in a positive bias in the DMC energy which, it is argued, falls off as the reciprocal of the target population. 18 This turns out to be a genuine problem in the case of Wigner crystals where we are able to work to extremely high precision. An example of the problem of population control biasing is shown in Fig. 1, where it can be seen that the bias is indeed positive and that it falls off roughly as the reciprocal of the target population. The simplest method for avoiding population control biasing is to use a large target population. Alternatively, the reweighting scheme developed by Umrigar et al. 18 can be used. In our tests we found this scheme to work well, provided the number of reweighting factors was about the same as the number of timesteps over which average local energies are correlated (between 100 and 1000 timesteps of 30 a.u.). However, with this many reweighting factors present, the total weight at each timestep fluctuated enormously and very long simulations were required in order to obtain meaningful statistics. We found it to be more efficient to use larger populations than to employ the reweighting scheme. For this reason, we did not use the reweighting scheme in our production runs.
Including more parameters in the Jastrow factor can lead to a significant reduction in population control biasing, as illustrated in Fig. 2. This shows that, when the reweighting scheme is not used, DMC energies obtained with a a poor Jastrow factor and a small target population (solid line) are too high, but when a larger population is used in conjunction with the same Jastrow factor (dashed line) the energies are similar to those obtained with a good Jastrow factor and a small population (dotted line). Improving the overall quality of the trial wave function reduces the population control bias because it reduces the fluctuations in the population. The results shown in Fig. 1 in which the Jastrow factor consisted only of the u 0 term therefore represent a worst-case scenario. configurations; dashed line: DMC energies for the same Jastrow factor, but with a target population of 800 configurations; dotted line: DMC energies obtained with a fully optimized Jastrow factor (with both u0 and S1 terms) and a target population of 100 configurations.

C. Timestep biasing
The variation of DMC energy with timestep is shown in Fig. 3 for three different Jastrow factors. It can be seen that the bias is always positive and that it grows linearly with timestep; therefore we can largely eliminate the bias by carrying out simulations at a number of different timesteps and performing a linear extrapolation to zero timestep.
The differences between the DMC energies in Fig. 3 are due to population control and timestep bias. The solid and dashed curves were obtained using the same target population of 100 configurations, but with different Jastrow factors. The population control bias at fixed timestep is clearly smaller for the Jastrow factor with A = 438.389 (solid line) than for A = 563.157 (dashed line). The dotted line was obtained with a target population of 800 configurations, which essentially removes the population control bias. Because the solid and dotted curves are approximately parallel we deduce that the population and timestep errors are approximately independent of one another. Furthermore, it is clear that altering the Jastrow factor has a considerably greater effect on the population control bias than on the timestep bias.
In Fig. 4 we show that timestep bias remains a problem even with a well-optimized Jastrow factor and a large target population which essentially removes the population bias. However, this figure again shows that a linear fit is appropriate when extrapolating to zero timestep.

D. Parameters for the production runs
The final production runs were characterized as follows.
1. The Gaussian exponents in the orbitals were optimized by minimizing the DMC energy, the A parameters by minimizing the VMC energy, and the other parameters in the Jastrow factors by minimizing the variance of the 3. At each density, DMC calculations were performed using between four and six different timesteps and the energy was extrapolated linearly to zero timestep. 4. A variety of system sizes were used (see Table II), and the energies were extrapolated to infinite system size using Eq. (9).

A. DMC energies obtained using DMC-optimized Gaussian orbitals
The values of the exponents obtained by optimizing the DMC energy are shown in Table I, along with the final DMC energies. These were found by using the DMC results shown in Table II in conjunction with the finite size extrapolation formula of Eq. (9). Using the results for r s = 100 and r s = 125, we obtained a good fit, giving b = 1.26(3).   In Fig. 5 we plot the electronic charge densities for a r s = 100 crystal, calculated using HF theory orbitals (see Sec. VII), from the DMC-optimized orbitals but without a Jastrow factor, and within DMC using an optimized Jastrow factor.
The HF theory orbitals are very localized whereas the orbitals optimized within DMC are much more diffuse. However, the inclusion of a Jastrow factor results in the charge density from the DMC-optimized orbitals becoming more localized on lattice sites, although not to the same extent as the HF orbitals. The VMC charge density obtained with the optimized Slater-Jastrow wave function is very close to the DMC density shown in Fig. 5. The peak difference between the DMC and VMC charge densities is 2.6 arb. units, and therefore an extrapolated estimation of the QMC charge density, ρ(r) ≈ 2ρ DMC − ρ VMC , would be very similar to the DMC density.
The Jastrow factor in a Wigner crystal serves to further localize the electrons on their lattice sites. Therefore, whereas HF theory gives very localized Gaussian orbitals, in an optimized Slater-Jastrow wave function we find the orbitals to be much more diffuse, with the localization being caused by correlation effects from the Jastrow factor. Correlation effects allow electrons to invade each other's space to some extent without incurring a high potential energy penalty; hence the overall charge density is less localized than in HF theory, and the kinetic energy is lower as a result.

VI. LOCATING THE FLUID-TO-CRYSTAL TRANSITION DENSITY
In order to locate the density at which the transition from the Fermi fluid to the crystal occurs, we fit the DMC energy data for these phases to interpolating functions. The correlation energy is defined as where E is the total energy, E HF is the HF ground state energy, ζ is the polarization. Ceperley 17 proposed a fitting form for the correlation energy of a Fermi fluid as a function of r s , where γ ζ , β ζ 1 and β ζ 2 are fitting parameters. We make use of the highly accurate DMC energies of Zong et al. 6 for the ferromagnetic fluid phase at low densities, which used trial wave functions including "backflow" effects. We found an excellent fit to Eq. (11) giving γ 1 = −0.09399, β 1 1 = 1.5268 and β 1 2 = 0.28882. We also tried fitting the fluid data to the form proposed by Perdew and Zunger 19 , which is based upon Eq. (11), but includes an assumed dependence on polarization, so that the partially polarized and unpolarized data of Zong et al. could be used. The χ 2 value of this fit was not so good, however.
At low densities the total energy of a crystal phase can be expanded as where the {f } are constants. 20 The first of these is taken to be f 0 = −0.89593 in order to give the Madelung energy in the low density limit. We found that our DMC data fitted Eq. (12) very well, giving f 1 = 1.3379 and f 2 = −0.55270. These values are in reasonable agreement with the parameters found using a completely different method by Carr 20 and Carr, Coldwell-Horsfall and Fein 21 , who have calculated the zero-point lattice-vibrational energy of a Wigner crystal in order to give an analytical result of f 1 = 1.325. Furthermore, they use perturbation theory to obtain the approximate result f 2 = −0.365. This phonon model is in good agreement with our Wigner crystal energies at large r s , but it gives energies which are too high at smaller r s . The energies of the ferromagnetic fluid and BCC crystalline phases at low densities calculated by different authors are shown in Fig. 6. We found the transition from the fluid to crystalline phases to occur at r s = 106 ± 1, in agreement with the original result of Ceperley and Alder. 3 Note, however, that the transition density predicted using the fluid data of Zong et al. 6 in conjunction with the crystal data of Ceperley and Alder would be somewhat lower, at about r s = 127. Our Wigner crystal energies are slightly lower than those of Ceperley and Alder, even though they studied a Bose crystal, which must have a lower energy than the corresponding fermion one. We believe this difference must be due to some small bias in the results of Ceperley and Alder. 3 Fitting Eq. (12) to the crystal data of Ceperley and Alder, we find that f 1 = 1.4309 and f 2 = −1.1058. The discrepancy with the analytical result of Carr for f 1 is consistent with the presence of a small, positive, systematic bias in the crystal results of Ceperley and Alder.
Our transition density is considerably lower than the value of r s = 65 ± 10 obtained by Ortiz et al. 5,22 . The statistical error bars on their data are much larger than on the fluid data of Zong et al. 6 or on our Wigner crystal data, which hampers detailed comparisons. However, it appears that the main reason for the discrepancy is that  21 Where error bars on the DMC data cannot be seen, it is because they are smaller than the symbols showing the data points.
For the crystal phase, we may estimate the fixed node error resulting from the use of the HF orbitals by taking the difference between DMC energies calculated using the HF and DMC-optimized orbitals. At r s = 100 we find this difference to be 8.9(1) × 10 −5 a.u. per electron. Zong et al. 23 have calculated the fixed node errors resulting from the use of plane-wave orbitals for the fluid phases as 1.87(5) × 10 −5 a.u. per electron for the unpolarized fluid and 0.84(7) × 10 −5 a.u. per electron for the fully polarized fluid. Therefore, although the correlation energy of crystals is smaller than fluids at the same density, the fixed node errors resulting from the use of HF orbitals are considerably larger in crystals than in fluids.

VII. HF AND OTHER SIMPLE THEORIES
HF theory is described as "restricted" when the spin-orbitals are products of space and spin parts which are occupied in pairs with identical space parts, and "unrestricted" when the space parts are different or when they are not occupied in pairs. Within the quantum chemistry community the standard definition of the correlation energy is the difference between the exact and restricted HF energies, but in a Wigner crystal the electrons are localized in space individually and a description within restricted HF theory is not possible. We therefore define the correlation energy as the difference between the exact and unrestricted HF energies.
A recent self-consistent unrestricted HF study of electrons in a uniform potential gave stable Wigner crystal solutions for r s ≥ 4.5 in three dimensions. 9 Here we develop simplified versions of the HF model which almost exactly reproduce the results of the fully self-consistent HF studies 9 at low densities and compare the results with our DMC ones.
In the low-density regime, the overlap between orbitals centered on neighboring lattice sites is small and therefore we expect the Hartree and HF energies to be similar. Let us take the wave function of the crystalline phase to be a Hartree product of normalized Gaussian orbitals, centered on lattice sites. The Gaussian exponent C is to be determined variationally. The resulting kinetic energy per electron is easily evaluated as T = 3C/2. The spatial charge density is simply the superposition of the Gaussian charge densities due to each orbital. The electrostatic energy E H of this charge distribution may be readily evaluated, but we must subtract off the self-energies E 0 of each Gaussian. The total energy per electron is therefore given by where the G are the reciprocal lattice vectors and Ω = 4πr 3 s /3 is the volume of the primitive unit cell. Differentiating this energy with respect to the Gaussian exponent and approximating the sum by an integral, we find that The approximation is valid for large r s , where the density of reciprocal lattice vectors is large. Hence we find the optimal value of C to be C = 1/2r 3/2 s , which is precisely the result obtained by Wigner 2 using a spherical approximation.
Integrating Eq. (15) with respect to C we obtain the energy E = 3C/2 + π/2ΩC + f (r s ), where the function f (r s ) is the "constant" of integration. Inserting the optimal value of C and making use of the fact that, in the limit of low densities, E must tend to the Madelung energy of the crystal lattice, we find where M is the Madelung constant of the lattice. The simple Hartree model of Eq. (14) and the further simplifying approximation of Eq. (16) give energies which are very close to the full HF results for r s > 50; for example, the energy of Eq. (14) agrees with that of the full HF result to within 0.006%, whereas the energies at r s = 100 agree to within 0.001%. The agreement between the fully self-consistent HF data and the simpler approximations is just as good for a face centered lattice. These results demonstrate that exchange energies between orbitals on different sites are extremely small. The exchange interactions are only significant between nearest neighbor Gaussian orbitals and the exchange energy per electron is given by E X = −(N n /2) C/π exp(−CR 2 ), where N n is the number of nearest neighbors and R is the nearest neighbor distance. This expression gives extremely small energies for the low densities studied here.
We have also calculated DMC energies using the orbitals obtained from HF theory 9 , and Jastrow factors optimized using energy variance minimization. The results are summarized in Table III. The extrapolation of the DMC results to infinite system size was carried out using Ceperley's extrapolation scheme with his value of b = 1.5. 17 Comparing the data in Table III with that in Tables I or II, we see that the DMC energy obtained at r s = 100 with the HF orbitals is significantly higher than that obtained with DMC optimized orbitals. Some population control bias is present in the results in Table III, but we estimate this to be about 2% of the difference between the r s = 100 DMC energies obtained with the HF and DMC-optimized orbitals.  The strength of correlations in a system may be measured by the ratio of the correlation energy to the total energy, E c /E. The DMC results of Zong et al. 6 indicate that in the fluid phases E c /E tends to a positive constant as r s → ∞, while for the Wigner crystal our results show that E c /E tends to zero as r s → ∞. In this sense one may think of Wigner crystals as being weakly correlated systems at low densities.

VIII. MAGNETIC BEHAVIOR OF THE CRYSTALLINE PHASES
The tiny energy differences between ferromagnetic and antiferromagnetic crystals proved to be too small to resolve in our QMC calculations. It might be possible to resolve them using a correlated sampling approach within VMC. Such an approach should provide an accurate value for the energy difference between two systems, 1 and 2, if |Ψ 1 | 2 ≃ |Ψ 2 | 2 throughout the configuration space, which should hold for ferromagnetic and antiferromagnetic crystalline phases at sufficiently low densities. HF theory predicts ferromagnetic behavior in the low density limit but according to the theory of Thouless 24 such a system should be antiferromagnetic. In their path integral QMC calculations, Candido, Bernu, and Ceperley 25 have indeed found antiferromagnetic behavior for BCC Wigner crystals at low densities, although the energy differences are much smaller than our statistical error bars.

IX. FLOATING WIGNER CRYSTALS
The Hamiltonian of the uniform electron gas is invariant under the simultaneous translation or rotation of the electron positions. However, our trial Wigner crystal wave functions break these symmetries and, for example, the resulting charge densities are inhomogeneous, see Fig. 5. These wave functions represent Wigner crystals which are "pinned" by some small external influence. Pinning of Wigner crystals may arise from the presence of impurities or boundaries, and therefore the broken symmetry solutions are physically relevant, but Wigner crystals may also be mobile, in which case it is better to describe them as floating crystals.
One way of restoring the homogeneous and isotropic nature of the trial function is to consider a linear combination of displaced and rotated copies of the fixed wave function Ψ. This gives a "floating" wave function, where R(Ω) represents a rotation of all of the electron positions and the integrals are over all possible solid angles Ω and displacements ∆. Ψ F gives rise to a homogeneous and isotropic charge density n(r) = N/Ω, where N is the number of electrons and Ω is the volume of the system. Wave functions for floating Wigner crystals in extended systems have been discussed briefly by, for example, Bishop and Lührmann 26 , and by Mikhailov and Ziegler. 27 Rather more attention has been devoted to restoring the (rotational) symmetry in two-dimensional models of quantum dots. 27,28 In finite systems the energy gain per electron from restoring the symmetry can be substantial but for an infinite system it turns out to be negligible.
Using a trial function where the single particle orbitals in the Slater determinant are φ = e −Cr 2 , we have obtained analytical results at the variational level for the case when the translational symmetry is restored. We found that the difference in total energy is equal to the kinetic energy of the center of mass of the fixed crystal (3C/2), making the energy difference per electron negligible (this result also holds when a translationally invariant Jastrow factor is included). We also found that the expectation value of any operator that only depends on relative electron coordinates is the same for Ψ and Ψ F at the variational level. We expect the same conclusions to be true in DMC as well; in particular, we found that under open boundary conditions the nodal surface of the fixed and the translationally averaged trial functions are identical.
Despite the similarities in their energies, there is an important qualitative difference between the fixed and floating Wigner crystal wave functions. The fixed wave function can be written as a sum of disconnected partial wave functions, in the form Ψ = M ψ M , where the overlap between the ψ M tend exponentially to zero as N → ∞. From this it follows that it represents an insulator. 29,30 On the other hand, the same is not true for the floating wave function, resulting in conducting behavior.

A. Definitions
The spin-dependent PCF is defined as It would be very costly to evaluate the full six-dimensional function g σ,σ ′ (r, r ′ ) for a Wigner crystal within QMC. However, for a homogeneous and isotropic system g depends only on the separation between electrons, r = |r − r ′ |, so that where we have used n σ (r) = N σ /Ω. This one-dimensional function is much less costly to evaluate accurately than Eq. (18), hence we can obtain the PCF for a floating Wigner crystal. Since the expectation value of any operator that depends only on the difference between electron coordinates is the same for fixed and floating crystals, see Sec. IX, we can evaluate g σ,σ ′ (r) using our fixed trial function. Furthermore, the expectation value of the potential energy of the system can be written in terms of the spinindependent PCF, g(r) = 1 4 σ,σ ′ g σ,σ ′ (r): which holds even for an inhomogeneous system such as the fixed Wigner crystal. Studying the PCF of a floating crystal therefore provides insight into the physics of the fixed crystal as well.

B. Discussion and Results for the PCFs
We evaluated Eq. (19) within VMC and DMC by accumulating g σ,σ ′ (r) in radial bins. Our best estimates of g were obtained using the extrapolated estimator g(r) ≈ 2g DMC − g VMC . We calculated g σ,σ ′ (r) at r s = 110 for a ferromagnetic and an antiferromagnetic Wigner crystal using our optimized trial wave functions. For comparison, we have also calculated g σ,σ ′ (r) for the unpolarized fluid phase at r s = 110, using a trial wave function consisting of determinants of plane waves and a Jastrow factor optimized using energy variance minimization.
All of the biasing effects that apply to the DMC energy may also affect the PCFs. We found that this was indeed the case, as we have obtained results that showed very small but statistically significant differences when, for example, different timesteps were used or when different Jastrow factors were used in conjunction with small population sizes. Unlike the energy, it was not possible to use an extrapolation scheme to remove the timestep bias, as it showed no clear pattern. We found that finite size effects in the PCFs were very small, however, as results obtained for 64 electrons were not significantly different from those obtained with 512 electrons.
A further source of bias, which does not apply to the energy, arises from the use of the extrapolated estimator. Fig. 7 shows that the extrapolation can make a significant difference, and to check the reliability of this method we have evaluated the PCF using different quality Jastrow factors. As might be expected, the VMC and DMC results varied significantly, but the final extrapolated results were almost independent of the Jastrow factor, the differences being only slightly larger than the statistical error. This source of bias is therefore small and on a comparable level to the others.
The final PCF calculations were carried out in a 512 (518) electron system for the Wigner crystal (fluid), using a timestep of 30 a.u. and a target population of 960. Fig. 8 shows the extrapolated spin-independent PCF for the antiferromagnetic crystal and unpolarized fluid phases at r s = 110. For comparison, the HF result is also shown, together with the result obtained from the DMC-optimized orbitals, but without the Jastrow factor. The spinindependent PCF for the ferromagnetic crystal is not included as it was found to be almost identical to that of the antiferromagnetic crystal. The HF orbitals are very strongly localized and give the most rapidly varying PCF, whereas the PCF from the more diffuse DMC-optimized orbitals is much smoother.
It is interesting that the extrapolated PCF of the floating crystal shows strong oscillations at distances much larger than r s . This can be understood in terms of Eq. (20). The potential energies of the fixed and floating crystals are the same, but their charge densities and PCFs are quite different. For the floating crystal n(r) is constant while g(r, r ′ ) oscillates strongly, whereas for a fixed crystal the charge density n(r) oscillates and g(r, r ′ ) is expected to be much smoother.
Another interesting point is the large difference between the extrapolated parallel/antiparallel-spin PCFs of the crystal and the fluid, see Fig. 9. For the crystal, the PCF strongly reflects the underlying crystal structure composed of alternating spins, whereas for the fluid phase the two are almost identical. The energies of the two systems are, however, almost the same.

XI. CONCLUSIONS
We have carried out a careful DMC study of the BCC Wigner crystal in the density range 100 ≤ r s ≤ 150. We have experimented with several types of orbital in the trial wave functions but have been unable to improve upon the Gaussian form used in previous work. We have, however, optimized the Gaussian exponent by directly minimizing the DMC energy, which reduces the fixed node errors. We have also taken care to eliminate other biases in our DMC simulations, particularly those from timestep errors, population control bias, and finite size effects. We estimate that the uniform electron gas undergoes a transition from a ferromagnetic fluid to a BCC crystal at r s = 106 ± 1.
We have used Slater-Jastrow-type trial wave functions for our studies. Multiplying the Slater determinant by a pairwise repulsive Jastrow factor makes the charge density more inhomogeneous because the electrons in Wigner crystals are localized in space individually. This behavior contrasts with that found in many other systems where a pairwise repulsive Jastrow factor tends to smooth out the charge density.
The results of HF theory and Hartree theory are very similar because the exchange interaction between orbitals on different sites is small. The orbitals obtained within unrestricted HF theory (and Hartree theory) are very strongly localized and the kinetic energy within HF theory is larger than in our DMC calculations with the fully optimized trial wave functions. We have defined the correlation energy to be the difference between the exact and unrestricted HF energies. With this definition, and in the low density limit, Wigner crystals are weakly correlated systems. The inclusion of correlation in a Wigner crystal wave function beyond the unrestricted HF level results in the electronic charge density spreading out from the lattice sites. In this sense correlation delocalizes the electrons. Although the correlation energies of the crystal phases are smaller than those of the fluid phases at the same density, the use of HF orbitals within the trial wave functions results in larger fixed node errors for the crystal phases.
The variational energy for a floating Wigner crystal is lower than that of the fixed crystal by the kinetic energy of the center of mass, which is a negligible energy per particle for large systems. The expectation value of any operator that depends only on the difference between electron coordinates is the same for the fixed and floating crystals. We can therefore obtain the PCFs of a floating Wigner crystal rather simply from calculations on the fixed crystal.