High statistics lattice study of stress tensor correlators in pure $SU(3)$ gauge theory

We compute the Euclidean correlators of the stress tensor in pure $SU(3)$ Yang-Mills theory at finite temperature at zero and finite spatial momenta with lattice simulations. We perform continuum extrapolations using $N_\tau=10,12,16,20$ lattices with renormalized anisotropy 2. We use these correlators to estimate the shear viscosity of the gluon plasma in the deconfined phase. For $T=1.5T_c$ we obtain $\eta/s=0.17(2)$.


I. INTRODUCTION
Since relativistic hydrodynamics is quite successful in the interpretation of heavy ion experiments [1][2][3][4][5], it would be of great interest to calculate the shear viscosity of the quark gluon plasma from first principles.
In classical transport theory, the shear viscosity to entropy density ratio for a dilute gas at temperature T is η=s ∼ Tl mfpv ∼ Tv nσ , where l mfp is the mean free path,v is the mean speed, n is the particle number density and σ the cross section. For a weakly interacting system, σ is small, and η=s is expected to be large. In particular, for a free gas, η=s is infinite. On the other hand, for a strongly interacting system, η=s is expected to be small [6]. As heavy ion phenomenology points to a rather small viscosity [1][2][3][4][5], a nonperturbative calculation of the shear viscosity would be a great success.
One possible route to determine the viscosity is through the Kubo formula, relating transport coefficients to the zero-frequency behavior of spectral functions. The relevant Kubo formula for the shear viscosity is where ρ ijij ðω; k; TÞ is the spectral function corresponding to the energy momentum tensor at the specified spatial indices i ≠ j. The direction of the momentum is j. In this paper, we will assume without any loss of generality, that the external momentum is in the 3rd direction, while the zeroth direction is the (Euclidean) time. By choosing a matching i index we will consider the component ρ 1313 ðω; k; TÞ.
In general, the correlator of the energy momentum tensor T μν is given in Euclidean spacetime as Both early [7][8][9] and more recent [10] lattice studies of the viscosity used the Kubo formula (1). In this approach the integral transform (3) has to be inverted. For T ≫ ω the kernel behaves like e −ωτ , i.e., our task is similar to inverting a Laplace transform numerically. It is well known that such an approach is bound to face great difficulties. There are two interrelated problems: (1) Equation (3) is a Fredholm equation of the first kind, which for most kernels very ill-posed. The difficulty can intuitively be compared to the process of deblurring an image. Also, in particular, both the Laplace kernel, and our kernel Kðω; τ; TÞ are known to lead to a very ill-conditioned inverse problem [11]. (2) For the particular case of the viscosity, the signal in the stress-energy tensor is strongly dominated by the high frequency part of the spectral function [12]. This makes reconstruction even harder as the blurring character of the integral transform (point 1) mixes the contributions from the high and low ω part of the spectral function in the measured Euclidean correlator. To see how bad a particular inversion problem is, it is very instructive to look at the spectral function for the free theory, as this will correspond to the asymptotic behavior of the spectral function in the continuum theory, because of asymptotic freedom. To get the asymptotic behavior up to a constant, one only has to perform simple dimensional analysis. For ρ 1313 this leads to an asymptotic ω 4 behavior, making the UV contamination especially severe.
To see an honest illustration of these problems for ρ 1313 ∼ ω 4 , look at Fig. 1, where we illustrate how insensitive the Euclidean correlator is to the IR features of the spectral function. There we show two different spectral functions, with a factor of 10 difference in the viscosity, that nevertheless lead to subpercent differences in the corresponding Euclidean correlators. The two mock spectral functions in Fig. 1 are actually both physically motivated. The featureless spectral function (#1 in Fig. 1) is reminiscent of the one obtained from calculations in N ¼ 4 SYM theory, with AdS=CFT methods [13], while the spectral function exhibiting a Lorentzian peak at ω ¼ 0 is reminiscent of the kind of results one obtains from leading log kinetic theory calculations in QCD itself [14]. Since the AdS=CFT calculation is a strong coupling calculation in the wrong theory, while the kinetic theory calculation is a calculation in the wrong regime of QCD, we do not know a priori which type of spectral function we can expect for QCD in the phenomenologically relevant temperature range, so a fully controlled calculation of the viscosity from the Kubo formula would necessarily need to distinguish between these two scenarios.
Note that for Fig. 1 we assumed that the asymptotic behavior of the spectral function is known completely accurately and there are no features of the spectral function at intermediate frequencies (e.g., no glueballs or remnants of melted glueballs). Though there has been progress in perturbative calculations of the UV part of the spectral functions [15][16][17], these assumptions are optimistic. Still, the fact that under these assumptions an order of magnitude difference in the viscosity leads to less than 1% difference in the Euclidean correlators nicely illustrates our point.
The bottom line of this discussion is that, for a credible lattice estimate of the viscosity, a high level of precision is necessary for the Euclidean correlators, especially if we want to use Eq. (1), as was done in Refs. [7][8][9][10]18].
The situation is much better for the correlators of conserved charges, appearing in the electric conductivity [19][20][21][22] and heavy quark diffusion [23][24][25][26] calculations. In those cases, the spectral function at large ω only grows like ω 2 , making the UV contamination problem less severe. It was an important realization of Refs. [11,27] that even for the case of the shear viscosity, the asymptotic ω 4 behavior can be made better, only ω 2 , by utilizing the following Ward identity, and using the ρ 0101 correlator, instead of the ρ 1313 . This would make the asymptotic behavior of the shear viscosity spectral function only as bad as that of the electric conductivity. But there are crucial differences as well. In the continuum, we have the thermodynamic identity [28] hT 01 T 01 iðτ; q ¼ 0Þ=T 5 ¼ s=T 3 . This means that we need nonzero momenta to obtain information about the viscosity from this correlator. Even with this knowledge, the calculation of the viscosity is still much more difficult than that of the electric conductivity. The source of the difficulty is the fact that the stress-energy tensor correlators CðτÞ have a quickly degrading signal as τ is increased beyond a few lattice spacings. Usually, the width of the distribution for these observables in a Monte Carlo simulation is much larger than the value, at least near the middle point τT ¼ 1=2, the very point where the correlator has its highest sensitivity to transport. Thus, the physically most relevant quantity is evaluated as an average of wildly fluctuating contributions (with fluctuating sign), which is typically the characteristic of a sign problem.  For the quenched case this problem can be ameliorated by using the multilevel algorithm [29,30]. This algorithm depends crucially on the locality of the action, and therefore it proved to be hard to generalize for dynamical fermions. Some progress in this regard has been made recently in [31,32]. Nevertheless, at least in the quenched case, high statistical precision can be achieved via the multilevel algorithm.
The study of cutoff effects of these correlators is rather limited in the literature. The tree level improvement coefficients for the plaquette action and two different discretizations of T μν where calculated in [33]. So far, no calculations of these correlators are available with three lattice spacings in the scaling regime.
In this paper, we take steps towards achieving the high precision necessary for the calculation of the shear viscosity, by investigating several technical aspects of such a calculation, namely, (i) Utilizing a different gauge action, the tree level Symanzik-improved action, as opposed to the plaquette action used in previous studies. (ii) Studing the continuum limit behavior by simulating at different values of the lattice spacing N t ¼ 10, 12, 16 and 20. (iii) Calculating the w 0 scale with high precision. (iv) Using the Wilson flow for anisotropy tuning, as advertised in [34]. (v) Using shifted boundary conditions for the renormalization of the energy momentum tensor, a technique that was worked out for the isotropic case in [35]. Here, we utilize it for an anisotropic lattice. (vi) Calculating the tree-level improvement coefficients for the Symanzik-improved gauge action. Throughout this paper, we will mostly focus on the calculation of the energy-momentum tensor, and not the inversion method for reconstructing the spectral function. We believe this to be an important first step. Before the inversion can be done, one needs to have reliable results for the correlator itself. Nevertheless, in the end we give an estimate of the viscosity, using a similar hydrodynamics motivated fit ansatz as some previous studies [11].

