Static Response of the Strongly Coupled Electron Gas from Diffusion Monte Carlo

Introduction

The response properties of the homogeneous electron gas (HEG) have been of considerable theoretical interest for over half a century [1, 2, 3, 4]. They describe how the gas rearranges in response to an external potential of a given wavelength [5]. For a static external potential, this information is contained in the static linear density response function $\chi(q)$, and the local field factor $G(q)$ isolates the part of it due to exchange and correlation [3]. Both enter directly into descriptions of dielectric screening [6], effective electron-electron interactions [7, 8], and the exchange-correlation (xc) kernel of time-dependent density functional theory [9]. Early studies of the static response were based on analytical and diagrammatic schemes that satisfied the compressibility sum rule and the large-$q$ limit [1, 3, 6, 10, 11, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] [12].

The first quantum Monte Carlo (QMC) calculations of the static response function $\chi(q)$ and the local field factor $G(q)$ of the gas at zero temperature were carried out in parallel by two groups in 1992. Sugiyama, Bowen, and Alder demonstrated the approach on the charged Bose gas [29], while Moroni, Ceperley, and Senatore applied it to the two-dimensional electron gas and liquid ${}^{4}\mathrm{He}$ [30]. The first application to the three-dimensional electron gas was made by Bowen et al. [31], using variational Monte Carlo, fixed-node diffusion Monte Carlo (DMC), and released-node Green’s function Monte Carlo to confirm that the fixed-node results are exact within statistical error bars. The definitive zero-temperature results were later provided by Moroni, Ceperley, and Senatore in their benchmark study [32], which yielded accurate DMC data for $\chi(q)$ and $G(q)$ across the metallic density range $r_s = 2$–10, alongside an analytical interpolation of $G(q)$ that satisfies the known asymptotic limits. Their data, calculated by adding an external perturbation to the Hamiltonian, showed that $G(q)$ varies smoothly between these limits with a crossover near $2k_F$. These calculations have served as the primary benchmark for describing exchange-correlation effects in materials for nearly three decades.

At finite temperature, there have been significant efforts, especially using the path integral Monte Carlo (PIMC) approach. Dornheim et al. have computed $\chi(q)$ and $G(q)$ for the strongly coupled, i.e., low-density, electron liquid at reduced temperatures $\Theta= T/T_F = 0.5$–4 for a wide range of $r_s$ values including $r_s = 20$, $r_s = 50$, $r_s = 100$, and $r_s = 200$ via imaginary-time correlation functions [33, 34, 35, 36]. Their findings show indications of structure in $G(q)$ at intermediate wave vectors; however, with thermal and correlation effects inextricably entangled at finite $T$, a clear origin of this feature could not be established.

Following the calculations of Moroni et al., several efforts to parametrize $G(q)$ have been made. Corradini et al. fitted an analytical form to the Moroni et al. data [37], which made $f_{\mathrm{xc}}(q)$ conveniently available for use in ab initio calculations. This parametrization has been widely adopted ever since [9, 35, 38, 39, 40, 41, 42]. More recently, Kukkonen and Chen provided new results for $G(q)$ at $r_{\mathrm{s}}=1$ and $r_{\mathrm{s}}=2$ using diagrammatic QMC [8]. These results (along with the Moroni et al. data) were then used by Kaplan and Kukkonen to construct an improved parametrization [43].

Crucially, despite these efforts on parametrizing $G(q)$, the underlying QMC data for the unpolarized gas do not extend beyond $r_{\mathrm{s}}=10$. Two distinct factors have obstructed the extension of the perturbative DMC approach of Moroni et al. to lower densities. The first is cost. Each point of a $\chi(q)$ curve requires DMC calculations at several perturbation amplitudes, each preceded by an independent optimization of the trial wave function, which converges more slowly as the density decreases; a single curve at a single density can easily amount to $\sim10^{3}$ separate QMC calculations. This cost is met by modern implementations and computing resources [44, 45], which were instrumental for the present work. The second factor is not relieved by computing resources: several systematic errors that are controlled at metallic densities remain uncharacterized at $r_{\mathrm{s}}\gtrsim10$. The usable range of perturbation strength to ensure linear response has not been established, which increases computational cost if one must repeat calculations at each density and several perturbation strengths to verify the validity of linear response. The effect of the nodal surface of the trial wave function on the extracted $\chi(q)$ has not been tested, and the finite-size effects in $\chi(q)$ have not been investigated at large $r_{\mathrm{s}}$.

In this work, we present DMC calculations of the static response function $\chi(q)$ and the local field factor $G(q)$ of the three-dimensional electron gas at zero temperature, covering the entire liquid phase up to $r_{\mathrm{s}}=100$. The physical implications of this new regime are discussed in a companion paper [46], where we show the emergence of a nonmonotonic structure in $G(q)$ at intermediate wave vectors and how it manifests as short-range order in the electron liquid, together with a new parametrization of $G(q)$ that can be used for many-body calculations. In this paper, we focus on the methodology that made those findings possible, addressing in turn each of the obstacles above. We determine the optimal variation parameters of the trial wave function and reduce them to a closed form valid across the whole density range, so that the optimization does not have to be repeated for every calculation; we establish how far the perturbation amplitude can be increased while remaining within the linear response regime as the density decreases; we examine whether more sophisticated forms of the trial wave function meaningfully

alter the extracted susceptibility, and we investigate the finite-size effects at low density.

The rest of the paper is organized as follows. Section II gives a complete description of the methodology, including how the response function is extracted, how the regime of linear response is determined, how the trial wave function is constructed and optimized, how finite-size and time-step errors are treated, and how $G(q)$ is parametrized. Section III presents the results for $\chi(q)$ and $G(q)$ for $r_{\mathrm{s}}=2$–100 and highlights the features that emerge at large $r_{\mathrm{s}}$. Subsection III B relates these features to the structure of the strongly coupled electron liquid, as seen in its pair correlation function and static structure factor, and to its evolution toward the Wigner crystal. We conclude in Sec. IV.

Method

Static response from a cosine perturbation

