The properties of hydrogen at extreme temperatures, densities, and pressures are of paramount importance for a wealth of applications.1,2 An excellent example is given by inertial confinement fusion (ICF),3 where the fuel capsule (most commonly a deuterium–tritium fuel) has to traverse such warm dense matter (WDM) conditions on its pathway toward ignition.4 In several recent breakthrough experiments, it has been demonstrated at the National Ignition Facility (NIF)5 that it is indeed possible to ignite the fuel6 and even produce a net energy gain with respect to the laser energy deposited in the capsule.7,8 While undoubtedly a number of challenges remain on the way to an operational fusion power plant,9 there is hope that ICF might emerge as a future source of safe and sustainable energy. A second class of applications is given by astrophysical objects, such as giant planet interiors1,10–12 (e.g., Jupiter in our solar system, but also exoplanets) and brown dwarfs.13,14 Here, the properties of hydrogen are of key relevance to describing the evolution of these objects, and to understanding observations such as the gravitational moments of Jupiter.15–17

Despite its apparent simplicity, a rigorous theoretical description of hydrogen remains notoriously elusive over substantial parts of the relevant phase diagram.18–20 A case in point is given by the insulator-to-metal phase transition, which is highly controversial from both a theoretical18,21–23 and an experimental24–26 perspective. The situation is hardly better within the WDM regime that is of particular relevance for both ICF and astrophysical applications.1,27

Specifically, the WDM regime is usually defined in terms of the Wigner–Seitz radius rs=(3/4πne)1/3 (with n_e = N_e/Ω the electron density and Ω the volume) and the degeneracy temperature Θ = k_BT/E_F (with E_F the Fermi energy of an electron gas at equal density28 and T the temperature), both of the order of unity,29 i.e., r_s ∼ Θ ∼ 1. This means that in the Hamiltonian of the system, both kinetic and interaction terms are of the same order.28 Thus, there is an absence of small parameters, necessitating a full treatment of the complex interplay between strong thermal excitations (generally ruling out electronic ground-state methods), Coulomb coupling between both electrons and ions (often ruling out weak-coupling expansions such as Green functions30), and quantum effects such as Pauli blocking and diffraction [precluding semiclassical schemes such as molecular dynamics (MD) simulations with effective quantum potentials31]. Indeed, there exists no single method that is reliable over the entire WDM regime.27 In practice, a widely used method is given by a combination of molecular dynamics (MD) simulations of the ions with a thermal density functional theory (DFT)32 treatment of the degenerate electrons based on the Born–Oppenheimer approximation. On the one hand, such DFT-MD simulations are often computationally manageable and, in principle, provide access to various material properties, including the equation of state,33–36 linear response functions,37–39 and the electrical conductivity.40–43 On the other hand, the accuracy of DFT depends strongly on the employed exchange–correlation (XC) functional,18,44 which has to be supplied as an external input. While the accuracy of many functionals is reasonably well understood in the ground state,45 the development of novel, thermal XC functionals that are suitable for application at extreme temperatures has only started very recently.43,46–49 Moreover, many calculations require additional input, such as the XC kernel for linear-response time-dependent DFT (LR-TDDFT) calculations of WDM and beyond.50

Indeed, the linear density response of a system to an external perturbation19 constitutes a highly important class of material properties in the context of WDM research. Such linear-response theory (LRT) properties are probed, for example, in X-ray Thomson scattering (XRTS) experiments,51,52 which have emerged as a key diagnostic of WDM.53,54 In principle, the measured XRTS intensity gives one access to the equation-of-state properties of the probed material,55–57 which is of paramount importance for ICF applications and the modeling of astrophysical objects. Unfortunately, in practice, the interpretation of the XRTS intensity is usually based on uncontrolled approximations such as the widely used Chihara decomposition58 into effectively bound and free electrons. Thus, the quality of inferred parameters such as temperature, density, and ionization remains generally unclear.

LRT properties are also ubiquitous throughout WDM theory and central to the calculation of stopping power and electronic friction properties,59,60 electrical and thermal conductivities,61,62 opacity,63,64 ionization potential depression,65 and effective ion–ion potentials,66,67 as well as the construction of advanced nonlocal XC functionals for thermal DFT simulations.68,69

Owing to its key role, the electronic linear density response has been extensively studied for the simplified uniform electron gas (UEG) model,28,70 where ions are treated as a rigid homogeneous charge-neutralizing background; see Ref. 19 and references therein. These efforts have culminated in different parametrizations of the static local field factor,71–73 which is formally equivalent to the aforementioned XC kernel—a basic input for LR-TDDFT simulations and other applications. Very recently, Böhme et al.74,75 have extended these efforts and obtained quasi-exact ab initio path integral Monte Carlo (PIMC)76–78 results for the density response of an electron gas in the external potential of a fixed proton snapshot. On the one hand, these results are directly comparable with thermal DFT simulations of hydrogen, which has given important insights into the accuracy of various XC functionals.50,79,80 On the other hand, these calculations are computationally very inefficient—requiring a large set of independent simulations to estimate the density response at a single wavenumber for a given snapshot—and, more importantly, miss the important dynamic interplay between electrons and protons. The latter is crucial to capture effects such as electronic localization around the ions in observables such as the electronic density response function; indeed, simulation results are not comparable with experimental observations without it.

In the present work, we overcome these fundamental limitations by presenting the first ab initio PIMC results for the linear density response of full two-component warm dense hydrogen. By treating the electrons and protons dynamically on the same level (i.e., without invoking the Born–Oppenheimer approximation), we access all components of the density response, including the electron–electron, electron–proton, and proton–proton local field factors (i.e., XC kernels). Indeed, a consistent treatment of the interaction between the electrons and ions is crucial to capture the effects of electronic localization around the protons, which has important implications for the interpretation of XRTS experiments. These effects are particularly important for solid state densities (r_s = 3.23), where they are significant over the entire relevant wavenumber range. They are also substantial for metallic densities (r_s = 2) and even manifest at strong compression (r_s = 1) for small wave numbers.

Our simulations are quasi-exact; no fixed-node approximation81 is imposed. To deal with the fermion sign problem82,83—an exponential computational bottleneck in quantum Monte Carlo simulations of degenerate Fermi systems—we average over a large number of Monte Carlo samples, making our simulations computationally very expensive; the cost of this study is estimated to be O(107) CPUh. Additionally, we employ the recently introduced ξ-extrapolation method84–88 to access larger system sizes. We find that finite-size effects are generally negligible under these conditions, which is consistent with previous results for the UEG model.70,89,90

