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.

Warm dense matter (WDM) is an extreme state that is characterized by high temperatures (T103108 K), high pressures (P1104 Gbar), and densities in the vicinity of and even partially exceeding solid state densities (n10211028 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 

FIG. 1.

The warm-dense matter (WDM, orange bubble) regime is schematically illustrated in the density-temperature plane. Characteristic parameters are given by the Wigner–Seitz radius rs (vertical blue) and the degeneracy temperature θ=kBT/EF (diagonal green), also shown is the effective coupling parameter of the electrons Γ (solid red), see Ref. 26. The gray bubbles indicate conditions of various astrophysical objects taken from Ref. 1, and the black line sows the implosion path of a DT fuel capsule taken from Hu et al.27 Adapted from Refs. 28 and 29.

FIG. 1.

The warm-dense matter (WDM, orange bubble) regime is schematically illustrated in the density-temperature plane. Characteristic parameters are given by the Wigner–Seitz radius rs (vertical blue) and the degeneracy temperature θ=kBT/EF (diagonal green), also shown is the effective coupling parameter of the electrons Γ (solid red), see Ref. 26. The gray bubbles indicate conditions of various astrophysical objects taken from Ref. 1, and the black line sows the implosion path of a DT fuel capsule taken from Hu et al.27 Adapted from Refs. 28 and 29.

Close modal

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 materials17,30–34 and hot-electron chemistry.35–37 A particularly important application is given by inertial confinement fusion38,39 as it is realized at the National Ignition Facility (NIF).40 On its pathway toward ignition, the fuel capsule traverses the WDM regime27 (see the black line in Fig. 1 that indicates the implosion path of a deuterium–tritium [DT] capsule27).

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 pressure17,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 Facility50–52 and the OMEGA laser53 in the USA, SACLA54 in Japan, and the European XFEL55 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 θ=kBT/EF measures the thermal energy in units of the electronic Fermi energy EF,59 with θ1 and θ1 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 Γ1 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 hydrogen77,78 over substantial parts of the WDM regime. Moreover, path-integral Monte Carlo (PIMC) simulations within the fixed-node approximation79—also known as restricted PIMC (RPIMC) in the literature—have allowed Militzer and co-workers80–82 to go to heavier elements and material mixtures.83 These simulations constitute the basis for an extensive equation-of-state (EOS) database84 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 UEG28,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 models119–121 have matured into production-level codes that accelerate DFT calculations for realistic materials at finite temperature and pressure122,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 properties141 such as stopping power109,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 temperature146 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,ω), 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 

S(q,ω)=dteiωtF(q,t),
(1)

which is defined as

F(q,t)=n̂(q,t)n̂(q,0)
(2)

for homogeneous systems in equilibrium. For sufficiently weak probes in the linear response regime, the DSF is related to the linear density response function χ(q,ω) through the fluctuation-dissipation theorem59 (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,ω) in terms of the same momentum and energy transfer as those parametrizing S(q,ω) and χ(q,ω), 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(ω),

I(q,ω)=S(q,ω)R(ω).
(3)

Nevertheless with careful characterization of R(ω), information about S(q,ω) and χ(q,ω) 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 χ(q,ω) as153 

χ(q,ω)=χ0(q,ω)14πq2[1G(q,ω)]χ0(q,ω),
(4)

with χ0(q,ω) 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 χ(q,ω) 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,ω). Indeed, setting G(q,ω)0 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,ω)=4π/q2G(q,ω). Indeed, replacing χ0(q,ω) by the dynamic Kohn–Sham response function χS(q,ω) 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 effects168–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 distances139,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 regime77 and constitutes a de facto uncontrolled approximation in practice. Therefore, previous EOS measurements14,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 mechanics193 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.

We consider an N-particle system that is described by the sum of an unperturbed Hamiltonian Ĥ0 and an external single-particle perturbation that we give in the most general form,

Ĥ=Ĥ0+V̂q,ω,A,Ĥ0=K̂+Ŵ+V̂pot,
(5)
V̂q,ω,A=2Af(t)l=1Ncos(q·r̂lγωt).
(6)

Here, K̂,Ŵ, and V̂pot are the kinetic, interaction, and external Coulomb energy (e.g., due to a configuration of ions) of the electrons, respectively. Moreover, V̂q,ω,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̂q,ω,A and V̂pot contribute to the total external single-electron potential energy V̂ext=V̂pot+V̂q,ω,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.

TABLE I.

Overview of possible excitation scenarios contained in the Hamiltonian Eq. (5). More details on linear-response theory (cases 1 and 3) can be found in Sec. II B; static nonlinear response theory (case 4) is discussed in Sec. II D, and its generalization to the dynamic nonlinear regime (case 2) has been presented in Ref. 175. The theoretical formalism behind cases 5 and 6 is discussed in Sec. II A 4.