We extract the response properties by subjecting an otherwise homogeneous electron gas to a weak static cosine perturbation,

where $v_{\mathbf{q}}$ is the amplitude of the perturbation. The response function $\chi(q)$ can then be extracted by expanding the perturbed energy (per particle) in powers of $v_{\mathbf{q}}$. Within linear response, the ground-state energy is

where $E_0$ is the unperturbed ground-state energy per particle (i.e., $v_{\mathbf{q}}=0$), and $\chi(q)$ is the static density-density response function. Higher-order terms carry information about the higher-order response, which we do not explore in this work. By fitting the energy to this form, one can extract $\chi(q)$ from the coefficient of the quadratic term.

A crucial quantity that can be obtained directly from $\chi(q)$ is the local field factor $G(q)$, which encodes the exchange and correlation effects beyond the random-phase approximation (RPA). It is defined via the Dyson equation,

where $f_{\mathrm{xc}}(q)$ is the exchange-correlation kernel, $v_{\mathrm{c}}(q)=4\pi/q^2$ is the Coulomb interaction, and $\chi_0(q)$ is the static response function of the noninteracting system [3]. The noninteracting and RPA limits are reached by setting $G(q)=1$ and $G(q)=0$, respectively. Since $\chi(q)$ will be obtained from stochastic estimates of the energy $E(v_{\mathbf{q}})$, it carries a statistical error, which propagates into $G(q)$ with a prefactor proportional to $(q/\chi(q))^2$ that inevitably diverges in both the low- and high-$q$ limits. However, this is of no practical concern since $G(q)$ is known analytically in both limits [25]. The low-$q$ limit of the local field factor is given by $\lim_{q\to0}G(q)=A(q/k_F)^2$ whereas the high-$q$ limit is given by $\lim_{q\to\infty}G(q)=C(q/k_F)^2+B$. All quantities $A$, $B$, $C$ are known expressions and parametrizations as a function of $r_s$ [25, 26, 32, 37]. We provide the explicit expressions for $A(r_s)$, $B(r_s)$, and $C(r_s)$ in Appendix A.

To compute the perturbed energies $E(v_q)$, we employ the state-of-the-art fixed-node DMC method [32, 47, 48], implemented in QMCPACK [44, 45]. DMC is the gold standard for computing ground-state properties of quantum many-body systems, and it has been extensively used as a benchmark for the homogeneous electron gas [48, 49]. Modifications were made to include the cosine perturba-

tion to the Hamiltonian. Once $E(v_q)$ is obtained, $\chi(q)$ is extracted from the energy by performing a least-squares fitting procedure to the form of Eq. (2). The procedure is summarized in Figure 1 for two different densities, $r_s=5$ and $r_s=70$, and for three different wave vectors. We carefully select the range of $v_q$ used to fit $E(v_q)$ and extract $\chi(q)$ so that it remains in the linear-response regime, where the signal-to-noise ratio is high. The resulting maximum perturbation is parametrized by the fit $v_q^{\max}=0.1212r_s^{-1.43}$ Ha, where the exponent is determined by the noninteracting problem and the prefactor is fitted to QMC data. More details of this selection process can be found in Appendix B.

Trial wave function

We construct our trial wave function using both Slater–Jastrow (SJ) and Slater–Jastrow-backflow (SJ-BF) forms [50, 51, 52, 53, 54]. The former is the standard form used in many QMC calculations in which the nodal surface is fixed by the Slater determinant of single-particle orbitals, while the latter allows for modifications of the single-particle coordinates in the Slater determinant by introducing a displacement for the $j$th electron, which depends on the positions of all other electrons. Such a transformation introduces additional variational parameters that can be optimized to further minimize the energy. To compare both forms and verify that our choice of trial wave function does not bias the extraction of $\chi(q)$, we performed a subset of calculations with and without a backflow transformation. We find that while the backflow transformation significantly improves the energies obtained from DMC, it does not affect the extraction of $\chi(q)$, which is related to the curvature of the energy with respect to the perturbation amplitude $v_q$, and not the absolute value of the energy itself. Figure 1 demonstrates this point for (a) $r_s=5$ with $N_e=162$, and for (b) $r_s=70$ with $N_e=38$, at three wave vectors. The resulting $\chi(q)$ values obtained with and without backflow transformation are in agreement within error bars, even though the energies themselves are significantly different. Thus, throughout this work, we will mainly use the Slater–Jastrow trial wave function without backflow transformation (unless stated otherwise) to extract $\chi(q)$, which is computationally more efficient and shown not to bias our results. The SJ trial wave function is constructed as $\Psi_{\mathrm{T}}(\mathbf{R})=e^{J(\mathbf{R})}\det[\phi_i(\mathbf{r}_j)]$, where $J(\mathbf{R})$ is the Jastrow factor that captures electron-electron correlations, and $\det[\phi_i(\mathbf{r}_j)]$ is the Slater determinant of noninteracting single-particle wave functions. For the latter, we will use the one of noninteracting single-electron [55] perturbed by an effective external potential of the same form,

By construction, $v_{\mathrm{ext}}^{\mathrm{eff}}$ does not have to be equal to the amplitude of the external potential $v_q$, so $\alpha$ must be treated as a variational parameter that yields the minimum energy. Therefore, the optimization is to be achieved in two steps. First, one performs variational Monte Carlo (VMC) calculations to minimize the energy computed from the Hamiltonian perturbed by Eq. (1), with respect to the parameters in the Jastrow factor, for a given $\alpha$ in Eq. (4). This procedure is then repeated for different $\alpha$ values to obtain the VMC energy as a function of $\alpha$, i.e., $E_v(\alpha)$, and extract the optimal $\alpha$. In this work, we find that the $\alpha$ that minimizes the variational energy is different for each value of $q$ and $r_s$ and needs to be optimized for each case separately. We emphasize that the $\alpha$-optimization is crucial to obtain reliable results for $\chi(q)$ for all $r_s$ values. This is demonstrated in Figure 2(a) for $r_s=50$ and $q=1.2k_F$. A poorly chosen $\alpha$ changes not only the energy but also the curvature of $E(v_q)$, whose sign can even be reversed, and gives rise to an unphysical $\chi(q)>0$. As $\alpha$ is reduced (for this particular wave vector and density), the curvature of $E(v_q)$ converges to a stable value, and yields the correct $\chi(q)$ at $\alpha\simeq0.15$. Figure 2(b) shows that the variance of the local energy decreases accordingly as $\alpha$ approaches the optimal value, i.e., the minimum of $E_v(\alpha)$ (vertical dashed line). Performing over 700 VMC calculations using the above approach, we found that the optimal $\alpha$ is reproduced by the following form for $r_s=2$–$100$:

