Physical insights from imaginary-time density–density correlation functions

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.¹ 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.⁷ 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.⁸ In addition, the quantum Monte Carlo (QMC) paradigm,^9–11 which is computationally even more expensive, can even give exact results in many situations,¹² both at finite temperature¹³ and in the ground state.

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.⁴⁶ 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, τ).

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),¹⁹ 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.¹⁹ The paper concludes with a summary and outlook in Sec. IV.

From a physical perspective, the UEG can be characterized fully by three dimensionless parameters.⁵⁷ (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-classical⁶² (Θ ≫ 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.

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.

Other concepts such as approximate expressions for the plasmon location at finite wave numbers¹¹⁹ 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 ω₀ 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.

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;²⁸ 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¹²³ by Schörner et al.¹²⁴ 

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

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

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.¹⁹ 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 DSF¹⁹ 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,¹²⁸ 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²² and inertial confinement fusion applications.¹³¹ In this context, we also note the advent of modern free-electron x-ray laser facilities such as the new European XFEL,¹³² 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³⁰ is the absence of reliable external input for the dynamic XC kernel K_xc(q, ω)¹³³ beyond the local density approximation or generalized-gradient expansions.¹³⁴ 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.¹⁴ 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¹³⁴ 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.¹⁴ Another interesting route is an explicit time-dependent DFT propagation along imaginary time,¹³⁹ 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.⁸² 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.¹⁴⁰ 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.¹⁴⁴ 

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.

The authors have no conflicts to disclose.

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

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

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.