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

In this paper we study the SU(3)-gluodynamics shear viscosity temperature dependence on the lattice. To do so, we measure the correlation functions of the energy-momentum tensor in the range of temperatures T /Tc ∈ [0.9, 1.5]. To extract the shear viscosity we used two approaches. The first one is to fit the lattice data with a physically motivated ansatz for the spectral function with unknown parameters and then determine the shear viscosity. The second approach is to apply the Backus-Gilbert method allowing to extract the 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 range T /Tc ∈ [0.9, 1.5] the SU(3)-gluodynamics reveals the properties of a strongly interacting system, which cannot be described perturbatively, and has the ratio η/s close to the value 1/4π of the N = 4 Supersymmetric Yang-Mills theory.


Introduction
Modern heavy ion collisions experiments such as RHIC and LHC are aimed at studying the quark-gluon plasma (QGP). The hydrodynamic description of QGP evolution proved to be efficient in understanding the experimental results [1,2]. Despite this success, the hydrodynamics is an effective theory which correctly represents the dynamics of infrared degrees of freedom only. The parameters of this effective theory, such as shear viscosity, bulk viscosity, conductivity, etc. cannot be calculated within the hydrodynamics itself but must be determined either from experiment or from a calculation based on first principles.
The measurement of the elliptic flow [3,4] allows to estimate the value of the QGP shear viscosity. In particular, the hydrodynamic approximation describes the experimental data if the ratio of shear viscosity η to entropy density s lies within the range η/s = (1 − 2.5) × 1/4π [5]. This value is close to the result of the N = 4 Supersymmetric Yang-Mills (SYM) theory at the strong coupling η/s = 1/4π [6] and deviates from the calculation at the weak coupling η/s ∼ const/(g 4 log(1/g)) ∼ 1 [7,8].
From this consideration one can conclude that the small value of the ratio η/s is governed by the nonperturbative dynamics. Attempts of the nonperturbative calculation of shear viscosity were undertaken in [9][10][11][12]. Unfortunately it is rather difficult to estimate the systematic uncertainty of these approaches, so today the analytical approach that fully accounts for the nonperturbative dynamics of QGP based on the first principles is absent. For this reason the only way to calculate the shear viscosity of QGP is the lattice simulation of QCD.
Despite a considerable progress in the lattice study of the QCD properties today it is not possible to calculate the shear viscosity of QGP with dynamical quarks. Even shear viscosity study within the pure gluodynamics is an extremely complicated task. There are only a few attempts to calculate the shear viscosity of the SU(3)-gluodynamics undertaken in [13][14][15][16][17] and the SU(2)-gluodynamics undertaken in [18,19]. In this paper we are going

JHEP04(2017)101
to study the temperature dependence of the shear viscosity of the SU(3)-gluodynamics in the vicinity of the confinement/deconfinement phase transition T /T c ∈ [0.9, 1.5].
The paper is organized as follows. In the next section we describe the details of the calculation. Section 3 is devoted to the calculation of the shear viscosity from the lattice measurements of the energy-momentum tensor correlation function. In the last section our results are discussed and the conclusion is drawn.

Details of the calculation
Shear viscosity is related to the Euclidean correlation function of the energy-momentum tensor T µν = 1 4 δ µν F a αβ F a αβ − F a µα F a να (here for simplicity we omitted the trace anomaly): where T is the temperature of the system. The correlation function (2.1) can be written in terms of the spectral function ρ(ω) as follows (2. 2) The spectral function contains valuable information about the properties of a medium. To find the shear viscosity from the spectral function one uses the Kubo formula [20] The lattice calculation of the shear viscosity can be divided into two parts. The first part is the measurement of the correlation function C(x 0 ) with a sufficient accuracy. This part of the calculation requires large computational resources but for gluodynamics the accuracy of the correlator can be dramatically improved with the help of the two-level algorithm [21]. 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 a continuous spectral function ρ(ω) from the integral equation (2.2) for the set of O(10) values of the function C(x 0 ) measured within the lattice simulation. Below we will need the properties of the spectral function. First we recall the very general properties: positivity of the spectral function ρ(ω)/ω 0 and oddness: ρ(−ω) = −ρ(ω). We also assume that the spectral function is a smooth function over frequency-scales of order of temperature. At large frequencies one expects that the asymptotic freedom manifests itself in the approach of the real spectral function to the one calculated at the weak coupling. For this reason it is also important to write the expression for the spectral function in the leading-order approximation in the strong coupling constant [22] ρ LO (ω) = 1 10 where d A = N 2 c − 1 = 8 for the SU(3)-gluodynamics.