II. LATTICE CALCULATION OF THE CORRELATORS
Our calculation uses the tree-level Symanzik-improved gauge-action, where β ¼ 2N c g 2 .μ andν are the complementer indices for μ and ν, such thatμ <ν and the four indices μ; ν;μ;ν are a permutation of 0,1,2,3. Here W μν ðn; a; bÞ are Wilson loops around rectangular a × b paths. Finally, c 0 ¼ 5 3 and with ξ 0 the bare anisotropy. For our study we use anisotropic lattices with renormalized anisotropy ξ R ¼ 2. For anisotropy tuning we use the Wilson flow technique introduced in [34]. The procedure for anisotropy tuning will be detailed later.
We use the clover discretization of the energy momentum tensor, mainly because the center of the operator is always located on a site, therefore the separation of the operators is always an integer in lattice units. If one were to use the plaquette discretization, there would be a component that is defined for integer separations and one that is defined for half integer separations, and one would need an interpolation to add them together. This would lead to the appearance of a systematic error coming from the interpolation, that we want to avoid.
Following the line of previous studies, we use the twolevel algorithm, but now with a tree-level Symanzik improvement. [37] Thus we have thick layers (having a width of a full temporal lattice spacing) between the blocks in the inner update.

A. Statistics
As was explained in the Introduction, the correlators are needed to a very high precision, if one wants to have useful information on transport. In the pure SUð3Þ theory this can be achieved by using a multilevel algorithm [29] and high statistics. Earlier lattice studies of the viscosity [7][8][9][10]18] also use a multilevel algorithm.
The main difference is-apart from the higher statistical precision-that we are working with four different lattice spacings, which allows us to study the correlators in the continuum limit. Our statistics are summarized in Table I.   TABLE I. Number of measurements (millions) of the energymomentum tensor correlators at the simulation points. Between every measurement, there are 100 regular updates and 500 inner multilevel updates.
The relative statistical precision of our lattice data is shown in Table II.