where $\alpha_0(r_s)=0.7276r_s^{-0.6050}$ and $\alpha_1(r_s)=1.1603r_s^{0.3513}$. Figure 2(c) shows $\alpha(q,r_s)$ as a function of $q/k_F$ for six different values of $r_s$. Once the optimal $\alpha$ is obtained, we perform DMC calculations to obtain the perturbed energies $E(v_q)$, which are then used to extract $\chi(q)$.

Size effects and time-step error

It has been shown that the local field factor $G(q)$ is much less sensitive to finite-size effects than the energy [32] for $r_s\leq10$. This is due to the fact that exchange and correlation effects are short-ranged. Therefore, one can assume $G(\mathbf{q},N)\simeq G(\mathbf{q},\infty)$ to mitigate finite-size effects in our calculations as follows

where $\Delta\chi_0^{-1}(\mathbf{q})=\chi_0^{-1}(\mathbf{q},\infty)-\chi_0^{-1}(\mathbf{q},N)$ represents the finite-size correction in the noninteracting response function. Throughout this work, we use Eq. (6) to correct the finite-size effects in $\chi(q)$. At low-$r_s$ values, Eq. (6) is known to hold quite well [32]. However, at high-$r_s$ values, such a claim is questionable and merits further investigation.

In this subsection, we investigate the validity of this hypothesis at large $r_s$ by testing two finite-size features of the Fermi sphere: its resolution, set by $N_e$, and its isotropy at fixed $N_e$. Figure 3(a) shows the local field factor $G(q)$ at $r_s=50$ for different system sizes: $N_e=38$ (squares), $N_e=54$ (triangles), and $N_e=162$ (circles). The inset zooms into the region around $2k_F$. Overall, the DMC data points for different system sizes collapse onto a single curve, confirming the insensitivity of $G(q)$ to finite-size effects, even at large $r_s$ values. Although the statistical uncertainty of $G(q)$ at $q>4k_F$ and $q\lesssim1.6k_F$ is too

large to resolve differences between system sizes, both the small- and large-$q$ limits are analytically known (given in Appendix A). Figure 3(b) shows the difference between two $\chi(q)$-values that have the same $|\mathbf{q}|$ but different directions, i.e., $\delta\chi_{\hat{\mathbf{q}}}=\chi\left(\mathbf{q}=(2,2,1)\frac{2\pi}{L}\right)-\chi\left(\mathbf{q}=(3,0,0)\frac{2\pi}{L}\right)$ as a function of $r_s$. Since in the thermodynamic limit, $\delta\chi_{\hat{\mathbf{q}}}$ should vanish, any nonzero value of $\delta\chi_{\hat{\mathbf{q}}}$ indicates the presence of finite-size effects. This is shown by the triangles in Figure 3(b), which correspond to the raw $\chi(q)$ values. If Eq. (6) captures the dominant finite-size error, the corrected $\chi(q)$ should give $\delta\chi_{\hat{\mathbf{q}}}=0$. This is shown by the circles in Figure 3(b), where the difference is effectively reduced to zero within error bars. The anisotropy of $\chi(q,N)$ at fixed $N_e$ is therefore entirely due to $\chi_0(q,N)$, and it is removed by the correction in Eq. (6).

The remaining question is how large this correction is, and how it changes with increasing coupling. For this, we show in Figure 4(a) the finite-size correction to the static response, $\Delta\chi\equiv\chi(q,\infty)-\chi(q,N)$, obtained from Eq. (6) and given in units of $N_F=k_F/\pi^2$. For all densities, $\Delta\chi$ is localized in $q$ and vanishes at small and large $q$. The correction is largest at high density around $q\simeq1.5k_F$, and it becomes smaller at low density, where it shifts to $q\simeq2k_F$. It is predominantly negative, i.e., the finite system underestimates $|\chi(q)|$ where the correction is largest, and its magnitude always remains below $0.1N_F$ (where $N_F=-\chi_0(q\to0)$ sets the natural scale of the response). Its root-mean-square [Figure 4(b)] over all computed wave vectors decreases from $\approx4\times10^{-2}$ at $r_s=2$ to $\approx1\times10^{-2}$ at $r_s\gtrsim70$ for $N_e=162$, and approaches the typical statistical uncertainty of $\chi(q)$ (around $0.5\times10^{-2}$ at $r_s=90$). $N_e=54$ shows the same decreasing trend. The finite-size correction to $\chi(q)$ from

Eq. (6) is therefore modest at all densities, and it decreases further as the coupling increases. Twist averaging is an alternative route to removing one-body finite-size errors, but is more expensive, since it requires an independent calculation for each twist. In Appendix C, we show that at fixed $N_e$ (canonical twist averaging), it approaches our corrected results to within 2–3%, with a small residual bias due to keeping the number of electrons fixed.

Besides finite-size effects, the DMC energies carry a time-step error. Since $\chi(q)$ is extracted from small energy differences between perturbed states, any time-step bias that depends on $q$ or $v_q$ would directly bias the curvature of $E(v_q)$. Figure 5 therefore shows the time-step extrapolation for the perturbed gas at two wave vectors and two densities. The DMC energies exhibit a linear dependence on the time step $\tau$ as it approaches zero. For $r_s=5$, the energy is converged within error bars for $\tau\leq0.5\ \mathrm{Ha}^{-1}$, while for $r_s=70$, the energy is converged for $\tau\leq14\ \mathrm{Ha}^{-1}$. This is consistent with the fact that for such a large $r_s$ value, the system is in a strongly correlated regime where the energy scales shrink, and the DMC algorithm can take larger time steps without losing accuracy. In our production runs, we have used $\tau=r_s/40\ \mathrm{Ha}^{-1}$ (vertical dashed lines), which is well within the converged regime for all $r_s$ values considered in this work.

