*Ab initio* quantum Monte Carlo methods, in principle, allow for the calculation of exact properties of correlated many-electron systems but are, in general, limited to the simulation of a finite number of electrons *N* under periodic boundary conditions. Therefore, an accurate theory of finite-size effects is indispensable to bridge the gap to realistic applications in the thermodynamic limit. In this work, we revisit the uniform electron gas at finite temperature, as it is relevant to contemporary research, e.g., in the field of warm dense matter. In particular, we present a new scheme to eliminate finite-size effects both in the static structure factor *S*(*q*) and in the interaction energy $v$, which is based on the density response formalism. We demonstrate that this method often allows us to obtain $v$ in the thermodynamic limit within a relative accuracy of ∼0.2% from as few as *N* = 4 electrons without any empirical choices or knowledge of results for other values of *N*. Finally, we evaluate the applicability of our method upon increasing the density parameter *r*_{s} and decreasing the temperature *T*.

## I. INTRODUCTION

The accurate numerical simulation of interacting many-electron systems constitutes a central challenge in many domains of theoretical physics, quantum chemistry, material science, and related disciplines.^{1} Yet, while the equations governing these systems have been known for approximately a century, there still does not exist a practical tool to solve them in many situations and one has to rely on approximations.

A particularly important method for this purpose is density functional theory, which has emerged as the de facto work horse of quantum chemistry and boasts a remarkable success regarding the description of both bulk materials and molecules.^{2,3} More specifically, density functional theory simulations are computationally efficient, as the complicated *N*-body problem of interest is mapped on an effective single-body problem that can actually be solved numerically. Yet, this mapping comes at the expense of self-sufficiency, as information about electronic correlations must be supplied externally, in general, in the form of the unknown exchange–correlation functional.

In contrast, computationally more expensive quantum Monte Carlo (QMC) methods are, at least in principle, capable of providing an exact solution to the *N*-electron problem^{4,5} without any *a priori* information or empirical input. In particular, it was only the seminal QMC study of the uniform electron gas (UEG) by Ceperley and Alder,^{6} which allowed for the subsequent construction of accurate electron exchange–correlation functionals,^{7–10} that made the success of density functional theory possible in the first place. At the same time, most QMC methods are, by definition, only applicable for a finite number of electrons *N*, while most physical systems are given in the thermodynamic limit, i.e., in the limit of *N* → *∞* with the density *n* = *N*/*V* remaining constant. This problem is substantially exacerbated by the notorious fermion sign problem,^{11–13} which leads to an exponential increase in the computation time with *N* (see Ref. 13 for a topical review article). In fact, the sign problem has been revealed as *NP*-hard for a certain type of Hamiltonian by Troyer and Wiese,^{12} and a general solution appears to be improbable.

For these reasons, the solid understanding of finite-size effects in QMC data and their subsequent estimation in the form of a finite-size correction (FSC) is of paramount importance for the QMC community and beyond and constitutes a highly active topic of research.^{14–33} Moreover, while most works have been devoted to the study of electrons under ambient conditions, i.e., in the ground state,^{34} there has recently emerged a growing interest in the properties of matter at extreme densities and temperatures. Of particular importance is the regime of the so-called warm dense matter,^{27,35,36} which is characterized by the simultaneous importance of thermal excitations, Coulomb correlations, and fermionic exchange effects and naturally occurs in astrophysical objects such as giant planet interiors^{37} and brown dwarfs.^{38,39} In addition, warm dense matter has been predicted to occur on the pathway toward inertial confinement fusion^{40} and is routinely realized experimentally in large research centers around the globe.^{41}

This interest has sparked a series of new developments in the field of electronic QMC simulations at finite temperature,^{26,32,33,42–57} which, in turn, has caused the need to understand the impact of thermal excitations on finite-size effects. In general, this problem can be restated as the search for a short-range property that can be accurately inferred from a QMC simulation of a finite system and, in combination with a readily available theory such as the random phase approximation (RPA) [see Eq. (1) in the following], yields the full description of a system in the thermodynamic limit.

In the ground state, Chiesa *et al.*^{18} have proposed using the static structure factor *S*(*q*) for this purpose, which indeed has been shown empirically to only weakly depend on *N* at both zero and finite temperature.^{26,27} This finding has allowed us to introduce a simple analytical first-order correction to the finite-size error of the interaction energy per particle $v$, which was subsequently generalized by Brown *et al.*^{32} to arbitrary temperatures. Yet, this correction breaks down for both high temperature and density and, thus, is inapplicable over substantial parts of the warm-dense matter regime. The full estimation of the main contribution to the finite-size error of $v$ then requires the estimation of a discretization error [see Eq. (20) below] using a suitable trial function such as the static structure factor within RPA or a more sophisticated dielectric theory such as the approximate scheme by Singwi, Tosi, Land, and Sjölander^{58–60} (STLS). This was demonstrated by Dornheim, Groth, and co-workers,^{26,27,44} who showed that this procedure reduces the finite-size error in the QMC data by approximately two orders of magnitude, allowing for a reliable subsequent extrapolation to the thermodynamic limit.^{26,61}

In the present work, we go one step further and address the source of the residual finite-size errors after the discretization error has been eliminated. To this end, we employ the density response formalism and propose that in many cases, the static local field correction constitutes a more suitable choice for a short-range exchange–correlation function that can be accurately estimated from a QMC simulation of as few as *N* = 4 particles. More specifically, our new scheme allows for a direct estimation of the finite-size error in *S*(*q*) itself, which already constitutes an important finding in itself. Moreover, this FSC for the static structure factor can subsequently be used to eliminate the residual finite-size error in $v$, and we find that often *N* = 4 particles are sufficient to estimate $v$ with a relative accuracy of ∼0.2% without any empirical input, extrapolation, or knowledge about QMC results for multiple *N*.

This paper is organized as follows: In Sec. II, we introduce the required theoretical background starting with the density response formalism (Sec. II A) and the fluctuation–dissipation theorem (Sec. II B), followed by our new theory of finite-size effects in the static structure factor (Sec. II C), and the implications for the FSC of the interaction energy $v$ (Sec. II D). The application of our scheme is demonstrated in Sec. III, starting with an extensive analysis at extreme density and temperature in Sec. III A. In addition, we analyze the applicability of the method for the important cases of increasing coupling strength (Sec. III B) and low temperature (Sec. III C). This paper concludes with a brief summary and outlook in Sec. IV.

## II. THEORY

Throughout this work, we restrict ourselves to an unpolarized uniform electron gas (UEG; see Ref. 27 for details). In addition, the UEG is characterized by two parameters: (i) the density parameter $rs=r\xaf/aB$, where $r\xaf$ and *a*_{B} denote the average particle distance and first Bohr radius, and (ii) the degeneracy temperature *θ* = *k*_{B}*T*/*E*_{F}, where *E*_{F} is the usual Fermi energy.^{1,62} Recalling the typical definition of a coupling parameter Γ as the ratio of interaction ($W\u223crs\u22121$) to kinetic energy ($K\u223crs\u22122$), it is easy to see that Γ ∼ *r*_{s}. Therefore, *r*_{s} constitutes a direct measure for the coupling strength of the system.

