Ab initio simulation of warm dense matter

Warm dense matter (WDM) -- an exotic state of highly compressed matter -- has attracted high interest in recent years in astrophysics and for dense laboratory systems. At the same time, this state is extremely difficult to treat theoretically. This is due to the simultaneous appearance of quantum degeneracy, Coulomb correlations and thermal effects, as well as the overlap of plasma and condensed phases. Recent breakthroughs are due to the successful application of density functional theory (DFT) methods which, however, often lack the necessary accuracy and predictive capability for WDM applications. The situation has changed with the availability of the first \textit{ab initio} data for the exchange-correlation free energy of the warm dense uniform electron gas (UEG) that were obtained by quantum Monte Carlo (QMC) simulations, for recent reviews, see Dornheim \textit{et al.}, Phys. Plasmas \textbf{24}, 056303 (2017) and Phys. Rep. \textbf{744}, 1-86 (2018). In the present article we review recent further progress in QMC simulations of the warm dense UEG: namely, \textit{ab initio} results for the static local field correction $G(q)$ and for the dynamic structure factor $S(q,\omega)$. These data are of key relevance for the comparison with x-ray scattering experiments at free electron laser facilities and for the improvement of theoretical models. In the second part of this paper we discuss simulations of WDM out of equilibrium. The theoretical approaches include Born-Oppenheimer molecular dynamics, quantum kinetic theory, time-dependent DFT and hydrodynamics. Here we analyze strengths and limitations of these methods and argue that progress in WDM simulations will require a suitable combination of all methods. A particular role might be played by quantum hydrodynamics, and we concentrate on problems, recent progress, and possible improvements of this method.


I. INTRODUCTION
Warm dense matter has become a mature research field on the boarder of plasma physics and condensed matter physics, e.g. [1][2][3][4]. There are many examples in astrophysics such as the plasma-like matter in brown and white dwarf stars [5][6][7], giant planets, e.g. [8][9][10][11][12][13] and the outer crust of neutron stars [14,15]. Warm dense matter is also thought to exist in the interior of our Earth [16]. 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 Lawrence Livermore National Laboratory [18,19], the Z-machine at Sandia National Laboratory [20,21], the Omega laser at the University of Rochester [22], 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][19][20] where electronic quantum effects are important during the initial phase. Aside from dense plasmas, also 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, among others, by electronic quantum effects, moderate to strong Coulomb correlations and finite temperature effects. Quantum effects of electrons are of relevance at low temperature and/or if matter is very highly compressed, such that the temperature is of the order of (or lower than) the Fermi temperature (for the relevant parameter range, see Fig. 1 and, for the parameter definitions, see Sec. II).
A recent breakthrough occurred with the application of Kohn-Sham density functional theory (DFT) simulations because they, for the first time, enabled the selfconsistent simulation of realistic warm dense matter, that includes both, plasma and condensed matter phases, e.g. [60][61][62]. Further developments include orbital-free DFT methods (OF-DFT) and time-dependent DFT (TD-DFT), e.g. [63] 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 [4] or against electron-ion quantum Monte Carlo [64][65][66]. Also, the use of finite-temperature functionals was shown to be important [67,68] when the XCcontribution is comparable to the thermal energy, see Ref. [69] for a topical discussion. One goal of this paper is to present an overview of these results and discuss future research directions.
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. A possible candidate 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 [78,79]. 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 Ref. [80,81], and a systematic derivation of the QHD equations from the time-dependent Kohn-Sham equations was given in Ref. [82].
The goal of this paper is to present a summary of some of the recent in 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 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.