JHEP04(2017)101
The next-to-leading order expression for the spectral function at large ω is also known [23]: It should be noted here that at the large ω the spectral function scales as ρ(ω) ∼ ω 4 , what leads to the large perturbative contribution to the correlation function for all values of the Euclidean time x 0 . Calculation shows that even at x 0 = 1/(2T ) the tree level contribution is ∼ 80 − 90% of the total value of the correlation function. Note also that the large ω behavior of the spectral function leads to a fast decrease of the correlation function C(x 0 ) ∼ 1/x 5 0 for the small x 0 . For this reason the signal/noise ratio for the C(x 0 ) is small at x 0 a and the lattice measurement of the correlation function at x 0 ∼ 1/(2T ) becomes computationally very expensive.
In the numerical simulation we use the Wilson gauge action for the SU(3)-gluodynamics where U µ,ν (x) is the product of link variables along the elementary rectangular (µ, ν), which starts at x. For the tensor F µν we use the clover discretization scheme: (2.7) One can easily define the energy-momentum tensor on the lattice using its continuum expression and the lattice discretization (2.7) for the F µν tensor. Note also that instead of the correlation function T 12 (x)T 12 (y) in this paper 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 [13]. Meanwhile the renormalization properties of the diagonal components of the energy-momentum tensor T µν are known (see below).
It has become conventional to present the value of shear viscosity as the ratio viscosityto-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 [24].
The energy-momentum tensor in the continuum theory is a set of the Noether currents which are related to the translation invariance of the action. In the lattice formulation of field theory continuum rotational invariance does not exist and the renormalization for energy-momentum tensor is required. For the correlation function considered in this paper the renormalization is multiplicative [25]. The renormalization factors depend on the discretization scheme. For instance, for the diagonal component of the T µν (when µ = ν) and the plaquette-based discretization of the T µν : lated to the anisotropy coefficients [24,26]: where c σ and c τ are defined in [24].
Using the renormalization factors for the plaquette-based discretization of the 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 (plaq) T We measured the correlation functions C(x 0 ) on the lattice 16 × 32 3 with the following set of the β-parameter: β = 6.491, 6.512, 6.532, 6.575, 6.647, 6.712, 6.811, 6.897, which corresponds to the temperatures T /T c 0.9, 0.925, 0.95, 1.0, 1.1, 1.2, 1.35, 1.5. Application of the two-level algorithm allowed us to get the uncertainties not larger than ∼ 2-3% at the distance T x 0 = 0.5 for all temperatures under consideration. For the other points the accuracy is even better. In figure 1 we plot the correlation functions for various temperatures. The next step in the calculation of shear viscosity is the extraction of the spectral function from the integral equation (2.2). In this section we are going to use the physically motivated ansätze for the spectral function with the unknown parameters. These parameters will be determined through the fitting procedure with full covariance matrix of our data.
The first ansatz which will be used in the calculation is motivated by the QCD sum rules [27]. In order to build a tentative spectral function we join the first-order hydrody-

