A test on analytic continuation of thermal imaginary-time data

Some time ago, Cuniberti et al have proposed a novel method for analytically continuing thermal imaginary-time correlators to real time, which requires no model input and should be applicable with finite-precision data as well. Given that these assertions go against common wisdom, we report on a naive test of the method with an idealized example. We do encounter two problems, which we spell out in detail; this implies that systematic errors are difficult to quantify. On a more positive note, the method is simple to implement and allows for an empirical recipe by which a reasonable qualitative estimate for some transport coefficient may be obtained, if statistical errors of an ultraviolet-subtracted imaginary-time measurement can be reduced to roughly below the per mille level.


Introduction
It is a perennial problem that a reliable study of strong interactions at temperatures around a few hundred MeV requires non-perturbative methods, yet standard Monte Carlo techniques only work in Euclidean signature. This implies that many of the physically most interesting observables, such as transport coefficients or particle production rates (for reviews see e.g. refs. [1,2]), which are inherently Minkowskian in nature, are difficult to address.
A few years ago, a possible solution to this problem was proposed [3], through the use of a "Maximum Entropy Method" [4]. Unfortunately, despite many implementations and various related attempts (see e.g. refs. [5]- [11]; similar discussions continue also in the context of condensed matter physics problems, see e.g. ref. [12] and references therein), it remains a problem that there appears to be uncontrolled dependence on model input in these results. Therefore, it is perhaps worthwhile to look for alternatives as well.
Some time ago, Cuniberti et al [13] put forward a concrete suggestion which could help in this respect. An algorithm was provided which was shown to yield a correct analytic continuation in a specific limit; furthermore, it was suggested that the method might work even if there is a finite amount of data and the data points have error bars. Despite these attractive features, we are not aware of a previous numerical implementation of the algorithm. For the record, we wish to report one in this short note, even if quantitative success is somewhat marginal. Our hope is that some of our theoretical remarks are nevertheless of interest and that, on the qualitative level, the algorithm might turn out to be a useful addition to the existing tool kit.

Basic idea
The general philosophy of the method can be summarized as follows. We consider a periodic ("bosonic") Euclidean correlator, G(τ, ·), with G(τ + kβ, ·) = G(τ, ·) where k ∈ , β denotes the inverse temperature, and · denotes suppressed variables such as spatial momentum. The correlator should be analytic everywhere except at Re τ = 0+kβ; there it should still be continuous. We also assume the further property G(β − τ, ·) = G(τ, ·), 0 < τ < β, which was not imposed in ref. [13] but is satisfied by typical gauge-invariant current-current correlators measured in lattice QCD. Given the finiteness of β, the Fourier representation involves a discrete set of Matsubara frequencies, G(ω n , ·) ≡ β 0 dτ e iωnτ G(τ, ·), ω n = 2πnT , n ∈ , T ≡ β −1 . Making use of a general hypothesis about the asymptotic behaviour of various correlators in the Minkowskian time domain, namely that they show at most powerlike growth at infinity (physically relevant correlators are actually expected to vanish at infinity, see e.g. ref. [14]), a unique analytic continuation to a part of the complex plane can be shown to exist and can be constructed explicitly.
(v) Given J + , the spectral correlator reads , whereJ + denotes a complex conjugate. Noting that J + as produced by eq. (2) is real under our assumptions, the spectral function can finally be obtained as

Practical implementation
To apply the previous algorithm to a situation where only a finite number of data points are available, we adopt the following steps, with the same numbering as in section 3.
(ii) Defining the coefficients φ ℓ,n ≡ (−1) ℓ n! 2 F 1 (−ℓ, n + 1; 1; 2), n ≥ 0, which can be constructed from the recurrence relation the a ℓ can be obtained from The precise upper limit (be it N − 2 or e.g. N/2) has little importance, because φ ℓ,n turn out to decrease rapidly for n > n max , where in general n max ≪ N/2 (cf. below). (iii) According to ref. [13], the |a ℓ | decrease only up to some ℓ max , meaning that ℓmax ℓ=0 |a ℓ | 2 shows a plateau as a function of ℓ max , but eventually it diverges. Ref. [13] suggests choosing ℓ max from this plateau. We find, however, that in practice it is more useful to monitor the sum ℓmax ℓ=0 a ℓ ; the reason is discussed in connection with eq. (9) below. (iv) Choosing ℓ max according to some criterion, J + can be approximated through where, as usual, the L ℓ can be constructed from In general the asymptotic value, does not vanish. In contrast, physically relevant current-current correlators should vanish at infinite time separation [14]. It turns out that this supplementary information can be used to choose values of ℓ max , namely those for which J + (∞, ·) vanishes approximately, offering "windows of opportunity", in which the algorithm appears to perform reasonably well even with non-ideal data. (v) The spectral function can be obtained as where we have subtracted by hand any possible "remnant" J + (∞, ·). (Otherwise ρ(ω, ·) would diverge as ∼ 1/ω at small frequencies.) As the proofs in ref. [13] show, in the limit N → ∞ and vanishing errors these steps do yield the correct spectral function for current-current correlators of the considered type.