Analytical fit

All solid curves labeled “This work” that we will present in Sec. III are obtained from a single parametrization of $G(q)$, introduced in the companion paper [46]. It interpolates between a low- and a high-$r_s$ form,

where both terms take the form of Kaplan and Kukkonen [43],

with $x=q/k_{\mathrm{F}}$, the exact coefficients $A$, $B$, and $C$ given in Appendix A, and the crossover function

The function $h$ decays from 1 to 0, with position set by $a_4$ and width by $a_3$. Since $h \to1$ for $x \to0$ and the intermediate term vanishes faster than $x^2$ for $n>0$, Eq. (8) satisfies the compressibility sum rule and the large-$q$ limit exactly, at every $r_s$. The low-$r_s$ form $G_{\mathrm{lo}}(q)$ is the original parametrization of Ref. [43], with $n=4$ and $\theta=a_0+a_1e^{-a_2r_s}$. Kaplan and Kukkonen reported $a_0$ and $a_2$ with sizable uncertainties, due to the growing error bars of the diagrammatic QMC data above $2k_{\mathrm{F}}$ [see orange shaded area in Figure 6(d)]. Refitting these two parameters to include our DMC data, we obtain $a_0=-0.0046277$, within the original error bar, and $a_2=0.622624$, at the edge of the original uncertainty, and keep the remaining parameters at their published values, $a_1=0.0155766$, $a_3=3.516054$, and $a_4=1.015830$.

