Matter at extreme temperatures and pressures—commonly known as warm dense matter (WDM)—is ubiquitous throughout our Universe and occurs in astrophysical objects such as giant planet interiors and brown dwarfs. Moreover, WDM is very important for technological applications such as inertial confinement fusion and is realized in the laboratory using different techniques. A particularly important property for the understanding of WDM is given by its electronic density response to an external perturbation. Such response properties are probed in x-ray Thomson scattering (XRTS) experiments and are central for the theoretical description of WDM. In this work, we give an overview of a number of recent developments in this field. To this end, we summarize the relevant theoretical background, covering the regime of linear response theory and nonlinear effects, the fully dynamic response and its static, time-independent limit, and the connection between density response properties and imaginary-time correlation functions (ITCF). In addition, we introduce the most important numerical simulation techniques, including path-integral Monte Carlo simulations and different thermal density functional theory (DFT) approaches. From a practical perspective, we present a variety of simulation results for different density response properties, covering the archetypal model of the uniform electron gas and realistic WDM systems such as hydrogen. Moreover, we show how the concept of ITCFs can be used to infer the temperature from XRTS measurements of arbitrary complex systems without the need for any models or approximations. Finally, we outline a strategy for future developments based on the close interplay between simulations and experiments.

## I. INTRODUCTION

Warm dense matter (WDM) is an extreme state that is characterized by high temperatures ($T\u2208103\u2212108$ K), high pressures ($P\u223c1\u2212104$ Gbar), and densities in the vicinity of and even partially exceeding solid state densities ($n\u223c1021\u22121028$ cm^{–3}). These conditions are ubiquitous throughout Nature and occur in a host of astrophysical objects,^{1,2} see the overview in Fig. 1. Prominent examples are not only giant planet interiors,^{3,4} such as Jupiter in our solar system,^{5–9} but also exoplanets.^{10,11} Other natural realizations of WDM include brown dwarfs,^{6,12} white dwarf atmospheres,^{13,14} neutron star crusts,^{15,16} and the remnants of meteor impacts.^{17,18} Another prominent example of WDM is the core of Earth. It consists primarily of iron at a temperature of ∼6000 K and pressure of ∼300 GPa, which is considered at the cold end of the WDM parameter space. Understanding the response properties, such as the electrical and thermal conductivity in warm dense iron,^{19,20} is tied to geophysical dynamics within Earth's interior that generate its magnetic field.^{21} Experimental measurements and modeling are the subjects of very active investigations.^{22–25}

In addition to its fundamental importance for stellar objects, WDM is also highly relevant for a number of technological applications, such as the discovery of novel materials^{17,30–34} and hot-electron chemistry.^{35–37} A particularly important application is given by inertial confinement fusion^{38,39} as it is realized at the National Ignition Facility (NIF).^{40} On its pathway toward ignition, the fuel capsule traverses the WDM regime^{27} (see the black line in Fig. 1 that indicates the implosion path of a deuterium–tritium [DT] capsule^{27}).

Due to the high current interest in WDM, such extreme states are realized in large experimental facilities using different techniques, see the topical overview by Falk.^{41} Indeed, many ground-breaking results have been reported over the last few years, including the formation of nanodiamonds at extreme pressure^{17,30,42} and the study of the liquid-liquid phase transition in hot dense hydrogen.^{43–47} In addition, it can be expected that emerging capabilities at facilities such as NIF,^{40,48} LCLS at SLAC,^{18,49} the Sandia Z Pulsed Power Facility^{50–52} and the OMEGA laser^{53} in the USA, SACLA^{54} in Japan, and the European XFEL^{55} in Germany will open up new avenues for the study of matter at increasingly extreme densities and temperatures, with direct relevance to technological applications and laboratory astrophysics.^{13,56,57}

At the same time, it is important to note that a rigorous theoretical description of WDM is notoriously difficult.^{11,26,28} More specifically, WDM can be conveniently characterized in terms of three dimensionless parameters that are of the order of unity.^{28,58} First, the Wigner-Seitz radius $rs=d/aB$ is given by the ratio of the average interparticle distance *d* to the in Bohr radius $aB$, see the vertical blue lines in Fig. 1. Second, the degeneracy temperature $\theta =kBT/EF$ measures the thermal energy in units of the electronic Fermi energy $EF$,^{59} with $\theta \u226a1$ and $\theta \u226b1$ corresponding to the fully degenerate and semi-classical regime, respectively. Different values of *θ* are included as the green lines in Fig. 1. A third useful parameter is the classical coupling parameter of the electrons Γ that measures the ratio of the interaction to the kinetic energy on the level of a mean-field description;^{58} it is depicted by the solid red line in Fig. 1, and we also find $\Gamma \u223c1$ in the WDM regime.

Consequently, there are no small parameters that can serve as the basis for a suitable expansion.^{11,26} Nevertheless, the strong demand for theoretical models and simulations has sparked a surge of developments in the field. For example, *ab initio* quantum Monte Carlo (QMC)^{60} methods are capable of providing exact results for the static properties of the uniform electron gas (UEG)^{61–75}—the quantum version of the classical one-component plasma (OCP)^{28,59,76}—and light elements such as hydrogen^{77,78} over substantial parts of the WDM regime. Moreover, path-integral Monte Carlo (PIMC) simulations within the fixed-node approximation^{79}—also known as restricted PIMC (RPIMC) in the literature—have allowed Militzer and co-workers^{80–82} to go to heavier elements and material mixtures.^{83} These simulations constitute the basis for an extensive equation-of-state (EOS) database^{84} that can be used for a host of practical applications.

A second important tool for the *ab initio* simulation of WDM is the thermal density functional theory (DFT),^{85–89} which combines a computational cost that is generally less demanding compared to QMC methods with an often acceptable accuracy.^{90–93} In particular, first accurate parametrizations of the exchange-correlation (XC) free energy $fxc$ of the warm dense UEG^{28,69,94,95} allow for thermal DFT simulations of WDM on the level of the local density approximation (LDA) that take consistently into account the dependence of the XC functional on the temperature parameter *θ*. Such temperature effects can have a substantial impact on the DFT results for different observables,^{96–98} and the development of improved XC functionals for WDM calculations constitutes an active area of research.^{99–102} This is complemented by the development of efficient algorithms that allow one to extend DFT to higher temperatures,^{103–111} which are unfeasible for standard Kohn–Sham implementations.^{112} Moreover, there has been a surge of research activities in utilizing machine-learning techniques in the context of DFT.^{113} The most promising approaches for improving calculations in the WDM regime include machine-learning interatomic potentials that enable large-scale *ab initio* molecular dynamics calculations and actually replacing DFT calculations with machine-learning surrogate models.^{114} The development of accurate machine-learning interatomic potentials has become a broad field of research.^{115} They enable covering a large parameter space in temperature and pressure than is feasible with conventional *ab initio* molecular dynamics as has recently been demonstrated in several works.^{42,116–118} Replacing DFT all along in terms of machine-learning models is an avenue of research with great potential. Pioneering efforts on models^{119–121} have matured into production-level codes that accelerate DFT calculations for realistic materials at finite temperature and pressure^{122,123} and have recently enabled simulations at unprecedented system size.^{124}

In combination with complementary methods such as quantum hydrodynamics,^{125–130} quantum scattering theory,^{131,132} and kinetic models,^{133,134} these developments have given new insights into the physics of WDM, including the EOS,^{84,135,136} effective potentials,^{137–140} and a number of transport properties^{141} such as stopping power^{109,142–144} and electrical conductivity.^{133,145}

Many of the properties of WDM are encoded in its response to an external perturbation. Of particular interest are perturbations that couple to the electronic charge density, e.g., an incident electromagnetic field or the Coulomb field of an incident charged particle. The response to these probes contains information about a wide range of excitations, and many of which are sufficiently sensitive to the thermodynamic state of the WDM that they can be used to infer conditions like temperature^{146} or density.^{136,147,148} X-ray Thomson Scattering (XRTS)^{147,149–151} is a widely used diagnostic that enables these inferences (among others) through measurements of the intensity of hard x-rays that are inelastically scattered from WDM samples. The fact that hard x-rays are required is a reflection of the fact that the densities of WDM are sufficiently high that they are opaque to lower energy probes. Those same hard x-rays can be used to probe the state of these systems, while isochorically heat it on femtosecond time scales.^{145,152}

More specifically, XRTS experiments rely on measurements of the angle- and energy-resolved scattering intensity to probe the dynamic structure factor (DSF), $S(q,\omega )$, where **q** and *ω* are a wave vector and frequency, respectively. The DSF itself can be conveniently defined as the Fourier transform of the intermediate scattering function,^{147}

which is defined as

for homogeneous systems in equilibrium. For sufficiently weak probes in the linear response regime, the DSF is related to the linear density response function $\chi (q,\omega )$ through the fluctuation-dissipation theorem^{59} (FDT), Eq. (31). Kinematic constraints from the first Born approximation dictate that the scattering angle is related to the momentum transferred **q** to or from the sample in the scattering process. Similarly, the shift in energy upon scattering *ω* is the quantity of interest, rather than the absolute energy of the detected light, as this encodes information about the excitation or de-excitation of the sample. Thus, it is typical to parametrize the XRTS scattering intensity $I(q,\omega )$ in terms of the same momentum and energy transfer as those parametrizing $S(q,\omega )$ and $\chi (q,\omega )$, even though the salient experimental quantities (detection angle and photon energy) are different.

In fact, because the XRTS probe beam is only approximately monochromatic and the detector has energy-dependent sensitivity, the actually detected intensity is a convolution of the DSF and the combined source-instrument function $R(\omega )$,

Nevertheless with careful characterization of $R(\omega )$, information about $S(q,\omega )$ and $\chi (q,\omega )$ can be obtained from the measured intensity. Thus in addition to its practical utility as a multi-purpose diagnostic, XRTS is ultimately a powerful technique to directly access the density response function.

From a theoretical perspective, it is very convenient to express $\chi (q,\omega )$ as^{153}

with $\chi 0(q,\omega )$ is a known reference function, such as the density response of an ideal Fermi gas (also known as Lindhard function in the literature) in the case that $\chi (q,\omega )$ describes the density response of the UEG. In this case, the complete wave-vector- and frequency-resolved information about electronic XC effects is contained in the local field correction (LFC) $G(q,\omega )$. Indeed, setting $G(q,\omega )\u22610$ in Eq. (4) corresponds to the well-known random phase approximation (RPA), which describes the density response on a mean-field level.^{154} Consequently, the LFC constitutes key input for a host of applications, such as the construction of advanced XC functionals for DFT via the adiabatic-connection formula,^{99,155,156} the incorporation of electron–electron correlations into quantum hydrodynamics,^{125,126} or the interpretation of XRTS experiments.^{157–160} In addition, the LFC is formally equivalent to the XC kernel from linear response time-dependent DFT (LR-TDDFT);^{161–163} they are related by $Kxc(q,\omega )=\u22124\pi /q2G(q,\omega )$. Indeed, replacing $\chi 0(q,\omega )$ by the dynamic Kohn–Sham response function $\chi S(q,\omega )$ in Eq. (4) constitutes the basis for the LR-TDDFT simulation of real materials,^{162,164–166} see Sec. II A 2.

While the assumption of a linear response is fairly ubiquitous throughout quantum many-body theory, it holds, strictly speaking, only in the limit of an infinitesimal perturbation amplitude *A*. Recently, Dornheim *et al.*^{167} have presented the first *ab initio* PIMC simulation results for the nonlinear density response of the warm dense UEG. Such nonlinear effects^{168–177} become important, for example, for XRTS experiments with ultra-high intensities that can be realized at x-ray free-electron laser (XFEL) facilities via the novel seeding technique.^{178} Moreover, they are an important ingredient to the construction of effective pair potentials at small distances^{139,179,180} and stopping power calculations for slow and/or heavy particles.^{181–183} Finally, it has been shown that nonlinear effects depend much more sensitively on system parameters like the temperature *T*,^{167,175} which makes their experimental observation a potentially interesting new method of diagnostics.^{184} Describing the dynamics of a quantum many-body system due to a time-dependent external perturbation beyond the linear response regime is achieved by the real-time formulation of time-dependent DFT (RT-TDDFT),^{185} which will be introduced later.

Finally, we note that even the interpretation of XRTS experiments with comparably low scattering intensity such that any nonlinear effects can safely be neglected generally relies on a number of model assumptions. More specifically, the most widely used model is the Chihara decomposition,^{147,148,151,186} which is a chemical model assuming a clear distinction between bound and free electrons. This assumption is questionable in the WDM regime^{77} and constitutes a de facto uncontrolled approximation in practice. Therefore, previous EOS measurements^{14,187,188} can substantially depend on the employed model and do not necessarily constitute a reliable benchmark for the development of new theoretical approaches such as TDDFT,^{189} or reliable input for other applications. To overcome this unsatisfactory situation, Dornheim *et al.*^{146,190–192} have very recently suggested switching from the usual frequency-representation to the imaginary-time domain. This gives one direct access to the physical properties of a given system of interest and allows one to extract, for example, the temperature without any models, approximations, or simulations. In addition, such imaginary-time correlation functions (ITCF) naturally emerge within Feynman's imaginary-time path-integral picture of statistical mechanics^{193} and constitute an important link between simulations,^{194,195} physical insights,^{190,192} and experimental measurements.^{146,191}

In this work, we attempt to give a coherent picture of the current state-of-the-art with respect to our understanding of the electronic density response of WDM. While we put more emphasis on a number of recent promising developments with respect to theoretical models and different simulation techniques, we also indicate the respective connection to relevant experiments. Indeed, although the notion of evaluating experimental data in the imaginary-time domain might seem to be somewhere between impractical to outright outlandish, we demonstrate the simplicity of the idea using an actual XRTS dataset.

The paper is organized as follows. Sec. II contains the relevant theoretical background, starting with an introduction of the most relevant numerical methods in Sec. II A, namely, path-integral Monte Carlo (Sec. II A 1), thermal DFT (Sec. II A 2), and real-time time-dependent DFT (RT-TDDFT) (Sec. II A 3). Moreover, we briefly touch upon nonequilibrium Green's functions (NEGF) in Sec. II A 4. In addition, we discuss important relations within linear response theory (Sec. II B), its relationship to imaginary-time correlation functions (Sec. II C), and finally some basic concepts for a nonlinear theory of the electronic density response (Sec. II D). Section III is devoted to the discussion of a gamut of simulation results both for a uniform electron gas and for real WDM systems. In particular, we start our discussion with results for the linear density response in the static limit in Sec. III A, including a discussion of PIMC simulations of the UEG (Sec. III A 1) and recent PIMC and thermal DFT simulations of warm dense hydrogen (Sec. III A 2). In Sec. III B, we extend these considerations to the fully dynamic case, again starting with a PIMC-based investigation of the UEG (Sec. III B 1) followed by TDDFT calculations for hydrogen, iron, and aluminum in Sec. III B 2. The discussion of simulation results is completed in Sec. III C, where we present new results for the nonlinear density response of the warm dense UEG. In Sec. IV, we connect our theoretical framework to experimental observations by re-interpreting the XRTS measurement of plasmons in warm dense beryllium by Glenzer *et al.*^{202} in the imaginary-time domain.^{146,191} This paper is concluded in Sec. V, where we discuss how our understanding of the electronic density response of WDM can be improved in future works by combining new developments in *ab initio* simulations with new experimental setups.

## II. THEORY

We consider an *N*-particle system that is described by the sum of an unperturbed Hamiltonian $H\u03020$ and an external single-particle perturbation that we give in the most general form,

Here, $K\u0302,\u2009W\u0302$, and $V\u0302pot$ are the kinetic, interaction, and external Coulomb energy (e.g., due to a configuration of ions) of the electrons, respectively. Moreover, $V\u0302q,\omega ,A$ corresponds to an additional external perturbation with **q** and *ω* being the corresponding wave vector and frequency, and $A(t)=Af(t)$ is a real perturbation amplitude that can explicitly depend on the time.^{203} We note that both $V\u0302q,\omega ,A$ and $V\u0302pot$ contribute to the total external single-electron potential energy $V\u0302ext=V\u0302pot+V\u0302q,\omega ,A$. In Eq. (6), we introduced a parameter *γ* that is either 0 or 1 to distinguish between static and monochromatic excitation scenarios. An overview of the different physical situations that have been treated with established analytical and computational approaches is presented in Table I. Note that we assume Hartree atomic units throughout unless indicated otherwise.