f(t)A0γRegimeApproach
Small Dynamic linear response LR-TDDFT (Sec. II A 2), PIMC (Sec. II A 1) + analytic continuation (Sec. II C
Large Dynamic nonlinear response Dynamic nonlinear response theory (Ref. 175
Small Static linear response PIMC (Sec II A 1), KS-DFT (Sec. II A 2), dielectric theories (Refs. 163 and 196–201
Large Static nonlinear response PIMC (Sec. II A 1), KS-DFT (Sec. II A 2
A0δ(t) Small Weak kick Real-time TDDFT (RT-TDDFT) (Sec. II A 3), NEGF (Sec. II A 4
A0δ(t) Large Strong kick RT-TDDFT (Sec. II A 3), NEGF (Sec. II A 4
f(t)A0γRegimeApproach
Small Dynamic linear response LR-TDDFT (Sec. II A 2), PIMC (Sec. II A 1) + analytic continuation (Sec. II C
Large Dynamic nonlinear response Dynamic nonlinear response theory (Ref. 175
Small Static linear response PIMC (Sec II A 1), KS-DFT (Sec. II A 2), dielectric theories (Refs. 163 and 196–201
Large Static nonlinear response PIMC (Sec. II A 1), KS-DFT (Sec. II A 2
A0δ(t) Small Weak kick Real-time TDDFT (RT-TDDFT) (Sec. II A 3), NEGF (Sec. II A 4
A0δ(t) Large 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 χ(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=r1r2 and t=t1t2, one usually studies χ(q,ω).

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)A 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 models204–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

In Sec. II A, we focus on four particularly important methods, namely, path-integral Monte Carlo (Sec. II A 1), thermal DFT (Sec. II A 2), real-time time-dependent DFT (Sec. II A 3), and nonequilibrium Green functions (Sec. II A 4).

1. Path-integral Monte Carlo

The ab initio path-integral Monte Carlo (PIMC) approach221–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 β=1/kBT, volume V=L3, and number density n=N/V are constant) with N and N spin-up and spin-down electrons, is given by

Zβ,N,V=1N!N!σSNσSNsgn(σ,σ)VdRR|eβĤ|π̂σπ̂σR.
(7)

Here, the summation over all elements σi of the respective permutation group SNi, with π̂σ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̂) and potential (V̂=Ŵ+V̂ext) contributions to the full Hamiltonian Ĥ=K̂+V̂, which renders the matrix elements of the density operator ρ̂=eβĤ 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 formula228 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=iτ over a discrete time step τ=β/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 implementations72,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 literature235,236), the sign function in Eq. (7) is always positive, and the PIMC method allows for the quasi-exact simulation of up to N104 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 (θ0.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 Ceperley79 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,τ), see Sec. II C.

2. Thermal density functional theory

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

Ĥ=K̂+Ŵ+V̂ext,
(8)

with K̂ being the kinetic part, Ŵ containing the electron–electron interaction, and V̂ext an external single-particle potential. In particular, V̂ext 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̂ext can also include the monochromatic external perturbation that is used to obtain the electronic density response, cf. Eq. (6) in the limit of ω0 and f(t)1. 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 generalization85 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,

Ω[n]=FT[n]+Vext[n],
(9)

where

FT[n]=minΓ̂n(K[Γ̂]+S[Γ̂]+W[Γ̂])
(10)

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 Γ̂ 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

Ω=minnΩ[n]=minn(FT[n]+Vext[n]).
(11)

Minimizing Ω[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,

(22+vS[n](r))ϕα,k(r)=ϵα,kϕα,k(r),
(12)

where ϕα,k(r) and ϵα,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

Ω[n]=KS[n]+TSS[n]+U[n]+ΩXC[n]+Vext[n],
(13)

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

n(r)=α,kfα,k(β,μ)|ϕα,k(r)|2
(14)

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 {ϕα,k} are populated according to the Fermi–Dirac distribution

fα,k(β,μ)=[1+eβ(ϵα,kμ)]1,
(15)

with β=1/kBT and μ is the chemical potential.59 In combination, Eqs. (12) and (15) give one direct access to KS[n]=1/2α,kdrϕα,k*(r)2ϕα,k(r) and SS[n]=α,k(fα,kln(fα,k)+(1− fα,k)ln(1fα,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)+ vH[n](r)+vXC[n](r) in terms of the external potential vext(r), the Hartree potential vH[n](r)=δU[n]/δn(r), and the XC potential vXC[n](r)=δΩXC[n]/δn(r). In particular, the XC contribution Ω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, Ω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 functionals90–93,255 and the construction of novel, more sophisticated approximations to Ω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 {ϕα} 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 function162,248 given by

χGGS(q,ω)=1Vk,α,αfαkfαk+qω+ϵαkϵαk+q+iη×ϕαk|ei(q+G)·r|ϕαk+qϕαk|ei(q+G)·r|ϕαk+q.
(16)

Here, the parameter 0<η1 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 are the reciprocal lattice vectors. We are here mainly interested in the macroscopic response functions, i.e., setting G=G=0.

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

χ(q,ω)=χS(q,ω)1[v(q)+Kxc(q,ω)]χS(q,ω),
(17)

with Kxc(q,ω) 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,ω) strongly depends on the particular form of the KS response function χS(q,ω). 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 {ϕα}. 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,ω)0 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 χ(q,ω) gives one direct access to a number of material properties, such as the dynamic structure factor S(q,ω), 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,

itϕα,k(r,t)={22+vS[n](r,t)}ϕα,k(r,t),
(18)

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 ĤS(t)=2/2+vS(r,t). This is typically posed as an initial value problem in which ϕα,k(r,t) is known for some t = ti, typically from the solution to a static Mermin–Kohn–Sham DFT problem. The time-dependent density is

n(r,t)=α,kfα,k(β,μ)|ϕα,k(r,t)|2,
(19)

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 vS, vext describes the attractive electron–ion interaction as well as impressed perturbations that drive the dynamics, vH is the time-dependent Hartree potential, and vxc 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.

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 δv0f(t), where δ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 effects258,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,

U(tf,ti)=Texp[ititfdτĤS(τ)],
(20a)
whereϕα,k(r,tf)=U(tf,ti)ϕα,k(r,ti),
(20b)

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 Δt=(tfti)/N, over which the effective Hamiltonian is approximated as a constant at the midpoint of each step, eliminating the need for time ordering,

U(tf,ti)k=0N1Uap[ti+(k+1)Δt,ti+kΔt],
(21a)
whereUap(t+Δt,t)=exp[iΔtĤS(t+Δt2)].
(21b)

While conceptually straightforward, the practical implementation of this approach requires the numerical approximation of a matrix exponential and, thus, the diagonalization of ĤS 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 ĤS 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 scalable259,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)=v0exp(iq·r)f(t) can be used to compute the DSF at non-zero q,189,225 the calculation of the conductivity in the q0 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)=(1/c)(A/t), 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 σ(ω). 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

Ji(ω)=σij(ω)Ej(ω).
(22)

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

ϵij(ω)=δij+i4πσij(ω)ω.
(23)

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.

FIG. 2.

(a) Induced current density in the z direction vs time; (b) frequency-dependent conductivity component σzz.

FIG. 2.

(a) Induced current density in the z direction vs time; (b) frequency-dependent conductivity component σzz.

Close modal

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, {ϕi} and the associated creation and annihilation operators, ĉi and ĉj. The central quantity is the single-particle Green's function, which is an ensemble average of the Heisenberg form of the field operators,

gij(t,t)=iTCĉi(t)ĉj(t),
(24)

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 C275 that allows the extension of the power of Feynman diagrams to arbitrary nonequilibrium situations (for details, see the text books276,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)=igij(t,t+ϵ), 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)ψ(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),

n(1)=ig(1,1+ϵ).
(25)

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

k[itδikhikσHF(t)]gkjσ(t,t)=δC(tt)δij+kCdsΣikσcor(t,s)gkjσ(s,t),
(26)

where hHF contains kinetic, potential, and mean-field (Hartree–Fock) energy contributions, whereas correlation effects are included in the self-energy Σ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 Σcor, which depends on g rather than on the density and, therefore, depends on two times. Thus, it includes memory effects. (4) Similar to vxc in the case of DFT, Σ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 vxc and Σ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.

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)cos(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 δn(q,t) follows directly from the perturbation of the NEGF, via Eq. (25),

pδg(p,t,t)=iδn(q,t)=dt¯χ(q,t,t¯)U0(t¯),
(27)

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

FIG. 3.

Real-time NEGF simulation of the density response of the warm dense UEG following an excitation with a potential Vext(r,t)=U0(t)cos(qr), cf. dash-dotted line. Even after the excitation has vanished, the electron density exhibits plasma oscillations where frequency and damping depend on the self-energy Σcor used in the solution of the KBE (26). HF (HF2): no correlation self-energy starting from an noninteracting (interacting) initial state—this corresponds to the RPA result of the equilibrium theory; HF+ second Born: Σcor within the second Born approximation (2BA) showing correlation induced red-shift and collisional damping in addition to Landau damping. Diagrams at the bottom show the relation to equilibrium approach: the use of the 2BA for Σcor is equivalent to solving the Bethe–Salpetr equation with a two-particle kernel Ξ containing the indicated diagrams. Figure modified from Ref. 154.

FIG. 3.

Real-time NEGF simulation of the density response of the warm dense UEG following an excitation with a potential Vext(r,t)=U0(t)cos(qr), cf. dash-dotted line. Even after the excitation has vanished, the electron density exhibits plasma oscillations where frequency and damping depend on the self-energy Σcor used in the solution of the KBE (26). HF (HF2): no correlation self-energy starting from an noninteracting (interacting) initial state—this corresponds to the RPA result of the equilibrium theory; HF+ second Born: Σcor within the second Born approximation (2BA) showing correlation induced red-shift and collisional damping in addition to Landau damping. Diagrams at the bottom show the relation to equilibrium approach: the use of the 2BA for Σcor is equivalent to solving the Bethe–Salpetr equation with a two-particle kernel Ξ containing the indicated diagrams. Figure modified from Ref. 154.

Close modal

A remarkable feature of real-time NEGF computations of the density response is that already fairly simple approximations for the self-energy Σ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 Bonitz203 and Bonitz et al.,280 

Ξ(1,3;4,2)=±δΣ(1;3)/δg(4;2).
(28)

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.

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), which leads to the harmonically perturbed Hamiltonian,167,253,281

ĤĤq,ω,A=Ĥ0+2Al=1Ncos(q·r̂lωt).
(29)

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,

δn(q,ω,A)=Aχ(q,ω)ϕq,ω,
(30)

where ϕq,ω 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 harmonics175 or mode-coupling between multiple perturbations177 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 χ(q,ω) 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 

S(q,ω)=Imχ(q,ω)πn(1eβω),
(31)

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,ω). In fact, S(q,ω) fully determines χ(q,ω), 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,ω) is given by the convolution of the DSF with the combined source and instrument function R(ω), see Eq. (3), and the wave vector q is determined by the scattering angle. Therefore, the availability of an accurate model for either χ(q,ω) [subsequently evaluating Eq. (31)] or directly for S(q,ω) 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 models147,148,151,186 that decompose the full DSF into contributions from bound and free electrons and their respective transitions directly work with S(q,ω). In practice, these models can easily be convolved with R(ω) 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

S(q)=dωS(q,ω).
(32)

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 

g(r)=1+1ndq(2π)3(S(q)1)eiq·r.
(33)

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

W=12Vdr1Vdr2n(2)(r1,r2)|r1r2|.
(34)

where n(2)(r1,r2)=n(r1)n(r2)g(r1r2) 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)n2(r1,r2)n(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 

fxc=01dλŴλλ.
(35)

In other words, integrating over the interaction energy of a system where the interaction contribution to the Hamiltonian has been re-scaled by λ[0,1],Ĥλ=Ĥ0+λŴ gives access to the free energy. From the λ 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

fxc(rs,θ)=1rs20rsrsW(rs,θ)drs,quantumOCP,
(36a)
fc(Γ)=0ΓW(Γ)ΓdΓ,classicalOCP.
(36b)

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,ω), which give access to χ(q,ω), 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 liquid288 with rs10, with notable progress reported by Tanaka200 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 scheme201 simplifies to the former200 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,ω)G(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 feature252,292 in the dilute electron gas.

Finally, we note that Eq. (35) has been used by a number of groups28,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 ω0. In the ground state, the stiffness theorem59 then gives a straightforward relation between the induced change in the total energy and the static density response function χ(q)χ(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

ρ̂kq,A=1Vl=1Neik·r̂lq,A,
(37)

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

ρ̂kq,A=δk,qχ(q)A.
(38)

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,θ), 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,τ) [Eq. (39)], which is discussed in Secs. II C.

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