Problems
We now turn to two problems that limit the usefulness of the algorithm specified above. For simplicity the discussion will be carried out in the combined continuum and infinite-volume limit; a finite cutoff alleviates problem (i) but adds more structure to the spectral function and thereby renders problem (ii) worse. A finite volume can in principle also change the behaviour of Euclidean correlators in a physically interesting way (cf. ref. [16] and references therein), however we wish to ignore these effects for now. (Formally a smooth infinite-volume type shape could be obtained e.g. by considering a suitable Gaussian smoothing of ρ(ω, ·)/ω.) (i) In order for the recipe to apply, we must be able to compute the Fourier coefficientsG(ω n , ·); this implies that G(τ, ·) must be integrable around τ = 0 mod β. In fact, as mentioned above, in ref. [13] G(τ, ·) was even assumed to be continuous (and therefore finite) at τ = 0 mod β. In contrast, the correlation functions of composite operators that are relevant for the determination of transport coefficients or particle production rates in QCD diverge at small τ in the continuum limit.
In principle, a possible way around the problem might be to consider the temperature derivative of a correlator, rather than a correlator as such. This could work if the derivative is taken in fixed physical units. Sometimes, such differences are rather taken in scaled units, i.e. by subtracting values ofĜ(τ , ·) ≡ β p G(τ β, ·), 0 <τ < 1, at different β's, but such differences are continuous only if there are no scaling violations at small τ . In terms of a spectral function, the finiteness of G(0 + , ·) necessitates ρ(ω, ·) ≤ C/(ω ln 2 ω) at ω ≫ T , cf. eq. (12). The asymptotics of various spectral functions in this regime have been analyzed in ref. [17], and for vector current correlators the decay of the thermal part is indeed fast enough to satisfy the bound. (ii) As can be seen from eq. (6), the factor (−1) n implies a cancellation between Fourier modes in the construction of the a ℓ 's; when the a ℓ 's multiply the Laguerre polynomials in eq. (7), further cancellations take place, particularly at larget (cf. eq. (9)). This implies a substantial significance loss. Furthermore, as can be seen from eq. (5), for a fixed ℓ the coefficients φ ℓ,n grow very fast at small n, reaching a maximal value at n max ≈ 2(ℓ + 1). The maximal value, φ ℓ,nmax , grows with ℓ, exceeding 10 4 at ℓ = 21 (for n max = 6), 10 6 at ℓ = 40 (for n max = 9) and 10 8 at ℓ = 64 (for n max = 11). So, the Fourier coefficients around n ∼ n max would need to be determined with the corresponding relative accuracy in order to obtain a meaningful signal even after the cancellations. This is obviously a formidable challenge, particularly considering that problem (i) already requires a subtraction and inflicts an associated significance loss.

