The second order hydrodynamic transport coefficient $\kappa$ for the gluon plasma from the lattice

The quark gluon plasma produced in heavy ion collisions behaves like an almost ideal fluid described by viscous hydrodynamics with a number of transport coefficients. The second order coefficient $\kappa$ is related to a Euclidean correlator of the energy-momentum tensor at vanishing frequency and low momentum. This allows for a lattice determination without maximum entropy methods or modelling, but the required lattice sizes represent a formidable challenge. We calculate $\kappa$ in leading order lattice perturbation theory and simulations on $120^3\times 6,8$ lattices with $a<0.1$ fm. In the temperature range $2T_\mathrm{c}-10T_\mathrm{c}$ we find $\kappa=0.36(15)T^2$. The error covers both a suitably rescaled AdS/CFT prediction as well as, remarkably, the result of leading order perturbation theory. This suggests that appropriate noise reduction methods on the lattice and NLO perturbative calculations could provide an accurate QCD prediction in the near future.


Introduction
One of the major findings of the experimental heavy ion programme [1][2][3][4] is that QCD matter at high temperatures and low densities behaves as a nearly ideal fluid with very low viscosity. This conclusion is based on the fact that experimental data are excellently described by relativistic hydrodynamics, with transport coefficients fitted to the data [5][6][7][8][9][10]. Unfortunately, theoretical predictions of transport coefficients from the fundamental theory QCD remain very difficult [11]. Up to a few times the transition temperature to the quark gluon plasma, the QCD coupling is not weak enough for perturbative methods to apply, which predict a less ideal fluid [12,13].
On the other hand, results in the opposite strong coupling limit can be obtained by AdS/CFT duality methods in certain supersymmetric models [14,15], but these do not correspond to QCD directly.
Unfortunately, lattice simulations of real time quantities are in general severely limited by the need for analytic continuation. Calculations of spectral functions on the lattice based on maximum entropy methods [16,17] or a model ansatz [18][19][20] require both functional input and high accuracy data to sufficiently constrain the results. An exception to this conceptual difficulty are the three second-order hydrodynamic coefficients κ, λ 3 , λ 4 [21][22][23][24][25], which can be related to Euclidean correlation functions through Kubo formulae. They have recently been computed to leading order in a weak coupling expansion in [26], where also possibilities for a lattice determination were discussed. The coefficients λ 3 , λ 4 are related to three-point functions, which are still too costly to numerically evaluate.
Here we present a first attempt to determine κ from the momentum expansion of a suitable two-point function in a lattice simulation. In order to approach the zero momentum limit, very large lattices are required, demanding an enormous numerical effort already in pure gauge theory. While the errors on our result are thus still too large to be satisfactory, our work demonstrates that the determination of the second order coefficients is possible without conceptual difficulties and should be improved in the future with appropriate noise reduction methods. Interestingly and in contrast to the first order transport coefficients, the lattice result for κ is within error bars compatible with the perturbative weak coupling result. It is also compatible with a suitably rescaled AdS/CFT result.
In section 2 we briefly summarise the relation between the transport coefficient κ and a Euclidean correlator of the energy-momentum tensor, in section 3 this is carried over to the lattice formulation, including a leading order perturbative evaluation and a discussion of renormalisation. Section 4 contains the numerical results of our simulations.

The transport coefficient κ
The definition of transport coefficients is based on a gradient expansion of the energymomentum tensor in relativistic hydrodynamics, but their respective values have to be determined from experiment or an underlying theory. In the case of the quarkgluon plasma this underlying theory is QCD. In this section we review the connection between the transport coefficient κ and a Euclidean correlator in QCD, which allows for a direct computation of κ without resort to maximum entropy methods or functional input for the spectral function.

