Thermodynamics of the uniform electron gas: Fermionic path integral Monte Carlo simulations in the restricted grand canonical ensemble

The uniform electron gas (UEG) is one of the key models for the understanding of warm dense matter—an exotic, highly compressed state of matter between solid and plasma phases. The difficulty in modelling the UEG arises from the need to simultaneously account for Coulomb correlations, quantum effects, and exchange effects, as well as finite temperature. The most accurate results so far were obtained from quantum Monte Carlo (QMC) simulations with a variety of representations. However, QMC for electrons is hampered by the fermion sign problem. Here, we present results from a novel fermionic propagator path integral Monte Carlo in the restricted grand canonical ensemble. The ab initio simulation results for the spin‐resolved pair distribution functions and static structure factor are reported for two isotherms θ=1,2 (T in the units of the Fermi temperature). Furthermore, we combine the results from the linear response theory in the Singwi‐Tosi‐Land‐Sjölander scheme with the QMC data to remove finite‐size errors in the interaction energy. We present a new corrected parametrization for the interaction energy v(rs,θ) and the exchange–correlation free energy fxc(rs,θ) in the thermodynamic limit, and benchmark our results against the restricted path integral Monte Carlo by Brown et al. [Phys. Rev. Lett. 110, 146405 (2013)] and configuration path integral Monte Carlo/permutation‐blocking path integral Monte Carlo by Dornheim et al. [Phys. Rev. Lett. 117, 115701 (2016)].

Characteristic of all these diverse systems is the governing role of electronic quantum effects, moderate to strong Coulomb correlations, and finite-temperature effects. Quantum effects of electrons are of relevance at low temperature and/or if matter is very highly compressed, such that the temperature is of the order of (or lower than) the Fermi temperature (for a recent overview, see ref. [22]). Due to the simultaneous relevance of these effects, a theoretical description of WDM is difficult. Among the actively used approaches are quantum kinetic theory, [23][24][25][26] quantum hydrodynamics, [27][28][29] and density functional theory (DFT) simulations. The latter, for the first time, enabled the self-consistent simulation of realistic WDM, that includes both plasma and condensed matter phases. [30][31][32] The most accurate results for WDM, so far, were obtained via first-principles computer simulations such as quantum Monte Carlo (QMC). [4,[33][34][35][36][37][38][39][40][41] The first QMC results for dense quantum plasmas and WDM were obtained by Ceperley and Militzer et al. [42,34] and Filinov et al. [35,43] However, the simulations of the electronic part of WDM are strongly hampered by the fermion sign problem (FSP), see ref. [44] for a recent overview, which makes simulations at strong electronic degeneracy impossible. One approach to relieve this problem is the fixed-node approximation, [42] which works well in the ground state. However, the associated systematic error at finite temperatures is largely unclear. On the other hand, in the direct fermionic QMC simulations of Filinov et al., the sign problem was reduced by optimized Monte Carlo updates [45] but the simulations remained restricted to moderate quantum degeneracy.
Recently, a number of breakthroughs in QMC simulations occurred. The first was the idea to "decompose" the complex WDM system into subsystems and start by developing improved simulations for the most challenging component-the partially degenerate electrons. This can be done by considering the model or the uniform electron gas (UEG or "jellium") that has proven highly successful in condensed matter physics and many-body theory for many decades, where the focus is on the ground-state properties and finite but low temperatures. [46] Brown et al. [47,48] presented the first restricted path integral Monte Carlo (RPIMC) simulations for the jellium model at finite temperature, 0.0625 ≤ ≤ 8, where = k B T∕E F (E F denotes the Fermi energy), for which ground-state QMC approaches are not applicable. Even though RPIMC simulations have no sign problem, they could not be extended into the degenerate regime, r s ≲ 1, where r s is the ratio of the mean interparticle distance to the Bohr radius, r s = r∕a B , and even for moderate r s , in the range between 1 and 5, the potential energy showed an unexpected non-monotonic behaviour. This behaviour was clarified by novel configuration PIMC (CPIMC) simulations-that is, PIMC formulated in Fock space-that were developed at Kiel University a decade ago. [37] Remarkably, in contrast to standard PIMC, CPIMC has no sign problem at strong degeneracy (r s ≲ 0.7). In combination with an improved configuration space approach-permutation blocking PIMC (PB-PIMC) [40,49,50] -and confirmed by independent density matrix QMC (DM-QMC), [51] it became clear that RPIMC in this parameter range exhibits large systematic errors. On the other hand, the combination of CPIMC and PB-PIMC made it possible to obtain thermodynamic results for the warm dense UEG, which are free of systematic errors, over the entire density range, for ≳ 0.5. These results were first obtained for finite particle numbers (N ≲ 100) and subsequently extrapolated to the thermodynamic limit by Dornheim et al. [52] These data were then connected to existing ground-state results [53] by Groth et al. [54] This led to a first-principles based accurate analytical parametrization of the exchange-correlation (XC) free energy of the UEG. [4,54] The comparisons with earlier parametrizations are discussed by Dornheim et al. [55] Despite the overall very good quality of the parametrization presented by Dornheim et al. [52] for some parameter combinations, inaccuracies were observed, in particular, when computing derivatives with respect to temperature or density. [56] For this reason, independent benchmarks are highly desirable.
In this article, we present our first results from a new fermionic PIMC simulations, which is called fermionic propagator PIMC (FP-PIMC). This new approach is formulated in the grand canonical ensemble (GCE) and has a number of important advantages. First, it allows for a direct access to a one-particle DM and, correspondingly, to a momentum distribution (not discussed in the present article). In contrast to a standard PIMC in the canonical ensemble (CE), it uses as input the chemical potential (for each particle species) and natively includes physically relevant fluctuations in the charge density and spin polarization. The isothermal compressibility can be studied directly as a function of temperature and density.
In present work, however, we use the so-called restricted GCE (R-GCE) where the particle number fluctuations are artificially suppressed to perform a more direct comparison with the CE. In spite of the different ensembles used, for quantities extrapolated to the thermodynamic limit this should not be relevant. A detailed analysis of two temperature isotherms, = 1, 2, is reported which covers the broad density range, 0.1 ≤ r s ≤ 10, studied previously [52] with the combination of two techniques-CPIMC and PB-PIMC. For each combination of parameters [ , r s , { }], the finite-size effects are analysed for the particle number in the range: 20 ≤ N ≤ 100.
First, we benchmark our method against the known analytical results for the ideal Fermi gas and then apply it to the interacting UEG. Furthermore, novel finite-size corrections (FSCs) are derived and used to extrapolate the potential energy and the XC free energy to the thermodynamic limit. Comparing our results with the parametrization by Dornheim et al., [52] we observe several deviations that are investigated in detail. Furthermore, we present results for the spin-resolved static structure factor (SSF) and pair distribution functions (PDFs) and compare to spin-resolved linear response results within the Singwi-Tosi-Land-Sjölander (STLS) scheme.
The article is organized as follows. Section 2 outlines the theoretical and simulation methods. Section 3 presents the simulation results, and the article is concluded by a discussion in section 4.

THEORETICAL METHODS
Here, we describe the theoretical methods that are used in this article to simulate finite-temperature jellium. The first is a new implementation of fermionic PIMC that is formulated in the GCE and will be called FP-PIMC. The second is linear response theory based on the spin-revolved STLS approximation that will be used for comparison and for the derivation of FSCs in order to extrapolate the FP-PIMC data to the thermodynamic limit.

Path integral representation of the DM
In the framework of the path integral formulation of quantum mechanics introduced by Feynman, the many-body DM is presented via an ensemble of the imaginary-time particle trajectories. In the following, we will consider Fermi statistics and treat the fermions with different spin projections as different species. The total radius vector for particles of the same type will be denoted as follows: where the upper index denotes the spin state s = {↑, ↓}, the first lower index counts a number of particles with the same spin (1 … N s ), and the second lower index denotes the imaginary time argument, p = p , with = ∕P and 0 ≤ p ≤ P. To further shorten the notations, we introduce the space-time variable, , which defines a system microstate-a specific microscopic configuration of particle trajectories. To account for Fermi statistics in the partition function Z F originally expressed via distinguishable particle trajectories, the sum over N s ! permutations is introduced. Each permutation contributes with a positive or negative sign, Sgn( s ) = ±1, depending on the permutation parity in each spin-subspace: This expression is a starting point for the standard PIMC scheme. [57] 2.1.2 Sign problem in fermionic QMC The PIMC method remains problematic for solving many-fermion systems due to the so-called FSP. Not all terms in the sum (1) are positive-definite. The standard Monte Carlo scheme performs the important sampling of each term in the series based on the absolute value, while in the evaluation of the partition function and thermodynamic observables a } should be weighed with the overall sign. As a result, the expectation value, A, is evaluated with the statistical error A, which scales inversely proportional to the average sign s The average sign can be expressed as the ratio of Fermi and Bose partition functions for a given system and decays exponentially with the particle number N, the inverse temperature , and the difference of the corresponding free energies, Δf . As an illustration of the inefficiency of this straightforward approach, in Figure 1, we consider the simulation of an ideal quantum gas. The PDFs are compared for three statistics, that is, the particles are treated as boltzmannons, bosons, and fermions (spin polarized). The temperature is set to = T∕T F = 1. The quantum exchange effects for bosons and fermions are clearly observed. In the bosonic case, the exchange leads to an effective attraction-the prerequisite for Bose condensation. For fermions, Pauli blocking leads to the formation of an exchange hole at the origin-the particles in the same spin state experience effective repulsion at close spatial separation. This physical effect is completely missing in the standard representation of the fermionic DM as the sign-alternating series with the (positive-definite) bosonic DM and distinguishable particle trajectories, c.f. Equation (1). The spatial distribution of trajectories for a Fermi gas will be sampled similar to the bosonic case with a tendency to collapse to a high density phase. To recover the exchange hole near the origin would require a strong cancellation between different terms in (1) from different permutation classes s ∈ ℘. The correct behaviour of PDF at small distances, in this case, could be only recovered by the attenuating sign: lim r→0 g Fermi (r) ≈ ⟨s⟩ ⋅ g Bose (r), with ⟨s⟩ ∼ 10 −5 . This would lead to large statistical uncertainties in the thermodynamic averages (2).

Reducing the FSP with determinant propagators
To recover the effect of XC hole already on the level of Monte Carlo sampling scheme, one needs an alternative representation of the fermionic partition function. One solution is to use the antisymmetric propagators (in the form of many-body Slater determinants) to sample propagation of particle trajectories in the imaginary time. This requires to write a modified partition function where the summation over different permutation classes is performed analytically in the kinetic energy part of the DM (or in the N-particle free propagator). In this representation, the N-particle DM contains the Slater determinants, M s (p − 1, p), between each successive imaginary times (or time slices) Now the fermionic DM explicitly contains a new effective action S A (p − 1, p) where the Slater determinants introduce the additional exchange interaction contribution W ex along with the standard potential energy term W The antisymmetric (fermionic) free-particle propagators (denoted in the following as "FP") between two adjacent time-slices are given by where ⃗ M s is the N s × N s diffusion matrix This representation has several advantages. First, the fermionic DM is antisymmetric with respect to any pair exchange of spin-like fermions. Second, the probability of microstates is now proportional to the value of Slater determinants and the degeneracy of the latter at small particle separations correctly recovers the Pauli blocking effect. Third, the fermionic propagators partially reduce the sign problem. The combined contribution of different permutations, ∑ ↑ ( ↓ ) (sum of the terms in the sign-alternating series), now reduces to the composite sign of the determinants M ↑(↓) , while their absolute value enters explicitly in the probability density.
The efficiency of similar representation of the DM in the PIMC has been demonstrated by several authors: Takahashi and Imada, [58] Filinov, [35] and Lyubartsev. [59] Chin [60] used the determinant propagators to simulate relatively large ensembles of electrons in 3D quantum dots. A similar approach termed permutation-blocking PIMC (PB-PIMC) was recently developed by Dornheim et al. [40] and applied to the UEG at WDM conditions. [61] The present FP-PIMC method in some aspects (discussed later) is similar to those carried out by Dorheim et al. and Chin [40,60] but is different in the implemented Monte Carlo algorithm formulated for the sampling of the fermionic DM in the GCE (a detailed description will be given in a separate paper). It is based on Equations (4)-(10) applied to the grand canonical partition function (e.g., two-component system) In contrast to the CE simulations cited earlier, we specify values of the chemical potential for each particle species, { 1 , 2 , … }. The particle number in the simulation volume fluctuates around the average values, ⟨N 1 ⟩, ⟨N 2 ⟩, and the fluctuations are self-consistently controlled by the changes in the DM once we add/remove a particle trajectory as a special Monte Carlo update. The main complication, compared to classical simulations in the GCE, arises due to the indistinguishability of particle trajectories once they form the Slater determinants. Hence, the add/remove update in a quantum system has a non-local effect due to the exchange interaction contribution (7) and leads to the change in the sign-prefactors, Sgn p , in Equation (4).
Next, in the CE, the thermodynamic observables (such as the internal energy and the pressure) are related to the diagonal elements of the DM or partial derivatives of the partition function Z. The grand canonical formulation can provide important additional information and natively includes new physical effects: • direct account of density, charge, or spin fluctuations in the thermodynamic observables; • information on the chemical potential, the temperature, and the density dependence, = (T, V, ⟨N⟩); • information on the particle number distribution, the average density, and the isothermal compressibility, In this article, we report on the first application of our novel approach. As a model system, we use the jellium model of the warm dense UEG [62] (electrons on a uniform neutralizing ionic background with Ewald summation applied) and T . Different curves correspond to different temperatures (in 10 3 K) and different numbers of high-temperature propagators P used in Equation (14): , and P12(P = 12)]. Thermodynamic averages delivered by quantum Monte Carlo with a statistical error ≤ 1% remain feasible for ⟨S⟩ ≥ 10 −3 benchmark our results against the most accurate finite-temperature data by Dornheim et al. [52] that were obtained in the CE. To minimize possible discrepancies due to different ensembles, here we implement the so-called restricted-GCE (R-GCE) by adding to the original DM the slope Gaussian function, which suppresses the particle fluctuations around the mean value, with the typical choice i = 2. In the following, we sketch the main ideas shared by the present approach and those by Dorheim et al. and Chin [40,60] to explain why fermionic PIMC simulations with the determinant propagators become feasible and efficient. The effect of exchange interaction captured by the FP depends on the number of factorization factors (P) in the DM. The largest effect is achieved once the effective wavelength, ( = ∕P), in each propagator (10) is comparable with the interparticle distance, and the particles "feel" their fermionic nature. The Pauli exclusion principle prevents "condensation" of particles to a high-density phase as observed for the bosons (see Figure 1). In the simulations, this results in the increase of both the average interparticle distance (for the spin-like electrons) and in the average sign, which in its turn leads to reduction of the statistical error (2).
As an example, in Figure 2, we demonstrate how the average sign ⟨S⟩ decays with the degeneracy parameter = n 3 T . The set of upper curves corresponds to different number of propagators P in Equation (14). As shown by the upper set of curves (P = 3 and T = 31 … 250 ⋅ 10 3 K), the sign scaling is given solely by the degeneracy factor, ⟨S⟩ ∝ e − (P)⋅ . There is an additional prefactor, (P), which increases with the number of propagators, making the sign problem more severe.
This behaviour motivates to use of more accurate short-time propagators in the DM. [40] The resulting larger time step in the high-temperature factorization extends the applicability range of fermionic simulations to a higher degeneracy, when the average sign stays above a small but acceptable value, for example, ⟨S⟩ ≥ 10 −3 . In the present work, we take advantage of the fourth-order factorization scheme proposed by Chin and Chen [63] and Sakkos et al., [64] with the choice = ∕P, (2t

Ideal Fermi gas: Equation of state
To validate the present approach and demonstrate its efficiency, we have performed extensive tests against the known analytical results for an ideal Fermi gas. The equation of state (EoS), the electron density, and the chemical potential are mutually related via the Fermi-Dirac integrals as follows As a result, the reduced pressure (scaled in the units of the classical pressure) presented in Figure 3 can be expressed as The ratio of both integrals can be expressed via the Padé approximation in the form As a final result, to present the analytical form for the compressibility factor (17), one needs as input only the electron density expressed via the r s parameter, n = (4∕3 r 3 s ) −1 , and the thermal wavelength (T). In Figure 3, the set of dashed lines corresponds to the analytical result (18) for the three isotherms along with the FP-PIMC data for the pressure (expressed via the kinetic energy). It can be observed that the simulation results quite well agree with the expected analytical behaviour for a broad range of densities. The characteristic statistical errors in the FP-PIMC simulations are shown for a single Monte Carlo run: 1 day−single thread CPU.
We again refer to Figure 1 and remark that this type of simulations using the direct fermionic PIMC, based on the distinguishable "bosonic" trajectories (to account the sum over permutations), would not be possible, except for either high temperatures, T∕T F > 4, corresponding to low degeneracy or strongly non-ideal gas, when the Coulomb repulsion effectively substitutes the effect of Pauli blocking. In contrast, the use of the FP-PIMC with the anti-symmetrized DM allows the investigation of the electron gas at moderate degeneracy when the deviations from the classical pressure can reach 20-30%.

STLS theory for the electron gas resolving the two spin components
To analyse spatial correlations in the UEG with two spin components, we compare our simulations to the linear response formalism. Within this framework, to probe the spin-resolved correlations, one studies the response of the electron density to an external weak perturbation (in the form of a weak space and time-dependent electric or magnetic field). As a result, the finite-temperature spin-density linear response function can be expressed as [65] (we consider a uniform and isotropic system, q = |⃗ q|), where 0 (q, ; T, ) is the free electron gas polarizability. The value of the chemical potential at a given temperature T is defined by the normalization condition of the noninteracting system: with the dimensionless kinetic energy, t = k ∕E F . The dynamic local field correction (LFC), G(q, ), accounts for the XC effects beyond the random-phase approximation (RPA). Different approximation methods have been developed to evaluate the LFC. In the following, we consider the static approximation, G(q) = G(q, = 0), proposed by Pines, [66] which is the key quantity in the STLS scheme. [67] Then the dielectric function theory will be validated by the PIMC simulations. By assuming that the two-particle distribution function can be expressed in terms of the single-particle distribution function, the static LFC can expressed as where the SSF S(k) is expressed via the density response function using the fluctuation-dissipation theorem For the numerical implementation, it is more convenient, to switch from the -integration to a sum over the Matsubara frequencies [67] S(q, with l = 2 lk B T∕ℏ. In the STLS scheme, one obtains the self-consistent solution of two coupled equations (21) and (23). Once the SSF is known, the interaction energy in the system is evaluated as In a multicomponent system, one generalizes the above mentioned equations by considering the effective potential (⃗ q, ), which acts on an -species particle produced by the density fluctuation n (⃗ q, ) in the subsystem of the -type particles where all microscopic correlation effects are included in the LFC. Next, the density perturbations in the linear response approximation are introduced according to where V ext (q, ) represents the space-time Fourier transform for an external field. Finally, the density-density response functions, (k, ), between the and the particle species are defined by the relation By comparison of Equations (25)- (27) and omitting the arguments, one can write the final result for a two-component system as with v k = 4 ∕k 2 . This system of equations is solved within the STLS scheme introduced earlier and involves the self-consistent solution for S 11 , S 22 , S 12 (21) and G 11 , G 22 , G 12(21) , where the lower indices denote the spin components ′ ( , ′ =↑, ↓).

Spin-resolved SSF: comparison of QMC and STLS results
The spin-resolved SSF obtained from the STLS theory presented in section 2.2 is plotted in Figures 4 and 5 (see the dashed lines) for r s = 1, 2, 4, 8 and temperatures = 1, 2, along with the RPIMC results by Brown et al. [47] and the present FP-PIMC data. In all plots, the wave vector is scaled with the density parameter: qa B ⋅ r s . First, we notice a very good agreement between the two sets of QMC data. The discretized wave vectors correspond to the system of N = 66 unpolarized electrons in both cases, for = 2 (r s = 2 − 4), and, for ⟨N⟩ = 20(50) (FP-PIMC data), for = 1 (r s = 2 − 4). The error bars show the characteristic statistical error in the FP-PIMC simulations for a single Monte Carlo run (1 day−single thread CPU), while in the final results for the energies and the SSF (presented later), the average over 10-40 Monte Carlo runs was used (if necessary) to reduce the statistical uncertainty. In general, we can confirm a nearly perfect agreement for the static correlations for the spin-unlike electrons (S ↑↓ ), with some deviations observed for the spin-like electrons (S ↑↑ ) at high densities ( = 1, r s = 1, 2) and the smallest wave vector q 0 = 2 ∕L. Here, we observe that due to finite density fluctuations in the R-GCE, the finite-size effects in the SSF are, typically, larger than that in CE for the long wavelengths (in finite size systems at q 0 ). This effect can be explicitly observed in Figures 4 and 5 by the comparison the SSF for ⟨N⟩ = 20 and ⟨N⟩ = 66: only the first point, S(q 0 ), is noticeably shifted upwards for the smaller system size. Although in the present comparison the FP-PIMC results look almost identical with the previous RPIMC data, the deviations due to the approximate treatment of spin statistics by RPIMC will become well pronounced in the spin-resolved PDFs presented in section 3.2. This should lead to a slight shift in the S(q) values within some finite range of momenta, which seems to be difficult to resolve on the scale of Figures 4 and 5.
Now we turn to comparison with the STLS theory. In general, the agreement is quite satisfactory. The deviations are reduced (increased) with an increase (decrease) of the wave vector. Also note that the deviations to the QMC data in the ↑↑ and ↑↓ cases are of the opposite sign: we plot S ↑↓ as −S 12 , and the QMC data for S 11 (S 12 ) compared to the STLS are slightly more positive (negative) for qa B ≲ 3. This is the main reason why in the total SSF (see Figure 6) the agreement with the STLS looks much better as the deviations nearly cancel. This in its turn also results in a better agreement in the interaction energy v STLS(QMC) expressed via the integral (37).
With the last note, we want to underline the importance of the STLS theory for the present analysis of the finite-size effects presented in section 3.3. Due to the absence of QMC data for small wave vectors (qa B r s ≲ 1), the use of STLS remains an excellent option to provide reliable predictions for the total SSF, S q = 0.5 ⋅ (S 11 + S 22 + 2S 12 ), in the long-wavelength limit and, thereby, include quite accurately both exchange and correlation effects. [52] From Figure 6, we can confirm quite a good match between the STLS (dashed curve) and the QMC data at the smallest momentum, q 0 = 2 ∕L. This allows us to construct a combined solution (the model SSF is presented by the solid line), S m q , which accurately interpolates between the ab initio QMC for q > q 0 and includes the exact screening sum rule [68] satisfied by the STLS theory in the q → 0 limit with the plasma frequency p = √ 3∕r 3 s . In practice, the part of the S m q curve that goes through the QMC data points was obtained as the Fourier transform of the radial PDF g(r). The later one was evaluated quite accurately as the average over 10-40 independent QMC runs  N), , kT}. The interpolation between the S STLS and S QMC was done as the weighted average in the momentum interval [q 0 , q 1 ], where q 1 is the next smallest discrete momentum value after q 0 Finally, the reconstructed model SSF was used to accurately evaluate the finite-size effects in the interaction energy via Equations (37) and (39).
By a close inspection of the deviations between the STLS result and the model solution, S m q , in Figure 6, we can notice the following general trend that holds for both temperatures, = 1, 2. At high density (r s = 1), the STLS underestimates the true SSF (based on the QMC data) for q > q 0 . Next, by lowering the density (increasing the r s value) the S m q start to reduce (increase) its value systematically for qa B r s ≲ 3.5 (qa B r s ≥ 3.5), and, occasionally, quite a good coincidence with the STLS is achieved around r s = 2. By further lowering the density (r s = 4 and 8), the relative position of two results mutually swap: now the STLS curve goes slightly above at least up to qa B r s ∼ 3.5, and then stays always below the S m q -curve and demonstrates a much slower convergence to 1.
This peculiar behaviour has a direct consequence in the observed density dependence of the energy isotherms presented in Figure 12 (see section 3.3). Compared to the QMC data (solid and dashed lines with symbols), the STLS prediction, v STLS , shows more negative values in the high-density regime (r s ≲ 2), that is, it overestimates the XC hole at the origin (see also the PDF in Figure 7); next, the STLS symbols get close to the QMC data around r s = 2-3 (exactly when S STLS and S m q nearly match) and, finally, for r s ≥ 4 the STLS data points go above the QMC results. This last trend can be explained by the analysis of the PDF for the spin-unlike electrons in Figure 8: the STLS theory significantly underestimates the Coulomb correlation hole for the electrons with different spins in the low-density phase, see panels with r s = 2, 4, 8.

Spin-resolved PDFs
In the next section, we perform a comparison of the spin-resolved PDF evaluated with three different techniques. First, the self-consistent solution of the STLS equations for the UEG, for different combinations of r s and T, allows us to reconstruct where the lower indices denote the spin projections , ′ =↑ (↓). These results can be directly compared with the RPIMC data of Brown et al. [47] and the PDF obtained from the R-GCE: In Figure 8, we plot the PDF for spin-unlike electrons for three temperature values, = 1, 2, and 4. Each panel corresponds to a chosen density parameter r s in the regime of moderate (r s = 1) to strong (r s = 8) coupling.
First, we note that the interplay between correlation and thermal effects in g ↑↓ , as predicted by our FP-PIMC simulations, well agrees with the previous RPIMC results with the free-particle nodes. Here, the quantum exchange effects do not play a dominant role, since for the electrons with different spin projections, the approximation in treating Fermi statistics by RPIMC is not relevant.
Let us now consider the comparison of the STLS results. As expected, for the moderately interacting UEG (r s = 1 and r s = 2), both QMC-results are found to be in a very good agreement with the STLS predictions. Both the value g(0) and the slope of the PDF (at least for r∕a B ≤ 1.5) are reproduced by the STLS quite accurately. This also validates the correct implementation and solution of the two-component STLS equations (see section 2.2), and the self-consistent iterative solution for the local-field correction factor G . For larger coupling (r s > 4), we observe systematic deviations between the STLS and the QMC results. At the same time, the Coulomb correlation hole develops at the origin and g(0) goes to zero. The r s dependence of g(0) has a direct relation with the asymptotic behaviour of the momentum distribution function, as it was recently demonstrated by Hunger et al. [69] To explore the role of exchange effects, as a next step, we analyse the PDF for spin-like electrons g ↑↑(↓↓) . Here, the XC hole (g(0) = 0) is present at all densities being the combined result of the Pauli exclusion principle and Coulomb repulsion. By comparing both g ↑↑ and g ↑↓ , we observe that the XC hole is initially driven by quantum statistics (for r s ≲ 4) and, for larger r s , by Coulomb correlations (the slopes of g ↑↑ and g ↑↓ become similar for r s = 8).
Here, for all temperature-density combinations, the three approaches deviate from each other. For small particle separations, the STLS theory always predicts negative values of g(r). This is a well-known drawback of STLS-for strongly coupled systems, it becomes less accurate. In this case, the dynamical correlations going beyond the static LFC approximation (used in G (q, T)) become important. Here we refer to recent ab initio QMC simulation results for the dynamic LFC of the UEG. [70] Moreover, the static LFC in a neural net representation has been obtained by Dornheim et al. [71]

Interaction energy: extrapolation of QMC results to the thermodynamic limit
The interaction energy (per particle) in the TDL follows from the model SSF On the other hand, the interaction energy in a finite N-particle system (independent from QMC simulations) can be estimated as a sum over the discrete momenta G using the same model SSF: where M is the Madelung constant. The corresponding curves (dark green) are presented in Figures 9-11. The FSCs to the QMC data, v N , are estimated as the difference of the interaction energy per particle in the infinite system and the N-particle system with the imposed periodic boundary conditions. [52] which are shown by the brown curves (with star) in Figures 9-11. We compare this to the FSCs within the RPA proposed by Chiesa et al. [72] Δ RPA = p 4N coth shown by the dashed-blue curves (with open circle) in these figures.
In Figure 12, we compare our results for the potential energy in the thermodynamic limit (with FSC applied) over a broad range of densities, 0.1 ≤ r s ≤ 10, with earlier QMC data: the RPIMC results by Brown et al. [48] and CPIMC/PB-PIMC results by Dornheim et al. [52] Two isotherms are shown for = 1 and = 2. Figure 12c shows the relative deviations (in percent) from the present FP-PIMC results.

F I G U R E 9
Finite-size scaling of the potential energy per particle for = 1 and r s = {1, 2, 4, 8}: "v N "-the uncorrected (raw) quantum Monte Carlo data; "Δ N " with the removed discretization error, Equation (39), using the reconstructed (in the thermodynamic limit) static structure factor S m (k); "Δ RPA "-the finite-size corrections using the random-phase approximation, Equation (40); "v N (s k )"-the expected potential energy (per particle) with the finite-size effects obtained using S m (k) and summation over the reciprocal lattice vectors for given simulation cell size (Equation 38). The dotted lines show the result of extrapolation to the thermodynamic limit. The red crosses show the results obtained by Dornheim et al. [52] for the same set of parameters (the value in the thermodynamic limit (TDL) [the cross at the origin] and v N=66 ) F I G U R E 10 Same as in Figure 9, but for higher densities, r s = {0.1, 0.3, 0.5, 0.7}. The red cross at the origin shows the results obtained by Dornheim et al. [52] for the same set of parameters. In the high-density regime (r s = 0.1 … 0.3), the random-phase approximation-finite-size correction (open blue circles) converges much slower with the particle number compared with Δ N (Equation (39)) F I G U R E 11 Same as in Figure 9, but for = 2 and r s = {1,2,4,6}. The red cross at the origin shows the results of ref. [52] for the same set of parameters F I G U R E 12 Potential energy of the macroscopic uniform electron gas for two isotherms, = 1 and = 2. (a, b) Comparison of three sets of quantum Monte Carlo data and the Singwi-Tosi-Land-Sjölander (STLS) theory: (1) Brown et al., [48] (2) Dornheim et al., [52]  As a general trend, we observe that the isotherms from FP-PIMC (Label 3) and CPIMC/PB-PIMC (Label 2) converge to each other at high densities, r s ≤ 1, see Figure 12a. In particular, for = 2 and r s < 2, the relative deviations do not exceed 0.25%. However, for lower densities, the deviations become quite visible and reach a local maximum at the intermediate densities, that is, at r s ≈ 4 ( = 2) and r s ≈ 2 ( = 1). In contrast, the deviations to the RPIMC data are systematically larger for all densities and temperatures and can reach up to 1.7%.
A possible origin of the discrepancies observed between FP-PIMC and CPIMC/PB-PIMC at r s = 2 … 4 is not clear at the moment. Both isotherms (in the TDL) have been obtained by extrapolation of the finite-size results using the TA B L E 1 Potential energy per particle (in Hartree units) obtained using several extrapolation schemes STLS-PIMC combined with SSF, S m q , and the extrapolation procedure seems to be not the source of the error. The discrepancies are already present on the level of the finite-size QMC data (similar to RPIMC to PBPIMC comparison) and can be seen in Figures 9-11: compare the energies for N = 66 electrons denoted for CPIMC/PB-PIMC by the red cross. The slope of scaling with the particle number N is mainly given by v N (s k )-the interaction energy estimated from Equation (38) (see, e.g., the green curves with symbols in Figure 11), which in all cases lies slightly above (and more close to) the CPIMC/PB-PIMC energies. This estimate involves the angular-averaged value of the total SSF, S m , and, therefore, miss some information on the angular-resolved correlation effects. They are expected to be less relevant at high densities when the particle trajectories strongly overlap: see Figure 10 the cases with r s = 0.1 … 0.7), where we, indeed, find a good agreement between Equation (38) and the raw QMC data points. However, by lowering the density (r s ≥ 2), the QMC curve systematically shifts to more negative values; thus, the expectation from Equation (38) becomes less accurate. This trend F I G U R E 13 (a) Finite-size scaling of the potential energy in grand canonical ensemble (GCE) simulations (labelled "gce") versus the restricted GCE (R-GCE) defined by the particle number-variance parameter = 2,4,8. Similar to Figure 9 extrapolation to the TDL is demonstrated (for "gce" and = 2 case). Simulation parameters: r s = 8 and = 1. In all cases, the chemical potential is kept fixed ( = −0.077). The red crosses at the origin and 1∕N = 66 −1 show the results obtained by Dornheim et al. [52] (b,c) For ⟨N⟩ = 66 electrons: (1) the particle number distribution in the GCE versus the R-GCE for each spin component P(N ↑(↓) ) (indistinguishable on the non-log scale), and (2) the particle number imbalance distribution function, P(ΔN), with Δ N = (N ↑ − N ↓ ) characterizing the spin fluctuations in the ⟨N⟩ = 66 system as an example. In a macroscopic uniform electron gas, the spin susceptibility retains a constant value and leads to ⟨ΔN 2 ⟩ ∝ ⟨N⟩kT, characterizing the scaling of fluctuations in a finite-size system. In the case of the R-GCE, the fluctuations are additionally suppressed by the slope function, see Equation (13) is also confirmed for the PB-PIMC data (the N = 66 point); however, the effect is less pronounced than that for the present FP-PIMC.
Another source of discrepancy lies in the use of different thermodynamic ensembles. In the CE, all microstates used in the statistical average correspond to the unpolarized electron gas (the number of electrons is fixed by N ↑ = N ↓ ), whereas in the present R-GCE the spin polarization only on average remains zero (⟨N ↑ ⟩ = ⟨N ↓ ⟩). It is allowed that a large fraction of microstates has a non-zero spin polarization. This results in an effectively higher electron density in one of the spin components and should lead to more negative values of the interaction energy, v(r s , ), as can be seen from the density dependence presented in Table 1. A more detailed comparison between different thermodynamic ensembles (canonical and [restricted] grand canonical) is left for future analysis. As a first preliminary result in Figure 13 (first panel), we compare the FSC of the potential energy in the grand canonical simulation with the R-GCE ensemble (with = 2). Both simulations are performed for the same value of the chemical potential. The use of the slope function in Equation (13) slightly changes the average density (or the r s -value) from r s = 7.983 (GCE) to r s = 7.994 (R-GCE). This fact can explain that we have obtained a slightly different value, v(r s , ), in the TDL being more negative for the GCE case due to a higher density. Still as expected, both simulations provide nearly the same value in the TDL, while the energy in a finite-size system can be substantially different in both ensembles: compare the uncorrected QMC data v N for ⟨N⟩ ≲ 50.
For a complete comparison, the CE results by Dornheim et al. [52] are also indicated in Figure 13a by the red crosses. They are shifted upwards, but the value of the shift, as shown in the analysis, cannot be explained by the choice of a different ensemble (however, for any finite N, the results can, indeed, be different due to the fluctuations in the spin polarization). Thus, the origin of observed discrepancy with the present results remains unresolved.
In Figure 13b,c, we provide additional information on the particle number distribution and the relative fluctuations of both spin components.