. | f(t)
. | A_{0}
. | γ
. | Regime . | Approach . |
---|---|---|---|---|---|

1 | 1 | Small | 1 | Dynamic linear response | LR-TDDFT (Sec. II A 2), PIMC (Sec. II A 1) + analytic continuation (Sec. II C) |

2 | 1 | Large | 1 | Dynamic nonlinear response | Dynamic nonlinear response theory (Ref. 175) |

3 | 1 | Small | 0 | Static linear response | PIMC (Sec II A 1), KS-DFT (Sec. II A 2), dielectric theories (Refs. 163 and 196–201) |

4 | 1 | Large | 0 | Static nonlinear response | PIMC (Sec. II A 1), KS-DFT (Sec. II A 2) |

5 | $A0\delta (t)$ | Small | 0 | Weak kick | Real-time TDDFT (RT-TDDFT) (Sec. II A 3), NEGF (Sec. II A 4) |

6 | $A0\delta (t)$ | Large | 0 | Strong kick | RT-TDDFT (Sec. II A 3), NEGF (Sec. II A 4) |

. | f(t)
. | A_{0}
. | γ
. | Regime . | Approach . |
---|---|---|---|---|---|

1 | 1 | Small | 1 | Dynamic linear response | LR-TDDFT (Sec. II A 2), PIMC (Sec. II A 1) + analytic continuation (Sec. II C) |

2 | 1 | Large | 1 | Dynamic nonlinear response | Dynamic nonlinear response theory (Ref. 175) |

3 | 1 | Small | 0 | Static linear response | PIMC (Sec II A 1), KS-DFT (Sec. II A 2), dielectric theories (Refs. 163 and 196–201) |

4 | 1 | Large | 0 | Static nonlinear response | PIMC (Sec. II A 1), KS-DFT (Sec. II A 2) |

5 | $A0\delta (t)$ | Small | 0 | Weak kick | Real-time TDDFT (RT-TDDFT) (Sec. II A 3), NEGF (Sec. II A 4) |

6 | $A0\delta (t)$ | Large | 0 | Strong kick | RT-TDDFT (Sec. II A 3), NEGF (Sec. II A 4) |

We are here interested in response functions whose most general form explicitly depends on two spatial locations and two times $\chi (r1,t1;r2,t2)$ for inhomogeneous and nonequilibrium situations. However, for sufficiently homogeneous equilibrium systems, only the difference variables remain important, and after the Fourier transform of $r=r1\u2212r2$ and $t=t1\u2212t2$, one usually studies $\chi (q,\omega )$.

In the limit of small *A*, linear response theory (LRT)^{59} becomes valid, which leads to a number of simplifications and useful relations that are discussed in Sec. II B. A more general nonlinear theory of the electronic density response is discussed in Sec. II D. While it is often convenient to restrict oneself to a static perturbation amplitude $A(t)\u2261A$ in Eq. (5), the more general form is useful to study the response of a given system to a finite perturbation pulse with real-time methods, which is discussed in more detail in Sec. II A 4.

From a practical perspective, we note that a gamut of numerical methods for the description of WDM systems has been presented in the literature, including molecular dynamics, classical Monte Carlo simulations, average-atom models^{204–209} as well as integral equation theory approaches within the hypernetted-chain approximation for effective quantum potentials,^{210–213} dielectric formalism schemes,^{163,196–200,214–219} and quantum hydrodynamics.^{125,126,129,130,220}

### A. Numerical methods

#### 1. Path-integral Monte Carlo

The *ab initio* path-integral Monte Carlo (PIMC) approach^{221–224} is one of the most successful methods for the description of quantum degenerate many-body systems. It is based on an evaluation of the partition function in coordinate representation, which, in the case of an electron gas in the canonical ensemble (inverse temperature $\beta =1/kBT$, volume $V=L3$, and number density $n=N/V$ are constant) with $N\u2191$ and $N\u2193$ spin-up and spin-down electrons, is given by

Here, the summation over all elements *σ ^{i}* of the respective permutation group $SNi$, with $\pi \u0302\sigma i$ being the corresponding permutation operator, ensures the exact anti-symmetry with respect to the exchange of coordinates of identical fermions. Since a detailed introduction to the PIMC method has been presented elsewhere,

^{221}we here restrict ourselves to a brief sketch of the main idea. The central obstacle with respect to a direct evaluation of Eq. (7) is the non-commutability of the kinetic ($K\u0302$) and potential ($V\u0302=W\u0302+V\u0302ext$) contributions to the full Hamiltonian $H\u0302=K\u0302+V\u0302$, which renders the matrix elements of the density operator $\rho \u0302=e\u2212\beta H\u0302$ generally unfeasible. To overcome this issue, one can exploit the well-known semi-group property of the density operator, which eventually leads to a summation over

*P*particle coordinates $Ri$. Crucially, each of the corresponding

*P*density matrices has to be evaluated at

*P*-times the original temperature, which allows for the introduction of suitable factorization schemes;

^{225–227}the associated factorization error can be estimated from the Baker–Campbell–Hausdorff formula

^{228}and vanishes in the limit of large

*P*.

In addition, we note that each high-temperature factor can be interpreted as a propagation in the imaginary time $t=\u2212i\u210f\tau $ over a discrete time step $\tau =\beta /P$. Consequently, each quantum particle is represented as a path along the imaginary time, which can be mapped onto an ensemble of classical ringpolymers; this is the origin of the celebrated classical isomorphism by Chandler and Wolynes.^{229} Moreover, PIMC gives one direct access to different imaginary-time correlation functions, such as the imaginary-time version of the intermediate scattering function defined in Eq. (39). The basic idea of the PIMC method is to stochastically generate all possible path configurations according to their respective contribution to the partition function *Z* based on the Metropolis algorithm.^{230} For indistinguishable quantum particles—bosons or fermions—this also requires the sampling of all possible permutation topologies,^{231} which can be achieved efficiently using different implementations^{72,232} of the worm algorithm by Boninsegni *et al.*^{233,234}

In the case of bosons (or hypothetical distinguishable quantum particles, which are sometimes referred to as Boltzmannons in the literature^{235,236}), the sign function in Eq. (7) is always positive, and the PIMC method allows for the quasi-exact simulation of up to $N\u223c104$ particles, which has given important insights into phenomena such as superfluidity.^{237–239} For fermions, such as the electrons that are a key constituent of WDM, the sign function changes its sign for every pair exchange. The resulting cancelation of positive and negative terms is the root cause of the notorious fermion sign problem,^{240–243} which leads to an exponential increase in the required compute time upon increasing the system size *N* or decreasing the temperature *T*. Therefore, the direct PIMC method that is used in the present work is limited to specific parts of the WDM parameter space ($\theta \u22730.5$) and to light elements such as hydrogen.^{78}

For completeness, we note that the sign problem can formally be avoided by imposing restrictions on the nodal structure of the thermal density matrix; this restricted PIMC method developed by Ceperley^{79} is exact if the true nodal structure of the quantum many-body system of interest was known. In practice, however, one has to introduce approximations, and often the nodal structure of an ideal Fermi gas at the same conditions is used.^{79,244} First, this fixed-node approximation introduces an uncontrolled error that can be of the order of 10% in the XC energy of a UEG at high density and low temperature,^{68} whereas it is less pronounced in the momentum distribution function around *θ* = 1.^{72,73} Second, the nodal restrictions break the imaginary-time translation invariance, which prevents the direct estimation of imaginary-time correlation functions such as $F(q,\tau )$, see Sec. II C.

#### 2. Thermal density functional theory

Let us consider an (time-independent) electronic Hamiltonian of the form

with $K\u0302$ being the kinetic part, $W\u0302$ containing the electron–electron interaction, and $V\u0302ext$ an external single-particle potential. In particular, $V\u0302ext$ takes into account the electron–ion interaction when a snapshot of ions is considered, which is a common practice when DFT calculations are coupled to molecular dynamics (MD) simulations of classical ions within the Born–Oppenheimer approximation.^{11} Furthermore, $V\u0302ext$ can also include the monochromatic external perturbation that is used to obtain the electronic density response, cf. Eq. (6) in the limit of $\omega \u21920$ and $f(t)\u22611$. The basic idea of DFT is to express the ground-state energy of a many-electron system defined by Eq. (8) in terms of the many-electron density $n(r)$. At finite temperature, this is achieved by Mermin's generalization^{85} of the Hohenberg–Kohn theorems.^{245} This enables a thermodynamic description of many-electron systems within the context of DFT by the grand-canonical potential as a density functional,

where

is the universal functional that includes contributions from the kinetic energy, entropy $S$, and electron–electron interaction. Note that $FT[n]$ itself is minimized over the statistical operator $\Gamma \u0302$ defined by the many-particle states of the Hamiltonian in Eq. (8). Then, the grand-canonical potential of a many-electron system in thermal equilibrium is found by minimization

Minimizing $\Omega [n]$ over trial densities requires knowledge of all terms in Eq. (9) as a functional of the density. This is, however, not straightforward because explicit density functionals of these terms (besides the trivial term $Vext[n]$) are not known.

A practical solution to this problem is the KS approach.^{246} Here, one defines an auxiliary system of noninteracting particles described by a set of noninteracting single-particle Schrödinger equations,

where $\varphi \alpha ,k(r)$ and $\u03f5\alpha ,k$ denote the KS orbitals and corresponding eigenvalues with band index *α* and Bloch wave number $k$, respectively. The KS approach allows writing the grand-canonical potential as

where $KS[n]$ denotes the KS kinetic energy, $SS[n]$ the KS entropy, $U[n]$ the Hartree energy, and $\Omega XC$ the XC free energy. The KS system is constructed such that the density it yields

is identical to the many-electron density of the corresponding interacting many-electron system defined in Eq. (8). Due to their nature as effective single-particle states, the KS orbitals ${\varphi \alpha ,k}$ are populated according to the Fermi–Dirac distribution

with $\beta =1/kBT$ and *μ* is the chemical potential.^{59} In combination, Eqs. (12) and (15) give one direct access to $KS[n]=\u22121/2\u2211\alpha ,k\u222bdr\u2009\varphi \alpha ,k*(r)\u22072\varphi \alpha ,k(r)$ and $SS[n]=\u2212\u2211\alpha ,k(f\alpha ,k\u2009ln\u2009(f\alpha ,k)+(1\u2212\u2009f\alpha ,k)\u2009ln\u2009(1\u2212f\alpha ,k))$, whereas all non-ideal contributions to the full kinetic energy *K* and entropy $S$ are contained in the XC functional. The equality of the density in Eq. (14) with the true many-electron density is achieved by the KS potential defined as $vS[n]=vext(r)+\u2009vH[n](r)+vXC[n](r)$ in terms of the external potential $vext(r)$, the Hartree potential $vH[n](r)=\delta U[n]/\delta n(r)$, and the XC potential $vXC[n](r)=\delta \Omega XC[n]/\delta n(r)$. In particular, the XC contribution $\Omega XC[n]$ contains the full information about many-body correlations and, therefore, would require the exact solution of the original many-electron problem, which is not feasible. In practice, $\Omega XC[n]$, therefore, has to be approximated. Specifically, the particular choice of the XC functional substantially influences the accuracy of a DFT calculation,^{247} which makes both the benchmarking of existing functionals^{90–93,255} and the construction of novel, more sophisticated approximations to $\Omega XC[n]$ indispensable.

We note that the KS orbitals by themselves should be viewed not as physical quantities but as auxiliary properties that are connected to the energy and density of a system. At the same time, the ${\varphi \alpha}$ constitute an important ingredient to a number of other applications, such as the incorporation of nonlocality into quantum hydrodynamics via an *ab initio* Bohm potential term.^{129,220} A particularly important application of the KS orbitals that have been obtained from an equilibrium DFT calculation is given by linear response time-dependent DFT (LR-TDDFT), which is based on the KS response function^{162,248} given by

Here, the parameter $0<\eta \u226a1$ ensures the retardation, the sum runs over the different eigenvalues *ϵ* and Fermi functions *f* at different momenta $k$ and $k+q$, *ω* is the energy of the excited mode, and $G,G\u2032$ are the reciprocal lattice vectors. We are here mainly interested in the macroscopic response functions, i.e., setting $G=G\u2032=0$.

The physically meaningful dynamic electronic linear density response of a system of interest is then given by

with $Kxc(q,\omega )$ being the XC kernel mentioned in the discussion of Eq. (4) above. Since the LHS of Eq. (17) is a well-defined physical observable, the kernel $Kxc(q,\omega )$ strongly depends on the particular form of the KS response function $\chi S(q,\omega )$. In other words, the true XC kernel has to depend both on the material and on the XC functional that has been used for the computation of the ${\varphi \alpha}$. From a theoretical perspective, the XC kernel is readily defined in terms of a double functional derivative of the XC free energy;^{163} this is difficult to evaluate in practice and, to our knowledge, only possible for simple XC functionals in the case of bulk systems. Consequently, accurate and dynamic XC kernels exist only for model systems such as the UEG.^{185,249–252}

An alternative approach to the computation of the XC kernel from the second order functional derivatives is based on the perturbation of the system by an external harmonic field to measure the density response and from that to extract the XC kernel by inverting Eq. (17). This method was used by Moroni *et al.*^{253} to calculate the static density response function and the LFC of the UEG in the ground state from quantum Monte Carlo simulations. Recently, it was used within KS-DFT by Moldabekov *et al.* to compute the static XC kernel of the UEG and warm dense hydrogen on the basis of various XC functionals across Jacob's Ladder.^{254,255} Some illustrative results from these studies are presented in Sec. III A 2.

We note that setting $Kxc(q,\omega )\u22610$ in Eq. (17) is often denoted as RPA in the DFT literature but does not correspond to a mean-field description as the KS orbitals automatically contain a certain amount of information about electronic XC effects due to the employed XC functional. Finally, having the LR-TDDFT result for $\chi (q,\omega )$ gives one direct access to a number of material properties, such as the dynamic structure factor $S(q,\omega )$, see Eq. (31).

#### 3. Real-time time-dependent DFT (RT-TDDFT)

The real-time approach to TDDFT (RT-TDDFT)^{185,256,257} is a computationally efficient method for studying nonequilibrium electronic dynamics in first principles models of materials with an explicit treatment of the electron–ion interaction (usually within the Born–Oppenheimer approximation or sometimes using Ehrenfest molecular dynamics) and large basis sets. The equations of motion in RT-TDDFT are the time-dependent Kohn–Sham equations,

where we have $vS[n](r,t)=vext(r,t)+vH[n](r,t)+vxc[n](r,t)$. In what follows, we shall abbreviate the effective Hamiltonian as $H\u0302S(t)=\u2212\u22072/2+vS(r,t)$. This is typically posed as an initial value problem in which $\varphi \alpha ,k(r,t)$ is known for some *t* = *t _{i}*, typically from the solution to a static Mermin–Kohn–Sham DFT problem. The time-dependent density is

and these occupation factors are typically consistent with the initial condition and held constant over the course of the dynamics. Among the three terms in *v _{S}*, $vext$ describes the attractive electron–ion interaction as well as impressed perturbations that drive the dynamics,

*v*is the time-dependent Hartree potential, and

_{H}*v*is the time-dependent XC potential that describes electronic exchange and correlation. The creation of approximations to the exact time-dependent Kohn–Sham potential that are both accurate and computationally efficient is one of the central problems in RT-TDDFT. The practical trade-off between the cost and accuracy of available and prospective approximations is the central factor in determining whether RT-TDDFT is the right method for a particular calculation.

_{xc}While RT-TDDFT is capable of approximating these dynamics owing to arbitrary excitations [i.e., forms of $vext(r,t)$], it can also capture linear response functions in the limit that $vext(r,t)$ is proportional to $\delta v0f(t)$, where $\delta v0$ is a constant that determines the strength of the perturbation and *f*(*t*) is an arbitrary smooth function whose Fourier transform has support on the energies over which the linear response function is desired. As in LR-TDDFT, the response function computed using RT-TDDFT captures collective effects^{258,259} that are absent in approaches, in which the Kubo–Greenwood (KG) formula is directly applied to the Kohn–Sham orbitals.^{260–262} Furthermore, the computational details of these formulations differ in such a way that RT-TDDFT can be made asymptotically less expensive than LR-TDDFT. To understand when and how RT-TDDFT is most efficient, we first need to describe typical approaches to the numerical solution of Eq. (18).