F(q,τ)=n̂(q,0)n̂(q,τ).
(39)

From a theoretical perspective, the estimation of Eq. (39) within a PIMC simulation is straightforward190,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,τ) and the static density response function χ(q), which is known as the imaginary-time version of the fluctuation–dissipation theorem,190,293

χ(q)=n0βdτF(q,τ).
(40)

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,τ) and the DSF, which simply comprises a two-sided Laplace transform,

F(q,τ)=L[S(q,ω)]=dωeωτS(q,ω).
(41)

Traditionally, Eq. (41) constitutes the starting point for a so-called analytic continuation,308 i.e., a numerical inversion to obtain S(q,ω) 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,ω) 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,ω) 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,ω), 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 DSF59 

S(q,ω)=m,lPmnml(q)2δ(ωωlm).
(42)

In other words, S(q,ω) 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̂(q), with Pm being the probability to occupy the initial state m, nml(q) is the corresponding transition matrix element, and ωlm=(ElEm)/ denotes the energy difference. Inserting Eq. (42) into Eq. (41) gives the analogous representation in the τ-domain,

F(q,τ)=m,lPmnml(q)2eτωlm.
(43)

Evidently, the τ-decay of F(q,τ) for a given wave vector q is shaped by the characteristic frequencies in the corresponding S(q,ω). In particular, Eq. (43) directly implies that energetically low-lying excitations such as the roton feature in the dilute UEG292 will directly manifest as a less pronounced decay with τ. Conversely, it holds that ω(q)q2 in the non-collective single-particle regime with qqF (with qF being the usual Fermi wave number59), 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:

ΔFτ(q)=F(q,0)F(q,τ)F(q,0),
(44)

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,τ), 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,ω), 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,ω) of real materials in the WDM regime.

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

L[S(q,ω)]=L[S(q,ω)R(ω)]L[R(ω)]=L[I(q,ω)]L[R(ω)].
(45)

Indeed, the LHS of Eq. (45) is the ITCF F(q,τ), 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

S(q,ω)=S(q,ω)eβω.
(46)

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χ(q,ω). Inserting Eq. (46) into Eq. (41) gives the important symmetry relation,146,191

F(q,τ)=0dωS(q,ω){eωτ+eω(βτ)}=F(q,βτ).
(47)

In particular, Eq. (47) implies that F(q,τ) is symmetric around τ=β/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 as59,153,320

MαS=ωαS=dωS(q,ω)ωα.
(48)

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

M1S=ω1S=q22.
(49)

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

MαS=(1)ααταF(q,τ)|τ=0.
(50)

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,ω) based on static structural properties.323–326 

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 δn(r,t) expansion,

δn(1)=l=2Nw(l1)({l})i=2lV(i)d(i),
(51)

where 1={r1,t1},{l}=(1,2,,l),d1={dr1,dt1}. After writing the first three terms explicitly and setting χ(1,2)w(1)(1,2), we have171,181

δn(1)=d2χ(1,2)V(2)+d2d3w(2)(1,2,3)V(2)V(3)+d2d3d4w(3)(1,2,3,4)V(2)V(3)V(4).
(52)

Here, all time integrations run from to t1, 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)=iT{ψ+(1)ψ(1)} with respect to the external perturbation. They are special cases of more general higher-order correlation functions. For the linear response function, we have

χ(1,2)=L(12,12)|1=1+2=2=±iδg(11)δV(22)|1=1+2=2.
(53)

Here, t1+=t1+ε. 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

w(2)(1,2,3)=(±i)2δ2g(11+)δV(33)δV(22).
(54)

Again, this is a special case of the three-particle correlator Y(123,123). 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

χ(1,2)=d3Π(13)K(23),
(55)

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

w(2)(1,2,3)=d4d5Π(1,4,5)K(35)K(24)+d4Π(14)K(3,2,4).
(56)