B. Anisotropy tuning and scale setting
To fix the anisotropy we use the method introduced in [34]. The bare anisotropy ξ 0 ðβÞ is tuned so that ξ R ≡ 2. For the tuning, we define a spatial and a temporal w 0 scale, To tune the anisotropy, we use the following procedure: (1) We simulate the SUð3Þ theory at fixed β and several bare anisotropies around our estimate [34], targeting ξ R ¼ 2.
(2) We calculate the gradient flow using ξ R ¼ 2 and monitor w 0;x =w 0;t as a function of ξ 0 (see Fig. 2). The correct tuning of the anisotropy is achieved when w 0;x =w 0;t ¼ 1.
(3) The ξ 0 ðβÞ data set is fitted with a Padé formula. The fitted curve is plotted also in Fig. 2. Our parametrization reads We likewise fit w 0 ðβÞ, with a parametrization that interpolates smoothly between the two-loop running of the coupling and the lattice data, So far, we expressed the scale using w 0 . In order to be able to translate to the T c scale, we have to determine the combination w 0 T c in the continuum limit. We did this using four different (isotropic) actions (Wilson, tree-level Symanzik, Iwasaki and DBW2). With the exception of the last one, we found similar results using lattices up to N τ ¼ 12 and a continuum limit using an N 2 τ as well as an N 4 τ term. The uncontrolled systematics of the DBW2 result is no surprise, this action is known to poorly sample topological sectors, see for e.g., [39].
In all cases, T c was defined by the peak of the Polyakov loop susceptibility. The summary plot for this study is shown in Fig. 3 For each action, the resulting jackknife error was very small. Therefore we use the spread between the Wilson, Symanzik and Iwasaki results as an error estimate instead, and use the Symanzik result (which lies central between the others) as mean. We conclude that w 0 T c ¼ 0.2535ð2Þ.

