Temperature dependence of shear viscosity of $SU(3)$--gluodynamics within lattice simulation

In this paper we study the shear viscosity temperature dependence of $SU(3)$--gluodynamics within lattice simulation. To do so, we measure the correlation functions of energy-momentum tensor in the range of temperatures $T/T_c\in [0.9, 1.5]$. To extract the values of shear viscosity we used two approaches. The first one is to fit the lattice data with some physically motivated ansatz for the spectral function with unknown parameters and then determine shear viscosity. The second approach is to apply the Backus-Gilbert method which allows to extract shear viscosity from the lattice data nonparametrically. The results obtained within both approaches agree with each other. Our results allow us to conclude that within the temperature range $T/T_c \in [0.9, 1.5]$ SU(3)--gluodynamics reveals the properties of a strongly interacting system, which cannot be described perturbatively, and has the ratio $\eta/s$ close to the value ${1}/{4\pi}$ in $N = 4$ Supersymmetric Yang-Mills theory.


Introduction
One of the main purposes of heavy ion collision experiments at RHIC and LHC is the experimental study of Quark-gluon plasma (QGP). It was shown that experimental results can be described within hydrodynamical approach [1,2]. In hydrodynamics transport coefficients such as shear viscosity, bulk viscosity, conductivity, etc. determine the dynamics of conserved charges. These transport coefficients can be determined from the experiment or calculated theoretically. For instance, the elliptic flow measurement [3,4] allows to determine the shear viscosity of QGP. It turns out that the ratio of shear viscosity η to entropy density s lies in range η/s = (1 − 3) × 1/4π [7], which is close to N = 4 Super Yang Mills (SYM) theory prediction at strong coupling η/s = 1/4π [8] and differs from a weak coupling result η/s ∼ c/(g 4 log(1/g)) ∼ 1 [5,6]. Thus a ratio η/s is very small and cannot be described within perturbation theory. It would be interesting to study the value of the shear viscosity within such nonperturbative approach as lattice simulations of QCD.