The time integration of Eq. (18) requires the implementation of a numerical approximation to the exact time-ordered unitary propagator that translates any KS orbital in time,

and $T$ is the time-ordering operator. One approach to approximating this is first-order Trotterization of Eq. (20a) into *N* time steps of width $\Delta t=(tf\u2212ti)/N$, over which the effective Hamiltonian is approximated as a constant at the midpoint of each step, eliminating the need for time ordering,

While conceptually straightforward, the practical implementation of this approach requires the numerical approximation of a matrix exponential and, thus, the diagonalization of $H\u0302S$ at each time step or some other more efficient method for numerically effecting this operation.

Naive approaches to these operations have a cost that is cubic in the number of KS orbitals, which grows rapidly with temperature in and beyond the warm dense regime. The parallel scalability of these approaches is also limited by a conflict between the optimal data distribution for effecting the product of $H\u0302S$ and a trial orbital, and orthogonalization of those orbitals, i.e., the optimal distribution for one is suboptimal for the other. It turns out that it is possible to avoid both of these costs in RT-TDDFT, and, thus, it can be made to be more efficient and scalable^{259,263–265} than implementations that rely directly on dense linear algebra, including many implementations of LR-TDDFT. The precise conditions under which an RT-TDDFT linear response calculation is more efficient than an LR-TDDFT calculation are sensitive to many details of the calculation, but a large number of thermally occupied orbitals common in WDM calculations tend to favor RT-TDDFT as temperature increases.

While a scalar perturbation of the form $vext(r,t)=v0\u2009exp\u2009(iq\xb7r)f(t)$ can be used to compute the DSF at non-zero $q$,^{189,225} the calculation of the conductivity in the $q\u21920$ limit requires the extension of RT-TDDFT to vector perturbations. More specifically, these vector perturbations are used to simulate a spatially uniform electric field, $E(t)=\u2212(1/c)(\u2202A/\u2202t)$, impressed over the entire supercell.^{266,267} While the strength of this perturbation need not be weak in RT-TDDFT, as with the DSF, one can implement an LR-TDDFT calculation with a sufficiently weak $E(t)$ to compute the conductivity $\sigma (\omega )$. Then, rather than the time-dependent electronic charge density, the spatially uniform component of the time-dependent current density is related to the response function of interest. Calculations of the conductivity are realized using a microscopic version of Ohm's law,^{23,24}

Here, the current density (*J*) and electric field (*E*) are vectors, while the conductivity (*σ*) is a tensor. The dielectric tensor is related to the conductivity tensor via

Figure 2 illustrates the methodology starting with a weak sigmoidal pulse applied in the **z** direction. The induced current density along **z** is shown in the left panel, and the frequency-dependent conductivity obtained using the Fourier transform of the induced current density is shown in the right panel. Note that the values for $jz(t)$ in Fig. 2 are in atomic units.

#### 4. Nonequilibrium Green's functions

The nonequilibrium Green's functions (NEGF) theory is a powerful approach to treat electronic correlations, quantum and spin effects as well as dynamical screening in non-ideal quantum systems, e.g., Refs. 154 and 268. There exist many applications to the warm dense uniform electron gas and to dense fully ionized plasmas.^{268–270} Being computationally very demanding, many studies of real materials have been restricted to ground-state properties, e.g., Refs. 271 and 272, providing important benchmarks for DFT, where QMC results are not available. On the other hand, an accurate analysis of the behavior of electrons out of equilibrium, including strong excitation and the short-time dynamics, is possible but has been performed mostly for model systems, such as the UEG or lattice systems.^{273,274} Thus, the NEGF theory is complementary to both PIMC and DFT.

The NEGF theory is formulated in second quantization (Fock space), starting from a complete set of single-particle orbitals, ${\varphi i}$ and the associated creation and annihilation operators, $c\u0302i\u2020$ and $c\u0302j$. The central quantity is the single-particle Green's function, which is an ensemble average of the Heisenberg form of the field operators,

where $TC$ is the time-ordering operator that is analogous to the one in Eq. (20a). Here, the subscript *C* indicates the use of the Keldysh time contour $C$^{275} that allows the extension of the power of Feynman diagrams to arbitrary nonequilibrium situations (for details, see the text books^{276,277}) The angular brackets denote the expectation value that is computed either with the *N*-particle wave function (pure state) or N-particle density operator (mixed state) of the non-excited system.

The Green's function (24) gives access to all one-particle observables, such as the density matrix, which follows from the time-diagonal Green's function, $nji(t)=\u210figij(t,t+\u03f5)$, where *ϵ* is an infinitesimal positive constant assuring the proper ordering of the two field operators. A particular advantage of the NEGF theory arises from the two-time structure of the function *g*: its values for different time arguments give direct access to the spectral function, the density of states, and the interaction energy of the system.

A special case is the coordinate representation, which will be used in Sec. II D. There, a special notation for the field operators is used, $ci(t)\u2192\psi (1)$, where $1={r,t1}$ (the spin index is suppressed). Then, the single-particle space and time-dependent density is obtained from the space and time-diagonal elements of Green's function (24),

The equations of motion for the NEGF are the Keldysh–Kadanoff–Baym equation (KBE, to be supplemented by the adjoint equation),

where $hHF$ contains kinetic, potential, and mean-field (Hartree–Fock) energy contributions, whereas correlation effects are included in the self-energy $\Sigma cor[g]$, which is a functional of the one-particle NEGF.

The KBE is analogous to the time-dependent Kohn–Sham equations of TDDFT, Eq. (18). The main differences are: (1) the KBE are equations for a density matrix that is similar to a product of two KS orbitals. (2) The energy terms are grouped slightly differently than in TDDFT—the Hartree mean field and exact exchange are contained in $hHF$. (3) The exchange-correlation potential $vxc[n]$ is replaced by the correlation self-energy $\Sigma cor$, which depends on *g* rather than on the density and, therefore, depends on two times. Thus, it includes memory effects. (4) Similar to *v _{xc}* in the case of DFT, $\Sigma cor$ is the only approximation of the NEGF theory. Many systematic approaches exist for the self-energy, most importantly Feynman diagram methods that include the perturbation theory and partial resummations that allow the systematic inclusion of dynamical screening effects and bound states, e.g., Refs. 278 and 279. A comparison of the known approximations for

*v*and $\Sigma cor$ is a topic of active research as it allows for benchmarks between the two approaches and offers the derivation of novel approximations for both, DFT and NEGF.

_{xc}For the electronic density response of warm dense matter, the NEGF theory offers a variety of approaches. First, it is the derivation of equilibrium density response functions, which are obtained by functional derivation of *g*, as will be demonstrated in Sec. II D. This method allows one to systematically derive not only the linear density response but also nonlinear generalizations. Second, there exists a straightforward real-time approach to the linear and nonlinear density response that was developed by Kwong and Bonitz.^{203} To this end, an external potential is added to the general Hamiltonian (5) that is of the form $Vext(r,t)=U0(t)\u2009cos\u2009(qr)$. This potential imposes a density modulation to the system, and the short temporal duration of the pulse (“kick”) translates into a broad range of excitation energies. The solution of the KBE (26) for a UEG excited with this short pulse excites a monochromatic density perturbation that decays with the plasmon frequency for that *q* as shown in Fig. 3. The density fluctuation $\delta n(q,t)$ follows directly from the perturbation of the NEGF, via Eq. (25),

and is, to the lowest order (linear response), proportional to the excitation. In equilibrium, the density response function is stationary, $\chi (q,t,t\u2032)=\chi (q,t\u2212t\u2032)$, and Eq. (27) can be solved in Fourier space for *χ*: $\chi (q,\omega )=\delta n(q,\omega )/U0(\omega )$, which gives access to all linear response functions, including the dielectric function and the dynamic structure factor, as discussed in Sec. II B.

A remarkable feature of real-time NEGF computations of the density response is that already fairly simple approximations for the self-energy $\Sigma cor$ yield high level results for the density response functions. The relation to the first (equilibrium) approach which requires to solve the Bethe–Salpeter equation (BSE) for the two-particle Green's function with a kernel Ξ^{268} is given by Kwong and Bonitz^{203} and Bonitz *et al.*,^{280}

The diagrams contributing to Ξ, if the KBE is solved on the level of second Born self-energies, are shown in the bottom of Fig. 3. Thereby, use of conserving self-energies automatically yields sum rule preserving approximations for Ξ. These attractive properties of NEGF should also be fulfilled for other time-dependent approaches, including real-time TDDFT offering interesting opportunities for the treatment of correlation effects in WDM.

Finally, we note that the real-time NEGF approach can be straightforwardly extended to the nonlinear response of the UEG as demonstrated in Ref. 203.

### B. Linear response theory

At this point, let us set the perturbation amplitude in Eq. (5) to be a constant (except for an infinitesimal damping that assures the vanishing of the perturbation for $t\u2192\u2212\u221e$), which leads to the harmonically perturbed Hamiltonian,^{167,253,281}

In the limit of small *A*, the linear response theory (LRT)^{59} becomes valid, and the induced density is fully described by the simple linear relation,

where $\varphi q,\omega $ is the Fourier transform of the cosine perturbation, given by a sum of delta functions at positive and negative values of the wave vector **q** and the frequency *ω*. In other words, the system only exhibits a non-zero response at the wave vector and frequency of the original perturbation within LRT for homogeneous systems. For completeness, we note that this remains approximately true for weakly inhomogeneous systems, see Ref. 255 for a recent discussion. Effects such as the excitation of higher harmonics^{175} or mode-coupling between multiple perturbations^{177} are then exclusively nonlinear effects and are discussed in Sec. II D.

For now, we shall postpone the discussion of the theoretical estimation of the linear density response function $\chi (q,\omega )$ and instead focus on its utility for the description of WDM and beyond. In this regard, a key relation is given by the fluctuation-dissipation theorem (FDT),^{59}

where *n* denotes the average density. Equation (31) gives a direct and unique relation between the density response and the dynamic structure factor $S(q,\omega )$. In fact, $S(q,\omega )$ fully determines $\chi (q,\omega )$, as the real and imaginary parts of the latter are connected by the well-known Kramers–Kronig relations.^{59}

From a practical perspective, the DSF constitutes the central quantity in XRTS experiments with WDM;^{147,149,151} the measured scattering intensity $I(q,\omega )$ is given by the convolution of the DSF with the combined source and instrument function $R(\omega )$, see Eq. (3), and the wave vector **q** is determined by the scattering angle. Therefore, the availability of an accurate model for either $\chi (q,\omega )$ [subsequently evaluating Eq. (31)] or directly for $S(q,\omega )$ allows one to compare theoretical calculations with experimental measurements. Here, LR-TDDFT (Sec. II A 2) is an example for the first route, whereas approximate Chihara models^{147,148,151,186} that decompose the full DSF into contributions from bound and free electrons and their respective transitions directly work with $S(q,\omega )$. In practice, these models can easily be convolved with $R(\omega )$ and are used to infer a-priori unknown system parameters such as the temperature. Evidently, the quality of a, thus, inferred EOS significantly depends on the level of accuracy of the employed approximation. Very recently, Dornheim *et al.*^{146,190} have proposed to avoid this limitation by switching to the imaginary-time domain, which is discussed in Sec. II C.

A second highly important application of LRT stems from the combination of the FDT [Eq. (31)] with the adiabatic-connection formula, which we describe in the following. As an initial step, we note that the DSF gives straightforward access to the static structure factor (SSF) $S(q)$ via

For completeness, we point out that the SSF contains the same information as the pair correlation function $g(r)$, and they are related by Hansen and McDonald,^{282}

As the next step, we can use either $S(q)$ or $g(r)$ to estimate the interaction energy *W*, for example,^{282}

where $n(2)(r1,r2)=n(r1)n(r2)g(r1\u2212r2)$ denotes the two-particle density. Note that, for the classical and quantum OCP, owing to the rigid charge neutralizing background, the substitution $n2(r1,r2)\u2192n2(r1,r2)\u2212n(r1)n(r2)$ is implied, which is necessary for convergence.^{283} Finally, the adiabatic-connection formula relates Eq. (34) with the XC contribution $fxc$ to the free energy,^{59}

In other words, integrating over the interaction energy of a system where the interaction contribution to the Hamiltonian has been re-scaled by $\lambda \u2208[0,1],\u2009H\u0302\lambda =H\u03020+\lambda W\u0302$ gives access to the free energy. From the $\lambda \u2212$ parametric form, it is rather evident that the adiabatic-connection formula constitutes the finite temperature version of the Hellmann–Feynmann theorem,^{284} which is a very useful tool at zero temperature.^{59} In the context of the classical and quantum OCP, the coupling parameter forms of the adiabatic-connection formula are nearly exclusively employed,^{28,283,285}

Consequently, knowledge of the response of a given system to all possible external perturbations, in the limit of a weak perturbation amplitude, gives access to all thermodynamic properties and, thus, also to the free energy.

In practice, Eqs. (34) and (35) constitute the prime motivations for the field of the self-consistent dielectric formalism.^{163,196–201,214–219,286,287} Here, the basic idea is to develop approximate expressions for the LFC $G(q,\omega )$, which give access to $\chi (q,\omega )$, see Eq. (4) that can be employed via Eqs. (31), (32), (34), and (35) to compute a host of dynamic, static, and thermodynamic properties. It is noted that, in general, the LFC is a complicated functional of the DSF. In practice though, the LFC is a complicated functional of the SSF only, which still translates to a complex nonlinear system of integral (in wavenumber and frequency) equations. Recent advances in this field have mainly focused on the accurate description of the strongly coupled electron liquid^{288} with $rs\u227310$, with notable progress reported by Tanaka^{200} and Tolias *et al.*^{201} These schemes are based on the combination of the non-perturbative integral equation theory of classical liquids with the quantum dielectric formalism. The latter scheme^{201} simplifies to the former^{200} in well-defined limits and benefits from the availability of accurate OCP bridge functions extracted from MD simulations.^{289,290} Moreover, we stress that further development of dielectric theories remains an important goal also in the WDM regime, as the best available LFCs are frequency averaged, i.e., they neglect the dependence on *ω*, $G(q,\omega )\u2261G(q)$. Having a reliable and consistently frequency-dependent model for the LFC would allow a number of interesting investigations such as the determination of the effective mass within the concept of Fermi liquid theory,^{291} or the further analysis of a possible experimental detection of the roton feature^{252,292} in the dilute electron gas.

Finally, we note that Eq. (35) has been used by a number of groups^{28,69,94,163,197} to construct representations of the XC free energy of the UEG based on different input datasets for the interaction energy *W*.

As an alternative route to obtain insights into the local field correction without any dielectric approximations, Moroni *et al.*^{253} have suggested to carry out QMC calculations based on the modified Hamiltonian of Eq. (29), but in the static limit of $\omega \u21920$. In the ground state, the stiffness theorem^{59} then gives a straightforward relation between the induced change in the total energy and the static density response function $\chi (q)\u2261\chi (q,0)$. This approach has subsequently been used by different groups to estimate the static density response of a number of systems, including the UEG,^{281,293} a charged Bose gas,^{294} and neutron matter.^{295}

At finite temperature, it is more convenient to directly work with induced densities in reciprocal space,^{296,297}

with $\u27e8\u2026\u27e9q,A$ indicating that the expectation value is to be taken with respect to Eq. (29). The sought-after static linear density response function $\chi (q)$ can then be obtained via the relation

Furthermore, it is straightforward to invert Eq. (4) to obtain the static local field correction $G(q)$. This has allowed Moroni *et al.*^{281} to compute the first accurate results for $G(q)$ of the UEG in the zero-temperature limit, which have subsequently been used as input for the parameterization by Corradini *et al.*^{298}

While being formally exact, the approach that is delineated by Eqs. (37) and (38) requires, for any given combination of $(rs,\theta )$, the execution of multiple independent QMC simulations for each individual **q** and for different perturbation amplitudes *A* (also to ensure nonlinear effects are negligible). An elegant and computationally less demanding alternative is given by the PIMC estimation of the imaginary-time density-density correlation function $F(q,\tau )$ [Eq. (39)], which is discussed in Secs. II C.

### C. Imaginary-time density-density correlation function

The imaginary-time density–density correlation function (ITCF) $F(q,\tau )$ is given by the intermediate scattering function defined in Eq. (2) above but evaluated at the imaginary time $t=\u2212i\u210f\tau $, with $0\u2264\tau \u2264\beta $; in the following, we will simply denote it as

From a theoretical perspective, the estimation of Eq. (39) within a PIMC simulation is straightforward^{190,193} and requires the correlated evaluation of two density operators in reciprocal space at two different imaginary-time slices. For completeness, note that QMC results for the ITCF have been reported for a gamut of systems beyond WDM, including ultracold helium,^{299–303} exotic supersolids,^{304} and the UEG in the zero-temperature limit.^{305}

A first important relation is given by the connection between $F(q,\tau )$ and the static density response function $\chi (q)$, which is known as the imaginary-time version of the fluctuation–dissipation theorem,^{190,293}

In practice, Eq. (40) therefore allows one to estimate the full wave vector dependence of the static linear density response function from a single simulation of the unperturbed system without worrying about nonlinear effects. This approach was extensively used by Dornheim and co-workers over different parameter regimes.^{29,194,288,306,307}

A second important relation is given by the connection between $F(q,\tau )$ and the DSF, which simply comprises a two-sided Laplace transform,

Traditionally, Eq. (41) constitutes the starting point for a so-called analytic continuation,^{308} i.e., a numerical inversion to obtain $S(q,\omega )$ from the QMC data for the ITCF. This is a well-known, though notoriously difficult problem. In fact, an inverse Laplace transform is ill-posed, and the inevitable statistical error bars in the QMC input data lead to different instabilities in practice.^{308} As a consequence, a host of different strategies to deal with the analytic continuation have been suggested in the literature,^{252,299,300,306,309–316} but the quality of the, thus, reconstructed $S(q,\omega )$ generally remains unclear. A notable exception is given by the warm dense UEG, for which a number of exact analytical constraints allow one to sufficiently reduce the space of possible $S(q,\omega )$ to render the numerical inversion of Eq. (41) tractable.^{252,306,317,318} This has allowed Dornheim *et al.*^{252} to present the first accurate results for $S(q,\omega )$, which has given important insights into the nature of the XC induced red-shift in the dispersion of the UEG, cf. Sec. III B 1.

On the other hand, it is well understood that the two-sided Laplace transform constitutes a unique mathematical transformation. Therefore, the ITCF contains by definition the same amount of information as the DSF, albeit in an, at first, unfamiliar representation.^{146,190,192} For example, consider the exact spectral representation of the DSF^{59}

In other words, $S(q,\omega )$ is given by the sum over all possible transitions between the eigenstates *m* and *l* of the full *N*-body Hamiltonian induced by the density operator $n\u0302(q)$, with *P _{m}* being the probability to occupy the initial state

*m*, $\Vert nml(q)\Vert $ is the corresponding transition matrix element, and $\omega lm=(El\u2212Em)/\u210f$ denotes the energy difference. Inserting Eq. (42) into Eq. (41) gives the analogous representation in the

*τ*-domain,

Evidently, the *τ*-decay of $F(q,\tau )$ for a given wave vector **q** is shaped by the characteristic frequencies in the corresponding $S(q,\omega )$. In particular, Eq. (43) directly implies that energetically low-lying excitations such as the roton feature in the dilute UEG^{292} will directly manifest as a less pronounced decay with *τ*. Conversely, it holds that $\omega (q)\u223cq2$ in the non-collective single-particle regime with $q\u226bqF$ (with $qF$ being the usual Fermi wave number^{59}), which will lead to a steeper decay of the ITCF. Dornheim *et al.*^{190} have suggested characterizing this with a relative decay measure of the following form:

and the corresponding results are shown in Fig. 10.

In addition to being of great value for the interpretation of simulation results for $F(q,\tau )$, these considerations open up new avenues for the interpretation of XRTS experiments with WDM.^{146} In particular, the numerical deconvolution of the measured XRTS signal [Eq. (3) above] is generally prevented by the inevitable experimental noise. Therefore, XRTS experiments do not give direct access to the DSF, which contains the relevant physical information about the given system of interest. Instead, one has to take the aforementioned detour over a model description of $S(q,\omega )$, which then has to be inserted into Eq. (3) for a comparison to the experimental observation. Unfortunately, this forward-fitting approach introduces a bias to the interpretation of the experiment that depends on the employed approximation. To our knowledge, no exact simulation or theory is available for the description of $S(q,\omega )$ of real materials in the WDM regime.

In stark contrast, switching to the Laplace domain, i.e., inserting both the XRTS intensity $I(q,\omega )$ and the combined source and instrument function $R(\omega )$ into Eq. (41), makes the deconvolution trivial; this is due to the convolution theorem of $L[\u2026]$, which states that

Indeed, the LHS of Eq. (45) is the ITCF $F(q,\tau )$, which contains the same physical information as the DSF. As we shall see, this allows for the straightforward interpretation of XRTS experiments without any models, simulations, or approximations.

For a uniform system in thermodynamic equilibrium, the DSF fulfills the well-known detailed balance between positive and negative frequencies,^{59,147,319}

The detailed balance condition is a direct consequence of the exact spectral representation of the DSF, see Eq. (42), and is central to the derivation of the FDT, see Eq. (31). In essence, detailed balance is reflected by the FDT, given the odd frequency-parity of $Im\chi (q,\omega )$. Inserting Eq. (46) into Eq. (41) gives the important symmetry relation,^{146,191}

In particular, Eq. (47) implies that $F(q,\tau )$ is symmetric around $\tau =\beta /2$. In practice, one can, thus, directly diagnose the temperature of a given system from the Laplace transform of the XRTS signal, cf. Fig. 16 in Sec. IV; no forward-modeling or simulations are required.

As a final useful property of the ITCF, we mention its relation to the frequency moments of the DSF, which are defined as^{59,153,320}

For example, the first moment (i.e., *α* = 1) is given by the exact f-sum rule,^{59,321}

The same information is encoded into the ITCF via its derivatives with respect to *τ* around *τ* = 0,^{190,322}

In addition to being interesting in their own right, the frequency moments constitute exact constraints that assist in the analytic continuation that is necessary for the reconstruction of the DSF from the QMC ITCF data, see the inversion of the Laplace transform of Eq. (41).^{252,306,317,318} They also constitute the key input for the non-perturbative method of moments that is used by Tkachenko and co-workers to estimate $S(q,\omega )$ based on static structural properties.^{323–326}

### D. Nonlinear density response theory

Even though the idea of LRT is prevalent and very useful in almost all areas of physics, an incomparably richer treasure trove of physics can be obtained by considering stronger deviations from equilibrium, i.e., the nonlinear response. The response of the investigated system to stronger external perturbations not only can be described but also gives direct access to higher-order correlation functions.^{176}

The general definition of higher-order response functions is given by an induced density $\delta n(r\u2192,t)$ expansion,

where $1={r\u21921,t1},\u2009{l}=(1,2,\u2026,l),\u2009d1={dr\u21921,dt1}$. After writing the first three terms explicitly and setting $\chi (1,2)\u2261w(1)(1,2)$, we have^{171,181}

Here, all time integrations run from $\u2212\u221e$ to *t*_{1}, meaning that we need the retarded response functions. Still assuming a relatively weak perturbation, the series can be truncated after the second (quadratic) or third (cubic) term. Thus, the quadratic response function $w(2)$ and the cubic response function $w(3)$ are introduced.

Mathematically, these response functions can be defined as functional derivatives of the one-particle Green's function $g(11\u2032)=\u2212i\u27e8T{\psi +(1\u2032)\psi (1)}\u27e9$ with respect to the external perturbation. They are special cases of more general higher-order correlation functions. For the linear response function, we have

Here, $t1+=t1+\epsilon $. The function *L* is, in general, a true four-point function describing, among other properties, the correlations of two density fluctuations and the scattering of two particles in a medium. The definition of the quadratic response function is

Again, this is a special case of the three-particle correlator $Y(123,1\u20322\u20323\u2032)$. Similar considerations are valid for the cubic response function. Equations (53) and (54) can be used to derive exact equations of motion for the linear and quadratic response functions and also to obtain RPA and other approximations for any of the higher-order response functions.

It may sometimes be useful to introduce effective quantities like the polarization function and dielectric functions to describe in-matter effects. In the case of linear response, it holds

where the polarization function Π and the generalized linear dielectric function *K* were introduced. In such a way, the effective response described by the polarization function is connected to the total response described by *χ*. A similar relation for the quadratic response function reads

This defines the quadratic polarization $\Pi (123)$ and the quadratic dielectric function *K*(123). Naturally, this can be extended to the cubic case.

For simplicity, we restrict ourselves to the quadratic response; the higher-order response has been investigated elsewhere.^{175,195} We also assume the system to be homogeneous and the external perturbation to be of harmonic nature. Note that the latter does not restrict the generality of the considerations. We shall only treat the static *ω* = 0 case herein. The induced density contains various combinations of higher harmonics of the perturbing potential. In wavenumber space, the result is, up to the second order,

Here, $\chi (q)$ denotes the usual static (*ω* = 0) limit of the linear response function (53)^{59,327} and $w(2)(q1,q2)$ denotes the static quadratic response function. A more in-depth derivation up to the third order is given in Ref. 175.

In addition, we note that the quadratic density response function (see Ref. 195 for details) is connected to the imaginary-time structure of the system by the relation

The density operator

is not normalized, and $rl,\tau $ denotes the position of particle *l* at an imaginary time $\tau \u2208[0,\beta ]$. Equation (58), thus, directly implies that all quadratic terms of the nonlinear density response (including mode-coupling effects, see below) can be obtained from a single simulation of the unperturbed UEG.

One can determine all contributions of linear or nonlinear origin to the induced density at any wavenumber or frequency from a higher-order generalization of Eq. (57). For instance, the induced density at the first harmonic, i.e., original perturbing wave vector, is given as a sum of linear and cubic contributions,

where we have introduced $\chi (1,3)(q)=w(3)(\u2212q,q,q)+w(3)(q,\u2212q,q)+w(3)(q,q,\u2212q)$. The signal at the second harmonic is given as a sum of quadratic and quartic contributions,

The diagonal response function is again a special case of the general quadratic response function $\chi (2,2)(q)=w(2)(q,q)$ and can, therefore, be estimated from PIMC results for the quadratic ITCF defined by

cf. Equation (58). At the third harmonic, the signal is given as a sum of cubic and quintic contributions,

The cubic response again determines the nonlinear signal, but the cubic response at the third harmonic is a different realization of the general third order response function, $\chi (3,3)(q)=w(3)(q,q,q)$, compared to the cubic response at the first harmonic.

The ideal quadratic response function can be derived analogously to the (finite-T) Lindhard formula known from LRT. Two terms are contributing. At the second harmonic, a recursion relation is known expressing the ideal quadratic response in terms of the ideal linear response evaluated at the first and second harmonics. This result is due to Mikhailov,^{173,174}

The exact equation of motion for the quadratic response allows for an RPA-like approximation with two linear response functions in the denominator,

There is a square of the Lindhard dielectric function at the first harmonic and a factor comprised of the Lindhard dielectric function at the second harmonic.

The general structure of the equations of motion for the linear and quadratic responses is very similar. The idea of introducing LFCs into the RPA formula for the quadratic response is, thus, obvious

## III. SIMULATION RESULTS

Throughout this work, we limit ourselves to the discussion of simulation results for spin-unpolarized systems with an equal number of spin-up and spin-down electrons, $N\u2191=N\u2193=N/2$. We note that spin effects are expected to play a particularly important role for WDM in magnetic fields and have been studied extensively throughout the literature, e.g., Refs. 28, 69, 73, 94, 199, 218, 244, 328, and 329.

In addition, we note that both QMC and DFT simulations are restricted to a finite number of particles *N* in a finite simulation cell. The corresponding finite-size effects have been analyzed in detail both for ambient conditions^{330–334} and in the WDM regime^{64,317,335} and are not covered in the present work.

### A. Static linear density response

#### 1. Uniform electron gas

Let us begin the discussion of simulation results with an analysis of the static linear density response of the warm dense UEG.^{28,95} In particular, the UEG constitutes one of the most fundamental model systems in physics, quantum chemistry, and related disciplines. Often considered as the archetypal system of interacting electrons, the availability of accurate QMC simulations results^{336–339} in the zero-temperature limit that were subsequently used as input for various parametrizations^{340–343} has been of pivotal importance for the success of DFT at ambient conditions.^{112} The first ground-state calculations for the static linear density response of the UEG have been presented by Moroni *et al.*^{253,281} They were parametrized by Corradini *et al.*^{298} and have been confirmed in the recent study by Chen and Haule.^{344}

At finite temperature, a host of approximate results for the static linear density response has been reported based on, e.g., dielectric theories^{163,196,197,199–201,215,216,218} and classical mappings.^{212,345,346} To our knowledge, the first accurate QMC results for the linear response of the UEG at finite temperature have been presented in Refs. 28, 296, and 297, using the permutation blocking PIMC (PB-PIMC)^{347–349} and the configuration PIMC (CPIMC)^{61,68,70,350} methods. More specifically, these results have been based on a direct simulation of a harmonically perturbed system, see Eqs. (29), (37), and (38). While formally being exact, these efforts have required to perform a considerable number of independent QMC simulations to acquire the density response for a single combination of **q** with $(rs,\theta )$. Therefore, they have been limited to a few density-temperature combinations.

A more efficient strategy to investigate the static linear density response of any given system is offered by the imaginary-time version of the fluctuation-dissipation theorem [Eq. (40)], see the discussion in Sec. II C. In Fig. 4, we show a corresponding investigation of the warm dense UEG at the electronic Fermi temperature *θ* = 1 and the metallic density of *r _{s}* = 4, which is close to sodium.

^{351,352}Panel (a) shows the raw PIMC results for the ITCF $F(q,\tau )$ [Eq. (39)] in the

*τ*-

*q*-plane. For completeness, we note that the

*q*-grid is defined by the system size,

^{64,330,332,335,353}which limits us to discrete values with $q\u22652\pi /L$. The

*τ*-grid, on the other hand, directly follows from the number of high-temperature factors

*P*in our PIMC simulations (see Sec. II A 1), and can, in principle, be made arbitrary fine. In addition, we note that it is fully sufficient to restrict ourselves to the half-interval of $0\u2264\tau \u2264\beta /2$ due to the symmetry relation derived in Eq. (47). The $F(q,\tau )$ symmetry can also be discerned in panel (b), where we show the full

*τ*-dependence of the ITCF for a particular wave vector. Returning to panel (a), it becomes evident that the ITCF approaches the SSF $S(q)$ in the limit of $\tau \u21920,\u2009F(q,0)=S(q)$; this is indicated by the bold-dashed green curve.

While the SSF does not exhibit any maxima due to the absence of spatial order at these conditions, its thermal analogue given by $F(q,\beta /2)$ (light blue dashed curve) exhibits a distinct structure with a pronounced peak around $q\u223c2qF$. We note that the physical behavior of $F(q,\tau )$ has often been described as featureless in the literature^{221,299} and has been treated as a means for the estimation of other properties such as the DSF. Very recently, we have shown that $F(q,\tau )$ directly gives a number of physical insights into phenomena such as the roton feature of the UEG and the related XC induced red-shift of $S(q,\omega )$ without any analytic continuation,^{190,192} see also the discussion of Fig. 10. For example, the single-particle dispersion $\omega (q)\u223cq2$ directly manifests in $F(q,\tau )$ as the steep exponential decay for large *q*; in this regime, the ITCF can, indeed, be modeled very accurately based on a simple, semi-analytical single-particle model for the imaginary-time diffusion process, see Ref. 192 for details. Moreover, the well-known sum rules of the DSF manifest in the ITCF as derivatives with respect to *τ* around *τ* = 0 [Eq. (50)].

This is demonstrated in Fig. 4(b), where we illustrate $F(q,\tau )$ for $q=1.88qF$. More specifically, the solid green curve shows our exact PIMC results, and the dash-double-dotted yellow line has been obtained from the exact f-sum rule [Eq. (49)] that describes the slope of $F(q,\tau )$ at the origin; we find perfect agreement between the two curves, as it is expected. The availability of exact PIMC benchmark data for the case of the UEG also allows us to unambiguously assess the accuracy of different approximations. The dash-dotted gray curve corresponds to the ubiquitous RPA, which provides the least accurate description of all curves shown in Fig. 4(b). Evidently, the mean-field description of the dynamic density response function [i.e., setting $G(q,\omega )\u22610$ in Eq. (4)] on which RPA is based is not appropriate in this regime, where electronic XC effects play an important role. For completeness, we note that the RPA still exactly fulfills the f-sum rule, which only describes the derivative at *τ* = 0 but is agnostic with respect to the particular value of $F(q,0)$.

The dashed blue curve has been computed within the static approximation, i.e., by setting $G(q,\omega )\u2261G(q,0)$ in Eq. (4). In particular, Dornheim *et al.*^{252} have reported that the static approximation gives highly accurate results for the DSF $S(q,\omega )$ of the UEG for metallic densities $rs\u22724$. Evidently, the same trend manifests in the ITCF, which is in very good agreement with the exact PIMC results over the entire *τ*-range. For completeness, we note that the static approximation induces a slight overestimation of $F(q,0)=S(q)$ for $q\u22732qF$.^{159,354} While almost being negligible for an individual *q*, this systematic error accumulates within integrated properties such as the interaction energy *W* [Eq. (34)]. This problem can be solved by combining the static local field correction of the UEG with a consistent limit for large wave numbers $q=|q|\u226bqF$ that is connected to the so-called on-top pair correlation function, $g(r=0)$. The corresponding effective static approximation^{159} has been analytically parametrized in Ref. 354.

Finally, the dotted red curve in Fig. 4(b) has been obtained by combining the SSF with the semi-analytical single-particle imaginary-time diffusion model presented in the recent Ref. 192. Remarkably, this simple model, too, fulfills the exact f-sum rule and accurately captures the *τ*-dependence of the ITCF over the entire *τ*-range.

Let us next get to the task at hand, which is the estimation of the static linear density response function $\chi (q)$ from the PIMC data for $F(q,\tau )$ via Eq. (40). The results are shown in Fig. 4(c), where we show the, thus, obtained PIMC results for $\chi (q)$ as the green crosses. In the limit of small *q*, the exact density response is given by Kugler^{320} (dotted yellow curve),

this is a direct consequence of the perfect screening in the UEG and is not reproduced by the density response of the ideal Fermi gas (solid blue curve).

The density response of the UEG to an external static perturbation attains a maximum around $q\u223c2qF$. From a physical perspective, this can be understood intuitively from the following considerations: for intermediate *q*, the wavelength $\lambda =2\pi /q$ of the perturbation is comparable to the average interparticle distance *d*. Consequently, such a perturbation induces a spatial pattern that automatically minimizes the interaction energy between the electrons, which maximizes the reaction of the system.^{355} A more technical explanation can be gleamed from the ITCF depicted in panel (a). In particular, Eq. (40) states that $\chi (q)$ is proportional to the area under $F(q,\tau )$ for a given value of **q**. Therefore, the peak in $\chi (q)$ is directly related to the maximum in $F(q,\beta /2)$. In other words, the existence of a hypothetical spatial pattern that minimizes the interaction energy—but which is not actually manifested in the unperturbed UEG as there is no maximum in $S(q)$ at these conditions—manifests itself as a reduced decay along *τ* for intermediate *q*; this is directly translated to the static density response function shown in Fig. 4.

The dash-dotted black curve shows the static linear density response function within the RPA, which systematically underestimates the true response of the UEG for $q\u2273qF$. This, too, can be attributed to the behavior of the ITCF shown in Fig. 4(b), which is too small over the entire *τ*-range. At the same time, the RPA reproduces the correct perfect screening described by Eq. (67).

Finally, we have carried out extensive new PIMC simulations for the harmonically perturbed electron gas; the results for the dependence of the induced density [Eq. (37)] on the perturbation amplitude *A* are shown in Fig. 4(d) for integer multiples of the minimum wave number $qmin=2\pi /L$. Specifically, the data points show our raw PIMC simulation data (divided by *A*, such that the response attains a constant value in the LRT limit), and the corresponding dotted curves have been obtained from cubic fits according to Eq. (60). Evidently, the functional form nicely reproduces the PIMC results everywhere, as it is expected.^{167,175} The horizontal bars show the LRT limit computed from the ITCF, which is consistent both to the fits and to the PIMC data on which they are based. The linear coefficient in Eq. (60) directly corresponds to the linear response function $\chi (q)$, which are included in Fig. 4(c) as the red dots. For completeness, we note that the cubic coefficients correspond to the cubic response at the first harmonic.^{167,175} Evidently, these independent data for $\chi (q)$ that are based on the actual response of the UEG are in excellent agreement with the ITCF-based evaluation of Eq. (40) over the entire *q*-range. This demonstrates the high consistency of both our PIMC simulations as well as of the underlying theoretical framework.

Based on the described estimation of $\chi (q)$ from the ITCF, Dornheim *et al.*^{194} have carried out extensive PIMC calculations for $\u223c50$ density-temperature combinations; the corresponding parameters are depicted as the red crosses in the *r _{s}*-

*θ*-plane shown in Fig. 5(a). In addition, the black squares and dashed blue line show the ground-state QMC calculations by Moroni

*et al.*

^{281}(MCS) and the corresponding parameterization by Corradini

*et al.*

^{298}(CDOP). The task at hand was to obtain an accurate representation of the static local field correction $G(q)$, covering the entire relevant parameter range, i.e., the shaded green area in Fig. 5(a).

On the one hand, the inversion of Eq. (4) to compute $G(q)$ from PIMC data for $\chi (q)$ is straightforward. On the other hand, $G(q;rs,\theta )$ exhibits a nontrivial dependence on **q** based on different combination of *r _{s}* and

*θ*. This is illustrated in Fig. 5(b), where we show the static local field correction in the

*r*-

_{s}*q*-plane at the electronic Fermi temperature,

*θ*= 1. Specifically, the red crosses show the

*q*-dependence of our PIMC results for $G(q)$, and we observe the following trend. For low densities (i.e., large

*r*), $G(q)$ exhibits a pronounced increase in the limit of large

_{s}*q*. This increase becomes less pronounced with decreasing

*r*and actually becomes negative for high densities. From a physical perspective, this interesting behavior can be traced back to the XC contribution to the kinetic energy,

_{s}^{356,357}which is negative for some parameters at finite temperatures.

^{306,358}To construct a reliable representation of $G(q;rs,\theta )$ that is capable to capture these interesting trends, Dornheim

*et al.*

^{194}have trained a deep neural network to learn the appropriate functional dependence. The results are shown as the green surface in Fig. 5. Evidently, the neural net nicely reproduces the PIMC data where they are available and smoothly interpolates in between. The high quality of this surrogate model

^{359}was subsequently confirmed by the validation against independent data that had not been included into the training procedure, and by the very recent study of electronic exchange-correlation effects by Hou

*et al.*

^{357}

Since its publication in 2019, the neural-net representation of $G(q;rs,\theta )$ has been used for a number of applications, including the modeling of XRTS experiments,^{159} the estimation of ionization-potential depression,^{160} and as input for a theory of nonlinear effects in the UEG,^{175} which is covered in more detail in Sec. III C.

#### 2. Warm dense hydrogen

Let us next investigate the static linear density response of a real system. In Fig. 6, we show the static density response function for hydrogen at the electronic Fermi temperature, i.e., *θ* = 1. Panel (a) corresponds to the metallic density of *r _{s}* = 2, and the green circles have been obtained by Böhme

*et al.*

^{77}by exactly solving the electronic problem in a fixed ion potential with PIMC. In particular, they have applied an external harmonic potential and subsequently measured the response in the electronic density; see also Ref. 78 for a more detailed description of the simulation setup. The dotted red line shows results for the density response of the UEG at the same conditions and has been computed based on the neural-net representation of $G(q;rs,\theta )$.

^{194}Evidently, the electronic density response of hydrogen at these conditions strongly resembles the behavior of a free electron gas, as most of the electrons are delocalized here.

^{360}

Let us for now ignore the black stars and proceed to Fig. 6(b), where we present the same analysis for a lower density at *r _{s}* = 4. Such conditions are expected to exhibit a number of interesting physical effects and can be realized experimentally, for example, in hydrogen jets.

^{361}Low-density systems constitute an interesting benchmark for the assessment of electronic XC effects,

^{263,362}which are more pronounced at larger values of the quantum coupling parameter

*r*.

_{s}^{58}In fact, they constitute a realistic option for an experimental observation of the roton feature in the dispersion of the warm dense electron gas that has been reported and analyzed in earlier works.

^{190,252,292}In addition, these conditions are characterized by a partial localization of the electrons around the ions,

^{360}which means that widely used concepts such as the decomposition into bound and free electrons

^{148,186}are expected to break down. Indeed, the analysis by Böhme

*et al.*

^{77}has revealed that the actual electronic density response of hydrogen (green) is substantially lowered compared to the UEG model at the same conditions (dotted red). The dash-dotted yellow curve in Fig. 6(b) has also been computed based on the UEG, but with modified parameters $rs*$ and $\theta *$ that correspond to an effective free-electronic fraction of $\alpha =0.6$, i.e., 60%; this is consistent with the value of $\alpha =0.54$ reported by Militzer and Ceperley.

^{360}While this effective density response overall exhibits the correct magnitude around $q\u223c2qF$, it deviates from the PIMC reference data in particular, for large

*q*. In fact, the true electronic density response of hydrogen even exceeds the red UEG curve for $q\u22734qF$, which has been attributed to an isotropy breaking due to the presence of the ions at small length scales $\lambda q=2\pi /q$. This is a strong indication that ionization models do not universally work over all relevant length scales.

While being important in their own right, the availability of highly accurate PIMC benchmark results also allows one to benchmark a number of previously used approximations. Of particular interest is the assessment of the accuracy of thermal DFT, which constitutes the workhorse of WDM theory.^{11,89,363} In Fig. 6, we include DFT results that have been obtained by Moldabekov *et al.*^{255} based on the PBE XC functional^{364} for the electronic density response of the same ion snapshot as the black stars. As mentioned above, hydrogen basically behaves as a free electron gas at *r _{s}* = 2, and the DFT results are in good agreement with the PIMC reference data. For

*r*= 4, the agreement is less good, and the DFT calculations overestimate the actual response over the entire depicted

_{s}*q*-range. This can be traced back to an underestimation of the localization of the electrons around the ions due to the PBE XC functional, which is known to be afflicted with self-interaction errors.

^{365}At the same time, we note that the overall agreement between PIMC and PBE-based DFT calculations is better compared to the effective UEG model with $\alpha =0.6$. In particular, the DFT simulations are capable of accurately reproducing the effect of the isotropy breaking at large

*q*; this has important implications, as we shall discuss in the following.

Having found $\chi (q)$, e.g., from a DFT calculation, it is straightforward to invert Eq. (17) to obtain the static XC kernel $Kxc(q)$ that is fully consistent with any KS-reference function $\chi S(q,\omega )$. This procedure is outlined in more detail in Fig. 17 and allows one to generate the material-specific XC kernel for any XC functional; no evaluation of a functional derivative is required. Indeed, the latter point is notoriously difficult in practice and had hitherto to our knowledge prevented the computation of XC kernels of extended systems beyond adiabatic LDA (ALDA) or adiabatic generalized-gradient approximations (AGGA) for extended systems.^{366}

In Fig. 7, we show the corresponding results for the XC kernel of hydrogen for *r _{s}* = 4 and

*θ*= 1. As a reference function to obtain the kernel from the physical density response function $\chi (q)$, Moldabekov

*et al.*

^{255}have used the same mean-field response function $\chi 0(q)$ that has been obtained without any XC functional for all datasets. Specifically, the black squares have been computed from the PIMC data for $\chi (q)$ by Böhme

*et al.*

^{77}and constitute a highly accurate benchmark for the other results. The dashed red curve corresponds to the popular ALDA kernel, which is computed from the parabolic compressibility sum rule expansion around

*q*= 0 for a UEG,

at the same density and temperature. Clearly, it does not capture the true behavior of the kernel and substantially overestimates it over the entire depicted *q*-range. The solid dashed line has been obtained from the accurate neural-net representation of $G(q;rs,\theta )$ of the UEG;^{194} it, by definition, reproduces the ALDA expansion for small *q* and does not reproduce the PIMC results for larger *q* either. In fact, using the RPA expression for $\chi (q)$ in Eq. (17), which corresponds to setting $Kxc(q)\u22610$, constitutes a superior approximation compared to either ALDA or the full UEG model, which lead to an actual deterioration of the attained accuracy. This has important implications for the LR-TDDFT simulation of WDM, as we discuss in more detail in Sec. III B.

In stark contrast to the ALDA model, the new DFT-based results for $Kxc(q)$ from Ref. 255 that have been computed within the LDA (blue circles, using the functional by Perdew and Zunger^{341}) PBE^{364} (orange triangles) and SCAN^{367} (green diamonds) qualitatively capture the *q*-dependence of the PIMC reference data over the entire depicted *q*-range. In particular, SCAN exhibits the best accuracy for intermediate wave numbers with $q\u223c2qF$, although there remain significant differences to the exact PIMC reference data.

To summarize, the recent approach for the computation of $Kxc(q)$ constitutes a promising route to study electron–electron correlations of real materials within the framework of KS-DFT. In combination with accurate reference results such as the PIMC data for hydrogen by Böhme *et al.*,^{77} it provides an ideally suited tool for the assessment of the accuracy of different XC functionals. Future efforts in this direction will include the rigorous benchmark of orbital-dependent hybrid functionals, which have already been successfully applied to the study of WDM,^{145,368} and the development and benchmark of new functionals that are explicitly designed to meet the challenges of capturing the electronic density response of matter at extreme conditions.

### B. Dynamic linear density response

In Sec. III A, we have given an overview of some recent promising developments regarding the description of the electronic density response in the static limit of $\omega \u21920$. While being an important step in the right direction, many important applications such as the modeling of XRTS experiments^{147,164,166,189} or the construction of advanced XC functionals within the adiabatic-connection fluctuation–dissipation formulation of DFT^{99,156} require, as input, some information about the full frequency-dependence of the electronic density response of a given system.

Unfortunately, the accurate estimation of the dynamic electronic properties of WDM is even more difficult than in the static case, and only a handful of methods are serious contenders. Among these, the nonequilibrium Green's functions (NEGF) method^{276,278} deserves a special place, as it is, as the name suggests, capable of treating real nonequilibrium conditions as they occur, for example, during the stopping of a projectile in a medium.^{142,369–371} While some remarkable methodological improvements have been reported over the last few years,^{372,373} it remains unclear if NEGF calculations of highly excited states as they occur in WDM are computationally feasible, and if common approximations to the self-energy are sufficient to capture the impact of electron–electron interactions. To our knowledge, no NEGF calculations of real WDM systems have been presented in the literature.

A second route toward dynamic electronic properties of WDM is RT-TDDFT.^{162,259,374–376} RT-TDDFT is a computationally efficient approach that allows for the use of large supercells and basis sets along with an accurate description of the electron–ion interaction. Among its key capabilities for WDM purposes are its ability to model XRTS experiments without the ubiquitous Chihara decomposition^{189} and the possibility to estimate both nonequilibrium^{377} and nonlinear effects.^{267} For example, RT-TDDFT allows for a description of nonadiabatic electron–ion dynamics such as those that are relevant to electronic stopping power,^{109,142,378–380} which is critical to self-heating in inertial fusion applications.^{381,382} On the downside, systematically improvable approximations to the true dynamic XC potential^{383,384}—which replaces the standard XC functional of equilibrium DFT—have been difficult to develop, and most calculations are done with adiabatic approximations to standard functionals. Nevertheless, RT-TDDFT has performed well when compared to experiments, and stopping power calculations have so far demonstrated reasonable agreement with WDM experiments^{106,385} and databases of experimental results in less extreme conditions.^{386–389}

Currently, the arguably most widely used *ab initio* method for the description of electron dynamics in WDM is given by linear response TDDFT (LR-TDDFT).^{161,162,164,166} To be more specific, LR-TDDFT is based on a standard equilibrium KS-DFT simulation of a given system, which gives access to a set of KS orbitals ${\varphi \alpha}$. The latter are then used to compute the dynamic KS-response function $\chi S(q,\omega )$ via Eq. (16); the key ingredient is then given by the XC kernel $Kxc(q,\omega )$, which relates $\chi S(q,\omega )$ to the actual physical dynamic linear density response function $\chi (q,\omega )$ of the system of interest via Eq. (17). Formally, LR-TDDFT and RT-TDDFT should give equal results in the limit of weak perturbations; a practical investigation of this assumption for the density response of WDM constitutes an important topic for future research. The main bottleneck of LR-TDDFT is given by its dependence on a particular XC kernel, which has to be material-specific and must be consistent with the XC functional that has been used to determine $\chi S(q,\omega )$. In practice, both conditions are generally violated by the most commonly used model kernels,^{366,390} and the impact of such inconsistencies on the calculated observables remains as of yet poorly understood. Here, we stress that this is not a shortcoming of the general formulation of the LR-TDFT, but of the inconsistency in the *ad hoc* approximations used for the XC kernel. As we have elaborated in Sec. III A 2, Moldabekov *et al.*^{255} have recently introduced a formally exact framework to estimate the consistent static kernel for any given XC functional and material, and its utilization for LR-TDDFT calculations within the *static approximation* where $Kxc(q,\omega )\u2261Kxc(q)$ in Eq. (17) is discussed in Sec. III B 2.

As the final method for the computation of the dynamic density response, we mention PIMC. While, from a conceptual perspective, the real-time propagation of any observable in Feynman's path-integral picture of quantum mechanics^{228} is straightforward, it leads to an oscillating complex exponential function in practice; this is the origin of the infamous dynamic phase problem.^{391–393} Still, PIMC methods give direct access to the imaginary-time dynamics of a given system^{190,193,195} in thermodynamic equilibrium; exact PIMC results for the ITCF $F(q,\tau )$ of the warm dense UEG are shown in Fig. 4. Indeed, the ITCF, in principle, contains the same information as the DSF $S(q,\omega )$—albeit in an unfamiliar representation^{190,192}—and is directly useful for the interpretation of XRTS experiments,^{146} see Sec. IV. Yet, the analytic continuation^{308} from the *τ*-domain to real frequencies is notoriously difficult, and information about specific features might be lost in the process. It is, however, possible for the case of the UEG^{252,306,318} as we shall explain in Secs. III B 1.

#### 1. Uniform electron gas

The accurate estimation of the dynamic density response and dynamic structure factor of the UEG has been an active topic of investigation over many decades.^{161,203,214,249–251,268,287,318,326,394–400} In the ground state, Takada and collaborators^{396,397} have presented accurate data for $S(q,\omega )$ over a broad range of densities, with important implications for scattering experiments with aluminum.^{352,396} Only very recently, LeBlanc *et al.*^{400} have presented accurate data for the UEG density response based on novel diagrammatic Monte Carlo calculations at low temperatures, $\theta =0.1$; it remains to be seen if these promising efforts can be extended to the WDM regime in future works.

The first highly accurate results for the dynamic electronic density response of the warm dense UEG have been presented by Dornheim *et al.*^{252} based on exact PIMC results for the ITCF $F(q,\tau )$. To render the numerical inversion of Eq. (41) tractable, the original problem has been re-cast into the reconstruction of the dynamic local field correction $G(q,\omega )$. More specifically, trial solutions for $G(q,\omega )$ give straightforward access to $\chi (q,\omega )$ [Eq. (4)], and via the fluctuation-dissipation theorem [Eq. (31)] also to the corresponding DSF $S(q,\omega )$; finally, the latter is inserted into Eq. (41) and compared to the PIMC reference data for $F(q,\tau )$. While formally being equivalent to the original problem of the analytic continuation,^{308} this approach allows one to incorporate a number of exact properties of $G(q,\omega )$,^{161,252,287,306} which sufficiently constrains the space of possible trial solutions for $S(q,\omega )$.

In Fig. 8, we show corresponding results for both the real (top panel) and the imaginary part (bottom panel) of $G(q,\omega )$ for the UEG at *θ* = 1, *r _{s}* = 6

^{318}with the shaded gray area indicating the given uncertainty interval. For the real part, the static $\omega \u21920$ limit can be inferred from $\chi (q)$, while the high frequency $\omega \u2192\u221e$ limit follows from the DSF cubic sum rule, which allows expressing $M3S$ [cf. Eq. (48)] in terms of static properties such as the $S(q)$.

^{59,401,402}It is interesting that the real part exhibits a nontrivial and, in fact, non-monotonic behavior between these limits, with a pronounced maximum around $\omega \u22482\omega p$, with $\omega p=3/rs3$ being the usual plasma frequency.

^{59}The corresponding imaginary part of the dynamic local field correction vanishes in both the limits of $\omega \u21920$ and $\omega \u2192\u221e$ and exhibits a single broad peak around $\omega \u22485\omega p$.

From a theoretical perspective, accurate knowledge of $G(q,\omega )$ gives one direct access to the dynamic density response function $\chi (q,\omega )$ [cf. Eq. (4)] and, in this way, a gamut of other properties such as the dynamic dielectric function $\u03f5(q,\omega )$ or the electric conductivity $\sigma (q,\omega )$ at finite wave numbers. This has been explored in detail in recent works by Hamann *et al.*^{318,399}

Due to its central role in the interpretation and modeling of XRTS experiments^{147,149} [cf. Eq. (3)], one of the most interesting properties in this regard is given by $S(q,\omega )$, which has been studied extensively in Refs. 252, 292, 306, and 317. In particular, Dornheim *et al.*^{252} have found an XC induced red-shift in the position of the maximum of the true $S(q,\omega ),\u2009\omega (q)$, compared to the RPA for metallic densities. This is shown in Fig. 9(a) for *r _{s}* = 4 at the electronic Fermi temperature,

*θ*= 1. Specifically, the solid green curve depicts the exact $\omega (q)$ based on the full PIMC results for the dynamic local field correction $G(q,\omega )$, which are down-shifted compared to the RPA mean-field description that is depicted as the dash-dotted black curve. In addition, the dashed blue curve in Fig. 9(a) shows the static approximation,

^{252}where the dynamic local field correction in Eq. (4) is approximated by its exact static limit, i.e., $G(q,\omega )\u2261G(q)$. The results are in excellent agreement with the full reference PIMC data. In other words, the combination of a dynamic description on the level of the RPA [where the frequency-dependence comes from the dynamic Lindhard function $\chi 0(q,\omega )$] with a static (frequency averaged) treatment of electron–electron XC effects via $G(q)$ provides a highly accurate description of the DSF of the UEG at metallic densities. Note that this is consistent with previous observations in the ground-state limit.

^{161,403}

Let us postpone the discussion of the dotted red curve for now and proceed with Fig. 9(b), where we show the same analysis for a lower density with *r _{s}* = 10. Owing to the role of the Wigner–Seitz radius as the quantum coupling parameter of the UEG,

^{58,59}these conditions are located at the margin of the strongly coupled electron liquid regime where electron–electron correlation effects play the dominant role.

^{201,288}In this case, the exact PIMC reference data for $\omega (q)$ exhibit a nonmonotonous behavior, with a minimum at intermediate wave numbers, $q\u223c2qF$. In fact, this phenomenologically resembles the roton feature that is well-known from quantum liquids such as ultracold helium.

^{301,303,404–408}It is worth noting that a phonon-roton spectrum is an unambiguous feature of collective excitations in classical liquids and supercritical fluids.

^{408,409}This includes not only model liquids (hard-sphere,

^{410}Lennard-Jones,

^{411}Yukawa,

^{412}and inverse power law

^{413}) but also real liquids (water,

^{414,415}elemental noble liquids,

^{416}and liquid alkali metals

^{417}), for which the phonon-roton structure has been predicted by

*ab initio*molecular dynamics simulations and observed by inelastic x-ray and neutron scattering. Nevertheless, despite some numerical hints,

^{418}a plasmon-roton spectrum has not been univocally observed in MD simulation studies of the classical OCP,

^{419,420}which is characterized by the onset of a negative long-wavelength limit dispersion, $d\omega (q)/dq<0$ as $q\u21920$, at intermediate coupling $\Gamma \u223c10$.

^{421–425}It is interesting that both quantum and classical cases have been elucidated by the onset of spatial structure:

^{409}through the form of the Feynman ansatz

^{426}[$S(q)$ connection] for quantum liquids and through the quasi-localized charge approximation

^{427}[$g(r)$ connection] or the static local field corrected dielectric formalism

^{428}[$S(q)$ connection] for classical liquids. Yet, no signatures of spatial structure can be found in either the static structure factor $S(q)$ or the pair correlation function $g(r)$ for the UEG at the conditions of Fig. 9(b).

To explain the mysterious roton feature in the UEG, Dornheim *et al.*^{292} have recently introduced the pair alignment model (PA) that interprets the minimum in $\omega (q)$ in terms of the relative length scales upon applying an external harmonic perturbation. In particular, the fluctuation-dissipation theorem [Eq. (31)] directly implies that the behavior of $S(q,\omega )$ can be fully explained by understanding the response of any given system to a monochromatic external perturbation of the same wave vector **q** and frequency *ω*. In the vicinity of the roton minimum in $\omega (q)$ around $q\u223c2qF$, the wavelength of the cosinuoidal perturbation is comparable to the average interparticle distance, $\lambda q=2\pi /q\u223cd$. Therefore, the alignment of the electrons to the respective potential minima is associated with a reduction in the interaction energy *W*.

The PA model introduced in Ref. 292 postulates that the XC induced down-shift of the true $\omega (q)$ compared to the RPA is due to a deficiency of the latter to capture the true change in the interaction energy *W*. The latter can be accurately quantified using a suitable effective potential,^{140,429,430} and the difference between the true change in the interaction and the respective RPA result is indicated by the vertical red arrows labeled as $\Delta Wxc$ in Fig. 9; the corresponding dotted red curves have been computed as $\omega PA(q)=\omega RPA(q)+\Delta Wxc(q)$. For the metallic density of *r _{s}* = 4, the PA model nicely captures the true down-shift of $\omega (q)$ and, therefore, is directly relevant for the modeling of XRTS experiments of WDM.

^{157}For the strongly coupled case of

*r*= 10, the PA model follows the dashed blue static approximation curve rather than the actual PIMC results and, consequently, underestimates the true depth of the roton minimum at these conditions. This is not a coincidence, and it can be attributed to the construction of the PA model as an average energy shift. Empirically, the static approximation exhibits the same behavior and merges the sharp roton peak at low

_{s}*ω*with the shoulder at $\omega RPA(q)$, see Ref. 292 for a more detailed explanation.

We note that the pair alignment that causes the XC induced down-shift of $\omega (q)$ is the same mechanism that has been evoked to explain the peak in the static linear density response function $\chi (q)$ in Fig. 4. In fact, $S(q,\omega )$ and $\chi (q)$ are related via the inverse-frequency sum rule, which states that

Here, $M\u22121S$ denotes the inverse-frequency moment of the DSF defined in Eq. (48). It is evident from Eq. (69) that a shift of spectral weight in $S(q,\omega )$ to smaller frequencies *ω* is associated with an increase in the magnitude of $\chi (q)$. This explains the emergence of an increasingly sharp peak in $\chi (q)$ around $q=2qF$, which has been reported in previous studies of the electron liquid both at finite temperature^{201,288} and in the ground state.^{431}

While having exact reference results for $S(q,\omega )$ gives one valuable insights into a number of physical effects, the required analytic continuation^{308} is not expected to be feasible for arbitrary complex systems due to the ill-posed nature of the inverse Laplace transform. On the other hand, it is clear that, from a mathematical perspective, the ITCF $F(q,\tau )$ contains exactly the same information as the DSF, although in an unfamiliar representation. The task at hand is, thus, to understand the manifestation of a given physical effect of interest in the imaginary-time domain. For example, it is obvious from Eq. (47) that $F(q,\tau )$ is symmetric around $\tau =\beta /2$ and constitutes a direct and highly sensitive measure for the temperature of any system in thermodynamic equilibrium.^{146,190} This has important implications for the interpretation of XRTS experiments and is discussed in Sec. IV.

In addition, understanding the physics encoded into $F(q,\tau )$^{192} is also directly helpful for the interpretation of simulation results. Very recently, Dornheim *et al.*^{190} have demonstrated that the roton feature in the dispersion $\omega (q)$ of the DSF of the UEG manifests as a reduced decay with respect to *τ* of the ITCF for a given wave vector **q**, see the spectral representation of $F(q,\tau )$ in Eq. (43). In other words, the existence of energetically low-lying density excitations manifests as the stability of electron–electron correlations along the diffusion through the imaginary time *τ*; this can be quantified using the relative decay measure $\Delta F\beta /2(q,\tau )$ defined in Eq. (44).

The results are shown in Fig. 10 for *θ* = 1 with green and red datasets corresponding to *r _{s}* = 4 and

*r*= 10, respectively. More specifically, the symbols show the PIMC reference data for $\omega (q)$ that are also depicted in Fig. 9, and the corresponding solid lines show the relative decay measure $\Delta F\beta /2(q,\tau )$ that has been directly computed from the ITCF without the need for an analytic continuation. We note that all curves have been re-normalized with respect to their respective plasmon limit for $q\u21920$, which is given by the usual plasma frequency $\omega p$ in the case of $\omega (q)$. In addition, all curves exhibit a characteristic single-particle limit for $q\u226bqF$. For the DSF, it is given by the well-known single-particle dispersion $\omega (q)\u223cq2$. For the ITCF, one finds $\Delta F\beta /2(q\u2192\u221e,\tau )=1$, which means that the corresponding re-normalized curves attain a constant value in the short wavelength limit. This is due to the increasingly steep decay of $F(q,\tau )$ with $0\u2264\tau \u2264\beta /2$ in this regime, which can be explained with a simple semi-analytical Gaussian imaginary-time diffusion model; see Ref. 192 for an extensive discussion of the physical origin of the

_{s}*τ*-dependence of the ITCF.

For intermediate *q*, we find a close correspondence between the respective $\omega (q)$ and $\Delta F\beta /2(q,\tau )$ for the same value of the density parameter *r _{s}*. Indeed, the relative decay measure even becomes slightly negative for

*r*= 10. We note that the phenomenological similarity between the roton minimum in the DSF and the reduced

_{s}*τ*-decay of the ITCF becomes even more apparent at

*r*= 20, where the alignment of pairs of electrons

_{s}^{292}has a more pronounced impact on the physical observables of the system.

Despite these works, we emphasize that the thorough development of a dynamic quantum many-body theory in the imaginary-time domain as of yet remains in its infancy. Future works may start with the investigation of other systems such as quantum liquids^{301,303,404} and, of course, real WDM applications that include both electrons and ions. It is evident from Eq. (41) that every theoretical model such as not only the Chihara decomposition^{148,151,186} but also TDDFT simulations can be straightforwardly translated into a description of the ITCF. At the same time, we also stress that working in the *τ*-domain might allow for the development of entirely new concepts that naturally emerge from Feynman's imaginary-time path-integral representation of statistical mechanics but may not have an obvious analogue in the traditional frequency-domain.

#### 2. Warm dense matter

Let us next extend our analysis of the dynamic electronic density response of WDM beyond the UEG model. As we have elaborated above, no exact method for the simulation of frequency-dependent properties of real WDM systems is currently available to our knowledge. Being motivated by the remarkable performance of the static approximation $Kxc(q,\omega )\u2261Kxc(q)$ for the UEG for metallic densities $rs\u22724$, we, here, explore the combination of LR-TDDFT (cf. Sec. II A 2) with a consistent static XC kernel for the example of warm dense hydrogen. The results are shown in Fig. 11 for (a) *r _{s}* = 2 and (b)

*r*= 4 at the electronic Fermi temperature,

_{s}*θ*= 1. More specifically, all curves have been computed via Eq. (17) using the KS-response function $\chi S(q,\omega )$ based on a set of KS orbitals from an equilibrium DFT simulation with the LDA functional.

Let us start our analysis for *r _{s}* = 2, where hydrogen is mostly ionized, as it has been explained in Sec. III A 2. In this case, the depicted LR-TDDFT results closely resemble the DSF of the UEG.

^{252}In particular, the solid black curve has been computed by setting $Kxc(q,\omega )\u22610$ in Eq. (17), which is commonly being referred to as RPA in the DFT literature.