C. Renormalization
The translational symmetry is broken on the lattice. As a result, renormalization factors appear between the lattice definition of the energy momentum tensor T μν and the physical quantity. This factor depends on the action, the discretization scheme in the T μν observable and the lattice spacing (or the β parameter). Moreover, these factors are not the same for each component, since the off-diagonal (sextet), the diagonal (triplet), and the trace (singlet) correspond to different representations of the fourdimensional rotation group. On an isotropic lattice, one has three factors, and there is no summation over μ and ν in the above formulas.
We use the clover discretization of F a μν and define our correlators from the sextet (off-diagonal) components. In the presence of anisotropy, the renormalization constant Z 6 splits into three different renormalization constants: In our renormalization procedure, we get Z ts 6 from the thermodynamic identity (20), and we get the ratios Z ss 6 =Z ts 6 and Z tt 6 =Z ts 6 from shifted boundary conditions. For an isotropic gauge action, the renormalization constants have been worked out with shifted boundary conditions in [35]. Using shifted boundary conditions with shift vector ⃗ ξ ¼ ðξ 1 ; ξ 2 ; ξ 3 Þ ¼ ð1; 1; 1Þ, the off-diagonal T 0i components develop a nonvanishing expectation value. Since with this particular choice of the shift, the three spatial directions are equivalent, we have T 01 ¼ T 02 ¼ T 03 . We note that this renormalization condition is very similar to the one used in [40], where it was used for the secondorder hydrodynamic coefficient κ. Imposing this condition gives 2Z tt Therefore, the ratios Z ss 6 =Z ts 6 and Z tt 6 =Z ts 6 can be calculated from a single simulation with L −1 Thus, e.g., to renormalize T μν in a N τ ¼ 12 simulation with ξ R ¼ 2, we make an auxiliary run on a 48 × 96 × 48 × 3 lattice with the same bare parameters. The resulting factors will depend on β and N τ . The method requires that N τ =4 is an integer. We observe a 1=N 2 τ scaling of both Z ss 6 =Z ts 6 and Z tt 6 =Z ts 6 . For the renormalization of N τ ¼ 10 we can therefore use an interpolation in N τ . Our simulated results on the renormalization factors can be seen in Fig. 4.
The overall constant Z ts 6 can be determined from the following thermodynamic identity [41]: This can be used for renormalization by requiring that the value of C 0101 at τT ¼ 0.5 equals the continuum value of the entropy determined in [42]. We used the values s=T 3 ¼ 5.02 and 5.57 for 1.5T c and 2T c respectively. This identity also provides a way to estimate the order of magnitude of the discretization errors in C 0101 . Since in the continuum this correlator is independent of τ, the τ dependence of the correlator gives a very direct way to see discretization errors already on the finite N τ data (for the case of N t ¼ 16, see Fig. 6).

III. RESULTS ON THE CORRELATORS
A. Results at finite N t The 13 channel correlators can be seen in Fig. 5, while the results for the 01 channel can be seen in Fig. 6.
We note here, that trying a simple featureless ansatz ρ 1313 ¼ Aω þ Bω 4 on our N t ¼ 20 lattices (T ¼ 1.5T c ) for τ ≥ 6=10 and q ¼ 0 we get χ 2 =n dof ¼ 3.8=3, while an ansatz with an infinitely sharp transport peak ρ 1313 ¼ AδðωÞω þ Bω 4 gives χ 2 =n dof ¼ 5.7=3, giving a slight preference to the featureless scenario that we will later use in our more detailed analysis, while not yet ruling the sharp transport peak scenario out.
For the 01 channel, in the continuum, the correlator for zero spatial momentum should be a constant, equal to the entropy. The renormalization condition we used for this correlator is simply that at the middle point, τT ¼ 1=2 it should equal the continuum value of −s=T 3 . How different the correlators value is for τT ≠ 1=2 is some kind of measure of the cutoff effects. As we already discussed, we expect the 01 channel to have smaller cutoff errors and also to be more sensitive to transport, so this is the more important of the two correlators.

