An accurate theoretical description of nonideal (i.e., interacting) quantum many-body systems is centrally important in numerous research fields in physics, quantum chemistry, material science, and related disciplines. While the basic equations governing quantum mechanics have been known for around a century, they are usually too complex to be solved in practice, even in the case of a few particles. The first attempts to circumvent this obstacle were based on uncontrolled approximations such as the Hartree–Fock approach, where the difficult treatment of correlations is abandoned in favor of a substantially simplified mean-field picture.1 Nevertheless, such mean-field methods have given important qualitative insights into various phenomena, such as collective plasmon excitations2–4 and Bose–Einstein condensation.5,6

Over past decades, the exponential increase in available computing time has sparked a remarkable surge of activity in fields such as computational physics and computational chemistry. In particular, state-of-the-art numerical methods often allow one to drastically reduce or even completely avoid approximations and simplifications. In this regard, a case in point is the density functional theory (DFT) approach.7 More specifically, DFT combines an often sufficient level of accuracy with manageable computational cost, which arguably makes it the most successful electronic structure method available. Indeed, a sizeable fraction of the world’s supercomputing time is being spent on DFT calculations, and the number of DFT-based scientific publications has grown exponentially in recent years.8 In addition, the quantum Monte Carlo (QMC) paradigm,9–11 which is computationally even more expensive, can even give exact results in many situations,12 both at finite temperature13 and in the ground state.

From a practical perspective, a particularly important property of quantum many-body systems is their response to an external perturbation.14 Such linear-response properties are readily measured in scattering experiments15 and can (in principle) give access to the full thermodynamic properties of a system. In this context, the key quantity is the dynamic structure factor (DSF)4S(q,ω)=∫−∞∞dteiωtF(q,t),(1)where the intermediate scattering function for wave vector q is defined as the correlation of two density operators in reciprocal space, i.e.,F(q,t)=⟨n̂(q,t)n̂(−q,0)⟩.(2)For example, neutron scattering experiments with ultracold helium have given invaluable insights into the collective excitations of superfluids5 and normal quantum liquids,16–18 including the distinct roton feature19,20 present at intermediate wave numbers. A second important example is the diagnosis of warm dense matter (WDM), an extreme state21–23 that occurs naturally in astrophysical objects such as the interiors of giant planets24 and is important for technological applications such as inertial confinement fusion.25 Here, x-ray Thomson scattering (XRTS) experiments26,27 are a highly important diagnostic method and give insights into important system parameters such as the temperature T28,29 or the charge state Z.

Unfortunately, theoretical modeling of the dynamic properties of correlated quantum many-body systems, such as S(q, ω), is notoriously difficult. In practice, DFT cannot give direct access to many-particle density correlation functions such as Eq. (2), and time-dependent formulations30 require additional approximations, such as the dynamic exchange–correlation (XC) kernel or a time-dependent XC potential. Furthermore, exact QMC methods are usually restricted to the imaginary-time domain. For example, the ab initio path-integral Monte Carlo (PIMC) approach9,13,31 gives straightforward access to F(q, τ), corresponding to the intermediate scattering function evaluated at the imaginary time t = −iτ with τ ∈ [0, β] and β = 1/(k_bT). Note that we assume Hartree atomic units throughout this work. The connection to the sought-after dynamic structure factor is then given by a two-sided Laplace transform, i.e.,F(q,τ)=LS(q,ω)=∫−∞∞dωe−ωτS(q,ω).(3)Equation (3) is the basis for an analytic continuation (AC),32 i.e., the numerical inversion of LS(q,ω) to compute S(q, ω). This problem is well known but also notoriously difficult, the latter being because any noise in the QMC data for F(q, τ) leads to instabilities and ambiguity in the DSF. Despite these formidable difficulties, AC is still one of the most promising routes to capturing rigorously the complex interplay of nonideality, quantum degeneracy, and possibly thermal excitations. Consequently, there are many AC methods33–44 with different strengths and weaknesses. Unfortunately, different methods must often be benchmarked against each other,45 and so the accuracy of the reconstructed spectra generally remains unclear.

The present work is aimed at partially circumventing such challenges in describing the dynamics of correlated quantum many-body systems. More specifically, we argue that because of the uniqueness of the two-sided Laplace transform, the imaginary-time density–density correlation function (ITCF) F(q, τ) contains the same information as does S(q, ω) itself, only in a form that might be unfamiliar at first glance.46 While the traditional way of doing physics in the frequency domain emerges naturally (e.g., from scattering experiments detecting an energy-resolved signal), it is not the only—or necessarily preferred—option in practice because literally all physical concepts known from the ω domain of S(q, ω) have an analog in the τ domain of F(q, τ). Indeed, many features such as sharp quasi-particle (QP) peaks in S(q, ω) can (in principle) also be identified in the τ domain. Moreover, we stress that imaginary time is directly connected to the physical concept of quantum-mechanical delocalization, which means that F(q, τ) gives straightforward insights that can only be observed indirectly in the DSF. Consequently, our aim is not to find good approximations to concepts derived in the ω domain but instead to highlight how physical effects are manifested directly in F(q, τ).

