Improved real-time dynamics from imaginary frequency lattice simulations

The computation of real-time properties, such as transport coefficients or bound state spectra of strongly interacting quantum fields in thermal equilibrium is a pressing matter. Since the sign problem prevents a direct evaluation of these quantities, lattice data needs to be analytically continued from the Euclidean domain of the simulation to Minkowski time, in general an ill-posed inverse problem. Here we report on a novel approach to improve the determination of real-time information in the form of spectral functions by setting up a simulation prescription in imaginary frequencies. By carefully distinguishing between initial conditions and quantum dynamics one obtains access to correlation functions also outside the conventional Matsubara frequencies. In particular the range between $\omega_0$ and $\omega_1=2\pi T$, which is most relevant for the inverse problem may be more highly resolved. In combination with the fact that in imaginary frequencies the kernel of the inverse problem is not an exponential but only a rational function we observe significant improvements in the reconstruction of spectral functions, demonstrated in a simple 0+1 dimensional scalar field theory toy model.

The underlying difficulties are related to the fact that lattice computations are carried out in an artificial Euclidean time. This in turn necessitates an analytic continuation of simulated correlation functions, before genuine real-time information can be extracted. It is in principle possible to carry out this task, since both Euclidean and Minkowski correlators are described by the same real-valued spectral function, to which they are connected via integral relations. In turn the real-time quantities of interest may often be expressed directly in term of these spectral functions, in case of the shear viscosity η = − 1 π lim μ→0 + lim p→0 + ρ ππ (μ,p) μ with ρ ππ (μ) denoting the spectral function of the spatial part of the traceless energy momentum tensor correlator in real-time frequencies μ. In practice the task of extracting spectral functions is highly ill-posed, as a continuous functions needs to be unfolded from a noisy and finite set of simulated data. Bayesian inference has been found to provide a highly efficient tool for carrying out such unfolding tasks and in this study we will deploy the recently developed BR method [14] whenever spectra are reconstructed.The standard relations for bosonic fields read Two concrete challenges manifest themselves in the above equations. First, the inverse problem from Euclidean data will suffer from exponential information loss induced by the strongly damping cosh/sinh kernel. This issue was well illustrated in [6], in that the Euclidean correlator is only marginally sensitive to sizable changes in the spectral function, especially at small frequencies. One possible step towards resolving this problem would be to go to imaginary frequencies, where the Källén-Lehmann kernel is only a rational function and thus the correlator exhibits much more sensitivity to spectral features, also at small μ. This however is only a partial solution since another more fundamental issue persists. In the standard Euclidean formulation of thermal field theory, the imaginary time axis is compactified and its physical extent related to the inverse temperature of the system. Thus the imaginary frequency correlator is resolved only on the Matsubara frequencies ω n = 2πT n, even though the right expression of eq.(1) indicates that G(iω) exists also at intermediate frequencies.
The physics of transport coefficients is indeed most prominently encoded in G(iω) between ω 0 = 0 and ω 1 = 2πT and already at ω 2 only contributes marginally to the correlator. Resolving the Euclidean domain more and more finely, i.e. going to the continuum limit, only gives access to higher Matsubara frequencies but does not resolve the regime relevant in practice any better. I.e. from the point of view of the inverse problem it would be most advantageous to resolve this low imaginary frequency regime, which then significantly improves the chance of a successful reconstruction. Of course if all Matsubara frequencies were known exactly, analyticity permits us to extract the full real-time information already from the discrete Matsubara information. In practice however both the number and precision of the data is limited. In order to make progress in the determination of transport properties, we have recently devised a novel approach to thermal quantum fields [15], whose ability to dramatically improve the reconstruction of spectral functions is demonstrated in a scalar toy-model in (0+1) dimensions.