B. Continuum limit extrapolation
For the purpose of this paper we focus our discussion of the continuum limit extrapolation to the middle point of the correlators τT ¼ 1=2. We choose this approach for several reasons: (i) This is the most IR sensitive part of the correlators, therefore the most interesting part for studying transport. (ii) This is the part of the correlator with the least amount of cutoff effects, therefore one has to control the continuum extrapolation of this first, before attempting to go to smaller separations in imaginary time. FIG. 5. The renormalized shear correlator C 1313 at different lattice spacings and different spatial momenta. We also present a continuum estimate, that was produced by performing a spline interpolation of the finite N t data. We only present the continuum estimate in the range where the χ 2 was acceptable.
FIG. 6. The renormalized shear correlator C 0101 at N t ¼ 16 and for different spatial momenta for the temperature T ¼ 1.5T c .
Notice, that the C 1313 correlator is closely related to the τ derivative of the C 0101 correlator: as can be seen from a differentiation of the sum rule (3) and application of the Ward identity (5) respectively. Thus, taking the C 1313 ðτ; qÞ and C 0101 ðτ; qÞ correlators only at the value τT ¼ 1=2 already contains the leading τ dependence of C 0101 . Thus, we may continue with the extrapolation at τT ¼ 1=2. We will attempt a continuum limit extrapolation both with and without tree-level improvement. The tree-level improvement coefficients are the result of a tedious, but straightforward computation. The numerical values of the improvement coefficients are summarized in Appendix B. We will also attempt both linear and quadratic fits for the continuum limit extrapolation. Attempting a continuum limit extrapolation from our N t ¼ 10, 12, 16, 20 data yields the following behavior: (i) In the 0101 channel, since one applies the renormalization condition (11) after the tree-level improvement, the continuum extrapolation is quite flat, regardless of whether one uses tree-level improvement or not. (ii) A linear fit to the N t ¼ 10, 12, 16, 20 lattices in the 1313 channel with and without tree-level improvement does not always yield consistent results within 1σ for the continuum limit extrapolation.
(iii) A quadratic fit to the N t ¼ 10, 12, 16, 20 lattices in the 1313 channel with and without tree-level improvement does yield consistent results, but then we have one degree of freedom less, so the error on the continuum is larger, roughly on the 2%-3% level. (iv) Linear versus quadratic fits to the data obtained without tree-level improvement are not consistent within 1σ for the 1313 channel. (v) Linear versus quadratic fits to the tree-level improved data are closer, but still not consistent within 1σ for the 1313 channel. This behavior can be visually observed in Figs. 7-9 where the linear and quadratic extrapolations are shown.
From this analysis, we conclude that from our present data, the continuum extrapolation has error bars on the few  Table III.

C. Finite volume effects at tree level
From the tree-level calculation, we can estimate the finite volume effects on the UV contribution to the correlators. For the volumes used for our simulations, i.e., L x T ¼ L y T ¼ 2 and L z T ¼ 8 we calculated the tree-level (UV) contribution of the spectral function in Appendix B. The relative deviation from the infinite volume contribution is shown in Table IV. Thus, the tree-level finite volume effect on the observables considered here is on the 10% level. While 10% error on the final viscosity is probably harmless at this point, it may shift the relative weight of the UV and IR contributions. It is important to note, that the finite volume correction depends very weakly on q and for C 0101 it weakly depends on τ, too. Thus, in the viscosity fits the volume dependence approximately factorizes, and it affects only the value of the c parameter in Eq. (24). The values we quote later will correspond to the raw fit, which is expected to be roughly 10% below the infinite volume value.