All formulas and results are given in Hartree atomic units.

### A. Density response and local field correction

The density response of an electron gas to an external harmonic perturbation^{53} of wave number *q* and frequency *ω* is—within linear response theory—fully described by the dynamic density response function,^{1,63}

Here, *χ*_{0}(*q*, *ω*) denotes the density response function of the ideal Fermi gas and the information about exchange–correlation effects is contained in the dynamic local field correction *G*(*q*, *ω*).

Let us next consider the static limit, i.e.,

In this limit, accurate data for Eq. (1) have been presented by Dornheim *et al.*^{64–66} based on the relation^{67}

with the imaginary-time density–density correlation function being defined as

In addition, it is possible to obtain *χ*(*q*) from a simulation of the harmonically perturbed electron gas, which has been done both in the ground state^{28,67,68} and at finite temperature.^{33,51,53} In principle, it is then straightforward to use *χ*(*q*) to solve Eq. (1) for the static local field correction,

For the present paper, it is important to explicitly indicate the dependence of different estimated properties on the system size, the purpose for which we shall henceforth use the superscript *N*. In particular, the simplest way to write Eq. (5) is then given by

Yet, we empirically know that *χ*^{N}(*q*) substantially depends on *N*, and this finite-size error is propagated into *G*^{N}(*q*) as defined in Eq. (6). Fortunately, we also know that this error is almost exclusively due to the inconsistency between $\chi 0TDL(q)$ and *χ*^{N}(*q*), as the consistent determination of *G*(*q*) would require us to instead use the density response function of the ideal system at the same system size, $\chi 0N(q)$. Thus, the finite-size corrected static local field correction can be computed as^{28,68}

Moreover, we can readily define a FSC for *G*(*q*),

such that

### B. Fluctuation–dissipation theorem

The fluctuation–dissipation theorem^{1}

relates Eq. (1) to the dynamic structure factor *S*(*q*, *ω*) and, thus, directly connects the local field correction to different material properties. In particular, the static structure factor is defined as the normalization of the DSF,

and thus entails an averaging over the full frequency range. We stress that this is in contrast to the static density response function *χ*(*q*) introduced in Sec. II A, which is defined as the limit of *ω* → 0. The static structure factor, in turn, gives direct access to the interaction energy of the system, and for a uniform system, it holds that^{27}

Finally, we mention the adiabatic connection formula,^{27,61,69}

which implies that the free energy *f* (also, equivalently, the partition function *Z*) can be inferred from either *χ*(*q*, *ω*), *S*(*q*, *ω*), or *S*(*q*).

### C. Finite-size correction of *S*(*q*)

Since the full frequency dependence of *G*(*q*, *ω*) remains unknown in most cases, one might neglect dynamic effects and simply substitute *G*(*q*) in Eq. (1). This leads to the dynamic density response function within the *static approximation*,^{70–73}

which entails the frequency dependence on an RPA level, but exchange–correlation effects are incorporated statically. Indeed, it has been shown that Eq. (14) allows for an accurate (although not exact) description of many material properties.^{70,71,73,74} This finding constitutes the motivation for the present scheme to overcome finite-size effects in *S*(*q*) and related quantities.

Very recently, Dornheim and Moldabekov^{75} have introduced the concept of an effectively frequency-averaged local field correction $G\xaf(q)$, i.e., a static local field correction that is not defined in the limit of *ω* → 0 but effectively incorporates the dependence on *ω*. For example, we can define such a local field correction in a way that, when it is being inserted into Eqs. (10), (11), and (14), it exactly reproduces the QMC data for *S*^{N}(*q*),