Details of the calculation
Shear viscosity is related to the Euclidean correlation function of the energy-momentum tensor T µν (here we omitted for simplicity trace anomaly): where T is the temperature of the system. The correlation function (2) can be written in terms of the spectral function ρ(ω) as follows In order to find shear viscosity from the spectral function one uses the Kubo formula [16]: Lattice calculation of shear viscosity can be divided into two parts. The first part is the measurement of the correlation function C(x 0 ) with sufficient accuracy. This part of the calculation requires large computational resources but for the gluodynamics the accuracy of the correlator can be dramatically improved with the help of the two-level algorithm [17]. The second part is the determination of the spectral function ρ(ω) from the correlation function C(x 0 ). The last part of the calculation is probably the most complicated, since one should determine continuous spectral function ρ(ω) from integral equation (3) for the set of O(10) values of the function C(x 0 ) measured in lattice simulation.
First we recall very general properties: the positivity of the spectral function ρ(ω)/ω ≥ 0 and oddness: ρ(−ω) = −ρ(ω). At large frequencies one expects asymptotic freedom behaviour. For this reason it is also important to write the expression for the spectral function in the leading-order approximation in strong coupling constant expansion [18]: where d A = N 2 c − 1 = 8 for the SU(3)-gluodynamics. One also knows next-to-leading order expression for the spectral function at large ω [19]: It should be noted here that at large ω the spectral function scales as ρ(ω) ∼ ω 4 , what leads to a large perturbative contribution to the correlation function for all values of Euclidean time x 0 . Calculation shows that even at the x 0 = 1/(2T ) the tree level contribution is ∼ 80 − 90% of the total value of the correlation function. Note also that large ω behaviour of the spectral function leads to a fast decrease of the correlation function C(x 0 ) ∼ 1/x 5 0 for small x 0 . For this reason the signal/noise ratio for the C(x 0 ) is small at x 0 a and lattice measurement of the correlation function at x 0 ∼ 1/(2T ) becomes computationally very expensive.
In numerical simulation we use Wilson gauge action for the SU(3)-gluodynamics where U µ,ν (x) is the product of the link variables along elementary rectangular (µ, ν), which starts at x.
For the tensor F µν we use the clover discretization scheme: It causes no difficulties to build energy-momentum tensor from Eq. (1) having expression for the tensor F µν .
To calculate shear viscosity one should measure the correlation function (2) using two-level algorithm [17] to decrease computation time. Note also that instead of the correlation function T 12 (x)T 12 (y) we measure the correlation function 1 2 ( T 11 (x)T 11 (y) − T 11 (x)T 22 (y) ). Both correlation functions are equal in the continuum limit [9].
It is conventional to present the value of shear viscosity as a the ratio viscosity-to-entropy density η/s. For homogeneous systems the entropy density s can be expressed as s = +p T , where is the energy density and p is the pressure. These thermodynamic quantities were measured with the method described in [20].
Energy-momentum tensor in continuum theory is a set of Noether currents which are related to the translation invariance of the action. In lattice formulation of field theory continuum translation invariance does not exist and renormalization for energy-momentum tensor is required. For the correlation function considered in this paper the renormalization is multiplicative [21]. Renormalization factors depend on the discretization scheme. For instance, for the diagonal component of T µν (when µ = ν) and the plaquette-based discretization of T µν : the renormalization factors are related to the anisotropy coefficients [22,23]: where c σ and c τ are defined in [20] and for S U(2) are computed on the lattice in [24].
Using the renormalization factors for the plaquette-based discretization of T 00 , we can find the renormalization factors for the clover discretization simply by fitting the vacuum expectation values of the renormalized T 00 : Z  Application of two-level algorithm allowed us to get the uncertainty ∼ 2 − 3% at the distance T x 0 = 1/2 for all lattices. For the other points the accuracy is better. In Fig. 1 we plot the correlation functions for various temperatures.
The next part of the calculation is the extraction of the spectral function from the integral equation (3). In this section we are going to use physically motivated ansatz for the spectral function with unknown parameters that can be determined through the fitting procedure.
In the calculation of viscosity we use ansatz motivated by QCD sum rules [25]. We join hydrodynamic behaviour at small frequencies with asymptotic freedom at large frequencies 1 In the last formula ρ lat (ω) is a tree level lattice expression for the spectral function calculated for the correlation function ∼ 1 2 ( T 11 (x)T 11 (y) − T 11 (x)T 22 (y) ) with clover discretization of the tensor F µν at lattice with fixed L t and L s → ∞. The function ρ lat (ω) was calculated in paper [15].
The correlation function C(x 0 ) is very sensitive to the ultraviolet properties of the spectral function, so we use ρ lat (ω) instead of the continuum expression (5). Application of the function ρ lat (ω) considerably improves the quality of the fit.
Evidently the inclusion of interactions modifies tree level formula, so we multiplied the ρ lat (ω) by some factor A which effectively takes into account radiative corrections at large frequencies. Introduction of the factor A is very important for successful description of lattice data.
The obtained values ω 0 ∼ 3 GeV are physically well motivated. The value of the strong coupling constant is α s (ω 0 ) ∼ 0.2 − 0.3. This allows us to expect that perturbative expression for the spectral function is applicable for ω > ω 0 . The values of the factor A for all temperatures, which takes into account the effect of interactions, are close to unity what confirms applicability of the asymptotic freedom at high frequencies.
Low frequency part of the spectral function (9) is given by the first-order hydrodynamic expression ∼ ω. Comparison of the spectral functions of the energy-momentum tensor correlation functions obtained in N = 4 SYM [26] and the first-order hydrodynamic expressions allows us to expect that this approximation works well up to ω ≤ πT 1 GeV [27].
On the other hand high frequency perturbative expression for the spectral function is rigidly fixed and works well up to ω ≥ ω 0 ∼ 3 GeV. The form of the spectral function in the region 1 GeV ≤ ω ≤ 3 GeV is not clear. This is the main source of the uncertainty of the calculation based on the fitting procedure. In order to get rid of these uncertainties, one should apply non-parametrical estimation techniques.
The method can be formulated as follows 3 . One needs to solve equation (3). To do this we rewrite it in the following form where f (ω) is some function with the property f (t)| t→0 ∼ t and K(x, ω) is a rescaled kernel of the integral equation Our aim is to determine the ρ(ω). In the BG method instead of the true function ρ(ω) we consider the estimator of this functionρ(ω) which can be written as where the function δ(ω, ω) is called the resolution function. This function has a peak around the point ω and it is normalized as ∞ 0 dωδ(ω, ω) = 1. The sharper the peak, the closerρ is to ρ. The function is expanded over the K(x i , ω) as If a resolution functon is chosen this way, the estimator is a linear conbination of the values of the correlation functionρ The above formulae are typical for all non-parametrical estimation procedures. The specifics of the BG-method are the following: one minimizes the sum of functionals A + B, where A ensures that the estimator reproduces correlation function, the B functional punishes estimator for being too dependent on the data (prevents overfittng).
To get a better approximation for the ρ(ω) by the estimatorρ(ω) one should sharpen the δ(ω, ω) function. This minimizes the A functional, and the BG-method uses the second moment of the resolution function as the measure of its width: Minimization requirement gives the following values of the coefficients The matrix W i j (ω) is poorly conditioned and there are large cancelations between different terms in equation (14). To make extimator less dependent on the data, one should take the B functional int account. This regularizes the matrix W i j (ω) by the covariance matrix of our data S i j The λ parameter represents the trade-off between stability and data fitting. The exact λ value should be chosen depending on the matrix W condition number.
To determine the spectral density at small frequencies ω Λ QCD it is reasonable to choose the function f (ω) = ω. In this case the ratio ρ(ω)/ω| ω→0 = η π . To calculate shear viscosity, we choose λ = 0.002. For this value of the λ the uncertainty in the restored spectral density at zero frequency is less than 1% for all temperatures.
As already noted, we expect that hydrodynamic approximation works well up to ω ∼ πT [27]. At the same time from the width of the resolution function at λ = 0.002 is ∆ω ∼ 4T . We conclude that the convolution of the spectral density with the resolution function at λ = 0.002 averages over reasonable region and gives good appoximation for the shear viscosity.
According to the results of previous section asymptotic freedom starts to work for frequencies ω/T ≥ 7 − 8. The resolution function is considerably suppressed in this region. However, the spectral function at large frequencies rises very quickly ρ(ω) ∼ ω 4 , so this convolution contributes to the answer. So, to get reliable estimate one should subtract the asymptotic freedom.
The simplest way to do this is to subtract the convolution ρ uv (0) = δ(0,ω)ρ uv (ω)dω, where the ρ uv (ω) = Aθ(ω − ω 0 )ρ lat (ω). The determination of A and ω 0 parameters can be done by many procedures, the simplest way is to fit colleration function C(τ) for low τ, where it is completely governed by the UV-contribution. Despite the fact that this method of subtraction makes assumptions about the UV spectral function form, it still leaves the intermediate and hydrodynamic regions untouched and the usage of BG method remains sensible.