An accurate theoretical description of the dynamic properties of correlated quantum many-body systems, such as the dynamic structure factor *S*(**q**, *ω*), is important in many fields. Unfortunately, highly accurate quantum Monte Carlo methods are usually restricted to the imaginary time domain, and the analytic continuation of the imaginary-time density–density correlation function *F*(**q**, *τ*) to real frequencies is a notoriously hard problem. Here, it is argued that often no such analytic continuation is required because by definition, *F*(**q**, *τ*) contains the same physical information as does *S*(**q**, *ω*), only represented unfamiliarly. Specifically, it is shown how one can directly extract key information such as the temperature or quasi-particle excitation energies from the *τ* domain, which is highly relevant for equation-of-state measurements of matter under extreme conditions [T. Dornheim *et al.*, Nat. Commun. **13**, 7911 (2022)]. As a practical example, *ab initio* path-integral Monte Carlo results for the uniform electron gas (UEG) are considered, and it is shown that even nontrivial processes such as the roton feature of the UEG at low density [T. Dornheim *et al.*, Commun. Phys. **5**, 304 (2022)] are manifested straightforwardly in *F*(**q**, *τ*). A comprehensive overview is given of various useful properties of *F*(**q**, *τ*) and how it relates to the usual dynamic structure factor. In fact, working directly in the *τ* domain is advantageous for many reasons and opens up multiple avenues for future applications.

## I. INTRODUCTION

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 excitations^{2–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 temperature^{13} and in the ground state.

^{14}Such

*linear-response*properties are readily measured in scattering experiments

^{15}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)

^{4}

^{,}

**q**is defined as the correlation of two density operators in reciprocal space, i.e.,

^{5}and normal quantum liquids,

^{16–18}including the distinct roton feature

^{19,20}present at intermediate wave numbers. A second important example is the diagnosis of warm dense matter (WDM), an extreme state

^{21–23}that occurs naturally in astrophysical objects such as the interiors of giant planets

^{24}and is important for technological applications such as inertial confinement fusion.

^{25}Here, x-ray Thomson scattering (XRTS) experiments

^{26,27}are a highly important diagnostic method and give insights into important system parameters such as the temperature

*T*

^{28,29}or the charge state

*Z*.

*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 formulations

^{30}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) approach

^{9,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*

_{b}

*T*). 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.,

*analytic continuation*(AC),

^{32}i.e., the numerical inversion of $LS(q,\omega )$ 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 methods

^{33–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**, *τ*).

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

*S*(

**q**,

*ω*) is generally impossible because of numerical instabilities.

*τ*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 find

*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 decomposition

