Warm dense matter is one of the most active frontiers in plasma physics due to its relevance for dense astrophysical objects and for novel laboratory experiments in which matter is being strongly compressed, e.g., by high-power lasers. Its description is theoretically very challenging as it contains correlated quantum electrons at finite temperature—a system that cannot be accurately modeled by standard analytical or ground state approaches. Recently, several breakthroughs have been achieved in the field of fermionic quantum Monte Carlo simulations. First, it was shown that exact simulations of a finite model system ($30\u2026100$ electrons) are possible, which avoid any simplifying approximations such as fixed nodes [Schoof *et al.*, Phys. Rev. Lett. **115**, 130402 (2015)]. Second, a novel way to accurately extrapolate these results to the thermodynamic limit was reported by Dornheim *et al.* [Phys. Rev. Lett. **117**, 156403 (2016)]. As a result, now thermodynamic results for the warm dense electron gas are available, which have an unprecedented accuracy on the order of 0.1%. Here, we present an overview on these results and discuss limitations and future directions.

## I. INTRODUCTION

The uniform electron gas (UEG) (i.e., electrons in a uniform positive background) is inarguably one of the most fundamental systems in condensed matter physics and quantum chemistry.^{1} Most notably, the availability of accurate quantum Monte Carlo (QMC) data for its ground state properties^{2,3} has been pivotal for the success of the Kohn-Sham density functional theory (DFT).^{4,5}

Over the past few years, interest in the study of matter under extreme conditions has grown rapidly. Experiments with not only inertial confinement fusion targets^{6–8} and laser-excited solids^{9} but also astrophysical applications such as planet cores and white dwarf atmospheres^{10,11} require a fundamental understanding of the warm dense matter (WDM) regime, a problem now at the forefront of plasma physics and materials science. In this peculiar state of matter, both the dimensionless Wigner-Seitz radius $rs=r\xaf/a0$ (with the mean interparticle distance $r\xaf$ and Bohr radius *a*_{0}) and the reduced temperature $\theta =kBT/EF$ ($EF$ being the Fermi energy) are of order unity, implying a complicated interplay of quantum degeneracy, coupling effects, and thermal excitations. Because of the importance of thermal excitation, the usual ground-state version of DFT does not provide an appropriate description of WDM. An explicitly thermodynamic generalization of DFT^{12} has long been known in principle but requires an accurate parametrization of the exchange-correlation free energy (*f _{xc}*) of the UEG over the entire warm dense regime as an input.

^{13–17}

This seemingly manageable task turns out to be a major obstacle. The absence of a small parameter prevents a low-temperature or perturbation expansion, and consequently, Green function techniques in the Montroll-Ward and *e*^{4} approximations^{18,19} break down. Further, the linear response theory within the random phase approximation^{20,21} (RPA) and also with the additional inclusion of static local field corrections as suggested by, e.g., Singwi, Tosi, Land, and Sjölander^{22–24} (STLS) and Vashista and Singwi^{24,25} (VS), is not reliable. Quantum classical mappings^{26,27} are exact in some known limiting cases but constitute an uncontrolled approximation in the WDM regime.

The difficulty of constructing a quantitatively accurate theory of WDM leaves thermodynamic QMC simulations as the only available tool to accomplish the task at hand. However, the widely used path integral Monte Carlo^{28} (PIMC) approach is severely hampered by the notorious fermion sign problem^{29,30} (FSP), which limits simulations to high temperatures and low densities and precludes applications to the warm dense regime. In their pioneering work, Brown *et al.*^{31} circumvented the FSP by using the fixed-node approximation^{32} (an approach hereafter referred to as restricted PIMC, RPIMC), which allowed them to present the first comprehensive results for the UEG over a wide temperature range for $rs\u22651$.

Although these data have subsequently been used to construct the parametrization of *f _{xc}* required for thermodynamic DFT (see Refs. 24, 33, and 34), their quality has been called into question. First, RPIMC constitutes an uncontrolled approximation,

^{35–38}which means that the accuracy of the results for the finite model system studied by Brown

*et al.*

^{31}was unclear. This unsatisfactory situation has sparked remarkable recent progress in the field of fermionic QMC.

^{39–50}In particular, a combination of two complementary QMC approaches

