Ab initio Quantum Monte Carlo simulation of the warm dense electron gas

Ab initio Quantum Monte Carlo simulation of the warm dense electron gas

Tobias Dornheim dornheim@theo-physik.uni-kiel.de Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, D-24098 Kiel, Germany    Simon Groth [ Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, D-24098 Kiel, Germany    Fionn D. Malone Department of Physics, Imperial College London, Exhibition Road, London SW7 2AZ, United Kingdom    Tim Schoof Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, D-24098 Kiel, Germany    Travis Sjostrom Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA    W.M.C. Foulkes Department of Physics, Imperial College London, Exhibition Road, London SW7 2AZ, United Kingdom    Michael Bonitz Institut für Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel, D-24098 Kiel, Germany
July 15, 2019

Warm dense matter is one of the most active frontiers in plasma physics due to its relevance for dense astrophysical objects as well as 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 ( electrons) is possible that 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 that have an unprecedented accuracy on the order of . Here we present an overview on these results and discuss limitations and future directions.

Warm dense matter, Quantum Monte Carlo, uniform electron gas.
Valid PACS appear here
preprint: AIP/123-QED

]TD and SG contributed equally to this work.

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 chemistryloos (). Most notably, the availability of accurate quantum Monte Carlo (QMC) data for its ground state propertiesalder (); perdew () has been pivotal for the success of Kohn-Sham density functional theory (DFT) ks (); dft_review ().