IV. ESTIMATING THE VISCOSITY
To get an educated guess on the viscosity, one needs to assume an ansatz. Here, we assume a very simple hydrodynamics plus tree-level ansatz for the spectral function, corresponding to the featureless scenario in Fig. 1: The hydrodynamic predictions for the spectral function are written out in Appendix A. The tree-level spectral function is taken to be the one at infinite volume and in the continuum. Notice that the integral in ω for this UV part is cut off at ω ¼ jqj in the IR. The part at lower ω is responsible for the hydrodynamical behaviour, which is taken into account by the other term. Actually, the free gas formula would give an infinite contribution to the viscosity. The formulas for the continuum tree-level spectral function can be found in Ref. [27] and are also summarized in Appendix B.
We introduce a constant c in front of the spectral function. We do so to account in a simple way for higher-order and also finite volume corrections. The assumption that all of these effects can be put into a single constant is a very strong one. One hint that it might be a good estimate is given by Table II, where we have the finite  volume correction factors for different correlators. This model has two free parameters [43], the factor c in front of the tree-level correlator, and the shear viscosity to entropy ratio η=s may be contaminated by higher-order contributions, as well. Clearly, our estimate of the viscosity is correct only in the range of validity of this simple model for the spectral function. In the following, we will work out how data constrain the model parameters. Assuming that our system is within the model's range of validity, we can make a quantitative statement on the shear viscosity.

A. Sensitivity to model parameters
Before describing the fitting procedure let us show how the different model parameters influence the observables considered here. Since it is the more interesting quantity, our discussion here will focus on C 0101 . Figures 10 and 11 concentrate on changing one of the parameters, c or η=s, respectively. From these pictures the conclusion one can draw is that C 0101 , while certainly sensitive to the hydrodynamic parameter η=s, is also sensitive to the UV parameter c. This is not surprising, but it is a big advantage compared to the C 1313 correlator, where the sensitivity to η=s is smaller. Still, one has to acknowledge that while this is a pretty useful quantity, because of the sensitivity to both parameters, it is not enough to constrain the value of c and η=s together. To do that one has to consider in addition the τ dependence of C 0101 , or equivalently, the correlator C 1313 at τ ¼ 0, which, as we have already shown, corresponds to the second τ derivative of C 0101 .

B. Fits at finite N t
We will present two different fits for the parameter c and η=s. First we study only C 0101 as a function of τ and q at N τ ¼ 16. This approach is similar to what was used in earlier publications, where only data at one finite N t were available. Our choice of the channel C 0101 is motivated by the smaller cutoff errors compared to C 1313 , as well as its higher sensitivity to the transport part of the spectral function. We constrain our fits to the range τT ∈ ½0.3; 0.5. This range is motivated by two observations: i) the quantity C 0101 ðτ; q ¼ 0Þ, which is, in principle, independent of τ in the continuum, is indeed constant for N t ¼ 16 in this range (see Fig. 6); ii) the finite volume corrections at tree level are also τ independent in this range. The latter fact motivates the assumption that most finite volume effects can be captured by a modified value of the c parameter.
T η=s c 1.5T c 0.178 (15) 0.60(6) 2.0T c 0.157 (13) 0.63 (7) Here the error includes statistical errors, as well as systematic errors coming from the choice of τ min ¼ 5.=16 or 6.=16 and q max ¼ 2π=4 or q max ¼ 3π=4. These numbers are, of course, only valid once we assume our hydrodynamic ansatz. The viscosity appears to be temperature independent from our analysis. Here we have to mention a serious drawback of our fit ansatz. It assumes that the hydrodynamic prediction for the spectral function, strictly valid only for ω ≪ T, is also a good approximation for higher frequencies. This is true for N ¼ 4 SYM theory, where AdS=CFT can be used to calculate the spectral function [13]. Our ansatz cannot produce a quasiparticle peak, that would appear in weak coupling treatements of QCD, like kinetic theory [14,44,45]. This means that the physical mechanism that makes the viscosity diverge for T → ∞, namely the sharpening of the peak in ρ 1313 near ω ¼ 0 is missing from our ansatz. This implies that our ansatz can certainly not be used at very high temperatures, where the weak coupling calculation is trustworthy, and even at intermediate temperatures we might underestimate the viscosity, effectively smearing out the transport peak by enforcing the ansatz in the data analysis. This is a weakness shared by all previous lattice estimates of the shear viscosity, since they either use a very similar hydrodynamic ansatz, or the Backus-Gilbert method, in which case the base functions used for the reconstruction similarly prefer a featureless behavior, as opposed to a transport peak [11].