Relativistic hydrodynamics
The basic quantity in relativistic hydrodynamics is the energy momentum tensor (for a review, see [27]), which can be decomposed into an ideal part T µν 0 and a dissipative part Π µν The ideal part is determined by the hydrodynamic degrees of freedom, wich are the energy density , pressure p, the fluid's four velocity u µ and the metric tensor g µν . Lorentz symmetry and the identifications T 00 (0) = , T 0i (0) = 0 and T ij (0) = p δ ij in the local rest frame restrict its form to The dissipative contribution consists of a traceless part π µν and a part with nonvanishing trace Π The former has been specified for a non-conformal fluid in a second order gradient expansion within N = 4 Super-Yang-Mills theory [25] π µν = −ησ µν + ητ π Dσ µν + ∇ · u 3 σ µν + κ R µν − 2u α u β R α µν β + . . . . (2.4) Besides the shear viscosity η and the relaxation time τ π , to second order the transport coefficient κ enters the expansion and couples to the symmetrized Riemann curvature tensor R, its contractions and the fluid's four velocity u µ . For explanations of ∇, σ µν , D and further terms in the expansion we refer to [25]. Note that even in flat spacetime the transport coefficient κ has a non-vanishing value [26,28].

Thermal field theory
For the computation of the transport coefficient κ from QCD a relation between its definition in relativistic hydrodynamics and thermal field theory is necessary. This can be achieved by considering the fluid's linear response to a metric perturbation [24] and establishes a connection between the transport coefficient κ and the retarded thermal correlator of the energy momentum tensor T ij in momentum space, The transport coefficient κ is identified as the leading low momentum coefficient at zero frequency with momentum aligned in z-direction, q = (0, 0, q 3 ) [24,28], While the retarded correlator is a real time quantity, it is related by analytic continuation to the Euclidean correlator with the discrete Matsubara frequencies ω n = n2πT , n ∈ Z. This can be seen by writing both correlators in their spectral representation Appropriate boundary conditions render the analytic continuation unique [29], For vanishing frequency ω = 0 this can be written [11] G R (ω = 0, q) = G E (ω = 0, q) + B. (2.13) The contact term B arises from the missing commutator in the definition of the Euclidean correlator (2.8) compared to its retarded analogue (2.5) and corresponds to the correlator evaluated at equal spacetime points, ∼ T 12 (0)T 12 (0). An investigation of this contact term B by an operator product expansions shows that it is momentum independent [30]. Hence equation (2.7) can be rewritten where we have absorbed the constant G(0) and the contact term B into a new constant The transport coefficient κ can now be obtained as the slope of the low momentum correlator G E (q 2 ), which provides a possibility for a direct determination using lattice QCD. This is in contrast to computations of the shear viscosity [19] or heat conductivity [17,20]. These are true dynamical quantities which cannot be related to Euclidean correlators without non-trivial analytic continuation. Their determination by lattice calculations thus requires additional input, e.g. an ansatz for the spectral function or the maximum entropy method.
So far the discussion was completely general. We now specify to Yang-Mills theory and its energy momentum tensor [11] T µν = θ µν + 1 4 δ µν θ, (2.15a) where F a µν corresponds to the field strength tensor. The term θ = β(g)/(2g)F a αβ F a αβ includes the renormalisation group function β(g) and corresponds to the trace anomaly, which is caused by breaking of scale invariance. Since the transport coefficient κ is defined in the shear channel, T 12 T 12 , it does not enter the computation. Equation (2.14) has been evaluated perturbatively in pure gluodynamics in the ideal gas limit, i.e. at vanishing coupling, with the result [26,28]

Computation of κ in lattice QCD
In this section we describe the discretisation of the action and the energy-momentum tensor and explain the need for renormalisation. Furthermore, we calculate κ in lattice perturbation theory and compare with the continuum result (2.16).

Lattice framework
We employ Wilson's Yang-Mills action on an anisotropic lattice with different lattice spacings in temporal and spatial direction, a σ and a τ , respectively, with lattice coupling β = 2N c /g 2 and plaquette variables U µν . The bare anisotropy ξ 0 gets renormalised to the actual anisotropy ξ = a σ /a τ , .
We take the numerical evaluation of the renormalisation factor from [31]. The scale is set for a specific value of the anisotropy, ξ = 2, by comparison of the string tension from the lattice √ σ L [32] to its experimental value √ σ exp = 440 MeV [33]. The spatial lattice spacing follows from As will be discussed in section 3.4, the discretised energy-momentum tensor requires multiplicative renormalisation due to the reduced translational invariance on the lattice. For this purpose it is favourable to express the correlator (2.8) in terms of diagonal elements instead of nondiagonal ones. This is achieved by rotating the lattice by π/4 in the plane of the corresponding channel, i.e. the (1, 2)-plane for q = (0, 0, q 3 ). As shown in appendix A the trace anomaly θ does not enter the transformed correlator, although it includes diagonal elementes of the energymomentum tensor [34], Additionally, temporal and spatial elements of the energy-momentum tensor require separate renormalisation factors Z τ and Z σ on an anisotropic lattice. The diagonal energy-momentum tensor elements in the clover discretisation read where the bare elements are given by The clover plaquette [35] consists of four ordinary plaquettes (see figure 1) and is given by In contrast to an implementation with simple plaquette terms [34] the clover version has reduced discretisation errors and an improved signal-to-noise ratio [36], cf. figure 2.

Relation of κ to the lattice correlator
In order to extract κ numerically from equation (2.14), we compute the Euclidean correlator of the energy-momentum tensor within the lattice framework and perform a Fourier transform to momentum space with vanishing frequency. The determination requires the momenta to be aligned orthogonally to the studied channel of the energy-momentum tensor, i.e. q = (0, 0, q 3 ) for T 12 . This is also the case for the corresponding Kubo formula [11]. Thus the correlator in momentum space is given by Additionally, we include the channels T 13 and T 23 with corresponding momenta in our analysis, since rotational invariance allows to average over all three channels. We need small momenta compared to temperature, which sets the relevant scale, i.e. q i /T < 1. With the discretised versions of temperature and momenta we have for the ratio on the lattice The temporal lattice extent N τ is fixed by the temperature and lattice spacing. In order to fit the transport coefficient κ to equation (2.14), we need at least three different momenta satisfying this constraint (3.11). Thus the simulation requires large spatial lattice extents N σ , which makes the calculation costly. This can be partly moderated by working with anisotropic lattices ξ > 1.

Lattice perturbation theory
In order to estimate lattice artefacts and check our numerics, we first compute the transport coefficient κ in lattice perturbation theory on a lattice with anisotropy ξ in the case of vanishing coupling (g = 0). Definitions of relevant quantities and intermediate results can be found in appendix B, for an overview see e.g. [37].
In the case of vanishing coupling the field strength tensor simplifies to where we replace the differential operator by the central difference and define a lattice spacing a µ , which is excluded from Einstein's sum convention (3.14) In lattice perturbation theory the dynamical variables are the gauge fields A µ and we can plug the energy-momentum tensor from equation (2.15b) together with the field strength tensor (3.12) into the correlator (2.8). Then sixteen terms of the generalised form have to be transformed to momentum space. After transforming the individual gauge fields A µ (x) to momentum space by (B.1b), we apply Wick's theorem using the free gauge field propagator (B.3). Because of translational invariance it is sufficient to consider y = 0 or C i 1 i 2 j 1 j 2 l 1 l 2 m 1 m 2 (x, 0), and we obtain with the lattice momenta q, k as defined in appendix B. Evaluating the correlator (2.8) and aligning the outer momentum to q = (0, 0, 0, q 3 ) we find for its Fourier transform We perform the finite Matsubara sums by the residue theorem using the formula [38] 1 and list the results for the individual terms in appendix C. As will be described in section 3.4 we subtract the temperature independent vacuum part to avoid ultraviolet divergences. The three-momentum integration can be performed after expanding the integrals around the continuum limit. This step extends the integration measure to infinite volume [− π /a, π /a] 3 → R 3 and produces correction terms in small lattice spacings a σ . Together with the expansion in small momenta q 3 the remaining integrals can be solved analytically and one finds for the different terms For fixed temperature T = (a τ N τ ) −1 we can rewrite the dependence on lattice spacings a τ and a σ as a dependence on the temporal lattice extent N τ and the anisotropy ξ = a σ /a τ . Combining the results of (3.19) we obtain the following expression for the dimensionless energy-momentum tensor correlator in momentum space from which we identify the dimensionless transport coefficient κ/T 2 as At fixed temperature the continuum limit a µ → 0 is performed by taking N τ → ∞, where we reproduce the result of equation (2.16).
Although the computation has been performed in the ideal gas limit and thus lacks corrections in the coupling g, it may serve as a check of our numerics at high temperatures and helps to estimate the size of cut-off effects. The computed correction in the inverse temporal lattice extent suggests an anisotropy of ξ ≈ 2.33 in order to eliminate leading order lattice artefacts. In the case of other values for the anisotropy we can determine the required temporal lattice extent to decrease the leading discretisation error below a desired treshold. As stated in section 4 we use ξ = 2 in order to use previous results for the scale setting. Thus a temporal lattice extent of N τ ≥ 6 is required in order to reduce the leading lattice artefacts below 10% in the ideal gas limit. Note that an anisotropy larger than ξ > 2.33 causes a quadratic increase of the lattice artefacts, though it would milden the constraint (3.11).