We are convinced that our results will open up multiple opportunities for future research. First, the quasi-exact nature of our results makes them a rigorous benchmark for computationally less costly but approximate approaches such as thermal DFT or PIMC within the fixed-node approximation. Moreover, the obtained LRT properties constitute direct predictions for upcoming XRTS experiments with hydrogen jets91,92 and ICF plasmas.5 This makes them also ideally suited to benchmark Chihara models54,58,93 that are commonly used to diagnose XRTS measurements and infer an equation of state. Of particular value are our results for the various partial local field factors, which can be used as input for innumerable applications such as transport property estimates or the construction of novel XC functionals for thermal DFT simulations. Finally, we note that the study of the density response of warm dense hydrogen presented here is interesting in its own right, and gives us new insights into the interplay of electronic localization, quantum effects, and electron–ion correlations in the WDM regime.

The remainder of the paper is organized as follows. In Sec. II, we provide the relevant theoretical background: the full hydrogen Hamiltonian (Sec. II A), a brief introduction to the PIMC method and its estimate of imaginary-time correlation functions (ITCFs) (Sec. II B), the linear density response theory (Sec. II C), and its relation to XRTS experiments (Sec. II D). In Sec. III, we present our extensive novel simulation results covering the cases of metallic density (Sec. III A), solid state density (Sec. III B), and strong compression (Sec. III C). The paper is concluded by a summary and outlook in Sec. IV. Additional technical details are provided in the Appendixes A and B.

We assume Hartree atomic units (i.e., 4πϵ₀ = ℏ = m_e = e = 1) throughout this work, unless otherwise specified. The full Hamiltonian governing the behavior of hydrogen is given byĤH=−12∑l=1N∇l,e2−12mp∑l=1N∇l,p2+∑k=1N∑l=1Nl<kWE(r̂l,r̂k)+∑k=1N∑l=1Nl<kWE(Îl,Îk)−∑k=1N∑l=1NWE(Îl,r̂k),(1)where the first and second terms correspond to the kinetic energy of the electrons and protons, respectively, and m_p ≈ 1836.15 is the proton mass. The pair interaction W_E is given by the usual Ewald summation, and we follow the conventions introduced in Ref. 94; r̂ and Î denote the position operators of electrons and protons.

The ab initio PIMC method76,77,95 constitutes one of the most successful tools for the description of interacting, quantum degenerate many-body systems at finite temperature. The central property is given by the canonical (i.e., the inverse temperature β = 1/T, volume Ω = L³, simulation box length L, and numbers of both electrons and protons are fixed) partition function evaluated in the coordinate representationZN,Ω,β=1N↑!N↓!∑σN↑∈SN↑∑σN↓∈SN↓ξNpp×∫dR〈R|e−βĤH|π̂σN↑π̂σN↓R〉,(2)where the meta-variable R=(r1,…,rN,I1,…,IN)T contains the coordinates of both electrons (r_l) and protons (I_l). Throughout this work, we consider the fully unpolarized case with an equal number of spin-up and spin-down electrons, N_↑ = N_↓ = N/2. Equation (2) contains a summation over all possible permutations σNi of electrons of the same spin orientation i ∈ {↑, ↓}, which are realized by the corresponding permutation operator π̂σNi; N_pp denotes the total number of pair permutations of both spin-up and spin-down electrons. Note that we treat the heavier protons as distinguishable quantum particles (sometimes referred to as boltzmannons in the literature), which is exact under the investigated conditions. In fact, even protonic quantum delocalization effects are only of the order of 0.1% here; cf. Figs. 3, 9, and 13 below. The variable ξ determines the type of quantum statistics of the electrons, with ξ = −1, ξ = 0, and ξ = 1 corresponding to Fermi–Dirac, Boltzmann–Maxwell, and Bose–Einstein statistics, respectively. Only ξ = −1 has distinct physical meaning for warm dense hydrogen, although other values can give valuable insights into the importance of quantum degeneracy effects for different observables.87,88

A detailed derivation of the PIMC method is beyond our scope and has been presented elsewhere.75,77,78 In essence, the complicated quantum many-body problem of interest, as defined by Eq. (2), is mapped onto an effectively classical system of interacting ring polymers with P segments each. This is the celebrated classical isomorphism.96 A schematic illustration of this idea is presented in Fig. 1, where we show a configuration of N = 3 hydrogen atoms in the τ–x plane, with τ ∈ [0, β] and t = −iτ the imaginary time, and x is a spatial coordinate. In the path integral picture, each particle is represented by an entire path of particle coordinates along the imaginary time τ with P imaginary-time slices. The filled and empty circles correspond to the coordinates of electrons and ions, respectively, and beads of the same particle on adjacent time slices are connected by a harmonic spring potential. The latter connections are depicted by the red and green lines for electrons and protons. The extension of these paths corresponds to the thermal wavelength λTa=2πβ/ma, which is smaller by a factor of 1/mp≈0.023 for the protons. The basic idea of the PIMC method is to use the Metropolis algorithm97 to generate a Markov chain of such path configurations X, where the meta-variable X contains all coordinates and also the complete information about the sampled permutation structure. A first obstacle is given by the diverging Coulomb attraction, which prevents the straightforward application of the Trotter formula.98 To avoid the associated path collapse, we employ the well-known pair approximation,75,77,99 which effectively leads to an off-diagonal quantum potential that remains finite even at zero distance. An additional difficulty arises from the factor (−1)Npp in the case of Fermi statistics. For every change in the permutation structure, the sign of the configuration weight changes. For example, we show a configuration with N_pp = 1 pair exchange in Fig. 1 (i.e., the two electrons on the right), which has a negative configuration weight. In practice, the resulting cancellation of positive and negative terms leads to an exponential increase in the computational time toward low temperatures and large numbers of electrons N; this bottleneck is known as the fermion sign problem in the literature;82,83 see also Appendix A for additional details. While different, often approximate, strategies to mitigate the sign problem have been suggested,81,100–104 here, we carry out exact direct PIMC solutions that are subject to the full amount of cancellations. The comparably large noise-to-signal level is reduced by averaging over a large number of independent Monte Carlo samples, and our results are exact within the given error bars.

An additional advantage of the direct PIMC method concerns its straightforward access to various imaginary-time correlation functions (ITCFs).78,105–107 In the context of the present work, the focus lies on the density–density ITCF (hereinafter referred to simply as ITCF)Fab(q,τ)=⟨n̂a(q,0)n̂b(−q,τ)⟩,(3)which measures the decay of correlations between particles of species a and b along the imaginary-time diffusion process; see the yellow lines and symbols in Fig. 1. As we shall see in Sec. II C, the ITCF gives direct access to the static density response of warm dense hydrogen from a single simulation of the unperturbed system. The ITCF is interesting in its own right108,109 and can be measured in XRTS experiments;53,110,111 see Sec. II D.

