Ab Initio Path Integral Monte Carlo Simulations of the Uniform Electron Gas on Large Length Scales

The accurate description of non-ideal quantum many-body systems is of prime importance for a host of applications within physics, quantum chemistry, materials science, and related disciplines. At finite temperatures, the gold standard is given by ab initio path integral Monte Carlo (PIMC) simulations, which do not require any empirical input but exhibit an exponential increase in the required computation time for Fermionic systems with an increase in system size N. Very recently, computing Fermionic properties without this bottleneck based on PIMC simulations of fictitious identical particles has been suggested. In our work, we use this technique to perform very large (N ≤ 1000) PIMC simulations of the warm dense electron gas and demonstrate that it is capable of providing a highly accurate description of the investigated properties, i.e., the static structure factor, the static density response function, and the local field correction, over the entire range of length scales.

(N ≤ 1000) PIMC simulations of the warm dense electron gas and demonstrate that it is capable of providing a highly accurate description of investigated properties, i.e., the static structure factor, the static density response function, and local field correction, over the entire range of length scales.q/q F this work state-of-the-art

Table of Contents (TOC)/Abstract graphic.
The accurate description of fermionic quantum many-body systems is a paramount task within physics, quantum chemistry, material science, and related fields.An important subcategory is given by thermal simulations that describe quantum systems at finite temperatures.For example, warm dense matter (WDM), 1 an extreme state combining high temperatures (T ∼ 10 4 − 10 8 K), densities (n ∼ 10 22 − 10 27 cm −3 ) and pressures (P ∼ 1 − 10 4 GBar), is ubiquitous throughout our universe and occurs in a host of astrophysical objects such as giant planet interiors 2 and brown dwarfs. 3,4In addition, the fuel capsule in inertial confinement fusion applications has to traverse the WDM regime 5 in a controlled way to reach ignition. 6Consequently, WDM constitutes a highly active topic and is routinely realized in experiments at large research facilities using different techniques. 7[15][16] From a theoretical perspective, the rigorous description of such systems constitutes a difficult challenge as it must capture the complex interplay of effects such as non-ideality, quantum degeneracy and diffraction, as well as strong thermal excitations. 1,17In this context, the ab initio path integral Monte Carlo (PIMC) method 18 is a promising candidate as it is, in principle, capable to deliver an exact solution to the full quantum many-body problem of interest without the need for empirical input such as the exchange-correlation functional in density functional theory 19 or the self-energy in Green function approaches. 20,21fortunately, the PIMC simulation of fermions is afflicted with the notorious fermion sign problem, 22 which leads to an exponential increase in the required compute time for example upon increasing the system size N . 23,24ile a complete solution of the sign problem remains unlikely, the pressing need for an accurate description of quantum Fermi systems has ignited a number of promising developments in the field of quantum Monte Carlo simulations over the last decade, e.g.[27][28][29][30][31][32][33][34][35][36][37][38][39][40] On the one hand, these methods allow for a very accurate description of a given N -body system.In combination with appropriate finite-size corrections, 31,41,42 this has led to the first reliable parametrizations of the uniform electron gas (UEG) covering the entire WDM parameter space, 32,40,43,44 allowing for thermal DFT calculations of real WDM systems on the level of the local density approximation [45][46][47] and beyond. 48,49On the other hand, these tools by themselves are not capable to capture long-range phenomena 50 that manifest on length scales larger than the given box length L which, for a given density, is determined by the number of simulated particles N .This precludes, for example, the description of X-ray Thomson scattering (XRTS) experiments, a key diagnostic for WDM applications, [51][52][53][54] in a forward scattering geometry where the system is probed at a small momentum transfer q and, consequently, a long wavelength λ = 2π/q.Moreover, electrical 55 and thermal 56 conductivities are necessarily defined in the optical limit of q → 0, which makes the simulation of large systems even more important.ing to the physically relevant cases of Bose-, Boltzmann-, and Fermi-statistics, respectively.Subsequently, Dornheim et al. 58 have implemented the same idea into PIMC and found that it works remarkably well for weakly to moderately quantum degenerate systems, including the warm dense UEG.In a nutshell, these findings imply the intriguing possibility of PIMC simulations of large Fermi systems without the exponential bottleneck due to the sign problem.