JHEP04(2017)101
namic behavior at the small frequencies with the asymptotic freedom at large frequencies 1 In the last formula ρ lat (ω) is the 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 the clover discretization of the tensor F µν at the lattice with the fixed L t and L s → ∞. The function ρ lat (ω) was calculated in [19].
It was noted above that the correlation function C(x 0 ) is very sensitive to the ultraviolet properties of the spectral function. To get an accurate description of the lattice data we used the lattice spectral function at large frequencies ρ lat (ω) instead of the continuum expression (2.4). Application of the function ρ lat (ω) takes the discretization effects in the temporal direction into account and, as the result, considerably improves the quality of the fit.
Evidently the inclusion of the interactions modifies the tree level formula. However, due to the asymptotic freedom at large frequencies one can expect that this modification is not dramatic. In particular, it is seen from (2.5) that the one-loop radiative corrections modify the spectral function by a factor close to unity. Inspired by this observation, in expression (3.1) we multiplied the ρ lat (ω) by the factor A which effectively accounts the radiative corrections at large frequencies. Our results show that the introduction of the factor A is crucial for the successful description of the lattice data.
A drawback of ansatz (3.1) is that this function has a discontinuity at the point ω = ω 0 . It causes no difficulties to build the spectral function with the properties similar to ansatz (3.1) but free from this drawback Notice that in this ansatz we have introduced the parameter γ controlling the width of the transition from the infrared hydrodynamic regime to the ultraviolet regime of the asymptotic freedom in the spectral function. Evidently the width of the transition is ∼ 1/γ. In ansätze (3.1), (3.2) the first-order hydrodynamic regime directly continues to the regime of the asymptotic freedom. However, it is reasonable to assume that there is a region of frequencies where the spectral function deviates from the first-order hydrodynamics. In order to study this deviation in addition to spectral functions (3.1), (3.2) we use the following ansatz In the last formula we introduced a correction to the first-order hydrodynamic approximation which is controlled by the parameter C. We did not introduce the next-to-leading order correction to the first-order hydrodynamics ∼ ω 2 since the spectral function is an odd function of frequency. Notice that in the hard-thermal-loop framework the hydrodynamic behavior at small frequencies is replaced by the transport peak of a final width ω → ω/(1 + Γω 2 ) [28]. So, ansatz (3.3) can be considered as a first term of the expansion of the transport peak with Γ = −C.
We fit the lattice data (14 x 0 /a 2) for each temperature by formula (2.2) with spectral functions (3.1), (3.2), (3.3). In table 1 we show the parameters of functions (3.1), (3.2), (3.3) obtained through this fit. Instead of the parameter B in the fifth column of table 1 we show the ratio η/s which is related to the B as η/s = πB/s.

Now few comments are in order
• It is seen from table 1 that functions (3.1), (3.2), (3.3) fit the lattice data for various temperatures quite well. It is also seen that for all ansätze the ratio η/s quickly drops when temperature approaches the critical point T c and then either slowly rises after the T c or stays constant. This behavior was seen in the various models aimed at the calculation of shear viscosity in QCD.
• The values of η/s, A and ω 0 obtained through the fitting of the data by the various ansätze at the same temperature are in agreement with each other within the uncertainty of the calculation. However, the uncertainties for these parameters are different for the various ansätze.