Test
Despite the problems of the previous section, we now wish to demonstrate that in principle the method may still work on the qualitative level. In order to achieve this, we assume that a suitable ultraviolet subtraction has been carried out, and that our data are quite precise. As a model, we take inspiration from a spectral function related to an "electric field" correlator yielding the momentum-diffusion coefficient of a heavy quark [18,19,20]. We assume a substantial positive intercept of ρ(ω, ·)/ω at zero frequency (cf. ref. [21]) and decreasing behaviour at large frequency, just fast enough to yield a continuous G(τ, ·): ,ω ≡ ω 2πT .
This spectral function is not positive-definite in order to reflect the fact that a suitable ultraviolet subtraction has been carried out, and because a negative ρ ∼ −T 4 /(ω ln 2 ω) is precisely the qualitative asymptotic behaviour found in perturbation theory [22] (the logarithm squared assumes that the gauge coupling is let to run with ω). 3 The coefficient C appears linearly in all steps so that, without loss of generality, we set C ≡ 4T 3 in the following, thereby normalizing ρ test /(ωT 3 ) to unity atω → 0. The Euclidean correlator is subsequently integrated numerically from for τ = kβ/N ≤ β/2. At each τ we add a random error from a Gaussian distribution of relative variance σ 2 , and then mirror G(τ, ·) to the whole interval through the symmetry in τ → β − τ . On this "data" the steps of section 4 are applied. Small variants of eq. (11) bring along little change, but the situation deteriorates rapidly if the structure is more peaked either in the ultraviolet (ω ≫ 2πT ) or in the infrared (ω ≪ 2πT ).
The behaviour of the coefficients a ℓ , in particular the "diagnostic" sum of eq. (9), is illustrated in fig. 2; the quality of recovering the spectral function using windows of opportunity deduced from fig. 2 is shown in fig. 3. Fig. 4 demonstrates the effect of statistical noise on the final results. The lessons we draw are the following: (i) The recovery is reasonable only when ℓmax ℓ=0 a ℓ ≈ 0. For good accuracy, the sum ℓmax ℓ=0 a ℓ shows a nearzero minimum; then ℓ max should be chosen close to the minimum (the function ρ(ω, ·)/ω is also extremal there). For poor accuracy it could happen that no clear minimum is seen, cf. σ = 10 −3 in fig. 2(left); then ℓ max should be chosen close to a point where ℓmax ℓ=0 a ℓ crosses zero. (ii) Other values of ℓ max lead in general to nonsensical results. (iii) Within a given "robust" near-zero minimum, the dependence on N and σ is quite mild. (iv) The recovery can be qualitatively improved only by increasing the accuracy so much that a second window opens up (cf. σ = 10 −8 in fig. 2). This is unlikely to be reached in practice and, in any case, the improvement is not that overwhelming (cf. fig. 3). (v) The statements above apply to any single random configuration. With a sample of them, ℓ max could be separately fixed for each configuration. Within our toy model, it requires an accuracy σ < ∼ 10 −4 to find a useful minimum for almost every configuration (cf. fig. 4). For σ = 10 −3 , typical configurations show no near-zero minimum (cf. fig. 2), but rare ones do and if it is possible to restrict the statistics to those and to carry out the averaging on the final ρ(ω, ·)/ω, then a rough estimate can still be obtained.
Although it goes beyond the scope of the present note to carry out a detailed investigation of issues related to  Fig. 3.ρ test /ω ≡ ρ test /(ωT 3 ) as a function ofω = ω/(2πT ), for various relative accuracies σ, as well as ℓmax chosen from "windows of opportunity" in which ℓmax ℓ=0â ℓ approximately vanishes, cf. fig. 2 (results are shown for one random configuration). The thick solid line is the correct (input) result. The case σ = 10 −3 with N = 48 shows a "failed" example: there is a minimum in ℓmax ℓ=0â ℓ but it is not near zero (cf. fig. 2).
statistical analysis, we note that, in general, the output function depends non-linearly on input data, because the value of ℓ max varies and affects significantly the result. Error estimation should therefore be carried out with e.g. jackknife or bootstrap methods, perhaps with blocked configurations (the effect of blocking has been shown to be beneficial in connection with the Maximum Entropy Method, see e.g. ref. [23]).

Conclusions
The algorithm of ref. [13] possesses a number of attractive features: it can be fully specified in a small number of explicit steps; it requires no priors; it does not necessitate a positive-definite spectral function; and it projects out the Matsubara zero-mode contribution whose handling has been considered a problem in certain contexts. Unfortunately, from a practical point of view, the algorithm of ref. [13] cannot be guaranteed to yield a quantitatively accurate analytic continuation of thermal imaginary-time data. In some sense, the situation is akin to the sign problem hampering simulations of QCD with a finite baryon number density: there are significant cancellations taking place, particularly if a spectral function at a small frequency ω ≪ 2πT needs to be determined. Also, short-distance divergences need to be subtracted from the Euclidean correlator G(τ, ·), which constitutes a significance loss of its own.
Nevertheless, we have demonstrated that in a lucky case with a structureless spectral function and precise data (with relative errors < 0.1% after the ultraviolet subtraction), already N > ∼ 20 data points may yield a qualitative reproduction of a transport coefficient (zero-frequency intercept of ρ(ω, ·)/ω). In general, it is difficult to estimate systematic errors, but if a clear near-zero minimum in ℓmax ℓ=0 a ℓ is found as a function of ℓ max , then it appears that a < ∼ 50% uncertainty can be expected. This could already be useful, given that current model-independent determinations of transport coefficients might contain errors of more than 100% [21].