In practice, we determine $G\xafinvertN(q)$ numerically by scanning the corresponding results for *S*(*q*) over a dense grid of $G\xaf$ for each particular wave number *q*. In fact, it has been shown^{74} that *G*^{N}(*q*) is almost exactly equal to $G\xafinvertN(q)$ for wave numbers *q* ≲ 3*q*_{F}. It is therefore reasonable to approximate the (*a priori* unknown) FSC for $G\xafinvertN(q)$ by the FSC for the static limit Δ*G*^{N}(*q*) defined in Eq. (9) above,

The resulting corrected value $G\xafinvertFSC(q)$ is then used to compute the finite-size corrected value of the static structure factor $SNFSC(q)$, again via Eqs. (10), (11), and (14). The corresponding value of the FSC for *S*^{N}(*q*) is then given by

We stress that our approach is inherently superior to the *static approximation* introduced in Refs. 70, 71, and 74, as we only use it to estimate the finite-size error in $G\xafinvertN(q)$, which is comparably small. Thus, the full frequency dependence of *G*(*q*, *ω*) is consistently incorporated into the definition of $G\xafinvertN(q)$, whereas we only neglect it for the FSC in Eq. (16).

### D. Finite-size correction of the interaction energy

The interaction energy in the thermodynamic limit is defined by evaluating the integral in Eq. (12) above using the static structure factor in the thermodynamic limit, *S*^{TDL}(*q*),

In contrast, the corresponding energy per particle for a finite system is defined as a sum over the reciprocal vectors **G** of the simulation cell,

where *ξ*_{M} is the usual Madelung constant taking into account the interaction of a charge with its own background and periodic array of images.^{14} We note that the different power in *q* in Eq. (18) compared to Eq. (19) is a direct result of the transformation in spherical coordinates and the resulting one-dimensional integral over the modulus of **q**.

Following a short analysis that has been presented elsewhere,^{18,24,26,27} it can be seen that there are two potential sources for the difference between Eqs. (18) and (19): (i) the discretization error $\Delta vdN$ due to the approximation of a continuous integral by a discrete sum and (ii) the error intrinsic in *S*^{N}(*q*) ≠ *S*^{TDL}(*q*), leading to $\Delta viN$.

It is well known that (i) is the dominant effect (at least for a sufficiently large system-size *N*), and it can be accurately estimated using a decent trial function for *S*(*q*), e.g., using STLS,^{26,59,60} or even RPA,

In addition, we mention that the first-order contribution to Eq. (20) is a direct consequence of the missing **G** = **0** term in Eq. (19) and can be exactly evaluated in RPA, which gives^{18,27,32}

where $\omega p=rs3/3$ denotes the usual plasma frequency.

Yet, often (ii), too, is not negligible, which is particularly true either at high temperature/small *r*_{s} or in the low-temperature regime where *S*^{N}(*q*) is afflicted with momentum shell effects.^{31} Moreover, the fermion sign problem^{12,13} still prevents QMC simulations except for very small *N* for substantial parts of the WDM regime, which makes the accurate estimation of the intrinsic error highly desirable.

We thus make use of the new FSC for *S*(*q*) given in Eq. (17) above and write the intrinsic finite-size error of $v$ as

The fully finite-size corrected estimate for the interaction energy is then given by

### E. Momentum shell effects

At a low temperature, when the momentum distribution approaches a step function, momentum shell effects constitute an additional source of finite-size errors. More specifically, momentum shell effects are most pronounced when the number of electrons of a certain spin does not allow to exactly fill a given shell in **q**-space. Therefore, one typically considers only particle numbers that are commensurate to the respective cubic lattice, i.e., *N* = 14, 38, 54, 66, … for an unpolarized system. Empirically, it was found^{33} that such effects start to manifest for *θ* ≲ 0.5. In contrast, these effects are smoothed out at higher temperatures due to thermal broadening of the occupation numbers.

A common way to remove such momentum shell effects is to apply the well-known twist-averaging procedure.^{17,20,21,76,77} The basic idea is to use generalized boundary conditions, where the momentum grid is shifted by a finite twist angle $t\u2208R3$, so that the quantized wave vectors are given by **q** = 2*πL*^{−1}**l** + **t** with $l\u2208Z3$ and to subsequently average over different **t**.

Yet, as the present manuscript is devoted to relatively large temperatures, all of our QMC results have been obtained without any twist-averaging. Furthermore, we mention that the momentum-shell effects are fully present in $\chi 0N(q)$ too (see Ref. 33 for a corresponding analysis), and the FSC presented in Secs. II A–II D thus potentially constitutes an alternative to the twist-averaging procedure. This possibility is addressed in Sec. III C.

## III. RESULTS

All results presented in this work have been obtained for the unpolarized UEG, i.e., with the number of spin-up and spin-down electrons being equal.

### A. High densities and temperatures