The non-compact Euclidean time setup
Starting point for the novel approach is to formulate the thermal field theory as genuine initial value problem on the Schwinger-Keldysh contour, with a forward propagating branch populated by fields ϕ + (t), a backward propagating branch ϕ − (t) and a compact Euclidean domain of ϕ E (τ), which encodes the thermal initial conditions. The corresponding partition function reads where implicitly the fields at the latest real-time step t f are identified ϕ + (t f ) = ϕ − (t f ). Note that here the Euclidean timeτ is nothing but a mathematical tool introduced to efficiently sample the thermal initial conditions for ϕ + and ϕ − provided by the density matrix ρ = exp[−βH]. What makes thermal equilibrium special is the fact that the spectral function of the system can be computed directly from the fields ϕ + in form of the correlator G ++ (t) = �ϕ + (t)ϕ + (0)�, with the following momentum space decomposition before and after analytic continuation to imaginary frequencies I.e. in contrast to non-equilibrium settings no cross correlations between ϕ + and ϕ − are needed. In addition due to time translational invariance we may move the initial time point t 0 =τ 0 to negative infinity and extend the contour to positive infinity. At this stage we argue that, due to causality, no field ϕ + at finite time t may know about the implicit identification of ϕ + and ϕ − , now residing at t → ∞, since all correlations damp away. In turn we consider the two real-time paths effectively cut open as indicated by the dots on the right in Fig.1. The second step is to also change the way how the initial conditions are handled. Since the physics is time translational invariant, we can measure any real-time correlation functions of interest by inserting operators on the forward contour close to positive infinity, which again due to causality cannot influence the behavior at the starting point. Therefore we treat the Euclidean contour of ϕ E (τ) as genuinely compact and identify ϕ E (0) = ϕ E (β), which for a finite extent of the real-time path would of course not be permissible.
The final step consists of analytically continuing the time variable t of ϕ + and ϕ − onto a new infinitely extended imaginary time axis. As we regard the Euclidean timeτ as a mathematical tool only, we do still have the freedom to continue to a new τ = it. Due to the decoupling of the two fields at infinity, from now on we consider only the ϕ + and ϕ E degrees of freedom. With the extent of the imaginary time τ not being constrained, we may actually carry out a genuine Fourier transform and consider the system directly in imaginary frequencies, which now are not restricted to the conventional Matsubara ones.
Since several unconventional steps were taken in constructing this approach, let us have a look at the outcome in a non-trivial low-dimensional scalar toy-model, i.e. ϕ 4 theory in (0+1) dimensions. This theory is equivalent to the quantum mechanical anharmonic oscillator and it may be solved by computing an appropriate Schrödinger equation, in turn providing us with independent data for correlation-and spectral functions.
The computation is implemented via stochastic quantization on two finite grids with the same spacing Δτ, concurrently updating both the fields ϕ E , encoding temperature, as well as the fields ϕ + , from which the correlation function and spectrum is eventually determined. Based on the Euclidean we set up real-valued Langevin updates for both the compact and non-compact time domains, denote Monte-Carlo time as s and use independent Gaussian noises η and η � Note that of course in a finite box we are not able to represent an infinitely extended imaginary time axis. I.e. we will incur a finite volume artifact from having to specify an artificial boundary condition at the final τ. This issue already affects the free theory, due to the presence of time derivatives in S E . Putting a boundary by hand (periodic, Dirichlet, etc.) on a finite imaginary time path for ϕ + yields incorrect values for G ++ . They however can be made to converge to the correct result by extending the τ domain further and further, as we have explicitly tested. I.e. even by simulating directly in the new imaginary time one can systematically approach the correct result, paying as a price the numerical costs related to a large imaginary time extent. The convergence to the infinite volume result may be significantly improved however by formulating the action for ϕ + directly in continuous imaginary frequencies iω and subsequently discretizing it there. The free part of the action then is diagonal in imaginary frequencies and the correct (temperature independent) noninteracting G ++ is obtained trivially when computed via (still real-valued) Langevin stochastic quantization. No boundary conditions in iω need to be specified in the free case. Now in the presence of interactions, the ϕ 4 term is not diagonal in imaginary frequencies but is diagonal in imaginary time. Therefore we decide to implement the Langevin update in two steps, first computing the drift from the kinetic term of the action in imaginary frequencies, where it is local, then transforming the fields via discrete Fourier transform to non-compact imaginary time, where ϕ 3 is evaluated locally before transforming back and evaluating the two drift contributions together as in �η � (iω, s)� = 0 with �η � (iω)η � (ω � )� = 2δ(ω + ω � ) denotes Gaussian noise. While no explicit boundary conditions need to be specified in this mixed representation setting, the discrete Fourier transform of course pre-supposes periodic boundary conditions. This is why there still remains a residual finite volume artifact, which however is significantly smaller as compared to computing both the kinetic and interacting drift term directly in imaginary time.
We still have to answer how information about temperature is relayed from ϕ E to ϕ + . Since we know that in the free theory G ++ is temperature independent we cannot simply replace ϕ + (t = 0) = ϕ E (τ = 0) at each step from the concurrent update of ϕ E . Information about temperature only enters into G ++ in the interacting theory, which is why we implement the above mentioned replacement of ϕ + when evaluating the interacting drift term. This can be seen as an additional constraint on the drift, which is provided by the finite temperature initial density matrix.