While both domains are formally equivalent, they tend to emphasize different aspects. Consequently, when viable, one should analyze both S(q, ω) and F(q, τ) to obtain a more complete physical picture. In practice, working with F(q, τ) instead of S(q, ω) has a number of important benefits that are summarized in Fig. 1 for the case of XRTS experiments. Traditionally, one would construct an approximate model for the DSF and then convolve it with the instrument function R(ω) that characterizes the laser beam and detector for comparison with the experimentally measured intensity signal I(q, ω),47 i.e.,I(q,ω)=S(q,ω)⊛R(ω).(4)A straightforward deconvolution of the experimental signal to obtain an experimental S(q, ω) is generally impossible because of numerical instabilities.

In stark contrast, these obstacles are completely absent in the τ domain. First, note that computing the two-sided Laplace transform of the experimental signal is both easy and numerically well behaved.28,29 Moreover, the deconvolution is trivial in the τ domain because here the convolution works as a simple multiplication; thus we findLS(q,ω)=LS(q,ω)⊛R(ω)LR(ω),(5)where both the numerator and denominator on the RHS of Eq. (5) are easily evaluated in practice. Thus, it is possible to obtain experimental results directly for F(q, τ), which in turn give model-free access to physical parameters such as the temperature.28,29 In addition, QMC methods such as PIMC can give exact theoretical results for F(q, τ), which opens up the enticing possibility for unprecedented accurate comparisons between theory and experiment. For completeness, note that it is also possible to translate well-known approximate models for S(q, ω) such as the widely used Chihara decomposition48 into the τ domain, where they can be benchmarked against more-accurate simulation data.

The present aim is to provide the first systematic overview of a number of recent developments focusing on different aspects of the ITCF.28,29,46,49,50 In addition, we connect the imaginary-time domain with other concepts such as the roton feature in the strongly coupled uniform electron gas (UEG),19 and we introduce the conceptual basis for various new developments, including the extraction of QP energies and the possibility of comparing experimental measurements of F(q, τ) both with exact and approximate simulations.

This paper is organized as follows. In Sec. II, we introduce the relevant theoretical background starting with the UEG model51,52 that we consider throughout this investigation. However, we emphasize that all conclusions drawn about the physical insights from the ITCF are completely general and in no way particular to the UEG. Section II B is devoted to the ab initio PIMC method9,13,53 and how it gives straightforward access to F(q, τ). In addition, in Sec. II C we discuss the connection between the DSF and ITCF in linear-response theory. This is followed in Sec. II D by a comprehensive overview of some general properties of F(q, τ). In Sec. III, we present our new results, starting with a qualitative investigation of some general trends based on synthetic trial spectra S(ω) in Sec. III A. This is followed by a detailed investigation of exact PIMC results for F(q, τ), which give important insights into a number of physical processes including the connection between temperature and imaginary-time diffusion, as well as the nontrivial roton feature in the DSF.19 The paper concludes with a summary and outlook in Sec. IV.

The UEG51,52 (also known as the jellium in the literature) is the quantum version of the classical one-component plasma.54,55 More specifically, the Hamiltonian of the UEG with N electrons is given byĤ=−12∑l=1N∇l2+∑l<kNϕE(r̂l,r̂k)+N2ξM,(6)with ξ_M being the usual Madelung constant, and the Ewald pair potential ϕE(r̂l,r̂k) taking into account both the interaction between two electrons (and their infinite periodic array of images) and that with their respective positive neutralizing homogeneous background. See Fraser et al.56 for an extensive and accessible discussion of the Ewald potential for the UEG.

From a physical perspective, the UEG can be characterized fully by three dimensionless parameters.57 (1) The density parameter (also known as the Wigner–Seitz radius in the literature) rs=r̄/aB, with r̄ being the average particle separation and a_B being the first Bohr radius, plays the role of the quantum coupling parameter. Consequently, the UEG attains the limit of a noninteracting Fermi gas for r_s → 0 and forms a Wigner crystal for r_s ≳ 100.58–61 (2) The degeneracy temperature Θ = T/T_F, with T_F being the Fermi temperature,4,52,57 determines whether the UEG is fully quantum degenerate (Θ ≪ 1) or semi-classical62 (Θ ≫ 1). (3) The spin-polarization parameter ξ = (N − N)/N, where N and N are the numbers of spin-up and spin-down electrons, respectively; here we restrict ourselves to the fully unpolarized (i.e., paramagnetic) case of ξ = 0. The WDM regime is typically defined by the condition r_s ∼Θ ∼ 1.

The UEG is one of the most fundamental model systems in physics, quantum chemistry, and related fields. Indeed, accurate parameterization of various UEG properties63–68 based on QMC simulations both in the ground state69–71 and at finite temperature72–76 has been fundamentally important for a wide spectrum of applications, most notably as input for DFT simulations of real materials.