^{48}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 model^{51,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 method^{9,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.

## II. THEORY

### A. Uniform electron gas

^{51,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

*ξ*

_{M}being the usual Madelung constant, and the Ewald pair potential $\varphi E(r\u0302l,r\u0302k)$ 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\u0304/aB$, with $r\u0304$ 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-classical^{62} (Θ ≫ 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 properties^{63–68} based on QMC simulations both in the ground state^{69–71} and at finite temperature^{72–76} has been fundamentally important for a wide spectrum of applications, most notably as input for DFT simulations of real materials.

### B. Path-integral Monte Carlo and imaginary-time correlation functions

^{4}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 gives

*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}(

*i*∈ {

*↑*,

*↓*}) of the respective permutation group $SNi$, and the corresponding permutation operators $\pi \u0302\sigma i$. As can be seen, the density operator $\rho \u0302=e\u2212\beta H\u0302$ in Eq. (7) can be interpreted straightforwardly as a propagation in imaginary time by an interval of

*t*= −

*iβ*.

*ϵ*=

*β*/

*P*, leading to

*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

*primitive*

*factorization*

^{80}

^{,}

^{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.,

*P*imaginary-time slices, and it holds that

**R**

_{0}=

**R**

_{P}. For simplicity, we assume that the summation over all permutations as well as the combinatorial normalization factors are contained implicitly in d

**X**. The weight function

*W*(

**X**) can now be evaluated explicitly; see the second line in Eq. (11). More specifically, the potential term is given by

*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

*l*th particle on imaginary-time slice

*α*. The kinetic term is off-diagonal and corresponds to the density matrix of a noninteracting system, i.e.,

*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 $\sigma \u03f5=\lambda \u03f5/2\pi $, with

*λ*

_{β}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 version

^{85}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.

^{33,82,96}In particular, the imaginary-time version of the intermediate scattering function is defined as the ITCF, i.e.,

*F*(

**q**,

*τ*) at integer multiples of the factorization step

*ϵ*,

^{82}i.e.,

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*^{−1}**n** and **n** ≠ **0** being an integer vector.

### C. Connection to dynamic structure factor and linear-response theory

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.

^{4}

^{,}

*χ*(

**q**,

*ω*) known from linear-response theory. It is convenient to express the latter as

^{100}

^{,}

*χ*

_{0}(

**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**,

*τ*).

^{111}see also the Appendix herein for a short derivation.

### D. Properties of imaginary-time intermediate scattering function

*detailed balance*relation of the DSF,

^{4,26}i.e.,

^{28}

^{,}

*F*(

**q**,

*τ*) is symmetric about

*τ*=

*β*/2. Therefore, any knowledge about the ITCF—e.g., from an XRTS measurement

^{26–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.

^{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 relation

^{4}

^{,}

*τ*, which gives

*F*(

**q**,

*τ*) at the origin give straightforward access to both odd and even frequency moments,

^{50,113}i.e.,

*α*≥ 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(\alpha )$, 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).

^{4},

*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\u0302(q)$. In addition, $Pm=e\u2212Em\beta /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.,

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

## III. RESULTS

### A. General trends: Synthetic spectra

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

*σ*. The results are shown in Fig. 3, with the top (resp. bottom) panel showing the results for the DSF with

*β*= 1 and

*ω*

_{0}= 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.,

*ω*

_{QP}being the QP energy, because Eq. (29) is a nascent delta-function. It is straightforward to show that the corresponding ITCF is given by

*σ*= 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 numbers^{119} that are sometimes used to diagnose plasma parameters^{26,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 *ω*_{0} 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.

*τ*≤

*β*/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 UEG

^{19}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.,

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 beryllium^{123} by Schörner *et al.*^{124}

### B. Uniform electron gas

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 factor* *F*(**q**, *β*/2), where *F*(**q**, *τ*) attains its minimum value for all *q*. All three curves exhibit a qualitatively similar peak around *q* = 2*q*_{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}.

*τ*decay as

*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*= 2

*r*

_{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* ∼ 2*q*_{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.

## IV. CONCLUSIONS

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 DSF^{19} 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 temperature^{28,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 objects^{22} 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 DFT^{30} 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 UEG^{135,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 calculations^{134} 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}

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Tobias Dornheim**: Conceptualization (equal); Writing – original draft (equal). **Zhandos A. Moldabekov**: Conceptualization (equal); Writing – original draft (equal). **Panagiotis Tolias**: Conceptualization (equal); Writing – original draft (equal). **Maximilian Böhme**: Conceptualization (equal); Writing – original draft (equal). **Jan Vorberger**: Conceptualization (equal); Writing – original draft (equal).

## DATA AVAILABILITY

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

### APPENDIX: IMAGINARY-TIME VERSION OF FLUCTUATION–DISSIPATION THEOREM

Even though the imaginary-time FDT has found applications in the literature, to our knowledge no derivation has ever been reported. The key is selecting an appropriate imaginary-time integration interval that will allow the application of the Kramers–Kronig relation.

*τ*∈ [0,

*β*] exposes the

*ω*= 0 pole, i.e.,

## REFERENCES

*Many-Particle Physics*

*Physics of Solids and Liquids*

*Quantum Theory of the Electron Liquid*

*Bose-Einstein Condensation*

*Bose-Einstein Condensation*

*Quantum Monte Carlo: Origins, Development, Applications*

*Fundamentals of Many-Body Physics: Principles and Methods*

^{4}He

^{3}He without fixed nodes

*High-Energy-Density Physics: Foundation of Inertial Fusion and Experimental Astrophysics*

*Graduate Texts in Physics*

*Frontiers and Challenges in Warm Dense Matter*

*Fundamentals of Time-Dependent Density Functional Theory*

*Lecture Notes in Physics*

^{4}He. Path integrals and maximum entropy

*Ab initio*low-energy dynamics of superfluid and solid

^{4}He

^{4}He from quantum Monte Carlo: Maximum entropy revisited

*Ab initio*path integral Monte Carlo results for the dynamic structure factor of correlated electrons: From the electron liquid to warm dense matter

*Ab initio*path integral Monte Carlo approach to the static and dynamic density response of the uniform electron gas

*Plasma Scattering of Electromagnetic Radiation: Theory and Measurement Techniques*

*Ab initio*path integral Monte Carlo simulations and analytical theory

*g*(

*q*) and the exchange-correlation kernel

*K*

_{xc}(

*r*) of the homogeneous electron gas

*Ab initio*exchange-correlation free energy of the uniform electron gas at warm dense matter conditions

*ab initio*path integral Monte Carlo study and machine learning representation

*Ab initio*quantum Monte Carlo simulation of the warm dense electron gas in the thermodynamic limit

*Ab initio*thermodynamic results for the degenerate electron gas at finite temperature

^{4}

*Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets*

*EBL-Schweitzer*

*Ab initio*path integral Monte Carlo simulations of the warm dense electron gas

*Ab initio*path integral Monte Carlo approach to the momentum distribution of the uniform electron gas at finite temperature without fixed nodes

^{4}

*B*

_{1g}Raman spectrum of the two-dimensional Heisenberg model

*The Method of Moments and Its Applications in Plasma Physics*

*Interacting Electrons*

*Ab initio*path integral Monte Carlo simulations and dielectric theories

*ab initio*path integral Monte Carlo simulations

*Ab initio*path integral Monte Carlo simulations

*Ab initio*dielectric response function of diamond and other relevant high pressure phases of carbon