Very recently, Xiong and Xiong
In the present work, we rigorously test this hypothesis by carrying out unprecedented large-scale PIMC simulations (with N ≤ 10 3 , see Fig. 1) of the UEG both at WDM parameters, and in the strongly coupled electron liquid regime. 59,60In this way, we unambiguously demonstrate the capability of the ξ-extrapolation method to accurately describe the system on all length scales, including the difficult long-wavelength limit of q → 0.
These findings open up a gamut of new possibilities to study WDM, ultracold atoms, and a host of other Fermi systems.
Path integral Monte Carlo and the ξ-extrapolation method.Let us consider a system of N = N ↑ + N ↓ unpolarized electrons in a cubic simulation cell of volume V = L 3 at an inverse temperature β = 1/k B T .Writing the partition function in coordinate representation then gives where the meta-variable R = (r 1 , . . ., r N ) T contains the coordinates of both spin-up and spindown electrons.Due to the indistinguishable nature of electrons of equal spin-orientation, we have to sum over all possible permutations of particle coordinates that are realized by the corresponding permutation operators πσ N i with i ∈ {↑, ↓}.A special role is played by the aforementioned continuous spin-variable ξ, which effectively controls the likelihood of pair exchanges (with N pp being the particular number of pair exchanges needed to realize a particular permutation) in the PIMC simulation. 58,61For ξ ≥ 0, all contributions to Z are positive, and no sign problem occurs.In contrast, positive and negative contributions increasingly cancel for ξ < 0, which is most severe in the fermionic limit of ξ = −1.To avoid the associated exponential decrease in the accuracy, Xiong and Xiong 57 have proposed to extrapolate to the latter by fitting a quadratic polynomial of the form to PIMC results for an observable Ô in the sign-problem free domain, i.e., ξ ≥ 0. Subsequently, Dornheim et al. 58 have found that this approach works remarkably well for a host of observables, including the energy E, static structure factor S(q), and even the imaginarytime density-density correlation function F (q, τ ). 54,62ditional details regarding the derivation of the imaginary-time PIMC method 18 as well as the fermion sign problem 23,24 have been presented in the literature and need not be repeated here.For completeness, we note that we use a canonical adaption of the worm algorithm by Boninsegni et al. 63,64 based on the extended ensemble idea presented in Ref. 65 All results have been obtained for P = 200 imaginary-time propagators with a primitive factorization scheme, and the convergence with P has been carefully checked.
Results.Let us start our investigation by considering the UEG at a Wigner-Seitz radius of r s = 2 and degeneracy temperature Θ = k B T /E F = 1 (where E F is the usual Fermi energy 68 ).This corresponds to a metallic density that can be realized for example in experiments with aluminium 69 at the electronic Fermi temperature of T = 12.53 eV.In the left panel of Fig. 2, we show our new results for the static structure factor S(q).More specifically, the dashed blue curve corresponds to the well-known random phase approximation (RPA), where the electronic density response to an external perturbation is treated on a mean-field level. 54,68For the case of the UEG, the RPA becomes exact in the long-wavelength limit of q → 0, but systematically underestimates the true SSF around the Fermi wavenumber q ∼ q F .These shortcomings have been corrected in the semi-empirical effective static approximation (ESA), 66,67 which is shown as the solid black line, and which is expected to give a highly accurate description of the UEG over the entire q-range at these conditions.This is indeed verified by the green crosses depicting exact PIMC results for S(q) that have been obtained based on simulations with N = 34 electrons without the ξ-extrapolation method.
Alas, these data are only available above a minimum wavenumber of q min = 0.63q F which is defined by the size of the simulation cell. 31 overcome this limitation, we have carried out extensive new calculations with N = 200 and N = 1000 electrons using the method explored in Refs. 57,58In Fig. 1, we show a snapshot from a corresponding PIMC calculation with N = 1000.A particular strength of the method is the full treatment of quantum diffraction effects over the entire length scale, as each electron is represented by an entire path of P = 200 positions along the imaginary-time domain; these paths would collapse to point particles in the classical limit of β → 0.
Coming back to the SSF shown in Fig. 2, we find that the new results for N = 200 (red circles) and N = 1000 (yellow circles) that have been obtained based on Eq. (2) using input data from the sign-problem free domain of ξ ≥ 0 (shaded grey area) are in excellent agreement with the ESA results and the exact PIMC results for N = 34, where they are available.In particular, we unambiguously demonstrate that this approach is capable to recover the q → 0 limit, which is known exactly in the case of the UEG; see also the inset showing a magnified segment around this regime.The somewhats larger spread of the data points for the N = 200 and N = 1000 cases is due to the computational cost increasing with ∼ N 2 and thus the significantly reduced number of snapshots available for averaging with manageable computational cost.
In the right panel of Fig. 2, we repeat this analysis for the static linear density response function χ(q), which can be computed from PIMC results for the ITCF F (q, τ ) based on the Overall, we find the same general trends as for S(q): the RPA underestimates the true magnitude of the density response around the Fermi wavenumber; the ESA is quasi-exact for all q and nicely agrees with the direct PIMC results; most importantly, our new, largescale PIMC simulations are capable to predict the correct q → 0 limit, see also the inset.
To further highlight the high quality of these results, we consider the exact expression where χ 0 (q) is the temperature-dependent Lindhard function describing the density response of an ideal Fermi gas 68 and the complete, wavenumber resolved, information about electronic exchange-correlation effects is included in the static local field correction (LFC) G(q). 71,72 Therefore, the LFC constitutes key input for a number of applications, such as the computa-tion of electrical conductivities, 73 the interpretation of XRTS experiments, 74 the development of thermal XC functionals, 75-78 and as the exchange-correlation kernel in time-dependent DFT simulations. 79,80In Fig. 3, we show the LFC for the same conditions as in Fig. 2.
The solid blue line corresponds to the recent analytical parametrization of the ESA 67 that includes the correct q → 0 limit given by the compressibility sum-rule (CSR), where f xc is the exchange-correlation energy of the UEG. 32,40,43,44In practice, we evaluate Eq. ( 5) using the parametrization by Groth et al., 32 and the results are included as the dash-dotted yellow curve.Furthermore, the dashed black line shows the PIMC based neural net representation of G(q; r s , Θ) from Ref. 71 The ESA and ML curves are very similar in the depicted q-range and only noticeably differ in the single-particle limit of q ≫ q F .As usual, the green crosses show exact PIMC results for N = 34 that have been obtained by inverting Eq. ( 4), and the red circles the corresponding ξ-extrapolated PIMC results for N = 200.
Again, we find that the ξ-extrapolation method very accurately describes this sophisticated exchange-correlation property over all length scales.This raises the intriguing possibility to study long-range correlations in other Fermi systems such as warm dense hydrogen, where the true q → 0 limit is a subject of active investigations. 50,55,81 a second example, we consider the UEG at r s = 10 and Θ = 0.5.3][84] In Fig. 4, we show the corresponding SSF and static linear response function with the usual color code.To assess the accuracy of the ξ-extrapolation at these conditions, we have carried out new, exact PIMC simulations for N = 14 electrons (green crosses).We find an average sign of S = 0.02, making them computationally involved, but still feasible. 23Considering the static structure factor, the