Since its original inception for describing ultracold ⁴He in the 1960s,77 the ab initio PIMC method has become one of the most successful simulation tools in statistical physics and related fields.9,13,53 The basic idea is to express the canonical partition function (i.e., volume V, number density n = N/V, and inverse temperature β are fixed) in coordinate space, which for an unpolarized electron gas givesZβ,N,V=1N↑!N↓!∑σ↑∈SN↑∑σ↓∈SN↓sgn(σ↑,σ↓)×∫VdR〈R|e−βĤ|π̂σ↑π̂σ↓R〉.(7)Here, the metavariable R=(r1,…,rN)T contains the coordinates of all N = N + N electrons in the system. Note that the antisymmetry of the fermionic electrons with identical spin orientation is taken into account by the sums over all possible permutation elements σ (i ∈ {↑, ↓}) of the respective permutation group SNi, and the corresponding permutation operators π̂σi. As can be seen, the density operator ρ̂=e−βĤ in Eq. (7) can be interpreted straightforwardly as a propagation in imaginary time by an interval of t = −iβ.

Unfortunately, direct evaluation of Eq. (7) is precluded by the absence of knowledge about the matrix elements 〈R|e−βĤ|R′〉, given that the kinetic and potential contributions to the full Hamiltonian Ĥ=K̂+V̂ do not commute. To circumvent this obstacle, we use the exact semigroup propertye−βĤ=e−ϵĤP(8)with ϵ = β/P, leading toZβ,N,V=1N↑!N↓!∑σ↑∈SN∑σ↓∈SNsgn(σ↑,σ↓)×∫VdR0…dRP−1〈R0|e−ϵĤ|R1〉×〈R1|…〈RP−1|e−ϵĤ|π̂σ↑π̂σ↓R0〉.(9)

As can be seen, the single imaginary-time propagation has been replaced by P shorter steps of length ϵ. From a practical perspective, the main point is that each matrix element must be evaluated at P times the original temperature. In the limit of large P, one can introduce a suitable high-temperature factorization, and the error can be made arbitrarily small. A detailed discussion of different factorization schemes is beyond the present scope but is given in Refs. 78 and 79. For the parameters of interest herein, it is sufficient to restrict ourselves to the simple primitivefactorization80e−ϵĤ=e−ϵK̂e−ϵV̂+Oϵ−2.(10)Note that the convergence of Eq. (10) is ensured by the Trotter formula,81 and that the convergence with P is checked carefully in practice. Inserting Eq. (10) into Eq. (9) then gives the final expression for the partition function, i.e.,Z=∫VdXW(X)=∫VdX∏α=0P−1WV(Rα)WK(Rα,Rα+1),(11)where the new metavariable X=(R0,…,RP−1)T contains the coordinates of all particles on all P imaginary-time slices, and it holds that R₀ = R_P. For simplicity, we assume that the summation over all permutations as well as the combinatorial normalization factors are contained implicitly in dX. The weight function W(X) can now be evaluated explicitly; see the second line in Eq. (11). More specifically, the potential term is given byWV(Rα)=exp−ϵ2∑k≠lNW(rl,α,rk,α),(12)with W(r, s) being the pair interaction between two particles [combining both the Ewald pair potential and the Madelung constant from the UEG Hamiltonian of Eq. (6)], and r_l,α being the coordinate of the lth particle on imaginary-time slice α. The kinetic term is off-diagonal and corresponds to the density matrix of a noninteracting system, i.e.,WK(Rα,Rα+1)=exp−∑l=1N(rl,α−rl,α+1)22σϵ22πσϵ3N/2.(13)

A graphical illustration of Eq. (11) is presented in Fig. 2, where we show a fictitious configuration of N = 4 particles in the τ–x plane. First and foremost, note that each particle is now represented by an entire closed path along the imaginary-time direction. This is the well-known classical isomorphism,83 given that we have effectively mapped the complicated quantum system of interest onto a system of classical ring polymers. The beads of each particle interact via the usual pair potential; see Eq. (12). Furthermore, the beads of the same particle but on adjacent imaginary-time slices interact via the harmonic spring potential introduced in Eq. (13). In particular, the characteristic width of the corresponding Gaussian is directly proportional to the thermal wavelength σϵ=λϵ/2π, withλϵ=λβPandλβ=2πβ.(14)In practice, this spring potential determines the diffusion throughout imaginary time; see the yellow Gaussian curve at the RHS of Fig. 2. For high temperature, λ_β is small and the diffusion is severely restricted. In the classical limit of β → 0, the Gaussian becomes infinitely narrow, and the paths collapse to point particles. In contrast, the paths become more extended for low temperatures, which is a direct manifestation of quantum delocalization. The nontrivial interplay of the potential and kinetic terms in Eq. (11) then gives an exact description of quantum diffraction. The basic idea of the PIMC method is to randomly generate all possible paths (including all different possible permutations;31 e.g., see the permutation cycle including two identical particles in the center of Fig. 2) using a particular implementation of the Metropolis algorithm.84 In practice, we use an off-diagonal canonical version85 of the worm algorithm idea by Boninsegni et al.86,87 for all calculations in this work. An additional obstacle is given by the antisymmetry of identical fermions under the exchange of particle coordinates [cf. Eq. (7)], which implies that the configuration weight W(X) can be both positive and negative. This is the origin of the notorious fermion sign problem,88–90 which leads to an exponential increase in computing time with system parameters such as N or β. While there are various approximate frameworks for dealing with the sign problem,74,91–95 here we carry out exact unrestricted PIMC simulations, so our simulations are computationally demanding but exact within the given Monte Carlo error bars. See Ref. 90 for a detailed discussion of the sign problem.