Let us start our investigation with an example at high density (*r*_{s} = 0.3) and extreme temperature (*θ* = 2), which is close to the conditions in the solar core.^{78} In Fig. 1, we show the corresponding dependence of the interaction energy per particle *V*^{N}/*N* on the system size *N*. Let us first consider the bottom panel, where we show the entire relevant energy range. In particular, the green crosses depict the raw uncorrected QMC data [configuration path integral Monte Carlo (CPIMC)^{79,80} data partly taken from Ref. 26 and a few data points obtained with the standard path integral Monte Carlo (PIMC)^{81–83}] that have been computed from Eq. (19). Evidently, these data exhibit a strong dependence on *N* reaching ∼100% for *N* = 4. Even for the largest depicted system size, *N* = 300, finite size errors have not completely vanished. The yellow triangles have been obtained by adding to the raw data the FSC for the discretization error using *S*^{STLS}(*q*) as a trial function [see Eq. (20) above]. This leads to a remarkable reduction of the finite-size errors in the dataset, and the residual error appears to be of the order of ∼1% even for as few as *N* = 4 particles. This validates the previous conclusion from Ref. 26 that $\Delta vdN$ does indeed constitute the dominant contribution to the total finite-size error Δ$v$^{N}. For completeness, we also show the first-order correction $\Delta vd,0N$ [Eq. (21)], which is depicted by the dark gray stars. Yet, this first-order correction is only reasonable for very large numbers of electrons *N* and even leads to an increase in finite-size errors for *N* ≲ 20. Finally, the blue dotted line depicts a linear fit to the yellow triangles for *N* ≥ 20, which has been used to extrapolate the residual error of these data in Refs. 26, 27, 44, and 61.

Let us postpone the discussion of the red circles in Fig. 1 and instead focus on the source of the remaining *N*-dependence after the discretization error has been eliminated. To this end, we show QMC results for the static structure factor in Fig. 2 under the same conditions for three values of *N*. More specifically, the colored symbols depict the raw, uncorrected QMC results for *S*^{N}(*q*) for *N* = 4 (blue diamonds), *N* = 20 (red circles), and *N* = 66 (green crosses). First and foremost, we note that the dependence of *S*^{N}(*q*) on *N* is indeed small when compared to the substantially more pronounced finite-size error of *V*^{N}/*N* itself. Still, there do appear significant differences between the data for different *N*, which is particularly pronounced for *N* = 4, as expected.

Adding the FSC Δ*S*^{N}(*q*) from Eq. (17) to *S*^{N}(*q*) leads to the black symbols in Fig. 2. Evidently, this correction works well under these conditions, and the black diamonds, squares, and crosses can hardly be distinguished with the naked eye. This, in turn, means that as few as *N* = 4 particles are sufficient to obtain a wave-number resolved description of electronic correlations in the UEG, despite the apparently large finite-size effects at such an extreme density.

Having thus verified our theory for finite-size errors in the static structure factor, we can readily compute the second source of finite-size errors in the interaction energy, namely, the intrinsic error $\Delta viN$ [Eq. (22)]. The results of removing both intrinsic and discretization errors are depicted by the red circles in Fig. 1. In particular, our new procedure leads to a striking reduction in finite-size errors even for *N* = 4 electrons without any empirical input or knowledge about the behavior for other *N*. Furthermore, the corrected data for all *N* are within ∼0.2% of the value that was obtained from the extrapolation of the yellow triangles. For comparison, we have also included the prediction for $v$ from STLS^{59,60} (horizontal gray line) that is obtained by inserting *S*^{STLS}(*q*) in Eq. (18). Even though STLS and related dielectric theories are often considered to be rather accurate under these conditions,^{66} we can clearly resolve a systematic overestimation of $v$ by ∼1%. Yet, this deficiency has no impact on the estimation of $\Delta dN[SSTLS(q)]$, since the systematic error in *S*^{STLS}(*q*) cancels in Eq. (20).

To further illustrate the idea behind our new FSC scheme, we show the density response function of the corresponding non-interacting (ideal) Fermi system in Fig. 3. In particular, the black solid curve shows $\chi 0TDL(q)$, and the blue diamonds, red circles, and green crosses depict $\chi 0N(q)$ for *N* = 4, *N* = 20, and *N* = 66, respectively. Evidently, there appear significant deviations between $\chi 0TDL$ and $\chi 0N$ that closely resemble the trend observed in *S*^{N}(*q*) shown in Fig. 2 above. In addition, we have also included *χ*^{STLS}(*q*) in the figure (dashed-dotted yellow curve), which serves as a guide for the eye.

To conclude our analysis for the present example, we show the actual effectively static local field correction $G\xafinvertN(q)$ in Fig. 4 for the same particle numbers. The top panel shows the uncorrected raw data, and we see substantial finite-size effects for all datasets exceeding 100% for *N* = 4. The bottom panel shows the corresponding results for $G\xafFSC(q)$ [see Eq. (16)], and we find that the correction Δ*G*^{N}(*q*) defined in Eq. (8) removes nearly the entire dependence on *N*. This finding can be interpreted in the following way: the local field corrections *G*(*q*, 0) and $G\xafinvert(q)$ have been defined as the difference between the exact density response/static structure factor of the system and the prediction of mean-field theory, i.e., RPA. Therefore, it is a short-range function taking into account exchange–correlation effects, which are hardly affected by the system size at these conditions. Correspondingly, we can accurately extract either *G* or $G\xafinvert$ from a single simulation of the UEG even for *N* = 4. Finally, this information about the short-range behavior can be combined with the mean-field description, thus yielding the full wave number resolved description of the system in the thermodynamic limit.

To discern the limits of our new FSC scheme, we next consider a substantially more extreme case, where even thousands of particles are not sufficient to resemble the thermodynamic limit. In particular, we show the interaction energy per particle in Fig. 5 for *r*_{s} = 0.1 and *θ* = 8. Let us first consider the bottom panel showing the entire relevant energy range of the uncorrected data, which again have been partly taken from Ref. 26 and computed with CPIMC and partly obtained using standard PIMC in the context of the present work. The raw QMC results (green crosses) naturally exhibit severe finite-size effects that exceed 500% for *N* = 4. Moreover, the dependence of *V*^{N}/*N* is substantial even for *N* = 2000 (the largest *N* shown in Fig. 5) and does not follow any obvious functional form. The yellow triangles have again been obtained by removing the discretization error $\Delta vdN$ estimated using *S*^{STLS}(*q*), which reduces the dependence on *N* by at least two orders of magnitude. Yet, there still remains a bias of several percent for the smallest depicted values of *N*. The first-order correction $\Delta vd,0N$, on the other hand, is not appropriate even for *N* ∼ 10^{3} and only becomes accurate for even larger *N*. Finally, the red circles have been computed by subtracting from the yellow triangles the finite-size error due to the intrinsic dependence on *N* of *S*^{N}(*q*), $\Delta viN$, computed with our new scheme. Even for this most extreme case, where thousands of particles do not suffice for an appropriate description, we can predict the interaction energy in the thermodynamic limit with a relative accuracy of ∼0.2% from just *N* = 4 electrons.

### B. Increasing the coupling strength

To further assess the remarkable performance of our new FSC scheme demonstrated at high densities in Sec. III A, we will now investigate the manifestation of finite-size effects upon increasing the density parameter *r*_{s}. For a quantum system such as the UEG, this is equivalent to increasing the coupling strength (see Refs. 1 and 62). A first example is shown in Fig. 6, where we plot the interaction energy per particle for *r*_{s} = 1 and *θ* = 2. First and foremost, we find that finite-size errors are overall reduced in magnitude compared to *r*_{s} = 0.3 shown in Fig. 1, as expected. More specifically, it is known empirically that the system-size dependence increases with the density when *N* is being kept constant^{27} and, conversely, decreases with the increase in *r*_{s}. In addition, the discretization error again accounts for the bulk of finite-size effects, as the yellow triangles constitute a substantial improvement over the green crosses. Furthermore, the purple squares taken from Ref. 26 have been obtained using the Permutation Blocking PIMC (PB-PIMC) method^{45,50,84,85} and, too, were corrected using $\Delta vdN[SSTLS(q)]$. These data are in excellent agreement with the independent standard PIMC data obtained for the present work, which further corroborates the high quality of the existing description of the UEG.^{27,61} At the same time, the first-order correction $\Delta vd,0N$ is not sufficiently accurate to replace the full evaluation of Eq. (20). Finally, the new FSC scheme for the intrinsic error $\Delta viN$ results in the red circles, where the residual error does not exceed ∼0.3% even for *N* = 4.

Let us next further increase the density parameter to *r*_{s} = 2, which is shown in Fig. 7. This is a typical density for the conduction band electrons in metals, and thus, such a case is of paramount importance. In warm dense matter research, such states are created in many experiments, e.g., using aluminum.^{86–89} In this case, finite-size effects are again reduced in relative importance, and the first-order correction from Eq. (21) becomes accurate for *N* ∼ 34. Yet, there still remain significant residual errors in *V*/*N* after the discretization error is removed, which are of the order of a few percent for *N* = 4 and decrease to ∼0.5% for *N* = 20. While it is certainly feasible to carry out a controlled extrapolation of the remaining *N*-dependence of the data under these conditions, it is still worth investigating if our new scheme makes this effort superfluous.

