Electronic density response of warm dense matter

Warm dense matter (WDM) is an extreme state that is characterized by high temperatures (⁠$T∈103−108$ K), high pressures (⁠$P∼1−104$ Gbar), and densities in the vicinity of and even partially exceeding solid state densities (⁠$n∼1021−1028$ 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.²¹ Experimental measurements and modeling are the subjects of very active investigations.^22–25 

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

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.⁴¹ Indeed, many ground-breaking results have been reported over the last few years, including the formation of nanodiamonds at extreme pressure^17,30,42 and the study of the liquid-liquid phase transition in hot dense hydrogen.^43–47 In addition, it can be expected that emerging capabilities at facilities such as NIF,^40,48 LCLS at SLAC,^18,49 the Sandia Z Pulsed Power Facility^50–52 and the OMEGA laser⁵³ in the USA, SACLA⁵⁴ in Japan, and the European XFEL⁵⁵ 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$⁠,⁵⁹ 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;⁵⁸ 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)⁶⁰ methods are capable of providing exact results for the static properties of the uniform electron gas (UEG)^61–75—the quantum version of the classical one-component plasma (OCP)^28,59,76—and light elements such as hydrogen^77,78 over substantial parts of the WDM regime. Moreover, path-integral Monte Carlo (PIMC) simulations within the fixed-node approximation⁷⁹—also known as restricted PIMC (RPIMC) in the literature—have allowed Militzer and co-workers^80–82 to go to heavier elements and material mixtures.⁸³ These simulations constitute the basis for an extensive equation-of-state (EOS) database⁸⁴ that can be used for a host of practical applications.

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

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

Many of the properties of WDM are encoded in its response to an external perturbation. Of particular interest are perturbations that couple to the electronic charge density, e.g., an incident electromagnetic field or the Coulomb field of an incident charged particle. The response to these probes contains information about a wide range of excitations, and many of which are sufficiently sensitive to the thermodynamic state of the WDM that they can be used to infer conditions like temperature¹⁴⁶ 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,¹⁴⁷ 

which is defined as

for homogeneous systems in equilibrium. For sufficiently weak probes in the linear response regime, the DSF is related to the linear density response function $χ(q,ω)$ through the fluctuation-dissipation theorem⁵⁹ (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(ω)$⁠,

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,ω)$ as¹⁵³ 

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.¹⁵⁴ 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.¹⁶⁷ have presented the first ab initio PIMC simulation results for the nonlinear density response of the warm dense UEG. Such nonlinear effects^168–177 become important, for example, for XRTS experiments with ultra-high intensities that can be realized at x-ray free-electron laser (XFEL) facilities via the novel seeding technique.¹⁷⁸ Moreover, they are an important ingredient to the construction of effective pair potentials at small distances^139,179,180 and stopping power calculations for slow and/or heavy particles.^181–183 Finally, it has been shown that nonlinear effects depend much more sensitively on system parameters like the temperature T,^167,175 which makes their experimental observation a potentially interesting new method of diagnostics.¹⁸⁴ 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),¹⁸⁵ which will be introduced later.

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

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.²⁰³ 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.

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=r1−r2$ and $t=t1−t2$⁠, one usually studies $χ(q,ω)$⁠.

In the limit of small A, linear response theory (LRT)⁵⁹ 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 models^204–209 as well as integral equation theory approaches within the hypernetted-chain approximation for effective quantum potentials,^210–213 dielectric formalism schemes,^{163,196–200,214–219} and quantum hydrodynamics.^{125,126,129,130,220}

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

