Ab initio simulation of warm dense matter

Warm dense matter (WDM) has become a mature research field on the boarder of plasma physics and condensed matter physics, e.g., Refs. 1–4. There are many examples in astrophysics such as the plasma-like matter in brown and white dwarf stars,^5–7 giant planets, e.g., Refs. 8–13, and the outer crust of neutron stars.^14,15 Warm dense matter is also thought to exist in the interior of Earth.¹⁶ In the laboratory, WDM is being routinely produced via laser or ion beam compression or with Z-pinches, see Ref. 17 for a recent review article. Among the facilities, we mention the National Ignition Facility at the Lawrence Livermore National Laboratory;^18,19 the Z-machine at the Sandia National Laboratory;^20,21 the Omega laser at the University of Rochester;²² the Linac Coherent Light Source (LCLS) in Stanford;^23,24 the European free electron laser facilities FLASH and X-FEL in Hamburg, Germany;^25,26 and the upcoming FAIR facility at GSI Darmstadt, Germany.^27,28 A particularly exciting application is inertial confinement fusion^18–20 where electronic quantum effects are important during the initial phase. Aside from dense plasmas, many condensed matter systems exhibit WDM behavior – if they are subject to strong excitation, e.g., by lasers or free electron lasers.^29,30

The behavior of all these very diverse systems is characterized by, among others, electronic quantum effects, moderate to strong Coulomb correlations, and finite temperature (FT) effects. Quantum effects of electrons are of relevance at a low temperature and/or if matter is very highly compressed, such that the temperature is of the order of (or lower than) the Fermi temperature (for the relevant parameter range, see Fig. 1 and, for the parameter definitions, see Sec. II).

An important role in the theoretical description of quantum plasmas is being played by the quantum kinetic theory.^31–38 During the last 25 years, improved and generalized quantum kinetic equations have been derived starting from reduced density operators, e.g., Refs. 39 and 40, or nonequilibrium Green functions (NEGFs);^41–44 for text books, see Refs. 40 and 45–47 and references therein. Another direction in quantum plasma theory is first principles computer simulations such as quantum Monte Carlo (QMC),^4,48–55 semiclassical molecular dynamics (SC-MD) with quantum potentials, e.g., Ref. 56, electronic force fields,^57,58 and various variants of quantum MD, e.g., Refs. 59–63.

A recent breakthrough occurred with the application of Kohn–Sham density functional theory (DFT) simulations because they, for the first time, enabled the self-consistent simulation of realistic warm dense matter that includes both plasma and condensed matter phases, e.g., Refs. 64–66. Further developments include orbital-free DFT (OF-DFT) methods, e.g., Refs. 67 and 68, and time-dependent DFT (TD-DFT), e.g., Ref. 69. In DFT simulations, however, a bottleneck is the exchange–correlation (XC) functional for which a variety of options exist, the accuracy of which is often poorly known, what limits the predictive power of the method. This requires tests against independent methods such as quantum Monte Carlo simulations for the electron component⁴ or against electron-ion quantum Monte Carlo.^70–72 Also, the use of finite-temperature functionals was shown to be important^73,74 when the XC-contribution is comparable to the thermal energy, see Ref. 75 for a topical discussion and Ref. 76 for an extensive investigation of hydrogen. One goal of this paper is to present an overview of these results and discuss future research directions.

Motivated by time-resolved experiments, e.g., Ref. 77, the theoretical description of the nonequilibrium dynamics of warm dense matter is attracting increasing interest, e.g., Ref. 78. Time-dependent x-ray Thomson scattering was modeled in Refs. 79 and 80. Here, the powerful methods are quantum kinetic equations^81,82 and nonequilibrium Green functions, e.g., Refs. 83 and 84.

All of the above-mentioned simulation approaches are complex and require substantial amounts of computer time. At the same time, the above-mentioned simulations are currently only feasible for small length scales and simulation durations. Therefore, simplified models that would allow to reach larger length and time scales are highly desirable. Possible candidates are fluid models for quantum plasmas that are obtained via a suitable coarse graining procedure, as in the case of classical plasmas. Quantum hydrodynamics (QHD) models for dense plasmas have experienced high activity since the work of Manfredi and Haas.^85,86 However, their version of QHD involved several assumptions, the validity of which remained open for a long time. Corrections of the coefficients in the QHD equations were recently obtained in Refs. 87 and 88, and a systematic derivation of the QHD equations from the time-dependent Kohn-Sham equations is given in Ref. 89. We also mention a recent alternative approach that is based on the computation of semiclassical Bohm trajectories.⁹⁰ 

The goal of this paper is to present a summary of some of the recent ab initio simulations of the electron gas under warm dense matter conditions, including thermodynamic functions and local field corrections developments. Furthermore, we summarize recent progress in the field of QHD for quantum plasmas. In addition to an overview of recent developments, we present new results for (a) the application of the finite-temperature exchange correlation free energy in DFT simulations of dense hydrogen and carbon (Sec. IV); (b) for the dynamic density response function, $χ(ω,q)$ (Sec. III C); (c) for the screened potential of an ion in a correlated plasma, based on the ab initio QMC input for the local field correction (Sec. V F); and (d) on the dispersion of ion-acoustic modes in a correlated quantum plasma (Sec. V G).

This paper is organized as follows: in Sec. II, we recall the main parameters of warm dense matter and the relevant temperature and density range. Section III presents an overview on recent quantum Monte Carlo simulations followed by finite-temperature DFT results in Sec. IV. WDM out of equilibrium and its treatment via a QHD model is discussed in Sec. V.

Let us recall the basic parameters of warm dense matter:^40,89 the first are the electron degeneracy parameters $θ=kBT/EF$ and $χ=nΛ3$⁠, where $Λ=h/2πmkBT$⁠, is the thermal DeBroglie wavelength, and the Fermi energy of electrons (in 3D) is

where n is the electron density and T is the electron temperature. The ion degeneracy parameter is a factor $(mT/miTi)3/2$ smaller than the one of the electrons and typically negligible for WDM. Second, is the classical coupling parameter of ions $Γi=Qi2/(aikBTi)$⁠, where Q_i is the ion charge and a_i is the mean inter-ionic distance. Furthermore, the quantum coupling parameter (Brueckner parameter) of electrons in the low-temperature limit is

where $a=(4/3πn)−1/3$ denotes the mean distance between two electrons, a_B is the Bohr radius, and $mr=mmi/(m+mi)$ϵ_b are the reduced mass and background dielectric constants, respectively, for hydrogen $mr≈m, ϵb=1$⁠, and $aB=0.529Å$⁠. Note that another way to measure the coupling strength in the degenerate limit, that is directly related to r_s, is via

We can introduce an effective coupling parameter that interpolates between the classical and strongly degenerate limits,

and a simple estimate for the boundary between ideal and nonideal plasmas is the line $Γeff=0.1$ that has been included in Fig. 1. Finally, the degree of ionization of the plasma—the ratio of the number of free electrons to the total (free plus bound) electron number, $αion=n/ntot$⁠, determines how relevant plasma properties are compared to neutral gas or fluid effects.