This defines the quadratic polarization Π(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,

δn(k)=Aχ(k){δ(kq)+δ(k+q)}+A2{w(2)(kq,q)[δ(k2q)+δ(k)]+w(2)(k+q,q)[δ(k+2q)+δ(k)]}.
(57)

Here, χ(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

w(2)(q1,q2)=12V0βdτ10βdτ2ρ̃(q1+q2,0)ρ̃(q1,τ1)ρ̃(q2,τ2).
(58)

The density operator

ρ̃(q,τ)=l=1Nexp(iq·r̂l,τ)
(59)

is not normalized, and rl,τ denotes the position of particle l at an imaginary time τ[0,β]. 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,

ρ̂qq,A=χ(q)A+χ(1,3)(q)A3+,
(60)

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

ρ̂2qq,A=χ(2,2)(q)A2+χ(2,4)(q)A4+.
(61)

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

F(2)(q,τ1,τ2)=ρ̃(2q,0)ρ̃(q,τ1)ρ̃(q,τ2),
(62)

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

ρ̂3qq,A=χ(3,3)(q)A3+χ(3,5)(q)A5+.
(63)

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, χ(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

χ0(2,2)(q)=2q2[χ0(1,1)(2q)χ0(1,1)(q)].
(64)

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

χRPA(2,2)(q)=χ0(2,2)(q)[1v(q)χ0(1,1)(q)]2[1v(2q)χ0(1,1)(2q)].
(65)

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

χLFC(2,2)(q)=χ0(2,2)(q)[1v(q)[1G(q)]χ0(1,1)(q)]2×[1v(2q)[1G(2q)]χ0(1,1)(2q)]1.
(66)

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=N=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 conditions330–334 and in the WDM regime64,317,335 and are not covered in the present work.

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 results336–339 in the zero-temperature limit that were subsequently used as input for various parametrizations340–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 theories163,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,θ). 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 rs = 4, which is close to sodium.351,352 Panel (a) shows the raw PIMC results for the ITCF F(q,τ) [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 q2π/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τβ/2 due to the symmetry relation derived in Eq. (47). The F(q,τ) 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 τ0,F(q,0)=S(q); this is indicated by the bold-dashed green curve.

FIG. 4.

PIMC results for the ITCF and the static linear density response of the UEG for N = 34, rs = 4, and θ = 1. (a) ITCF in the qτ-plane in the interval 0τβ/2; the dashed green line indicates the limit of τ = 0, where the ITCF becomes equal to the static structure factor, F(q,0)=S(q), and the light blue line indicates the thermal structure factor F(q,β/2). (b) Full τ-dependence of F(q,τ) for a particular wave vector with q=1.88qF; solid green: PIMC data from Ref. 190; dotted red: single-particle decay model from Ref. 192; dash-dotted black: RPA, dashed blue: static approximation,252G(q,ω)G(q,0). (c) Static linear density response function χ(q), green crosses: ITCF-based evaluation of Eq. (40); red dots: estimation from PIMC simulations of the harmonically perturbed system [cf. Eq. (29)] via Eq. (60); solid blue: ideal Lindhard function, χ0(q); dash-dotted black: RPA; dotted yellow: long-wavelength expansion Eq. (67). (d) Symbols: raw PIMC results for the induced density, Eq. (37) for different q as a function of the perturbation amplitude A; dotted lines: corresponding fits according to Eq. (60); solid bars on the y-axis: LRT limit computed from the ITCF.

FIG. 4.

PIMC results for the ITCF and the static linear density response of the UEG for N = 34, rs = 4, and θ = 1. (a) ITCF in the qτ-plane in the interval 0τβ/2; the dashed green line indicates the limit of τ = 0, where the ITCF becomes equal to the static structure factor, F(q,0)=S(q), and the light blue line indicates the thermal structure factor F(q,β/2). (b) Full τ-dependence of F(q,τ) for a particular wave vector with q=1.88qF; solid green: PIMC data from Ref. 190; dotted red: single-particle decay model from Ref. 192; dash-dotted black: RPA, dashed blue: static approximation,252G(q,ω)G(q,0). (c) Static linear density response function χ(q), green crosses: ITCF-based evaluation of Eq. (40); red dots: estimation from PIMC simulations of the harmonically perturbed system [cf. Eq. (29)] via Eq. (60); solid blue: ideal Lindhard function, χ0(q); dash-dotted black: RPA; dotted yellow: long-wavelength expansion Eq. (67). (d) Symbols: raw PIMC results for the induced density, Eq. (37) for different q as a function of the perturbation amplitude A; dotted lines: corresponding fits according to Eq. (60); solid bars on the y-axis: LRT limit computed from the ITCF.

Close modal

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,β/2) (light blue dashed curve) exhibits a distinct structure with a pronounced peak around q2qF. We note that the physical behavior of F(q,τ) has often been described as featureless in the literature221,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,τ) 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,ω) without any analytic continuation,190,192 see also the discussion of Fig. 10. For example, the single-particle dispersion ω(q)q2 directly manifests in F(q,τ) 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,τ) 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,τ) 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,ω)0 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,ω)G(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,ω) of the UEG for metallic densities rs4. 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 q2qF.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|qF that is connected to the so-called on-top pair correlation function, g(r=0). The corresponding effective static approximation159 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 χ(q) from the PIMC data for F(q,τ) via Eq. (40). The results are shown in Fig. 4(c), where we show the, thus, obtained PIMC results for χ(q) as the green crosses. In the limit of small q, the exact density response is given by Kugler320 (dotted yellow curve),

limq0χ(q)=q24π;
(67)

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 q2qF. From a physical perspective, this can be understood intuitively from the following considerations: for intermediate q, the wavelength λ=2π/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 χ(q) is proportional to the area under F(q,τ) for a given value of q. Therefore, the peak in χ(q) is directly related to the maximum in F(q,β/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 qqF. 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π/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 χ(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 χ(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 χ(q) from the ITCF, Dornheim et al.194 have carried out extensive PIMC calculations for 50 density-temperature combinations; the corresponding parameters are depicted as the red crosses in the rs-θ-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).

FIG. 5.

Neural network representation of the static local field correction G(q;rs,θ). (a) Available training data in the rs-θ-plane: red crossed: PIMC simulations based on the ITCF via Eq. (40); black squares (MCS): ground-state QMC calculations by Moroni et al.281 based on simulations of the harmonically perturbed system, cf. Eq. (29); dashed blue (CDOP): ground-state parameterization based on the MCS data by Corradini et al.298 The shaded green area indicates the validity range of the neural network. (b) Static local field correction at θ = 1 in the rs-q-plane. The red crosses show PIMC results, and the green surface is the prediction by the neural network. Reprinted with permission from Dornheim et al.,194 J. Chem. Phys. 151, 194104 (2019). Copyright 2019 AIP Publishing.

FIG. 5.

Neural network representation of the static local field correction G(q;rs,θ). (a) Available training data in the rs-θ-plane: red crossed: PIMC simulations based on the ITCF via Eq. (40); black squares (MCS): ground-state QMC calculations by Moroni et al.281 based on simulations of the harmonically perturbed system, cf. Eq. (29); dashed blue (CDOP): ground-state parameterization based on the MCS data by Corradini et al.298 The shaded green area indicates the validity range of the neural network. (b) Static local field correction at θ = 1 in the rs-q-plane. The red crosses show PIMC results, and the green surface is the prediction by the neural network. Reprinted with permission from Dornheim et al.,194 J. Chem. Phys. 151, 194104 (2019). Copyright 2019 AIP Publishing.

Close modal

On the one hand, the inversion of Eq. (4) to compute G(q) from PIMC data for χ(q) is straightforward. On the other hand, G(q;rs,θ) exhibits a nontrivial dependence on q based on different combination of rs and θ. This is illustrated in Fig. 5(b), where we show the static local field correction in the rs-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 rs), G(q) exhibits a pronounced increase in the limit of large q. This increase becomes less pronounced with decreasing rs 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,356,357 which is negative for some parameters at finite temperatures.306,358 To construct a reliable representation of G(q;rs,θ) 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 model359 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,θ) 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 rs = 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,θ).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 

FIG. 6.

Static density response of hydrogen at the electronic Fermi temperature, θ = 1, for (a) rs = 2 and (b) rs = 4. Green circles: PIMC reference data by Böhme et al.;77 black stars: DFT results using the PBE XC functional by Moldabekov et al.;255 dotted red: density response of the UEG at the same conditions evaluated using the neural-net representation from Ref. 194; dash-dotted yellow [only panel (b)]: density response of a UEG with modified rs* and θ* corresponding to a free-electron fraction of α=0.6.

FIG. 6.

Static density response of hydrogen at the electronic Fermi temperature, θ = 1, for (a) rs = 2 and (b) rs = 4. Green circles: PIMC reference data by Böhme et al.;77 black stars: DFT results using the PBE XC functional by Moldabekov et al.;255 dotted red: density response of the UEG at the same conditions evaluated using the neural-net representation from Ref. 194; dash-dotted yellow [only panel (b)]: density response of a UEG with modified rs* and θ* corresponding to a free-electron fraction of α=0.6.