II. WARM DENSE MATTER PARAMETERS
Let us recall the basic parameters of warm dense matter [40,82]: the first are the electron degeneracy parameters θ = k B T /E F and χ = nΛ 3 , where Λ = h/ √ 2πmk B T , is the thermal DeBroglie wave length, and the Fermi energy of electrons (in 3D) is where n is the electron density and T the electron temperature. The ion degeneracy parameter is a factor (mT /m i T i ) 3/2 smaller than the one of the electrons and typically negligible for WDM. Second is the classical coupling parameter of ions Γ i = Q 2 i /(a i k B T i ), where Q i is the ion charge, and a i is the mean inter-ionic distance. Further, 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 m r = mm i /(m + m i ) and b are the reduced mass and background dielectric constant, respectively, for hydrogen m r ≈ m , b = 1, and a B = 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/n tot , Figure 1. Density-temperature plain with examples of plasmas and characteristic plasma parameters. ICF and WDM denote inertial confinement fusion and warm dense matter, respectively. Metals (semiconductors) refers to the electron gas in metals (electron-hole plasma in semiconductors). Weak electronic coupling is found outside the line Γ eff = 0.1, cf. Eq. (4). Electronic (ionic) quantum effects are observed to the right of the line χ = 1 (χp = 1). The coupling strength of quantum electrons increases with rs (with decreasing density). Atomic ionization due to thermal effects (due to pressure ionization) is dominant above (to the right of) the red line, for the case of an equilibrium hydrogen plasma [83]). The values of χp and rs refer to the case of hydrogen. Figure modified from Ref. [82].
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 density whereas the quantum coupling parameter decreases with the density n. We underline that the parameters a B , E F , Λ, θ contain the density of free electrons and Γ i , χ i the density of free ions. This means the lines of constant Γ i , χ, χ i , r s 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 −|E b |/k B T , 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.  The uniform electron gas (UEG) is one of the most fundamental model systems in physics [4,84,85]. In particular, the accurate parametrization [86,87] of the ground-state exchange-correlation energy e xc (r s ), based on ab inito quantum Monte Carlo simulations [88], has been essential for the striking success of 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 [89,90] simulation of warm dense matter, see Sec. IV, requires a parametrization of the exchange-correlation free energy f xc (r s , θ), which explicitly depends both on density and on the temperature [69,91].
While being an important mile stone, these data have been obtained on the basis of the uncontrolled fixed node approximation [108], which has recently been revealed to be surprisingly inaccurate with systematic errors in the exchange-correlation energy exceeding 10% at high density and low temperature [54].
This unsatisfactory situation has sparked a surge of new developments in the field of fermionic QMC simulations at finite temperature [55,[109][110][111][112][113][114][115][116]. In particular, Groth, Dornheim, and co-workers have introduced a combination of two complementary QMC methodspermutation blocking PIMC (PB-PIMC) and configuration PIMC (CPIMC)-that allows 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 [117], the same authors presented the first ab initio parametrization of f xc 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. [67] (GDB), 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. [103] (BCDC, blue circles) and a corresponding parametrization (KSDT, solid blue), HNC data and a parametrization thereof by Tanaka [96] (HNC, green squares and line), and the STLS-based parametrization by Ichimaru et al. [93]. At high density (r s = 0.1, panel a), both dielectric theories agree relatively well with the GDB 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,117,118]. 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 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 elevated temperature.
In the mean time, the availability of the accurate GDB benchmark data has led to a revised version of the KSDT parametrization (denoted as corrKSDT in Refs. [68,119]), which basically reproduces GDB over the entire WDM regime [68]. First and foremost, we note that both GDB and corrKSDT are suitable to be used as an exchange-correlation functional on the level of the local density approximation [the GDB 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 [119]. 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. [68] have found that there occur some potentially unphysical oscillations in quantities, that are derived from f xc , 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 GDB and corrKSDT.
Let us conclude this section by proposing a few possible solutions to the remaining open questions: i) the unphysical behaviour in the second-and higher-order derivatives are most likely a consequence of the functional form of the GDB and corrKSDT parametrizations [68]. 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 GDB parametrization to 0 ≤ r s ≤ 20 can be significantly extended by incorporating the recent ab initio PIMC results for the electron liquid regime (20 ≤ r s ≤ 100) by Dornheim et al. [120]. iii) new ab intio QMC results at low temperature, θ ∼ 10 −1 , could help to more accurately resolve open thermodynamic questions like the effective mass enhancement [121], but are difficult to obtain due to the notorious fermion sign problem [122]. 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.
B. Summary of ab initio results for the static local field correction One important step in going beyond local approximations, such as LDA or 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 = 4e 2 π/q 2 [frequently atomic units are used, then this becomes 4π/q 2 ] is the Fourier transform of the Coulomb potential. Further, χ 0 (q, ω) denotes the usual Lindhard function that describes the density response of the ideal (i.e., noninter-  [93]. The bottom panel shows the relative deviation in v with respect to the GDB reference data. Taken from Ref. [4] with the permission of the authors. acting system) [85], and the local field correction G(q, ω) entails the full frequency-and wave number resolved information about exchange-correlation effects on χ [123]. 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.
Naturally, there have been many attempts to find suitable approximations for G(q, ω), most commonly within the purview of dielectric theories [147][148][149][150][151]. The first accurate benchmark data for the LFC have been obtained by Moroni et al. [152,153] 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 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 input for the parametrization of G(q, 0) by Corradini et al. [154] (CDOP), which, in turn, has been used for many applications, e.g., Refs. [3,128,131,132,155,156]. Recently, Dornheim, Groth and co-workers [157,158] 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 at WDM conditions. While being conceptually valid and interesting in their own right, these simulations suffer from a prohibitive computational cost: the full characterisation of G(q, 0; r s , θ) requires a dense grid of densities, temperatures, and wave num-bers, (r s , θ, 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 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; r s , θ) covering the entire WDM regime.
A different, more convenient route is given by the imaginary-time version of the fluctuation-dissipation theorem [157,159], 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 [160,161], which means that the entire wavenumber 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. [162], This strategy-in combination with the efficient finitesize correction introduced in Ref. [158]-was recently used by Dornheim et al. [163] to obtain extensive new ab initio PIMC results for G(q, 0) for N param ∼ 50 densitytemperature 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, r s , θ) and predicts as output the corresponding LFC G(q, 0; r s , θ) in the range of 0 ≤ θ ≤ 4, 0.7 ≤ r s ≤ 20, and 0 ≤ q ≤ 5q F . 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. [153], and the dashed blue line the corresponding CDOP parametrization, which incorporates both, the compressibility sum-rule [95], for q → 0, and the exact large-q limit found by Holas [164,165]. 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. [163]. Evidently, the machine-learning representation nicely reproduces the input data where they are available, and . The static local field correction G(q, 0; rs, θ) in the wave number-temperature plain at rs = 5. The black squares and dashed blue line depicts the ground-state QMC data from Ref. [153] (MCS) and an accurate parametrization thereof [154] (CDOP), and the red crosses correspond to the new finite-temperature PIMC data from Ref. [163]. The green surface shows the prediction by the neural net, that is available at continuous q, rs, and θ. Taken from Ref. [163] with the permission of the authors.
smoothly interpolates in between. Moreover, its capability to predict G(q, 0; r s , θ) has been validated against independent benchmark data that had not been included in the training procedure, see Ref. [163] 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 3q F , but exhibits a significant drop and an apparent saddle point, for large wave numbers. At the Fermi temperature (blue dashes), G(q, 0; r s , θ) exhibits an even more interesting behavior: while it is approximately equal to the θ = 0 curve, for q 2q F , there appears a complicated shape with a maximum around q ≈ 2.7q F , a subsequent minimum at q ≈ 4.5q F , 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 [163].
For completeness, we mention that new ab initio results for G(q, 0; r s , θ) at strong coupling beyond the WDM regime (r s ≥ 20) have recently been presented . Wave-number dependence of the static LFC at metallic density, rs = 4. All curves have been obtained by using the machine-learning representation of G(q, 0; rs, θ) from Ref. [163], which approaches CDOP [154] in the lowtemperature limit.