C. Fits on continuum data
For the second fit we consider the q 3 =ðπT=4Þ ¼ 0, 1, 2, 3 dependence of C 0101 ðτT ¼ 0.5Þ and C 1313 ðτT ¼ 0.5Þ in the continuum. Our results are Here, the error is statistical only. The systematic error coming from the choice of the τ range does not exist here, since we always use just τ ¼ 1=2. This fit uses the same hydrodynamic ansatz as discussed above. This is the first estimate of η=s using continuum extrapolated data. It is consistent with earlier estimates using a single lattice spacing [9,46].

V. SUMMARY
In this paper we studied the continuum behavior of the energy-momentum tensor correlators in pure SUð3Þ gauge theory. We found cutoff errors at a few percent level for N t ¼ 16. For some quantities, namely C 1313 and C 0101 at several spatial momenta, and τT ¼ 1=2 continuum extrapolation was possible. Out of these quantities C 0101 is actually sensitive to the transport part of the spectral function.
The achieved percent level precision of the data does not yet allow us to distinguish different scenarios for the spectral function. The statistical precision was already boosted by the multilevel algorithm on an anisotropic lattice. Despite recent efforts of using the gradient flow to measure energy-momentum tensor correlator [47,48] and the promising extension of the multilevel algorithm to full QCD [31] it is hardly possible to achieve even this precision with dynamical quarks. Thus, for the testing of the hydrodynamical model one has to seek for alternative methods.
Once, however, a model is postulated, it is possible to give a model-dependent estimate of the shear viscosity to entropy ratio. We gave the first estimate of this phenomenologically important quantity from continuum extrapolated lattice data. Our estimate is in the same ballpark as earlier estimates based on finite N t lattices. The combination of linearized relativistic hydrodynamics and linear response theory allows for a derivation of the low frequency behavior of the energy-momentum tensor correlators. For a nice derivation of these formulas, see the Appendix of [13]. Here, we just collect the relevant formulas in our notation, for easy reference. We assume the spatial momentum is in the z direction k ¼ ð0; 0; kÞ. In this case the spectral functions in the shear channel are:

ACKNOWLEDGMENTS
where s is the entropy, η is the shear viscosity and T is the temperature. We use these formulas for both our finite N τ and continuum data fits. The zero spatial momentum limit of ρ 1313 =ω is a constant equal to η=π, while the zero spatial momentum limit of the ρ 0101 =ω is a delta function at the origin: as can be easily shown using Eq. (A1). This is the hydrodynamic identity we utilize for our renormalization procedure. Formulas (A1) and (A2) are also the basis for the derivation of the Kubo formulas, like Eq. (1).

APPENDIX B: TREE-LEVEL SPECTRAL FUNCTION IN THE CONTINUUM
The leading-order perturbative result for the spectral function at high frequency is [27] −ρ The tree-level result in the ρ pert 1313 channel follows trivially using the Ward identity (5).
In our analysis, we only take the first part ω > q. The reason for that is that the ω < q part describes the transport properties of a free gas of gluons, and corresponds to an infinite viscosity. We therefore drop this term and substitute it with the ansatz from hydrodynamics (which describes a strongly coupled system).

APPENDIX C: TREE-LEVEL IMPROVEMENT COEFFICIENTS
The formulas for the tree-level improvement are the result of a tedious but straightforward leading-order calculation. The resulting formulas still contain Matsubara sums, that can be easily evaluated numerically. For reference, we include here the numerical values of the tree-level improvement coefficients relevant for our study.