The high-$r_s$ form $G_{\mathrm{hi}}(q)$ keeps the same structure, but lets the exponent evolve with density, $n=a_5+a_6e^{-a_7r_s}$, while $\theta=a'_0+a'_1e^{-a'_2r_s}$ keeps the original form with refitted parameters. The crossover $h$ is moved to larger $x$ through $a'_4$, with a broader width $a'_3$. At large $r_s$, $\theta$ becomes negative, and the intermediate term allows for a possible nonmonotonic behavior of $G(q)$, which we will discuss in Sec. III. The values of the primed parameters are given in Ref. [46]. The switching function $s$ is negligible at $r_s=1$, reaches $s\simeq0.06$ at $r_s=2$, and is close to unity for $r_s\geq5$, so that $G_{\mathrm{lo}}(q)$ is recovered in the regime of metallic densities, and $G_{\mathrm{hi}}(q)$ describes the strongly coupled regime.

Results

We have performed DMC calculations for the unpolarized three-dimensional electron gas between $2\leq r_s\leq100$, with $N_e=38$, 54, and 162 electrons. The values of $r_s$ are chosen to cover the entire liquid phase of the electron gas, and perhaps even beyond the Wigner crystallization transition, since current estimates for the transition vary between $r_s\simeq87$ and $r_s\simeq113$ [48, 56, 57, 58]. Note, however, that our ansatz is limited to represent the liquid phase only, even if it is metastable. Our results are compared to the existing Corradini et al. parametrization, which becomes an extrapolation at $r_s>10$, in order to highlight new features emerging in the strongly coupled regime.

Static response function and local field factor from DMC

We present our results for the static density response function $\chi(q)$ and the local field factor $G(q)$, using Eq. (3), for $r_s=2$, 5, 10 in the upper and lower panels of Figure 6, respectively. The DMC data are shown as circles, while the dash-dotted curve is the analytical parametrization of Ref. [37] designed specifically for these $r_s$ values to fit the data from Ref. [32]. As a first demonstration, our results for $\chi(q)$ are in excellent agreement with the existing parametrizations fitted at these same densities. The response of the system vanishes as $q\to0$ quadratically (governed by the compressibility limit) due to the long-range nature of the Coulomb interaction, and peaks at $q\simeq2k_{\mathrm{F}}$ before approaching the RPA limit at large $q$. The peak becomes more pronounced with decreasing density, while its departure from RPA is an indication of growing exchange-correlation effects. Importantly, our simulations provide the first DMC data for $q<k_{\mathrm{F}}$, and they are consistent with the compressibility limit. This is important for numerous applications [8, 13, 59, 60, 61].

On the other hand, the case of $G(q)$ is more delicate, as it isolates the exchange-correlation contribution that may not be immediately apparent in $\chi(q)$. We show $G(q)$ in Figure 6(d)–(Figure 6f), where the small- and large-$q$ limits are given as gray dashed lines. For $r_s=2$, we present additional data (orange shaded region) obtained from diagrammatic QMC [8] which show a hump in $G(q)$ near $2k_{\mathrm{F}}$ that is not present in previous parametrizations. This feature is included in the parametrization of Kaplan and Kukkonen [43], i.e., the intermediate term in Eq. (8). However, the amplitude of the hump is uncertain due to increasingly large error bars in the QMC data for $q>2k_{\mathrm{F}}$. Our DMC data confirm the presence of such a feature, but we find the enhancement to be quite small. Moreover, our significantly smaller error bars (obtained using a backflow transformation for $q>1.75k_{\mathrm{F}}$ and $r_s=2$ to reduce the variance) resolve the discrepancy between DMC and diagrammatic QMC reported in Ref. [8], as both methods now agree within error bars. Our DMC results, combined with the diagrammatic QMC data, therefore allow us to determine the parameters of Kaplan and Kukkonen more precisely (Sec. II D).

With the smaller error bars, the large-wave-vector regime becomes reliably accessible. Starting from $r_s=5$, we find that $G(q)$ exhibits a softening above $2k_{\mathrm{F}}$ and deviates from the monotonic behavior fitted by the existing parametrizations, before approaching the high-$q$ limit. This is because the intermediate-$q$ regime was essentially unexplored by previous DMC data (e.g., at the lowest density, $r_s=10$, the data of Moroni et al., shown as red squares in Figure 6(f), do not extend beyond $q\simeq3k_{\mathrm{F}}$, but the feature is discernible in their data as well), so the parametrizations of Moroni et al. and Corradini et al. were designed as smooth crossovers between the known low- and high-$q$ limits [62]. Our parametrization, designed to capture the nonmonotonic behavior of $G(q)$ in the intermediate-$q$ region, accurately reproduces the DMC data throughout the whole range.

This feature is a precursor of a more pronounced nonmonotonicity due to growing exchange-correlation effects, as we show in Figure 7 at $r_s=20,30,50$. At stronger coupling, we see that $G(q)$ develops a clear structure that gives rise to a local minimum. In contrast to the existing parametrization proposed by Corradini et al. (dash-dotted line, which is an extrapolation in this density range) that assumed monotonic behavior, the amplitude of this structure grows systematically with $r_s$ and $G(q)$ approaches 1 at its local minimum. The position of the minimum shifts only slightly with density, moving from $q \simeq3.3k_F$ to $3.2k_F$ as $r_s$ increases. The suppression toward $G(q) \simeq1$ implies that the response function $\chi(q)$ tends toward the noninteracting response around the same $q$-value. This is an interesting feature, since it implies that the effective interaction $v_c(q)[1-G(q)]$ between electrons nearly vanishes (see Eq. (3)). The exchange-correlation effects cancel the Coulomb interaction, and the gas responds to a perturbation of this wave vector like a noninteracting system despite being strongly correlated. Our parametrization (Sec. II D) is constructed to capture this nonmonotonic behavior. Similar observations in the Bose liquid, the one-component plasma, and the 2D electron gas have been reported in the literature, with various possibilities attributed to strong liquid [63] or Coulomb interaction [64], scattering of electron pairs of opposite spin [65], the $2k_F$-singularity in the response function responsible for Friedel oscillations [37], etc. In our companion paper [46], as well as in Sec. III B, we show how this feature in $G(q)$ relates to short-range order in the electron liquid as it approaches the Wigner crystallization transition.

Beyond the local minimum, we see in Figure 7 that $G(q)$ approaches its large-$q$ limit. Whether it does so with further, smaller oscillations cannot be resolved with our current data. Similar nonmonotonic structure in $G(q)$ has been reported in PIMC calculations [33, 34], but whether it is a ground-state property or a finite-temperature effect has remained an open question, since the fermion sign problem prevents PIMC from reaching low temperatures [35, 66]. Therefore, for the high-$r_s$ panels of Figure 7, we also show the finite-temperature PIMC results of Ref. [34] for $\Theta=0.5$, which allow a direct comparison between ground-state and finite-temperature results. As panels (d)–(f) show, the PIMC data are in quantitative agreement with our $T=0$ results, demonstrating that this structure was inherited from the ground state of the strongly coupled electron liquid rather than being a thermal effect.[67] Moreover, via Eq. (3), the agreement in $G(q)$ implies that the thermal changes in $\chi^{-1}(q)$ and $\chi_0^{-1}(q)$ largely cancel, i.e., that the temperature dependence of $\chi(q)$ at $\Theta=0.5$ is carried predominantly by $\chi_0(q)$. This is consistent with panels (a)–(c) of Figure 7, where $\chi(q)$ itself, in which no such cancellation occurs, agrees less well between PIMC and DMC. This suggests that $G(q)$ determined at moderate temperatures can be used to obtain $\chi(q)$ at $T\to0$, and vice versa.

Finally, we show in Figure 8 the results for $r_s=100$, the lowest density we consider. The inset shows the ratio $\chi(q)/\chi_0(q)$. At such low density, whether the ground state of the electron gas is liquid or crystalline is still not precisely known, and our data show interesting features that could be related to the onset of such a transition. First, we see that the suppression of $\chi(q)$ toward the noninteracting response (dashed gray line) is the strongest among all densities considered, with $\chi(q_{\min})/\chi_0(q_{\min}) = 1.00 \pm0.02$ at $q \simeq3.2k_F$ (see inset). Second, we also see a secondary peak in $\chi/\chi_0$ at about twice the wave vector of the first peak, which leaves a direct imprint on the excitation spectrum of the electron liquid [46].

Relation to short-range order

The structure in $G(q)$ emerges in the same range of $r_s$ in which $g(r)$ and $S(q)$ develop pronounced correlation features, suggesting that the two are related. Here, we show that the feature in $G(q)$ is a manifestation of the same physical phenomenon that shapes $g(r)$ and $S(q)$, commonly referred to as short-range order.

Electron liquid with short-range order

Before turning to our analysis, we characterize how the electron liquid evolves as the coupling increases until it eventually crystallizes into a bcc Wigner crystal at $r_s \simeq87$–$113$ [48, 56, 57, 58]. Figure 9 shows the pair correlation function $g(r)$ from weak coupling up to near crystallization ($r_s = 90$), obtained from separate DMC simulations of the unperturbed electron gas with a liquid trial wave function. For comparison, we also show $g(r)$ of the bcc Wigner crystal at $r_s = 110$ (solid black line), taken from the DMC calculations of Ref. [57]. As the coupling increases, a peak develops in $g(r)$ and becomes pronounced for $r_s \gtrsim20$. The location of this peak is $r \simeq3.3/k_F$, which corresponds to the typical separation between an electron and its nearest neighbors. Interestingly, this peak coincides with the nearest-neighbor distance in the bcc Wigner crystal, $r_{\rm nn} = 3.38/k_F$ [68]. While the oscillations of $g(r)$ decay in the liquid, they persist to long range in the crystal phase.

The emergence of such a peak at the nearest-neighbor distance has been identified as short-range order in the literature [69, 70], and indicates that the liquid is approaching the Wigner crystallization transition [71, 72]. Apart from the peak and oscillations of $g(r)$, short-range order also shows up in the steepness of $g(r)$ on its rising edge, around $k_F r = 2$, close to the radius of the sphere available to each electron, i.e., the Wigner–Seitz radius $r_s \simeq1.92/k_F$. Compared with the weakly coupled liquid at $r_s = 5$, $g(r)$ drops much more steeply toward $r \to0$ when short-range order is present. This depletion is more easily seen in the xc hole, $n_0[g(r)-1]$, which measures the depletion of electrons around a reference electron [3, 73]. The inset of Figure 9 shows the xc hole as a two-dimensional cut through the reference electron for $r_s=5$ and $r_s=90$. The liquid with short-range order has a noticeably more localized xc hole than the weakly coupled liquid. This can be quantified by the fraction of the missing electron within the Wigner–Seitz radius, $-n_0\int_0^{r_s}[g(r)-1]4\pi r^2\,dr$, which would equal unity for a completely empty Wigner–Seitz sphere. This fraction grows from 57% at $r_s=5$ to 95% at $r_s=90$ (98% in the crystal). The localized xc hole and the nearest-neighbor peak thus appear together, suggesting that both are signatures of short-range order: the region around each electron is nearly depleted up to its Wigner–Seitz radius, while its neighbors accumulate at the nearest-neighbor distance, reminiscent of the local structure of the bcc Wigner crystal.

In reciprocal space, short-range order mainly corresponds to a peak in the static structure factor $S(q)$, arising from the oscillations in $g(r)$, followed by further oscillations, arising from the steepness of the xc hole, that are damped out at large wave vectors. Just as for $g(r)$, this peak in $S(q)$ coincides with the first Bragg peak of the incipient Wigner crystal, located at $q\simeq2.28k_F$. The location of the peak therefore tracks the evolution from short-range order in the liquid to long-range order in the crystal. In the following, we relate these signatures of short-range order, the nearest-neighbor peak in $g(r)$ and $S(q)$ and the steepness of the xc hole, to the local field factor $G(q)$.

Calculations of $S(q)$ and $g(r)$ from $G(q)$

The static structure factor $S(q)$ can be obtained directly from DMC. Alternatively, it follows from linear response theory via $S(q)=-\frac{1}{\pi n_0}\int_0^\infty\operatorname{Im}\chi(q,\omega)\,d\omega$, which, however, requires the dynamic density response function $\chi(q,\omega)$. The latter is determined by the dynamic local field factor $G(q,\omega)$, the frequency-dependent generalization of $G(q)$. Here, we use the adiabatic (static) approximation $G(q,\omega)\approx G(q)$, which is exact at $\omega=0$ and neglects the frequency dependence. Our parametrization of $G(q)$ then yields $\chi(q,\omega)$ via the Dyson equation,

where $\chi_0(q,\omega)$ is the noninteracting density response function [3]. From here, one can obtain the dynamic structure factor as $S(q,\omega)=-\frac{1}{\pi n_0}\operatorname{Im}\chi(q,\omega)$.

Figures 10 and 11 show the static structure factor $S(q)$ and the pair correlation function $g(r)$ obtained from DMC calculations (circles, calculated independently from the unperturbed gas) and from the adiabatic approximation (lines) using Eq. (10) for $r_s=20,30,50,80,90$. Since our goal is to investigate the effect of the feature in $G(q)$ on $S(q)$, we also show the case where $G(q)$ is featureless for comparison, i.e., the parametrization of Corradini et al. (dash-dotted lines). In all cases, $g(r)$ is obtained from the corresponding $S(q)$ by Fourier transform. Since the adiabatic approximation neglects the frequency dependence of $G(q,\omega)$, the resulting large-$q$ behavior of $S(q)$ is incorrect, and the corresponding small-$r$ behavior of $g(r)$ is unreliable (gray curves in Figure 11). In the DMC data, a peak develops in $S(q)$ near $q\simeq2.2k_F$ as the liquid becomes more strongly coupled. As correlations grow, the peak becomes more pronounced and is followed by a dip, consistent with the short-range order discussed above.

The solid curves in Figure 10 show that the effect of the feature in the local field factor $G(q)$ is twofold. First, it reproduces the dip in $S(q)$ after the main peak as $r_s$ increases, which a featureless $G(q)$ completely misses (red arrows). The equivalent picture in real space is shown in Figure 11, where the red arrows mark the corresponding suppression of $g(r)$ on its rising edge, at $k_F r \simeq2.24$, due to the feature in $G(q)$. Since this suppression is tied to the localization of the xc hole (Figure 9), the feature in $G(q)$ is a signature of short-range order in the electron liquid. Consistently, the real-space xc kernel $f_{\rm xc}(r)$, i.e., the Fourier transform of $-v_{\rm c}(q)G(q)$, has a sharp positive peak at $k_F r \simeq2.2$, on the rising edge of the xc hole (dashed circles in the inset of Figure 9), at the same distance as the red arrows in Figure 11. This peak implies that density modulations cost exchange-correlation energy at this distance. Both this peak and the suppression marked by the red arrows disappear when the feature is removed from $G(q)$. For details, see the companion paper [46].

Second, the feature also reproduces the correct location of the main peak in $S(q)$. Without it, the main peak position is overestimated at all densities, which in real space (Figure 11) corresponds to an underestimation of the typical separation between an electron and its nearest neighbors (black arrows in Figure 10 and Figure 11). These observations indicate that the features in $S(q)$, $g(r)$, and $G(q)$ are related and emerge together with short-range order. This is consistent with studies of the strongly coupled one-component plasma, the classical analog of our system, where the relation between $G(q)$ and $S(q)$ is well established [74]. Although the adiabatic approximation is quantitatively inaccurate at large $q$ and small $r$, as discussed above, the qualitative agreement with DMC is clear. The same agreement also holds for other approximations, such as the single-mode approximation (SMA), in which the dynamic structure factor is approximated by a single sharp collective mode, $S(q,\omega) \approx S(q)\delta(\omega-\omega_q)$ [46].

Conclusions and Outlook

We have presented diffusion Monte Carlo calculations of the static density response function $\chi(q)$ and the local field factor $G(q)$ of the three-dimensional electron gas at zero temperature. Our results extend the zero-temperature QMC data beyond the metallic regime and across the liquid phase up to the vicinity of Wigner crystallization. These calculations were made possible by a careful analysis and control of the systematic errors that arise in the strongly coupled regime. First, we established a semi-empirical scaling of the linear-response window with density, which guides the choice of $v_q$. Second, we demonstrated that a backflow transformation, while lowering the energy substantially, leaves the extracted $\chi(q)$ unchanged, so that Slater–Jastrow trial wave functions suffice for the response. Third, we showed that the scaling $\alpha$ of the perturbing potential entering the Slater determinant must be optimized, since a poor choice alters the curvature of $E(v_q)$ and can even reverse its sign and thereby lead to an unphysical response. We provided a closed-form expression for the optimal $\alpha(q,r_s)$, valid for $2 \le r_s \le100$, which removes the need for a separate optimization in future calculations.

Finally, we showed that the local field factor remains insensitive to system size even at very large $r_s$ values wherever our data resolve it, and that the correction based on the noninteracting response alone removes the dominant finite-size error in $\chi(q)$. The time step was chosen small enough to ensure that the time-step bias is negligible. Together, these results show that the strongly coupled regime can be treated with standard Slater–Jastrow wave functions and modest system sizes, paving the way for future calculations of the response of the electron gas and related systems.

At metallic densities, our results agree with previous DMC data and parametrizations, extend them to $q < k_F$, where they are consistent with the compressibility sum rule, and resolve the discrepancy between DMC and diagrammatic QMC at intermediate wave vectors. Beyond this regime, $G(q)$ develops a pronounced nonmonotonic structure at intermediate wave vectors, already discernible at $r_s = 5$, whose amplitude grows with $r_s$ and which is absent from all existing parametrizations. The quantitative agreement with finite-temperature PIMC data at $r_s = 20$–$50$ establishes this feature as a ground-state property of the electron liquid. Finally, we have shown that the inclusion of this feature in $G(q)$ reproduces the correct physical behaviors of the pair correlation function $g(r)$ and the static structure factor $S(q)$, which are missed by a featureless $G(q)$, and that this feature is related to the emergence of short-range order in the electron liquid, as seen in the calculations of $g(r)$ and $S(q)$.

The parametrization of $G(q)$ used throughout this work, given in the companion paper [46], reproduces our data over the full liquid phase and is directly usable for numerous applications. For instance, it can be used to construct exchange-correlation functionals in time-dependent density functional theory [75, 76, 77], which is widely used for electronic structure calculations of materials. It can also be used to describe dynamical and spectral properties in the electron gas [9, 78] and in real materials [9, 40, 79, 80, 81, 82], to construct effective electron–electron interactions [7, 83, 84], and to extend response calculations of warm dense matter to lower temperatures [35], where current PIMC calculations are limited by the fermion sign problem. Additionally, our results can serve as benchmarks for future theoretical and computational studies of the electron gas and related systems.

The methodology developed here is not specific to the density response of the unpolarized three-dimensional gas. In fact, it can be applied directly to the spin response, through a spin-dependent perturbation, to partially polarized systems, and to the two-dimensional electron gas, where crystallization occurs at much higher density and where the interplay between charge and spin correlations at intermediate wave vectors remains to be established.

Data Availability

The data that support the findings of this study, i.e., the local field factor $G(q)$ from which the static density response function $\chi(q)$ can be obtained, the static structure factor $S(q)$, and the pair correlation function $g(r)$, are openly available in the Zenodo repository [85].

Acknowledgments

The authors would like to thank Saverio Moroni, Carl Kukkonen, and Carlos Rodriguez Perez for helpful discussions. The Flatiron Institute is a division of the Simons Foundation. This work has received state funding managed by the French National Research Agency (ANR) under the France 2030 program, reference “Grant No. ANR-22-2EXES-0013”. YY acknowledges support by the National Science Foundation (NSF) under grant number DMR-2532734. This work used Delta through allocation PHY250273 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. This research used the Delta advanced computing and data resource which is supported by the National Science Foundation (award OAC 2005572) and the State of Illinois. Delta is a joint effort of the University of Illinois Urbana-Champaign and its National Center for Supercomputing Applications.

Low- and high-$q$ limits of $G(q)$

We use Hartree atomic units, $k_F=(9\pi/4)^{1/3}/r_s$ and $n_0=3/(4\pi r_s^3)$. The local field factor satisfies

The coefficient $A$ follows from the compressibility sum rule [32, 37],

where $\mu_c=d(n\epsilon_c)/dn$, primes denote derivatives with respect to $r_s$, and the exchange contribution is the constant $1/4$.

The coefficient $C$ is set by the correlation contribution to the kinetic energy [25, 32],

where $t_c=-\epsilon_c-r_s\epsilon_c'$ is the correlation kinetic energy per electron given by the virial theorem. For the correlation energy per electron $\epsilon_c(r_s)$ we use the Perdew–Wang parametrization [86] of the unpolarized gas,

where $D=2a_0^{\mathrm{PW}}\left(\beta_1r_s^{1/2}+\beta_2r_s+\beta_3r_s^{3/2}+\beta_4r_s^2\right)$, with $a_0^{\mathrm{PW}}=0.031091$, $\alpha_1^{\mathrm{PW}}=0.21370$, $\beta_1=7.5957$, $\beta_2=3.5876$, $\beta_3=1.6382$, and $\beta_4=0.49294$. The same $\epsilon_c$ enters both $A$ and $C$, and we use it for every model of $G(q)$ shown in this work, including those of Refs. [32, 37]. Evaluating $A$ and $C$ with the Perdew–Zunger form [87] instead changes $A$ by less than 0.5% and $C$ by 1–4%. The resulting change in $\chi(q)$ is below 0.3% for $r_s\leq20$, and reaches about one statistical error bar at $r_s=90$.

The constant $B$ is taken from the parametrization of Ref. [32], which is based on the DMC calculations of Ref. [26],

with $a_1^B=2.15$, $a_2^B=0.435$, $b_1^B=1.57$, and $b_2^B=0.409$. (15) was fitted for $r_s\leq10$, so for $r_s>10$ it is an extrapolation, and is not analyzed in this work since the regime it affects is either not covered by our DMC data or the error bars are too large to be sensitive to it.

Determining the perturbation range

In this appendix, we describe our procedure for selecting the range of perturbation amplitudes $v_q$ used to fit the DMC energies $E(v_q)$ and extract $\chi(q)$, which is crucial for obtaining reliable results. The main challenge is to select a range that is sufficiently small to remain in the linear-response regime, but sufficiently large to ensure a good signal-to-noise ratio (SNR) in the DMC energies. At small $v_q$, the energy shift $E(v_q)-E_0\propto v_q^2$ becomes comparable to the statistical error of the DMC energy, so the curvature is poorly determined, which leads to very large error bars in $\chi(q)$. At large $v_q$, the higher-order terms in Eq. (2) become significant and bias the quadratic coefficient, and hence $\chi(q)$. Since each amplitude requires an independent optimization and DMC calculation, a poorly chosen range is also costly. We have found that this range depends on both the density and the wave vector, and we have developed a systematic procedure to determine it. Over the densities considered here, the appropriate amplitudes span almost two orders of magnitude.

Figure 12 demonstrates our strategy for selecting the range of $v_q$ used to fit $E(v_q)$ and extract $\chi(q)$, for a single density and wave vector. The main panel shows the DMC energies (symbols) with the quadratic fit of Eq. (2) (solid). As one can see, one needs to select a range of $v_q$ that is neither too small, where the SNR is low, nor too large, where the response is no longer linear. We indicate these two regions with hatched areas in the main panel. The lower panel shows the SNR (circles) and reduced $\chi_\nu^2$ of the quadratic fit (triangles). The best extraction of $\chi(q)$ is obtained in the linear-response regime, where the SNR is sufficiently high and the reduced $\chi_\nu^2$ of the quadratic fit is sufficiently low. Beyond the linear regime, $\chi_\nu^2$ rises by one to three orders of magnitude within a few $v_q$ points (see lower panel), so any threshold of order unity selects the same range. We have used the criteria of SNR $>10$ and $\chi_\nu^2<5$ to select the range of $v_q$ used for the quadratic fit.

The demonstration in Fig. 12 is done for a single point. Selecting the valid range of $v_q$, on the other hand, has a nontrivial dependence on both $r_s$ and $q$. As we have seen in Figure 1(a) and Figure 1(b), the linear-response regime shrinks significantly as one goes from $r_s=5$ to $r_s=70$. Determining this dependency by brute force is prohibitively expensive. For this, we have solved the noninteracting problem perturbed by the same cosine potential, and found that the maximum $v_q$ allowed to remain in the quadratic regime scales as $r_s^{-p}$ with $p\simeq1.4$–1.5. Since we have found the interacting linear-response window to be smaller than the noninteracting one, we then determined the interacting $v_q^{\max}$ at a single density ($r_s=5$), by sampling $v_q$ between zero and the noninteracting maximum, and we transferred it to the other densities assuming the same $r_s^{-p}$ scaling. In each case, $E(v_q)$ was verified to remain in the quadratic regime (by the criteria defined above). The resulting amplitudes are shown in the right panel of

Figure 12, together with the fit $v_q^{\max} = 0.1212r_s^{-1.43}\,\mathrm{Ha}$, whose exponent is taken from the noninteracting problem and whose prefactor is fitted to the data points. This avoids scanning $v_q$ at every density, and it provides a guide for choosing the perturbation amplitude in future calculations. Note that the window also depends on $q$: around $q \sim2k_F$, the linear-response regime for higher $r_s$ values is significantly smaller compared to other wave vectors, as shown by the dash-dotted lines in Figure 1(a) and Figure 1(b). We handle this by successively discarding the largest perturbation amplitude until the criteria defined above, i.e., SNR $> 10$ and $\chi_\nu^2 < 5$, hold. All results that we present in this work are obtained using this procedure and criteria.

Effects of twist averaging on $\chi(q)$

Twist-averaged boundary conditions [88] are a common alternative to analytic corrections for one-body finite-size errors. We test them at $r_s = 5$ for $N_e = 54$ with $\mathbf{q} = (3,0,0)2\pi/L$ and $N_e = 162$ with $\mathbf{q} = (3,2,2)2\pi/L$. For each twist, drawn from a Halton sequence starting at $\Gamma$, we perform DMC calculations at fixed $N_e$ and at the production perturbation amplitudes, average $E(v_q)$ over twists, and extract $\chi(q)$ from its curvature. Figure 13 shows the result as a function of the number of twists, together with the $\Gamma$-point value corrected with Eq. (6). With increasing $N_{\mathrm{tw}}$, the twist average moves away from the uncorrected $\Gamma$ value ($N_{\mathrm{tw}} = 1$) toward the corrected one, recovering most of the correction. For $N_e = 162$, $-\chi/n_0$ goes from 10.84 at $\Gamma$ to $11.48 \pm0.13$ after 32 twists, compared with $11.76 \pm0.05$ for the corrected $\Gamma$ value. For $N_e = 54$, the uncorrected $\Gamma$ value is already closer to the corrected one, and the twist average remains within $2\sigma$ of it for all $N_{\mathrm{tw}} \ge2$, reaching $12.64 \pm0.16$ after 96 twists, compared with $12.29 \pm0.16$. In both cases, the twist average agrees with the corrected $\Gamma$ value to within 2–3%, i.e., within about $2\sigma$. The remaining difference is expected for fixed-$N_e$ twist averaging. For a size-independent $G(q)$ [Figure 3(a)], the twist dependence of $\chi(q)$ enters only through the noninteracting response at each twist, so the procedure can reach the thermodynamic limit only if it recovers the Lindhard function for noninteracting electrons. At fixed $N_e$ it does not, because each twist fills its own $N_e$ lowest levels and the averaged Fermi sphere is distorted [88]. Fixing the chemical potential instead avoids this at the cost of a fluctuating particle number (grand-canonical twist averaging), which could in principle lead to a more accurate size correction but was not tested here. Throughout this work, we keep $N_e$ fixed for all twists (canonical twist averaging) and use Eq. (6), which inserts the thermodynamic-limit $\chi_0$ directly from a single twist.

References