C. Ab initio dynamic results
In the 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 [141] is given by the dynamic structure factor S(q, ω), which is defined as the Fourier transform of the intermediate scattering function Naturally, the straightforward evaluation of F (q, t) requires real time-dependent simulations [166][167][168], for which, presently an exact simultaneous treatment of exchange-correlation, thermal, and degeneracy effects is not possible. Therefore, previous results [147,166,[169][170][171][172][173] 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 the Sections before, are effectively rendered unfeasible regarding time-dependent simulations due to an additional dynamical sign problem [174,175]. An alternative to simulations in real time is given by the method of analytic continuation, with the imaginarytime 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 S trial (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 S trial (q, ω) might fulfill Eq. (11) within the given statistical uncertainty [176,177]. A first step to further restrict the space of possible S trial (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 [162,178]. 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 [179][180][181][182], 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 [85] which gives a straightforward relation between the 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 G trial (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 G trial (q, ω) such that a significant number of exact properties are build 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 ≈ 2q F at r s = 6 and θ = 1. First and foremost, we note that the stochastic sampling of G trial (q, ω) still allows for nontrivial structures in S(q, ω), and even solutions with two peaks are possible. The colored curves correspond to those S trial (q, ω) that are not consistent with our PIMC Figure 6. Stochastic sampling method for the dynamic structure factor. Shown is the frequency dependence of S(q, ω) for the unpolarized UEG at rs = 6 and θ = 1 for a fixed wave number q ≈ 2qF. The black curves depict trial solutions, S trial (q, ω), that are consistent with our PIMC data [cf. Eqs. (11) and (12)], whereas the colored curves are being discarded. Taken from Ref. [162] with the permission of the authors.
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 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.8q F , which is not captured by RPA at all. This feature had previously been reported by Takada  15) (static LFC, SLFC), and the usual random phase approximation (RPA), respectively. The corresponding shaded areas illustrate the full width at half maximum, and the grey areas enclosed by two parabolas the zero temperature electron-hole pair continuum [85]. Taken from Ref. [178] with the permission of the authors.
and Yasuhara [172,173] based on 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 [141], where the DSF of the UEG is used to describe the free electronic part [145].
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 [63] 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. [81], 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 [135][136][137], the dynamic conductivity [138,139], 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 [164,166,168,169], which are commonly used for WDM research. For example, Dornheim et al. [178] 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 [145]. 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 [163], 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 [93][94][95], 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 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 XCfunctional in DFT, or the Chihara decomposition [146], which presupposes a potentially unrealistic distinction between bound and free electrons [63].

A. Kohn-Sham-Mermin DFT
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 [183]. For the electrons, finite temperature DFT developed by Mermin [89] for the Kohn-Sham scheme [184] is applied, which minimizes the grand potential, 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 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, n(r) is the charge density of the electrons. Further, 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 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 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) [185] and Perdew-Burke-Ernzerhof (PBE) [186]. 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 [69], and, instead of E xc , an accurate exchangecorrelation free energy, F xc , has to be used. In this section, we quantitatively examine the importance of the related finite temperature effects.

B. Equation of state of warm dense hydrogen
The DFT-MD simulations for dense hydrogen have been performed using the CP2K code [187]. 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) = R i (r)Y li,mi (θ, φ), with R i (r) denoting the radial part and Y l,m denoting the angular part. Goedecker-Teter-Hutter pseudopotentials (GTH) of LDA (Pade) 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 [67,188]. The GDSMFB functional is accessed in the CP2K code using the library of exchange-correlation functionals (LIBXC) commonly supported by DFT codes [189,190].
Due to the huge computational cost at high densities and the large temperatures considered in these cases, we choose 256 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 [191]. For N < 256, 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. [192].
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 [65,193]. The first order phase transition (LLPT) is well characterized using ( ∂P ∂V )| T > 0 in the EOS visibly more prominent in the QMC results compared to DFT-MD which requires more sampling near the transition region [194]. 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,[195][196][197]. 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 r s = 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% [198]. With increasing temperature, the pressure obtained using GDSMFB converges towards the path integral Monte Carlo results obtained by Hu et al. [199]. The relative dif-ference 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. [69].

C. Equation of state of warm dense carbon
The DFT-MD simulations for dense carbon are performed using the recently developed ext-FPMD method [200] implemented in the Quantum-Espresso code [201], which combines the analytical treatment of high-energy electrons as plane waves and the numerical treatment of the remaining electrons within 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 [202]. The interaction between the carbon ions and the electrons is described by an all-electron PAW potential [203]. 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 [204]. The ion temperature is controlled by an Andersen thermostat [205]. 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 3 , corresponding to r s = 1.4755, is shown in Fig. 10. For the lowest temperatures, i.e. at 1eV and 5eV, we find that LDA [185] and the finite temperature GDSMFB results [67] are close to each other. Both deviate from the PBE result [186], which shows that the gradient correction of the exchange-correlation energy is more important than the finite-temperature effects, for 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 importance role. The relative deviation for the electronic part of the pressure between zero-temperature exchangecorrelation functionals and their finite-temperature counterpart reach a maximum of ∼4% at 10 5 K for PW and 2 × 10 5 K for PBE, where Θ = 0.374 and Θ = 0.749, respectively.
We note that a further increase of 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 hightemperature region around 10 7 K. PIMC results [206] are PW: zero-temperature LDA (PW-functional) [185]; PBE: zero-temperature PBE-functional [186]); Danel-OF: orbitalfree MD results [207]; Benedict-PIMC: restricted PIMC results [206]. Interpolation is applied to align the different data grid when necessary.
only available for high temperatures, and they are within 1% to our GDSMFB results shown in this figure. The deviations are comparable to the data accuracy, due to the statistical errors. OFMD simulations [207] 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.