^{162}Including the appropriate PIMC results for the static XC kernel $Kxc(q)$ by Böhme

*et al.*

^{77}leads to the dotted red curve; it exhibits a similar XC induced red-shift as we have observed in the case of the UEG in Sec. II B 1. Finally, the dashed green curve has been obtained by replacing the true static kernel with the ALDA, which is based on the $q\u21920$ limit of the UEG at the same parameters, cf. Eq. (68). Evidently, the ALDA is in good agreement with the dotted red curve. This is expected, as the electrons of hydrogen basically behave like a UEG at these conditions due to the high degree of ionization.

^{77,360}

In stark contrast, the electrons are strongly influenced by the ions at *r _{s}* = 4, which is shown in Fig. 11(b). In particular, we find a diffusive peak around

*ω*= 0, which is a direct consequence of the partial localization of the electrons around the protons. For these parameters, the RPA and the static approximation are in close agreement, as the exact PIMC results for the static XC kernel nearly vanish, cf. Fig. 7. The corresponding ALDA kernel fails to capture this trend, and the dashed green line substantially deviates from the dotted red reference curve. In other words, using the ALDA kernel that is based on the UEG in combination with the dynamic KS-response function $\chi S(q,\omega )$ that has been computed from the KS orbitals for a specific model system leads to an actual deterioration of the quality of the results compared to the RPA. Therefore, we view the construction or utilization of universal model kernels as more or less futile, since a good kernel must (1) be consistent with the employed XC functional that has been used to compute $\chi S(q,\omega )$ and (2) depend on the physical behavior of the given system of interest.

A promising route is given by the recent framework by Moldabekov *et al.*,^{255} which allows one to compute the appropriate static XC kernel for any given system, and for arbitrary XC functionals across Jacob's Ladder^{432} without the need for an explicit evaluation of functional derivatives. The accuracy of a corresponding LR-TDDFT calculation will then hinge on the static approximation, which has to be studied in more detail in future works. In this regard, we mention the advent of direct PIMC simulations of real WDM systems without the fixed-node approximation,^{77,78} which will give us the first exact results for the ITCF $F(q,\tau )$ of hydrogen in the near future. In addition to being interesting in their own right, such results will constitute an unassailable benchmark for the development of TDDFT simulations of WDM and will guide the development of improved approximations that can be applied to more complex systems that are beyond the scope of PIMC.

Based on the methodology introduced in Sec. II A 3, RT-TDDFT is a powerful method for evaluating experimentally relevant features such as the DSF, and the electrical conductivity (in the optical limit $q\u21920$). Figure 12 shows the frequency-dependent electrical conductivity^{23} of iron at earth-core conditions ($P\u223c320$ GPa or $rs\u223c2.2,\u2009T\u223c0.55$ eV). The red curve shows the tensor averaged result of several ionic configurations of size N = 16 with an energy resolution $\Delta \omega \u223c$ 0.11 eV, which is proportional to the inverse of the total propagation time ($t\u223c1500$ a.u.). We compare our calculations with prior results (black dashed and dotted curves) obtained using the KG formula based on static DFT (up to 20 eV).^{433,434} In general, KG results are affected by finite-size effects and are highly dependent on the density and location of the KS eigenvalues.^{435,436} The red dotted line shows the Drude fit to RT-TDDFT data to obtain the DC conductivity, which is in the range observed in previous *ab initio*^{433,434,437} methods.

Next, in Fig. 13, we consider a RT-TDDFT calculation of the DSF of isochorically heated aluminum ($Te=0.3$ eV) at a momentum transfer of $q=0.67$ a.u., compared to both experiment and a DFT-based theory for these same conditions first studied in Ref. 145. The RT-TDDFT calculations make use of in-house extensions to the Vienna *ab initio* simulation package (VASP)^{438–440} first reported in Ref. 189, and we consider adiabatic PBE and SCAN exchange-correlation potentials. The plasmon for both RT-TDDFT functionals is notably red-shifted relative to the experiment, similar to the LR-TDDFT ALDA results reported in Ref. 167, but inconsistent with LR-TDDFT results for cold solid aluminum reported in Ref. 441], which also indicated better agreement (and a small blue shift) between ALDA and experiment. The DFT-based theory to which we compare from Ref. 145 uses separate models for the bound-free edge (the $q\u21920$ limit of the optical conductivity) and the plasmon (the Mermin approximation, with collision frequencies taken from Kubo–Greenwood calculations^{442}), but both are based on thermal DFT in a limit that is equivalent to an LR-TDDFT calculation in which *f _{Hxc}* = 0 and both make use of the HSE exchange-correlation functional in the attendant ground-state calculations.

^{443}To capture the bound-free edge, we used an 11-electron projector augmented-wave pseudization

^{444}of the electron–ion interaction. Relative to a three-electron pseudization, the additional spectral width adversely impacts the condition number of the linear system of equations in the time propagation scheme, making HSE impractical relative to simply verifying expected systematics for much cheaper PBE and SCAN calculations. We expect that the only impact of using HSE would be to more strongly bind the 2p states and to shift the bound-free edge a bit further than SCAN, consistent with prior DFT-MD calculations.

^{445}This is consistent with the systematics that are observable in switching between PBE and SCAN in RT-TDDFT, in which no change in the shape of the plasmon is observable.

The RT-TDDFT DSF for a range of plausible ionic structures is convolved with a fit to the reported instrument function, and while the simulated scattering intensity is in reasonable agreement with experiment, the discrepancy between RT-TDDFT and the methods based strictly on DFT is notable. The bound-free edge is approximately consistent with Kubo–Greenwood at $q\u21920$, which suggests that this region of the DSF should be well treated within a ground-state framework. However, the plasmon is notably distorted relative to the *ab initio* formulation of the Mermin approximation. We note that this has previously been observed,^{189} and more severe discrepancies between Kubo–Greenwood and RT-TDDFT have been observed in the context of simulations of bound-bound features in XRTS.^{258} In as far as RT-TDDFT is ostensibly a higher level of theory, this warrants further investigation to reach a better understanding of the relative strengths/weaknesses of the two approaches as well as ultimate consistency.

### C. Nonlinear density response

The concept of linear density response is ubiquitous throughout physics,^{327} in general, and WDM description, in particular. From a theoretical perspective, neglecting all nonlinear terms, e.g., in Eq. (52), leads to substantial simplifications, which often render calculations feasible in the first place. Indeed, the assumption of a linear response is well justified for many applications,^{167} such as XRTS experiments with low to moderate intensities as they are employed for WDM diagnostics.^{147} At the same time, nonlinear effects are known to play an important role in high-Z stopping power calculations^{181–183} and have been suggested to be important in the construction of effective pair potentials.^{139,179,180} To be more specific, the Bethe stopping power formula, together with its many corrections, remains the cornerstone of our understanding of electron ionization-excitation energy losses of charged heavy particles in matter.^{446} One of the primary limitations of the bare (uncorrected) Bethe formula is reflected on its origin from the first-order quantum mechanical perturbation theory.^{447} In fact, for high-Z particles, it can be expected that higher-order terms of the Born series become significant.^{448,449} Translating from the language of quantum scattering theory to the language of density response theory, the stopping power of the high-Z particle is the force that it experiences from its own induced field in the medium.^{450} Thus, it is straightforward to deduce that the next high-order correction will roughly be $\u221dZ3$ and will emerge from the quadratic density response.^{182,451,452} This is the famous Barkas effect and corresponding Barkas correction that also accounts for the stopping power differences between equal mass projectiles of opposite charge sign.^{453,454}

In addition, it can be shown^{176} that nonlinear response functions are directly connected to higher-order correlation functions. Therefore, the experimental investigation of nonlinear observables might allow one to gain new insights into many-body effects in a given system. Finally, we note that the nonlinear density response has been shown^{167,175} to depend much more sensitively on system parameters like the temperature compared to the usually considered linear response, which might make them useful as a new tool for WDM diagnostics.^{184}

To our knowledge, the first rigorous analysis of the static nonlinear density response of WDM has been presented in Ref. 167 based on exact PIMC simulations of a harmonically perturbed UEG, cf. Eq. (29). The basic idea can be seen in Fig. 4(d), where we show the induced density $\u27e8\rho \u0302k\u27e9q,A$ [Eq. (37)] as a function of the perturbation amplitude *A*. More specifically, the symbols show the raw PIMC results, and the dotted curves fit functions based on the expansion given in Eq. (60). Note that the induced density has been divided by *A*, so that the datasets attain a constant value in the linear response limit. For the smallest depicted value of the wave number *q* (blue diamonds), the density response of the UEG is comparably small, and nonlinear effects are basically negligible over the entire depicted *A*-range. In stark contrast, the induced density for $q=3qmin$ is substantially larger in magnitude and visibly deviates from the constant LRT limit for $A\u22730.02$. This is reflected by a large value in the cubic coefficient in Eq. (60), which directly corresponds to the cubic density response at the first harmonic, i.e., at the third power of the perturbation amplitude but exactly at the perturbation wave vector (and, in general, frequency). Detailed investigations of the latter have been presented in Refs. 167 and 175.

While being formally exact, the approach based on the direct PIMC simulation of the harmonically perturbed system is computationally very demanding as it requires one to carry out a set of independent simulations over a sufficient interval of perturbation amplitudes *A* to extract the nonlinear density response for a single combination of *r _{s}*,

*θ*, and

*q*. A promising route to circumvent this issue has been suggested by Moldabekov

*et al.*,

^{455}who have shown that the nonlinear density response can be studied very accurately based on KS-DFT simulations. In particular, this will allow one to transcend the current limitations of PIMC and to study nonlinear effects in real materials. From the perspective of PIMC itself, Dornheim

*et al.*

^{195}have presented relations between generalized nonlinear response functions and high-order imaginary-time correlation functions, which can be viewed as generalizations of the imaginary-time version of the fluctuation-dissipation theorem given in Eq. (40). In the present work, we restrict ourselves to quadratic orders in the perturbation strength

*A*, and the corresponding relation between the quadratic density response function at the second harmonic and the respective imaginary-time three-body correlation function $F(2)(q;\tau 1,\tau 2)$ is given by Eq. (58).

In Fig. 14, we show new *ab initio* PIMC results for $F(2)(q;\tau 1,\tau 2)$ in the *τ*_{1}-*τ*_{2}-plane for two different values of the wave number *q* for the UEG at *r _{s}* = 4 and

*θ*= 1, i.e., for the same conditions as in Fig. 4. Overall, we find that $F(2)(q;\tau 1,\tau 2)$ exhibits a qualitatively similar behavior to the two-body ITCF $F(q,\tau )$. Specifically, we observe an increasingly steep decay with respect to both

*τ*

_{1}and

*τ*

_{2}with increasing

*q*, whereas the three-body ITCF becomes flatter in the long-wavelength limit. At the same time, we note that the physical behavior of $F(2)(q;\tau 1,\tau 2)$ remains very poorly understood, and its improved understanding along similar lines as it has been achieved for $F(q,\tau )$ in Refs.190,192 remains an important task for future works.

In the present context, the main utility of $F(2)(q;\tau 1,\tau 2)$ is given by its direct connection to the quadratic static density response function of the second harmonic, $\chi (2,2)(q)$, via Eq. (58), which is explored in Fig. 15 for the same conditions as in Fig. 14. In particular, the green crosses have been obtained from a single simulation of the unperturbed UEG in this way. Evidently, $\chi (2,2)(q)$ exhibits the opposite sign compared to the linear static density response function $\chi (q)$ shown in Fig. 4, which, together with the positive sign of the cubic response at the first harmonic, leads to a saturation of the true response compared to LRT for medium to large perturbation amplitudes *A*. In addition, we find that the PIMC-based evaluation of Eq. (58) is in excellent agreement with the five red dots, which have been obtained from a multitude of PIMC simulations of the harmonically perturbed system by fitting Eq. (61) to the induced density at the second harmonic, i.e., $\u27e8\rho \u03022q\u27e9q,A$. In practice, the imaginary-time framework introduced in Ref. 195, thus, gives superior results for the quadratic density response of the second harmonic for all relevant wave numbers with a fraction of the computation cost. For completeness, we note that similar relations also exist for the nonlinear interaction between multiple external perturbations, which has been explored in more detail in Ref. 177.

In addition to their direct value for the understanding of nonlinear effects in WDM, the availability of exact PIMC reference data is also pivotal to guide the development of new theoretical frameworks. In this regard, a highly useful result is given by the recursion formula for the quadratic density response of the ideal Fermi gas by Mikhailov,^{173,174} cf. Eq. (64). The results are shown as the solid yellow curve in Fig. 15 and exhibit the correct asymptotic behavior for $q\u22732qF$. Screening effects for small *q*, on the other hand, are not included. In Ref. 175, five of us have introduced an RPA-like expression for $\chi (2,2)(q)$—and also for high-order response functions which are beyond the scope of the present overview—that is truncated at the linear level with regard to screening effects, see Eq. (65). The results have been included as the dashed black curve. They exhibit the correct trends in the limits of $q\u21920$ and $q\u2192\u221e$ but substantially underestimate the true quadratic response around its maximum at $q\u223c2qF$. Including the static local field correction $G(q;rs,\theta )$ evaluated from the neural-net representation from Ref. 194 via Eq. (66) gives the dotted blue curve, which is in perfect agreement with the PIMC datasets over the entire *q*-range. This is a strong empirical verification for the small effect of the screening truncation employed in Eqs. (65) and (66). Moreover, it confirms the more pronounced importance of electronic XC effects for the description of nonlinear effects, since the RPA curve is less accurate for the description of $\chi (2,2)(q)$ compared to $\chi (q)$; this is consistent with previous investigations presented in Refs. 175 and 176.

Future investigations of the nonlinear density response of the UEG will include the study of dynamic effects, which is possible in multiple ways. For example, the RPA and LFC-based expressions for the quadratic density response Eqs. (65) and (66) remain the same for all frequencies. As a first step, one might, therefore, study dynamic nonlinear effects based on the static approximation, which is straightforward given the easy availability of $G(q;rs,\theta )$ based on different representations.^{72,194} Moreover, these efforts can be improved by including the full, frequency-dependent results for $G(q,\omega )$ that are available for certain parameters of the UEG based on the analytic continuation^{252} discussed in Sec. III B 1. Such an approach will allow one to gain both qualitative and quantitative insights into a number of interesting phenomena, such as double plasmons^{250,456,457} and dynamic three-body correlations.^{458–460}

In addition, upcoming direct PIMC simulations (i.e., without the fixed-node approximation) of real WDM systems such as hydrogen^{77,78} will allow one to directly estimate $F(2)(q;\tau 1,\tau 2)$ as well as higher-order ITCFs that are connected to the cubic density response and so on.^{195} First and foremost, this will facilitate investigations of the static nonlinear density response of WDM, without any assumptions or approximations. In particular, the presence of both electrons and ions leads to various cross terms, which deserve a closer examination. Finally, we envision the systematic development of the study of dynamic many-body effects in the imaginary-time domain similar to the efforts that have been developed in Refs. 146, 190, and 192, for $F(q,\tau )$, which have been touched upon throughout the present work.

## IV. INTERPRETATION OF XRTS EXPERIMENTS

One of the most important practical applications of the electronic density response of WDM is given by the interpretation of XRTS experiments. In the traditional way, one constructs a theoretical model for the DSF $S(q,\omega )$ with a-priori unknown system variables such as the temperature *T*, the density *n*, or the charge state *Z* being the free parameters.^{147} Convolving the DSF with the combined source and instrument function $R(\omega )$, cf. Eq. (3), then allows one to compare these results with the experimentally measured scattering intensity $I(q,\omega )$, and the best fit between the model and the experiment then gives one access to the system parameters of interest. We note that the numerical deconvolution of Eq. (3) is generally rendered highly unstable by the inevitable experimental noise. Therefore, XRTS does not give one direct access to the DSF, which could have been used, e.g., to extract the temperature via the detailed balance relation Eq. (46).^{319}

To our knowledge, the most widely used theoretical model for the interpretation of XRTS experiments with WDM is given by the Chihara decomposition,^{147,148,151,186} which is based on a separation into bound and free electrons. Yet, such a chemical picture is expected to break down for significant parts of the WDM regime, where the electrons are partially localized around the ions; in this case, they are neither bound nor free, see also the discussion of warm dense hydrogen in Sec. III A 2.

We note that overcoming these difficulties with more advanced simulation methods such as LR-TDDFT is also no trivial matter for numerous reasons. First and foremost, such *ab initio* simulations are computationally very costly and are unclear if it will be feasible to carry them out on a sufficient grid of free parameters to infer thermodynamic variables such as the temperature from XRTS measurements. In this regard, modern machine-learning based interpolation methods^{114} will likely be helpful, but their practical application to this problem remains an important task for future research. Second, the accuracy of LR-TDDFT simulations is determined by the XC kernel, and the development of consistent useful approximations for real WDM systems is still in its infancy.

In any case, it is safe to say that previous approaches for the interpretation of XRTS measurements are based on approximations, and their accuracy remains largely unclear. Therefore, EOS measurements^{187,188} are actually model-dependent and do not necessarily constitute a reliable baseline either for practical applications or to guide the development of improved theoretical methods.

Very recently, Dornheim *et al.*^{146,190,191} have suggested that this problem can be circumvented by switching from the usual frequency-domain that is centered around $S(q,\omega )$ to the imaginary time *τ*; here, the main property is the ITCF $F(q,\tau )$, which contains exactly the same information as the DSF. In particular, the convolution theorem Eq. (45) implies that the deconvolution is trivial in the Laplace domain. In other words, XRTS gives us direct access to $F(q,\tau )$, which can subsequently be used to infer information about the given system without any intermediate models or simulations. For example, the symmetry relation Eq. (47) indicates that the ITCF is symmetric around $\tau =\beta /2=1/2T$. In this sense, one could say that XRTS measurements effectively function as a thermometer as the relation between $F(q,\tau )$ and *T* is trivial.

In Fig. 16, we demonstrate this idea for the pioneering observation of plasmons in warm dense beryllium by Glenzer *et al.*^{202} Panel (a) shows the actual experimental scattering intensity (green) and the corresponding combined source and instrument function (blue). In particular, the elastic feature around *ω* = 0 is dominated by the latter, and both the up- and down-shifted plasmons are nicely visible around $\omega \u2248\xb125$ eV. From a practical perspective, it is clear that the evaluation of the Laplace transform Eq. (41)—and also the proof of the convolution theorem, Eq. (45)—requires an integration over an infinite frequency-range, whereas the available *ω*-range is clearly limited in any given experiment. In practice, we, thus, define a symmetrically truncated Laplace transform,

with the corresponding truncated ITCF

It is easy to see that the truncation error vanishes in the limit of large *x*

and the convergence with *x* has to be carefully checked.

This is demonstrated in Fig. 16(b), where we plot the position of the minimum in $Fx(q,\tau )$ as a function of the integration boundary. We note that the shaded green area constitutes a measure of the respective uncertainty due to the experimental noise, and the utilized procedure is explained in detail in Ref. 191. Evidently, we observe a convergence of the, thus, inferred temperature around $\omega \u227330$ eV, i.e., upon capturing the main plasmon peak in the integration range. As our final estimate for the temperature, we find $T=14.8\xb12$ eV. This is comparably close to the nominal temperature of $T=12$ eV, which was found in the original publication^{202} based on a simple Mermin model.^{461}

Let us conclude this reassessment of the beryllium XRTS data by examining the actual result for $F(q,\tau )$, shown in Fig. 16 for the converged value of $x=40$ eV. Specifically, the green line shows the deconvolved ITCF (up to the unknown normalization factor), and the shaded gray area is the corresponding uncertainty interval due to the experimental noise. The vertical yellow line indicates the location of the minimum, which has been used to infer the temperature via Eq. (47). Moreover, the dotted blue line shows the Laplace transform of the instrument function $R(\omega )$, i.e., the denominator of the RHS of Eq. (45). Without this correction, the minimum of the ITCF would be shifted to smaller values of *τ*, leading to an overestimation of the temperature *T.*^{146,191} Finally, the dashed red curve has been obtained by mirroring the ITCF around $\tau =\beta /2$. Evidently, the red curve is in excellent agreement with the green curve, which means that the ITCF computed from the XRTS data shown in Fig. 16(a) is really symmetric to a remarkable degree given the comparably high noise level in the latter. The stability of the method with respect to experimental noise is analyzed in detail in Ref. 191.

In summary, the formally exact method proposed in Ref. 146 allows one to extract very accurate results for the temperature of arbitrary complex systems in thermodynamic equilibrium from an XRTS measurement without any models, approximations, or simulations. From a practical perspective, a decisive factor is given by the accurate quantification of the combined source and instrument function $R(\omega )$. This is relatively straightforward using source monitors at modern XFEL facilities,^{18,55} whereas the characterization of backlighter emission spectra^{462} at facilities like NIF is more challenging. On the other hand, we note that the high temperatures typical for NIF implosion experiments make the influence of $R(\omega )$ on the extracted temperature less pronounced in any case. Furthermore, it has been shown^{191} that the characteristic width of the instrument function also determines the minimum temperature that can be inferred from the ITCF method, which becomes apparent both from considerations of the respective DSF and the corresponding ITCF. Specifically, the detailed balance relation Eq. (46) of the DSF directly implies that the signal at negative frequencies is exponentially damped with increasing *β*, i.e., upon lowering the temperature. At some point, the signal for $\omega <0$ will vanish within the given experimental noise, and no inference of a temperature will be possible. In the *τ*-domain, large *β* means that one has to resolve the ITCF over a correspondingly large imaginary-time interval. With increasing *τ*, the effect of the convolution becomes exponentially more pronounced. In this case, every small uncertainty in the determination of $R(\omega )$ will have increasing consequences for the evaluation of Eq. (45) and, thus, make the localization of the minimum around $\tau =\beta /2$ less accurate.

Additional challenges for future works include the investigation of the effect of inhomogeneity onto the ITCF, which is particularly important for spatially extended systems such as the fuel capsule in implosion experiments at the NIF.^{463,464} Moreover, the symmetry relation Eq. (47) naturally only holds in thermodynamic equilibrium, and the influence even of small nonequilibrium effects^{270,465} should be quantified in the future. As a final point, we mention the potential impact of the (small) frequency-dependence of the wave vector, $q=q(\omega )$, even though a first investigation in Ref. 191 has led to the conclusion that the effect can likely be neglected for XRTS experiments and is of no consequence for the interpretation of the beryllium data shown in Fig. 16.

## V. SUMMARY AND OUTLOOK

In this work, we have given an overview of our current understanding of the electronic density response of WDM. More specifically, we have attempted to form a coherent picture that elucidates the connections between various theoretical approaches and simulation methods such as PIMC and different flavors of DFT, and between different representations such as the DSF $S(q,\omega )$ and the ITCF $F(q,\tau )$. In addition, we have discussed both the ubiquitous weak-perturbation limit described by linear response theory and its extension to include nonlinear effects with respect to the perturbation amplitude.

From a theoretical perspective, density response theory constitutes a powerful framework for the investigation of a host of observables. Indeed, it forms the basis for the widely used LR-TDDFT method and helps to connect approaches such as dielectric theories with PIMC and DFT. In addition, the density response theory is pivotal for the understanding and modeling of experiments, with XRTS measurements of WDM being a prime example. In this regard, special attention has been given to the imaginary-time domain, which combines a number of key advantages. First, obtaining the ITCF $F(q,\tau )$ from XRTS experiments is straightforward; this allows for exact and model-free diagnostics of the temperature^{146} and can likely be extended to other observables in future works.^{190} Second, the ITCF contains the same information as the DSF, which makes it a useful tool for the investigation of dynamic effects from exact PIMC simulations.^{190} Finally, we argue that the ITCF is at the heart of density response theory as it is directly connected to $\chi (q)$ and $S(q,\omega )$ in different ways. In fact, we view the imaginary-time domain as a complementary representation of dynamic quantum many-body theory. As such, it is reasonable to consider it even in cases where frequency-dependent properties such as the DSF are exactly known, since different representations tend to emphasize different aspects of the same information.

Let us conclude this work by outlining three particularly promising routes for future research.

### A. Improved real-time TDDFT and NEGF simulations

As we discussed in Secs. II A 3 and II A 4, real-time simulations offer a number of attractive features: direct access to nonlinear response properties, to non-adiabatic effects as well as to high-level electronic correlations. The observation that real-time NEGF simulations with rather simple self-energies allow one to reproduce linear response results that, otherwise, require to solve complicated Bethe–Salpeter-type equations indicates a promising route for future developments in NEGF and real-time TDDFT. This can benefit from recent developments in the NEGF theory that allowed to dramatically speed up the simulations achieving time-linear scaling^{372,373} also for advanced self-energies.^{279} These simulations can be used, among other things, to benchmark existing and derive improved exchange-correlation potentials for RT-TDDFT. In this context, we also mention the recently developed quantum fluctuations approach,^{466} which constitutes another promising concept for the real-time computation of the density response.

### B. Systematic improvement of LR-TDDFT

In Fig. 17, we illustrate the recent framework by Moldabekov *et al.*^{255} for the DFT-based estimation of the static XC kernel $Kxc(q)$. As the first step (top left), KS-DFT calculations are carried out with respect to the original electronic Hamiltonian $H\u03020$ and with respect to a modified Hamiltonian $H\u0302=H\u03020+V\u0302ext$ that is perturbed by a monochromatic static potential $V\u0302ext$, cf. Eq. (29). The difference between these two calculations gives one the induced density shown in the bottom left panel. In the limit of small *A*, one gets direct access to the physical static density response $\chi (q)$ and, in this way, to the static XC kernel $Kxc(q)$ of real materials (top right).

As the first future improvement, it will be indispensable to compare these DFT results for different conditions and for a variety of XC functionals across Jacob's Ladder to upcoming exact PIMC simulations of real WDM, starting with hydrogen. It can be expected that a comparison of the static XC kernel by itself will give invaluable insights into the performance of different functionals and may guide the development of improved approximations. In this context, we also mention the enticing possibility to actually measure the exact $Kxc(q)$ in XRTS experiments, which is discussed in more detail below. Second, we propose to utilize the DFT results for $Kxc(q)$ as input for LR-TDDFT simulations within the static approximation.^{252} This can be motivated by the high accuracy of this approach for the UEG at metallic densities and should work particularly well for highly ionized matter. The corresponding results for $S(q,\omega )$ can then a) be used to interpret and model XRTS experiments [Eq. (3)] and, after performing the Laplace transform Eq. (41), compared to exact PIMC reference data for $F(q,\tau )$, e.g., for hydrogen.

In fact, the last step with the DFT results can be viewed as a linear response imaginary-time-dependent DFT (LR-iTDDFT) calculation within the respective static approximation. From a theoretical perspective, iTDDFT has a crucial advantage over usual TDDFT methods that operate in frequency-space: exact PIMC calculations can be used to extract an exact dynamic XC kernel that contains the full dependence on the complex frequency *z*, $Kxc(q,z)$. Such information can give key insights into the limits of the static approximation. Moreover, it can form the basis for the construction of novel dynamic XC kernels for LR-iTDDFT calculations of real materials that transcend the current limitations of TDDFT. We re-iterate that such calculations will allow one to accurately predict the ITCF $F(q,\tau )$, which is directly accessible from XRTS experiments.

Finally, we propose to use either the aforementioned LR-TDDFT calculations within the static approximation or improved LR-iTDDFT calculations with a dynamic, z-dependent kernel as the basis for the estimation of electron–electron correlation functions such as the static structure factor $S(q)$ or the electron–electron pair correlation function $g(r)$. From a philosophical perspective, the key ingredient to this idea is given by the fluctuation-dissipation theorem [Eq. (31)] that provides a straightforward relation between the induced single-particle density—a property that is easily accessed within standard KS-DFT—and an electron–electron correlation function such as the DSF $S(q,\omega )$. In this way, it will be possible to reconstruct information about electronic correlations; we also note that this idea is not limited to pair correlations and can be extended to higher-order correlation functions^{176} by using KS-DFT to study the nonlinear density response of a given system of interest.^{455}

### C. Improved XRTS measurements at modern XFEL facilities

The proposed systematic improvements of our theoretical understanding of the electronic density response of WDM can be decisively aided by new emerging experimental capabilities. Indeed, the new framework for the interpretation of XRTS experiments with WDM in the imaginary-time domain allows us to suggest new experimental setups that exploit these ideas in an optimal way. While this applies to XRTS experiments in general, we find that modern XFEL facilities such as the European XFEL^{55} are particularly promising in this regard, as they offer a high degree of flexibility for new developments, and for the exploration of novel concepts and ideas. In practice, a key advantage of the European XFEL is given by the rep-rated high-power drive laser systems ReLaX and DiPOLE100X, which are made available by the HIBEF user consortium and match the macro-bunch frequency of the x-ray laser of 10 Hz.^{467} This increase in the data rate of up to several thousand times in comparison to similar installations at LCLS and SACLA leads to a revolutionary improvement in photon statistics,^{468} which can result in high-quality scattering spectra with a dynamic range over four orders of magnitude or more. This accuracy is particularly important for clear measurements of the blue-shifted part of the scattering spectrum, which is exponentially damped by the Boltzmann factor in detailed balance. Furthermore, the use of self-seeding, both with and without a monochromator, can substantially reduce the bandwidth of the instrument function while keeping the photon flux in a comparable realm as for pure SASE radiation or, indeed, even better.^{469} This reduces the influence of uncertainties of the instrument function for further analysis.

A clear example of a first experiment to test the proposed imaginary-time framework would be the measurement of XRTS of high energy density matter with unprecedented quality for an entire set of scattering angles. This would give us access to $I(q,\omega )$ over the full relevant *q*-range and will be beneficial for two reasons. First, this would allow the ITCF-based temperature diagnostics developed in Ref. 146 to be applied for all wave vectors. Since the extracted temperature has to be the same independent of a particular *q*, such an experiment can be used to further demonstrate the consistency of the idea and to quantify the underlying uncertainty. Second, recording the XTRS spectrum for multiple *q* will be indispensable to test theoretical models and to guide the development of methodological improvements. In Sec. III A 1, we have made extensive use of the imaginary-time version of the fluctuation-dissipation theorem, Eq. (40), which gives a straightforward relation between $F(q,\tau )$ and the static density response function $\chi (q)$. The proposed set of XRTS measurements will, thus, make it possible to measure the static density response function of WDM systems. In addition to being very interesting in their own right, these data can be used to invert Eq. (17) and, in this way, to obtain unassailable experimental reference data for the XC kernel of real materials.

## ACKNOWLEDGMENTS

This work was partially supported by the Center for Advanced Systems Understanding (CASUS), which is financed by Germany's Federal Ministry of Education and Research (BMBF) and by the Saxon state government out of the State budget approved by the Saxon State Parliament. M.B. acknowledges funding by the Deutsche Forschungsgemeinschaft via Project No. BO1366-15. A.D.B. acknowledges support from Sandia's Laboratory Directed Research and Development Program and U.S. Department of Energy Science Campaign 1. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA0003525. The work of Ti.D. and F.G. was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344.

The PIMC and DFT calculations were partly carried out at the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN) under Grant No. shp00026, on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden, and on the Hemera cluster at Helmholtz-Zentrum Dresden-Rossendorf (HZDR).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Tobias Dornheim:** Writing – original draft (lead); Writing – review & editing (lead). **Tilo Doeppner:** Writing – review & editing (equal). **Frank Reno Graziani:** Writing – review & editing (equal). **Michael Bonitz:** Writing – original draft (equal); Writing – review & editing (equal). **Attila Cangi:** Writing – original draft (equal); Writing – review & editing (equal). **Jan Vorberger:** Writing – original draft (equal); Writing – review & editing (equal). **Zhandos Moldabekov:** Writing – original draft (equal); Writing – review & editing (equal). **Kushal Ramakrishna:** Writing – original draft (equal); Writing – review & editing (equal). **Panagiotis Tolias:** Writing – original draft (equal); Writing – review & editing (equal). **Andrew Baczewski:** Writing – original draft (equal); Writing – review & editing (equal). **Dominik Kraus:** Writing – original draft (equal); Writing – review & editing (equal). **Thomas Robert Preston:** Writing – original draft (equal); Writing – review & editing (equal). **Dave A. Chapman:** Writing – review & editing (equal). **Maximilian Böhme:** Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.