The ab initio path-integral Monte Carlo (PIMC) approach^221–224 is one of the most successful methods for the description of quantum degenerate many-body systems. It is based on an evaluation of the partition function in coordinate representation, which, in the case of an electron gas in the canonical ensemble (inverse temperature $β=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

Here, the summation over all elements σⁱ 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,²²¹ 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 formula²²⁸ 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.²²⁹ 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.²³⁰ For indistinguishable quantum particles—bosons or fermions—this also requires the sampling of all possible permutation topologies,²³¹ which can be achieved efficiently using different implementations^72,232 of the worm algorithm by Boninsegni et al.^233,234

In the case of bosons (or hypothetical distinguishable quantum particles, which are sometimes referred to as Boltzmannons in the literature^235,236), the sign function in Eq. (7) is always positive, and the PIMC method allows for the quasi-exact simulation of up to $N∼104$ 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.⁷⁸ 

For completeness, we note that the sign problem can formally be avoided by imposing restrictions on the nodal structure of the thermal density matrix; this restricted PIMC method developed by Ceperley⁷⁹ 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,⁶⁸ 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.

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

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.¹¹ 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 generalization⁸⁵ of the Hohenberg–Kohn theorems.²⁴⁵ This enables a thermodynamic description of many-electron systems within the context of DFT by the grand-canonical potential as a density functional,

where

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

Minimizing $Ω[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.²⁴⁶ Here, one defines an auxiliary system of noninteracting particles described by a set of noninteracting single-particle Schrödinger equations,

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

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

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

with $β=1/kBT$ and μ is the chemical potential.⁵⁹ In combination, Eqs. (12) and (15) give one direct access to $KS[n]=−1/2∑α,k∫dr ϕα,k*(r)∇2ϕα,k(r)$ and $SS[n]=−∑α,k(fα,k ln (fα,k)+(1− fα,k) ln (1−fα,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,²⁴⁷ which makes both the benchmarking of existing functionals^90–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 function^162,248 given by

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

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;¹⁶³ 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.²⁵³ 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).

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

where we have $vS[n](r,t)=vext(r,t)+vH[n](r,t)+vxc[n](r,t)$⁠. In what follows, we shall abbreviate the effective Hamiltonian as $Ĥ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 = t_i, typically from the solution to a static Mermin–Kohn–Sham DFT problem. The time-dependent density is

and these occupation factors are typically consistent with the initial condition and held constant over the course of the dynamics. Among the three terms in v_S, $vext$ describes the attractive electron–ion interaction as well as impressed perturbations that drive the dynamics, v_H is the time-dependent Hartree potential, and v_xc 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 effects^258,259 that are absent in approaches, in which the Kubo–Greenwood (KG) formula is directly applied to the Kohn–Sham orbitals.^260–262 Furthermore, the computational details of these formulations differ in such a way that RT-TDDFT can be made asymptotically less expensive than LR-TDDFT. To understand when and how RT-TDDFT is most efficient, we first need to describe typical approaches to the numerical solution of Eq. (18).

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

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

While conceptually straightforward, the practical implementation of this approach requires the numerical approximation of a matrix exponential and, thus, the diagonalization of $Ĥ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 scalable^{259,263–265} than implementations that rely directly on dense linear algebra, including many implementations of LR-TDDFT. The precise conditions under which an RT-TDDFT linear response calculation is more efficient than an LR-TDDFT calculation are sensitive to many details of the calculation, but a large number of thermally occupied orbitals common in WDM calculations tend to favor RT-TDDFT as temperature increases.

While a scalar perturbation of the form $vext(r,t)=v0 exp (iq·r)f(t)$ can be used to compute the DSF at non-zero $q$⁠,^189,225 the calculation of the conductivity in the $q→0$ 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

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

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

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,

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

The Green's function (24) gives access to all one-particle observables, such as the density matrix, which follows from the time-diagonal Green's function, $nji(t)=ℏ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),

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

where $hHF$ contains kinetic, potential, and mean-field (Hartree–Fock) energy contributions, whereas correlation effects are included in the self-energy $Σ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 v_xc 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 v_xc 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.²⁰³ 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),

and is, to the lowest order (linear response), proportional to the excitation. In equilibrium, the density response function is stationary, $χ(q,t,t′)=χ(q,t−t′)$⁠, 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.

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 Ξ²⁶⁸ is given by Kwong and Bonitz²⁰³ and Bonitz et al.,²⁸⁰ 

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

In the limit of small A, the linear response theory (LRT)⁵⁹ becomes valid, and the induced density is fully described by the simple linear relation,

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 harmonics¹⁷⁵ or mode-coupling between multiple perturbations¹⁷⁷ 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),⁵⁹ 

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.⁵⁹ 

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 models^{147,148,151,186} that decompose the full DSF into contributions from bound and free electrons and their respective transitions directly work with $S(q,ω)$⁠. 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

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,²⁸² 

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