Close modal

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 rs = 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 rs.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 electrons148,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 θ* that correspond to an effective free-electronic fraction of α=0.6, i.e., 60%; this is consistent with the value of α=0.54 reported by Militzer and Ceperley.360 While this effective density response overall exhibits the correct magnitude around q2qF, 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 q4qF, which has been attributed to an isotropy breaking due to the presence of the ions at small length scales λq=2π/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 functional364 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 rs = 2, and the DFT results are in good agreement with the PIMC reference data. For rs = 4, the agreement is less good, and the DFT calculations overestimate the actual response over the entire depicted 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 α=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 χ(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 χS(q,ω). 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 rs = 4 and θ = 1. As a reference function to obtain the kernel from the physical density response function χ(q), Moldabekov et al.255 have used the same mean-field response function χ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 χ(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,

limq0G(q)=q24π2n2(nFxc),
(68)

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,θ) 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 χ(q) in Eq. (17), which corresponds to setting Kxc(q)0, 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.

FIG. 7.

Static XC kernel, plotted as the local field correction G(q)=Kxc(q)/v(q), of hydrogen at θ = 1 and rs = 4. The solid black and dashed red line show the full G(q) of the UEG at the same conditions evaluated from the neural-net representation from Ref. 194 and the corresponding exact long-wavelength expansion given by the compressibility sum rule (CSR) Eq. (68). The black squares have been computed from the exact PIMC reference data by Böhme et al.,77 and the blue circles, orange triangles, and green diamonds from DFT results for χ(q) (see Fig. 6) using the LDA,341 PBE,364 and SCAN367 XC functionals, respectively. Adapted from Moldabekov et al.255 

FIG. 7.

Static XC kernel, plotted as the local field correction G(q)=Kxc(q)/v(q), of hydrogen at θ = 1 and rs = 4. The solid black and dashed red line show the full G(q) of the UEG at the same conditions evaluated from the neural-net representation from Ref. 194 and the corresponding exact long-wavelength expansion given by the compressibility sum rule (CSR) Eq. (68). The black squares have been computed from the exact PIMC reference data by Böhme et al.,77 and the blue circles, orange triangles, and green diamonds from DFT results for χ(q) (see Fig. 6) using the LDA,341 PBE,364 and SCAN367 XC functionals, respectively. Adapted from Moldabekov et al.255 

Close modal

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 Zunger341) PBE364 (orange triangles) and SCAN367 (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 q2qF, 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.

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 ω0. While being an important step in the right direction, many important applications such as the modeling of XRTS experiments147,164,166,189 or the construction of advanced XC functionals within the adiabatic-connection fluctuation–dissipation formulation of DFT99,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) method276,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 decomposition189 and the possibility to estimate both nonequilibrium377 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 potential383,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 experiments106,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 {ϕα}. The latter are then used to compute the dynamic KS-response function χS(q,ω) via Eq. (16); the key ingredient is then given by the XC kernel Kxc(q,ω), which relates χS(q,ω) to the actual physical dynamic linear density response function χ(q,ω) 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 χS(q,ω). 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,ω)Kxc(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 mechanics228 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 system190,193,195 in thermodynamic equilibrium; exact PIMC results for the ITCF F(q,τ) of the warm dense UEG are shown in Fig. 4. Indeed, the ITCF, in principle, contains the same information as the DSF S(q,ω)—albeit in an unfamiliar representation190,192—and is directly useful for the interpretation of XRTS experiments,146 see Sec. IV. Yet, the analytic continuation308 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 UEG252,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 collaborators396,397 have presented accurate data for S(q,ω) 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, θ=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,τ). 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,ω). More specifically, trial solutions for G(q,ω) give straightforward access to χ(q,ω) [Eq. (4)], and via the fluctuation-dissipation theorem [Eq. (31)] also to the corresponding DSF S(q,ω); finally, the latter is inserted into Eq. (41) and compared to the PIMC reference data for F(q,τ). 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,ω),161,252,287,306 which sufficiently constrains the space of possible trial solutions for S(q,ω).

In Fig. 8, we show corresponding results for both the real (top panel) and the imaginary part (bottom panel) of G(q,ω) for the UEG at θ = 1, rs = 6318 with the shaded gray area indicating the given uncertainty interval. For the real part, the static ω0 limit can be inferred from χ(q), while the high frequency ω 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 ω2ωp, with ω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 ω0 and ω and exhibits a single broad peak around ω5ωp.

FIG. 8.

Ab initio PIMC results for the frequency-dependence of the dynamic local field correction G(q,ω) of the UEG with θ = 1 and rs = 6 for q=1.88qF. The shaded gray areas indicate the respective uncertainty intervals of the real [top] and imaginary [bottom] parts. Reprinted with permission from Hamann et al.,318 Phys. Rev. B 102, 125150 (2020). Copyright 2020 American Physical Society.

FIG. 8.

Ab initio PIMC results for the frequency-dependence of the dynamic local field correction G(q,ω) of the UEG with θ = 1 and rs = 6 for q=1.88qF. The shaded gray areas indicate the respective uncertainty intervals of the real [top] and imaginary [bottom] parts. Reprinted with permission from Hamann et al.,318 Phys. Rev. B 102, 125150 (2020). Copyright 2020 American Physical Society.

Close modal

From a theoretical perspective, accurate knowledge of G(q,ω) gives one direct access to the dynamic density response function χ(q,ω) [cf. Eq. (4)] and, in this way, a gamut of other properties such as the dynamic dielectric function ϵ(q,ω) or the electric conductivity σ(q,ω) 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 experiments147,149 [cf. Eq. (3)], one of the most interesting properties in this regard is given by S(q,ω), 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,ω),ω(q), compared to the RPA for metallic densities. This is shown in Fig. 9(a) for rs = 4 at the electronic Fermi temperature, θ = 1. Specifically, the solid green curve depicts the exact ω(q) based on the full PIMC results for the dynamic local field correction G(q,ω), 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,ω)G(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 χ0(q,ω)] 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

FIG. 9.

Position of the maximum of the DSF as a function of the wave number, ω(q), of the UEG at the electronic Fermi temperature, θ = 1, at (a) rs = 4 and (b) rs = 10. Solid green: exact PIMC results taken from Ref. 252; dashed blue: static approximation, i.e., setting G(q,ω)G(q,0) in Eq. (4); dash-dotted black: RPA; dotted red: pair alignment model introduced in Ref. 292. The vertical red arrows indicate the XC induced down-shift ΔWxc of the true ω(q) compared to the mean-field description within RPA.

FIG. 9.

Position of the maximum of the DSF as a function of the wave number, ω(q), of the UEG at the electronic Fermi temperature, θ = 1, at (a) rs = 4 and (b) rs = 10. Solid green: exact PIMC results taken from Ref. 252; dashed blue: static approximation, i.e., setting G(q,ω)G(q,0) in Eq. (4); dash-dotted black: RPA; dotted red: pair alignment model introduced in Ref. 292. The vertical red arrows indicate the XC induced down-shift ΔWxc of the true ω(q) compared to the mean-field description within RPA.

Close modal

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 rs = 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 ω(q) exhibit a nonmonotonous behavior, with a minimum at intermediate wave numbers, q2qF. 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 law413) but also real liquids (water,414,415 elemental noble liquids,416 and liquid alkali metals417), 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ω(q)/dq<0 as q0, at intermediate coupling Γ10.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 ansatz426 [S(q) connection] for quantum liquids and through the quasi-localized charge approximation427 [g(r) connection] or the static local field corrected dielectric formalism428 [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 ω(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,ω) 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 ω(q) around q2qF, the wavelength of the cosinuoidal perturbation is comparable to the average interparticle distance, λq=2π/qd. 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 ω(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 ΔWxc in Fig. 9; the corresponding dotted red curves have been computed as ωPA(q)=ωRPA(q)+ΔWxc(q). For the metallic density of rs = 4, the PA model nicely captures the true down-shift of ω(q) and, therefore, is directly relevant for the modeling of XRTS experiments of WDM.157 For the strongly coupled case of rs = 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 ω with the shoulder at ωRPA(q), see Ref. 292 for a more detailed explanation.

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

M1S=χ(q)2n.
(69)

Here, M1S 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,ω) to smaller frequencies ω is associated with an increase in the magnitude of χ(q). This explains the emergence of an increasingly sharp peak in χ(q) around q=2qF, which has been reported in previous studies of the electron liquid both at finite temperature201,288 and in the ground state.431 

While having exact reference results for S(q,ω) gives one valuable insights into a number of physical effects, the required analytic continuation308 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,τ) 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,τ) is symmetric around τ=β/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,τ)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 ω(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,τ) 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 ΔFβ/2(q,τ) defined in Eq. (44).