^{51,52}has recently been used to simulate the warm dense UEG without nodal restrictions,

^{42}revealing that the nodal constraints in RPIMC cause errors exceeding 10%. Second, Brown

*et al.*

^{31}extrapolated their QMC results for

*N*= 33 spin-polarized (

*N*= 66 unpolarized) electrons to the macroscopic limit by applying a finite-

*T*generalization of the simple first-order finite-size correction (FSC) introduced by Chiesa

*et al.*

^{53}for the ground state. As we have recently shown,

^{47}this is only appropriate for low temperature and strong coupling and, thus, introduces a second source of the systematic error.

In this paper, we give a concise overview of the current state of the art of quantum Monte Carlo simulations of the warm dense electron gas and present new results regarding the extrapolation to the thermodynamic limit (TDL). Further, we discuss the remaining open questions and challenges in the field.

After a brief introduction to the UEG model (II), we introduce various QMC techniques, starting with the standard path integral Monte Carlo approach (A) and a discussion of the origin of the FSP (B). The sign problem can be alleviated using the permutation blocking PIMC (PB-PIMC, C) method, the configuration PIMC algorithm (CPIMC, D), or the density matrix QMC (DMQMC, E) approach. In combination, these tools make it possible to obtain accurate results for a finite model system over almost the entire warm dense regime (IV). The second key issue is the extrapolation from the finite to the infinite system, i.e., the development of an appropriate finite-size correction,^{47,53–57} which is discussed in detail in Sec. V. Finally, we compare our QMC results for the infinite UEG to other data (2) and finish with some concluding remarks and a summary of open questions.

## II. THE UNIFORM ELECTRON GAS

### A. Coordinate representation of the Hamiltonian

Following Refs. 44 and 54, we express the Hamiltonian (using Hartree atomic units) for $N=N\u2191+N\u2193$ unpolarized electrons in coordinate space as

with the well-known Madelung constant $\xi M$ and the periodic Ewald pair interaction

Here, $R=n1L$ and $G=n2/L$ denote the real and reciprocal space lattice vectors, respectively, with $n1$ and $n2$ three-component vectors of integers, *L* the box length, $\Omega =L3$ the box volume, and *κ* the usual Ewald parameter.

### B. Hamiltonian in second quantization

In second quantized notation using a basis of spin-orbitals of plane waves, $\u3008r\sigma \u2009|ki\sigma i\u3009=1L3/2eiki\xb7r\delta \sigma ,\sigma i$, with $ki=2\pi Lmi,\u2009mi\u2208\mathbb{Z}3$ and $\sigma i\u2208{\u2191,\u2193}$, the Hamiltonian, Eq. (1), becomes

The antisymmetrized two-electron integrals take the form $wijkl\u2212=wijkl\u2212wijlk$, where

and the Kronecker deltas ensure both momentum and spin conservation. The first (second) term in the Hamiltonian, Eq. (3), describes the kinetic (interaction) energy. The operator $a\u0302i\u2020$ ($a\u0302i$) creates (annihilates) a particle in the spin-orbital $|ki\sigma i\u3009$.

## III. QUANTUM MONTE CARLO

### A. Path integral Monte Carlo

Let us consider *N* spinless distinguishable particles in the canonical ensemble, with the volume Ω, the inverse temperature $\beta =1/kBT$, and the density $N/\Omega $ being fixed. The partition function in coordinate representation is given by

where $R={r1,\u2026,rN}$ contains all 3*N* particle coordinates, and the Hamiltonian $H\u0302=K\u0302+V\u0302$ is given by the sum of a kinetic and a potential part, respectively. Since the low-temperature matrix elements of the density operator, $\rho \u0302=e\u2212\beta H\u0302$, are not readily known, we exploit the group property $e\u2212\beta H\u0302=(e\u2212\u03f5H\u0302)P$, with $\u03f5=\beta /P$ and positive integers *P*. Inserting *P* – 1 unities of the form $1\u0302=\u222bdR\alpha \u2009|R\u3009\alpha \u3008R|\alpha $ into Eq. (5) leads to