The parameters $χ,χi$⁠, and $Γi$ are shown in Fig. 1 where we indicate where these parameters equal one. Note that the classical coupling parameter increases with the density whereas the quantum coupling parameter decreases with the density n. We underline that the parameters $aB,EF,Λ,θ$ contain the density of free electrons and $Γi,χi$ the density of free ions. This means the lines of constant $Γi,χ,χi,rs$ shown in Fig. 1 refer to the free electron (ion) density. In cases when the plasma is only partially ionized, the free electron density has to be replaced by $n→αion×n$⁠. The degree of ionization decreases when the temperature is lowered, according to the Saha equation, $αion∼e−|Eb|/kBT$⁠, where E_b denotes the binding energy of the atom, and in Fig. 1 we indicate the line where a classical hydrogen plasma has a degree of ionization of 0.5. Qualitatively, a quantum plasma is found to the right of this line. Figure 2 shows a zoom into the warm dense matter range and also contains lines of constant r_s- and Θ-values.

The uniform electron gas (UEG) is one of the most fundamental model systems in physics.^4,92,93 In particular, the accurate parametrization^94,95 of the ground-state exchange–correlation energy $exc(rs)$⁠, based on ab initio quantum Monte Carlo simulations,⁹⁶ has been essential for the striking success of the density functional theory. While the influence of temperature on the electrons is negligible for most applications in, e.g., condensed matter or chemistry, the recent interest in matter under extreme conditions has led to new demands regarding our understanding of the UEG. More specifically, it has long been known that a thermal DFT^97,98 simulation of warm dense matter, see Sec. IV, requires a parametrization of the exchange–correlation free energy $fxc(rs,θ)$⁠, which explicitly depends both on the density and on the temperature.^75,99

Consequently, many such parametrizations have been presented over the last decades that are based on various approximations such as dielectric theories,^100–105 quantum–classical mappings,^106,107 and perturbative expansions;^108,109 see Refs. 4 and 110, for a topical overview. In addition, Brown et al.¹¹¹ presented the first quantum Monte Carlo results for the warm dense UEG using the restricted path integral Monte Carlo (RPIMC) method, which have subsequently been used as input for many applications,^112–114 most notably the parametrization of $fxc$ by Karasiev et al.¹¹⁵ (KSDT).

While being an important mile stone, these data have been obtained on the basis of the uncontrolled fixed node approximation,¹¹⁶ which has recently been revealed to be surprisingly inaccurate with systematic errors in the exchange–correlation energy exceeding 10%, at a high density and a low temperature.⁵⁴ This unsatisfactory situation has sparked a surge of new developments in the field of fermionic QMC simulations at a finite temperature.^55,117–124 In particular, Groth and co-workers have introduced a combination of two complementary QMC methods—permutation blocking PIMC (PB-PIMC) and configuration PIMC (CPIMC)—that allow for a highly accurate description of the UEG over a broad parameter range without the fixed node approximation. After developing a new finite-size correction scheme,¹²⁵ the same authors presented the first ab initio parametrization of $fxc$ with respect to density, temperature, and spin-polarization covering the entire WDM regime with an unprecedented accuracy of $∼0.3%$⁠.

This is illustrated in Fig. 3, where we show the temperature-dependence of the interaction energy of the spin-polarized UEG for two different densities and compare different parametrizations and data sets. The red solid line corresponds to the parametrization by Groth et al.⁷³ (GDB, referred to as GDSMFB hereafter), which is based on finite-T QMC data (red crosses) for $0.5≤θ≤8$ and temperature-corrected ground-state QMC data (red diamonds) for $θ≤0.25$⁠. The other curves depict the RPIMC data from Ref. 111 (BCDC, blue circles) and a corresponding parametrization (KSDT, solid blue), HNC data and a parametrization thereof by Tanaka¹⁰⁴ (HNC, green squares and line), and the Singwi–Tosi–Land–Sjölander (STLS)-based parametrization by Ichimaru et al.¹⁰¹ At a high density [$rs=0.1$⁠, panel (a)], both dielectric theories agree relatively well with the GDSMFB reference data (see also the bottom panel showing the relative deviation), while the KSDT curve exhibits deviations of up to $Δv/v∼15%$⁠. This is a direct consequence of the absence of RPIMC input data in this regime, and the insufficient finite-size correction of these data, where they are available.^4,125,126 In the WDM regime [r_s = 1, panel (b)], the situation somewhat changes as both, HNC and STLS, become less accurate and exhibit deviations of up to $Δv/v∼5%$ in the relevant temperature range. Moreover, the RPIMC data exhibit systematic deviations from the other curves as they are systematically too large, for a small θ, and too low, in the opposite case. This is due to a combination of the fixed-node approximation and the extrapolation to the thermodynamic limit; see Ref. 4 for an extensive discussion. Interestingly, the KSDT curve is remarkably accurate in the low-temperature limit and does not reproduce the biased RPIMC input data on which it is based. Still, there occur deviations of up to $Δv/v∼8%$⁠, at an elevated temperature.

In the meantime, the availability of the accurate GDSMFB benchmark data has led to a revised version of the KSDT parametrization (denoted as corrKSDT in Refs. 74 and 127) which basically reproduces GDSMFB over the entire WDM regime.⁷⁴ First and foremost, we note that both GDSMFB and corrKSDT are suitable to be used as an exchange–correlation functional on the level of the local density approximation⁷⁶ (the GDSMFB parametrization is available in the libxc library as “GDSMFB,” cf. Sec. IV), and as the basis for more sophisticated functionals such as a temperature-dependent generalized gradient approximation.¹²⁷ This opens up new avenues for DFT simulations of WDM systems without neglecting thermal effects in the XC-functional itself. On the other hand, Karasiev et al.⁷⁴ have found that there occur some potentially unphysical oscillations in quantities, that are derived from $fxc$⁠, such as the specific heat C_V. This is not surprising as C_V involves the second derivative of the fit with respect to the temperature which may contain a large error. In addition, the entropy was found to become negative at strong coupling and low temperature, which, however, is outside of the intended domains of application of both GDSMFB and corrKSDT.

Let us conclude this section by proposing a few possible solutions to the remaining open questions: (i) the unphysical behavior in the second- and higher-order derivatives are most likely a consequence of the functional form of the GDSMFB and corrKSDT parametrizations.⁷⁴ Therefore, addressing this issue would require a modification of the corresponding Padé approximation to automatically fulfill some additional constraints. (ii) The current validity domain of, e.g., the GDSMFBB parametrization to $0≤rs≤20$ can be significantly extended by incorporating the recent ab initio PIMC results for the electron liquid regime (⁠$20≤rs≤100$⁠) by Dornheim et al.¹²⁸ (iii) New ab initio QMC results at a low temperature, $θ∼10−1$⁠, could help to more accurately resolve open thermodynamic questions like the effective mass enhancement,¹²⁹ but are difficult to obtain due to the notorious fermion sign problem.¹³⁰ (iv) Neural networks are known to be valuable as universal function approximators (see also Sec. III B) and can be designed to fulfill all known constraints on the UEG. Moreover, they constitute a handy way to combine data for different quantities from different methods into a single, unified representation.

One important step in going beyond local approximations, such as the local density approximation (LDA) and generalized gradient approximations (GGAs), is to consider the response of an electron gas to an external harmonic perturbation [cf. Eq. (6)], which is described by the density response function

with q and ω being the corresponding wave number and frequency, respectively, and $ṽq=4πe2/q2$ (frequently atomic units are used, then this becomes $4π/q2$⁠) is the Fourier transform of the Coulomb potential. Furthermore, $χ0(q,ω)$ denotes the usual Lindhard function that describes the density response of the ideal (i.e., noninteracting system),⁹³ and the local field correction $G(q,ω)$ entails the full frequency- and wave number resolved information about exchange–correlation effects on χ.¹³¹ For example, setting $G(q,ω)=0$ in Eq. (5) leads to the widely used random phase approximation (RPA), which describes the density response of the electron gas on a mean field level.