JHEP04(2017)101
• The values of the threshold parameter ω 0 for all ansätze are physically well motivated. The value of the strong coupling constant at the threshold parameter ω 0 (ω 0 ∼ 2-3 GeV in the physical units) is α s (ω 0 ) ∼ 0.2-0.3. This allows us to expect that the perturbative expression for the spectral function is applicable for ω > ω 0 . The values of the factor A, which takes into account radiative corrections, are close to unity for all temperatures, this confirms the applicability of the asymptotic freedom at high frequencies.
• Notice also that contrary to the infrared part of the spectral function the parameters of the ultraviolet part for all ansätze are determined from the fit with a very good accuracy. This feature results from the fact that the dominant contribution to the correlation function is due to high frequencies.
• The ansatz ρ 3 (ω) allows to study the deviation from the first-order hydrodynamics. This deviation is controlled by the parameter C. From table 1 one sees that within the uncertainty of the calculation the values of the C are zero for all temperatures. This fact implies that our data does not allow to observe the deviation from the first-order hydrodynamics.
It is worth noting that we tried to fit our data by the spectral function similar to the ρ 1 (ω) but with the substitution ω → ω/(1 + Γω 2 ). This substitution accounts for the transport peak [28]. The result of this fit is very similar to that for ansatz (3.3). The parameters Γ equal zero within the uncertainty of the calculation.
Low frequency parts of spectral functions (3.1) and (3.2) are given by the first-order hydrodynamic expression ∼ ω. One can expect that the first-order hydrodynamic approximation works well up to ω πT 1 GeV [29]. From the other side the high frequency perturbative expression for the spectral function is fixed very accurately and it works well for ω ω 0 ∼ 3 GeV. The form of the spectral function in the region 1 GeV ω 3 GeV is not clear. We believe that the poor knowledge of the spectral function in the region 1 GeV ω 3 GeV is the main uncertainty source of the calculation based on the fitting of the lattice data by the functions ρ 1 (ω) and ρ 2 (ω). Notice that this source of uncertainty in the values of shear viscosity is not accounted in table 1.
The function ρ 3 (ω) modifies the first-order hydrodynamic expression in the intermediate region due to the term ∼ Cω 2 . Thus the function ρ 3 (ω) at least partly takes uncertainty in shear viscosity due to our poor knowledge of the spectral function in the intermediate region into account. For this reason we take the results for the ratios η/s obtained through the fitting by the function ρ 3 (ω) as the results of this section. In addition to the statistical uncertainties in the ratio η/s shown in table 1, there are uncertainties in the entropy density s and the renormalization coefficient of the clover discretized energy-momentum tensor (2.7). The former uncertainties are 4-6 % for all temperatures under consideration. The latter uncertainties are ∼ 3 % for the temperatures T /T c 1.0, ∼ 6 % for temperatures T /T c = 0.925, 0.95 and ∼ 12 % for the temperature T /T c = 0.9. The results for the ratios η/s obtained within the fitting procedure including all uncertainties are shown in the second column of table 2 and in figure 5.