We stress that Eq. (6) is still exact and constitutes an integral over *P* sets of particle coordinates ($dX=dR0\u2026dRP\u22121$), the integrand being a product of *P* density matrices, each at *P* times the original temperature *T*. Despite the significantly increased dimensionality of the integral, this recasting is advantageous as the high temperature matrix elements can easily be approximated, most simply with the primitive approximation, $e\u2212\u03f5H\u0302\u2248e\u2212\u03f5K\u0302e\u2212\u03f5V\u0302$, which becomes exact for $P\u2192\u221e$. In a nutshell, the basic idea of the path integral Monte Carlo method^{28} is to map the quantum system onto a classical ensemble of ring polymers.^{58} The resulting high dimensional integral is evaluated using the Metropolis algorithm,^{59} which allows one to sample the 3*PN*-dimensional configurations **X** of the ring polymer according to the corresponding configuration weight $W(X)$.

### B. The fermion sign problem

To simulate *N* spin-polarized fermions, the partition function from the previous Section III A has to be extended to include a sum over all $N!$ permutations of particles:

where $\pi \u0302s$ denotes the exchange operator corresponding to the element *s* from the permutation group *S _{N}*. Evidently, Eq. (7) constitutes a sum over both positive and negative terms, so that the configuration weight function $W(X)$ can no longer be interpreted as a probability distribution. To allow fermionic expectation values to be computed using the Metropolis Monte Carlo method, we introduce the modified partition function

and compute fermionic observables as

with averages taken over the modified probability distribution $W\u2032(X)=|W(X)|$ and $S=W(X)/|W(X)|$ denoting the sign. The average sign, i.e., the denominator in Eq. (9), is a measure for the cancellation of positive and negative contributions and exponentially decreases with inverse temperature and system size, $\u3008S\u3009\u2032\u221de\u2212\beta N(f\u2212f\u2032)$, with *f* and $f\u2032$ being the free energy per particle of the original and the modified system, respectively. The statistical error of the Monte Carlo average value $\Delta O$ is inversely proportional to $\u3008S\u3009\u2032$

The exponential increase in the statistical error with *β* and *N* evident in Eq. (10) can only be compensated by increasing the number of Monte Carlo samples, but the slow $1/NMC$ convergence soon makes this approach unfeasible. This is the notorious fermion sign problem,^{29,30} which renders standard PIMC unfeasible even for the simulation of small systems at moderate temperature.

### C. Permutation blocking path integral Monte Carlo

To alleviate the difficulties associated with the FSP, Dornheim *et al.*^{43,44,48} recently introduced a novel simulation scheme that significantly extends fermionic PIMC simulations towards lower temperature and higher density. This so-called permutation blocking PIMC (PB-PIMC) approach combines (i) the use of antisymmetrized density matrix elements, i.e., determinants;^{60–62} (ii) a fourth-order factorization scheme to obtain accurate approximate density matrices for relatively low temperatures (large imaginary-time steps);^{63–66} and (iii) an efficient Metropolis Monte Carlo sampling scheme based on the temporary construction of artificial trajectories.^{43}

where the $W\u0302$ operators denote a modified potential term, which combines the usual potential energy $V\u0302$ with double commutator terms of the form

and, thus, requires the evaluation of all forces in the system. Furthermore, for each high-temperature factor, there appear three imaginary time steps. The final result for the partition function is given by

where the determinants incorporate the three diffusion matrices for each of the *P* factors^{44}

The key problem of fermionic PIMC simulations is the sum over permutations, where each configuration can have a positive or a negative sign. By introducing determinants, we analytically combine both positive and negative contributions into a single configuration weight (hence the label “permutation blocking”). Therefore, parts of the cancellation are carried out beforehand, and the average sign of our simulations [Eq. (9)] is significantly increased. Since this effect diminishes with increasing *P*, we employ the fourth-order factorization, Eq. (11), to obtain sufficient (although limited,^{44} $|\Delta V|/V\u2009\u2272\u20090.1%$) accuracy with only a small number of high-temperature factors. PB-PIMC is a substantial improvement over regular PIMC, but the determinants can still be negative, which means that the FSP is not removed by the PB-PIMC approach. To illustrate this point, in Fig. 1, we show simulation results for the average sign (here denoted as *S*) as a function of the density parameter *r _{s}* for a UEG simulation cell containing

*N*= 33 spin-polarized electrons subject to periodic boundary conditions. The red, blue, and black curves correspond to PB-PIMC results for three isotherms and exhibit a qualitatively similar behavior. At high

*r*, fermionic exchange is suppressed by the strong Coulomb repulsion, which means that almost all configuration weights are positive and

_{s}*S*is large. With increasing density, the system becomes more ideal and the electron wave functions overlap, an effect that manifests itself in an increased number of negative determinants. Nevertheless, the value of

*S*remains significantly larger than zero, which means that, for the three depicted temperatures, PB-PIMC simulations are possible over the entire density range. In contrast, the green curve shows the density-dependent average sign for standard PIMC simulations

^{31}at

*θ*= 1 and exhibits a significantly steeper decrease with density, limiting simulations to $rs\u22654$.

### D. Configuration path integral Monte Carlo

For CPIMC,^{40,41} instead of performing the trace over the density operator in the coordinate representation [see Eq. (5)], we trace over Slater determinants of the form

where, in the case of the uniform electron gas, *n _{i}* denotes the fermionic occupation number ($ni\u2208{0,1}$) of the

*i*-th plane wave spin-orbital $|ki\sigma i\u3009$. To obtain an expression for the partition function suitable for Metropolis Monte Carlo, we split the Hamiltonian into diagonal and off-diagonal parts, $H\u0302=D\u0302+Y\u0302$ (with respect to the chosen plane wave basis, see Sec. II), and explore a perturbation expansion of the density operator with respect to $Y\u0302$

with $Y\u0302(\tau )=e\tau D\u0302Y\u0302e\u2212\tau D\u0302$. In this representation, the partition function becomes

The matrix elements of the diagonal and off-diagonal operators are given by the Slater-Condon rules

where the multi-index $si=(pqrs)$ defines the four orbitals in which ${n(i)}$ and ${n(i\u22121)}$ differ, and we note that *p* < *q* and *r* < *s*. As in standard PIMC, each contribution to the partition function (17) can be interpreted as a $\beta \u2212$periodic path in imaginary time, but the path is now in Fock space instead of coordinate space. Evidently, the weight corresponding to any given path (second line of Eq. (17)) can be positive or negative. Therefore, to apply the Metropolis algorithm, we have to proceed as explained in Sec. III B and use the modulus of the weight function as our probability density. In consequence, the CPIMC method is also afflicted with the FSP. However, as it turns out, the severity of the FSP as a function of the density parameter is complementary to that of standard PIMC, so that weakly interacting systems, which are the most challenging for PIMC, are easily tackled using CPIMC. For a detailed derivation of the CPIMC partition function and the Monte Carlo steps are required to sample it see, e.g., Refs. 40–42, and 51.

### E. Density matrix quantum Monte Carlo

Instead of sampling contributions to the partition function, as in path integral methods, DMQMC samples the (unnormalized) thermal density matrix directly by expanding it in a discrete basis of outer products of Slater determinants

where $\rho {n},{n\u2032}=\u3008{n}|e\u2212\beta H\u0302|{n\u2032}\u3009$. The density-matrix coefficients $\rho {n},{n\u2032}$ appearing in Eq. (21) are found by simulating the evolution of the Bloch equation

which may be finite-differenced as

Following Booth and coworkers,^{67} we note that Eq. (23) can be interpreted as a rate equation and can be solved by evolving a set of positive and negative walkers, which stochastically undergo birth and death processes that, on average, reproduce the full solution. The rules governing the evolution of the walkers, as derived from Eq. (23), can be found elsewhere.^{45,67} The form of $\rho \u0302$ is known exactly at infinite temperature (*β* = 0, $\rho \u0302=1\u0302$), providing an initial condition for Eq. (22). For the electron gas, however, it turns out that simulating a differential equation that evolves a mean-field density matrix at inverse temperature *β* to the exact density matrix at inverse temperature *β* is much more efficient than solving Eq. (22), an insight that leads to the “interaction picture” version of DMQMC^{39,46} used throughout this work.

The sign problem manifests itself in DMQMC as an exponential growth in the number of walkers required for the sampled density matrix to emerge from the statistical noise.^{67–70} Working in a discrete Hilbert space helps to reduce the noise by ensuring a more efficient cancellation of positive and negative contributions, enabling larger systems and lower temperatures to be treated than would otherwise be possible. Nevertheless, at some point, the walker numbers required become overwhelming and approximations are needed. Recently, we have applied the initiator approximation^{71–73} to DMQMC ($i$ − DMQMC). In principle, at least, this allows a systematic approach to the exact result with an increasing walker number. More details on the use of the initiator approximation in DMQMC and its limitations can be found in Ref. 39.

### F. Applicability of the QMC methods

To conclude the discussion of Quantum Monte Carlo, in Fig. 2, we give a schematic overview of the parameter combinations where the different methods can be used to obtain results in the thermodynamic limit (for a discussion of finite-size corrections, see Sec. V) with a relative accuracy of $\Delta V/V\u223c0.003$. Standard PIMC (black) is only useful for high temperatures and low densities where fermionic exchange does not play an important role and, therefore, does not give access to the WDM regime. PB-PIMC (green) significantly extends the possible parameter combinations to lower temperature (down to $\theta =0.5$ for $rs\u22651$) and is available over the entire density range for $\theta \u22732$. In contrast, both CPIMC (red) and DMQMC (blue) are feasible for all *θ* at small *r _{s}* and eventually break down with increasing

*r*due to coupling effects. Despite their apparent similar range of applicability, it turns out that CPIMC is significantly more efficient at higher temperature, while DMQMC is superior at low

_{s}*θ*.

## IV. SIMULATION RESULTS FOR THE FINITE SYSTEM

The first step towards obtaining QMC results for the warm dense electron gas in the thermodynamic limit is to carry out accurate simulations of a finite model system. In Fig. 3, we compare results for the density dependence of the exchange correlation energy *E _{xc}* of the UEG for

*N*= 33 spin-polarized electrons and two different temperatures. The first results, shown as blue squares, were obtained with RPIMC

^{31}for $rs\u22651$. Subsequently, Groth, Dornheim, and co-workers

^{44,51}showed that the combination of PB-PIMC (red crosses) and CPIMC (red circles) allows for an accurate description of this system for $\theta \u22650.5$. In addition, it was revealed that RPIMC is afflicted with a systematic nodal error for densities greater than the relatively low value at which

*r*= 6. Nevertheless, the FSP precludes the use of PB-PIMC at lower temperatures and, even at $\theta =0.5$ and

_{s}*r*= 2, the statistical uncertainty becomes large. The range of applicability of DMQMC is similar to that of CPIMC, and the DMQMC results (green diamonds) fully confirm the CPIMC results.

_{s}^{39,46}Further, the introduction of the initiator approximation (i-DMQMC) has made it possible to obtain results up to

*r*= 2 for $\theta =0.5$. Although i-DMQMC is, in principle, systematically improvable and controlled, the results suggest that the initiator approximation may introduce a small systematic shift at lower densities.

_{s}In summary, the recent progress in fermionic QMC methods has resulted in a consensus regarding the finite-*N* UEG for temperatures $\theta \u22650.5$. However, there remains a gap at $rs\u22482\u22126$ and $\theta <0.5$ where, at the moment, no reliable data are available.

## V. FINITE SIZE CORRECTIONS

In this section, we describe in detail the finite-size correction scheme introduced in Ref. 47 and subsequently present detailed results for two elucidating examples.

### A. Theory

As introduced above (see Eq. (1) in Sec. II A), the potential energy of the finite simulation cell is defined as the interaction energy of the *N* electrons with each other, the infinite periodic array of images, and the uniform positive background. To estimate the finite-size effects, it is more convenient to express the potential energy in *k*-space. For the finite simulation cell of *N* electrons, the expression obtained is a sum over the discrete reciprocal lattice vectors **G**

where $S(k)$ is the static structure factor. In the limit as the system size tends to infinity and $\xi M\u21920$, this yields the integral

The task at hand is to find an accurate estimate of the finite-size error from Eq. (26), which, when added to the QMC result for $VN/N$, gives the potential energy in the thermodynamic limit. As a first step, we note that the Madelung constant may be approximated by^{55}

an expression that becomes exact in the limit as $\u03f5\u21920$. The Madelung term thus cancels the minus unity contributions to both the sum and the integral in Eq. (27).

The remaining two possible sources of the finite-size error in Eq. (26) are (i) the substitution of the static structure factor of the infinite system *S*(*k*) by its finite-size equivalent $SN(k)$ and (ii) the approximation of the continuous integral by a discrete sum, resulting in a discretization error. As we will show in Sec. V B, $SN(k)$ exhibits a remarkably fast convergence with system size, which leaves explanation (ii). In particular, about a decade ago, Chiesa *et al.*^{53} suggested that the main contribution to Eq. (26) stems from the $G=0$ term that is completely missing from the discrete sum. To remedy this shortcoming, they made use of the random phase approximation (RPA) for the structure factor, which becomes exact in the limit $k\u21920$. The leading term in the expansion of $SRPA(k)$ around *k* = 0 is^{26}

with $\omega p=3/rs3$ being the plasma frequency. The finite-*T* generalization of the FSC introduced by Chiesa *et al.*, hereafter called the BCDC-FSC, is^{31,47}

Eq. (30) would be sufficient if (i) $S0RPA(k)$ were accurate for $k\u2009\u2272\u20092\pi /L$ and (ii) all contributions to Eq. (26) beyond the $G=0$ term were negligible. As is shown in Sec. V B, both conditions are strongly violated in parts of the warm dense regime.

To overcome the deficiencies of Eq. (30), we need a continuous model function $Smodel(k)$ to accurately estimate the discretization error from Eq. (27)

A natural choice would be to combine the QMC results for $k\u2265kmin$, which include all short-ranged correlations and exchange effects, with the STLS structure factor for smaller *k*, which is exact as $k\u21920$ and incorporates the long-ranged behavior that cannot be reproduced using QMC due to the limited size of the simulation cell. However, as we showed in Ref. 47, a simpler approach using $SSTLS(k)$ [or the full RPA structure factor $SRPA(k)$] for all *k* is sufficient to accurately estimate the discretization error.

### B. Results

#### 1. Particle number dependence

To illustrate the application of the different FSCs, Fig. 4 shows results for the unpolarized UEG at *θ* = 2 and *r _{s}* = 1. The green crosses in panel (b) correspond to the raw, uncorrected QMC results that, clearly, are not converged with system size

*N*. The raw data points appear to fall onto a straight line when plotted as a function of $1/N$. This agrees with the BCDC-FSC formula, Eq. (30), which also predicts a $1/N$ behavior, and suggests the use of a linear extrapolation (the green line). However, while the linear fit does indeed exhibit good agreement with the QMC results, the computed slope does not match Eq. (30). Further, the points that have been obtained by adding $\Delta VBCDC$ to the QMC results, i.e., the yellow asterisks, do not fall onto a horizontal line and do not agree with the prediction of the linear extrapolation (see the horizontal green line). To resolve this peculiar situation, we compute the improved finite-size correction [Eq. (31)] using both the static structure factor from STLS ($SSTLS$) and the combination of STLS with the QMC data ($Scomb$) as input. The resulting corrected potential energies are shown as black squares and red diamonds, respectively, and appear to exhibit almost no remaining dependence on system size. In panel (c), we show a segment of the corrected data, magnified in the vertical direction. Any residual finite-size errors [due to the QMC data for

*S*(

*k*) not being converged with respect to

*N*, see panel (d)] can hardly be resolved within the statistical uncertainty and are removed by an additional extrapolation. In particular, to compute the final result for

*V*/

*N*in the thermodynamic limit, we obtain a lower bound via a linear extrapolation of the corrected data (using $SSTLS$) and an upper bound by performing a horizontal fit to the last few points, all of which are converged to within the error bars. The dotted grey line in panel (b), which connects to the extrapolated result, shows clearly that the results of this procedure deviate from the results of a naive linear extrapolation.

Finally, in panel (d) of Fig. 4, we show results for the static structure factor *S*(*k*) for the same system. As explained in Sec. V A, momentum quantization limits the QMC results to discrete *k* values above a minimum value $kmin=2\pi /L$. Nevertheless, the *N* dependence of the *k* grid is the only apparent change of the QMC results for *S*(*k*) with system size, and no difference between the results for the three particle numbers studied can be resolved within the statistical uncertainty (see also the magnified segment in the inset). The STLS curve (red) is known to be exact in the limit $k\u21920$ and smoothly connects to the QMC data, although for larger *k* there appears an almost constant shift. The full RPA curve (grey) exhibits a similar behavior, albeit deviating more significantly at intermediate *k*. Finally, the RPA expansion around *k* = 0 [Eq. (29), light blue] only agrees with the STLS and full RPA curves at very small *k* and does not connect to the QMC data even for the largest system size simulated.