The results are shown in Fig. 10 for θ = 1 with green and red datasets corresponding to rs = 4 and rs = 10, respectively. More specifically, the symbols show the PIMC reference data for ω(q) that are also depicted in Fig. 9, and the corresponding solid lines show the relative decay measure ΔFβ/2(q,τ) 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 q0, which is given by the usual plasma frequency ωp in the case of ω(q). In addition, all curves exhibit a characteristic single-particle limit for qqF. For the DSF, it is given by the well-known single-particle dispersion ω(q)q2. For the ITCF, one finds ΔFβ/2(q,τ)=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,τ) with 0τβ/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 τ-dependence of the ITCF.

FIG. 10.

Dispersion of dynamic density fluctuations in the UEG at θ = 1 for rs = 4 (green) and rs = 10 (red). Symbols: PIMC-based results for the position of the maximum in the DSF ω(q), cf. Fig. 9; the corresponding curves depict the relative decay measure of the ITCF ΔFβ/2(q) defined in Eq. (44). All datasets have been re-normalized with respect to the plasmon limit for q0. Data taken from Refs. 192 and 252.

FIG. 10.

Dispersion of dynamic density fluctuations in the UEG at θ = 1 for rs = 4 (green) and rs = 10 (red). Symbols: PIMC-based results for the position of the maximum in the DSF ω(q), cf. Fig. 9; the corresponding curves depict the relative decay measure of the ITCF ΔFβ/2(q) defined in Eq. (44). All datasets have been re-normalized with respect to the plasmon limit for q0. Data taken from Refs. 192 and 252.

Close modal

For intermediate q, we find a close correspondence between the respective ω(q) and ΔFβ/2(q,τ) for the same value of the density parameter rs. Indeed, the relative decay measure even becomes slightly negative for rs = 10. We note that the phenomenological similarity between the roton minimum in the DSF and the reduced τ-decay of the ITCF becomes even more apparent at rs = 20, where the alignment of pairs of electrons292 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 liquids301,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 decomposition148,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,ω)Kxc(q) for the UEG for metallic densities rs4, 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) rs = 2 and (b) rs = 4 at the electronic Fermi temperature, θ = 1. More specifically, all curves have been computed via Eq. (17) using the KS-response function χS(q,ω) based on a set of KS orbitals from an equilibrium DFT simulation with the LDA functional.

FIG. 11.

LR-TDDFT results for hydrogen at the electronic Fermi temperature θ = 1 for (a) rs = 2 and (b) rs = 4 for a wave vector of q=1.69qF. Solid black: RPA, i.e., setting Kxc(q,ω)0 in Eq. (17); dashed green: adiabatic LDA (ALDA), cf. Eq. (68); dotted red: static approximation, Kxc(q,ω)Kxc(q) using exact PIMC data for the XC kernel in the static limit. Adapted from Ref. 77.

FIG. 11.

LR-TDDFT results for hydrogen at the electronic Fermi temperature θ = 1 for (a) rs = 2 and (b) rs = 4 for a wave vector of q=1.69qF. Solid black: RPA, i.e., setting Kxc(q,ω)0 in Eq. (17); dashed green: adiabatic LDA (ALDA), cf. Eq. (68); dotted red: static approximation, Kxc(q,ω)Kxc(q) using exact PIMC data for the XC kernel in the static limit. Adapted from Ref. 77.

Close modal

Let us start our analysis for rs = 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,ω)0 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 q0 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 rs = 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 χS(q,ω) 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 χS(q,ω) 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 Ladder432 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,τ) 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 q0). Figure 12 shows the frequency-dependent electrical conductivity23 of iron at earth-core conditions (P320 GPa or rs2.2,T0.55 eV). The red curve shows the tensor averaged result of several ionic configurations of size N = 16 with an energy resolution Δω 0.11 eV, which is proportional to the inverse of the total propagation time (t1500 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 initio433,434,437 methods.

FIG. 12.

Frequency-dependent electrical conductivity of iron under Earth-core conditions (T0.55 eV, rs2.2) from RT-TDDFT calculations (red curve). This is compared with previous works (black) using the KG formula.433,434 The Drude fit to RT-TDDFT data is indicated by the red dotted line.

FIG. 12.

Frequency-dependent electrical conductivity of iron under Earth-core conditions (T0.55 eV, rs2.2) from RT-TDDFT calculations (red curve). This is compared with previous works (black) using the KG formula.433,434 The Drude fit to RT-TDDFT data is indicated by the red dotted line.

Close modal

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 q0 limit of the optical conductivity) and the plasmon (the Mermin approximation, with collision frequencies taken from Kubo–Greenwood calculations442), but both are based on thermal DFT in a limit that is equivalent to an LR-TDDFT calculation in which fHxc = 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 pseudization444 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.

FIG. 13.

(Top) The RT-TDDFT density response to a perturbation of the form v0exp(iq·r)f(t), where q=0.67 a.u., and f(t) is a Gaussian envelope. Results for adiabatic PBE and SCAN are included, and the envelope for a postprocessing window with a 3-eV Lorentzian broadening is indicated in the background. (Bottom) The RT-TDDFT results postprocessed into simulated XRTS spectra through convolution with a 4-Gaussian least squares fit to the instrument function in Ref. 145. We simulate a range of ionic structures, spanned by the shaded regions, by subtracting the elastic peak from the RT-TDDFT and adding a range of model ionic peaks. We compare to digitizations of both the experimental spectra and a subset of the theoretical data by Witte et al.145 The theoretical data from Witte et al.145 are the bound-free edge from the Kubo–Greenwood optical conductivity and the rest of the spectrum from the Mermin approximation with collision frequencies determined using Kubo–Greenwood—both of which used HSE. Note that the SCAN results are on top of the PBE results, with the exception of a slight difference in the bound-free edge.

FIG. 13.

(Top) The RT-TDDFT density response to a perturbation of the form v0exp(iq·r)f(t), where q=0.67 a.u., and f(t) is a Gaussian envelope. Results for adiabatic PBE and SCAN are included, and the envelope for a postprocessing window with a 3-eV Lorentzian broadening is indicated in the background. (Bottom) The RT-TDDFT results postprocessed into simulated XRTS spectra through convolution with a 4-Gaussian least squares fit to the instrument function in Ref. 145. We simulate a range of ionic structures, spanned by the shaded regions, by subtracting the elastic peak from the RT-TDDFT and adding a range of model ionic peaks. We compare to digitizations of both the experimental spectra and a subset of the theoretical data by Witte et al.145 The theoretical data from Witte et al.145 are the bound-free edge from the Kubo–Greenwood optical conductivity and the rest of the spectrum from the Mermin approximation with collision frequencies determined using Kubo–Greenwood—both of which used HSE. Note that the SCAN results are on top of the PBE results, with the exception of a slight difference in the bound-free edge.

Close modal

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 q0, 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.

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 calculations181–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 Z3 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 shown176 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 shown167,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 ρ̂kq,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 A0.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 rs, θ, 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;τ1,τ2) is given by Eq. (58).

In Fig. 14, we show new ab initio PIMC results for F(2)(q;τ1,τ2) in the τ1-τ2-plane for two different values of the wave number q for the UEG at rs = 4 and θ = 1, i.e., for the same conditions as in Fig. 4. Overall, we find that F(2)(q;τ1,τ2) exhibits a qualitatively similar behavior to the two-body ITCF F(q,τ). 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;τ1,τ2) remains very poorly understood, and its improved understanding along similar lines as it has been achieved for F(q,τ) in Refs.190,192 remains an important task for future works.

FIG. 14.