Consequently, the local field correction (LFC) is of paramount importance to incorporate nonlocal exchange–correlation effects into other theories, like QHD,^87,88,132 effective potentials,^3,133–138 and the construction of advanced exchange–correlation functionals for DFT^139–142 and time-dependent DFT.⁶⁹ Moreover, it can directly be used to compute important material properties like the stopping power,^143–145 electrical and thermal conductivities,^146,147 and energy transfer rates.¹⁴⁸ Finally, we mention the interpretation of x-ray Thomson scattering (XRTS) experiments,^149–153 e.g., within the Chihara decomposition,¹⁵⁴ which is of paramount importance as a method of diagnostics.

Naturally, there have been many attempts to find suitable approximations for $G(q,ω)$⁠, most commonly within the purview of dielectric theories.^155–159 The first accurate benchmark data for the LFC have been obtained by Moroni et al.^160,161 by performing ground-state QMC simulations of a perturbed electron gas governed by the Hamiltonian

with $ĤUEG$ corresponding to the unperturbed UEG and A being the perturbation amplitude. More specifically, they have performed multiple simulations for a single wave number q to measure the response of the electron gas in dependence of A, which is linear for a small A with $χ(q,0)$ being the slope. While being limited to the static limit (i.e., ω = 0), these data have subsequently been used as an input for the parametrization of $G(q,0)$ by Corradini et al.¹⁶² (CDOP), which, in turn, has been used for many applications, e.g., Refs. 3, 136, 139, 140, 163, and 164.

Recently, Dornheim and co-workers^165,166 have extended the idea behind Eq. (6) to finite temperature using the novel permutation blocking PIMC and configuration PIMC methods, which has allowed to obtain the first ab initio results for the static density response of the UEG under WDM conditions. While being conceptually valid and interesting in their own right, these simulations suffer from a prohibitive computational cost: the full characterization of $G(q,0;rs,θ)$ requires a dense grid of densities, temperatures, and wave numbers, $(rs,θ,q)$⁠. Unfortunately, each such tuple would, in turn, require multiple simulations for different perturbation amplitudes A, and potentially also different N to eliminate possible finite-size effects. Therefore, the aforementioned strategy is valuable to produce an accurate benchmark data at specific points, but cannot be feasibly used to generate the bulk of input data needed for a full description of $G(q,0;rs,θ)$ covering the entire WDM regime.

A different, more convenient route is given by the imaginary-time version of the fluctuation–dissipation theorem,^165,167

with

