Electronic Density Response of Warm Dense Hydrogen: Ab initio Path Integral Monte Carlo Simulations

The properties of hydrogen under extreme conditions are important for many applications, including inertial confinement fusion and astrophysical models. A key quantity is given by the electronic density response to an external perturbation, which is probed in X-ray Thomson scattering (XRTS) experiments -- the state of the art diagnostics from which system parameters like the free electron density $n_e$, the electronic temperature $T_e$, and the charge state $Z$ can be inferred. In this work, we present highly accurate path integral Monte Carlo (PIMC) results for the electronic density response of hydrogen. We obtain the exchange-correlation (XC) kernel $K_{xc}$, which is of central relevance for many applications, such as time-dependent density functional theory (TD-DFT). This gives us a first unbiased look into the electronic density response of hydrogen in the warm-dense matter regime, thereby opening up a gamut of avenues for future research.

The properties of hydrogen under extreme conditions are important for many applications, including inertial confinement fusion and astrophysical models. A key quantity is given by the electronic density response to an external perturbation, which is probed in X-ray Thomson scattering (XRTS) experiments-the state of the art diagnostics from which system parameters like the free electron density ne, the electronic temperature Te, and the charge state Z can be inferred. In this work, we present highly accurate path integral Monte Carlo (PIMC) results for the electronic density response of hydrogen. We obtain the exchange-correlation (XC) kernel Kxc, which is of central relevance for many applications, such as time-dependent density functional theory (TD-DFT). This gives us a first unbiased look into the electronic density response of hydrogen in the warm-dense matter regime, thereby opening up a gamut of avenues for future research.
Matter at extreme densities (n ∼ 10 21−27 /cm −3 ) and temperatures (T ∼ 10 4−8 K) is ubiquitous throughout our universe [1]. This so-called warm dense matter [2] (WDM) naturally occurs in a number of astrophysical objects such as giant planet interiors [3,4], brown dwarfs [5], and the outer layers of neutron stars [6,7]. In addition, WDM is highly important for technological applications such as the discovery of novel materials [8,9] and hot-electron chemistry [10]. Moreover, a hydrogen fuel capsule has to traverse the WDM regime [11] on its compression path towards inertial confinement fusion (ICF) [12].Consequently, WDM is a highly active research area and a strongly increasing number of experiments [13] are performed at large research facilities such as the National Ignition Facility (NIF) at the Lawrence Livermore National Laboratory [14] or the European X-FEL in Hamburg, Germany [15]. These developments have led to a number of spectacular recent discoveries,including the high-precision measurement of the stopping power in WDM [16], the use of record peak-brightness free-electron lasers to study shock compressed aluminium [17], and the possible observation of the molecular-to-metallic transition in hydrogen at extreme pressure [18,19].
At the same time, we stress that the interpretation of such experiments decisively depends on a rigorous theoretical modelling. In fact, even basic system parameters like the density or electronic temperature are generally unknown and have to be inferred. In this regard, the Xray Thomson scattering (XRTS) technique [20,21] has emerged as the most promising method of diagnostics of WDM experiments.
Unfortunately, the theoretical description of WDM constitutes a most formidable challenge due to the highly nontrivial interplay of a number of physical effects [2,22,23]. WDM states are correlated due to the Coulomb in- * t.dornheim@hzdr.de teraction, the electrons (and sometimes the nuclei) are partially quantum degenerate, and WDM is a highly excited state, which rules out the well-stocked arsenal of ground-state quantum many-body methods [24,25]. This intricacy is usually expressed in terms of two characteristic parameters, that are both of the order of unity [26]: a) the density parameter r s = r/a B , where r and a B are the average electronic distance and first Bohr radius, and b) the degeneracy temperature θ = k B T /E F , with E F being the Fermi energy [27].
In this situation, thermal density functional theory (DFT) [29] has emerged as the work horse in WDM theory [2], as it often combines a manageable computation cost with a reasonable degree of accuracy. Yet, the results can substantially depend on the employed exchange-correlation (XC) functional [30], which cannot be obtained within DFT itself, and has to be supplied as a semi-empirical a-priori input; see Ref. [31] for a critical discussion focusing on hydrogen.
An additional challenge pertaining to the interpretation of XRTS experiments is the wave-number dependence of the measured signal. The ab initio estimation [32] of the corresponding dynamic structure factor S(q, ω) (beyond model assumptions such as the Chihara decomposition [33,34]) requires one to carry out timedependent DFT (TD-DFT) simulations [35], which need as input the frequency-and wave vector-dependent XCkernel K xc (q, ω). Yet, little is known about the actual K xc (q, ω) of a real WDM system, and previous studies [33,36] have been exclusively based on drastic simplifications such as the adiabatic local density approximation (ALDA). This is highly problematic, since these approximate kernels have been designed for the application at T = 0 and, in the case of ALDA, are valid only for small wave numbers q = |q| and break down for strong degrees of inhomogeneity [37].
In the present work, we aim to overcome these limitations by performing ab initio path integral Monte Carlo (PIMC) calculations [38,39] for the electronic density re- Perturbation strength dependence of the electronic density response of hydrogen at q = 2πL −1 (0, 0, 2) T (q ≈ 1.5qF) at the electronic Fermi temperature (θ = 1) for rs = 2 (left), rs = 4 (center), and rs = 6 (right). The green crosses show our new PIMC data for hydrogen, and the red circles corresponding results for a uniform electron gas partly taken from Ref. [28]. The dashed black line shows a cubic fit to the PIMC data, and the solid blue horizontal lines the linear-response limit.
sponse of warm dense hydrogen. This allows us to thoroughly address a number of fundamental questions, including i) the quantification of the impact of electron-ion correlations on electronic properties, ii) the assessment of the validity range of linear response theory and the importance of nonlinear effects [28,40,41], and iii) to check widespread model assumptions about the decomposition into bound and free electrons [21]. Most notably, we present the first, highly accurate results for the XCkernel of a realistic WDM system in the static limit (i.e., ω → 0). Our theoretical predictions are directly useful for upcoming experiments with hydrogen, for example at NIF or the European X-FEL. Results: We have carried out thermal DFT-MD simulations to obtain a set of ionic configurations, and use the PIMC method to obtain exact solutions to the electronic problem in the external ionic potential. This allows us to directly compare our PIMC results to DFT, and use the former as input for the latter. We note that we impose no nodal restrictions on the thermal density matrix [42]. Therefore, our PIMC simulations are without approximation, but computationally expensive due to the notorious fermion sign problem [43,44]. We estimate the full computation cost of the present study to be of the order of O 10 7 CPUh. To compute the electronic density response, we apply an external harmonic perturbation [28,[45][46][47] of wave vector q and perturbation amplitude A, which leads to the full electronic Hamiltonian (we assume Hartree atomic units throughout this work) cos (q ·r l ) . (1) Here N = N ↑ + N ↓ is the total number of electrons, and we restrict ourselves to the unpolarized case of N ↑ = N ↓ = N/2. The operatorsŴ ee andV ei denote the interaction between the electrons and the external potential due to the fixed ions, respectively. We use the PIMC method to compute the expectation value of the electronic density in Fourier space ρ(q, A) = ρ q q,A , and the sought-after density response is simply given by The results for the density response of warm dense hydrogen are shown in Fig. 1 for r s = 2 (left), r s = 4 (center), and r s = 6 (right) at the electronic Fermi temperature, i.e., θ = 1. Specifically, we show ∆ρ/A, and the green crosses (red circles) depict our new PIMC results for hydrogen (for a uniform electron gas (UEG) [28]). The main qualitative trends are as follows. For the metallic density of r s = 2, the electrons are only weakly localized around the ions, and the density response both qualitatively and quantitatively resembles the behaviour of a UEG at the same conditions. The central panel corresponds to a lower electronic density as it is realised experimentally for example within hydrogen jets [48]. In this case, a sizeable fraction of the electrons is assumed to be involved in bound states with the protons. Consequently, we observe a starkly reduced density response compared to the UEG, as only the unbound (free) electrons react to the external perturbation; see the discussion of Fig. 2 below. Furthermore, we observe that the data sets for hydrogen and the UEG converge in the limit of large A, since the external perturbation will eventually predominate over the ionic potential. Lastly, the right panel shows results for r s = 6. While being hard to probe in current experimental set-ups, these conditions are interesting as a challenging benchmark for theoretical methods due to the pronounced impact of electronic XC-effects [49,50]. In this case, we find that the electronic density response is reduced by an order of magnitude compared to the UEG at small A, due to the high probability of bound states. An additional advantage of the present approach is its capability to study nonlinear effects and, in this way, to unambiguously assess the validity range of linear response theory. Overall, we find similar trends for hydrogen as in previous investigations of the UEG [28,40], and nonlinear effects are even increased in hydrogen compared to the UEG at r s = 6.
The central task of the present work is the estimation of the linear electronic density response function. The density response at the first harmonic of the original perturbation (i.e., at the same q) can be expanded in powers of the perturbation amplitude [40] as ∆ρ(q, A) = χ(q)A + χ cubic (q)A 3 , and the dashed black curves depict corresponding fits [52] to the PIMC data for sufficiently small A. The horizontal blue lines depict the linear coefficient χ(q) := χ(q, 0), i.e., the static limit of the linear density response function [53] Here χ 0 (q, ω) can be either the Kohn-Sham response function [54], or the Lindhard function [27] in the case of a uniform system. In addition, G(q, ω) = − q 2 4π K xc (q, ω) denotes the local field correction (LFC), which contains the full wave-vector and frequency-resolved information about electronic XC-effects, and which is equivalent to the XC-kernel K xc usually employed in the context of TD-DFT [35]. Evidently, both G(q, ω) and K xc (q, ω) depend on the particular choice of χ 0 (q, ω) and, in the context of DFT, on the employed XC-functional.
The relevant linear electronic density response is shown in Fig. 2. The green symbols show our new PIMC results for different N that have been obtained by following the fitting procedure of ∆ρ(q, A) outlined above for different wave vectors q. Evidently, no systematic dependence on the system size can be resolved, which is consistent to previous findings for the UEG [51,55,56]. The density response of the UEG has been included as dotted red curves and is in good agreement to the green symbols at r s = 2; we only observe small deviations around q ∼ 1.5q F . In stark contrast, there appear pronounced differences between the UEG and hydrogen (exceeding 30%) at r s = 4, which is a direct consequence of the increased location of the electrons around the ions. It is often assumed that the total number of electrons can be decomposed into a bound and into an effectively free fraction; we denote the latter as α = N free /N in this work. This leads to the modified parameters r s (α) and θ(α). Empirically, we find that the choice of α = 0.6 leads to a good qualitative agreement between the UEG model and our PIMC results for hydrogen at r s = 4 for small to intermediate wave numbers, see the dash-dotted yellow curve in the bottom panel of Fig. 2. This is close to the result of α ≈ 0.54 by Militzer and Ceperley [57] based on a cluster analysis in real space. Interestingly, the agreement between the PIMC data for hydrogen and the effective free electron model deteriorates towards large q, where the response of hydrogen even slightly exceeds the response of the full UEG model. This is a direct consequence of the relevant length scales l = 2π/q. Excitations with small q correspond to large distances in coordinate space that only affect free electrons. Large q, on the other hand, are directly translated into small l, on which bound electrons, too, can react to the external perturbation.
The blue dots in Fig. 2 depict the Kohn-Sham response function χ 0 (q) that has been computed within the local density approximation (LDA) [58]. For r s = 2, these data closely resemble the Lindhard function [27] (dashed black) describing the uniform ideal Fermi gas. This is expected, as the Kohn-Sham orbitals resemble plane waves when localization is weak. The black stars show the DFT response function within the random phase approximation, i.e., by setting G(q, 0) = 0 in Eq. (2). Evidently, this leads to the correct asymptotes for large and small q, but becomes inaccurate around q ∼ 1.5q F where electronic XC-effects are known to be important [23]. For r s = 4, the Kohn-Sham response function exhibits a more interesting behaviour: for large q, it approaches the Lindhard function of the full electronic density, whereas it more closely resembles the reduced Lindhard function corresponding to α = 0.6 (solid grey) for q 2q F . In other words, some information about electron-ion correlations, in general, and about bound states, in particular, is included in the KS-response. Remarkably, the corresponding RPA data are in good agreement to the exact PIMC results; more so than for the less strongly coupled case of r s = 2. We also include the KS response function that has been obtained from a DFT simulation without any XC-effects for r s = 4 as the grey circles. Interestingly, this leads to substantial differences compared to the blue dots for small to medium q, and to a deterioration in the corresponding RPA response function (purple triangles).
Our new PIMC results for the electronic density response of hydrogen give us access to the XC-kernel of a realistic WDM system on a true ab initio level. The results are shown in Fig. 3 for the same conditions as in Fig. 2. The dotted red curves show results for the full UEG model [51] and the dashed black lines the ALDA, which is a parabolic expansion around q → 0. Evidently, the latter is only appropriate for q < 2q F even in the case of a UEG. Let us next consider the different results that have been obtained from our PIMC results for hydrogen for χ(q) by inverting Eq. (2) using as input different χ 0 (q). The blue (green) symbols have been obtained using the KS response within LDA (the Lindhard function). At r s = 2, the DFT kernel agrees well qualitatively with the LFC of the UEG. Interestingly, the kernel that has been computed in terms of the Lindhard function exhibits substantial deviations to the other data sets over the entire q-range as it has to balance the absence of electron-ion correlations from χ 0 (q).
For r s = 4, the situation is more complicated as none of the data sets computed from χ(q) resemble the UEG model. For q 2q F , the LDA-RPA response function is in good agreement with the PIMC results, and the thus extracted LFC is small in magnitude in this regime. This changes for intermediate q, where the LFC becomes comparable in magnitude to the UEG model, but with a negative sign. This is a direct consequence of the overestimation of the actual response by the RPA in the case of hydrogen, whereas it is well known [51,55] that the opposite holds for the UEG. Using the response function of a uniform ideal Fermi gas to compute the LFC enhances this trend, and the KS-response from a DFT simulation with no XC-effects is located in between. In addition, this means that using either the ALDA or full UEG model as the XC-kernel in a TD-DFT calculation of hydrogen leads to an actually worse estimation of the electronic density response compared to the bare RPA at r s = 4.
Conclusion: In this work, we have presented the first, highly accurate PIMC results for the electronic density response of hydrogen in the WDM regime. This has allowed us to unambiguously quantify the importance of nonlinear effects, which are even more important compared to the previously investigated case of a uniform electron gas [28,40,60]. Moreover, we have obtained extensive data for the exact linear-response limit of the static density response function χ(q, 0), which have been used to extract the first exact results for the static LFC G(q) and XC-kernel K xc (q). We stress that the latter in general depends on the employed χ 0 (q, ω), which can be either a Lindhard function or a Kohn-Sham response function that depends on the particular choice of the XC-functional. This has profound consequences in the presence of bound states, where the actual XC-kernel does not even qualitatively resemble the familiar UEGmodels. Consequently, the commonly employed ALDA is appropriate when most electrons are free and behave similarly to a UEG, but dramatically breaks down when the degree of electronic localization around the protons is substantial; see the Supplemental Material [54] for a corresponding TD-DFT study.
Our results indicate a strong need for the development and re-evaluation of presently used model XC-kernels such as ALDA. A particular advantage of our study is given by the direct possibility to compare our PIMC results for different properties to corresponding DFT results for the same ionic configuration. This, in turn, will give us invaluable lessons regarding the performance of different XC-functionals [37,61], and guide the development of new approaches [62][63][64][65][66]. Finally, we re-iterate the capability of our set-up to study nonlinear effects in warm dense hydrogen, which may give unprecedented insights into many-body correlation effects in WDM [41] and are known to sensitively depend on important system parameters such as the electronic temperature [67].

PATH INTEGRAL MONTE CARLO
We use the direct PIMC method [1][2][3] to simulate the electronic system governed by the Hamiltonian Eq. (1) of the main text. A typical snapshot is presented in Fig. 1, where we show a configuration of N = 20 electrons at θ = 1 and r s = 2. The green orbs represent the (fixed) ions, and the smaller blue-red paths the configuration of the electrons. We do not impose any nodal restrictions [4] on the structure of the thermal density matrix. Therefore, our simulations are exact within the given statistical error bars but computationally expensive.
Being based on a decomposition of the thermal density matrix, the PIMC method becomes exact only in the limit of an infinite number of imaginary-time steps P . In this work, we use the pair approximation that is based on the exact solution of the two-body problem between an electron and an ion following the procedure described by Militzer [6]. In Fig. 2, we show the convergence with P of the density in Fourier space ρ(q, A) [see the main text] at r s = 4 and q = 3.35q F . The top and bottom panels show results for the unperturbed case (A = 0), and a perturbation of A = 0.07. The green crosses correspond to our implementation of the pair approximation, and the propagator error is small even for as few as P = 50. To ensure that any residual convergence error is negligible, we use P = 500 throughout this work. As an additional benchmark, we have also implemented a decomposition of the density matrix based on the diagonal Kelbg potential [5]. The results are shown as the red circles in Fig. 2 and exhibit a substantially larger propagator error compared to the pair approximation, as it is expected. At the same time, they nicely converge towards the green crosses in the limit of P → ∞, which is a strong indication for the correctness of our implementation.

IMPACT OF THE PIMC XC-KERNEL ON TD-DFT SIMULATIONS
We have carried out linear-response time-dependent DFT (TD-DFT) calculations to assess the impact of the electronic XC-kernel on the dynamic structure factor (DSF) S(q, ω) of the electrons in hydrogen; all relevant DFT parameters are listed below. The results are shown in Fig. 3, and the top (bottom) panel shows the DSF at θ = 1 for q ≈ 1.69q F at r s = 2 (r s = 4). The dotted blue and green curves depict corresponding results for the UEG [8,9] that have been obtained within RPA, and by using the static LFC of the UEG, respectively. Evidently, the incorporation of the latter leads to a correlation induced red-shift in both cases, which is more pronounced at the lower density, as it is expected. The ω-dependence of G(q, ω) has a negligible impact on the DSF of the UEG at these conditions, and the dotted red curves are exact for a pure electron gas without the ions [8,9]. The solid black and red curves show TD-DFT results (averaged over N snap = 13 independent snapshots) based on LDA calculations on the RPA level and using the appropriate exact static XC-kernel based on PIMC results for hydrogen shown in Fig. 3 in the main text, respectively. To explicitly focus on the impact of K xc (q), we compute the electron-electron part of the full dynamic structure factor of hydrogen, which, however, intrinsically depends on the ionic structure. The full DSF includes an additional electron-ion term (leading to the well-known inelastic peak) which we do not consider in this work. For r s = 2, both curves are qualitatively similar to the corresponding UEG curves, as it is by now expected. In stark contrast, we observe substantial qualitative differences to the UEG at r s = 4. In particular, our TD-DFT results for S(q, ω) exhibit a diffusive peak induced by the localization of the electrons in the attractive potential of the ions around ω = 0 (see also Refs. [10,11] for observations of such a feature in other systems) and a second peak around twice the electronic plasma frequency, ω = 2ω p . Interestingly, the incorporation of the appropriate electronic LFC has a negligible impact onto the TD-DFT results, as the RPA is nearly exact in the static limit, see the main text. We re-iterate our earlier observation that this is the opposite trend compared to the UEG, where the impact of the LFC increases monotonically with r s . Finally, the yellow curves show TD-DFT results for S(q, ω) that have been obtained on the basis of the ALDA model. Evidently, this approximation is well justified for r s = 2, where hydrogen behaves similarly to an UEG. Contrary, ALDA leads to a significant deterioration of the quality of S(q, ω) at r s = 4. In order words, it can indeed be preferable to neglect the XC-kernel altogether, instead of employing models that have been derived in a different context, and under assumptions that can be strongly violated in realistic WDM applications.
The KS-DFT calculations have been performed using the GPAW code [14][15][16][17], which is a real-space implementation of the projector augmented-wave method. A Monkhorst-Pack [18] sampling of the Brillouin zone was used. The calculations are carried out using a plane-wave basis.
For the calculation of the static density response function χ 0 (q, 0), the following parameters have been used: For r s = 2, θ = 1 (T = 12.5280 eV), and the number of particles N = 14 (N = 20), the main simulation cell size is L = 4.1Å (L = 4.631Å) with a spacing of h = 0.1708Å (h = 0.193Å), the cut-off energy in the thermal equilibrium calculation is set to be E cut = 550 eV (E cut = 650 eV), the number of bands is given by 1500 (1100), and a k -point grid of 4 × 4 × 4. The wavefunctions from the thermal equilibrium calculation are used to compute the static density response function χ 0 (q, 0) using a plane-wave cut-off of 325 eV (150 eV), where the wave numbers q = jq min , with j = 1...5 (j = 1...4) and q min = 2π/L, have been considered. The data point for q min /12 is generated using 300 bands, with E cut = 550 eV, h = 0.185Å, and a 12 × 12 × 12 k -point grid.
For the calculation of the dynamic density response function χ 0 (q, ω) that has been used to compute S(q, ω) shown in Fig. 3, the following parameters are used: At r s = 2 (r s = 4), χ 0 (q, ω) is computed using the main simulation cell size L = 4.1Å (8.224Å ) with a spacing of h = 0.25Å, the cut-off energy in the calculation of the equilibrium state is set to be E cut = 350 eV, the number of bands 500 (400), and the k -point grid of 12×12×12 (4×4×4). The wavefunctions from the equilibrium state calculation are used to compute χ 0 (q, ω) with a plane-wave cut-off of 80 eV and the broadening parameter γ = 0.2 eV.