The basic idea of the density response formalism is to study a system that is subject to a weak external perturbation that couples linearly to the Hamiltonian:19Ĥq,ω;Ae,Ap=ĤH+2Ae∑l=1Ncos(q⋅r̂l−ωt)+2Ap∑l=1Ncos(q⋅Îl−ωt).(4)Here, ĤH denotes the unperturbed Hamiltonian of the full hydrogen system [Eq. (1)], and q and ω are the wavevector and frequency of the external monochromatic perturbation. It should be noted that we distinguish the perturbation amplitudes of the electrons (A_e) and protons (A_p), which allows us to study different species-resolved components of the density response, as we shall discuss in detail. In the limit of an infinitesimally small external perturbation strength, the induced density is a linear function of A_a (i.e., the amplitude of the perturbation acting on species a), and linear response theory is applicable:19,112δ⟨ρ̂a(q,ω)⟩=χaa(0)(q,ω)×Aa+∑bvab(q)[1−Gab(q,ω)]δ⟨ρ̂b(q,ω)⟩,(5)where the induced density perturbation δ⟨ρ̂a⟩ is evaluated in reciprocal space and where χaa(0)(q,ω) denotes the linear density response of a noninteracting system (of species a) under the same conditions; v_ab denotes the interaction between species a and b. Further, the complete information about species-resolved exchange–correlation effects is contained in the local field factors G_ab(q, ω). Equation (5) defines a coupled system of equations that explicitly depends on the perturbation A_a of all components.

The fundamental property in linear density response theory is the linear density response function χ_ab(q, ω) that describes the density response of the species a to a potential energy perturbation applied to the species b:δ⟨ρ̂a(q,ω)⟩=Abχab(q,ω).(6)This dynamic susceptibility constitutes a material property of the unperturbed system and is directly related to the dielectric function.113,114 Moreover, it is of central importance to LR-TDDFT simulations,115 and for the interpretation of XRTS experiments with WDM and beyond, as we explain in Sec. II D below.

Let us focus on the static limit χ_ab(q, 0) ≡ χ_ab(q) that describes the response to a time-independent cosinusoidal perturbation. Moroni et al.116,117 have presented the first accurate results for the linear static density response function and the local field factor of the ground state UEG based on diffusion quantum Monte Carlo simulations. This idea was subsequently adapted to PIMC and configuration PIMC (CPIMC) simulations of the UEG in the WDM regime.118,119 More specifically, the basic idea is to utilize the formal perturbation strength expansion of nonlinear density response theory, which readsδ⟨ρ̂a(q)⟩=χab(1,1)(q)Ab+χab(1,3)(q)Ab3+O(Ab5),(7)with χab(m,l)(q) the partial static density response of order l at a harmonic m, where χab(1,1)(q)≡χab(q).120 This is a particular case of a general result stating that the total induced density at odd (even) harmonics is given by an infinite series of all l ≥ m odd (even) powers of the perturbation strength with coefficients equal to χab(m,l)(q).120 Thus, the induced density is estimated for various A_b, and a polynomial expansion19,116,121δ⟨ρ̂a(q)⟩=c1Ab+c3Ab3+O(Ab5)(8)allows us to identify the linear and cubic density response functions at the first harmonic (i.e., at the wavenumber of the perturbation) via the correspondence c₁ = χ_ab(q) and c3=χab(1,3)(q). On the one hand, this direct perturbation method is formally exact, and can easily be adapted to other methods such as DFT.50,122,123 On the other hand, it is computationally very expensive, since independent simulations need to be carried out for multiple perturbation amplitudes A_b just to extract the response χ_ab(q) for a single wavevector at a given density and temperature combination. An elegant alternative is given by the imaginary-time version of the fluctuation–dissipation theorem,108 whose species-resolved version readsχab(q,0)=−NaNbΩ∫0βdτFab(q,τ),(9)which implies that the species-resolved density responses at all accessible wavevectors can be extracted from a single simulation of the unperturbed system by estimating the corresponding partial F_ab(q, τ). This method has been extensively applied to study the density response of the UEG at finite temperatures over a vast range of parameters.71,124–130 Here, we employ this second route for the bulk of our new results, while we use the direct perturbation route as an independent cross-check.

Let us next consider the local field factors G_ab(q, ω) in more detail. For the UEG, the polarization potential approach leads to the well-known linear density response expression28,70,131,132χee(q,ω)=χee(0)(q,ω)1−vee(q)[1−Gee(q,ω)]χee(0)(q,ω).(10)Setting G_ee(q, ω) ≡ 0 in Eq. (10) leads to the random phase approximation (RPA), which describes the electronic density response on a mean-field level. It should be pointed out that spin-resolved UEG generalizations have been presented in the literature,112,133 where the spin-up and spin-down electrons are treated as two distinct species.

For multicomponent systems, the coupling between the different species makes the situation considerably more complicated. We note that the terms partial and species-resolved are used interchangeably throughout this paper for the associated LRT quantities of multicomponent systems. Introducing the so-called vertex-corrected interaction28θab(q,ω)=vab(q)[1−Gab(q,ω)],(11)the polarization potential approach leads to the following expression for the linear density perturbation δn induced by the potential energy perturbation δU:134δna(q,ω)=χaa(0)(q,ω)δUa(q,ω)+∑cχaa(0)(q,ω)θac(q,ω)δnc(q,ω),(12)which clearly constitutes the generalized version of Eq. (5) for non-monochromatic perturbations. The functional derivative definition of the partial density response function χ_ab(q, ω) = δn_a(q, ω)/δU_b(q, ω) and the identity δU_a(q, ω)/δU_b(q, ω) = δ_ab yieldχab(q,ω)=χaa(0)(q,ω)δab+∑cχaa(0)(q,ω)θac(q,ω)χcb(q,ω).(13)The ideal density response functions χab(0) are nonzero only if a = b. Note the reciprocal connections, namely, χ_ab(q, ω) = χ_ba(q, ω), θ_ab(q, ω) = θ_ba(q, ω), and G_ab(q, ω) = G_ba(q, ω), which are a consequence of Newton’s third law v_ab(q) = v_ba(q). The above 3 × 3 set of linear equations can be explicitly solved for the partial density response functions in terms of their dependence of θ_ab(q, ω). For hydrogen, a, b, c = {e, p}, this leads toχee(q,ω)=1Δ(q,ω)χee(0)(q,ω)[1−θpp(q,ω)χpp(0)(q,ω)],χpp(q,ω)=1Δ(q,ω)χpp(0)(q,ω)[1−θee(q,ω)χee(0)(q,ω)],χep(q,ω)=1Δ(q,ω)[χee(0)(q,ω)θep(q,ω)χpp(0)(q,ω)],where the auxiliary response Δ(q, ω) is the determinant given byΔ(q,ω)=[1−θee(q,ω)χee(0)(q,ω)][1−θpp(q,ω)χpp(0)(q,ω)]−θep2(q,ω)χee(0)(q,ω)χpp(0)(q,ω).When first-principles results are available for χ_ab(q, ω), as in our case for the static limit of ω = 0, the above 3 × 3 set of linear equations can be explicitly solved for the hydrogen partial local field factors:θee(q,ω)=1χee(0)(q,ω)−χpp(q,ω)χee(q,ω)χpp(q,ω)−χep2(q,ω),(14)θpp(q,ω)=1χpp(0)(q,ω)−χee(q,ω)χee(q,ω)χpp(q,ω)−χep2(q,ω),(15)θep(q,ω)=χep(q,ω)χee(q,ω)χpp(q,ω)−χep2(q,ω).(16)The limit of weak electron–proton coupling is expressed as χee(q,ω)χpp(q,ω)≫χep2(q,ω), which gives the one-component expressionsθee(q,ω)=1χee(0)(q,ω)−1χee(q,ω),θpp(q,ω)=1χpp(0)(q,ω)−1χpp(q,ω).The opposite limit of strong electron–proton coupling is expressed as χee(q,ω)χpp(q,ω)≪χep2(q,ω), which again gives a one-component expression, but without a noninteracting contribution:θep(q,ω)=−1χep(q,ω).

The measured intensity in an XRTS experiment can be expressed as19,51,110I(q,ω)=See(q,ω)⊛R(ω),(17)where S_ee(q, ω) is the dynamic structure factor (DSF) and R(ω) the combined source and instrument function. In practice, R(ω) is often known with sufficient accuracy from source monitoring, or from the characterization of a utilized backlighter X-ray source.135 We note that a direct deconvolution to solve Eq. (17) is generally not stable, owing to noise in the experimental data; a model for S_ee(q, ω) is usually convolved with R(ω) and compared with I(q, ω) instead. The connection between S_ee(q, ω) and the dynamic density response function introduced in Sec. II C is given by the well-known fluctuation–dissipation theorem:28See(q,ω)=−Im{χee(q,ω)}πne(1−e−βω).(18)In combination with the Kramers–Kronig causality relation that connects Im{χ_ee(q, ω)} and Re{χ_ee(q, ω)},28 the fluctuation–dissipation theorem directly implies that χ_ee(q, ω) and S_ee(q, ω) contain exactly the same information. Moreover, the static density response function χ_ee(q) is related to the DSF via the inverse frequency moment sum rule:136χee(q)=−2ne∫−∞∞dωSee(q,ω)ω.(19)Therefore, accurate knowledge of χ_ee(q) gives direct insights into the low-frequency behavior of S_ee(q, ω), which is dominated by electron–proton coupling effects such as localization and ionization. However, direct evaluation of the right-hand side of Eq. (19) on the basis of experimental XRTS data is generally prevented by the convolution with R(ω); see Eq. (17). In practice, this problem can be circumvented elegantly by considering the relation between S_ee(q, ω) and the ITCF F_ee(q, τ):Fee(q,τ)=L[See(q,ω)]=∫−∞∞dωSee(q,ω)e−τω,(20)where L[⋯] denotes the two-sided Laplace transform operator. Making use of the convolution theoremL[See(q,ω)]=L[See(q,ω)⊛R(ω)]L[R(ω)],(21)one can thus deconvolve the measured scattering intensity in the imaginary-time domain. In fact, the evaluation of Eq. (21) turns out to be remarkably stable with respect to experimental noise53,110,111 and thus gives one direct access to the ITCF from an XRTS measurement. It is then straightforward to obtain the static electron–electron density response function χ_ee(q) from the experimentally measured ITCF via Eq. (9), which amounts to the imaginary-time analog of Eq. (19). Our new results for the static density response of full hydrogen thus constitute an unambiguous prediction for upcoming XRTS experiments with hydrogen jets, fusion plasmas, etc.

To gain additional insights into the physical implications of χ_ee(q), we consider the widely used Chihara decomposition of the DSF:54,58,93See(q,ω)=Sel(q,ω)+Sbf(q,ω)+Sff(q,ω)︸Sinel(q,ω).(22)The basic idea of Eq. (22) is to divide the electrons into effectively bound and free populations. The first contribution to the full DSF is given by the pseudo-elastic componentSel(q,ω)=WR(q)δ(ω),(23)where the Rayleigh weight contains both the atomic form factor of bound electrons and a screening cloud of free electrons.137 Thus, this term originates exclusively from the electronic localization around the protons. The second contribution to the full DSF contains all inelastic contributions, i.e., transitions between bound and free electronic states (and the reverse process, free–bound transitions93) described by S_bf(q, ω), as well as transitions between free states described by S_ff(q, ω). It is important to note that the decomposition into bound and free states is arbitrary in practice and breaks down at high compression, where even the orbitals of bound electrons overlap;138 the PIMC method that we use in the present work does not distinguish between bound and free electrons and consequently is not afflicted by these problems. Nevertheless, the simple picture of the Chihara model (22) gives important qualitative insights into the density response of the system. For example, Eq. (19) directly implies that χ_ee(q) is highly sensitive to electronic localization around the ions, which should result in an increased density response of hydrogen compared with the UEG model;71,72 this is indeed what we infer from our PIMC results, as we shall see in Sec. III.

We use the extended ensemble sampling scheme from Ref. 139, implemented in the Imaginary-Time Stochastic High-Performance Tool for Ab initio Research (ISHTAR) code,140 which is a canonical adaption of the worm algorithm by Boninsegni et al.78,141 All results are freely available online142 and can be used as input for other calculations, or as a rigorous benchmark for other methods. A short discussion of the ξ-extrapolation method,84,86 used to simulate larger systems, is given in Appendix A.