where $n(2)(r1,r2)=n(r1)n(r2)g(r1−r2)$ 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.²⁸³ Finally, the adiabatic-connection formula relates Eq. (34) with the XC contribution $fxc$ to the free energy,⁵⁹ 

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,²⁸⁴ which is a very useful tool at zero temperature.⁵⁹ In the context of the classical and quantum OCP, the coupling parameter forms of the adiabatic-connection formula are nearly exclusively employed,^28,283,285

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

In practice, Eqs. (34) and (35) constitute the prime motivations for the field of the self-consistent dielectric formalism.^{163,196–201,214–219,286,287} Here, the basic idea is to develop approximate expressions for the LFC $G(q,ω)$⁠, 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 liquid²⁸⁸ with $rs≳10$⁠, with notable progress reported by Tanaka²⁰⁰ and Tolias et al.²⁰¹ These schemes are based on the combination of the non-perturbative integral equation theory of classical liquids with the quantum dielectric formalism. The latter scheme²⁰¹ simplifies to the former²⁰⁰ 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,²⁹¹ or the further analysis of a possible experimental detection of the roton feature^252,292 in the dilute electron gas.

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

As an alternative route to obtain insights into the local field correction without any dielectric approximations, Moroni et al.²⁵³ 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 theorem⁵⁹ 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,²⁹⁴ and neutron matter.²⁹⁵ 

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

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

Furthermore, it is straightforward to invert Eq. (4) to obtain the static local field correction $G(q)$⁠. This has allowed Moroni et al.²⁸¹ 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.²⁹⁸ 

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

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

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

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,

Traditionally, Eq. (41) constitutes the starting point for a so-called analytic continuation,³⁰⁸ 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.³⁰⁸ 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.²⁵² 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 DSF⁵⁹ 

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 P_m being the probability to occupy the initial state m, $‖nml(q)‖$ is the corresponding transition matrix element, and $ωlm=(El−Em)/ℏ$ denotes the energy difference. Inserting Eq. (42) into Eq. (41) gives the analogous representation in the τ-domain,

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 UEG²⁹² will directly manifest as a less pronounced decay with τ. Conversely, it holds that $ω(q)∼q2$ in the non-collective single-particle regime with $q≫qF$ (with $qF$ being the usual Fermi wave number⁵⁹), which will lead to a steeper decay of the ITCF. Dornheim et al.¹⁹⁰ have suggested characterizing this with a relative decay measure of the following form:

and the corresponding results are shown in Fig. 10.

In addition to being of great value for the interpretation of simulation results for $F(q,τ)$⁠, these considerations open up new avenues for the interpretation of XRTS experiments with WDM.¹⁴⁶ 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

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

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

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 as^59,153,320

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

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

In addition to being interesting in their own right, the frequency moments constitute exact constraints that assist in the analytic continuation that is necessary for the reconstruction of the DSF from the QMC ITCF data, see the inversion of the Laplace transform of Eq. (41).^{252,306,317,318} They also constitute the key input for the non-perturbative method of moments that is used by Tkachenko and co-workers to estimate $S(q,ω)$ 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.¹⁷⁶ 

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

where $1={r→1,t1}, {l}=(1,2,…,l), d1={dr→1,dt1}$⁠. After writing the first three terms explicitly and setting $χ(1,2)≡w(1)(1,2)$⁠, we have^171,181

Here, all time integrations run from $−∞$ to t₁, 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′)=−i⟨T{ψ+(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

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

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

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

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

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

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

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