Next, we consider the estimation of imaginary-time correlation functions within the PIMC paradigm.33,82,96 In particular, the imaginary-time version of the intermediate scattering function is defined as the ITCF, i.e.,F(q,τ)=⟨n̂(q,0)n̂(−q,τ)⟩.(15)Naturally, the evaluation of Eq. (15) requires the correlated evaluation of the density operator in reciprocal space at two different imaginary-time arguments. Using the factorization introduced in Eq. (9), it can be seen that the PIMC method allows for a straightforward evaluation of F(q, τ) at integer multiples of the factorization step ϵ,82 i.e.,F(q,τj)=1Zβ,N,V1N↑!N↓!∑σ↑∈SN↑∑σ↓∈SN↓sgn(σ↑,σ↓)×∫dR0…dRP−1〈R0|n̂(q)e−ϵĤ|R1〉〈R1|e−ϵĤ|R2〉…〈Rj|n̂(−q)e−ϵĤ|Rj+1〉…〈RP−1|e−ϵĤ|π̂σ↑π̂σ↓R0〉.(16)

In practice, the τ grid can be made arbitrarily fine by increasing the number of high-temperature factors P, with a linear increase in the required computing time. In contrast, the wave-vector grid is determined by the system size,72,97–99 with q = 2πL⁻¹n and n ≠ 0 being an integer vector.

A central relation in the context of the present study is given by Eq. (3), which unambiguously connects the DSF S(q, ω) to the ITCF F(q, τ). By definition, both quantities thus contain exactly the same information because the two-sided Laplace transform is a unique transformation.

A second crucial relation is given by the well-known fluctuation–dissipation theorem (FDT)4S(q,ω)=−Imχ(q,ω)πn(1−e−βω),(17)which links the DSF to the dynamic density response function χ(q, ω) known from linear-response theory. It is convenient to express the latter as100χ(q,ω)=χ0(q,ω)1−4πq21−G(q,ω)χ0(q,ω),(18)with χ₀(q, ω) being the analytically known density response function of a noninteracting system, and G(q, ω) being the dynamic local field correction that contains the full information about XC effects. Consequently, setting G(q, ω) ≡ 0 in Eq. (18) leads to a description of the density response at the mean-field level, which is typically known as the random phase approximation. Evidently, it is straightforward to compute the ITCF from any dielectric theory for G(q, ω)101–110 by inserting the corresponding (static or dynamic) local field correction into Eq. (18), evaluating the FDT [Eq. (17)], and finally computing the two-sided Laplace transform [Eq. (3)] to obtain F(q, τ).

Finally, note that the ITCF is connected to the static limit of Eq. (18) via a simple one-dimensional integral, i.e.,χ(q,0)=−n∫0βdτF(q,τ),(19)which is sometimes referred to as the imaginary-time version of the FDT;111 see also the Appendix herein for a short derivation.

In this subsection, we summarize a number of properties of the ITCF. As a starting point, we consider the detailed balance relation of the DSF,4,26 i.e.,S(q,−ω)=S(q,ω)e−βω,(20)which holds for uniform systems in thermodynamic equilibrium. Inserting Eq. (20) into the two-sided Laplace transform [Eq. (3)] then directly yields the important symmetry relation28F(q,τ)=∫0∞dωS(q,ω)e−ωτ+e−ω(β−τ)=F(q,β−τ).(21)Evidently, Eq. (21) implies that F(q, τ) is symmetric about τ = β/2. Therefore, any knowledge about the ITCF—e.g., from an XRTS measurement26–29—allows one to read the temperature directly from the plot by locating its extremum, which is always a minimum; no theoretical model, simulation, or indirect inference is required.