Renormalisation
The correlator defined in (2.8) suffers from ultraviolet divergences. Although they become finite on the lattice, we have to correct the correlator by additive renormalisation. Therefore we subtract the vacuum part, which is defined as the correlator computed at vanishing temperature, from the measured correlator. We define a new vacuum corrected expectation value by where O T is an observable evaluated at a given temperature T and O Tvac its vacuum contribution, i.e. evaluated at vanishing temperature T vac = 0. The energy-momentum tensor is the Noether current corresponding to translational invariance. In the continuum it is protected from renormalisation by Wardidentities [39]. However, on the lattice translations only form a discrete symmetry group and thus multiplicative renormalisation becomes necessary. (The lattice perturbation theory computation in section 3.3 does not require multiplicative renormalisation because it is the non-interacting case).
For an isotropic lattice the finite renormalisation factor only depends on the lattice coupling β whereas on an anisotropic lattice it also depends on the anisotropy ξ. Additionally, temporal and spatial direction (3.6) require separate renormalisation factors Z σ (β, ξ) and Z τ (β, ξ). Then the renormalised energy-momentum tensor in the diagonal channel reads Applying the cubic symmetry (3.4) we rewrite the correlator (3.9) using the above notation and find with the newly defined bare correlators

25c)
and their vacuum subtracted versions Performing the renormalisation procedure we need the ratio Z σ (β, ξ 0 )/Z τ (β, ξ 0 ) and the absolute scale Z τ (β, ξ 0 ). The former can be obtained from renormalisation group invariant quantities [40]. To this end one introduces three differently sized lattices 27) and the renormalisation group invariant quantities Since the renormalisation factors do not depend on the temperature, all directions are symmetric and it follows where the expectation values are computed by lattice simulations of (3.28). We compute the ratio Z σ (β, ξ 0 )/Z τ (β, ξ 0 ) from all three equations in (3.29) and average the results. The simulations have to be performed for every lattice coupling β and anisotropy ξ. We obtain the absolute renormalisation factor by utilising the physical interpretation of the energy-momentum tensor, whose diagonal spatial elements are equivalent to the pressure θ ii = p. (3.31) The absolute renormalisation factor enters the energy-momentum tensor correlator quadratically. Therefore the renormalisation procedure is very sensitive to the exact value of the pressure and encourages us to use a highly precise value for it. For this reason we use the continuum extrapolated lattice data from [41]. Figure 3 illustrates the difference between the continuum value of the pressure and the not multiplicatively renormalised energy-momentum tensor. The difference between them at a given temperature corresponds to the absolute renormalisation factor. ii /T 4 for N τ = 6 and ξ = 2 to the continuum extrapolated pressure p/T 4 from the lattice [41], where the line is obtained by a cubic spline interpolation. The difference between them at a given temperature corresponds to the absolute renormalisation factor.