Thus, we have added our estimation of $\Delta viN$ to the yellow triangles, and the results are as usual depicted by the red circles. We find that our scheme significantly reduces the dependence on *N* for all numbers of electrons, although the actual error seems to be overestimated by Eq. (22). For example, the thus fully corrected value for the interaction energy per particle for *N* = 4 deviates by almost 1% from the extrapolated result (horizontal blue dotted line), and these residual deviations decrease with the increase in *N*. We thus conclude that our new scheme does indeed constitute an improvement over the previous state-of-the-art of FSCs, but the finite-size error cannot be completely removed for *N* = 4.

In order to understand this observed trend, it makes sense to go to even larger values of the coupling strength. To this end, we show the interaction energy per particle for *r*_{s} = 5 and *θ* = 2 in Fig. 8. We note that such a low density serves as a valuable laboratory for the impact of electronic exchange–correlation effects on material properties^{90,91} and can be realized experimentally, e.g., via hydrogen jets.^{92} First and foremost, we observe a further reduced relative magnitude of finite-size effects, as shown by the green crosses depicting the raw PIMC data for *N* ≥ 20. Moreover, both the full expression for $\Delta vdN$ and the first-order term $\Delta vd,0N$ give fairly similar results and are all but indistinguishable for *N* ≳ 10. Moreover, the residual error after removing the discretization error does not exceed 0.3% in this case. Finally, the red circles have been obtained by further adding to the yellow triangles the result for $\Delta viN$. Evidently, Eq. (22) substantially overestimates true residual error and, in fact, exhibits an even somewhat more pronounced and less trivial dependence on *N*. The rest of this section will attempt to explain this unexpected breakdown of our expression for $\Delta viN$ in some detail.

Let us begin this analysis by considering the central ingredient to the interaction energy, i.e., the static structure factor *S*(*q*) shown in Fig. 9. The top panel shows standard PIMC results for the raw uncorrected *S*^{N}(*q*) for *N* = 4 (blue diamonds), *N* = 14 (red circles), *N* = 20 (green crosses), and *N* = 66 (gray stars). In addition, the dashed-dotted yellow curve corresponds to *S*^{STLS}(*q*) and has been included as a reference. Interestingly, the only substantial intrinsic dependence of *S*^{N}(*q*) on *N* appears for *N* = 4, whereas the other data seem to fall on a continuous curve.

The bottom panel of the same figure shows the same information, but for the finite-size corrected version of the static structure factor, $SFSCN(q)$. In this case, we find that adding our expression for Δ*S*^{N}(*q*) from Eq. (17) onto *S*^{N}(*q*) leads to a substantial increase in finite-size effects and we find a significant difference even between *N* = 20 and *N* = 66 that was absent in the top panel. We have thus found that our new scheme for the FSC for $v$ breaks down because the involved expression for Δ*S*^{N}(*q*) becomes inappropriate.

Proceeding to the next deeper level of explanation, we show different local field corrections in Fig. 10 for the same conditions. The top panel corresponds to the effectively static local field correction $G\xafinvert(q)$ defined in Eq. (15) above, and the blue diamonds, red circles, and green crosses show the raw uncorrected data for *N* = 4, *N* = 14, and *N* = 20, respectively. In addition, the yellow dashed-dotted line shows the local field correction from STLS and has been included as a reference. Similar to our previous observation regarding *S*^{N}(*q*) above, here too, we only find small finite-size effects for *N* = 4, and no dependence on *N* is visible for *N* = 14 and *N* = 20. The black symbols in the same panel have been obtained by adding to $G\xafinvertN(q)$ the finite-size correction for the static local field correction *G*^{N}(*q*, 0), Δ*G*^{N}(*q*), given in Eq. (8). Again, the FSC exacerbates the actual *N*-dependence of the uncorrected data.

This can be potentially explained by two distinctly different effects: (i) the FSC for the static local field correction *G*^{N}(*q*, 0) could be fundamentally different from the real *a priori* unknown FSC for $G\xafinvertN(q)$, or (ii) our expression for Δ*G*^{N}(*q*) is inappropriate for both *G*^{N}(*q*, 0) and $G\xafinvertN(q)$. This question is resolved by the analysis presented in the bottom panel of Fig. 10, where we show the exact static limit of the local field correction for the same conditions as in the top panel. Evidently, the different data for *G*(*q*, 0) exhibit the same trend as $G\xafinvertN(q)$. In particular, the FSC Δ*G*^{N}(*q*) does not constitute an improvement over the uncorrected PIMC data, which means that explanation (ii) holds: we do not have an appropriate theory for finite-size effects in both *G*^{N}(*q*, 0) and $G\xafinvertN(q)$ under these conditions.

This can be understood more intuitively by examining Fig. 11, where we show data for different density response functions *χ*(*q*). More specifically, the colored symbols depict the results for $\chi 0N(q)$, and the dark gray solid curve depicts $\chi 0TDL(q)$. We note that there are substantial finite-size effects, particularly for *N* = 4, as expected. In contrast, the black symbols depict the static density response function of the UEG that was computed using standard PIMC via Eq. (3), and no dependence on *N* can be resolved with the bare eye. We have thus found that finite-size effects in the actual density response function *χ*^{N}(*q*) [also in *S*^{N}(*q*)] disappear with increasing *r*_{s}, which does not hold for $\chi 0N(q)$. Yet, if it is *χ*^{N}(*q*) ≈ *χ*^{TDL}(*q*), then using $\chi 0TDL(q)$ in Eq. (6) is actually consistent and directly leads to the static local field correction in the thermodynamic limit, *G*^{N}(*q*) ≈ *G*^{TDL}(*q*), just as we have observed in Fig. 10.

The spurious effects of the finite-size correction Δ^{N}*G*(*q*) can then finally be interpreted in the following way: The mean-field description that is directly based on $\chi 0N(q)$ predicts substantial finite-size effects. These, however, are absent from the actual density response function *χ*^{N}(*q*), which we computed via PIMC. Therefore, we can conclude that the system-size dependence is suppressed by the increased coupling strength, since the particles are considerably more localized due to the Coulomb repulsion. It is then only meaningful to define the local field correction as the difference between $\chi 0TDL(q)$ and *χ*^{N}(*q*), as the finite-size effects in $\chi 0N(q)$ do not manifest in the latter. If we instead would define *G*^{N}(*q*, 0) as the difference between *χ*^{N}(*q*) and $\chi 0N(q)$, this local field correction would have to balance the system-size dependence from $\chi 0N(q)$ so that the actual density response function remains independent of *N*. This then explains the unreasonable behavior of Δ*G*^{N}(*q*) shown in Fig. 10.

In a nutshell, our new FSC scheme implicitly assumes that the description of the system can be decomposed into a mean-field part that exhibits finite-size effects and a short-range local field correction that does not. At large values of *r*_{s}, this assumption breaks down, as the finite-size effects in $\chi 0N(q)$ are suppressed by the strong Coulomb repulsion. At the same time, we note that this is not catastrophic, as the dependence of QMC data on *N* is small for lower densities in the first place.^{26,27}

### C. Going to low temperature

The final part of the investigation in this work is devoted to the performance of our new FSC scheme for low temperatures. This is a particularly interesting regime, as accurate QMC data for the UEG are sparse under these conditions due to the fermion sign problem.^{13,44,55,57} In addition, the UEG has attracted renewed interest at high densities also in the ground state,^{93,94} where the full configuration interaction QMC (FCIQMC) method^{95} is capable of giving accurate results for finite *N*.