In this section, we investigate in detail the linear density response of hydrogen at a metallic density r_s = 2 and at the electronic Fermi temperature Θ = 1. We begin our analysis with a study of the ITCF F_ab(q, τ), which constitutes the basis for a significant part of the present work. In Figs. 2(a)–2(c), we illustrate F_ee(q, τ), F_pp(q, τ), and F_ep(q, τ), respectively, in the relevant q–τ plane. Note that the symmetry relation F_ab(q, τ) = F_ab(q, β − τ) holds (see also Fig. 3 below), as a consequence of the imaginary-time translation invariance at thermodynamic equilibrium or, equivalently, from the detailed balance relation S_ab(q, ω) = e⁻S_ab(q, −ω) for the DSF.53,108 Thus, we can restrict ourselves to discussion of the interval τ ∈ [0, β/2]. The partial electron–electron ITCF F_ee(q, τ) shown in Fig. 2(a) exhibits a rich structure that is the combined result of multiple physical effects. The τ = 0 limit corresponds to the static structure factor F_ab(q, 0) = S_ab(q), with S_ee(q) approaching unity for large wavenumbers, and a finite value for q → 0.88 In addition, F_ee(q, τ) exhibits an increasingly pronounced decay with τ for large q. From a physical perspective, this is due to the Gaussian imaginary-time diffusion process that governs the particle path in the path-integral picture.108,109 In essence, the electrons are quantum delocalized, with the extension of their imaginary-time paths being proportional to the thermal wavelength, as has been explained in the discussion of Fig. 1 above. With increasing q, one effectively measures correlations on increasingly small length scales λ = 2π/q. For small λ, any correlations completely decay along the imaginary-time diffusion, and F_ee(q, β/2) goes to zero. From a mathematical perspective, this increasing τ decay also follows from the f-sum rule, which states that19,108,111,143∂∂τFab(q,τ)τ=0=−δabq22ma.(24)These different trends can also be seen in Figs. 3(a) and 3(b), where we show the τ dependence of F_ee(q, τ) for q = 1.53 Å⁻¹ and q = 7.65 Å⁻¹, respectively, as the solid red curves. Additional insights come from a comparison with the corresponding UEG results under the same conditions (the double-dashed yellow curves). For the smaller q value, the DSF of the UEG is dominated by a single broadened plasmon peak in the vicinity of the plasma frequency,144 leading to a moderate decay with τ. For full hydrogen, we find a very similar τ dependence, but the entire curve appears to be shifted by a constant offset; this is a direct signal of the elastic feature [Eq. (23)] in S_ee(q, ω) that originates from the electronic localization around the protons. For q = 7.65 Å⁻¹, on the other hand, the UEG and full hydrogen give very similar results for F_ee(q, τ). In fact, they are indistinguishable around τ = 0, since S_ee(q) = 1 holds in the large-q limit, and both curves satisfy Eq. (24) exactly. Interestingly, we observe a slightly reduced τ decay in H compared with the UEG model around τ = β/2, which has small though significant implications for χ_ee(q), as we shall see below.

Let us next consider the proton–proton ITCF F_pp(q, τ) and the electron–proton ITCF F_ep(q, τ) which are shown in Figs. 2(b) and 2(c), respectively. First and foremost, both correlation functions appear to be independent of τ over the entire depicted q range. For F_ep(q, τ), this is indeed the case within the given MC error bars; see also the blue curves in Fig. 3. This is consistent with the f-sum rule [Eq. (24)], which states that F_ep(q, τ) is constant at least at first order in τ. The proton–proton ITCF exhibits a richer behavior on a magnified view; see the green curves in Fig. 3. Within our PIMC simulations, the protons are treated as delocalized quantum particles, just like the electrons, although their thermal wavelength is smaller by a factor of 1/m_p ≈ 0.02. This is reflected by the less extended paths along the imaginary-time diffusion process (cf. Fig. 1) and consequently by a strongly reduced τ decay of F_pp(q, τ) compared with F_ee(q, τ). This is also reflected by the mass in the denominator of the right-hand side of Eq. (24). Remarkably, for large q, we can resolve a small yet significant τ dependence of F_pp(q, τ), of the order of ∼0.1%; see the right inset in Fig. 3(b). Overall, it is still clear that the behavior of both F_pp(q, τ) and F_ep(q, τ) is essentially governed by the corresponding static structure factors S_pp(q) and S_ep(q), while the τ dependence is essentially negligible for the linear density response in this regime.

We now proceed to the central topic of the present work, namely, the investigation of the partial static density response of warm dense hydrogen. Figure 4 depicts the static density response function χ_ab(q) as a function of the wavenumber q. The results represented by the red, green, and blue symbols have been computed from the ITCF through Eq. (9) for the electron–electron, proton–proton, and electron–proton responses, respectively. The crosses and diamonds correspond to N = 14 and N = 32 atoms, and no dependence on the system size can be resolved. Before comparing the datasets with the other depicted models and calculations, we shall first summarize their main trends:(i)The electron–proton density response χ_ep(q) = χ_pe(q) has the same negative sign as χ_ee(q) and χ_pp(q), since the unperturbed protons will follow the induced density of the perturbed electrons, and vice versa.(ii)The electron–proton density response χ_ep(q) decays monotonically with q as the electron–proton coupling vanishes in the single-particle limit; the same trend also leads to the well-known decay of the elastic feature in S_ee(q, ω) quantified by the correspondingly vanishing Rayleigh weight W_R(q).(iii)The electron–electron density response χ_ee(q) is relatively flat for q ≲ 2 Å⁻¹ and decays monotonically for larger q. Equation (9) directly implies that this is a quantum delocalization effect. In practice, the static density response is proportional to the area under the corresponding ITCF. The latter vanishes increasingly rapidly with τ as q increases, as discussed above, leading to the observed reduction in χ_ee(q). Heuristically, this can be understood as follows. While the static density response of a single (or noninteracting) classical particle is independent of wavenumber, the response of a delocalized quantum particle is reduced when its thermal wavelength is comparable to the perturbation wavelength. In fact, a quantum particle will stop reacting altogether when its extension is much larger than the excited wavelength.(iv)The proton–proton density response χ_pp(q) increases with q and seemingly becomes constant for q ≳ 6 Å⁻¹. The reduction in χ_pp(q) for small q is a consequence of the proton–proton coupling (and its interplay with electrons), whereas its behavior for large q comes from the larger proton mass and the resulting strongly reduced quantum delocalization. For completeness, note that we can resolve small deviations from the classical limit of χppcl(q)=−npβ for large wavenumbers. We also note that the general property |χAAq(q)|≤|χAAcl(q)|, valid for static linear response functions associated with any Hermitian operator Â,²⁸ should be respected by all the species-resolved density responses.

