Thermal dynamics on the lattice with exponentially improved accuracy

We present a novel simulation prescription for thermal quantum fields on a lattice that operates directly in imaginary frequency space. By distinguishing initial conditions from quantum dynamics it provides access to correlation functions also outside of the conventional Matsubara frequencies $\omega_n=2\pi n T$. In particular it resolves their frequency dependence between $\omega=0$ and $\omega_1=2\pi T$, where the thermal physics $\omega\sim T$ of e.g.~transport phenomena is dominantly encoded. Real-time spectral functions are related to these correlators via an integral transform with rational kernel, so their unfolding is exponentially improved compared to Euclidean simulations. We demonstrate this improvement within a $0+1$-dimensional scalar field theory and show that spectral features inaccessible in standard Euclidean simulations are quantitatively captured.

We present a novel simulation prescription for thermal quantum fields on a lattice that operates directly in imaginary frequency space. By distinguishing initial conditions from quantum dynamics it provides access to correlation functions also outside of the conventional Matsubara frequencies ωn = 2πnT . In particular it resolves their frequency dependence between ω = 0 and ω1 = 2πT , where the thermal physics ω ∼ T of e.g. transport phenomena is dominantly encoded. Real-time spectral functions are related to these correlators via an integral transform with rational kernel, so their unfolding is exponentially improved compared to Euclidean simulations.
We demonstrate this improvement within a 0 + 1-dimensional scalar field theory and show that spectral features inaccessible in standard Euclidean simulations are quantitatively captured.
Two conceptual challenges plague the necessary reconstruction of real-time correlation functions from those simulated in imaginary time: the crucial one is the finite extent of the Euclidean time axis with τ ∈ [0 , 1/T ]. It leads to discrete Matsubara frequencies ω n = 2πnT , which limits the resolution of imaginary frequencies. In turn, increasing the number of temporal points N τ simply increases the maximum frequency available while keeping the spacing fixed. Thus, a large N τ does not improve the access to the relevant regime ω ∼ T , where thermal physics is dominantly encoded in G(ω). The signal to noise ratio for thermal contributions actually decays rapidly above ω ≈ T and indeed, already G(ω 1 = 2πT ) may only be marginally relevant for thermal physics.
In order to resolve this conceptual issue, we present a simulation algorithm, which operates in imaginary frequency space with arbitrary resolution, given by the number of frequency points N ω . In particular this approach resolves frequencies between ω 0 and ω 1 .
The second issue concerns the unfolding of spectral functions ρ(µ) from imaginary time correlators G(τ ), or imaginary frequency correlators G(ω). In Euclidean time the inverse problem reads (1) Eq. (1) leads to an exponentially hard ill-posed problem [35]. In this work we simulate directly in imaginary frequency space, where numerical data and spectral functions are related via the rational Källén-Lehmann kernel, (2) Its polynomial decay significantly improves the success of the spectral reconstruction already for Mastubara correlators with ω = ω n . Such a rational kernel is moreover an analytic implementation of the idea behind the Backus-Gilbert [36] or Sumudu [37] approaches.
In summary, it is the combination of the rational kernel in (2) and the access to arbitrary frequencies that leads to exponentially improved accuracy in the spectral reconstruction, see Fig. 4 for the results. Now we discuss the lattice setup that implements the above ideas. The starting point is a real-time thermal field theory, formulated on the Schwinger-Keldysh contour [38,39]. It amounts to an initial value problem with the following path integral representation for the partition function with S E being the Euclidean-and S M the Minkowski space action. ϕ lives on the forward and backward branch, ϕ + (t, x) lives on the forward branch from the inital time t = t 0 = −∞ to t = ∞, and ϕ − (t, x) lives on the backward branch, see Now we use that the real-time contour extends to infinity, and any correlations between fields at finite t and t = +∞ vanish. This allows us to cut open the Keldysh contour at t = +∞, and to rotate the real-time contour to an imaginary time axis, not to be confused with the time axis of the intial condition path integral. This rotation is implemented such that the path integrals for the rotated forward and backward branches have bounded statistical measures: the forward branch is rotated to the upper complex half plane while the backward branch to the lower one (see Fig. 1).
Furthermore, in equilibrium the correlator G in (2) can be derived as the analytic continuation of the real-time correlator of the fields on the forward branch, Thus it suffices to compute the correlation functions on the forward contour, and we do not consider the backward contour any further here. The forward branch may now be treated using stochastic quantization [40] with a Euclidean action S E . Its initial condition is given by the Euclidean field, ϕ + (τ = t 0 ) = ϕ E (0), and the correlation function G(ω) can be computed for all ω ∈ R.
In practice the simulation is carried out on a finite temporal lattice with periodic boundary condition which introduces finite volume effects. We have checked that these lattice artefacts disappear in the infinite volume limit, but the convergence is very slow. A qualitatively enhanced convergence of both, infinite volume and continuum limit, can be achieved as follows: the Fourier transform to frequency space translates the initial condition to the constraint l ϕ + (ω l ) = ϕ E (τ 0 ), where |l| = 0, ..., N ω /2 , and ω l = 2πl/(N ω ∆τ ). Here ∆τ is the temporal lattice spacing on the standard Euclidean lattice. This constraint can be rewritten as one for the field at the highest frequency, ϕ + (ω Nω/2 ) = ϕ E (τ 0 ) − l ϕ + (ω l ). The kinetic term 1/2 l ϕ + (ω l )P(ω l )ϕ + (ω l ) with lattice dispersion P is local in frequency space. Firstly, the constraint is irrelevant in the continuum limit, as its statistical weight is vanishing due to ϕ + (ω Nω/2 )P(ω Nω/2 )ϕ + (ω Nω/2 ) → ∞. Secondly, in the infinite volume we have an infinite number of momentum modes. This entails that the weight of one specific momentum mode tends towards zero. In summary, the dependence of the kinetic term on the inital condition is a well-understood lattice artefact. In the spirit of improved lattice actions we can lift the Matsubara constraint for all terms in the lattice action that are local in frequency space. Here we apply this reasoning to S 0 , the quadratic part of the action. The efficient convergence of this approach is confirmed for the propagator, see Fig. 3.
The full theory is then formulated in a mixed representation in τ and ω. While S 0 is simulated in ω-space, S int = S − S 0 is directly implemented via a discrete Fourier transform in τ -space, since interactions are local on the sites. Approximating the convolutions in the interaction term via DFT introduces a finite volume artifact, which vanishes as N ω → ∞. The thermal initial condition ϕ + (τ 0 ) = ϕ E (τ 0 ) solely enters via this interaction term. This concludes the formulation and discussion of our lattice setup.
To crosscheck our computations we have also performed an equivalent quantum mechanical (QM) computation of the thermal real-time correlation functions G QM (t) ≡ x(t)x(0) using the anharmonic oscillator in  the truncated Hilbert space of 32 energy eigenfunctions of the harmonic oscillator [38]. The spectral function is obtained from a Fourier transform of Im[ x(t)x(0) ] along N t = 3.2 × 10 4 steps of ∆t = 0.05, i.e. it is computed only over a finite real-time extent. Due to the restricted Hilbert space and integration this computation underestimates the height of delta-peak like structures, present e.g. in the free spectrum and contains ringing.
In Fig. 2 we show the Euclidean correlators computed for the thermal scalar system at β = 1 and m = 1 along the standard compactified Euclidean axis for the free λ = 0 (top points) and the interacting case λ = 24 (bottom points). The black crosses denote the outcome of the QM computation and open symbols denote lattice results. For the free case we only carry out simulations at Nτ = 16, while at λ = 24 we produce data at four different Nτ = 16 . . . 128. The inset shows the relative difference between the interacting correlator at Nτ = 16 computed from the lattice and via quantum mechanics. For our choice of ds = 10 −4 they agree within the statistical errors of the lattice simulation. Note that increasing Nτ , i.e. a smaller ∆τ , leads to a larger drift term, so that at Nτ = 128 the ds needs to be reduced by two orders of magnitude to maintain an accurate outcome.
We focus on simulations with Nτ = 16, and inspect the corresponding correlation functions in imaginary frequency space given in the top panel of Fig. 3 for both the free (gray) and interacting (colored) case. We show only the relevant low frequency regime up to ω/m = 15 while the simulations reach up to ω 8 /m ≈ 50. As ϕ + (τ ) ∈ R the correlators are real and symmetric around ω = 0.
The black filled symbols denote the correlators obtained from simply Fourier transforming the fields along the compact Euclidean domain and thus are located at  Open symbols arise from our imaginary frequency simulation. The inset shows the low frequency regime at λ = 24, in particular that finite volume artifacts are present at ω0 in the Matsubara correlator, which are removed going to finer resolution in ω. (bottom) Relative difference standard Matsubara correlators and those from our direct simulation. Note the excellent agreement except for ω0 at λ = 24, reflecting finite volume artifacts in the standard Matsubara value.
the Matsubara frequencies ω n . The open symbols are obtained from the fields simulated directly along the imaginary frequency domain using between N ω = 16 . . . 512 points. In the free λ = 0 case our novel simulation yields a smooth interpolation between the Matsubara frequencies and furthermore as can be seen from the gray points in the lower panel of Fig. 3 agrees excellently with the values obtained from the Fourier transform of the compact Euclidean domain. Note that at λ = 0 the two concurrent simulations (5) and (7) are completely decoupled.
In the interacting case denoted by colored open symbols in Fig. 3, we again find a smooth interpolation between Matsubara frequencies and very good agreement with the standard Matsubara propagators for ω n>0 . Let Free λ=0  us inspect the regime close to the origin, given by the inset on the top of Fig. 3. Already at N ω = 16 the new simulation (dark blue) leads to a value at ω = 0 that differs from the conventional Matsubara one (black). Increasing N ω quickly lets the correlator converge significantly above the standard estimate. Since in a lattice simulation ω = 0 reflects the sum over all corresponding imaginary time extent, we infer that the standard simulation along the compactified Euclidean domain is not able to capture all the relevant late Euclidean-time physics that however can be vital to thermal processes. Its purpose is only to correctly sample ϕ(τ 0 ), which it does. Now we extract spectral functions from the correlators, with kernelsK(µ,τ , T ) = cosh[µ(τ − /2T )]/sinh[µ/2T ] and K(µ, ω) = 2µ µ 2 +P(ω) . The latter is the lattice analogue of the continuum Källén-Lehmann kernel, and its use leads to a more efficient convergence of the spectral reconstruction. Here we deploy a recent Bayesian reconstruction (BR) method [41] and discretize positive real-time frequencies µ ∈ [10 −3 , 200] with N µ = 4000 points including a high frequency window of N HP µ = 1000 around the lowest lying peaks. For minimal bias we use a flat default model m(µ) = 1. In the free case we know that only a single delta-peak at the mass of the field exists in the spectrum. The quantum mechanical computation given by the thick gray line in the left panel of Fig. 4 while being able to capture the position of the peak, gives an artificial finite width due to the underlying finite realtime extent Fourier transform. We have checked that it can be reduced systematically and therefore for visualization purposes multiply the curve here by a factor of 256. The reconstruction from the Nτ = 16 compact Euclidean domain (dark gray dashed) both gives a too large width and its position is still slightly off. A Fourier transform of the Euclidean fields and reconstructing along Matsub-ara frequencies already improves the results significantly.
The light gray dashed line shows a much reduced width and is positioned very close to unity. Since here the standard Matsubara data already leads to a very good result, the improvement from using the correlators of the new simulation only leads to a further decrease in the reconstructed width. Note that the width depends on the error of the data and can be improved by more statistics.
Turning to the interacting case λ = 24, we first carry out the reconstruction from Eucldiean data, using different resolutions alongτ . In the center panel of Fig. 4 the results are shown for Nτ = 16 . . . 128. We find that the reconstruction does not capture the position of the lowest lying peak in the spectrum and even more importantly does not show a systematic improvement with Nτ . The reason is that with larger Nτ higher ω n become accessible but contribute only marginally to the relevant physics.
If on the other hand we use data in imaginary frequency space the situation is very different as shown in the right panel of Fig. 4. Again simply using the standard Matsubara correlators already improves the width of the reconstruction while the peak position is still off (gray dashed). However the access to the correlator at frequencies between ω 0 and ω 1 allows us to systematically and significantly improve the reconstruction of the lowest lying peak, with N ω = 512 reproducing the peak position with high accuracy.
In summary we have shown that by carefully distinguishing initial conditions from quantum dynamics it is possible to devise a simulation prescription for thermal quantum fields along general imaginary frequencies.
In particular this approach gives access to the thermal physics between ω 0 and ω 1 . Real-time spectral functions are related to imaginary time correlators via an integral transform with a rational kernel, and in combi-nation their extraction becomes exponentially improved. The success of this approach has been demonstrated here with scalar fields. It can be straightforwardly extended to Abelian and non-Abelian gauge theories which is work in progress. Evidently, computations of phenomenologically relevant quantities such as transport coefficients and heavy quark spectral functions will be improved significantly with our simulation prescription.
We thank N. Christiansen and S. Flörchinger for discussions. A.R. acknowledges fruitful exchanges with P. de Forcrand and A. Kurkela at CERN and during the ECT* workshop on advances in transport. This work is supported by EMMI, the grants ERC-AdG-290623, BMBF 05P12VHCTG. It is part of and supported by the DFG Collaborative Research Centre "SFB 1225 (ISO-QUANT)".