Over the past few years, interest in the study of matter under extreme conditions has grown rapidly. Experiments with inertial confinement fusion targets nora (); schmit (); hurricane3 () and laser-excited solids ernst (), but also astrophysical applications such as planet cores and white dwarf atmospheres knudson (); militzer (), 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 (with the mean interparticle distance and Bohr radius ) and the reduced temperature ( 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 DFTmermin () has long been known in principle, but requires an accurate parametrization of the exchange-correlation free energy () of the UEG over the entire warm dense regime as an inputkarasiev2 (); dharma (); gga (); burke (); burke2 ().

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 approximationskremp_springer (); vorberger_pre04 () break down. Further, linear response theory within the random phase approximationgupta (); pdw84 () (RPA) and also with the additional inclusion of static local field corrections as suggested by, e.g., Singwi, Tosi, Land, and Sjölanderstls_original (); stls (); stls2 () (STLS) and Vashista and Singwivs_original (); stls2 () (VS), are not reliable. Quantum classical mappingsrpa0 (); dwp () 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 Carlocep () (PIMC) approach is severely hampered by the notorious fermion sign problemloh (); troyer () (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.brown () circumvented the FSP by using the fixed-node approximationnode () (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 .

Although these data have subsequently been used to construct the parametrization of required for thermodynamic DFT (see Refs. brown3, ; karasiev, ; stls2, ), their quality has been called into question. Firstly, RPIMC constitutes an uncontrolled approximationvfil1 (); vfil2 (); hyd1 (); hyd2 (), which means that the accuracy of the results for the finite model system studied by Brown et al.brown () was unclear. This unsatisfactory situation has sparked remarkable recent progress in the field of fermionic QMCvfil4 (); blunt (); malone (); malone2 (); dornheim (); dornheim2 (); dornheim_prl (); dornheim_cpp (); tim_prl (); tim1 (); tim2 (); blunt2 (). In particular, a combination of two complementary QMC approachesdornheim3 (); groth () has recently been used to simulate the warm dense UEG without nodal restrictionstim_prl (), revealing that the nodal contraints in RPIMC cause errors exceeding . Secondly, Brown et al.brown () extrapolated their QMC results for spin-polarized ( unpolarized) electrons to the macroscopic limit by applying a finite- generalization of the simple first-order finite-size correction (FSC) introduced by Chiesa et al.chiesa () for the ground state. As we have recently showndornheim_prl (), this is only appropriate for low temperature and strong coupling and, thus, introduces a second source of 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. 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 (III.1) and a discussion of the origin of the FSP (III.2). The sign problem can be alleviated using either the permutation blocking PIMC (PB-PIMC, III.3) method, or the configuration PIMC algorithm (CPIMC, III.4), or the density matrix QMC (DMQMC, III.5) 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 correctionchiesa (); lin (); drummond (); fraser (); krakauer (); dornheim_prl (), which is discussed in detail in Sec. V. Finally, we compare our QMC results for the infinite UEG to other data (V.2.2) and finish with some concluding remarks and a summary of open questions.

Ii The Uniform Electron Gas

ii.1 Coordinate representation of the Hamiltonian

Following Refs. fraser, ; dornheim2, , we express the Hamiltonian (using Hartree atomic units) for unpolarized electrons in coordinate space as


with the well-known Madelung constant and the periodic Ewald pair interaction


Here and denote the real and reciprocal space lattice vectors, respectively, with and three-component vectors of integers, the box length, the box volume, and the usual Ewald parameter.

ii.2 Hamiltonian in second quantization

In second quantized notation using a basis of spin-orbitals of plane waves, , with , and , the Hamiltonian, Eq. (1), becomes


The antisymmetrized two-electron integrals take the form , 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 () creates (annihilates) a particle in the spin-orbital .

Iii Quantum Monte Carlo

iii.1 Path integral Monte Carlo

Let us consider spinless distinguishable particles in the canonical ensemble, with the volume , the inverse temperature , and the density being fixed. The partition function in coordinate representation is given by


where contains all particle coordinates, and the Hamiltonian is given by the sum of a kinetic and a potential part, respectively. Since the low-temperature matrix elements of the density operator, , are not readily known, we exploit the group property , with and positive integers . Inserting unities of the form into Eq. (5) leads to


We stress that Eq. (6) is still exact and constitutes an integral over sets of particle coordinates (), the integrand being a product of density matrices, each at times the original temperature . 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, , which becomes exact for . In a nutshell, the basic idea of the path integral Monte Carlo method cep () is to map the quantum system onto a classical ensemble of ring polymers chandler (). The resulting high dimensional integral is evaluated using the Metropolis algorithm metropolis (), which allows one to sample the -dimensional configurations of the ring polymer according to the corresponding configuration weight .

iii.2 The fermion sign problem

To simulate spin-polarized fermions, the partition function from the previous section has to be extended to include a sum over all permutations of particles:


where denotes the exchange operator corresponding to the element from the permutation group . Evidently, Eq. (7) constitutes a sum over both positive and negative terms, so tht the configuration weight function 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 and 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, , with and being the free energy per particle of the original and the modified system, respectively. The statistical error of the Monte Carlo average value is inversely proportional to ,


The exponential increase of the statistical error with and evident in Eq. (10) can only be compensated by increasing the number of Monte Carlo samples, but the slow convergence soon makes this approach unfeasible. This is the notorious fermion sign problem loh (); troyer (), which renders standard PIMC unfeasible even for the simulation of small systems at moderate temperature.

iii.3 Permutation blocking path integral Monte Carlo

To alleviate the difficulties associated with the FSP, Dornheim et al.dornheim (); dornheim2 (); dornheim_cpp () 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 det1 (); det2 (); det3 (); (ii) a fourth-order factorization scheme to obtain accurate approximate density matrices for relatively low temperatures (large imaginary-time steps) ho1 (); ho2 (); ho3 (); ho4 (); and (iii) an efficient Metropolis Monte Carlo sampling scheme based on the temporary construction of artificial trajectories dornheim ().

In particular, we use the factorization introduced in Refs. ho3, ; ho2,


where the operators denote a modified potential term, which combines the usual potential energy 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 factorsdornheim2 (),


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 , we employ the fourth-order factorization, Eq. (11), to obtain sufficient (although limited dornheim2 (), ) 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.

Figure 1: Density dependence of the average sign of a PB-PIMC simulation of the uniform electron gas. Also shown are standard PIMC data taken from Ref. brown, . The figure has been taken from Dornheim et al. dornheim2 ().

To illustrate this point, in Fig. 1 we show simulation results for the average sign (here denoted as ) as a function of the density parameter for a UEG simulation cell containing 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 , fermionic exchange is suppressed by the strong Coulomb repulsion, which means that almost all configuration weights are positive and 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 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 simulationsbrown () at and exhibits a significantly steeper decrease with density, limiting simulations to .

iii.4 Configuration path integral Monte Carlo

For CPIMCtim1 (); tim2 (), 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 case of the uniform electron gas, denotes the fermionic occupation number () of the -th plane wave spin-orbital . To obtain an expression for the partition function suitable for Metropolis Monte Carlo, we split the Hamiltonian into diagonal and off-diagonal parts, (with respect to the chosen plane wave basis, see Sec. II), and explore a perturbation expansion of the density operator with respect to :


with . 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 defines the four orbitals in which and differ and we note that and . As in standard PIMC, each contribution to the partition function (17) can be interpreted as a 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.2 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 required to sample it see, e.g., Refs. tim1, ; tim2, ; groth, ; tim_prl, .

iii.5 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 . The density-matrix coefficients appearing in Eq. (21) are found by simulating the evolution of the Bloch equation,


which may be finite-differenced as


The matrix elements of the Hamiltonian are as given as in Eqs. (18) and (19).

Following Booth and coworkersbooth (), 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 elsewherebooth (); blunt (). The form of is known exactly at infinite temperature (, ), 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 DMQMCmalone2 (); malone () 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 noisebooth (); spencer (); kolodrubetz (); shepherd3 (). 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 approximationcleland (); shepherd (); shepherd2 () to DMQMC (DMQMC). In principle, at least, this allows a systematic approach to the exact result with increasing walker number. More details on the use of the initiator approximation in DMQMC and its limitations can be found in Ref. malone2, .

Iv Simulation results for the finite system

Figure 2: Exchange-correlation energy of spin-polarized electrons as a function of the density parameter for two isotherms. Shown are results from CPIMC and PB-PIMC taken from Ref. groth, , restricted PIMC from Ref. brown, , and DMQMC from Ref. malone2, . For , all data has been shifted by Hartree. In the case of DMQMC, the initiator approximation is used.

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. 2, we compare results for the density dependence of the exchange correlation energy of the UEG for spin-polarized electrons and two different temperatures. The first results, shown as blue squares, were obtained with RPIMC brown () for . Subsequently, Groth, Dornheim and co-workersdornheim2 (); groth () showed that the combination of PB-PIMC (red crosses) and CPIMC (red circles) allows for an accurate description of this system for . In addition, it was revealed that RPIMC is afflicted with a systematic nodal error for densities greater than the relatively low value at which . Nevertheless, the FSP precludes the use of PB-PIMC at lower temperatures and, even at and , 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 resultsmalone (); malone2 (). Further, the introduction of the initiator approximation (i-DMQMC) has made it possible to obtain results up to for . Although i-DMQMC is, in principle, systematically improvable and controlled, the results suggest that the initiator approximation may introduce a small systematic shift at higher densities.

In summary, the recent progress in fermionic QMC methods has resulted in a consensus regarding the finite- UEG for temperatures . However, there remains a gap at and 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. dornheim_prl, and subsequently present detailed results for two elucidating examples.

v.1 Theory

As introduced above (see Eq. (1) in Sec. II.1), the potential energy of the finite simulation cell is defined as the interaction energy of the 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 -space. For the finite simulation cell of electrons, the expression obtained is a sum over the discrete reciprocal lattice vectors :


where is the static structure factor. In the limit as the system size tends to infinity and , this yields the integral


Combining Eqs. (24) and (25) yields the finite-size error for a given QMC simulation:


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 , gives the potential energy in the thermodynamic limit. As a first step, we note that the Madelung constant may be approximated bydrummond ()


an expression that becomes exact in the limit as . 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 by its finite-size equivalent ; and (ii) the approximation of the continuous integral by a discrete sum, resulting in a discretization error. As we will show in Sec. V.2, exhibits a remarkably fast convergence with system size, which leaves explanation (ii). In particular, about a decade ago, Chiesa et al. chiesa () suggested that the main contribution to Eq. (26) stems from the 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 . The leading term in the expansion of around isrpa0 ()


with being the plasma frequency. The finite- generalization of the Chiesa et al. FSC, hereafter called the BCDC-FSC, is brown (); dornheim_prl ():


Eq. (30) would be sufficient if (i) were accurate for , and (ii) all contributions to Eq. (26) beyond the term were negligible. As is shown in Sec. V.2, 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 to accurately estimate the discretization error from Eq. (27):


A natural choice would be to combine the QMC results for , which include all short-ranged correlations and exchange effects, with the STLS structure factor for smaller , which is exact as 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. dornheim_prl, , a simpler approach using [or the full RPA structure factor ] for all is sufficient to accurately estimate the discretization error.

v.2 Results

v.2.1 Particle number dependence

Figure 3: Finite-size correction for the UEG at and : a) dependence of the FSCs; b) potential energy per particle, ; the dotted grey line corresponds to the TDL value where the had been substracted; c) extrapolation of the residual finite-size error; and d) corresponding static structure factors from QMC (for ), STLS, RPA, and the RPA expansion around , Eq. (29).

To illustrate the application of the different FSCs, Fig. 3 shows results for the unpolarized UEG at and . The green crosses in panel b) correspond to the raw, uncorrected QMC results that, clearly, are not converged with system size . The raw data points appear to fall onto a straight line when plotted as a function of . This agrees with the BCDC-FSC formula, Eq. (30), which also predicts a 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 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 () and the combination of STLS with the QMC data () 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 not being converged with respect to , 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 in the thermodynamic limit, we obtain a lower bound via a linear extrapolation of the corrected data (using ) 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. 3, we show results for the static structure factor for the same system. As explained in Sec. V.1, momentum quantization limits the QMC results to discrete values above a minimum value . Nevertheless, the dependence of the grid is the only apparent change of the QMC results for 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 and smoothly connects to the QMC data, although for larger there appears an almost constant shift. The full RPA curve (grey) exhibits a similar behavior, albeit deviating more significantly at intermediate . Finally, the RPA expansion around [Eq. (29), light blue] only agrees with the STLS and full RPA curves at very small and does not connect to the QMC data even for the largest system size simulated.

Figure 4: Finite-size correction for the UEG at and : a) dependence of the FSCs; b) potential energy per particle, ; c) extrapolation of the residual finite-size error; and d) corresponding static structure factors from QMC (for ), STLS, RPA, and the RPA expansion around , Eq. (29).

To further stress the importance of our improved finite-size correction scheme, Fig. 4 shows results again for but at higher density, . In this regime, the CPIMC approach (and also DMQMC) is clearly superior to PB-PIMC and simulations of unpolarized electrons in basis functions are feasible. Due to the high density, the finite-size errors are drastically increased compared to the previous case and exceed for particles [see panels a) and b)]. Further, we note that the BCDC-FSC is completely inappropriate for the 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 and ) and approaches Eq. (30) only in the limit of very large systems [; see panel a)]. The small residual error is again extrapolated, as shown in panel c).

Finally, we show the corresponding static 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.

v.2.2 Comparison to other methods

Figure 5: Potential energy per particle of the uniform electron gas at –simulations versus analytical models. Squares: QMC results for particles dornheim3 (), (red) rhombs: finite-size corrected QMC data (TDL) dornheim_prl (), green (yellow) curves: RPA (STLS) data stls2 (), blue: results of the parametrization of Ref. karasiev () (KSDT). Lower Fig.: relative deviations of all curves from the fit to the thermodynamic QMC results.

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. 5 a), the potential energy per particle, , is shown as a function of for the isotherm with . 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 particles (see Ref. dornheim3, ) and the red diamonds to the finite-size corrected data from Ref. dornheim_prl, . As expected, the finite-size effects drastically increase with density from , at , to , at . 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 the accuracy quickly deteriorates and, already at moderate coupling, , the systematic error is of the order of . The yellow asterisks show the SLTS result which agrees well with the simulations (the systematic error does not exceed ) over the entire -range considered, i.e., up to . Finally, the blue curve has been obtained from the recent parametrization of by Karasiev et al.karasiev () (KSDT), for which RPIMC data have been used as an input. While there is a reasonable agreement with our new data for (with ), there are significant deviations at smaller , which only vanish for .

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 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, DMQMC). With this it is now possible to obtain highly accurate results for up to particles in the entire density range and for temperatures . 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 loosing accuracy dornheim_prl (). 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 . 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 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, Vahista-Singwi, 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. perdew_13 ()

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

  1. 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.

  2. 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, 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 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.

  3. The present ab initio results allow for an entirely new view on previous theoretical models. For the first time, a clear judgement about the accuracy becomes possible which more clearly maps out the sphere of applicability of the various approaches, e.g. 4ebeling (). 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.

  4. 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.

  5. A major challenge of Metropolois-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. filinov_pra_12 () and references therein which are encouraging also for applications to the uniform electron gas, in the near future.

  6. Finally, there is 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, , 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.