Numerical setup
We create the gauge field configurations using the standard heatbath algorithm [42][43][44] adapted to an anisotropic lattice. Our implementation is based on the library QDP++ [45].
In order to compute the vacuum part necessary for additive renormalisation, we run extra simulations with increased temporal lattice extent N τ . For our fine and spatially large lattices this is very costly. We therefore choose T vac ≈ 0.8 T c , with the critical temperature T c ≈ 260 MeV for Yang-Mills theory [46]. For our purposes this temperature is low enough since firstly the vacuum divergence is temperature independent, and secondly it is well known that the pressure or the deviation of screening masses from their vacuum values are exponentially small in the confined phase (see [41,[46][47][48] for numerical evidence and [49] for an analytic explanation).
The set of momenta has to fulfil the constraint (3.11), which basically dictates the simulation parameters. An anisotropy ξ > 1 benefits this constraint. As discussed in section 3.3 a value for the anisotropy of ξ ≈ 2.33 minimizes the first order lattice corrections. However, we choose an anisotropy of ξ = 2, which allows to set the scale by taking the lattice spacing as a function of the lattice coupling a = a(β) from [32]. Adjusting the temporal lattice extent to N τ ≥ 6 reduces the computed lattice errors in (3.21) below 10%. A numerical analysis of the relevant correlators in lattice QCD [36] even suggests values for the temporal lattice extent of N τ ≥ 8.
Extracting the transport coefficient κ from (2.14) by performing a linear fit in q 2 3 requires at least three different momenta q 3 , where the highest momentum still has to fulfil the constraint (3.11). More momenta would be favourable improving the fit's quality. Thus we choose for the temporal lattice extent N τ = 6 and for the spatial lattice extent N σ = 120 at a given anisotropy ξ = 2. All simulation parameters are listed in table 1. In the deconfined phase topological fluctuations are suppressed [50] and we expect no difficulties in using very fine lattices. Due to the large computational effort creating gauge fields on 120 3 × N τ lattices, we do not exclude any configurations but account for existing correlations by jackknife error sampling, see e.g. [51]. The multiplicative renormalisation procedure requires knowledge of the renormalisation factor ratio Z σ (β, ξ 0 )/Z τ (β, ξ 0 ). As described in section 3.4 we determine it from computing the quantities (3.28) on lattices (3.27) with L = 48. The simulations must be performed for every lattice coupling β of table 1. Intermediate results for the computation of the renormalisation factors are shown in table 6 and table 7 in appendix D with reference to run (i) of table 1.

Comparison to lattice perturbation theory
Our first simulation aims at making contact to lattice peturbation theory, section 3.3. The weak coupling regime is reached by increasing the temperature. Adopting the parameters from the previous section 4.1 we choose for the lattice coupling β = 7.1, corresponding to a temperature of T = 9.4 T c , and a spatial lattice spacing of a σ = 0.026 fm (see column (i) in table 1).   Figure 4 shows the correlator G E (q)/T 4 for five momenta compared to the result from lattice perturbation theory and table 2 the corresponding data points. The large errors of the correlator are almost entirely due to the additive renormalisation procedure. Table 4 lists the data of the bare correlators (3.25) regarding this simulation, whereas table 5 lists the data of the additively renormalised correlators (3.26). The vacuum subtraction causes a significant loss of accuracy. Computing the pressure by means of the interaction measure [46] suffers from the same phenomenon. Thus, we create a large amount of statistics (see table 1) to provide a significant signal for the correlators. In terms of error reduction it is highly favourable to perform the additive renormalisation before the multiplicative one. Otherwise, the propagated errors entering from the multiplicative renormalisation add to the described loss of precision.
Fitting the datapoints of the correlator to a line yields for the y-intercept G (0)/T 4 = 0.69(4) and for the transport coefficient κ/T 2 = 0.40 (26), which is consistent with the leading order lattice perturbation theory result κ LPT /T 2 = 0.47. Note that full agreement is not yet expected since at T = 9.4 T c there are still significant corrections due to interactions, i.e. we are still far from the ideal gas limit.

Temperature dependence
In principle, the temperature can be varied at fixed β and lattice spacing by changing N τ , where lower temperature implies larger N τ . However, due to the constraint on the momenta from equation (3.11) this would require a similar increase of the spatial volume and thus a drastical growth of the numerical effort. Hence the fixed scale approach is not practical for temperatures approaching the phase transition. We therefore investigate the temperature dependence of κ at fixed N τ /N σ by repeating the simulations at various lattice couplings β. In this case the different temperatures are evaluated at different lattice spacings, and consequently also different spatial volumes in physical units. However, since our lattice spacings are all a σ < 0.1 fm, we expect the lattice artefacts on the temperature dependence of the transport coefficient κ/T 2 to be negligible. As a consistency check for this, we also perform simulations at different temperatures but the same lattice spacings (simulations (i) and (ii) in table 1).
The results are shown in figure 5. The datapoint at T = 7.1 T c suffers from large errorbars since the spatial lattice extents have been kept fixed while increasing the temporal lattice extent N τ . This corresponds to less momenta fulfilling the constraint (3.11) and generates a loss of accuracy in the fit (4.1). Within the errorbars, the values of κ/T 2 at T = 9.4 T c and T = 7.1 T c agree (c.f. table 3) and thus justify the comparison at different lattice spacings and temperatures.
The numerical values for the transport coefficient κ are also summarised in table 3. Within errorbars, the temperature dependence of the transport coefficient is consistent with that of the ideal gas, κ ∼ T 2 , which is also the prediction of AdS/CFT [24] for the opposite strong coupling limit. Assuming this functional dependence, we may increase the accuracy by averaging the data points with N τ = 6 to give our final result, κ avr = 0.36(15)T 2 . (4. 2) The prediction from AdS/CFT correspondence for this coefficient is [27] κ where η is the shear viscosity and s the entropy density. The latter is proportional to the number of degrees of freedom of the theory, which is higher in the SUSY Yang-Mills used for the correspondence 1 . In order to compare with the QCD calculations, we thus use the AdS/CFT results η/s = 1/4π and eq. (4.3), but take the QCD entropy density from a lattice calculation [41]. The result is about a third of the perturbative prediction and also consistent with the simulation results.

Conclusions
We have calculated the second order hydrodynamic transport coefficient κ for the Yang-Mills plasma using lattice perturbation theory and Monte Carlo simulations.  Table 3: Lattice results for the transport coefficient κ/T 2 at different spatial lattice spacings a σ and temperatures T /T c . This is possible because the retarded correlator of the energy momentum tensor at zero frequency has a trivial analytic continuation to a corresponding Euclidean correlator. The transport coefficient parametrises the low momentum behaviour of this correlator, whose realisation requires large spatial lattice directions, making a numerical calculation very challenging and thus leaving large statistical errors. Their main source are the vacuum subtractions leading to similar problems in calculations of the equation of state at low temperatures. One might hope that alternative methods avoiding this step [52] may improve this situation.
In the investigated temperature range 2T c < T < 10T c our data are consistent with κ ∝ T 2 , as predicted both by weak and strong coupling methods. Because of still large errorbars, our result also quantitatively covers both the leading order perturbative as well as the AdS/CFT prediction rescaled by the QCD entropy. This would suggest that, besides improved simulation methods, next-to-leading order analytic calculations should be able to give a result with improved accuracy.

Acknowledgments
We thank H. Meyer and D. Rischke for their comments on the manuscript. This project is supported by the Helmholtz International Center for FAIR within the LOEWE program of the State of Hesse. C.S. acknowledges travel support by the Helmholtz Research School H-QM. The calculations have been performed on LOEWE-CSC at Goethe-University Frankfurt, we thank the team of administrators for support.

A Cubic symmetry of the energy momentum tensor
The correlator T 12 (x)T 12 (y) can be expressed in terms of diagonal energy-momentum tensor elements by exploiting rotation invariance on a spatially isotropic lattice (and medium) under rotations by α = π/4 about the z-direction. The transformation of a second rank tensor reads and the corresponding transformation matrix is given by For the energy-mometum tensor components of interest, this means where we used T 12 = T 21 . With the definition of the energy-momentum tensor (2.15a) we find for the correlator Note that the trace anomaly θ cancels completely. From rotational invariance follows

B Definitions in lattice perturbation theory
The Fourier transforms of the gauge field A µ to momentum space and back are defined by The shift to the center of the link variables x + a µ µ/2 in the Fourier transform simplifies the computation. The free gauge field propagator is given by where we use Feynman-'t Hooft gauge with ξ = 1. The momenta in lattice perturbation theory are given by with a 0 ≡ a τ and a i ≡ a σ . We do not imply a sum over the index µ.