Another useful property of a correlated quantum many-body system is given by the frequency moments of the DSF, which we define asMS(α)=⟨ωα⟩S=∫−∞∞dωS(q,ω)ωα.(22)In particular, the odd moments can in principle be obtained by evaluating nested commutator terms.112 In the quantum-mechanical case,100 the cases of α = −1, 1, 3 are known from such sum rules. For example, the well-known f-sum rule gives the important relation4MS(1)=⟨ω1⟩S=q22.(23)To connect Eq. (22) to the ITCF, we differentiate the latter with respect to τ, which givesdndτnF(q,τ)=(−1)n∫−∞∞dωωne−τωS(q,ω).(24)Therefore, the derivatives of F(q, τ) at the origin give straightforward access to both odd and even frequency moments,50,113 i.e.,MS(α)=−1α∂α∂ταF(q,τ)τ=0(25)for α ≥ 0. In practice, knowledge about such frequency moments is important for further constraining a potential AC from the imaginary-time domain to S(q, ω); e.g., see Refs. 114 and 115. In fact, the DSF is fully determined by MS(α), which is known as the Hamburger problem.116,117 Therefore, the frequency moments [and thus the derivatives of F(q, τ) around the origin contain important physical insights as discussed in Sec. III. Very recently, Eq. (25) was used in Ref. 49 to evaluate the f-sum rule in the imaginary-time domain, which gives access to the absolute normalization of the XRTS intensity and therefore to the electron–electron static structure factor S(q) = F(q, 0).

As a final useful relation, we consider the exact spectral representation of the DSF, which is given by4S(q,ω)=∑m,lPmnml(q)2δ(ω−ωlm).(26)Here, l and m denote the eigenstates of the Hamiltonian, ω_lm = (E_l − E_m)/ℏ the corresponding energy difference, and n_ml the transition element from state m to l due to the density operator n̂(q). In addition, Pm=e−Emβ/Z is the probability of occupying the initial state m. The two-sided Laplace transform of Eq. (26) then gives an analogous spectral representation of the ITCF, i.e.,F(q,τ)=∑m,lPmnml(q)2e−τωlm=∑Em<Elnml(q)2Pme−τωlm+Pleτωlm,(27)where in the second line we have ordered the eigenstates according to their energy. At finite temperature (i.e., T > 0, leading to a finite value of β), transitions between eigenstates occur in both directions. Specifically, excitations that increase the energy always lead to exponential decay with τ (see also Appendix D of Ref. 118), whereas conversely, excitations that reduce the energy lead to a corresponding exponential increase. Note that this is reflected directly in the symmetry relation of Eq. (21).

To get a general feeling for the relationship between S(q, ω) and F(q, τ), we start by considering synthetic DSFs of a simple Gaussian form that obey the detailed balance relation of Eq. (20), i.e.,S(ω)=G(ω;ω0,σ)if ω≥0,eβωG(−ω;ω0,σ)otherwise,(28)with the usual definition of a normalized Gaussian, i.e.,G(ω;ω0,σ)=exp(ω−ω0)22σ22πσ2.(29)The corresponding two-sided Laplace transform [Eq. (3)] is given byF(τ)=A(τ;ω0,σ)+A(β−τ;ω0,σ),(30)with the definitionA(τ;ω0,σ)=e−τω0eτ2σ2/21+erfω0−σ2τ2σ22.(31)Note that any nascent delta function that includes a parameter controlling the width is acceptable for the present analysis and leads to identical conclusions; however, Eq. (29) is particularly convenient because of its simple Laplace transform.

We start the investigation of such synthetic spectra by considering the impact of the peak width σ. The results are shown in Fig. 3, with the top (resp. bottom) panel showing the results for the DSF with β = 1 and ω₀ = 1 [resp. the corresponding F(τ)]. For σ = 1, the DSF consists of a single broad peak combining the positive and negative frequency range; note that the damping for ω < 0 is a quantum effect described by the detailed balance relation of Eq. (20). This manifests as a more pronounced decay with τ in the ITCF compared to the other depicted examples. For σ = 0.5, there remains some significant overlap in S(ω) around ω = 0, but overall we find a double-peak structure in the DSF. This pronounced change in the DSF is translated directly to F(τ) decaying less steeply with τ. Finally, the dotted blue curve is the result for σ = 0.1, which leads to two narrow peaks in S(ω). In practice, such features are often interpreted as distinct QP excitations. Indeed, our synthetic model approaches a delta-like QP excitation in the limit of σ → 0, i.e.,SQP(ω)=δ(ω−ωQP)+e−βωQPδ(ω+ωQP),(32)with ω_QP being the QP energy, because Eq. (29) is a nascent delta-function. It is straightforward to show that the corresponding ITCF is given byFQP(τ)=e−τωQP+e−(β−τ)ωQP.(33)The results for Eq. (33) are shown as the solid yellow curve in the bottom panel of Fig. 3 and agree excellently with the ITCF of the Gaussian trial DSF with σ = 0.1. This finding has important consequences for the practical interpretation of scattering experiments. For example, it is well known that the XRTS signals of a WDM probe exhibit a sharp plasmon excitation at the plasma frequency ω_p for small wave number q = |q|. In this way, the location of the plasmon gives direct insight into the density of the unbound electrons that take part in this collective excitation and therefore is an import diagnostic tool. In practice, this endeavor is complicated by the convolution with the instrument function, which in combination with other features in the DSF might mask the true location of the plasmon. In contrast, Eq. (33) implies that we can extract the plasma frequency directly from an exponential fit to F(q, τ). As explained in Fig. 1, the latter can be computed without the bias due to the instrument function because the convolution is just a multiplication in the τ domain. We thus conclude that working with F(q, τ) not only is advantageous for diagnosing the temperature of a sample (see the recent work by Dornheim et al.28 and the discussion of Fig. 5 below) but also can potentially give access to the free electronic density in a forward-scattering experiment where the scattering angle and thus the wave number are small enough to probe the regime of collective plasma excitations.

Other concepts such as approximate expressions for the plasmon location at finite wave numbers119 that are sometimes used to diagnose plasma parameters26,120,121 can also be translated directly into the τ domain. In general, however, we argue against the blind transformation of ω-native concepts to F(q, τ) because it makes more sense to work directly with the rich physics that is naturally inherent to F(q, τ).

We continue our investigation of synthetic spectra by investigating the manifestation of the peak position ω₀ shown in Fig. 4. In fact, the impact of the peak position can qualitatively be seen immediately from Eqs. (27) and (33). Specifically, transferring spectral weight from larger to smaller excitation energies ω leads to a less steep decay along τ (for τ < β/2). This is indeed what we find in the bottom panel of Fig. 4.

From a physical perspective, this finding has important implications, as we shall see in the discussion of real physical results for the UEG in Sec. II A. First, it means that the persistence of two-body correlations throughout the imaginary time 0 ≤ τ ≤ β/2—which we can observe directly in our PIMC simulations—implies a downshift in the dominant excitation energies in S(q, ω); see the spectral representation in Eq. (26) above. This is strongly related to the roton feature in the spectrum of density fluctuations in both the UEG19 and quantum liquids such as ultracold helium.5,16–18 In other words, quantifying the decay of correlations in F(q, τ) is a straightforward alternative to the usual pseudo dispersion relation ω(q) constructed from the position of the maximum in the DSF.122 In addition, the aforementioned effects lead to a maximum in the static linear density response function χ(q, 0). This can be seen either from the imaginary-time version of the FDT [Eq. (19)] or by recalling the relation between χ(q, 0) and the inverse frequency moment of the DSF,41,42 i.e.,MS(−1)=−χ(q,0)2n.(34)

We conclude this investigation of synthetic spectra by considering the impact of the temperature β. A corresponding analysis is shown in Fig. 5 for three slightly different values of β. Indeed, the main difference between the spectra is given by the varying ratio of the peaks at positive and negative frequency. In stark contrast, we observe a pronounced influence of the temperature on F(τ) shown in the bottom panel of Fig. 5. Even though from a mathematical perspective S(q, ω) and F(q, τ) are completely equivalent representations of the same information, both domains emphasize different aspects of it. For example, the peak width of the DSF is a concept of the ω domain and can be seen most clearly in S(ω); see Fig. 3. On the other hand, the τ domain is intimately connected to the temperature of a system, which manifests as the somewhat subtle detailed balance relation of Eq. (20) in the DSF. Meanwhile, the corresponding symmetry [Eq. (21)] of the ITCF is enhanced substantially. From a physical perspective, this is not surprising because the (inverse) temperature determines the characteristic variance of the imaginary-time diffusion process (cf. Fig. 2) and therefore decisively shapes the decay of correlations with τ. In practice, this means that the τ domain is the representation of choice for extracting the temperature from an XRTS measurement;28 see also the recent Ref. 29 for a more quantitative analysis. Last, we note that the utility of F(q, τ) has been confirmed by the independent reanalysis of an XRTS measurement of warm dense beryllium123 by Schörner et al.124

Next, we turn our attention to exact simulation results for a physical system. In Fig. 6, we show our ab initio PIMC results for F(q, τ) for the UEG at r_s = 10 and at different values of the degeneracy temperature Θ. Note that it is sufficient to consider the interval 0 ≤ τ ≤ β/2 because of the symmetry relation of Eq. (21). These parameters are located at the margin toward the strongly coupled electron liquid regime;4,41,125 still, we observe almost no correlation-induced features in the static structure factor S(q) = F(q, 0) for all depicted values of Θ. First and foremost, we observe a decay along the τ direction in the depicted interval for all values of the wave number q. This is expected because correlations can only become weaker—or in the extreme case remain unaffected—because of the imaginary-time diffusion process that is sampled in our PIMC simulations. In addition, this decay of correlations is clearly more pronounced at low temperatures. This can be explained directly by recalling the discussion of the path sampling in Fig. 2.

In Fig. 7, we show snapshots of PIMC simulations of N = 14 unpolarized electrons at r_s = 10 for Θ = 1 (left column) and Θ = 4 (right column). More specifically, the top row shows path configurations (with the red and blue lines corresponding to spin-up and spin-down electrons, respectively) for the two different values of Θ. At the lower temperature, the paths exhibit a more pronounced diffusion along the τ direction compared to Θ = 4. This is expected and a direct consequence of the definition of the thermal wave length λ_β in Eq. (14). Naturally, the larger displacements in coordinate space with increasing τ ≤ β/2 lead to a decreasing density–density correlation function. For illustration, we also show the same configurations as paths in the 3D simulation box in the bottom row of Fig. 7. In this representation, the more pronounced imaginary-time diffusion process in the case of Θ = 1 manifests as more-extended paths. With increasing temperature, the paths become less extended and attain the limit of classical point particles in the limit of β → 0.

An additional trend that can be seen from Fig. 6 is that the imaginary-time decay of F(q, τ) becomes increasingly steep in the limit of large q. This marks the transition from the collective to the single-particle regime, which can be reproduced with remarkable accuracy by a simple Gaussian imaginary-time diffusion model; see Ref. 46 for an extensive discussion. This trend is investigated in more detail in Fig. 8, where we show F(q, τ) along the τ direction for two selected values of q. Specifically, the solid colored curves correspond to the different values of Θ. First and foremost, we note that we have not rescaled the τ axis by the respective values of the inverse temperature β as in Fig. 6. Therefore, the minima in the different curves are located at different positions in descending order of Θ. In fact, this representation gives a direct insight into the observed less-pronounced decay of F(q, τ) along the τ direction: it is not the consequence of a less steep decay by itself, but rather of the reduced imaginary time that is available for the diffusion process.

This is illustrated further by the dashed black lines, which depict a linear expansion of F(q, τ) around τ = 0 evaluated from the exact f-sum rule of Eq. (23). In particular, the lines are parallel because of Eq. (23), but the initial points are different because of the different static structure factors S(q) = F(q, 0). In addition, Eqs. (23) and (25) clearly predict a parabolically increasing decay with respect to q along the τ direction around τ = 0, which is reflected by the observed steep decay in the limit of large q in Fig. 6.

From a practical perspective, we find that the localization of the position of the minimum in F(q, τ) with respect to τ becomes more difficult for low temperatures in the single-particle regime. However, this is not a fundamental obstacle because the temperature can always be inferred via F(q, 0) = F(q, β), which also follows directly from the symmetry relation of Eq. (21).29

A further interesting research question is given by the dependence of F(q, τ) on the coupling strength. This is investigated in Fig. 9, where we show our PIMC results for F(q, τ) at the electronic Fermi temperature Θ = 1 for r_s = 4 (top), r_s = 10 (center), and r_s = 20 (bottom). In this case, the most pronounced trend is the increased structure in S(q), i.e., in the limit of τ = 0. This can be seen particularly well in the top panel of Fig. 10, where we show this limit for all three considered values of r_s. The inset shows a magnified segment around the peaks at r_s = 20 (blue) and r_s = 10 (red); no peak can be resolved for r_s = 4 (green).

The bottom panel of Fig. 10 shows the same information for the thermal structure factorF(q, β/2), where F(q, τ) attains its minimum value for all q. All three curves exhibit a qualitatively similar peak around q = 2q_F, which is indicative of a reduced decay along the τ direction. In turn, this implies a shift in the spectral weight of the DSF S(q, ω) toward lower frequencies; cf. the discussion above of Fig. 4. Interestingly, the peak height exhibits nonmonotonic behavior with respect to r_s and is minimal for r_s = 10. On the other hand, the peak position is increasing monotonically (albeit weakly) with r_s.

To investigate this effect further, we define a measure of relative τ decay asΔFτ(q)=F(q,0)−F(q,τ)F(q,0),(35)where the values of zero and unity indicate no decay and full decay, respectively. The results for ΔF_β/2(q) are shown in the top panel of Fig. 11. First, we see that all curves converge for large q and eventually attain the limit of lim_q→∞ΔF_β/2(q) = 1. This is a direct consequence of the vanishing value of F(q, β/2) in the single-particle limit. Around twice the Fermi wave number, the three curves in the top panel of Fig. 11 start to deviate from each other in a highly nontrivial way. The horizontal lines indicate the q → 0 limit of ΔF_β/2(q), which is determined by the sharp plasmon excitation at the plasma frequency ω_p; see also Eq. (33). For r_s = 4, the curve is relatively featureless and interpolates smoothly between the q = 0 and q → ∞ limits. For r_s = 10 and particularly r_s = 20, the curves become nonmonotonic and exhibit a minimum at intermediate wave numbers. In other words, the decay along the τ direction is suppressed when the wave number is of the order of the average particle separation d = 2r_s.

We have already mentioned several times that such a reduced decay implies a shift of spectral weight in the DSF to lower frequencies ω. This is indeed the case for the UEG under the present conditions; see the original Ref. 41 for all technical details about the corresponding calculations. The results for the wave-number dependence of the position of the maximum in the DSF, ω(q), are shown in Fig. 12. We observe a monotonic curve for r_s = 4 and a distinct roton minimum for r_s = 10 and 20 around the same position as that of the nonmonotonic behavior of ΔF_β/2(q) in Fig. 11.

To illustrate further the direct physical correspondence between the reduced τ decay on the one hand and the roton feature in the DSF on the other hand, we normalize ΔF_β/2(q) by its q → 0 limit in the bottom panel of Fig. 11. Note that this is analogous to the normalization with respect to the plasmon frequency of the dispersion relation ω(q) in Fig. 12. The resulting curves resemble even more closely the dispersion of the DSF for small to intermediate wave numbers, and they exhibit the same qualitative trend. In particular, we find a comparable roton minimum for r_s = 10 and 20. This is unsurprising given that both F(q, τ) and S(q, ω) contain by definition the same physical information. Therefore, any physical process such as the roton feature has to manifest itself in both the ω and τ domains. This is also evident from comparing Eqs. (26) and (27).

We conclude this investigation with a brief discussion about the physical mechanism behind the roton feature in both S(q, ω) and F(q, τ). Very recently, Dornheim et al.19 showed that the energy reduction in the spectrum of density fluctuations can be explained by the spatial alignment of pairs of electrons on different length scales. More specifically, the roton appears when the wave length λ = 2π/q of a density fluctuation is comparable to the average interparticle distance, and this is indeed the case for q ∼ 2q_F. In this case, the alignment of two electrons at distance λ leads to a reduction in the interaction energy of the system, which manifests as (i) the downshift of ω(q) in the DSF and (ii) the stability of spatial correlations along the imaginary-time propagation. Note that this pair alignment model was subsequently substantiated further in Refs. 126 and 127.

The present work was devoted to investigating the dynamic properties of correlated quantum many-body systems in the imaginary-time domain. In particular, it was argued that the usual approach in terms of the DSF S(q, ω) in the frequency domain is not the only option because by definition the ITCF F(q, τ) in the imaginary-time domain contains exactly the same information, albeit in an unfamiliar representation. Some properties such as the peak width of the DSF are more easily accessible in the ω domain, whereas other physical effects such as the temperature or quantum-mechanical delocalization are more clearly emphasized in the τ domain. In fact, even nontrivial physical effects such as the roton feature in the DSF19 can easily be identified in the ITCF.

Instead of attempting the notoriously difficult analytic continuation from QMC results for F(q, τ) to S(q, ω), we find that it can be highly advantageous to pursue the opposite direction. Indeed, transforming an experimental signal from the ω domain to the τ domain is straightforward and well behaved with respect to the inevitable experimental noise.28,29 In addition, it allows for an easy deconvolution of the instrument function; this is an important point, for example in XRTS experiments.26,47

First and foremost, we have demonstrated that the availability of such experimental results for F(q, τ) directly allows for a number of physical insights. In particular, it allows for the model- and simulation-free extraction of important system parameters such as the temperature28,29 and QP excitation energies such as the plasmon frequency ω_p. In addition, experimental data for F(q, τ) can be compared straightforwardly to exact QMC results for the same quantity, which opens up the enticing possibility for unprecedented agreement between theory and experiment.

We are convinced that the proposed physical interpretation of the dynamic properties of correlated quantum systems in the τ domain will have a strong impact in a number of disciplines. A prime example is the interpretation of XRTS experiments with WDM,128 which can be used for the systematic construction of model-free equation-of-state databases.129,130 In turn, the latter are of paramount importance for the description of astrophysical objects22 and inertial confinement fusion applications.131 In this context, we also note the advent of modern free-electron x-ray laser facilities such as the new European XFEL,132 which with their high repetition rate will allow for high-quality results for I(q, ω) and in this way also F(q, τ); see the outlook of Ref. 14 for more details. This will facilitate measurements at a number of scattering angles and thus wave vectors, which will be important for resolving physical dispersion effects such as the roton feature. Moreover, we stress that the basic idea to use F(q, τ) directly instead of S(q, ω) is in no way limited to either XRTS or the study of WDM, and it can easily be applied to a wide range of other systems such as ultracold atoms.5,17,18,114

In addition to its considerable value for interpreting experiments, we note that the τ domain also opens up the way for new developments in the theoretical modeling of quantum many-body systems. For example, a central limitation of linear-response time-dependent DFT30 is the absence of reliable external input for the dynamic XC kernel K_xc(q, ω)133 beyond the local density approximation or generalized-gradient expansions.134 While a number of model kernels exist for either the UEG135,136 or more-generic systems,137,138 they are often restricted to the static limit. In this regard, exact QMC-based input data for the imaginary-time dependence of density–density correlations can help in multiple ways.14 First, we propose pursuing linear-response imaginary-time-dependent DFT simulations of real materials, which can utilize exact fully τ-dependent information about XC effects based on ab initio QMC simulations. More specifically, one might either (i) combine a static material-specific XC kernel based on DFT calculations134 with the τ dependence extracted from a PIMC simulation of either the UEG or a more realistic system or (ii) directly attempt a generalization of suitable PIMC simulation data. While it is certainly true that some features about the dynamic density response of a given system of interest might be unresolvable in the τ domain (e.g., the width of a peak in the DSF), often one is primarily interested in frequency-integrated properties such as the electron–electron static structure factor S_ee(q) or its inverse Fourier transform, which is given by the electron–electron pair correlation function. To this end, operating in the τ domain is not a disadvantage and might allow for the extraction of electron–electron correlations from DFT simulations with high fidelity.14 Another interesting route is an explicit time-dependent DFT propagation along imaginary time,139 although this will likely be more challenging in practice.

Finally, note that the analysis of ITCFs can be extended straightforwardly to higher-order density correlators.82 These higher-order ITCFs can easily be estimated in PIMC simulations and will give new insights into dynamic three-body and four-body correlation functions.140 Moreover, such higher-order ITCFs are directly related to nonlinear response properties.73,141–143 The latter are known to depend very sensitively on system parameters such as the temperature and therefore may constitute an additional diagnostic tool.144

Acknowledgment. This work was supported partially by the Center for Advanced Systems Understanding (CASUS), which is financed by Germany’s Federal Ministry of Education and Research (BMBF), and by the state government of Saxony from the State budget approved by the Saxon State Parliament. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2022 research and innovation program (Grant No. 101076233, “PREXTREME”). The PIMC calculations were carried out at the Norddeutscher Verbund für Hoch-und Höchstleistungsrechnen (HLRN) under Grant No. shp00026, and on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden.