Virial coefficients of the Uniform Electron Gas from Path Integral Monte Carlo Simulations

The properties of plasmas in the low-density limit are described by virial expansions. Analytical expressions are known from Green's function approaches only for the first three virial coefficients. Accurate path integral Monte Carlo (PIMC) simulations have recently been performed for the uniform electron gas, allowing the virial expansions to be analyzed and interpolation formulas to be derived. The exact expression for the second virial coefficient is used to test the accuracy of the PIMC simulations and the range of validity of the interpolation formula of Groth {\it et al.}~[Phys.~Rev.~Lett.~\textbf{119}, 135001 (2017)]. We discuss the fourth virial coefficient, which is of interest, e.g., for properties of solar plasmas, but has not yet been precisely known. Combining PIMC simulations with benchmarks from exact results of the virial expansion would allow us to obtain precise results for the equation of state (EoS) in a wide range of parameters.


I. INTRODUCTION
The thermodynamic properties of Coulomb systems in a wide region of density and temperature are of high interest with respect to various applications.A particularly important regime is given by so-called warm dense matter (WDM) [1], which naturally occurs in a gamut of astrophysical objects such as giant planet interiors [2] and brown dwarfs [3].Moreover, WDM plays an important role in technological applications such as the discovery and synthesis of materials [4][5][6], hot-electron chemistry [7], and inertial confinement fusion [8].As a result, WDM is actively realized in experiments at various research facilities such as the National Ignition Facility (NIF) [9], the Linac Coherent Light Source (LCLS) [10].and the Omega laser facility [11] in the USA, or the European XFEL in Germany [12]; a topical overview of different relevant experimental techniques has been presented by Falk [13].At the same time, we stress that a rigorous theoretical description of such extreme states of matter is indispensable to interpret experimental measurements [14,15], and to guide the development of new set-ups [16,17].
In recent years, new possibilities to obtain results for the thermodynamic properties beyond perturbation theory have arisen, applying numerical simulations to solve the basic expressions [1,16,18,19].Density-functional theory has been successfully applied to evaluate properties of warm dense matter [20][21][22][23], but a main deficit of it is that electron-electron interaction is treated in a certain approximation [24,25].Therefore, computationally more involved path integral Monte Carlo (PIMC) simulations [19,[26][27][28] are of growing interest since they allow the correct treatment of electron-electron interaction.
As a simple example, we consider the homogeneous electron gas (uniform electron gas, UEG [19,29,30]), where the electrons move over a positively charged background which is added to ensure charge neutrality.The electronic part of the Hamiltonian is given by where m and e are the electron mass and charge, ϵ 0 is the permittivity of the vacuum, and pi and ri denote the momentum and position operators of the i-th electron.
In thermodynamic equilibrium, the state of the plasma is determined by the temperature T in addition to the number density n = N/Ω (with Ω being the volume) or the corresponding chemical potential µ.Note that we consider the unpolarised UEG throughout, where n ↑ = n ↓ = n/2 so that both spin directions have the same density.The relationships between the various state variables such as internal energy U , free energy F , entropy S, pressure P , etc., are called equations of state (EoS).All thermodynamic properties can be derived from a thermodynamic potential; F (Ω, N, T ) as function of Ω, N, T constitutes an example.We note that EoS databases constitute key input for a host of practical applications such as the modeling of laser fusion [8] or the description of astrophysical objects [1,3,31].
Correlations appear for the plasma owing to the Coulomb interaction term (1) proportional to e 2 .No closed-form solutions are known, and we must perform approximations (or use numerical techniques [16]) to solve this many-body problem.We discuss two possibilities: (i) Perturbation expansion with respect to e 2 .We obtain analytic expressions for arbitrary orders of e 2 in terms of noninteracting equilibrium correlation functions, which can be easily evaluated using Wick's theorem.However, we have no proof of the convergence of this series expansion and no error estimate.In order to make this analytical approach more efficient, the methods of thermodynamic Green's functions and Feynman diagram technique were elaborated [32,33].The perturbation approach is improved by performing partial summations corresponding to special concepts such as the introduction of the quasiparticle picture (self-energy Σ), screening of the potential (polarization function Π), or formation of bound states (Bethe-Salpeter equation).This leads to useful results for the properties of the plasma in a wide range of T and n.However, as characteristic for perturbative approaches, exact results can be found only in some limiting cases.
(ii) In principle, an accurate evaluation of thermodynamic potentials is possible using path-integral Monte Carlo (PIMC) simulations, see Refs.[18,19,34] and references therein.The shortcomings of this approach include the relatively small number of particles (a few dozen up to one hundred at the present time [35,36]) and the sign problem for fermions [37,38].Over recent years, this emerging approach has been put forward together with improving computer facilities.At present, accurate calculations have been performed mainly for the UEG over a broad range of parameters [28,39].
The UEG is the simplest example.In a next step, calculations for the two-component hydrogen plasma would be of interest for both thermodynamics and transport properties [27,[40][41][42].There are some low-density results, see Militzer [27] or Filinov and Bonitz [41] and further references given in these works.However, highprecision PIMC simulations for hydrogen plasmas in the low-density region, which allow the extraction of higher order virial coefficients, are presently not available.
In this work, we investigate in detail exact virial expansions, which are of considerable value as a rigorous benchmark for numerical methods in certain limits, and as a useful constraint for (semi-)analytical EoS interpolations.To this end, we present new PIMC simulations for the UEG under extreme conditions, i.e. at very low densities and very high temperatures.We investigate the virial expansion and discuss higher order virial coefficients not considered in previous publications [28].In particular, we discuss the high temperature limit of the fourth virial coefficient.A comparison is made with the interpolation formula [43] and the limits of its applicability are shown.
The paper is organized as follows: A brief introduction to the virial expansion of the mean potential energy is given in Sec.II.Effective virial coefficients and virial plots are introduced.PIMC simulations at high temperatures and small densities are presented in Sec.III.An interpolation formula [43] is shown in Sec.IV, and the second virial coefficient is considered as a benchmark.The fourth virial coefficient is analysed in Sec.V, where a high-temperature approximation is given and compared with PIMC simulations.The exact temperature dependence of the fourth virial coefficient is not yet known, but remains a challenge for future PIMC simulations, as we conclude in Sec.VI.

II. VIRIAL COEFFICIENTS FROM ANALYTICAL APPROACHES
A. Virial expansions for the UEG Using the method of thermodynamic Green's functions from quantum statistics, the virial expansion of the free energy of the UEG is written as see Refs.[28,33] where expressions for the lowest virial coefficients F i are also given.The mean potential energy V is given by (for the relation to the internal energy see Ref. [44]).
From the virial expansion of F (T, Ω, N ), see Ref. [33], we get the following virial expansion of with the variables Here ζ(x) denotes the Riemann zeta function, and C = 0.57721 . . . is Euler's constant.We express this expansion in terms of T, n and introduce atomic units ℏ = m = e 2 /4πϵ 0 = 1 (see Appendix A), so that k B T is measured in Hartree (Ha) and n in electrons per The virial expansion of the specific mean potential en- No closed expression for v 4 (T ) is known.

B. The effective second virial coefficient
In Refs.[28,45], a method for extracting the virial coefficients from data was presented.We demonstrate this approach for the second virial coefficient v 2 (T ) for which the exact expression ( 7) is known.Since in the low-density limit the lowest virial coefficients dominate the function v(T, n) (6), we subtract the "trivial" contributions of v 0 (T ) (Debye term) and v 1 (T ).The remaining part is then dominated by v 2 (T ) in the low-density limit.
To extract the value of v 2 (T ) from numerical (or measured) results for v(T, n), we consider isotherms and calculate an effective, density-dependent second virial coefficient We have the result v 2 (T ) = lim n→0 v eff 2 (T, n).The density dependence of v eff 2 (T, n) in the low-density limit is given according to Eq. ( 6) as Ha ), the isotherms should meet the co-ordinate at v 2 (T ) and the slope is v 3 (T ).The linear pattern is violated when higher virial coefficients become relevant.Note that both v 2 (T ) and v 3 (T ) are known for the UEG according to Eq. ( 7).Corresponding plots for three isotherms are shown below in Figs.
2-3, T Ha = 0.589307, 0.294653, and 100.We will apply this method of the virial plot to PIMC simulations v PIMC (T, n) to obtain the values v eff,PIMC 2 (T, n) according to Eq. (8).It is clear that this method of extracting virial coefficients requires a high precision of the calculated data, since we are analysing the difference of large numbers.This is because the lower virial coefficients, such as the Debye term, dominate the low-density limit of the potential energy density v(T, n).

III. PIMC SIMULATIONS FOR THE UEG A. PIMC simulations at high temperatures and small densities
To compare with the virial expansion, we have to consider high temperatures T Ha ≥ 1 and n Ha ) ≪ 1 so that the contributions of higher order virial coefficients are sufficiently small.This requires us to go beyond the conditions that had been explored in Ref. [28].Clearly, the fermion sign problem [37] does not pose an obstacle as quantum degeneracy effects become negligible in the limit of high temperature.Instead, the biggest challenge is given by finitesize effects in the simulation data, which are substantial Extrapolating PIMC results for the interaction energy of the UEG to the thermodynamic limit (TDL) at rs = 2 and Θ = 217.204(THa = 100).Left: Raw PIMC results for the interaction energy per particle V /N (green crosses), and finite-size corrected values (red circles) as a function of system size.Right: Magnified segment around the finite-size corrected results; the solid black lines show empirical linear and constant fits, and the blue cross depicts the extrapolated result in the limit of in this regime.This is illustrated in the left panel of Fig. 1, where we show raw PIMC simulation results for the interaction energy per particle V /N at r s = 2 and Θ = 217.204(for r s , Θ, see Appendix A).Indeed, the dependence on the system size is of the order of 100%.To overcome this bottleneck, we employ the finite-size correction scheme developed in Ref. [35], which constitutes a finite-temperature version of the approach originally introduced by Chiesa et al. [46].We refer the interested reader to the overview [19] for a more detailed discussion.The thus corrected interaction energies are shown as the red circles in Fig. 1 and exhibit a drastically reduced dependence on the system size; the residual error is of the order of ∼ 0.1%.In the right panel, we show a magnified segment around these corrected values, and the solid black lines show empirical fits based on simple linear and constant functional forms.In practice, we give the final result based on the linear extrapolation, and the associated uncertainty is computed from the difference between the two solid lines, see the blue cross in Fig. 1.All PIMC results shown in this work have been obtained based on this procedure.

B. Calculations with PIMC simulation results
We show PIMC simulation results for the UEG as published in Ref. [28] in Tab.I.For our analysis, we need high accuracy data because we consider small differences of large numbers.We focus only on results for r s = 20 and Θ = 128 and Θ = 64.There, the virial coefficients are relatively large, and higher orders of the virial expansion are not too dominant.
The virial plots for the corresponding two isotherms are presented in Fig. 2, where Ha ).Isotherms are presented because the virial coefficients describe the expansion with respect to density n at fixed T .For comparison, the benchmarks v 2 (T )+v 3 (T )n Ha ) are also shown.There is a nice agreement.Deviations may be explained by the contribution of higher virial coefficients for the analytical results.In addition, the PIMC data have also uncertainties expressed by error bars.
To demonstrate the limiting behavior given in the virial plot by the linear relation (9), neglecting higher order terms O[n 1/2 ], more PIMC simulation data would be of interest.In this work, we performed high precision PIMC simulations for additional three parameter values.The results are shown in Tab.II.
In Tab.II, No. 1 and 2 belong to the isotherms of Tab.I, No. 1 and 2, but at lower density.As seen in Fig. 2, the PIMC simulations are consistent with the virial expansion.However, the error bars are quite large so that the extrapolation n → 0 to extract the second virial coefficient from PIMC simulations has also a large error.
In Tab.II, No. 3 belongs to the isotherm T Ha = 100 to study the high-temperature limit.As shown in Fig. 3, the PIMC simulation data are also consistent with the virial benchmark.However, the error bars are large, too.

IV. INTERPOLATION FORMULAS FOR THERMODYNAMIC PROPERTIES OF THE UEG
The PIMC simulations are computationally very expensive.Instead of performing time-consuming calcu- lations for each parameter value, interpolation formulas have been worked out which allow to reproduce the results for each parameter value within a given accuracy.
Because the limiting behavior of the free energy is known at low and at high density, Padé expressions can be used to good effect [43,47,48].We can also test these interpolation formulas with respect to their accuracy and the parameter range where they can be used.In particular, we study whether they can be used instead of PIMC simulations to extract a value for the virial coefficient such as v 2 (T ) or v 4 (T ).The GDSMFB interpolation formula for the XC free energy density of the spin-unpolarized UEG is [43] The coefficients a, b, c, d, e are again Padé formulae with respect to temperature given in the Supplemental material to [43] [see also Appendix B].
The exact relationship between the exchangecorrelation free energy and the potential energy is given as so that We find This expression will be used to calculate v eff,GDSMFB 2 (T, n) according to Eq. ( 8).The corresponding results are shown in Figs. 2 and 3 as blue lines.
It is obvious that the interpolation formula shows strong deviations in the low-density limit.The reason is that the Padé formula (13) is not constructed to reproduce the v 2 (T ) so that the analytical behavior of the Padé formula in the low-density limit is not consistent with the virial expansion.This discrepancy shows a limit of applicability of the interpolation formula.However, because in the low-density region the lowest order virial terms (e.g. the Debye shift) dominate, the error of the interpolation formula becomes small if these lowest orders are correctly included.
Is it possible to extract the virial coefficients from the Padé formula (13)?In Fig. 3 we show the values for the effective second virial coefficient v eff,GDSMFB 2 (T, n), Eq. ( 8), for the isotherm T Ha = 100.It is clearly shown that in the low-density limit (below n 1/2 B ln(4πn B /T 2 Ha ) = 1) the benchmark of the second virial coefficient is not reached.It is also shown that there the value of the interpolation formula lies outside the error bars of the PIMC simulation.We see that in the limit n → 0 the results for v eff,GDSMFB Ha ) > 2 gives a limit at n = 0 of -0.0128 which deviates from the exact value by about 10 %.
In conclusion, from Figs. 2 & 3 we see that PIMC simulations become difficult in the low-density region and the error bars become large.The low-order virial coefficients may be considered as a benchmark for the simulation.The interpolation formula (13) describes the virial plot for v eff 2 in a certain approximation only in an interme-diate parameter range.At high densities, higher orders of the virial expansion become important.At very low densities, the analytical behavior of the interpolation formula is not able to reproduce the exact virial coefficients.However, they may be estimated in certain approximation, if a linear behavior can be seen in the virial plot, see Fig. 3.

V. THE FOURTH VIRIAL COEFFICIENT
A. High-temperature limit and PIMC simulation data Analytical expressions for v 4 (T ) are not yet known.Approximations considering special classes of diagrams have been obtained within Green's function approaches.For instance, considering the diagrams of lowest order with respect to interaction, in Refs.[33,49], the following contribution to the fourth virial coefficient has been given [in atomic units, see Eqs. ( 6), (7)] This result leads to a high temperature behavior lim T →∞ v 4 (T ) ∝ T −2 if no other diagrams contribute to this limit.We follow the method explained for the second virial coefficient.At fixed T , in the low-density limit the lowest virial coefficient will dominate because of the analytical behavior near n = 0. Thus, subtracting the (k − 1) lowest virial coefficients from the thermodynamic quantity, the remaining part allows to determine the next virial coefficient v k (T ) [45].We use a virial plot for v eff 3 (T, n) (T, n), Eq. ( 15), plotted as a function of −1/ ln(4πnB/T 2 Ha ).The slope of v eff 3 (T, n) determines v4(100).The linear relation Eq. ( 15) with v3(100) = −8.352× 10 −7 is denoted as virial 3+4.The dashed curve corresponds to v GDSMFB 4 = 0.00038.The dash-dotted curve corresponds to v4 = 0.00129, Eq. ( 14).(b) The effective second virial coefficient v eff,GDSMFB 2 (T, n), Eq. ( 8), plotted as function of n 1/2 B .From analytical approaches (Eqs.( 7) and ( 14)) follows v eff 2 (100) ≈ −0.01147 + 0.00129n (shown as dash-dotted line).The dashed line represents a linear fit.(Atomic units used.)defined as In the low-density limit, the density dependence of v eff 3 (T, n) is given according to Eq. ( 6) as Thus, in the virial plot where v eff 3 (T, n) is shown as a function of 1/ ln(4πn B /T 2 Ha ), isotherms should meet the co-ordinate at v 3 (T ), and the slope is v 4 (T ).
We discuss here the high-temperature region because we expect that the odd virial coefficients become small in this limit T → ∞, as seen for the lowest virial coefficients (7), see also Tab.V in App. C. We expect a wider range of the linear relation ( 16) if the higher virial coefficients, in particular v 5 (T ), are small.As example, we consider T Ha = 100, see Fig. 4(a).The dash-dotted curve denotes the virial expansion with the exact value for v 3 (100), Eq. ( 7), and the approximation (14) for v 4 (100).In addition, a PIMC simulation is also shown (No. 3 from Tab. II).Within the error bars, the result agrees with the virial expansion.However, for this approach to extract virial coefficients from PIMC simulations, more data in the lowdensity region with higher accuracy are required which are not yet available.

B. Fourth virial coefficient from interpolation formulas
Because the full T dependence of v 4 (T ) is not yet known from the Green's function approach, it would be of interest to obtain results from simulations.As shown for v 2 (T ) in the previous section, see also [28], high-accurate PIMC simulations may be used to extract this quantity.However, they are not yet available.
To make some estimations with respect to v 3 (T ), we may use the interpolation formula (13) instead of the exact approach using PIMC simulations.Because this GDSMFB interpolation formula is only an approximation, significant deviations may occur.
Using v GDSMFB (T, n) (13) as input for v, we calculate v eff,GDSMFB 3 (T, n) according Eq. (15), see Appendix C.These values are plotted in Fig. 4(a) as a function of −1/ ln(4πn B /T 2 Ha ).As discussed above, in the high temperature region considered here, the odd virial coefficients give only small contributions.
From this curve, the linear extrapolation is possible for the values at larger densities.The value of v 4 (100) is estimated to be close to 0.00038.The values at smaller densities cannot be used for the extrapolation because the interpolation formula is an approximation, and deviations yield large effects for small densities as already seen for v eff,GDSMFB 2 (T, n) in Figs. 2 and 3.In Fig. 4(a), we also show the approximation (14) with the value v 4 (T Ha = 100) ≈ 0.001295.Compared with the PIMC value of Tab.II also shown in Fig. 4(a), we find that both results are consistent.The GDSMFB interpolation formula is not consistent with PIMC simulations in this parameter region.In addition, the extracted approx-imation for v 4 (100) is different from the approximation (14) which should be valid in the high temperature limit.