One of the main sources of finite-size errors at low temperature is given by momentum-shell effects, i.e., if the number of electrons of a certain spin does not exactly fill a corresponding shell in **q**-space (i.e., 1, 7, 19, 27, 33, …). For this reason, one typically selects a system size that is commensurate to the grid in momentum space, e.g., for an unpolarized UEG, *N* = 14, 38, 54, 66, …. An additional strategy is to employ periodic boundary conditions with a finite twist angle^{17,20,21,76,77} and perform a subsequent *twist-averaging* (see Sec. II E), which has been shown to greatly alleviate these effects. At the same time, we note that these momentum shell effects are already fully present in $\chi 0N(q)$, which makes our new FSC a potential alternative in this regime. This is illustrated in Fig. 12, where we show the density response function for *r*_{s} = 0.3 and *θ* = 0.5. Here, the blue diamonds, red circles, and green crosses show $\chi 0N(q)$ for *N* = 4, *N* = 6, and *N* = 34, respectively, and the dark gray solid curve corresponds to $\chi 0TDL(q)$. In addition, the dashed-dotted curve shows *χ*^{STLS}(*q*) and has been included as a reference. First, we note that *θ* = 0.5 only constitutes the boundary of the low-temperature regime, and thermal excitations still play an important role.^{27} Still, incipient momentum shell effects clearly manifest for $\chi 0N(q)$ in Fig. 12, particularly for *N* = 4. Similar results for lower temperatures have been presented by Groth *et al.*^{33}

The next step toward a finite-size correction for *S*^{N}(*q*) and *V*^{N}/*N* is given by the analysis of different local field corrections presented in Fig. 13 for the same conditions. Let us first consider the top panel where we show results for *N* = 4 (diamonds) and *N* = 6 (circles) obtained without any twist-averaging. The red data show $G\xafinvertN(q)$, which constitutes the basis of our FSC scheme for Δ*S*^{N}(*q*) and, in turn, $\Delta viN$. Apparently, the datasets for the two different particle numbers do not follow a smooth progression, and a finite-size correction is needed; similarly, *S*^{N}(*q*) exhibits finite-size effects (see Fig. 14 in the following). The blue data show results for the static limit of the local field correction *G*^{N}(*q*, 0) and have been obtained from PIMC simulations via Eqs. (3) and (6). On the one hand, these data, too, do exhibit substantial finite-size effects, which further underline the need for an FSC. On the other hand, the system-size dependence of *G*^{N}(*q*) does not resemble the observed behavior of $G\xafinvertN(q)$, which makes the applicability of our FSC scheme questionable.

Remarkably, adding the FSC Δ*G*^{FSC}(*q*) to the blue points gives the black curves, for which no dependence on *N* can be resolved within the given level of statistical uncertainty. This can be seen particularly well in the bottom panel where we show a magnified view around $GNFSC(q)$. Thus, QMC simulations of *N* = 4 and *N* = 6 electrons appear to give us access to *G*(*q*, 0) in the thermodynamic limit without any twist-averaging, but not $G\xafinvertTDL(q)$ as no theory for $\Delta G\xafinvertN(q)$ is at hand.

Let us next ponder the consequences of these findings for the FSC of the static structure factor *S*(*q*), which we show in Fig. 14 for the same conditions. Here, the blue diamonds, yellow triangles and green crosses show the raw uncorrected QMC data for *S*^{N}(*q*). We note that the observed dependence on *N* is qualitatively similar to the finite-size effects observed in $\chi 0N(q)$ in Fig. 12 above. In addition, the red diamonds have been obtained by using *G*^{N}(*q*, 0) to compute the static structure factor within the *static approximation* [Eq. (14)]. Evidently, these data do substantially disagree with the actual QMC data for *S*^{N}(*q*), which appear to indicate that the *static approximation* is not applicable under these conditions. This is not surprising and directly explains the difference between $G\xafinvertN(q)$ and *G*^{N}(*q*) in Fig. 13.

Yet, the full picture is even more subtle and can be deduced from the black diamonds in Fig. 14, which show the static structure factor $SstatFSC(q)$ obtained via the *static approximation* using as input the finite-size corrected local field correction $GNFSC(q)$ for *N* = 4. Remarkably, these data do not exhibit the unsmooth behavior of *S*^{N}(*q*) associated with the momentum shell effects but a smooth progression with *q*. Moreover, these data are much closer to the green crosses depicting the QMC data for *N* = 34 and the black crosses that have been obtained by applying the usual FSC Δ*S*^{N}(*q*) [Eq. (17)] to the latter.

These findings indicate that the *static approximation* does not break down under these conditions after all, since the black diamonds appear to be accurate. Instead, the large discrepancy between the blue and red diamonds follows from an inconsistent mixing of different response functions, which did not matter for *θ* = 2, but has a large impact in the low temperature regime. More specifically, a consistent *static approximation* for a finite number of electrons *N* should have the form

where the dynamic density response function of the ideal Fermi gas is computed for the same system size. Thus, combining $GNFSC(q)$ with $\chi 0TDL(q,\omega )$ gives reasonable results for *S*(*q*) (black diamonds), whereas mixing *G*^{N}(*q*) with $\chi 0TDL(q,\omega )$ does fit neither to the thermodynamic limit nor to the QMC data for *S*^{N}(*q*). Following this line of thinking, our definition of $G\xafinvertN(q)$ from Eq. (15) above, too, is inconsistent, as it must reproduce *S*^{N}(*q*) using the inappropriate thermodynamic limit result $\chi 0TDL(q,\omega )$. Instead, one should define a modified effective local field correction, where the functional $S[G\xaf(q)](q)$ is evaluated using the *static approximation* from Eq. (24). The implementation of this idea, however, is beyond the scope of the present work and constitutes a project for future research.

The final question to be discussed here is whether it is still possible to compute a reasonable estimation for the intrinsic contribution to the finite-size error of the interaction energy $\Delta viN$ without detailed knowledge of $\chi 0N(q,\omega )$. This is indeed possible, as we do have access to finite-size corrected data for *S*(*q*) in the form of the black diamonds in Fig. 14 that have been obtained from the *static approximation* using as the input $GNFSC(q)\u2248GTDL(q)$. From these, we can deduce an approximate expression for the finite-size effects in *S*^{N}(*q*) as

In particular, we stress that Eq. (25) cannot be exact, as $SstatFSC(q)$ entails the (small) systematic error from the *static approximation* [i.e., neglecting the frequency dependence of the full local field correction *G*(*q*, *ω*)], which is absent from *S*^{N}(*q*). Still, we stress that this systematic bias is small, especially at the high densities that are of interest to the present work,^{70,71} and that Eq. (25) is only used to estimate the small correction to the interaction energy per particle and not *V*^{N}/*N* itself. The corresponding expression for the intrinsic contribution to Δ$v$^{N} is then given by

Let us conclude our investigation of finite-size effects in electronic structure simulations under extreme conditions with an analysis of the interaction energy shown in Fig. 15. In particular, we only show a magnified segment around the finite-size corrected datasets, as the overall trend is quite similar to the case of *θ* = 2 shown in Fig. 1 above. Furthermore, the zero-order expression $\Delta vd,0N$ from Eq. (21) only becomes appropriate for large *N*, and the full estimation of $\Delta vdN$ is needed. Still, there remains a significant residual error that only vanishes within the given Monte Carlo error bars for *N* ≳ 66. Using Eq. (26) to estimate the intrinsic error in $v$ gives the purple squares, where the residual error is reduced by nearly an order of magnitude. We thus conclude that accurate knowledge of the static response functions *χ*^{N}(*q*) and $\chi 0N(q)$ can be used *in lieu* of the full frequency dependence of $\chi 0N(q,\omega )$ to compute a reasonable estimate of $\Delta viN$ within the *static approximation* at low temperatures.

## IV. SUMMARY AND OUTLOOK

### A. Summary

In summary, we have presented a FSC scheme for the static structure factor *S*(*q*) and, in this way, also the interaction energy per particle $v$. More specifically, we have employed the density response formalism and identified an effectively frequency-averaged static local field correction $G\xafinvert(q)$ as a short-ranged exchange–correlation property that can be extracted from a QMC simulation of only a few particles.

As a practical application, we have investigated the UEG under extreme conditions^{27} and demonstrated that often as few as *N* = 4 electrons are sufficient to obtain the interaction energy with a relative accuracy of ∼0.2%. We stress that our approach is completely non-empirical and does not require any external input or simulation data for different *N*.

In addition, we have analyzed the applicability of our scheme upon increasing the density parameter *r*_{s} (thus the coupling strength) and found that it eventually breaks down when correlation effects dominate. Finally, we have investigated the UEG at a low temperature, where momentum shell effects constitute an additional source of finite-size errors. Still, these errors, too, are present in the density response function of the non-interacting system $\chi 0N(q)$, and we have outlined two different strategies to mitigate the *N*-dependence in this case: (i) knowledge of the static interacting response function *χ*^{N}(*q*) allows us to obtain an accurate (though not exact) FSC for both *S*(*q*) and $v$, which reduces the residual error in $v$ after the elimination of the discretization error by an additional order of magnitude, and (ii) using $\chi 0N(q,\omega )$ would allow constructing a consistent *N*-dependent implementation of the *static approximation*, which has the potential to fully remove momentum shell effects and, in this way, remove the need for an additional twist-averaging procedure.

### B. Outlook

We believe that our findings open new avenues for additional topics of research to be pursued in future investigations. First and foremost, we call to mind the direct relation between *S*(*q*) and the pair distribution function *g*(*r*), which can be exploited to derive a FSC for the latter. In particular, this could be vital for the determination of *g*(0), which has most recently been investigated by Hunger *et al.*^{96} and is important for different applications (e.g., Refs. 74 and 75). In addition, the idea behind *S*(*q*) can be straightforwardly extended to the imaginary-time density–density correlation function *F*(*q*, *τ*) defined in Eq. (4) above. This function is of high value, as it allows for the computation of the dynamic structure factor *S*(*q*, *ω*) and related properties,^{70,71,73} which, in turn, are accessible in experiments, e.g., using the x-ray Thomson scattering technique.^{97,98} Furthermore, we mention the well-known thermodynamic relations between different energies, e.g., Eq. (13), which potentially allow for generalizing our correction for $v$ to other quantities such as the kinetic energy. A particularly interesting extension of our work is the generalization to multi-component systems such as electrons in the Coulomb potential of positively charged nuclei. It is well known that the concept of the local field correction can be used in this case as well,^{99} and we expect that the general concept of decomposing the system into a short-range exchange–correlation part estimated in a QMC simulation and a long-range part known from theory will be applicable. Finally, we point out that our method is particularly successful for high temperatures and densities, which indicates that it may apply to classical molecular dynamics simulations^{100–102} as well.

## ACKNOWLEDGMENTS

We are grateful to S. Groth for sending us various CPIMC datasets for *S*^{N}(*q*) and $\chi 0N(q)$. This work was partly funded by the Center for Advanced Systems Understanding (CASUS), which is financed by German Federal Ministry of Education and Research (BMBF), and the Saxon Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament. We gratefully acknowledge the CPU time at the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN) under Grant No. shp00026 and on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.