V. WDM OUT OF EQUILIBRIUM
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 [63,208], semiclassical kinetics [71], quantum kinetic theory and nonequilibrium Green functions [44,76,209,210], hydrodynamics [73,81,211] or rate equations [212]. Among the problems that were studied are the equilibration of the electron distribution by electron-electron collisions [76,209], non-thermal melting induced by fs xray pulses [208], the density response for nonequilibrium momentum distributions [72,210], density evolution following short-pulse laser excitation [71], collisional heating of quantum plasmas by a laser pulse [213,214], or ionization dynamics in a short laser pulse [73].
In the reminder of this section we discuss the quantum hydrodynamics approach and its relation to DFT [82] more in detail because the former is comparatively little 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.

A. Dynamics of N quantum particles. Wave function and density operator
We consider a non-relativistic quantum system of N electrons described by the spin-independent hamiltonian where R = (r 1 , σ 1 ; r 2 , σ 2 ; . . . r N , σ N ), r i are the particle coordinates and σ i their spin projections, and V is an external potential, e.g. due to the plasma ions. Assuming first as 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 d 3N R |Ψ(R, t)| 2 = N . For particles with spin s, there are g s = 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 [40], where the sum runs over projection operators on all solutions of the Schrödinger equation (23), and p a are real probabilities, 0 ≤ p a ≤ 1, with a p a = 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 p a =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. [36,40,215]. 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 Wigner function) [40]: The equation of motion forF 1 follows straightforwardly from Eq. (25), e.g. Refs. [40,82] an is given below, cf. Eq. (29).

B. Time-dependent Kohn-Sham equations for electrons in WDM
We now derive the equation of motion for the timedependent 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 V xc .
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 form where the mean occupation numbers in equilibrium are given by the Fermi function (18). The equation forF 1 in the mean field (Hartree) approximation has the form [40], 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 approximations the N -particle wave function obeying Eq. (23) is just their product. Inserting the 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 ) U H [n(r, t)] = g s dr 2 w(r − r 2 )n(r 2 , t; β, µ), where the Hartree mean field is the coordinate representation of the 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 U H ). This result can be directly extended beyond the mean field approximation by replacing and, as a consequence, Eqs. (32), (33) become the timedependent Kohn-Sham equations-the basic equations of time-dependent density functional theory (TD-DFT) [216]. A particular strength of this theory is its solid theoretical foundation on the Runge-Gross theorem [216] and the corresponding theorems for time-independent DFT [217]. The basic statement is that a system of N interacting fermions can be mapped exactly on a system of N non-interacting particles with the same density n(r, t) where all interactions are lumped into an effective singleparticle 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 exchangecorrelation 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 V xc 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 finite temperature, V xc carries a temperature dependence which was discussed in detail before, see Sec. III.

C. Microscopic quantum hydrodynamic equations for dense plasmas
Following Ref. [82], we now derive the microscopic QHD (MQHD) equations, starting from the timedependent Hartree equations (32) and (33). To this end we simply convert each orbital solution, φ i (r, t) = A i (r, t)e i Si(r,t) , into an individual pair of amplitude and phase equations [79,82], for i = 1 . . . N , where n i = A 2 i and p i = ∇S i , and we introduced a short notation for the total potential energy, U tot = V + U H . This system of MQHD equations is fully equivalent to TD-DFT and, for V xc → 0, it exactly coincides with the time-dependent nonlinear Hartree (quantum Vlasov) equations [40,218]. Moreover, it was shown in Ref. [82] that this approximation, in linear response, is exactly equivalent to the random phase approximation (RPA or 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).

D. Derivation of the QHD equations from MQHD
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. In Ref. [78] the authors assumed that all orbital amplitudes are equal whereas, in Ref. [79], it was assumed that one can substitute Finally, in Ref. [82] 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,37,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, n i → n (this is what was used in Refs. [78,79]), 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,37,38), and obtain the QHD equations that contain three correlation functions that we denote by J np , J pp and Q ∆ , The correlation function J pp is directly related to the pressure tensor [82], which is seen by rewriting in coordinates as (α, β = x, y, z, summation over repeated Greek indices is implied) The first (diagonal) term is nothing but the average kinetic energy per particle, arising from the momentum variance which is non-zero, even in a pure state. Note that this term contains classical (thermal) and quantum contributions simultaneously. Also, no approximate "closure" as in classical hydrodynamics is required because this contribution is fully determined by the mean kinetic energy of the quantum plasma and by the corresponding approximations. For example, for an ideal Fermi gas in the thermodynamic limit (N → ∞, L → ∞, n = N/L 3 =const) at T = 0, the diagonal term is directly related to the pressure of a three-dimensional Fermi gas via P id F 0 (n) = 2 5 nE F (n). The second term contains non-diagonal elements of the pressure tensor, i.e. the mean shear stress, which is important for strongly correlated electrons. Thus, for a 3D ideal Fermi gas, we obtain The results for 1D and 2D have been given in Ref. [82]. Obviously, the pressure term and the correlation functions J np and Q ∆ are only of relevance in a spatially inhomogeneous system.