To further stress the importance of our improved finite-size correction scheme, Fig. 5 shows results again for *θ* = 2 but at higher density, $rs=0.1$. In this regime, the CPIMC approach (and also DMQMC) is clearly superior to PB-PIMC and simulations of *N* = 700 unpolarized electrons in *N _{b}* = 189 234 basis functions are feasible. Due to the high density, the finite-size errors are drastically increased compared to the previous case and exceed 50% for

*N*= 38 particles [see panels (a) and (b)]. Further, we note that the BCDC-FSC is completely inappropriate for the

*N*values considered, as the yellow asterisks are clearly not converged and differ even more strongly from the correct TDL than the raw uncorrected QMC data.

Our improved FSC, on the other hand, reduces the finite-size errors by two orders of magnitude (both with $SSTLS$ and $Scomb$) and approaches Eq. (30) only in the limit of very large systems [$N\u2273104$; see panel (a)]. The small residual error is again extrapolated, as shown in panel (c).

Finally, we show the corresponding static structure factors in panel (d). The RPA expansion is again insufficient to model the QMC data, while the full RPA and STLS curves smoothly connect to the latter.

### 2. Comparison to other methods

To conclude this section, we use our finite-size corrected QMC data for the unpolarized UEG to analyze the accuracy of various other methods that are commonly used. In Fig. 6(a), the potential energy per particle, *V*/*N*, is shown as a function of *r _{s}* for the isotherm with

*θ*= 2. Although all four depicted curves exhibit qualitatively similar behavior, there are significant deviations between them [see panel (b), where we show the relative deviations from a fit to the QMC data in the TDL]. Let us start with the QMC results: the black squares correspond to the uncorrected raw QMC data for

*N*= 66 particles (see Ref. 52) and the red diamonds to the finite-size corrected data from Ref. 47. As expected, the finite-size effects drastically increase with density from $|\Delta V|/V\u22481%$, at

*r*= 10, to $|\Delta V|/V\u226550%$, at $rs=0.1$. This again illustrates the paramount importance of accurate finite-size corrections for QMC simulations in the warm dense matter regime. The RPA calculation (green curve) is accurate at high density and weak coupling. However, with increasing

_{s}*r*, the accuracy quickly deteriorates and, already at moderate coupling,

_{s}*r*= 1, the systematic error is of the order of 10%. The yellow asterisks show the SLTS result, which agrees well with the simulations (the systematic error does not exceed 3%) over the entire

_{s}*r*-range considered, i.e., up to

_{s}*r*= 10. Finally, the blue curve has been obtained from the recent parametrization of

_{s}*f*by Karasiev

_{xc}*et al.*

^{34}(KSDT), for which RPIMC data have been used as an input. While there is a reasonable agreement with our new data for $rs\u22731$ (with $|\Delta V|/V\u223c2%$), there are significant deviations at smaller

*r*, which only vanish for $rs<10\u22124$.

_{s}## VI. SUMMARY AND OPEN QUESTIONS

Let us summarize the status of *ab initio* thermodynamic data for the uniform electron gas at finite temperature. The present paper has given an overview of recent progress in *ab initio* finite temperature QMC simulations that avoid any additional simplifications such as fixed nodes. While these simulations do not “solve” the fermion sign problem, they provide a reasonable and efficient way on how to *avoid it*, in many practically relevant situations, by combining simulations that use different representations of the quantum many-body state: the coordinate representation (direct PIMC and PB-PIMC) and Fock states (CPIMC and DMQMC). With this, it is now possible to obtain highly accurate results for up to $N\u223c100$ particles in the entire density range and for temperatures $\theta \u22730.5$. As a second step, we demonstrated that these comparatively small simulation sizes are sufficient to predict results for the macroscopic uniform electron gas *not significantly losing accuracy*.^{47} This unexpected result is a consequence of a new highly accurate finite-size correction that was derived by invoking STLS results for the static structure factor.

With this procedure, it is now possible to obtain thermodynamic data for the uniform electron gas with an accuracy on the order of 0.1%. Even though pure electron gas results cannot be directly compared to warm dense matter experiments, they are of high value to benchmark and improve additional theoretical approaches. Most importantly, this concerns finite-temperature versions of the density functional theory (such as orbital-free DFT), which is the standard tool to model realistic materials and which will benefit from our results for the exchange-correlation free energy. Furthermore, we have also presented a few comparisons with earlier models such as RPA, STLS, or the recent fit of Karasiev *et al.* (KSDT), the accuracy and errors of which can now be unambiguously quantified. We found that among the tested models, the STLS is the most accurate one. We wish to underline that even though exchange-correlation effects are often small compared to the kinetic energy, their accurate treatment is important to capture the properties of real materials, see e.g., Ref. 74.

In the following, we summarize the open questions and outline future research directions.

Construction of an improved fit for the exchange-correlation free energy due to their key relevance as input for finite-temperature DFT. Such fits are straightforwardly generated from the current results but require a substantial extension of the simulations to arbitrary spin polarization. This work is currently in progress.

The presently available accurate data are limited to temperatures above half the Fermi energy, as a consequence of the fermion sign problem. A major challenge will be to advance to lower temperatures, $\Theta <0.5$, and to reliably connect the results to the known ground state data. This requires substantial new developments in the area of the three quantum Monte Carlo methods presented in this paper (CPIMC, PB-PIMC and DMQMC) and new ideas on how to combine them. Another idea could be to derive simplified versions of these methods that treat the FSP more efficiently but still have acceptable accuracy.

The present

*ab initio*results allow for an entirely new view on previous theoretical models. For the first time, a clear judgment about the accuracy becomes possible, which more clearly maps out the sphere of applicability of the various approaches, e.g., Ref. 75. Moreover, the availability of our data will allow for improvements of many of these approaches via adjustment of the relevant parameters to the QMC data. This could yield, e.g., improved static structure factors, dielectric functions or local field correlations.Similarly, our data may also help to improve alternative quantum Monte Carlo concepts. In particular, this concerns the nodes for Restricted PIMC simulations, which can be tested against our data. This might help to extend the range of validity of those simulations to higher density and lower temperature. Since this latter method does not have a sign problem, it may allow to reach parameters that are not accessible otherwise.

A major challenge of Metropolis-based QMC simulations that are highly efficient for thermodynamic and static properties is to extend them to dynamic quantities. This can, in principle, be done via analytical continuation from imaginary to real times (or frequencies). However, this is known to be an ill-posed problem. Recently, there has been significant progress by invoking stochastic reconstruction methods or genetic algorithms. For example, for Bose systems, accurate results for the spectral function and the dynamics structure factor could be obtained, e.g., Ref. 76 and references therein, which encourage also for applications to the uniform electron gas, in the near future.

Finally, there are a large number of additional applications of the presented

*ab initio*simulations. This includes the 2D warm dense UEG where thermodynamic results of similar accuracy should be straightforwardly accessible. Moreover, for the electron gas, at high density, $rs\u22720.1$, relativistic corrections should be taken into account. Among the presented simulations, CPIMC is perfectly suited to tackle this task and to provide*ab initio*data also for correlated matter at extreme densities.

## ACKNOWLEDGMENTS

This work was supported by the Deutsche Forschungsgemeinschaft via Project No. BO1366-10 and via SFB TR-24 Project No. A9 and Grant No. shp00015 for CPU time at the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN). T.S. acknowledges the support of the U.S. DOE/NNSA under Contract No. DE-AC52-06NA25396. F.D.M. was funded by an Imperial College PhD Scholarship. F.D.M. and W.M.C.F. used computing facilities provided by the High Performance Computing Service of Imperial College London, by the Swiss National Supercomputing Centre (CSCS) under Project ID s523, and by ARCHER, the UK National Supercomputing Service, under EPSRC Grant No. EP/K038141/1 and via a RAP award. F.D.M. and W.M.C.F. acknowledge the research environment provided by the Thomas Young Centre under Grant No. TYC-101.