Of equal interest as these observations about the density response of full two-component hydrogen are their comparisons with other models and calculations. The gray squares in Fig. 4 have been obtained from the direct perturbation approach, i.e., from independent PIMC simulations of hydrogen where one component has been perturbed by an external harmonic perturbation [cf. Eq. (4)]. This procedure is further illustrated in Fig. 5, where we show the induced density ρ_a(q) as a function of the perturbation amplitude. More specifically, the red crosses in Fig. 5(a) correspond to the electronic density of full two-component hydrogen induced by the electronic perturbation amplitude A_e (with A_p = 0) for a wavenumber q = 1.53 Å⁻¹. In the limit of A_e → 0, ρ_e(q)/A_e attains a finite value that is given by the static linear density response function χ_ee(q); the latter is shown as the horizontal dashed blue line, as computed from the ITCF F_ee(q, τ) via Eq. (9), and it is in excellent agreement with the red diamonds in the limit of small A_e. The solid green lines show cubic fits via Eq. (7). Evidently, the latter nicely reproduce the A_e dependence of the PIMC data for moderate perturbations, and the linear coefficient then corresponds to the linear density response function; the deviations between data and fits for A ≳ 10 eV can be taken into account by including higher-order terms in Eq. (7),120 which will be pursued in future work. The red crosses have been obtained from the same set of simulations (i.e., A_e > 0 and A_p = 0) and depict the corresponding induced density of the protons that is described by χ_ep(q). The unperturbed protons thus do indeed follow the perturbed electrons, although with a somewhat reduced magnitude. This can be discerned particularly well in Fig. 6, where we show the density in real space along the direction of the perturbation for a comparably small electronic perturbation amplitude of A_e = 1.36 eV. The solid red line shows the PIMC results for the electronic density and the dotted blue curve the linear-response theory estimate129ne(r)=n0+2Ae⁡cos(q⋅r)χee(q),(25)using the ITCF-based result for χ_ee(q). Both curves are in excellent agreement, which implies that LRT is accurate in this regime. The green curve shows the proton density from the same calculation. It exhibits the same cosinusoidal form, but with a reduced magnitude. Let us postpone a discussion of the other curves in Fig. 6 and return to Fig. 5(a). The yellow squares show results for the induced proton density ρ_p that have been obtained from a second, independent, set of PIMC calculations with a finite proton perturbation amplitude A_p > 0, but unperturbed electrons, A_e = 0. We find excellent agreement with the ITCF-based result for the LRT limit of χ_pp(q), and the cubic fit agrees well with the data. The protons react more strongly to an external static perturbation compared with the electrons; see the discussion of Fig. 4 above. Finally, the yellow crosses show the induced electron density ρ_e from the same calculation. We recover the expected symmetry χ_ep(q) = χ_pe(q) within the linear response limit. Interestingly, this breaks down for larger perturbation amplitudes, where the protons again react more strongly than the electrons. This nonlinear effect deserves further explanation, and will be investigated in detail in a dedicated future work.

In Appendix B, we present more results for the direct perturbation approach for different wavenumbers, which have been employed to obtain the linear density response functions that are shown as the empty squares in Fig. 4. We find perfect agreement with the ITCF-based datasets everywhere.

Let us now consider the solid yellow curve in Fig. 4, which shows the density response of the UEG model71 under the same conditions. We find good (though not perfect) agreement with the electronic density response of full hydrogen for q ≳ 4 Å⁻¹. In stark contrast, there appear to be substantial differences between the two datasets for χ_ee(q) in the limit of small q, where the density response of the UEG vanishes owing to perfect screening.145 This is not the case for hydrogen, owing to the electron–proton coupling, which can be understood intuitively in two different ways. From the perspective of density response theory, the fact that we study the response in the static limit of ω → 0 implies that the protons have sufficient time to follow the perturbed electrons, thereby effectively screening the Coulomb interaction between the latter. Ionic mobility breaks down the perfect screening relation of the UEG, and allows the electrons to react even to perturbations on very large length scales (i.e., q → 0), leading directly to the nonvanishing value of χ_ee(q) for large wavelengths. Additional insight comes from the relation of χ_ee(q) as the inverse frequency moment of the dynamic structure factor S_ee(q, ω); cf. Eq. (19) in Sec. II D above. For the UEG, S_ee(q, ω) consists simply of a sharp plasmon excitation around the (finite) plasma frequency for small q, and the weight of the plasmon vanishes quadratically in this regime.28 For full two-component hydrogen, on the other hand, S_ee(q, ω) contains additional contributions (i) from bound–free transitions and (ii) from the quasi-elastic feature, which is usually modeled as a sum of an atomic form factor and a screening cloud of free electrons.137 The latter feature increases with small q and, being located at very small frequencies, is manifested strongly in the inverse moment of S_ee(q, ω); the static electron–electron density response function is thus highly sensitive to electronic localization around the protons. We note that this has potentially important implications for the interpretation of XRTS experiments with WDM, since χ_ee(q) can be directly inferred from the measured intensity (cf. Sec. II D), and, from a theoretical perspective, it does not require dynamic simulations. We therefore suggest that after having inferred the temperature from the model-free ITCF thermometry method introduced in Refs. 53 and 110, one might calculate χ_ee(q) for a given wavenumber over a relevant interval of densities to infer the latter from the XRTS measurement. Such a strategy would completely circumvent the unphysical decomposition into bound and free electrons, while at the same time being very sensitive to a related, but well-defined, concept: electronic localization around the ions. For light elements such as H or Be,87 this might even be accomplished on the basis of quasi-exact PIMC results, offering a pathway for, in principle, approximation-free WDM diagnostics. At the same time, we should point out that the static density response might also be estimated with reasonable accuracy from computationally less demanding methods such as DFT or restricted PIMC (using the direct perturbation approach), since no dynamic information is required. A dedicated exploration of this idea thus constitutes an important route for future research. Let us conclude our comparison between the electronic density response of full two-component hydrogen and the UEG by inspecting the corresponding density profile in Fig. 6, which is shown as the solid yellow curve. As expected, the UEG reacts less strongly to an equal external perturbation. Note that the nearly perfect agreement with the green curve is purely coincidental and is a consequence of the intersection of the yellow and blue data in Fig. 8 for the considered wavenumber.

Finally, the data represented by the solid black dots in Fig. 4 have been taken from Böhme et al.74 and show χ_ee(q) of a nonuniform electron gas in the external potential of a fixed proton configuration. Evidently, keeping the protons fixed has a dramatic impact on the electronic density response. Those electrons that are located around a proton are substantially less likely to react to the external harmonic perturbation than the electrons in a free-electron gas,146 leading to an overall reduction of the density response for small to intermediate q. In particular, this snapshot-based separation of the electrons and protons completely misses the correct signal of the electronic localization around the protons that has been discussed in the previous paragraph. While being ideally suited to benchmark DFT simulations, and potentially to provide input to the latter, the physical content of these results is thus quite incomplete.

Let us conclude our analysis of warm dense hydrogen at r_s = 2 by considering the various local field factors θ_ab(q) [cf. Eqs. (14)–(16)] that are shown in Fig. 7. The solid yellow line corresponds to the machine learning (ML) representation of the UEG71 and has been included as a reference. The solid red line shows the electron–electron local field factor θ_ee(q). It attains the same limit of −1 for q → 0, but deviates substantially from the UEG result for all finite wavenumbers. The straightforward application of UEG models for the description of real hydrogen is thus questionable in this regime. This discrepancy becomes even more pronounced for solid state density, but vanishes for r_s = 1; cf. Secs. III B and III C below. The green curve shows the proton–proton local field factor θ_pp(q). Interestingly, it is basically indistinguishable from θ_ee(q) for q ≲ 4 Å⁻¹, but diverges from the latter for large wavenumbers. Finally, the blue curve shows the electron–proton local field factor θ_ep(q). Unsurprisingly, it is larger in magnitude than the other datasets for small to moderate wavenumbers, which reflects the importance of electron–proton coupling effects.

As a second example, we investigate hydrogen at Θ = 1 and r_s = 3.23, i.e., the density of solid state hydrogen. From a physical perspective, the low density is expected to lead to an increased impact of both electron–electron and electron–proton coupling effects,74,147,148 making these conditions a very challenging test case for simulations. We note that in contrast to the UEG, the lower density is also more challenging for our PIMC setup. This is a consequence of the incipient formation of H⁻ ions and molecules,88 leading to a more substantial degree of quantum degeneracy and consequently a more severe sign problem. Indeed, we find an average sign of S ≈ 0.05 for N = 14 hydrogen atoms, causing a factor of 1/S² = 400 in the required computational time. In addition, low-density hydrogen is expected to exhibit interesting physical effects, such as the roton-type feature149 in the dynamic structure factor S_ee(q, ω) for intermediate wavenumbers.150 Intriguingly, such conditions can be probed with XRTS measurements of optically pumped hydrogen jets,91,92 which makes our results directly relevant for upcoming experiments.

In Fig. 8, we show our new PIMC results for the different ITCFs F_ab(q, τ) in the relevant τ–q plane. Overall, we find the same qualitative trends as with r_s = 2 investigated above, but with two main differences: (i) the electron–proton ITCF F_ep(q, τ) attains larger values on the depicted q grid, indicating a higher degree of coupling between the two species; (ii) both F_ee(q, 0) and F_pp(q, 0) exhibit reduced decay for small q, for the same reason.

Additional insight comes from Figs. 9(a) and 9(b), where we show the various ITCFs for q = 0.95 Å⁻¹ and q = 4.73 Å⁻¹, respectively, along the τ axis. In Fig. 9(a), the main difference from the r_s = 2 case is the larger offset between the results for F_ee(q, τ) from the UEG and full two-component hydrogen; this is a direct consequence of the increased electronic localization around the protons and the correspondingly increased Rayleigh weight W_R(q). For the larger q value in Fig. 9(b), we again find that no dependence of F_ep(q, τ) can be resolved within the given confidence interval (shaded blue area). By contrast, we can clearly resolve protonic quantum effects; see the right inset. An additional interesting observation comes from a comparison of the solid red and double-dashed yellow curves corresponding to F_ee(q, τ) for hydrogen and the UEG. In the limit of τ → 0, the two curves are in perfect agreement, which is a consequence of the fact that S_ee(q) ≈ 1 and the f-sum rule [Eq. (24)] yielding the same slope for both datasets. For larger τ, on the other hand, we observe substantial differences between hydrogen and the UEG. Specifically, we find a reduced τ decay for the former system compared with the latter, which cannot simply be explained by a constant shift due to the Rayleigh weight. From the perspective of our PIMC simulations, this clearly indicates that electron–electron correlations are stabilized along the imaginary-time diffusion process by the presence of the protons. Equivalently, we can attribute this finding to a shift of spectral weight in S_inel(q, ω) to lower frequencies (see also Ref. 108 for a more detailed discussion), indicating a nontrivial structure of the full DSF. The presented data for F_ee(q, τ) thus constitute rigorous benchmarks for models (e.g., the Chihara decomposition) and simulations (e.g., LR-TDDFT), and a dedicated future comparative analysis will provide important insights into the ranges of validity of different methods.

In Fig. 10, we show the corresponding species-resolved static density response functions χ_ab(q), and again restrict ourselves to a discussion of the main differences from the r_s = 2 case:(i)χ_ee(q) increases monotonically with decreasing q over the entire range of q considered, and the electronic localization around the protons predominantly shapes its behavior for small q.(ii)In contrast to the UEG, χ_ee(q) does not converge toward the ideal density response χee(0)(q) for large q, which is a direct consequence of the reduced τ decay of F_ee(q, τ) depicted in Fig. 8(b).(iii)χ_ep(q) is substantially larger, with respect to χ_ee(q) and χ_pp(q), for r_s = 3.23 compared with r_s = 2.(iv)χ_pp(q) exhibits a reduced decay for small q compared with r_s = 2; this is due to the electronic screening of the proton–proton interaction, making the proton response more ideal.

Finally, we show the partial local field factors θ_ab(q) in Fig. 11. While the electron–electron local field factor of full two-component hydrogen (red curve) attains the same limit as the UEG model71 (yellow curve) for q → 0, it exhibits the opposite trend for intermediate q, followed by a steep increase for shorter wavelengths. This is consistent with previous findings by Böhme et al.74 for an electron gas in a fixed external proton potential under similar conditions (r_s = 4) and clearly indicates the breakdown of UEG models when electrons are strongly localized. The proton–proton local field factor is relatively featureless over the entire range of q, which might be due to the aforementioned electronic screening of the proton–proton interaction. Similar to r_s = 2, θ_ep(q) constitutes the largest local field factor for relevant wavenumbers.

As the final example, we investigate compressed hydrogen at Θ = 1 and r_s = 1. The corresponding PIMC results for the species-resolved ITCFs are shown in Fig. 12 in the τ–q plane and qualitatively closely resemble the case of r_s = 2 shown in Fig. 3. The main difference is the reduced magnitude of F_ep(q, τ), indicating a substantially weaker localization of the electrons around the protons, as expected. In Fig. 13, we show the ITCFs along the τ direction for q = 3.06 Å⁻¹ and q = 15.3 Å⁻¹, i.e., for the same values of q/q_F as in Figs. 3 and 9. We find very similar behavior as for r_s = 2, with a reduced electron–proton coupling. In fact, no deviations can be resolved between F_ee(q, τ) for the UEG model and full two-component hydrogen in Fig. 13(b). In other words, compressed hydrogen closely resembles a free-electron gas.

This observation is further substantiated in Fig. 14, where we show the corresponding partial static density response functions χ_ab(q). Evidently, χ_ee(q) is very similar to that in the UEG and converges toward the latter for q ≳ 5 Å⁻¹. Nevertheless, we can still clearly detect the signature of electronic localization around the protons as a somewhat increased response for smaller q. In fact, the convergence of the electronic density response of hydrogen toward the UEG model for high densities can be seen most clearly in Fig. 15, where we show the partial local field factors θ_ab(q). In stark contrast to r_s = 2 and mainly to r_s = 3.23, the UEG model is in very good agreement with the true θ_ee(q) of hydrogen at r_s = 1. This has important implications for laser fusion applications and clearly indicates that UEG-based models such as the adiabatic local density approximation are appropriate over substantial parts of the ICF compression path. Finally, we note that the electron–proton local field factor θ_ep(q) has the same magnitude as θ_ee(q) in this case, whereas θ_pp(q) is somewhat smaller.

In this work, we have presented the first ab initio results for the partial density response functions of warm dense hydrogen. This has been achieved on the basis of direct PIMC simulations that are computationally very expensive owing to the fermion sign problem, but are exact within the given Monte Carlo error bars. Moreover, we have employed the recently introduced ξ-extrapolation technique84–88 to access larger system sizes; no finite-size effects have been detected for the wavenumber-resolved properties, in agreement with previous results for the UEG model under the same conditions.89 A particular advantage of the direct PIMC method is that it allows us to estimate all ITCFs F_ab(q, τ). First and foremost, this gives us direct access to the full wavenumber dependence of the static density response χ_ab(q) and consequently the local field factor θ_ab(q) from a single simulation of the unperturbed system. As an additional cross-check, we have also carried out extensive simulations of hydrogen where either the electrons or the protons are subject to an external monochromatic perturbation. We find perfect agreement between the direct perturbation approach and the ITCF-based method in the linear response regime, as expected.

In addition to the anticipated impact on future investigations of WDM that is outlined below, the study presented in this paper is of great interest in its own right and has provided new insights into the complex interplay of the electrons and protons in different regimes. We repeat that both F_ee(q, τ) and χ_ee(q) can be obtained from XRTS measurements, and our results thus constitute unambiguous predictions for experiments with ICF plasmas and hydrogen jets. In particular, we have shown that χ_ee(q) is highly sensitive to the electronic localization around the protons; see Fig. 16, where we compare χ_ee(q) for hydrogen and the UEG model at Θ = 1 and three different densities. This effect is particularly pronounced for small wavenumbers, which can be probed in forward scattering geometries,51 and it is directly related to the important concept of effective ionization. At the same time, we stress that the latter is strictly speaking, ill-defined and ambiguous, whereas the reported impact of electronic localization on χ_ee(q) constitutes a well-defined physical observable both in experiments and in simulations. In terms of physical parameters, we have found that the electrons in hydrogen behave very differently from the UEG model for r_s = 2 and even more so for r_s = 3.23, while the UEG model is appropriate for strongly compressed hydrogen at r_s = 1. Finally, we have reported, to the best of our knowledge, the first reliable results for the local field factor θ_ab(q) that is directly related to the XC kernel from LR-TDDFT calculations.

We are convinced that our study opens up a wealth of opportunities for significant future work, and helps to further lift WDM theory onto the level of true predictive capability. We give a non-exhaustive list of particularly promising projects:1.The electron–electron density response function χ_ee(q) is ideally suited for the model-free interpretation of XRTS measurements of WDM. Specifically, we propose to first infer the temperature from the exact ITCF-based thermometry method introduced in Refs. 53 and 110 and to subsequently carry out a set of PIMC simulations for χ_ee(q) over a relevant set of densities. Matching the PIMC result to the experimental result for χ_ee(q), the inverse moment of S_ee(q, ω), then gives one model-free access to the density. This PIMC-based interpretation framework for XRTS experiments will provide new insights into the equation of state of WDM, with important implications for astrophysical models and laser fusion applications. Additionally, we note that computing χ_ee(q) does not require dynamical simulations or dynamic XC kernels, and therefore might be suitable for approximate methods such as DFT.2.The species-resolved local field factors θ_ab(q), and in particular the electron–electron XC kernel, constitute key input for a gamut of applications, including the estimation of thermal and electrical conductivities,⁶¹ the construction of effective potentials,^66,151,152 and the estimation of the ionization potential depression.⁶⁵ Two particularly enticing applications concern the construction of nonlocal XC functionals for DFT simulations based on the adiabatic connection formula and the fluctuation–dissipation theorem,⁶⁸ and LR-TDDFT simulations within the adiabatic approximation.^39,115 In fact, previous studies of the UEG^{72,73,124,144} have reported that with the use of a static XC kernel, it is possible to obtain highly accurate results for S_ee(q, ω) over substantial parts of the WDM regime. Extending these considerations for hydrogen and beyond are promising routes that will be explored in dedicated future works.3.Our quasi-exact PIMC results constitute a rigorous benchmark for computationally less expensive though approximate simulation methods, most importantly thermal DFT. Therefore, a dedicated comparative investigation for a real WDM system will give invaluable new insights into the range of applicability of available XC functionals, as well as guiding the development of new thermal functionals that are explicitly designed for application in the WDM regime.^43,46–49 In addition, the results presented for both F_ee(q, τ) and χ_ee(q) can be used to benchmark dynamic methods such as LR-TDDFT,^38,39 real-time TDDFT,¹⁵³ and indeed the popular but uncontrolled Chihara models.^51,54,58,934.A less straightforward, though highly rewarding endeavor is represented by the so-called analytic continuation of F_ee(q, τ), i.e., the numerical inversion of Eq. (20) to solve for the dynamic structure factor S_ee(q, ω). While such an inverse Laplace transform constitutes a notoriously difficult and, in fact, ill-posed problem,¹⁵⁴ this issue has been circumvented recently for the warm dense UEG based on stochastic sampling of the dynamic local field factor with rigorous constraints imposed on the trial solutions.^124,144 Finding similar constraints for warm dense hydrogen would open up the way for the first exact results for the DSF of real WDM systems, as well as for related properties such as the dynamic dielectric function, conductivity, and dynamic density response.^113,1145.Finally, we note that current direct PIMC capabilities allow for highly accurate simulations of elements up to Be.^87,88 It will be very interesting to see how the more complex behavior of such systems (e.g., double occupation of the K shell of Be for T = 100 eV^87,88) manifests in the different ITCFs and density response functions. Moreover, these considerations might be extended to material mixtures such as LiH, giving rise to additional cross-correlations that can be straightforwardly estimated by upcoming PIMC simulations.

Acknowledgment. This work was partially supported by the Center for Advanced Systems Understanding (CASUS), financed by Germany’s Federal Ministry of Education and Research (BMBF) and the Saxon State Government out of the State Budget approved by the Saxon State Parliament. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2022 Research and Innovation Program (Grant Agreement No. 101076233, “PREXTREME”). The views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. Computations were performed on a Bull Cluster at the Center for Information Services and High-Performance Computing (ZIH) at Technische Universität Dresden, at the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN) under Grant No. mvp00024, and on the HoreKa supercomputer funded by the Ministry of Science, Research and the Arts Baden-Württemberg and by the Federal Ministry of Education and Research.