E. Plasma oscillations in MQHD and QHD
An important test for the QHD and MQHD models is the result for electron plasma oscillations (Langmuir waves) in the limit of 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 f eq i reduce to unity, for i ≤ E F , and zero otherwise.
Considering first the MQHD equations, Eqs. (36,37), and (38), we apply a harmonic monochromatic excitation, φ 1 (r, t) ∼ e −iωt+iqr , and linearize n i (t) and p i (t) around the unperturbed solution. Finally, the density response is computed via orbital averaging and Fourier transformation, with the result given by Eq. (51). Second, we consider the QHD equations (45), (46) which already contain the orbital averaging. Linearization of these equations and Fourier transformation yields the plasmon dispersion for the case of D = 1, 2, 3 dimensions [219], given by Eq. (52).
where we introduced the Fermi velocity via E F = mv 2 F /2. Let us discuss this result. First, we observe that the MQHD-dispersion (51) 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. (51) are correct. However, comparing with the QMC results (red curves in Fig. 7), the MQHD dispersion at intermediate wave numbers, 1 q/q F 3, strongly overestimates the oscillation frequency. This is due to the neglect of interaction 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 constructionless 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 2 ) 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 2 term in the dispersion (52) the pressure has to be multiplied by a factorᾱ that was reported in Ref. [81]: for example for a 3D plasma, the Fermi pressure that appears in Eq. (50) has to be multiplied by a factorᾱ = 9/5. Note that this value applies only to small wave numbers whereas for large wavenumbers, q 2k F , it will approach the value 3/5 [81].
F. Screened ion potential from MQHD and QHD.

Influence of electronic correlations
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. [127,[220][221][222][223]. For example, in Ref. [220] 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 5meV) at about 6 Bohr radii (about 2.5 interparticle distances). Tests with DFT simulations that can be regarded as benchmarks [221,224,225] revealed that no such minimum exists. The reason for the unphysical predictions of the QHD model was clarified in Ref. [127]: 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 [127]. 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. Figure 12. Prefactor of the Bohm potential for low-frequency long-wavelength excitations, as a function of the degeneracy parameter. γ changes from 1/9, at T = 0, to 1/3, at high temperature. Dashed line: high-frequency limit of γ. Figure  from Ref. [127], published with the permission of the authors.
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 (53). 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. [163]. 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 [128]. 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 subplot in Fig. 13) [128].
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 [3] and, in the case of strongly coupled ions, deviations from the RPA screened potential may be quite large and have significant impact on the structural and dynamical properties of the ion component in a dense plasma [3,226]. 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 artefact, 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 of future work.

G. Ion-acoustic modes in a quantum plasma
Let us now turn to ion-acoustic oscillations in a quantum plasma that have been studied by many authors. For example, Haas and Mahmood [211] 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(γ, c s ) = γ 2 /(12m e m i c 2 s ) represents the quantum electron correction due to Bohm term, and c s = k B T /m i Li 3/2 (z)/Li 1/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 c s = 2/3E F /m i . The dispersion of ion-acoustic wave is sensitive to the properties of the surrounding electrons via screening of the ion potential, Eq. (53), and to ionic correlations, depending on the coupling parameter Γ i . This is illustrated in Fig. 14 where the QHD result (54) is compared to the data from molecular dynamics (MD) simulation of ions at Γ i = 15, r s = 1.5 and θ = 0.1. The MD data was obtained using the screened ion potential (53) with (q, ω = 0) in RPA [226]. 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 QHD result for the dispersion, Eq. (54), 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 wavenumbers is inserted into Eq. (54). The resulting dispersion (black dashed curve) agrees well with the MD data for all wavenumbers, ka 1.0. The failure at larger wavenumber is not surprising due to the standard limitation of fluid approaches that we discuss in Sec. V H. To summarize, our analysis reveals that the functional form (54) of the dispersion is reasonable which is promising for further applications of the QHD approach to twocomponent plasmas, taking into account electronic and ionic non-ideality effects.

H. Limitations and further improvement of QHD: nonlocal and exchange-correlation effects
The general validity limits of quantum hydrodynamic models have been discussed in a variety of papers, e.g. Refs. [79,221,227]. Most importantly, similar as in a classical plasma, a hydrodynamic description is valid only on length scales larger than the screening length. In a quantum plasma this is the Thomas-Fermi length, λ T F , which leads to the criterion In the weak coupling limit, Γ q ≤ 0.1, this leads to the restriction k 0.03k F . Of course, this limitation can be avoided if the coefficients in the QHD equations are adjusted using information from the MQHD equations (kinetic theory). For example, it was shown in Ref. [81] how to properly choose the QHD parameters for large wavenumber, k 2k F . Moreover, tests against MQHD (RPA) results, discussed in the preceding sections, allow us to correct the prefactors in the QHD equations for the high-frequency and low-frequency limits.
Thus, even if correlation effects are neglected (V xc → 0, J np → 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 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 finite temperature [81]. 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, lim q→0 Π RPA (q, ω). When this is converted into local hydrodynamic equations via orbital averaging, the result is different for high and low frequency, respectively. It is possible to avoid this problem by introducing a more general nonlocal expression for the QHD potential that was derived in Ref. [81]. Here we summarize these results starting from the RPA and, in addition, including exchange-correlation effects that we link to the dynamic local field correction for which we have obtained ab initio results via QMC simulations, cf. Sec. III B. The main idea behind non-local quantum hydrodynamics is to require [81] that the QHD polarisation function equals to the polarisation function that follows from kinetic theory or TD-DFT (MQHD), i.e Π QHD ≡ Π LR . The derivation of the QHD equations presented in Sec. V D gives a strict justification for this requirement.
Let us return to the QHD momentum equation (46), considering zero vorticity, and introduce the total potential µ[n(r, t)] that contains ideal and exchangecorrelation contributions (first and second terms), µ[n(r, t)] = µ id [n(r, t)] + µ xc [n(r, t)] + eφ(r, t), (59) whereas the last term is due to the total field -the sum of external field as well as mean field (Hartree) contributions which is the solution of Poisson's equation, i.e. eφ = U tot . The force field µ is defined by the functional derivative of the grand potential [81]: where Ω = Ω id +Ω xc , is the sum of an ideal and exchangecorrelation part that will be specified below.
In equilibrium (current-free case), Eqs. (58) and (60) reduce to the Euler-Lagrange equation δn 0 (r, t) where the subscript 0 indicates the equilibrium case. Assuming a weak perturbation [81,228], the force becomes where n 1 (r, t) = n(r, t) − n 0 (r), with |n 1 /n 0 | 1. Eqs. (58)- (62) and the requirement that the QHD density response agrees with the kinetic theory results in linear response, Π QHD ≡ Π LR , give [81]: where we introduced the modified linear response polarization, Π LR from which the long-wavelength limit at finite ω, Π 0 (ω) = mk 2 /ω 2 n 0 , is being subtracted, and F denotes the Fourier transform to frequency and wavenumber (k, ω) space. Considering first the case of ideal quantum electrons, where Π LR = Π RPA , we have: Eq. (65) yields a non-local potential, Eq. (62), and allows, among others, to find the correction factors α(ω, q; Θ, D) and γ(ω, q; Θ, D), Eqs. (57) and (56), in various limiting cases [81]. Consider now the general case of nonideal electrons, r s 1, Γ 1 1, cf. Sec. II. This requires to include the exchange-correlation contribution to the grand potential. This can be done by using approximations for the exchange-correlation potential, V xc , from DFT. A simple ground state expression in the local density approximation was introduced to QHD in Ref. [229]) and has been used in a number of subsequent papers, e.g. in Ref. [230]. A second strategy that is closer to the topics discussed in this paper is to directly use ab initio input from QMC. This is indeed possible via the dynamic local field correction G(q, ω), that was computed in Sec. III C, [81]: Equations (58)-(60), (65) and (66) represent a closed set of equations which is exact in the weak perturbation case, i.e. |n 1 /n 0 | 1, and can be summarized in one generalized nonlocal momentum balance equation [81]: where ϕ 1 (r, t) = ϕ(r, t) − ϕ 0 (r) .
This is a remarkable result that contains all relevant limiting cases and assures the highest accuracy possible via a link to quantum kinetic theory and ab initio input from QMC.
Let us discuss the limitations of the result (67). The main assumption in the derivations above is the validity of linear response, i.e. |n 1 /n 0 | 1. If this is not satisfied, equations (58)-(60), (65) and (66) can still be used, but the accuracy will be largely defined by the form of the functional Ω[n], e.g. see [231,232], and the results will be not reliable. This concerns, in particular, applications to nonlinear oscillations and waves in quantum plasmas. In this case, the linear response result, Π LR , has, in principle, to be replaced by solutions of a nonlinear kinetic equation. Another cases that is beyond the scope of Eq. (67) is very rapid external excitation. In this case the distribution function may be far from a Fermi function f eq which give rise to a strongly modified plasmon spectrum in quantum plasmas, e.g. [233,234], similar to the case of classical plasmas. The relevance of nonequilibrium plasmas under warm dense matter conditions was studied, e.g. by Gericke et al. [72]. For these situation, a nonequilbrium quantum kinetic theory is required that yields the time-dependent density response, eq (q, ω; t) [210] and local field corrections, G eq (q, ω; t) which assumes an equilibrium form in which the distribution function is replaced by a nonequilibrium function, f eq → f (t). However, even this approach maybe inappropriate if the excitation is on the scale of the plasma period or faster, t 2π/ω pl . In that case, the formation of the screening cloud and of the plasmon spectrum proceeds on the time scale of the excitation, and a true nonequilibrium theory for the density response functions is required, cf. Ref. [40] and references therein.