Ab initio PIMC results for the three-body ITCF F(2)(q;τ1,τ2) [Eq. (62)] for N = 34, rs = 4, and θ = 1 in the τ1-τ2-plane; (a) q=1.25qF and (b) q=2.51qF.

FIG. 14.

Ab initio PIMC results for the three-body ITCF F(2)(q;τ1,τ2) [Eq. (62)] for N = 34, rs = 4, and θ = 1 in the τ1-τ2-plane; (a) q=1.25qF and (b) q=2.51qF.

Close modal

In the present context, the main utility of F(2)(q;τ1,τ2) is given by its direct connection to the quadratic static density response function of the second harmonic, χ(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, χ(2,2)(q) exhibits the opposite sign compared to the linear static density response function χ(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., ρ̂2qq,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.

FIG. 15.

Quadratic density response at the second harmonic χ(2,2)(q,0) for N = 34, rs = 4, and θ = 1. Red circles: estimation from PIMC simulations of the harmonically perturbed system [cf. Eq. (29)] via Eq. (61); green crosses: integration of the three-body ITCF F(2)(q,τ1,τ2) shown in Fig. 14 via Eq. (58); solid yellow: quadratic density response of the ideal Fermi gas, Eq. (64); dashed black: RPA expression, Eq. (65); dotted blue: incorporation of XC effects via the static local field correction G(q) taken from Ref. 194 via Eq. (66).

FIG. 15.

Quadratic density response at the second harmonic χ(2,2)(q,0) for N = 34, rs = 4, and θ = 1. Red circles: estimation from PIMC simulations of the harmonically perturbed system [cf. Eq. (29)] via Eq. (61); green crosses: integration of the three-body ITCF F(2)(q,τ1,τ2) shown in Fig. 14 via Eq. (58); solid yellow: quadratic density response of the ideal Fermi gas, Eq. (64); dashed black: RPA expression, Eq. (65); dotted blue: incorporation of XC effects via the static local field correction G(q) taken from Ref. 194 via Eq. (66).

Close modal

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 q2qF. 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 χ(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 q0 and q but substantially underestimate the true quadratic response around its maximum at q2qF. Including the static local field correction G(q;rs,θ) 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 χ(2,2)(q) compared to χ(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,θ) based on different representations.72,194 Moreover, these efforts can be improved by including the full, frequency-dependent results for G(q,ω) that are available for certain parameters of the UEG based on the analytic continuation252 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 plasmons250,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 hydrogen77,78 will allow one to directly estimate F(2)(q;τ1,τ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,τ), which have been touched upon throughout the present work.

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,ω) 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(ω), cf. Eq. (3), then allows one to compare these results with the experimentally measured scattering intensity I(q,ω), 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 methods114 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 measurements187,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,ω) to the imaginary time τ; here, the main property is the ITCF F(q,τ), 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,τ), 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 τ=β/2=1/2T. In this sense, one could say that XRTS measurements effectively function as a thermometer as the relation between F(q,τ) 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 ω±25 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,

Lx[S(q,ω)R(ω)]=xxdωeτω{S(q,ω)R(ω)},
(70)

with the corresponding truncated ITCF

Fx(q,τ)=Lx[S(q,ω)R(ω)]L[R(ω)].
(71)
FIG. 16.

Temperature diagnosis of the XRTS measurement of warm dense beryllium by Glenzer et al.202 (a) XRTS signal (green) and combined source and instrument function (blue); (b) convergence of the inferred temperature with the integration boundary x [Eq. (71)] and the corresponding uncertainty interval due to the noise in the intensity I(q,ω); (c) extracted ITCF for x=40 eV (solid green); the shaded gray area illustrates the uncertainty in F(q,τ), and the dotted blue line shows the Laplace transform of the instrument function; the dotted red curve has been mirrored around τ=β/2, which illustrates that the symmetry relation Eq. (47) is fulfilled to a high degree. Taken from Ref. 191.

FIG. 16.

Temperature diagnosis of the XRTS measurement of warm dense beryllium by Glenzer et al.202 (a) XRTS signal (green) and combined source and instrument function (blue); (b) convergence of the inferred temperature with the integration boundary x [Eq. (71)] and the corresponding uncertainty interval due to the noise in the intensity I(q,ω); (c) extracted ITCF for x=40 eV (solid green); the shaded gray area illustrates the uncertainty in F(q,τ), and the dotted blue line shows the Laplace transform of the instrument function; the dotted red curve has been mirrored around τ=β/2, which illustrates that the symmetry relation Eq. (47) is fulfilled to a high degree. Taken from Ref. 191.

Close modal

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

limxFx(q,τ)=F(q,τ),
(72)

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,τ) 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 ω30 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±2 eV. This is comparably close to the nominal temperature of T=12 eV, which was found in the original publication202 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,τ), 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(ω), 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 τ=β/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(ω). This is relatively straightforward using source monitors at modern XFEL facilities,18,55 whereas the characterization of backlighter emission spectra462 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(ω) on the extracted temperature less pronounced in any case. Furthermore, it has been shown191 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 ω<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(ω) will have increasing consequences for the evaluation of Eq. (45) and, thus, make the localization of the minimum around τ=β/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 effects270,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(ω), 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.

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,ω) and the ITCF F(q,τ). 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,τ) from XRTS experiments is straightforward; this allows for exact and model-free diagnostics of the temperature146 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 χ(q) and S(q,ω) 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.

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 scaling372,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.

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 Ĥ0 and with respect to a modified Hamiltonian Ĥ=Ĥ0+V̂ext that is perturbed by a monochromatic static potential V̂ext, 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 χ(q) and, in this way, to the static XC kernel Kxc(q) of real materials (top right).

FIG. 17.

Schematic illustration of the DFT setup to the static XC kernel by Moldabekov et al.255 Top left: standard KS-DFT is used to compute the single-electron density n(r) a) with respect to the original electronic Hamiltonian of interest Ĥ0 (solid black) and (b) with respect to a modified Hamiltonian subject to an external monochromatic perturbation V̂ext (dashed blue). Bottom left: this gives access to the corresponding density modulation Δn(r) for different wave vectors q. Top right: in the linear-response regime (V̂extĤ0), the induced density gives access to the static density response function χ(q) and the corresponding XC kernel Kxc(q) (blue crosses). Bottom right: the FDT [Eq. (31)] provides a direct connection to electron–electron correlation functions such as the static structure factor See(q), which cannot be readily computed within standard DFT.

FIG. 17.

Schematic illustration of the DFT setup to the static XC kernel by Moldabekov et al.255 Top left: standard KS-DFT is used to compute the single-electron density n(r) a) with respect to the original electronic Hamiltonian of interest Ĥ0 (solid black) and (b) with respect to a modified Hamiltonian subject to an external monochromatic perturbation V̂ext (dashed blue). Bottom left: this gives access to the corresponding density modulation Δn(r) for different wave vectors q. Top right: in the linear-response regime (V̂extĤ0), the induced density gives access to the static density response function χ(q) and the corresponding XC kernel Kxc(q) (blue crosses). Bottom right: the FDT [Eq. (31)] provides a direct connection to electron–electron correlation functions such as the static structure factor See(q), which cannot be readily computed within standard DFT.

Close modal

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,ω) 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,τ), 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,τ), 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,ω). 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 functions176 by using KS-DFT to study the nonlinear density response of a given system of interest.455 

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 XFEL55 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,ω) 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,τ) and the static density response function χ(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.

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).

The authors have no conflicts to disclose.

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).

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