Numerical results
As a first benchmark we check how well spectral properties of the theory are reproduced within the standard compact Euclidean simulation approach using λ = 24, m = 1, β = 1 on N τ = 16 grid points. We collect the values of G E (τ) = �ϕ E (τ)ϕ E (0)� after each Δs = 1 and accumulate N conf = 10 6 samples, which leads to statistical errors of around ΔG/G ∼ 10 −3 . At the same time we carry out a quantum mechanical (QM) computation where the Schrödinger equation for the Euclidean G E (τ), as well as the real-time field correlator is solved in the truncated Hilbert-space of 32 energy eigenstates of the harmonic oscillator [16]. In Fig.2 we show in the top panel the correlator values from that quantum mechanical computation as black crosses. The simulated values from a standard Euclidean stochastic quantization are given by blue circles and the inset confirms that they agree within statistical errors. On the bottom, the QM spectral function from Fourier transforming the numerical solution of the real-time Schrödinger equation up to t/m = 1600 is shown as thick gray curve. The residual width due to the numerical Fourier transform is not relevant here, we focus on the peak position.
To extract the spectral functions from simulated correlators we deploy the BR method [14] with realtime frequencies, denoted here as μ ∈ [10 −3 , 200] discretized along N μ = 4000 points. Around the low lying peak structure a high resolution window with N HP μ = 1000 points is selected. The default model we use corresponds to the minimum bias choice of m(μ) = 1. For the Euclidean correlator we deploy the naive discretization K(μ,τ, T ) = cosh[μ(τ − /2T )]/sinh[μ/2T ] of the kernel introduced in eq.(1). The spectrum reconstructed in this way is shown in Fig.2 as the dark blue curve in the bottom panel (almost covered by the light blue curve). It is cleat that with relative errors of order 10 −3 and N τ = 16 the peak position is not yet well reproduced.
In standard simulations, since temperature fixes the Euclidean time extent, the only handle for improvement we have is to go to higher resolution in compact imaginary time, discretizing β = 1 with e.g. N τ = 32, 64 or 128 points. As we argued before this only gives access to higher lying Matsubara frequencies, which are not very sensitive to the physics of the lowest lying spectral peak. And indeed we find that as shown by the three lighter colored curves in the bottom panel of Fig.2, going up to N τ = 128 does not improve the situation. On the other hand, a higher resolution inτ means that the drift term in the Langevin update becomes larger and we have to reduce the Monte-Carlo step size Δs significantly to keep the numerics stable, making these finer computation more expensive with limited practical use.
We now turn to the novel simulations, where besides a standard compact Euclidean domain with N τ = 16, we also consider the fields ϕ + in imaginary frequencies resolved on N ω ≥ N τ points. The kinetic term of the discretized action reads As is shown in the bottom left panel of Fig.3, resolving G ++ at the same N ω = 16 imaginary frequencies (blue) as the Matsubara correlator (gray), leads to values, which agree among the two for all finite ω n > 0 within uncertainties.  Figure 4. BR method spectral reconstruction from standard Matsubara correlators (gray) and from the novel simulations (colored solid). While access only to Matsubara frequencies (gray, Euclidean) and (blue, new approach) does not yield satisfactory results, resolving between the ω n 's dramatically improves the reconstruction success.
Interestingly, close to iω = 0 G ++ shows a deviation, which, if N ω is increased beyond 32, settles down to an asymptotic value slightly above the Matsubara correlator. This difference turns out to be an important ingredient to a successful spectral reconstruction.
For the imaginary frequency spectral reconstruction in the following K(μ, ω) = 2μ μ 2 +P(ω) is used, which corresponds to the lattice analogue of the continuum kernel and leads to better convergence of the reconstruction algorithm than using the naive discretization of eq.(1). Fig.4 again shows the QM spectrum as thick gray line, and the result from the BR method reconstruction, based on the standard Matsubara correlator, as dashed light gray line. We see that similar to the result in the bottom panel of Fig.2, the peak position is not satisfactorily captured with N τ = 16 Matsubara frequencies. Since G ++ (iω) for N ω = 16 agrees with its Matsubara counterpart for all finite ω n the corresponding spectral reconstruction (blue solid) not surprisingly also misses the peak position. Interestingly, the small deviation at ω 0 has already nudged the peak in the direction of the correct result.
The strength of our novel simulation approach lies in its capability to resolve between the conventional Matsubara frequencies. As shown in the left panel of Fig.3 going from N ω = 16 to N ω = 512 eventually produces a collection of points, which appear to smoothly cover the whole imaginary frequency range. From the discussion of eq.(1) we expect that access to the correlator in between the first and second Matsubara frequency will significantly improve the spectral reconstruction success. And indeed as is shown by the colored curves in Fig.4, for this simulation parameter set, already twice the standard imaginary frequency resolution N ω = 32 moves the peak much closer to the correct position. Going to even larger numbers, such as N ω = 512 allows us to reproduce the peak position with a deviation of only around one per mille.
While these results are encouraging they have been obtained in a simple toy model only. Therefore it is paramount to extend the computations performed here to genuine field theories. We are currently working on the straight forward implementation of code both for real-valued scalar fields in (2+1) dimensions, as well as for complex scalar fields in (3+1) dimensions. While the focus in the former case lies in investigating bound state spectra, the latter will provide a useful and realistic testing ground for the computation of transport properties. Indeed for the charged scalar field the relevant energy momentum tensor correlator for the shear viscosity [17] can be written as (7) First preliminary results for a lightweight and purely exploratory setup of λ = 1, m 2 = 0.2 on 8 3 × N τ × N ω grids with N τ = 4 and N ω = 32 are shown in Fig.5. On the left we plot the field twopoint function D(p) = �ϕ + (p)ϕ + (p)� at vanishing spatial momentum. The outcome of the standard Euclidean simulation is denoted by the filled squares, while the novel imaginary frequency simulation produces the open circles. The correlator is not exactly symmetric, since the real-and imaginary part of the scalar field contribute independently to the positive and negative imaginary frequency sides. At least in the symmetric phase, in which the simulation is situated, symmetry emerges better and better with increased statistics. Interestingly, it seems that here the agreement between the different simulations is better, in particular we do not find a statistically significant deviation at ω 0 . There are Figure 5. Preliminary results for the (3+1)d complex scalar field in imaginary frequencies. We show (left) the field correlator D(iω, p = 0), as well as (right) the energy momentum tensor correlator relevant for shear viscosity. Note that the latter is not affected by an exponential decay of the signal to noise ratio but still shows significantly larger errors than the simple field two-point function. two possible explanations for this, either finite temperature effects are not as relevant for the chosen parameter set or due to the still sizable errorbars, we are just not yet able to resolve a difference that is nevertheless present. Both increasing statistics, as well as testing different parameter sets for higher temperatures is currently high priority work in progress.
As last item we take a look at the preliminary result for the energy momentum tensor on the right in Fig.5. This correlator is exactly symmetric, which is why only its positive frequencies are shown. As expected, this observable, while not suffering from an exponentially decaying signal in imaginary frequencies, still shows much more statistical variance than the simple field two-point functions. It is a pleasure to observe that even though the statistical errors are still sizable, also for this correlator we find agreement between the novel and the standard compact Euclidean simulations on the Matsubara frequencies. A more quantitative study of the correlator and a first spectral reconstruction will be attempted, once adequate statistics of the order of ΔD/D ∼ 10 −2 have been achieved.
In summary, we have presented a novel approach to treating thermal fields as an initial value problem in an additional non-compact Euclidean time. We have formulated a simulation prescription based on this picture, in which the fields ϕ + and ϕ E are concurrently sampled using a real-valued Langevin update. To mitigate the effects of finite volume artifacts related to the finite extent of the new imaginary time axis, we simulate ϕ + directly in imaginary frequencies, where the free theory is diagonalized. Temperature information is supplied to the the ϕ + fields by imprinting the information of ϕ E (τ = 0) onto the drift term arising from the interaction part of the action in the Langevin update. Our results in a (0+1) dimensional toy model show that resolving the correlator between the conventional Matsubara frequencies indeed dramatically improves the reconstruction of spectral features. It however will be important to further clarify the meaning of the deviations observed between the correlators produced in the novel and standard simulations around vanishing Matsubara frequencies, in particular since they appear to be important for obtaining the correct spectral reconstruction. Preliminary results for complex scalar fields in (3+1) dimensions are also encouraging and their maturation is work in progress. This work is supported by the DFG Collaborative Research Centre "SFB 1225 (ISOQUANT)"