This work was supported by the Deutsche Forschungsgemeinschaft via project BO1366-10 and via SFB TR-24 project A9 as well as grant shp00015 for CPU time at the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN). TS acknowledges the support of the US DOE/NNSA under Contract No. DE-AC52-06NA25396. FDM is funded by an Imperial College PhD Scholarship. FDM and WMCF 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 EP/K038141/1 and via a RAP award. FDM and WMCF acknowledge the research environment provided by the Thomas Young Centre under Grant No. TYC-101.


  • (1) P.-F. Loos and P.M.W. Gill, The uniform electron gas, WIREs Comput. Mol. Sci. 6, 410-429 (2016)
  • (2) J.P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23, 5048 (1981)
  • (3) D.M. Ceperley and B.J. Alder, Ground State of the Electron Gas by a Stochastic Method, Phys. Rev. Lett. 45, 566 (1980)
  • (4) W. Kohn, and L.J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys.Rev. 140, A1144 (1965)
  • (5) R.O. Jones, Density functional theory: Its origins, rise to prominence, and future, Rev. Mod. Phys. 87, 897-923 (2015)
  • (6) R. Nora et al., Gigabar Spherical Shock Generation on the OMEGA Laser Phys. Rev. Lett. 114, 045001 (2015)
  • (7) P.F. Schmit et al., Understanding Fuel Magnetization and Mix Using Secondary Nuclear Reactions in Magneto-Inertial Fusion Phys. Rev. Lett. 113, 155004 (2014)
  • (8) O.A. Hurricane et al., Inertially confined fusion plasmas dominated by alpha-particle self-heating, Nature Phys. 3720 (2016)
  • (9) R. Ernstorfer et al., The Formation of Warm Dense Matter: Experimental Evidence for Electronic Bond Hardening in Gold, Science 323, 5917 (2009)
  • (10) M.D. Knudson et al., Probing the Interiors of the Ice Giants: Shock Compression of Water to 700 GPa and , Phys. Rev. Lett. 108, 091102 (2012)
  • (11) B. Militzer et al., A Massive Core in Jupiter Predicted from First-Principles Simulations, Astrophys. J. 688, L45 (2008)
  • (12) N.D. Mermin, Thermal Properties of the Inhomogeneous Electron Gas, Phys. Rev. 137, A1441 (1965)
  • (13) T. Sjostrom and J. Daligault, Gradient corrections to the exchange-correlation free energy, Phys. Rev. B 90, 155109 (2014)
  • (14) K. Burke, J.C. Smith, P.E. Grabowski, and A. Pribram-Jones, Exact conditions on the temperature dependence of density functionals, Phys. Rev. B 93, 195132 (2016)
  • (15) M.W.C. Dharma-wardana, Current Issues in Finite-T Density-Functional Theory and Warm-Correlated Matter, Computation 4, 16 (2016)
  • (16) V.V. Karasiev, L. Calderin, and S.B. Trickey, The importance of finite-temperature exchange-correlation for warm dense matter calculations, Phys. Rev. E 93, 063207 (2016)
  • (17) A. Pribram-Jones, P.E. Grabowski, and K. Burke, Thermal Density Functional Theory: Time-Dependent Linear Response and Approximate Functionals from the Fluctuation-Dissipation Theorem, Phys. Rev. Lett. 116, 233001 (2016)
  • (18) D. Kremp, M. Schlanges, and W.D. Kraeft, Quantum Statistics of Nonideal Plasmas, Springer (2005).
  • (19) J. Vorberger, M. Schlanges, W.-D. Kraeft, Equation of state for weakly coupled quantum plasmas, Phys. Rev. E 69, 046407 (2004)
  • (20) U. Gupta and A.K. Rajagopal, Exchange-correlation potential for inhomogeneous electron systems at finite temperatures, Phys. Rev. A 22, 2792 (1980)
  • (21) F. Perrot and M.W.C. Dharma-wardana, Exchange and correlation potentials for electron-ion systems at finite temperatures, Phys. Rev. A 30, 2619 (1984)
  • (22) K.S. Singwi, M.P. Tosi, R.H. Land, and A. Sjölander, Electron Correlations at Metallic Densities, Phys. Rev. 176, 589-599 (1968)
  • (23) S. Tanaka and S. Ichimaru, Thermodynamics and Correlational Properties of Finite-Temperature Electron Liquids in the Singwi-Tosi-Land-Sjölander Approximation, J. Phys. Soc. Jpn. 55, 2278-2289 (1986)
  • (24) T. Sjostrom, and J. Dufty, Uniform Electron Gas at Finite Temperatures, Phys. Rev. B 88, 115123 (2013)
  • (25) P. Vashishta and K.S. Singwi, Electron Correlations at Metallic Densities. V, Phys. Rev. B 6, 875 (1972)
  • (26) S. Dutta and J. Dufty, Classical representation of a quantum system at equilibrium: Applications, Phys. Rev. B 87, 032102 (2013)
  • (27) M.W.C. Dharma-wardana and F. Perrot, Simple Classical Mapping of the Spin-Polarized Quantum Electron Gas: Distribution Functions and Local-Field Corrections, Phys. Rev. Lett. 84, 959-962 (2000)
  • (28) D.M. Ceperley, Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67, 279-355 (1995)
  • (29) E.Y. Loh, J.E. Gubernatis, R.T. Scalettar, S.R. White, D.J. Scalapino and R.L. Sugar, Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301-9307 (1990)
  • (30) M. Troyer and U.J. Wiese, Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations, Phys. Rev. Lett. 94, 170201 (2005)
  • (31) E.W. Brown, B.K. Clark, J.L. DuBois and D.M. Ceperley, Path-Integral Monte Carlo Simulation of the Warm Dense Homogeneous Electron Gas, Phys. Rev. Lett. 110, 146405 (2013)
  • (32) D.M. Ceperley, Fermion Nodes, J. Stat. Phys. 63, 1237-1267 (1991)
  • (33) E.W. Brown, J.L. DuBois, M. Holzmann and D.M. Ceperley, Exchange-correlation energy for the three-dimensional homogeneous electron gas at arbitrary temperature, Phys. Rev. B. 88, 081102(R) (2013)
  • (34) V.V. Karasiev, T. Sjostrom, J. Dufty and S.B. Trickey, Accurate Homogeneous Electron Gas Exchange-Correlation Free Energy for Local Spin-Density Calculations, Phys. Rev. Lett. 112, 076403 (2014)
  • (35) B. Militzer and E.L. Pollock, Variational density matrix method for warm, condensed matter: Application to dense hydrogen, Phys. Rev. E 61, 3470-3482 (2000)
  • (36) B. Militzer, Ph.D. dissertation, University of Illinois at Urbana-Champaign (2000)
  • (37) V.S. Filinov, Cluster expansion for ideal Fermi systems in the ‘fixed-node approximation’, J. Phys. A: Math. Gen. 34, 1665-1677 (2001)
  • (38) V.S. Filinov, Analytical contradictions of the fixed-node density matrix, High Temp. 52, 615-620 (2014)
  • (39) Fionn D. Malone, N.S. Blunt, Ethan W. Brown, D.K.K. Lee, J.S. Spencer, W.M.C. Foulkes, and James J. Shepherd, Accurate Exchange-Correlation Energies for the Warm Dense Electron Gas, Phys. Rev. Lett. 117, 115701 (2016)
  • (40) T. Schoof, M. Bonitz, A.V. Filinov, D. Hochstuhl and J.W. Dufty, Configuration Path Integral Monte Carlo, Contrib. Plasma Phys. 51, 687-697 (2011)
  • (41) T. Schoof, S. Groth and M. Bonitz, Towards ab Initio Thermodynamics of the Electron Gas at Strong Degeneracy, Contrib. Plasma Phys. 55, 136-143 (2015)
  • (42) T. Schoof, S. Groth, J. Vorberger and M. Bonitz, Ab Initio Thermodynamic Results for the Degenerate Electron Gas at Finite Temperature, Phys. Rev. Lett. 115, 130402 (2015)
  • (43) T. Dornheim, S. Groth, A. Filinov and M. Bonitz, Permutation blocking path integral Monte Carlo: a highly efficient approach to the simulation of strongly degenerate non-ideal fermions, New J. Phys. 17, 073017 (2015)
  • (44) T. Dornheim, T. Schoof, S. Groth, A. Filinov, and M. Bonitz, Permutation Blocking Path Integral Monte Carlo Approach to the Uniform Electron Gas at Finite Temperature, J. Chem. Phys. 143, 204101 (2015)
  • (45) N.S. Blunt, T.W. Rogers, J.S. Spencer and W.M. Foulkes, Density-matrix quantum Monte Carlo method, Phys. Rev. B 89, 245124 (2014)
  • (46) F.D. Malone et al., Interaction Picture Density Matrix Quantum Monte Carlo, J. Chem. Phys. 143, 044116 (2015)
  • (47) T. Dornheim, S. Groth, T. Sjostrom, F.D. Malone, W.M. Foulkes, and M. Bonitz, Ab Initio Quantum Monte Carlo Simulation of the Warm Dense Electron Gas in the Thermodynamic Limit, Phys. Rev. Lett. 117, 156403 (2016)
  • (48) T. Dornheim, H. Thomsen, P. Ludwig, A. Filinov, and M. Bonitz, Analyzing Quantum Correlations made simple, Contrib. Plasma Phys. 56, 371 (2016)
  • (49) V.S. Filinov, V.E. Fortov, M. Bonitz and Zh. Moldabekov, Fermionic path integral Monte Carlo results for the uniform electron gas at finite temperature, Phys. Rev. E 91, 033108 (2015)
  • (50) N.S. Blunt, Ali Alavi, and G.H. Booth, Krylov-Projected Quantum Monte Carlo Method, Phys. Rev. Lett. 115, 050603 (2015)
  • (51) S. Groth, T. Schoof, T. Dornheim, and M. Bonitz, Ab Initio Quantum Monte Carlo Simulations of the Uniform Electron Gas without Fixed Nodes, Phys. Rev. B 93, 085102 (2016)
  • (52) T. Dornheim, S. Groth, T. Schoof, C. Hann, and M. Bonitz, Ab initio quantum Monte Carlo simulations of the Uniform electron gas without fixed nodes: The unpolarized case, Phys. Rev. B 93, 205134 (2016)
  • (53) S. Chiesa, D.M. Ceperley, R.M. Martin, and M. Holzmann, Finite-Size Error in Many-Body Simulations with Long-Range Interactions, Phys. Rev. Lett. 97, 076404 (2006)
  • (54) L.M. Fraser et al., Finite-size effects and Coulomb interactions in quantum Monte Carlo calculations for homogeneous systems with periodic boundary conditions, Phys. Rev. B 53, 1814 (1996)
  • (55) N.D. Drummond, R.J. Needs, A. Sorouri and W.M.C. Foulkes, Finite-size errors in continuum quantum Monte Carlo calculations, Phys. Rev. B 78, 125106 (2008)
  • (56) C. Lin, F.H. Zong and D.M. Ceperley, Twist-averaged boundary conditions in continuum quantum Monte Carlo algorithms, Phys. Rev. E 64, 016702 (2001)
  • (57) H. Kwee, S. Zhang, and H. Krakauer, Finite-Size Correction in Many-Body Electronic Structure Calculations, Phys. Rev. Lett. 100, 126404 (2008)
  • (58) D. Chandler, and P.G. Wolynes, Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids, J. Chem. Phys. 74, 4078 (1981)
  • (59) N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller and E. Teller, Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087 (1953)
  • (60) M. Takahashi and M. Imada, Monte Carlo Calculation of Quantum Systems, J. Phys. Soc. Jpn. 53, 963-974 (1984)
  • (61) V.S. Filinov et al., Thermodynamic Properties and Plasma Phase Transition in dense Hydrogen, Contrib. Plasma Phys. 44, 388-394 (2004)
  • (62) A.P. Lyubartsev, Simulation of excited states and the sign problem in the path integral Monte Carlo method, J. Phys. A: Math. Gen. 38, 6659–6674 (2005)
  • (63) M. Takahashi and M. Imada, Monte Carlo of Quantum Systems. II. Higher Order Correction, J. Phys. Soc. Jpn. 53, 3765-3769 (1984)
  • (64) S.A. Chin and C.R. Chen, Gradient symplectic algorithms for solving the Schrödinger equation with time-dependent potentials, J. Chem. Phys. 117, 1409 (2002)
  • (65) K. Sakkos, J. Casulleras and J. Boronat, High order Chin actions in path integral Monte Carlo, J. Chem. Phys. 130, 204109 (2009)
  • (66) S.A. Chin, High-order Path Integral Monte Carlo methods for solving quantum dot problems, Phys. Rev. E 91, 031301(R) (2015)
  • (67) G.H. Booth, A.J. Thom and A. Alavi, Fermion Monte Carlo without fixed nodes: A game of life, death, and annihilation in Slater determinant space, J. Chem. Phys. 131, 054106 (2009)
  • (68) J.S. Spencer, N.S. Blunt and W.M.C. Foulkes, The sign problem and population dynamics in the full configuration interaction quantum Monte Carlo method, J. Chem. Phys. 136, 054110 (2012)
  • (69) M.H. Kolodrubetz, J.S. Spencer, B.K. Clark and W.M.C. Foulkes The effect of quantization on the full configuration interaction quantum Monte Carlo sign problem, J. Chem. Phys 138, 024110 (2013)
  • (70) J.J. Shepherd, G.E. Scuseria and J.S. Spencer, Sign problem in full configuration interaction quantum Monte Carlo: Linear and sublinear representation regimes for the exact wave function, Phys. Rev. B 90, 155130 (2014)
  • (71) D. Cleland, G.H. Booth and A. Alavi, Communications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum Monte Carlo, J. Chem. Phys 132, 041103 (2010)
  • (72) J.J. Shepherd, G.H. Booth and A. Alavi, Investigation of the full configuration interaction quantum Monte Carlo method using homogeneous electron gas models, J. Chem. Phys 136, 244101 (2012)
  • (73) J.J. Shepherd, G.H. Booth, A. Grüneis and A. Alavi, Full configuration interaction perspective on the homogeneous electron gas, Phys. Rev. B 85, 081103 (2012)
  • (74) John P. Perdew, Climbing the ladder of density functional approximations, MRS Bulletin 38 (9), 743-750 2013
  • (75) S. Groth, T. Dornheim, and M. Bonitz, Free Energy of the Uniform Electron Gas: Analytical Models vs. first principle simulations, to be published
  • (76) A. Filinov, and M. Bonitz, Collective and single-particle excitations in 2D dipolar Bose gases, Phys. Rev. A 86, 043628 (2012)
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description