1.
V. E.
Fortov
, “
Extreme states of matter on Earth and in space
,”
Phys.-Usp.
52
,
615
647
(
2009
).
2.
R. P.
Drake
,
High-Energy-Density Physics: Foundation of Inertial Fusion and Experimental Astrophysics, Graduate Texts in Physics
(
Springer International Publishing
,
2018
).
3.
A.
Benuzzi-Mounaix
,
S.
Mazevet
,
A.
Ravasio
,
T.
Vinci
,
A.
Denoeud
,
M.
Koenig
,
N.
Amadou
,
E.
Brambrink
,
F.
Festa
,
A.
Levy
,
M.
Harmand
,
S.
Brygoo
,
G.
Huser
,
V.
Recoules
,
J.
Bouchet
,
G.
Morard
,
F.
Guyot
,
T.
de Resseguier
,
K.
Myanishi
,
N.
Ozaki
,
F.
Dorchies
,
J.
Gaudin
,
P. M.
Leguay
,
O.
Peyrusse
,
O.
Henry
,
D.
Raffestin
,
S. L.
Pape
,
R.
Smith
, and
R.
Musella
, “
Progress in warm dense matter study with applications to planetology
,”
Phys. Scr.
2014
,
014060
.
4.
R. F.
Smith
,
J. H.
Eggert
,
R.
Jeanloz
,
T. S.
Duffy
,
D. G.
Braun
,
J. R.
Patterson
,
R. E.
Rudd
,
J.
Biener
,
A. E.
Lazicki
,
A. V.
Hamza
,
J.
Wang
,
T.
Braun
,
L. X.
Benedict
,
P. M.
Celliers
, and
G. W.
Collins
, “
Ramp compression of diamond to five terapascals
,”
Nature
511
,
330
333
(
2014
).
5.
J.
Vorberger
,
I.
Tamblyn
,
B.
Militzer
, and
S. A.
Bonev
, “
Hydrogen-helium mixtures in the interiors of giant planets
,”
Phys. Rev. B
75
,
024206
(
2007
).
6.
D.
Saumon
,
W. B.
Hubbard
,
G.
Chabrier
, and
H. M.
van Horn
, “
The role of the molecular-metallic transition of hydrogen in the evolution of Jupiter, Saturn, and Brown dwarfs
,”
Astrophys. J.
391
,
827
831
(
1992
).
7.
B.
Militzer
,
W. B.
Hubbard
,
J.
Vorberger
,
I.
Tamblyn
, and
S. A.
Bonev
, “
A massive core in Jupiter predicted from first-principles simulations
,”
Astrophys. J.
688
,
L45
L48
(
2008
).
8.
T.
Guillot
,
Y.
Miguel
,
B.
Militzer
,
W. B.
Hubbard
,
Y.
Kaspi
,
E.
Galanti
,
H.
Cao
,
R.
Helled
,
S. M.
Wahl
,
L.
Iess
,
W. M.
Folkner
,
D. J.
Stevenson
,
J. I.
Lunine
,
D. R.
Reese
,
A.
Biekman
,
M.
Parisi
,
D.
Durante
,
J. E. P.
Connerney
,
S. M.
Levin
, and
S. J.
Bolton
, “
A suppression of differential rotation in Jupiter's deep interior
,”
Nature
555
,
227
230
(
2018
).
9.
S.
Brygoo
,
P.
Loubeyre
,
M.
Millot
,
J. R.
Rygg
,
P. M.
Celliers
,
J. H.
Eggert
,
R.
Jeanloz
, and
G. W.
Collins
, “
Evidence of hydrogen-helium immiscibility at Jupiter-interior conditions
,”
Nature
593
,
517
521
(
2021
).
10.
J. J.
Fortney
,
M. S.
Marley
, and
J. W.
Barnes
, “
Planetary radii across five orders of magnitude in mass and stellar insolation: Application to transits
,”
Astrophys. J.
659
,
1661
(
2007
).
11.
Frontiers and Challenges in Warm Dense Matter
, edited by
F.
Graziani
,
M. P.
Desjarlais
,
R.
Redmer
, and
S. B.
Trickey
(
Springer, International Publishing
,
2014
).
12.
A.
Becker
,
W.
Lorenzen
,
J. J.
Fortney
,
N.
Nettelmann
,
M.
Schöttler
, and
R.
Redmer
, “
Ab initio equations of state for hydrogen (H-REOS.3) and helium (HE-REOS.3) and their implications for the interior of Brown Dwarfs
,”
Astrophys. J., Suppl. Ser.
215
,
21
(
2014
).
13.
S. V.
Lebedev
,
High Energy Density Laboratory Astrophysics, Astrophysics and Space Science
(
Springer
,
Netherlands
,
2007
).
14.
A. L.
Kritcher
,
D. C.
Swift
,
T.
Döppner
,
B.
Bachmann
,
L. X.
Benedict
,
G. W.
Collins
,
J. L.
DuBois
,
F.
Elsner
,
G.
Fontaine
,
J. A.
Gaffney
,
S.
Hamel
,
A.
Lazicki
,
W. R.
Johnson
,
N.
Kostinski
,
D.
Kraus
,
M. J.
MacDonald
,
B.
Maddox
,
M. E.
Martin
,
P.
Neumayer
,
A.
Nikroo
,
J.
Nilsen
,
B. A.
Remington
,
D.
Saumon
,
P. A.
Sterne
,
W.
Sweet
,
A. A.
Correa
,
H. D.
Whitley
,
R. W.
Falcone
, and
S. H.
Glenzer
, “
A measurement of the equation of state of carbon envelopes of white dwarfs
,”
Nature
584
,
51
54
(
2020
).
15.
E. H.
Gudmundsson
,
C. J.
Pethick
, and
R. I.
Epstein
, “
Structure of neutron star envelopes
,”
Astrophys. J.
272
,
286
300
(
1983
).
16.
P.
Haensel
,
A. Y.
Potekhin
, and
D. G.
Yakovlev
, eds., “
Equilibrium plasma properties. Outer envelopes
,” in
Neutron Stars 1
(
Springer
,
New York, New York
,
2007
), pp.
53
114
.
17.
D.
Kraus
,
A.
Ravasio
,
M.
Gauthier
,
D. O.
Gericke
,
J.
Vorberger
,
S.
Frydrych
,
J.
Helfrich
,
L. B.
Fletcher
,
G.
Schaumann
,
B.
Nagler
,
B.
Barbrel
,
B.
Bachmann
,
E. J.
Gamboa
,
S.
Göde
,
E.
Granados
,
G.
Gregori
,
H. J.
Lee
,
P.
Neumayer
,
W.
Schumaker
,
T.
Döppner
,
R. W.
Falcone
,
S. H.
Glenzer
, and
M.
Roth
, “
Nanosecond formation of diamond and lonsdaleite by shock compression of graphite
,”
Nat. Commun.
7
,
10970
(
2016
).
18.
S. H.
Glenzer
,
L. B.
Fletcher
,
E.
Galtier
,
B.
Nagler
,
R.
Alonso-Mori
,
B.
Barbrel
,
S. B.
Brown
,
D. A.
Chapman
,
Z.
Chen
,
C. B.
Curry
,
F.
Fiuza
,
E.
Gamboa
,
M.
Gauthier
,
D. O.
Gericke
,
A.
Gleason
,
S.
Goede
,
E.
Granados
,
P.
Heimann
,
J.
Kim
,
D.
Kraus
,
M. J.
MacDonald
,
A. J.
Mackinnon
,
R.
Mishra
,
A.
Ravasio
,
C.
Roedel
,
P.
Sperling
,
W.
Schumaker
,
Y. Y.
Tsui
,
J.
Vorberger
,
U.
Zastrau
,
A.
Fry
,
W. E.
White
,
J. B.
Hasting
, and
H. J.
Lee
, “
Matter under extreme conditions experiments at the Linac Coherent Light Source
,”
J. Phys. B
49
,
092001
(
2016
).
19.
K.
Ohta
,
Y.
Kuwayama
,
K.
Hirose
,
K.
Shimizu
, and
Y.
Ohishi
, “
Experimental determination of the electrical resistivity of iron at Earth's core conditions
,”
Nature
534
,
95
98
(
2016
).
20.
Z.
Konôpková
,
R. S.
McWilliams
,