being the density–density correlation function (also known as the intermediate scattering function) evaluated in the imaginary time $τ∈[0,β]$⁠. Equation (8) can be straightforwardly computed using the standard PIMC method,^168,169 which means that the entire wave-number dependence of $χ(q,0)$ can be obtained from a single simulation of the unperturbed UEG (i.e., setting A = 0 in Eq. (6). The corresponding results for $G(q,0)$ are then computed by solving Eq. (5) for G, i.e.,,¹⁷⁰ 

This strategy—in combination with the efficient finite-size correction introduced in Ref. 166—was recently used by Dornheim et al.¹⁷¹ to obtain extensive new ab initio PIMC results for $G(q,0)$ for $Nparam∼50$ density–temperature combinations covering a significant part of the relevant WDM regime. These data, together with the CDOP parametrization for θ = 0, were subsequently used as input to train a deep neural net, which takes as input a tuple of $(q,rs,θ)$ and predicts as output the corresponding LFC $G(q,0;rs,θ)$ in the range of $0≤θ≤4, 0.7≤rs≤20$⁠, and $0≤q≤5qF$⁠.

A typical result is shown in Fig. 4, where the static LFC is plotted in the wave number-temperature plain for a fixed value of the density parameter, r_s = 5. The black squares depict the ground-state QMC results by Moroni et al.,¹⁶¹ and the dashed blue line depicts the corresponding CDOP parametrization, which incorporates both, the compressibility sum-rule,¹⁰³ for $q→0$⁠, and the exact large-q limit found by Holas.^172,173 The red crosses show the new PIMC data computed from Eq. (7), which is available at a dense grid of wave numbers that is determined by the usual momentum quantization in a finite simulation cell. Finally, the green surface has been evaluated using the neural net published in Ref. 171. Evidently, the machine-learning representation nicely reproduces the input data where they are available, and smoothly interpolates in between. Moreover, its capability to predict $G(q,0;rs,θ)$ has been validated against independent benchmark data that had not been included in the training procedure, see Ref. 171 for details.

Let us conclude this section by explicitly investigating the impact of thermal excitations on the static LFC. To this end, we plot G at a metallic density, r_s = 4, for five different temperatures in Fig. 5. The solid green curve corresponds to the zero-temperature limit, where G is accurately represented by CDOP. Upon increasing the temperature to $θ=0.5$ (red dots), the LFC essentially remains unchanged, for $q ≲ 3qF$⁠, but exhibits a significant drop and an apparent saddle point, for large wave numbers. At the Fermi temperature (blue dashes), $G(q,0;rs,θ)$ exhibits an even more interesting behavior: while it is approximately equal to the θ = 0 curve, for $q ≲ 2qF$⁠, there appears a complicated shape with a maximum around $q≈2.7qF$⁠, a subsequent minimum at $q≈4.5qF$⁠, and a positive large wave number tail. Finally, further increasing the temperature to θ = 2 (dash-dotted black curve) and θ = 4 (long-dashed brown curve) leads to significant thermal effects, even for small q values, and G exhibits a pronounced maximum at intermediate wave numbers, followed by a tail with a negative slope.¹⁷¹ 

For completeness, we mention that new ab initio results for $G(q,0;rs,θ)$⁠, at strong coupling beyond the WDM regime (⁠$rs≥20$⁠), have recently been presented in Ref. 128.

In Secs. III A and III B, we have outlined the current state of the art regarding both, thermodynamics and the static density response of the UEG in the WDM regime. However, a direct comparison to experiments often requires the calculation of dynamic properties. For example, the central quantity in modern x-ray Thomson scattering (XRTS) experiments¹⁴⁹ is given by the dynamic structure factor $S(q,ω)$⁠, which is defined as the Fourier transform of the intermediate scattering function F(q, t) [cf. Eq. (8)],

Naturally, the straightforward evaluation of F(q, t) requires real time-dependent simulations,^174–176 for which, presently an exact simultaneous treatment of exchange–correlation, thermal, and degeneracy effects is not possible. Therefore, previous results^{155,174,177–181} were based on partly uncontrolled approximations, the quality of which had remained unclear. Moreover, ab initio QMC methods, which were pivotal for the accurate description of static properties, as was discussed in Secs. III A and III B, are effectively rendered unfeasible regarding time-dependent simulations due to an additional dynamical sign problem.^182,183

An alternative to simulations in real time is given by the method of analytic continuation, with the imaginary-time density–density correlation function $F(q,τ)$⁠, as defined in Eq. (8), being the starting point. Recall that the latter can be computed without any approximation from standard PIMC simulations, and it is connected to the dynamic structure factor via a Laplace transform,

The task at hand is then to find a trial solution $Strial(q,ω)$ that, when being inserted into Eq. (11), reproduces the PIMC data for $F(q,τ)$⁠, for all values of $τ∈[0,β]$⁠. Such an inverse Laplace transform is a well-known, but notoriously difficult task, as different $Strial(q,ω)$ might fulfill Eq. (11) within the given statistical uncertainty.^184,185 A first step to further restrict the space of possible $Strial(q,ω)$ are frequency moments, which are defined as

and the results for the cases $k=−1,0,1,3$ are known from different sum-rules.^170,186 While the combination of Eqs. (11) and (12) has often turned out to be sufficient to accurately reconstruct $S(q,ω)$ in the case of, e.g., ultracold bosonic atoms,^187–190 we have found that this does not hold for the UEG in the WDM regime, and additional information is indispensible.

To his end, we invoke the fluctuation–dissipation theorem⁹³ 

which gives a straightforward relation between the dynamic structure factor (DSF) and the imaginary part of the dynamic density response function introduced in Sec. III B. In this way, the original reconstruction problem has been recast into the quest for a suitable dynamic local field correction $Gtrial(q,ω)$⁠. This is extremely advantageous, as many additional exact constraints on G, such as the static and high-frequency limits, and the Kramers–Kronig relation between its real and imaginary part, are known.

In practice, we solve the inversion problem posed by Eq. (11) by stochastically sampling trial solutions $Gtrial(q,ω)$ such that a significant number of exact properties are built in by design. Subsequently, we use the corresponding $χtrial(q,ω)$ to generate trial solutions for the DSF via Eq. (13), which are finally plugged into Eqs. (11) and (12), and discarded if the deviation from our PIMC data is more than the Monte Carlo error bar. The final solution for $S(q,ω)$ is then computed as the average over a large number of such valid trial solutions, which also allows us to estimate the remaining variance around our estimate for $S(q,ω)$⁠.

This is illustrated in Fig. 6, where we show the frequency dependence of $S(q,ω)$ for a fixed wave number $q≈2qF$ at r_s = 6 and θ = 1. First and foremost, we note that the stochastic sampling of $Gtrial(q,ω)$ still allows for nontrivial structures in $S(q,ω)$⁠, and even solutions with two peaks are possible. The colored curves correspond to those $Strial(q,ω)$ that are not consistent with our PIMC data for $F(q,τ)$ and $⟨ωk⟩$⁠, whereas the black curves are valid solutions, and are included in the calculation of the average for the final solution for $S(q,ω)$⁠. Moreover, we note that all black curves fall within a narrow band around their average, and exhibit a single broad peak centered around twice the plasma frequency. Therefore, the remaining degree of uncertainty is small, and we have achieved an accurate reconstruction of the DSF of the warm dense UEG.

These new ab initio results for the dynamics of correlated electrons have opened up numerous new avenues for future research projects. Most importantly, we mention that the detailed investigation of $S(q,ω)$ is interesting in its own right and might, potentially, lead to the discovery of hitherto unobserved physical effects. This is illustrated in Fig. 7, where we show the dispersion relation of $S(q,ω)$ at the Fermi temperature for two different values of the density parameter r_s. Let us first focus on the red and green curves (shaded areas), which depict the peak positions (full width at half maximum) of the ab initio solution for $S(q,ω)$ computed from the stochastic sampling of the dynamic LFC (DLFC) and the random phase approximation (RPA). At the metallic density (r_s = 4, left panel), the exact DSF exhibits a significant red-shift and correlation induced broadening as compared to the mean-field solution, which are particularly pronounced for intermediate wave numbers. On the other hand, all curves converge for small and large q, as expected. The right panel shows the same information for stronger coupling strength, r_s = 10. Remarkably, we find a pronounced negative dispersion $ω(q)$ in the DLFC curve around $q≈1.8qF$⁠, which is not captured by the RPA at all. This feature had previously been reported by Takada and Yasuhara^180,181 based on the approximate results at zero temperature, and might indicate an incipient excitonic mode, which emerges in the electron liquid regime. A more detailed investigation of this effect, which includes a “phase diagram” of its appearance regarding r_s and θ, and a prediction of experimental conditions, where it can be measured, is currently in progress. As a second important application, we mention the interpretation of XRTS experiments,¹⁴⁹ where the DSF of the UEG is used to describe the free electronic part.¹⁵³ 

The third application of the new data for the dynamics of the warm dense UEG is their potential utility as input for other simulation methods. For example, the dynamic LFC is directly related to the exchange–correlation kernel of TD-DFT⁶⁹ via

which gives rise to the intriguing possibility to systematically go beyond the nearly ubiquitous adiabatic approximation for the XC-functional. In addition, such information can directly be used to further improve QHD simulations where a similar relation as (14) has been derived in Ref. 88, cf. Sec. V D.

The fourth application of our ab initio dynamic results is the computation of additional material properties, such as the stopping power,^143–145 the dynamic conductivity,^146,147 the dynamic dielectric function $ϵ(q,ω)$⁠, and the density response function $χ(q,ω)$⁠. As an example we show preliminary results for the dynamic density response function, Eq. (5) in Fig. 8. Again, we clearly see the effect of correlations, by comparing the dynamic results (red) to the RPA (green). While, for r_s = 2, the effect is relatively small and mainly seen in a redshift of the imaginary part; at r_s = 10, the RPA completely fails to describe the density response.

A fifth important application of the ab initio data for $S(q,ω)$ is that they unambiguously allow us to benchmark previous approximations,^{172,174,176,177} which are commonly used for WDM research. For example, Dornheim et al.¹⁸⁶ have reported that RPA exhibits significant inaccuracies even at relatively moderate coupling, r_s = 2 and θ = 1, where electronic correlation effects had often been assumed to play a minor role. This has potentially important consequences for the interpretation of WDM experiments, as the determination of plasma parameters, such as the electronic temperature and the degree of ionization, is sensitive to the exact dispersion relation of the DSF of free electrons.¹⁵³ Moreover, it has allowed us to introduce a significantly more accurate, yet computationally equally cheap alternative to the RPA. More specifically, we replace in Eq. (5) the dynamic LFC by its exact static limit,

that is conveniently available as a neural-net representation,¹⁷¹ see Sec. III B. The corresponding results for the dispersion relation of $S(q,ω)$ computed within this exact static approximation (as opposed to static dielectric theories like STLS,^101–103 where the results for $G(q,0)$ are approximate and systematically biased) are shown as the dashed black lines in Fig. 7. At warm dense matter conditions (r_s = 4 and θ = 1, left panel), the SLFC curve is basically indistinguishable from the exact data both with respect to the peak position and shape over the entire q-range. Upon approaching the electron liquid regime (r_s = 10, right panel), there do appear small yet significant deviations between the two curves, although the SLFC still captures both, the broadening and the negative dispersion. The same behavior is seen in the dynamic density response function, cf. Fig. 8. Thus, we conclude that the exact static approximation constitutes a distinct improvement over the RPA everywhere, without any additional effort.

Finally, an ambitious follow-up project regarding the ab initio calculation of dynamic properties is the extension of our simulations to real WDM systems, i.e., going beyond the UEG model and to also include ions. Although computationally challenging, this would allow for the first exact theoretical results for the dynamics of warm dense matter. More specifically, the combination of PIMC and the subsequent analytic continuation does not require any arbitrary external input, such as the XC-functional in DFT, or the Chihara decomposition,¹⁵⁴ which presupposes a potentially unrealistic distinction between bound and free electrons.⁶⁹ 

In this section, we explore the effect of the finite temperature exchange correlation functionals that were obtained by QMC simulations (cf. Sec. III A) in DFT simulations of dense plasmas. We present results for the equation of state (EOS) of dense hydrogen and carbon in Figs. 9 and 10.

The finite temperature DFT-MD method combines the quantum treatment of the fast moving electrons with the classical description of the slow ion dynamics.¹⁹¹ For the electrons, finite temperature DFT developed by Mermin⁹⁷ for the Kohn-Sham scheme¹⁹² is applied, which minimizes the grand potential, $Ω=E−TS−μN$⁠. Here, E is the total energy, T is the electron temperature, S is the entropy, μ is the chemical potential, and N is the number of electrons. For simplicity, E and S are expressed in the spin-averaged form as

and

where i is the index of the energy eigenvalues, E^H is the Hartree energy, E^xc is the exchange–correlation energy, V^ei is the ionic potential experienced by the electrons, and $n(r)$ is the charge density of the electrons. Furthermore,

represents the Fermi-Dirac equilibrium distribution. The energies ϵ_i, the wave functions ψ_i, the chemical potential μ, and the charge density $n(r)$ are self-consistently determined from the variational Kohn–Sham equation,

with

and μ is determined by the charge conservation equation

where the orthonormality of the orbitals has been assumed. When the Kohn–Sham–Mermin equations are solved self-consistently, the forces acting on each ion can be determined by the Hellman–Feynman theorem or its finite-temperature generalization. Then, classical Newton's equations are solved to compute the dynamics of the ions.

In the finite temperature DFT (FT-DFT) method, the many-body effects of the electrons beyond the Hartree mean field are accounted for by the exchange–correlation functional E^xc. In particular, if the exact functional would be used, one could reproduce the exact solution of the original many-body problem of interest. For practical applications, however, this term has to be approximated. Previous exchange–correlation functionals were limited to the case of zero temperature such as the expressions due to Perdew–Wang (PW)¹⁹³ and Perdew–Burke–Ernzerhof (PBE).¹⁹⁴ The latter are adequate for many condensed matter applications, but they become problematic when it comes to WDM. In this case, finite temperature and entropic effects in the exchange–correlation functional are becoming important,^75,76 and, instead of E^xc, an accurate exchange–correlation free energy, F^xc, has to be used. In this section, we quantitatively examine the importance of the related finite temperature effects.

The DFT-MD simulations for dense hydrogen have been performed using the CP2K code.¹⁹⁵ The Gaussian plane waves method is used to solve the Kohn–Sham equations with Gaussians as the basis set and additional plane waves as the auxiliary basis of the form $ψ(r)=Ri(r)Yli,mi(θ,ϕ)$⁠, with $Ri(r)$ denoting the radial part and $Yl,m$ denoting the angular part. Goedecker–Teter–Hutter (GTH) pseudopotentials of LDA (Padé) form are used for approximating the potential due to the usage of the LDA form of the xc-functional as the reference functional to compare with the parametrized LDA form incorporating finite temperatures, hereafter referred to as GDSMFB.^73,196 The GDSMFB functional is accessed in the CP2K code using the library of exchange–correlation (LIBXC) functionals commonly supported by DFT codes.^197,198

Due to the huge computational cost at high densities and the large temperatures considered in these cases, we choose 32 atoms in an hexagonal supercell for the simulations. On the same note, the sampling is performed only at the $Γ$-point. The system size and the sampling of k-points can improve the convergence of the EOS, especially near the liquid-liquid phase transition (LLPT) and in the low density limit.¹⁹⁹ We observe finite size effects resulting in a lower electronic pressure, at low temperatures, compared to the orbital-free MD results of Wang et al.²⁰⁰ 

At the low-temperature limit, especially near the LLPT (⁠$1000−2000$⁠) K, the pressure obtained using DFT-MD for these densities is dependent on a suite (Jacob's ladder, non-local, dispersion, $…$⁠) of xc-functionals ignoring other parameters such as system size and k-point sampling.^71,201 The first order phase transition (LLPT) is well characterized using $(∂P∂ρ)|T=0$ in the EOS visibly more prominent in the QMC results compared to DFT-MD which requires more sampling near the transition region.^202,203 The synthesis of metallic hydrogen is among the current key topics of high-pressure physics, and significant progress has been made over the last decade in the prediction of the transition using QMC and higher rungs of xc-functionals.^21,204–206 The inclusion of the finite-temperature component to the LDA can be simply ignored for characterizing the phase transition at temperatures (1000–5000) K, as this corresponds to $Θ<0.01$ and, instead, we focus on the improvement in the EOS results across a gamut of higher temperatures accessible using Kohn-Sham and orbital-free DFT.

The variation of electronic pressure with respect to temperature at two different densities for dense hydrogen is shown in Fig. 9. At $rs=1.4,$ the total pressure with the finite-temperature xc-functionals differs by less than 0.5% compared to LDA, in the temperature range (5000–10000) K, and accurate reptation quantum Monte Carlo (CEIMC) results show a deviation of less than 2%.²⁰⁷ With increasing temperature, the pressure obtained using GDSMFB converges toward the path integral Monte Carlo results obtained by Hu et al.^208,209 The relative difference in electronic pressure due to finite-temperature xc effects is more prominent at lower densities, with a maximum of 6%, observed at r_s = 1.4, for $Θ=0.33−0.42$⁠. The relative difference in total pressure is in good agreement with the Kohn-Sham DFT and orbital-free results obtained by Karasiev et al.⁷⁵ An analysis of finite temperature exchange–correlation effects on various optical and transport properties of deuterium was recently presented in Ref. 210, and an extensive topical investigation of hydrogen can be found in Ref. 76.

The DFT-MD simulations for dense carbon are performed using the recently developed ext-FPMD method²¹¹ implemented in the Quantum-Espresso code,²¹² which combines the analytical treatment of high-energy electrons as plane waves and the numerical treatment of the remaining electrons within the Kohn–Sham–Mermin scheme. This ext-FPMD method thus elevates the temperature limit of previous DFT-MD simulations and can be coherently applied from cold materials to hot dense plasmas.²¹³ The interaction between the carbon ions and the electrons is described by an all-electron projector augmented-wave (PAW) potential.²¹⁴ 32 carbon atoms are included in our simulation, which amounts to 192 electrons in total, in the cubic simulation box with periodic boundary conditions for all three directions. A shifted $2×2×2$ K mesh grid is used for all the simulations.²¹⁵ The ion temperature is controlled by an Andersen thermostat.²¹⁶ A sufficiently large number of time steps are applied to ensure that the system has reached equilibrium before data collection starts using the last 5000 time steps. The electronic part of pressure is converged to within 1% with respect to all parameters such as plane wave cutoff energy, K-mesh density, and finite size effects.

The variation of both, total pressure and electronic part of pressure of carbon, at a density of $ρ=10.0$ g/cm³, corresponding to $rs=1.4755$⁠, is shown in Fig. 10. For the lowest temperatures, i.e., at 1 and 5 eV, we find that LDA¹⁹³ and the finite temperature GDSMFB results⁷³ are close to each other. Both deviate from the PBE result,¹⁹⁴ which shows that the gradient correction of the exchange–correlation energy is more important than the finite-temperature effects, for a low temperature, as expected. As the temperature rises, LDA and PBE results get closer to each other, however, both deviate from the GDSMFB result. This shows that, in this region, finite-temperature effects play a more important role. The relative deviation for the electronic part of the pressure between zero-temperature exchange–correlation functionals and their finite-temperature counterpart reaches a maximum of ∼4% at 10⁵ K for PW and $2×105$ K for PBE, where $Θ=0.374$ and $Θ=0.749$⁠, respectively.

We note that a further increase in the temperature will eventually make the form of the exchange–correlation functionals less important, as the system approaches the hot classical plasma regime where many-body effects are less prominent. This is shown in Fig. 10 for the high-temperature region around 10⁷ K. PIMC results²¹⁷ are only available for high temperatures, and they are within 1% of our GDSMFB results shown in this figure. The deviations are comparable to the data accuracy, due to the statistical errors. OFMD simulations²¹⁸ struggle, in the low-temperature region, because they lack shell structure effects, and they are also found to be inaccurate in the high-temperature region, because they use the PBE exchange–correlation functional that does not account for finite-temperature effects.

The response of warm dense matter to an external excitation and the subsequent thermalization are of prime importance for many applications. This includes laser excitation and ionization of warm dense matter but also compression experiments including phase transitions and the path to inertial confinement fusion.

Theoretical methods for WDM out of equilibrium are even more challenging than equilibrium applications that were discussed above. They include time-dependent DFT,^69,219 semiclassical kinetics,⁷⁸ quantum kinetic theory and nonequilibrium Green functions,^{44,83,220,221} hydrodynamics^80,88,222 or rate equations.²²³ Among the problems that were studied are the equilibration of the electron distribution by electron–electron collisions,^83,220 non-thermal melting induced by fs x-ray pulses,²¹⁹ the density response for nonequilibrium momentum distributions,^79,221 density evolution following short-pulse laser excitation,⁷⁸ collisional heating of quantum plasmas by a laser pulse,^224,225 or ionization dynamics in a short laser pulse.⁸⁰ 

In the remainder of this section, we discuss the quantum hydrodynamics approach and its relation to DFT⁸⁹ more in detail because the former is comparatively less discussed for WDM applications, even though it appears to be filling a gap in the arsenal of simulation techniques, what we discuss in Sec. VI.

We consider a non-relativistic quantum system of N electrons described by the spin-independent hamiltonian

where $R=(r1,σ1;r2,σ2;…rN,σN), ri$ are the particle coordinates and σ_i their spin projections, and V is an external potential, e.g., due to the plasma ions. Assuming the first one as a pure state, the dynamics of the system are governed by the N-particle Schrödinger equation

that is supplemented by an initial condition and the normalization $∑σ1…σN∫d3NR |Ψ(R,t)|2=N$⁠. For particles with spin s, there are $gs=2s+1$ different spin projections, and each spin sum gives rise to a factor g_s (in the following, we will not write the spin arguments and spin sums explicitly).

If the many-body system (22) is coupled to the environment—as is typically the case in plasmas that we concentrate on in this section—a description in terms of wave functions and the Schrödinger equation (23) is no longer adequate. Instead, the system is described by an incoherent superposition of wave functions (“mixed state”). This can be taken into account, by replacing the N-particle wave function by the N-particle density operator,⁴⁰ 

where the sum runs over projection operators on all solutions of the Schrödinger equation (23), and p_a are real probabilities, $0≤pa≤1$⁠, with $∑apa=1$⁠. Here, we used a general representation-independent form of the quantum states. It is directly related to the wave functions if the coordinate representation is being applied: $⟨R|Ψa(t)⟩=Ψa(R,t)$ [$⟨R|$ are eigenstates of the coordinate operator in N-particle Hilbert space]. The previous case of a pure state is naturally included in definition (24) by setting p_k = 1 and all $pa≠k=0$⁠. The second relation (24) is the normalization condition where the trace denotes the sum over the diagonal matrix elements of $ρ̂$⁠, see below. The equation of motion of $ρ̂$ follows from the Schrödinger equation (23) and is the von Neumann equation supplemented by the initial condition,

The method of density operators is well established in quantum many-body and kinetic theory, and relevant representations are the coordinate representation, momentum, and Wigner representation, e.g., Refs. 36, 40, and 226. From the N-particle density operator, all time-dependent properties of a quantum system can be obtained. However, in many cases, simpler quantities are sufficient such as reduced s-particle density operators, including the single-particle density operator (which is related to the distribution function or the Wigner function),⁴⁰ 

The equation of motion for $F̂1$ follows straightforwardly from Eq. (25), e.g., Refs. 40 and 89, and is given below, cf. Eq. (29).

We now derive the equation of motion for the time-dependent single-particle orbitals of N interacting electrons. We follow the idea of DFT that the many-particle quantities are expressed in terms of single-particle quantities via a mean field description. Exchange and correlation effects are then taken into account a posteriori, by adding the potential $Vxc$⁠.

Considering an N-particle system in the grand canonical ensemble (specified by the inverse temperature β and chemical potential μ), the single-particle nonequilibrium density operator has the following form:

where the mean occupation numbers in equilibrium are given by the Fermi function (18). The equation for $F̂1$ in the mean field (Hartree) approximation has the following form:⁴⁰ 

Correlation effects would give rise to a collision integral on the r.h.s. of Eq. (29), for various approximations, see Ref. 40.

From Eq. (28), we obtain the density matrix by multiplying with coordinate eigenstates $⟨r′|$ and $|r″⟩$⁠:

The single-particle wave functions are the so-called “natural orbitals,” and in the mean field approximation, the N-particle wave function obeying Eq. (23) is just their product. By inserting Ansatz (31) into the coordinate representation of Eq. (29), it is easy to verify that the latter is solved when each orbital fulfills the following single-particle Schrödinger equation (⁠$i=1,…,N$⁠)

where the Hartree mean field is the coordinate representation of operator (30) and contains the densities of all occupied orbitals. Equations (32) and (33) are the time-dependent Hartree equations for weakly interacting fermions (interactions are taken into account only via the mean field $UH$⁠).

This result can be directly extended beyond the mean field approximation by replacing

and, as a consequence, Eqs. (32) and (33) become the time-dependent Kohn-Sham equations–the basic equations of the time-dependent density functional theory (TD-DFT).²²⁷ A particular strength of this theory is its solid theoretical foundation on the Runge-Gross theorem²²⁷ and the corresponding theorems for time-independent DFT.²²⁸ The basic statement is that a system of N interacting fermions can be mapped exactly onto a system of N non-interacting particles with the same density $n(r,t)$ where all interactions are lumped into an effective single-particle potential that is a direct generalization of the Hartree potential (33).

The first remarkable property of these equations is that, both, the mean field and the additional exchange–correlation potential do not explicitly depend on the individual orbital wave functions but only on the total density, so also the coordinate dependence is only implicit, via the functional $n(r,t)$⁠. However, the exact functional $Vxc$ does not only depend on the current density, $n(r,t)$⁠, but, in general, the dependence is also on the density profile at earlier times, $n(r,t̃), 0≤t̃≤t$⁠. At the same time, most current implementations neglect this “memory” effect and use an adiabatic approximation (e.g., adiabatic LDA, ALDA), $t̃→t$ which leads to systematic errors. In Eq. (35), we also indicated that, at a finite temperature, $Vxc$ carries a temperature dependence which was discussed in detail before, see Sec. III.

Following Ref. 89, we now derive the microscopic QHD (MQHD) equations, starting from the time-dependent Hartree equations (32) and (33). To this end, we simply convert each orbital solution, $ϕi(r,t)=Ai(r,t)eiℏSi(r,t)$⁠, into an individual pair of amplitude and phase equations,^86,89 for $i=1,…,N$⁠,

where $ni=Ai2$ and $pi=∇Si$⁠, and we introduced a short notation for the total potential energy, $Utot=V+UH$⁠. This system of MQHD equations is fully equivalent to TD-DFT and, for $Vxc→0$⁠, it exactly coincides with the time-dependent nonlinear Hartree (quantum Vlasov) equations.^40,229 Moreover, it was shown in Ref. 89 that this approximation, in linear response, is exactly equivalent to the random phase approximation (RPA or the linearized quantum Vlasov equation). In particular, the linearized MQHD equations then yield the correct plasmon spectrum and the correct screening of a test charge – in contrast to the standard QHD (see below).

To convert these microscopic equations into a single pair of density and momentum equations (QHD), a suitable averaging over the orbitals is necessary which we denote by a “bar,”

where the orbitals are weighted by the Fermi function. Here, $n¯$ is interpreted as the mean orbital probability density. In Ref. 85, the authors assumed that all orbital amplitudes are equal whereas, in Ref. 86, it was assumed that one can substitute

Finally, in Ref. 89, it was demonstrated how the QHD equations can be derived without uncontrolled assumptions, and here we briefly recall that approach.

In order to take the orbital average of the MQHD equations (36)–(38), we express each of the orbital quantities in terms of their averages and fluctuations,

and take into account that the average of products of two orbital quantities is given by

where, in addition to the product of averages, there appears a correlation function. Note that the averaging of the Bohm potentials, Eq. (41), has to be done with care because the orbital depending densities enter at two places and we have to apply relation (43), with the result

Here, the first term is just the Bohm potential from the single-particle case with the density replaced by the mean density, $ni→n¯$ (this is what was used in Refs. 85 and 86) corresponding to Eq. (42), and the second term is the deviation which is presented below, in Eq. (48). With this, we can perform the averaging of the MQHD equations (36)–(38), and obtain the QHD equations that contain three correlation functions that we denote by J_np, J_pp, and $QΔ$⁠,

The function J_pp contains the correlations of the fluctuations of the momentum field.

Another formulation of the QHD equations, that is closer to classical hydrodynamics, is obtained if, instead of the mean orbital density $n¯$ and the mean momentum $p¯$⁠, we consider the density $n(r,t)$ and the current density $j(r,t)$⁠, defined in terms of the orbital quantities n_i and $ji=ni vi$⁠, as^86,230

where a factor 2 was included to account for the two electron spin projections. Using the definitions for n and $j$⁠, we can further define a mean velocity field, $u(r,t)=j(r,t)/n(r,t)$ which differs from the mean velocity $v¯=p¯/m$ that was defined above, cf. Eq. (40).

The dynamical equations for n and $j$ follow from Eqs. (36)–(38). They read as (see also Ref. 230)

where $Π(r,t)$ appears as a “pressure” tensor,

The first contribution arises from fluctuations of the orbital velocity fields, while the second term is due to the microscopic Bohm potential, i.e., due to the curvature of the orbital amplitudes. Note that the correlation function J_np appearing in Eq. (45) is contained in the definition of the current density and does not appear explicitly in Eq. (51). Analogous equations have been discussed in Ref. 230, see also Refs. 231 and 232.

An important test for the QHD and MQHD models is the result for electron plasma oscillations (Langmuir waves) in the limit of the weak external field (linear response). Here, we consider the simplest case of a spatially homogeneous weakly non-ideal electron gas (interactions are included only via the Hartree mean field, whereas exchange–correlation effects are neglected) at zero temperature, i.e., the statistical weights $fieq$ reduce to unity, for $ϵi≤EF$⁠, and zero otherwise.

Considering first the MQHD equations, Eqs. (36)–(38), we apply a harmonic monochromatic excitation, $ϕ1(r,t)∼e−iωt+iq·r$⁠, and linearize $ni(t)$ and $pi(t)$ around the unperturbed solution. Finally, the density response is computed via orbital averaging and Fourier transformation, with the result given by Eq. (55). Second, we consider the QHD equations (51) and (52), which already contain the orbital averaging. Here, in order to make further progress, an approximation for $Π$ is required. If we approximate the tensor by the Fermi pressure and the Bohm term [with the density $n(r,t)$],

linearization of these equations and Fourier transformation yields the plasmon dispersion for the case of D = 1, 2, 3 dimensions,²³³ given by Eq. (56).

where we introduced the Fermi velocity via $EF=mvF2/2$⁠. Note that the plasma frequency depends on the dimensionality of the system.

Let us discuss this result. First, we observe that the MQHD-dispersion (55) exactly coincides with the zero temperature limit of the RPA result presented in Fig. 7. In particular, the small-q (large-q) limit given by the first (third) term in Eq. (55) is correct. However, comparing with the QMC results (red curves in Fig. 7), the MQHD dispersion at intermediate wave numbers, $1 ≲ q/qF ≲ 3$⁠, strongly overestimates the oscillation frequency. This is due to the neglect of correlation effects in the present calculation. These effects can be restored by including a proper expression for the exchange–correlation potential in the MQHD equations.

Due to the agreement of the MQHD dispersion with the RPA, we can use it to test the accuracy of the QHD model. Note that the QHD model is—by construction—less accurate than MQHD because it involves an orbital averaging and thus, a loss of resolution of small length and energy scales. First, we observe that the small-q limit is correct. The large-q limit, on the other hand, is correct for 1D and 3D systems, but it is incorrect for 2D quantum plasmas. Third, the behavior at intermediate wave numbers (second term proportional to q²) is correct in 1D. In 2D, the QHD yields a coefficient 1/2, whereas the correct one is 3/4. Similarly, in 3D, the QHD result (1/3) deviates from the MQHD result (3/5). Therefore, high-frequency electronic plasma oscillations (⁠$ω≥ωpl$⁠) are correctly reproduced only in one-dimensional quantum plasmas.

Going back to the QHD equations, it is easy to verify that the origin of this incorrect coefficient, as compared to the RPA (MQHD) result, is the Fermi pressure term. To recover the correct coefficient of the q² term in the dispersion (56), the pressure has to be multiplied by a factor $α¯$ that was reported in Ref. 88: for example, for a 3D plasma, the Fermi pressure that appears in Eq. (54) has to be multiplied by a factor $α¯=9/5$⁠. Note that this value applies only to small wave numbers, whereas for large wavenumbers, $q≫2kF$⁠, it will approach the value 3/5.⁸⁸ 

After having discussed the high-frequency plasma oscillations, we now turn to the frequency range $ω≪ωpl$⁠. This is of prime importance, e.g., for ion acoustic oscillations and for the screened potential in a quantum plasma. In fact, screening effects replace the Coulomb potential Q/r by the potential

where the dielectric function contains the quantum plasma properties and is taken in the static limit. The screened potential for quantum plasmas that improves the conventional Yukawa (or Thomas-Fermi) model has been actively studied in recent years, e.g., Refs. 135, 234–237. For example, in Ref. 234 the authors predicted, using a QHD model that, in a quantum degenerate plasma in thermodynamic equilibrium, the electrostatic potential of an ion would be attractive. Their result is shown in Fig. 11 by the blue line and exhibits a shallow minimum (depth approximately 5 meV) at about 6 Bohr radii (about 2.5 interparticle distances). Tests with DFT simulations that can be regarded as benchmarks^235,238,239 revealed that no such minimum exists. The reason for the unphysical predictions of the QHD model was clarified in Ref. 135: the coefficient of the Bohm term in the QHD equations (45) and (46) turns out to be incorrect for applications to low-frequency excitations. Compared to its value at high frequencies, $ω≥ωpl$ (which we use as a reference, γ = 1), it has to be reduced by almost an order of magnitude (⁠$γ→1/9$⁠) to reproduce the correct MQHD (RPA) result. Two screened potentials that contain the correct factor $γ=1/9$ are included in the plot as well and do not exhibit the minimum and show good agreement with the full (nonlocal) RPA screened potential.¹³⁵ The only difference is that these potentials cannot resolve the Friedel oscillations. Note that the correction factor γ depends on temperature which is shown in Fig. 12.

Let us now investigate the effect of electronic correlations (effects beyond RPA) on the screened potential. Here, we consider two approximations for the static local field correction G(q) that lead to a correlated dielectric function appearing in formula (57). The first is the standard STLS (Singwi–Tosi–Land–Sjölander) approximation and the second, the exact results for G(q) obtained from our QMC simulations, cf. Sec. III B, where we used the machine learning representation of Ref. 171. The two results are shown, together with the uncorrelated RPA result, in Fig. 13. There, we present the screened potential of a proton in a dense quantum plasma in the WDM regime, at r_s = 2 and $θ=0.5$ which is close to the parameters of Fig. 11, the main difference being the finite temperature. The first effect of correlations is a significantly stronger screening (more rapid decay of the potential), compared to the RPA case.¹³⁶ In addition, we observe that the STLS potential strongly deviates from the QMC result not only quantitatively but even qualitatively: it overestimates screening and develops an unphysical attractive part (negative minimum) at intermediate distances (see the subplot in Fig. 13).¹³⁶ 

As was explored in Ref. 3, this unphysical behavior of the STLS approximation leads to additional restrictions on its applicability for two-component plasmas with non-ideal ions. Note that, in general, the screening is not exponential³ and, in the case of strongly coupled ions, deviations from the RPA screened potential may be quite large and have a significant impact on the structural and dynamical properties of the ion component in a dense plasma.^3,240 Given the overall good accuracy of STLS for thermodynamic quantities, cf. Sec. III A, these problems for the screened potential are an unexpected result. This is a similar artifact, as was observed for the QHD screened potential with the wrong Bohm term in Fig. 11. This observation underlines the high importance of ab initio QMC results, in particular for the local field correction. The QMC potential presented in Fig. 13 is a preliminary result, and a systematic analysis for a broad range of parameters in the WDM range is an important task for future work.

Let us now turn to ion-acoustic oscillations in a quantum plasma that have been studied by many authors. For example, Schmidt et al. gave a hydrodynamic derivation of the dispersion relation and the dynamic structure factor.²⁴¹ Haas and Mahmood²²² used Euler fluid equations for the ions, neglecting the ionic pressure term, and the QHD equations for the electrons with the Bohm potential with the correct factor γ, in the low frequency long-wavelength limit (see Fig. 12), to study low-frequency waves in a two-component plasma. They derived the following dispersion for ion-acoustic waves:

Here, $A(γ,cs)=γℏ2/(12memics2)$ represents the quantum electron correction due to the Bohm term, and $cs=(kBT/mi Li3/2(z)/Li1/2(z))1/2$ is the ion sound speed of an ideal quantum plasma written in terms of the polylogarithmic function $Liν(z)$ of the ideal electron fugacity z. In the limit of strong electron degeneracy, the sound speed is given by $cs=2EF/(3mi)$⁠.

The dispersion of the ion-acoustic wave is sensitive to the properties of the surrounding electrons via screening of the ion potential, Eq. (57), and to ionic correlations, depending on the coupling parameter $Γi$⁠. This is illustrated in Fig. 14 where the QHD result (58) is compared to the data from the molecular dynamics (MD) simulation of ions at $Γi=15, rs=1.5$⁠, and $θ=0.1$⁠. The MD data was obtained using the screened ion potential (57)with $ϵ(q,ω=0)$ in RPA.²⁴⁰ As discussed above, in linear response, the RPA (MQHD) description of electrons is more accurate than QHD with the standard Bohm potential and serves as a benchmark. The first observation from Fig. 14 is that the QHD result for the dispersion, Eq. (58), strongly deviates from the MD data for all wave numbers, because the sound speed of an ideal plasma is being used. A much improved behavior is found if the sound speed fitted to the MD data at small wave-numbers is inserted into Eq. (58). The resulting dispersion (black dashed curve) agrees well with the MD data for all wavenumbers, $ka ≲ 1.0$⁠. The failure at a larger wavenumber is not surprising because the standard limitation of the fluid approaches that we will discuss in Sec. V H. To summarize, our analysis reveals that the functional form (58) of the dispersion is reasonable which is promising for further applications of the QHD approach to two-component plasmas, taking into account electronic and ionic non-ideality effects.

The general validity limits of quantum hydrodynamic models have been discussed in a variety of papers, e.g., Refs. 86, 235, and 242. Most importantly, similar to a classical plasma, a sufficient condition for the applicability of a hydrodynamic description is that the considered length scales are larger than the screening length, e.g., Ref. 243. In a quantum plasma. this is the Thomas-Fermi length, λ_TF, which leads to the criterion

In the weak coupling limit, $Γq≤0.1$⁠, this leads to the restriction $k ≲ 0.03kF$⁠. Previously, it was found that the above condition is not always necessary, i.e., the hydrodynamics in many cases applies also to length scales smaller than the screening length, for classical plasmas, see, e.g., Ref. 243. For correlated quantum plasmas, it was found that correlation effects reduce the screening length,^3,136,240 cf. also Fig. 13, which supports the same conclusion that QHD should be valid on scales below λ_TF.

Of course, a further extension of the validity range can be achieved if the coefficients in the QHD equations are adjusted using information from the MQHD equations (kinetic theory). For example, it was shown in Ref. 88 how to properly choose the QHD parameters for a large wavenumber, $k≳2kF$⁠. Moreover, tests against the MQHD (RPA) results, discussed in Sec. V E, allow us to correct the prefactors in the QHD equations for the high-frequency and low-frequency limits.

Thus, even if the correlation effects are neglected (⁠$Vxc→0,Jnp→0, QΔ→0$⁠), it is clear that there exists no universal result for the parameters in the QHD equations that would apply to arbitrary situations. Instead, depending on the frequency and wave number of the excitation, the coefficients in front of the Fermi pressure and of the Bohm term vary, and the results are also dependent on the temperature and the system dimensionality:

The values of $α¯$ and γ for the important limiting cases of high and low frequency, as well as high and low wave number are known analytically, even at a finite temperature.⁸⁸ Thus for these situations, reliable simulations are possible. The reason for this frequency and wave number dependence of the coefficients is the fact that the kinetic theory (RPA) polarization (density response) function $ΠRPA$ has different long wavelength limits for different frequencies, $limq→0ΠRPA(q,ω)$⁠. When this is converted into local hydrodynamic equations via orbital averaging, the result is different for high and low frequencies, respectively.