The Backus-Gilbert method for the spectral function
In this section we are going to determine the ratio η/s using the Backus-Gilbert(BG) method [30,31]. 2 This approach has a considerable advantage over the method based on the fitting procedure: one does not need to know the parametrical form of the spectral function to carry out the calculation.
The method can be formulated as follows. 3 One needs to solve the equation (2.2). To do this we rewrite it in the following form where f (ω) is some function with the property f (t)| t→0 ∼ t, x i are the lattice points where the calculation of the C(x i ) is carried out and the K(x, ω) is the rescaled kernel of the integral equation (3.5) Our aim is to determine the ρ(ω). In the BG method instead of the 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 Evidently to get a better approximation for the ρ(ω) by the estimatorρ(ω) one should minimize the width of the δ(ω, ω). However, a very narrow peak might build an estimator fitting the points themselves, but not the physics (generality) they present. This means that any method of this kind should be regularized. The Backus-Gilbert is aimed at minimization of the Backus-Gilbert functional H(ρ(ω)) = λA(ρ(ω)) + (1 − λ)B(ρ(ω)). The component A represents the width of the resolution function (the second moment of distribution): A = ∞ 0 dωδ(ω, ω)(ω −ω) 2 . In principle, it could be any other function with the same physical meaning. The advantage of the second moment is that it is quadratic in ω andω, making analytical minimization possible.
Putting everything together, the minimization of the H functional gives the following values of the coefficients If the λ is close to 1, the resolution function has the smallest width and the estimator gives the best approximation for the spectral function. However, the application of the Backus-Gilbert method with the λ ∼ 1 to the calculation of shear viscosity gives rise to large uncertainties. The result becomes very dependent on the data, the spectral function turns out to be noisy and unstable. Statistical uncertainties can be reduced at the expense of increasing the width of the resolution function through decreasing the value of the λ-parameter. Now let us discuss the choice of the function f (x). To determine the spectral density at small frequencies it is reasonable to choose the function (3.12) In this case the ratio ρ(ω)/f 1 (ω)| ω→0 = η/π. The main motivation for choosing this function is that it gives a rather small width of the resolution function at small frequencies.
To study the spectral function at large frequencies we choose the following function (3.13) One can expect that due to asymptotic freedom at large frequencies ω Λ QCD the ratio ρ(ω)/f 2 (ω) behaves like constant. The width of the resolution function with the f 2 (x) is larger than that with the f 1 (x). Now let us proceed to the shear viscosity calculation. To do this we use the function f 1 (x). As to the λ-parameter we choose the λ = 0.002. For this value of the λ the uncertainty in the restored spectral density at zero frequency is smaller than 1% for all temperatures under consideration. In figure 2 we plot the resolution function at the temperature T /T c = 1.35 andω = 0 for various values of the λ. The resolution functions for the other temperatures are very similar to those in figure 2 and we do not show them here.
From figure 2 one sees that the width of the resolution function at λ = 0.002 is ∆ω ∼ 4T . If for a while we forget about the ultraviolet contribution, the convolution of the spectral function ρ(ω) with resolution function (3.6) gives some average of the spectral function over the interval with the width ∼ 4T . One can expect that the first-order hydrodynamic approximation works well up to ω πT 1 GeV [29], what covers most of the interval (0, 4T ) Now let us discuss the ultraviolet contribution to convolution (3.6). According to the results of the previous section the ultraviolet part of the spectral function starts to work for frequencies ω/T ∼ 7-8. From figure 2 one sees that the resolution function is considerably suppressed in this region. However, it is not possible to disregard ultraviolet contribution since the spectral function at large frequencies rises very quickly ρ(ω) ∼ ω 4 . The calculation shows that for the most temperatures under consideration the contribution which results from the ultraviolet part of the spectral function is larger than the contribution of hydrodynamic part of the spectral function. So, to get a reliable result for shear viscosity one should subtract the ultraviolet contribution.
To study the spectral function at large frequencies we are going to use the BG method with the function f 2 (ω) and λ = 0.002. In figure 3 we plot the ratiosρ(ω)/f 2 (ω) as the functions ofωa for the temperatures T /T c = 0.9, 1.1, 1.35, 1.5. The figures are similar for the other temperatures under consideration and we do not show them here. The red curves correspond to the spectral functions restored by the BG method from our data. In order to compare the results of the BG method with the results obtained in the previous section we took the spectral function at large frequencies (3.14) and convoluted it with the resolution function. The values of the A and ω with the uncertainties were taken from the fitting procedure (see table 1). The results are presented as the blue curves. Finally we plotted the dashed lines corresponding to the values of the constants A with the uncertainties obtained within the fitting procedure. Now few comments are in order • It is seen from figure 3 that the red curves can be separated into two parts. The first part is the spectral function for small frequenciesωa 0.5 (ω/T 8). One can say that the spectral function in this region is in the infrared regime. Afterωa ∼ 0.5 (ωT ∼ 8) there is a transition to the second regime where the spectral function is close to the ultraviolet asymptotic which is given by the constant A.
• It is also clear from figure 3 that the behavior of the blue curves representing the ultraviolet contribution to the ratioρ(ω)/f 2 (ω) is similar to the red ones. In the ultraviolet regime the blue and the red curves are close. Transition from the ultraviolet to the infrared regime takes place within the same region inωa.
• In the infrared region the red curves are higher than the blue ones. The difference between them can be attributed to the contribution of the spectral functions at small frequencies. One sees that the smaller the temperature the smaller the difference. If we recall that shear viscosity is related to the spectral function at small frequencies one can state that the shear viscosity drops with temperature. Our results assume that the entropy density s drops with temperature more quickly than the shear viscosity. As a result the ratio η/s rises below T c .
• From figure 3 one can see that in the ultraviolet region the restored spectral function is not a constant but a slowly varying function of theωa. The deviation of this function from the asymptotic value A obtained within the fitting procedure is very small for all temperatures. For the most temperatures the deviation is only a few percent. The deviation of the restored spectral function from the asymptotic value A can be attributed to the radiative corrections to the tree level spectral function which is more complicated than the constant.
The study carried out in this section allows us to state that formula (3.14) describes the ultraviolet part of the spectral function quite well. For this reason below we are going to use (3.14) as the model for the ultraviolet part of the spectral function. The value of the constant A will be determined from the variation of the restored ratioρ(ω)/f 2 (ω) in the regionωa ∈ (1.5, 3). This interval is chosen since the contribution of the infrared part of the spectral function is this region is small and the ratiosρ(ω)/f 2 (ω) for all temperatures are in the ultraviolet regime. The values of the constants A determined in this way are in agreement with that obtained within the fitting procedure.
In addition to the constant A the ultraviolet part (3.14) depends on the threshold parameters ω 0 . Thus, if one subtracts the ultraviolet contribution in form (3.14), the ratio η/s obtained within the BG method will depend on the value of ω 0 . To study this dependence in figure 4 we plot the ratio η/s as a function of ω 0 for the temperatures T /T c = 0.9, 0.925, 0.95, 1.0, 1.1, 1.5. The curves for the temperatures T /T c = 1.2, 1.35 are very close to the curve at T /T c = 1.5. For this reason we do not show these temperatures on the figure. The parameter A is taken at the central value of the fit of the lattice data by the function ρ 3 (ω) (see table 1). From figure 4 one sees that the larger the temperature the larger the slope of the curves and the weaker the dependence of the η/s on ω 0 . The dependence of the η/s on ω 0 is weak for the temperatures T /T c 1.0 and it is stronger for the T /T c < 1.0. The strongest dependence of the η/s on ω 0 is for the temperature T /T c = 0.9. We believe that this property stems from the already mentioned fact: shear viscosity of gluodynamics drops with temperature and the extraction of viscosity from the observable containing large ultraviolet contribution becomes more and more complicated for lower temperatures.
Unfortunately it is not quite clear how the value of the threshold parameter ω 0 can be determined within the BG method. Note, however, that the positions of the transition from the infrared to ultraviolet regime (see figure 3) coincide for both restored spectral function and for function (3.14) with ω 0 obtained within the fitting procedure. Note also that the values of the parameter A obtained within the BG method and the fitting procedure agree quite well. For this reason one can expect that the fitting procedure gives a good approximation for the value of ω 0 and we will take it for the model of the ultraviolet contribution.
Subtracting the ultraviolet contribution from the ratio η/s in form (3.14) we get the results of this section. These results are shown in the third column of table 2 and in (see the previous section). From table 2 and figure 5 one sees that the results obtained within two approaches applied in this paper agree with each other. In addition to the results obtained in this paper in figure 5 we plot the lattice results obtained in [14,15,17]. It is seen that our results are in agreement with the previous lattice studies.
It is also interesting to draw the results of [12]. In this paper the authors calculated the shear viscosity in the Yang-Mills theory using the exact diagrammatic representation in terms of full propagators and vertices using gluon spectral functions as an external input. Our results are in agreement with the results of this paper.
In figure 5 we also plot the value of the ratio η/s for the N = 4 SYM theory at the strong coupling η/s = 1/4π [6] and the results of the perturbative calculation of the η/s. Perturbative results were obtained as follows. The scale Λ for the running coupling constant in gluodynamics was taken from [35]. The entropy density s was taken at oneloop accuracy [36]. We took the perturbative results for shear viscosity at the next-toleading-log approximation from [8]. In order to estimate the uncertainty of the perturbative results we varied the scale in the region from the first to the second Matsubara frequency µ ∈ (2πT, 4πT ). Comparing our results with the other approaches one can conclude that within the temperature range T /T c ∈ [0.9, 1.5] SU(3)-gluodynamics reveals the properties of a strongly interacting system, which cannot be described perturbatively, and has the ratio η/s close to the value 1/4π in N = 4 Supersymmetric Yang-Mills theory.
In our calculation we use lattices with an aspect ratio of LT = 2. The calculation on lattices with a larger aspect ratio requires much more computation power. In [19]  calculated shear viscosity at the temperature T /T c = 1.2 for the SU(2)-gluodynamics on the lattices with LT = 2 and LT 2.4. Within the uncertainty of the of the calculation the results do not depend on the aspect ratio. We believe that this conclusion is valid for the SU(3)-gluodynamics and for temperatures that are not too close to the critical one. For the temperature T /T c = 1 one might expect the finite volume effects due to a large correlation length. However, the uncertainties of our results are sufficiently large and the point T /T c = 1 does not deviate considerably from the other points. For this reason we believe that even at the temperature T /T c = 1 the finite volume effects will not change our results.