VI. CONCLUSIONS AND OUTLOOK
In this article we have presented an overview on recent simulation results for warm dense matter. First, we have presented thermodynamic results for the degenerate electron component, considering the warm dense uniform electron gas, that are based on ab initio quantum Monte Carlo simulations. The results include highly accurate parametrizations of the exchange-correlation free energy and results for the static local field correction G(q). Furthermore, we presented ab initio results for dynamic quantities, including the dynamic local field correction G(q, ω), the dynamic structure factor S(q, ω) and the dynamic density response χ(q, ω). These results can be further extended to other dynamic quantities that are of high relevance for current warm dense matter experiments. Moreover, the static and dynamic local field correction are of key importance as input for other models and simulations methods, including quantum hydrodynamics and density functional theory. DFT today is the only approach that is, eventually, capable to cover the entire WDM range including electron-ion plasmas as well as the condensed matter phase. These simulations give access to time-dependent properties: in Born-Oppenheimer MD the electrons adiabatically following the ions. In contrast, in time-dependent DFT a real (non-adiabatic) time-dependence of electrons and ions is described.
However, DFT is notoriously inaccurate in treating ionization energies, band gaps and electronic correlation effects, in particular under WDM conditions. Here, ab initio input from QMC can substantially improve the simulations. One such input is the exchange-correlation free energy that has been used to improve the local density approximation, cf. Sec. III A. While the finitetemperature effect on the LDA equation of state turned out to be just on the order of a few percent, this seems to bring the results significantly closer to QMC simulations of dense plasmas, cf. Sec. IV. Even more promising are the QMC data for the static and dynamic local field cor- rection that provide highly valuable input for improved exchange-correlation kernels of TD-DFT and QHD, cf. Secs. III C and V H.
At the same time, the present versions of DFT and TD-DFT are not able to describe thermalization and electronic correlation effects such as Auger processes and double excitation. These processes require a kinetic description in the frame of quantum kinetic theory or nonequilibrium Green functions (NEGF). However, both NEGF and TD-DFT are extremely computationally costly and, therefore, are currently limited to short length and time scales. Larger length and time scales are accessible with simpler approaches, such as DFT-MD [cf. Sec. IV] or the quantum Boltzmann equation [40]. A qualitative summary of the length and time scales that can be described by the different methods is given in Fig. 15. There we also include molecular dynamics with quantum potentials (semiclassical MD) which is applicable when the electronic quantum dynamics are not important. We also indicated the possibility to extend these simulations to larger length scales by means of parallelization, whereas longer time scales can only be achieved, in some cases, by acceleration concepts, e.g. by combination with analytical models, for a discussion see Refs. [235][236][237].
Thus, there is still a big gap between the length and time scales that are accessible by microscopic simulations and the scales that are relevant in experiments. Here, a reliable fluid theory for fermions that is similarly advanced and successful as in classical plasmas such as the QHD studied above, could serve as the missing link. As we have pointed out, the orbital averaging involved in deriving the QHD equations restricts this model to finite length and time scales. This means, fast processes, in particular related to the thermalization of the electron distribution, as well as small length scales on the order of the Bohr radius or the Thomas-Fermi screening length, cannot be resolved. For these effects more advanced methods, such as DFT and quantum kinetic theory have to be used (see below). At the same time, the relative simplicity of the QHD equations allows one, to extend them to large length scales and propagate them to long times where the aforementioned approaches cannot be appliied. Thus, QHD could be a very useful tool that is complementary to other methods. This situation is indicated qualitatively in Fig. 15.