XC free energy
Most of DFT calculations are performed with temperature-independent XC functionals developed for the ground state. [73][74][75][76] However, the thermal effects may play an important role in WDM conditions. [77] In particular, it has been TA B L E 2 Fit parameters used for the parametrization of the interaction energy, according to Equation (42) recently shown that the correct account of XC thermal effects significantly improves agreement with recent experimental measurements for transport and optical properties of deuterium. [78] The effect of finite-temperature XC on the EoS of hydrogen in the electron liquid regime has been studied, [79] and the significant corrections to the EoS were found in the WDM regime (T > 10 4 K).
Even these few examples demonstrate that the EC free energy is the essential quantity for finite-temperature DFT simulations, and its accurate extraction from high-precision numerical simulations is highly desirable. This problem gets complicated, as all QMC simulations are limited to finite-size systems and an accurate exclusion of the finite-size effects in the thermodynamic observables remains one of the crucial and important task.
In order to evaluate the XC part of the free energy (per particle) f xc , we use its definition as an integral over the average Coulomb interaction energy, v(r s , ) = V int (r s , )∕⟨N⟩, and the coupling constant r s , that is, Applying the partial derivative with respect to r s on both sides, we get the relation v(r s , ) = 2f xc (r s , ) + r s f xc (r s , ) r s .
As a next step, f xc is presented via the Padé approximation proposed by Karasiev et al. [80] F xc ∕N = f xc (r s , Substituting this parametrization in Equation (42), we obtain the interaction energy v(r s , ) as a function of density and a set of free fit parameters {a, b, c, d, e}( ), which are used to fit the finite-size corrected FP-PIMC data in section 3.3. Their optimal values for temperatures = 1 and = 2 are summarized in Table 2, along with the other methods, and represent analytical parametrizations for the XC free energy, f xc , of the UEG.
The corresponding comparison between three sets of QMC data and predictions by the STLS is demonstrated in Figure 14. The deviations with respect to current results, Δf xc ∕f xc , reach about 1.5-2% for the STLS and the RPIMC, [47,80] and are within 0.7% for the CPIMC/PBPIMC data.
The source of these deviations lies in the values of the interaction energy, v(r s , ), used as an input for the parametrization in Equation (42) and has been properly discussed in the earlier section.

DISCUSSION AND CONCLUSIONS
In this article, we report on accurate thermodynamic data for the unpolarized UEG. We used a novel simulation approach based on the PIMC representation and the antisymmetric free-particle propagators in the form of Slater determinants. The use of the later allows to significantly reduce the fermion sign problem and recover the Pauli blocking effect in the spatial (coordinate)-space representation of the N-particle DM as discussed in section 2.1.3. Our ab inito results cover a broad range of electron densities, 0.1 ≤ r s ≤ 10, and two temperature isotherms T = T F and 2T F . Both thermal and degeneracy effects are quite important at these parameters and should be accurately treated. To cross validate our results, we performed comparison with an ideal Fermi gas and the previous ab initio finite-temperature QMC simulations-the restricted PIMC by Brown et al. [47] and CPIMC/PBPIMC by Dornheim et al. [52] The comparison includes the spin-resolved SSF and the PDFs, the interaction energy, v(r s , ), and the XC free energy, f xc (r s , ), extrapolated to the TDL. For both quantities, we provide the fit analysis. Quite a good agreement is demonstrated with the CPIMC/PBPIMC at high densities (r s ≤ 1). This validates that both methods correctly treat the XC contribution to the thermodynamic properties at high degeneracy. In the low-density regime (r s > 2), however, we found the systematic deviations with the CPIMC/PBPIMC that can reach up to 0.7% in both Δv∕v and Δf xc ∕f xc . Similar deviations with the RPIMC data are systematically larger by factor 2-3. By the analysis of the PDF in section 3.2, we can validate that the use of the fixed-node approximation leads to underestimation of the XC hole in g ↑↑ (r) for the spin-like electrons, while we observe a nearly perfect agreement in g ↑↓ (r) where the spin statistics does not play a role.
Also, we present the detailed analysis of the spin-resolved correlations within the two-component (in spin) linear response method based on the finite-temperature STLS approximation. We characterize the accuracy of the STLS for different parameter ranges by analysing the systematic deviations from the present QMC data. The most important part of this analysis is the reconstruction of the total SSF with the long-wavelength behaviour provided by the STLS. The short-range correlation effects are included via the SSF from the QMC. For small wave vectors, the asymptotic behaviour reproduces the RPA. As a next step, the obtained SSF allows us to remove the finite-size errors in the potential energy related with the discretization of momentum and perform the extrapolation to the TDL.
In summary, by performing extensive simulations, we present accurate ab initio data for f xc (r s , ) and v(r s , ). A consistent and reliable fit of f xc is required for the construction of XC functionals in finite-temperature DFT, as well as in experimental applications in the WDM regime. For the parametrization of f xc , we apply the analytical form proposed by Karasiev et al., [80] which provides the accurate fit to the QMC data. Our data are a valuable benchmark of the previous simulation schemes. In particular, the non-negligible deviations with respect to the parametrization of f xc presented by Dornheim et al. [52] are demonstrated and need to be clarified in future analysis.
As a next prospective application of the developed approach, we view simulations for a broader temperature range, 0.5 ≤ ≤ 8, and in GCE to analyse the XC contributions in the chemical potential and pressure. Furthermore, a direct generalization to the two-component systems, such as hydrogen and helium plasmas, seems to be feasible.