Discussion and conclusion
This paper is aimed at studying the temperature dependence of the SU(3)-gluodynamics shear viscosity within lattice simulation. In particular, we measured the correlation functions of the energy-momentum tensor T 12 (x 0 )T 12 (0) at the lattice 16 × 32 3 for the temperatures T /T c ∈ [0.9, 1.5]. In order to get small uncertainties in our results we used the two-level algorithm which allowed to reach the accuracy not larger than ∼ 2-3% at the distance T x 0 = 0.5 for all the temperatures under consideration. For the other points the accuracy is better.  Using the lattice data for the correlation functions we calculated the ratios η/s for the temperatures under consideration. To do this we used the physically motivated ansätze for the spectral function with unknown parameters. These parameters were determined through the fitting procedure. All ansätze used in this paper are different ways of the interpolation between the hydrodynamic behavior at small frequencies and the asymptotic freedom at large frequencies. These ansätze fit the lattice data quite well for all temperatures. Another approach used to calculate the ratio η/s is the Backus-Gilbert method.

JHEP04(2017)101
In table 2 and figure 5 we plot the results obtained in this paper. From table 2 and figure 5 one sees that the results obtained within two approaches applied in this paper agree with each other.
In addition in figure 5 we plot the lattice results obtained in [14,15,17]. It is seen that our results are in agreement with the previous lattice studies. In figure 5 we also plot the value of the ratio η/s for the N = 4 SYM theory at strong coupling η/s = 1/4π and the results of perturbative calculation of the η/s. Comparing our results with other approaches one can conclude that the ratio η/s for gluodynamics is very close to the N = 4 SYM and cannot be described perturbatively.
It is also interesting to mention the results of [12], where the authors calculated the shear viscosity in the Yang-Mills theory using the exact diagrammatic representation in terms of full propagators and vertices using gluon spectral functions as an external input. Our results are in agreement with the results of this paper.
Currently it is not possible to carry out lattice calculation of shear viscosity in QCD with dynamical fermions. However, one can estimate the ratio η/s using the following formula The ratio η/s YM is calculated in this paper, while the ratio η s QCD / η s YM for N f = 3 quarks was estimated in [12]. For the η/s YM we took the result obtained using the JHEP04(2017)101  Figure 6. The ratio η/s in QCD for various temperatures. The red square points represent the estimation of the η/s in QCD done in this paper. The green curve is the result of the NJL model [10], the violet curve is the result of dynamical quasiparticle approach [11], blue curve is the result of [12] and black line is the result of the N = 4 SYM. In addition we plot the grey region which shows the experimental bound on the ratio η/s found from the experiment data [5].
BG method. Our results for the ratio η/s QCD are shown in table 2 and in figure 6.
In addition in figure 6 we plot the estimation of the η/s obtained within various models: NJL [10], dynamical quasiparticle approach [11], the result of [12] and N = 4 SYM. Finally in figure 6 we plot the gray region which shows the experimental bound on the ratio η/s found from experiment data [5]. It is seen that the results obtained in this paper agree with experiment.