C. Generalized virial plots
The virial plots use an abscissa which gives a linear relation for the next higher virial coefficient so that this next virial coefficient is extracted from the slope of the isotherms at zero density.It may happen that the virial expansion contains terms which are very small so that these terms are not relevant.For instance, at high temperatures, v 3 (T ) becomes very small because it behaves ∝ T −7/2 according Eq. ( 7).As shown in Fig. 3, the slope of virial 2+3 is near to zero, but the interpolation formula indicates a significant increase with density.This is the contribution of higher order virial coefficients.We assume that in the high-temperature limit, the virial term v 3 (T ) can be neglected so that the dominant contribution to the virial expansion in the low-density range follows from v 4 (T ).
To extract this leading virial coefficient v 4 (T ) from v eff 2 (T, n), Eq. ( 9), we introduce a generalized virial plot where the abscissa is n 1/2 .If we observe a linear behavior, the slope determines v 4 (T ).third virial coefficient gives a contribution only for very low densities, leading to an off-set of the linear extrapolation to n = 0, but may be neglected if v 3 (T ) is small.
In Fig. 4(b), we show this generalized virial plot for v eff 2 (T, n).The dashed line corresponds to v with the approximation (14) for v 4 (T ).Because v 3 (100) is very small, the off-set at very low densities is not seen.For comparison, the interpolation formula is also shown, and a linear behavior is seen.The extrapolation to n = 0 misses the exact value v 2 (100) as also discussed above.There it was argued that the interpolation formula does not contain this benchmark by construction.However, if we assume that the interpolation formula gives a reasonable approximation in a wide range of parameter values, the linear behavior in the generalized virial plot is clearly seen.The extracted slope v GDSMFB 4 (100) = 0.00135 is in good agreement with the value 0.00129 from the approximation (14).
In addition to the isotherm T Ha = 100, we studied also other isotherms ranging from T Ha = 50 to T Ha = 400.The extracted slope v GDSMFB 4 (T ) show the 1/T 2 behavior in accordance with Eq. ( 14).The investigation of the uniform electron gas is of interest not only for the discussion of the exchange-correlation term of the energy-density functional in DFT calculations, for which analytical formulae have been derived by Groth, Dornheim, and Bonitz [19,47].It is also a prerequisite for the treatment of the more interesting case of a two-component plasma, e.g., the hydrogen plasma.For instance, the equation of state at low densities is of interest in helioseismology [50] where the fourth virial coefficient v 4 (T ) is relevant [51].In this context, the high-temperature limit of v 2 (T → ∞) and the relation to v 4 (T ) has been discussed in Refs.[28,52].For a discussion of the fourth virial coefficient v 4 (T ) of the hydrogen plasma see also Alastuey and Ballenegger [53,54].The correct determination of the fourth virial coefficient v 4 (T ) of the UEG is an important prerequisite for finding expressions for the fourth virial coefficient F 4 (T ) in the free energy (2) associated with the density power n 5/2 .However, we leave the discussion of this question to future work.

VI. CONCLUSIONS
Quantum statistics gives us exact expressions for thermodynamic and transport properties of plasmas in terms of equilibrium correlation functions, but their evaluation is a complex problem in many-particle physics.Numerical simulations are becoming more accurate as computer capacity increases.However, they need to be checked for their limitations, such as size effects, but also for fundamental problems such as the correct description of electron-electron collisions in the framework of DFT or strategies to deal with the sign problem in PIMC simulations.PIMC simulations are expected to provide an adequate description of electron-electron interactions, but are currently unable to solve complex plasmas such as multiply charged ions at low temperatures.
The use of analytical results for the virial expansion of thermodynamic properties as a benchmark for PIMC calculations for the uniform electron gas is demonstrated.In particular, we show that high-precision PIMC simulations confirm the correct form of the virial expansion that has been recently discussed [28].It also seems possible to obtain numerical values for higher virial coefficients, in particular the interesting virial coefficient v 4 (T ) for the order n 5/2 of the free energy.These values can be considered as exact results in plasma physics.
Analytical theory gives us exact results in limiting cases as benchmarks.These can be used to obtain results for parameter ranges where numerical simulations are not efficient, e.g. in the range of low densities.Virial expansions are used to control theories and numerical simulations.They are of interest for the construction of interpolation formulas.
The UEG is a comparatively simple case where PIMC simulations are possible with high accuracy.It will be interesting to extend the present considerations to a twocomponent system like the hydrogen plasma [27,[40][41][42], the positronium plasma or the electron-hole plasma.

D
. The n 5/2 term

Table II .
New PIMC data, corresponding parameter values and effective second virial coefficients (atomic units, see Appendix A).