Χ(q)
q/q F -0.001 0 0 0.2 0.4 PIMC reference results are in very good agreement with the ESA curve; small systematic deviations are only visible in the vicinity of the peak where the latter slightly overestimates the true SSF.The red circles show the ξ-extrapolated results, which nicely agree with the green crosses everywhere.Interestingly, we find that the effect of quantum statistics on S(q) is small (see the shaded grey area indicating the sign-problem free domain of ξ ≥ 0), even though the degree of cancellations of positive and negative terms is substantial already for the relatively small system size of N = 14.A similar observation has been reported recently for the case of ultracold 3 He, 11 which, too, constitutes a strongly coupled quantum liquid.
Finally, the yellow circles show new, ξ-extrapolated results for N = 200 electrons, and we find the same good performance as for the case of r s = 2 investigated above.
In the right panel of Fig. 4, we show results for the static linear density response function χ(q).First, we observe that the impact of quantum statistics is considerably larger than Figure 5: ξ-dependence of the static linear density response function χ(q) for N = 14, r s = 10, and Θ = 0.5 at q ≈ 2q F .The blue diamonds depict the exact PIMC results and the various curves have been obtained by quadratic fits based on Eq. ( 2) using as input data in the interval of ξ ∈ [0, ξ max ].
for S(q) = F (q, 0).Indeed, Eq. ( 3) implies that χ(q) directly depends on the imaginarytime diffusion of the electrons, which is more sensitive to quantum effects compared to the static structure of the system.Second, we find that the ξ-extrapolation based on Eq. ( 2) becomes somewhat inaccurate when the full sign-problem free interval of ξ ∈ [0, 1] is used.This is investigated in more detail in Fig. 5, where we show the ξ-dependence of χ(q) at q ≈ 2q F .Specifically, the blue diamonds show exact PIMC results, which are available as a benchmark over the full ξ-range in this case, and the various curves correspond to quadratic extrapolations based on Eq. ( 2) using as input data within different intervals ξ ∈ [0, ξ max ].
The dashed black line corresponds to the usual choice of ξ max = 1 and underestimates the true density response in the fermionic limit of ξ = −1 by about 1%.We note that this is still a reasonable level of accuracy for the description of a strongly coupled Fermi liquid, but it is considerably worse compared to the WDM case investigated above.The dash-dotted yellow and dotted green curves correspond to ξ max = 0.7 and ξ max = 0.5, respectively, and exhibit a substantially improved agreement with the PIMC data for ξ < 0. Empirically, these values of ξ seem equally reasonable and we have used ξ = 0.5 to compute the red circles in the right panel of Fig. 4. Evidently, this truncated extrapolation scheme works very well over the entire q-range.Finally, the solid red curve in Fig. 5 corresponds to ξ max = 0.3.In this case, the ξ-extrapolation is based on only four data points, which leads to a larger uncertainty compared to the other analyzed fitting intervals.
Returning to the static linear density response function shown in Fig. 4, we find that the effect of quantum statistics is restricted to the range of 1.75q F ≲ q ≲ 4q F for these conditions.Unsurprisingly, the ξ-extrapolated results for N = 200 electrons (yellow circles) are in excellent agreement with the solid black ESA curve everywhere, and reproduce the correct q → 0 limit, see also the inset.We thus conclude that the ξ-extrapolation method is applicable beyond the WDM regime, and constitutes a valuable tool for the investigation of strongly coupled Fermi liquids.
Summary and outlook.In this work, we have presented extensive new ab initio PIMC simulations of the UEG with unprecedented system size (N ≤ 1000 electrons).In this way, we have unambiguously demonstrated that the ξ-extrapolation method 57,58 is capable to describe non-ideal Fermi systems over all length scales, including the long-range limit (q → 0) that is known exactly in the case of the UEG, but not for other systems such as warm dense two-component plasmas or ultracold atoms.From a physical perspective, we have considered both the WDM regime that is of high current interest e.g. for the description of laser fusion applications and astrophysical models, and the strongly coupled electron liquid regime, which exhibits phenomena such as the roton-type feature of the dynamic structure factor that is interesting in itself.We are convinced that these findings open up a host of new avenues for impactful future research in a great variety of research fields.
First, an intriguing application of the ξ-extrapolation method stems from its unique capability to accurately probe even the largest length scales of the UEG.Naturally, when approaching the long-wavelength limit, otherwise inaccessible information is unlocked that concerns specific thermodynamic quantities as well as transport coefficients.(a) The q → 0 limit of the static local field correction directly leads to the isothermal compressibility via the eponymous sum rule without the need for thermodynamic integration or differentia-tions. 68Such a thermodynamic route that is independent from the internal energy would provide a rigorous check for existing UEG equations of state 32,44 which are typically free energy representations whose derived thermodynamic properties are known to suffer from certain pathologies. 48(b) The frequency dependent electrical conductivity is connected to the long-wavelength limit of the current response function, while the thermal conductivity is connected to the hydrodynamic limit of the heat (or energy) current response function, see the respective Kubo and Green-Kubo formulas. 85,86(c) Within the current density version of linear response theory, it is known that the generalized visco-elastic coefficients are connected to long-wavelength limits that involve the respective longitudinal and transverse dynamic local field corrections. 87,88Moreover, the shear viscosity is connected to a hydrodynamic limit that involves the transverse local field correction. 68,87,88(d) Within the spin resolved density version of linear response theory, it is known that the static spin susceptibility is connected to the q → 0 limit of the static spin-antisymmetric local field correction. 89In addition, the spin diffusion coefficient/constant is connected to a hydrodynamic limit that involves the dynamic spin-antisymmetric local field correction. 851][92][93][94] This will facilitate the interpretation of XRTS experiments in a forward scattering geometry where the system is effectively probed on large length scales.6][97] In accordance with the above reasoning, a second potential application of the ξ-extrapolation method that is related to the study of WDM concerns the estimation of thermodynamic quantities and transport coefficients such as the thermal conductivity 56 and the electrical conductivity, 55 which may substantially depend on exchange-correlation effects. 50turning to the UEG itself, it might be interesting to study the spin-resolved density response, as well as the spin-resolved components of the static structure factor which, in contrast to the spin-averaged χ(q) and S(q), are not described properly by the RPA or by more sophisticated dielectric schemes [98][99][100][101][102][103][104] even in the long-range limit of q → 0. 105 Such a study might give new insights into the performance of different approximations to the LFC, and can provide useful input and benchmark data for other applications.
A further interesting line of research is given by the application of the ξ-extrapolation method to other systems.Based on our encouraging results for the strongly coupled electron liquid, we propose to utilize the method to ultracold atoms such as the short-range interacting 3 He, but also dipolar systems. 23,106For these systems, interesting long-range phenomena include the acoustic mode in the dynamic structure factor for small q, and the long-wavelength limit of the SSF that is given by the compressibility sum-rule. 107Furthermore, we note that the method has already been applied to electrons in quantum dots in previous studies. 57,58,108deed, it is well known that trapped Fermi systems exhibit interesting effects such as the formation of Wigner molecules 109,110 and a negative superfluid fraction. 16,106,111nally, we stress that, despite its impressive performance, the ξ-extrapolation method is a very new technique, and its further methodological improvement remains in its infancy compared to more established methods such as restricted PIMC. 25,112,113The observed improvement of the extrapolation due to the truncation of the ξ-interval in the sign-problem free domain for r s = 10 and Θ = 0.5 suggests that much can potentially be gained by developing better extrapolation schemes.In addition, the absence of the exponential bottleneck poses new challenges, which might be met by combining PIMC with other ideas such as the adaptive long-range potential scheme that has been suggested in the recent Ref., 114 or the quantum ring-polymer contraction method that has been introduced in the context of PIMD.

Figure 1 :
Figure 1: Snapshot from an ab initio PIMC simulation of the warm dense UEG at r s = 2 and Θ = 1 with N = 1000 unpolarized electrons.Each electron is represented by an entire path along the imaginary time with P = 200 time-steps (blue spheres); their extension corresponds to the thermal wavelength λ β = 2πℏ 2 β/m e and takes into account quantum effects such as diffraction.