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

In this paper the temperature dependence of the $SU(3)$--gluodynamics bulk viscosity is studied within lattice simulations. To carry out this study we measure the correlation function of the trace of the energy-momentum tensor for a set of temperatures within the range $T/T_c \in (0.9, 1.5)$. To extract the bulk viscosity from the correlation function we apply the Backus-Gilbert method and the Tikhonov regularization method. We show that the ratio $\zeta/s$ is small in the region $T/T_c \geqslant 1.1-1.2$ and in the vicinity of the transition $T/T_c \leqslant 1.1-1.2$ it quickly rises. Our results are in agreement with previous lattice studies and in a reasonable agreement with other phenomenological approaches. Obtained values of the bulk viscosity are significantly larger than perturbative results, what confirms that QGP is a strongly correlated system.


Introduction
Hydrodynamics is believed to describe time evolution of quark-gluon plasma (QGP) created in heavy ion collision experiments (such experiments are carried out at RHIC and LHC and planned in the future at FAIR and NICA). The basic object in hydrodynamics is the energy-momentum tensor built as an expansion in gradients [1,2]. The leading order of this expansion describes an ideal fluid. The next-to-leading order includes dissipation and can be parameterized by two coefficients: shear and bulk viscosities. Trying to describe the particle yield in heavy ion collisions one can determine the shear and bulk viscosities of QGP [3]. In particular, typical value of the shear viscosity to the entropy density ratio extracted from the hydrodynamic studies is η/s = (1 − 2.5) × 1/4π [4].
It is therefore important to calculate these observables based on our theoretical knowledge of the system. Since QGP is a strongly correlated system, one of the main ways to carry out the first principle study of its properties is the lattice simulation of QCD. Despite considerable success in lattice study of QGP, lattice calculations of the shear and bulk viscosities still remain challenging problems requiring huge statistics. For this reason nowadays it is not feasible to calculate viscosities in QCD with dynamical quarks. In the following we are going to address viscosities in gluodynamics.
Lattice calculations of the gluodynamics shear viscosity were carried out in [5][6][7][8][9][10][11][12]. They are in agreement with each other and with the experimental data [4] within the uncertainties. The lattice results are also close to η/s = 1/4π obtained within the N = 4 supersymmetric Yang-Mills theory at strong coupling [13]. For many years there was a disagreement between lattice results and the perturbative calculations [14,15]. However, recent calculations in the next-to-leading order [16] give much smaller values of the shear viscosity, which are consistent with the lattice data. It is also worth to mention the paper [17], where the authors determined the value of the shear viscosity using the diagrammatic representation (their results are also very small and close to the lattice and the experimental data).
Another important transport coefficient is the QGP bulk viscosity ζ. There are only few rather old papers devoted to lattice calculation of the bulk viscosity [6,18]. Taking into account the rapid improvement of supercomputers and theoretical developments it is reasonable to conduct an up-to-date study of the bulk viscosity.
Let us now consider what is known about the bulk viscosity. One can expect that at very large temperature the results of perturbative calculation are applicable. It gives a very small value ζ/s ∼ 0.02α 2 s for N f = 0 [19]. Unfortunately, this result cannot be applied for temperatures ∼ few × the critical temperature T c , which is of interest for heavy ion collision experiments. The bulk viscosity in this region was studied in [20,21]. The authors applied the low energy theorems of QCD and derived the formula relating the spectral function of the energy-momentum tensor trace correlator to the energy density and pressure of hot matter. Using a physically motivated ansatz of the spectral function the authors found ζ. At the critical temperature the bulk viscosity has a peak with the height ζ/s ∼ 1. For temperatures T > T c the bulk viscosity quickly drops becoming very small ζ/s < 0.1 already for T > 1.1T c . There are a lot of phenomenological studies of the bulk viscosity [22][23][24][25][26][27][28] confirming the existence of peak at the critical temperature.
In these paper we study the temperature dependence of the bulk viscosity in SU (3)-gluodynamics within lattice simulation. We calculate the correlation functions of the energy-momentum tensor trace for the set of temperatures. To extract the spectral function and the bulk viscosity from the correlator we use two model-independent estimation techniques: the Backus-Gilbert method [29,30] and the Tikhonov regularization approach [31].
This paper is organized as follows. In the next section we describe the details of the lattice measurements of the correlation functions under study. The results of this measurement are presented in Section III. In the last section we discuss our results and draw the conclusion.

Details of the calculation
Bulk viscosity is related to the Euclidean correlation function of the trace of the energy-momentum tensor: where θ = β(g) 2g F a µν F a µν , β(g) is the β-function of gluodynamics, T is the temperature. The correlator (1) can be expressed in terms of the spectral function ρ(ω) via the integral equation The spectral function contains valuable information about medium properties. In particular, one can find the bulk viscosity if the spectral function is known [32] Lattice calculation of bulk viscosity can be divided into two parts. First, one measures the correlation function C(τ ) with sufficient accuracy. Second, one determines the spectral function ρ(ω) from C(τ ). Although the first step is very complicated, sufficient accuracy can be achieved due to the multi-level algorithm [33]. The second step is the well known ill-posed problem which is difficult to solve. Important properties of the spectral function are positivity ρ(ω) 0, ω > 0 and oddness ρ(−ω) = −ρ(ω). It is also important that at the leading order approximation in the strong coupling constant the spectral function can be written as where d A = N 2 c − 1 = 8 for the SU (3)-gluodynamics. We expect that due to the asymptotic freedom the leading order expression (4) is a good approximation for the spectral function at large frequency. To account for the discretization errors in temporal direction instead of (4) we are going to use the tree level lattice expression ρ lat (ω) calculated within the approximation: L t is fixed and L s → ∞ [34]. The resulting expression for ρ lat (ω) is cumbersome, for this reason we do not show it here.
At small frequencies the spectral function is in the hydrodynamic regime. The first order hydrodynamic behavior for the spectral function reads In numerical simulation we use the Wilson gauge action for the SU (3)-gluodynamics. For F µν the clover discretization scheme is used. Similar to the shear viscosity the bulk viscosity is presented as the viscosity-to-entropy-density ratio ζ/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 [35].
The energy-momentum tensor in the continuum theory is a set of the Noether currents which are related to the Lorentz invariance of the action. In the lattice formulation of field theory continuum rotational invariance does not exist and the renormalization of energy-momentum tensor is required. For the trace of the energy-momentum tensor the renormalization is multiplicative. The renormalization factors depend on the discretization scheme [36]. For the plaquette-based discretization scheme the renormalization is defined by the β-function [36]. Using the renormalization factors for the plaquette-based discretization of θ, one can easily find the renormalization factors for the clover discretization by fitting the expectation value of the trace anomaly:

Correlation functions and their properties
We measured the correlation functions C(τ ) on the lattice 16 × 32 3 with the parameters listed in Table I. β 6.491 6.512 6.532 6.552 6.575 6.61 6.647 6.682 6.712 6.765 6.811 6.855 6.897 Table I: The set of parameters used in the calculation of C(τ ).
The two-level algorithm allowed us to reach relative errors of 2-3% at the middle point τ T = 1/2 for all the temperatures. For other values of Euclidean time τ the relative errors are even smaller. In Fig. 1 we show correlation functions (1) for the temperatures T /T c = 0.90, 1.10, 1.35, 1.5. In order to estimate the finite volume effects, we measured the correlation functions (1) on the larger lattice 16 × 48 3 for the temperatures T /T c = 0.9, 0.975, 1.0, 1.05, 1.5. We found that for the temperatures T /T c = 0.9, 0.975 the deviations of the correlators measured on two lattices are less than 2σ. For the temperatures T /T c = 1.05, 1.5 the deviations are less than σ. For the temperature T /T c = 1.0 the deviation of the correlators in the middle point vicinity is as large as 4σ. The deviations for other points are smaller. Thus, we expect finite volume effects to be small for all the temperatures except T = T c . Finite volume effects for the temperature T = T c might be important.

Calculation of the bulk viscosity using the Backus-Gilbert method
In this section we determine the ratio ζ/s using the Backus-Gilbert (BG) method [29,30]. This is a model independent approach estimating the spectral function 1 . Instead of ρ(ω) one reconstructs the estimatorρ(ω) expressed asρ where the f (x) is an arbitrary function and the δ(ω, ω) is called the resolution function. This function has a peak aroundω and normalized as ∞ 0 dωδ(ω, ω) = 1. The BG resolution function is taken in the form For this resolution function the estimator is a linear combination of the correlation function values For better approximation of ρ(ω) with the estimatorρ(ω) one needs to minimize the width of δ(ω, ω). 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.
If λ is close to 1, the resolution function has the smallest width. However, the BG method with λ ∼ 1 leads to large uncertainties. The result becomes very dependent on the data and the spectral function turns out to be noisy and unstable. Statistical uncertainties are reduced at the expense of increasing the width of the resolution function through decreasing of λ.
The minimization of H gives In [11] it was seen that the BG method should be modified to study the shear viscosity. The estimator of the spectral function (6) is the convolution of the real spectral function ρ(ω) and the resolution function δ(ω,ω), which has an ultraviolet tail. The ultraviolet behavior of ρ(ω) is ∼ ω 4 and it is convolved with the tail of the resolution function. For this reason the shear viscosity calculated within the BG method acquires large ultraviolet contribution which is non-physical. To get rid of this problem it was proposed to determine the ultraviolet tail of the spectral function and then subtract it from the estimator. Then the result is the convolution of only the infrared part of the real spectral function and the resolution function.
To subtract the ultraviolet tail we are applied the approach proposed in [11]. The spectral function at large frequency is determined using the rescaling function We used the running coupling constant α s at one-loop level with Λ = 237 MeV [38]. In the Backus-Gilbert method one reconstructs the ratio ρ(ω)/f uv (ω) which is divergent at ω = Λ. To get rid of this divergence we assume that α(ω) = α(1GeV) for ω 1 GeV. The result is not sensitive to the modification of the α s -running since we study the large-frequency behavior of ρ(ω), which is not affected by this modification.
Finally, one has to fix the value of λ. In the BG method larger λ leads to larger uncertainties of the calculation. On the other hand small λ increases the width of the resolution function and one needs to find a compromise between these two tendencies. We found that this compromise is satisfied at λ = 0.01 for our study of the ultraviolet properties of the spectral function. For the infrared study (see below) we used λ = 0.1.
Within the BG method with (12) we reconstruct the ratioρ(ω)/f uv (ω). In Fig. 2 we plot the reconstructed ratios for few temperatures. We see that at large frequencies the ratioρ(ω)/f uv (ω) reaches the plateau. The role of the α 2 s (ω) factor in (12) should be emphasized: if it were not for the running coupling, the ratioρ(ω)/f uv (ω) would not have a plateau. This means that the function ρ lat itself does not catch the essential behavior at large frequencies and the account of α 2 s (ω) is necessary. We conclude that at large frequencies the spectral function behaves as ρ(ω) = Af uv (x). Based on this finding we propose the following form of the ultraviolet spectral function where A is the value ofρ(ω)/f uv (ω) on the plateau. In the calculation we determine the value and the uncertainties of A from the plateau in the regionωa ∈ (1.5, 3). Another parameter of the ultraviolet tail is ω 0 frequency threshold from which the spectral function is given by the ultraviolet form (13). This parameter will discussed below.
Having determined the ultraviolet behavior of the spectral function, we proceed to the calculation of the bulk viscosity. In order to calculate the ζ, we found the estimator for ρ(ω)/ω atω = 0. We calculate it using the BG method with f (x) = x.
The resolution function for the temperature T /T c = 1.05 is shown in Fig. 3 (solid line). The resolution functions at the other temperatures are close to that for the T /T c = 1.05. It is seen that the width of the resolution function is ∼ 4T . Thus the spectral function ρ(ω)/ω is averaged over the region ∼ 4T .
With the resolution function we calculate the estimator for ρ(ω)/ω. Then we subtract the ultraviolet contribution given by the convolution of the spectral function (13)  not be determined within the BG analysis. To account its uncertainty, we vary ω 0 within the region ω 0 /T c ∈ (5, 10) ω 0 ∼ (1.4, 2.8) GeV . We believe that this region is sufficiently safe to estimate the uncertainty due to the unknown value of ω 0 . Our results for the ratios ζ/T 3 and ζ/s are shown in Fig. 4. The uncertainties shown in Fig. 4 are due to statistical errors and the uncertainties in the A and ω 0 .

Calculation of the bulk viscosity using the Tikhonov regularization
Finally we consider another approach to the bulk viscosity estimation, called the Tikhonov regularization (TR) 2 [31]. The TR method allows to make the resolution function narrower as compared to the BG method. In the calculation with the TR we follow the formulas (11). The difference of the TR as compared to the BG method is in the regularization of the W ij matrix. In the BG method one adds the covariance matrix λS ij to the matrix (1 − λ)W ij . In the TR method the matrix W ij is regularized as follows. Let us consider the SVD decomposition of the W −1 : We substitute the matrix D by the matrix D = diag (σ 1 + γ) −1 , (σ 2 + γ) −1 , . . . , (σ n + γ) −1 , where the γ is the regularization parameter of the TR method. The TR method thus smoothly cuts-off small singular values σ, making the results more stable.
The application of the TR method is the same as the application of the BG, but with the W ij regularized in the other way. We first choose the regularization parameter γ = 0.1 which is a compromise between the width of the resolution function and the uncertainty of the calculation. The resolution function in this case is shown in Fig. 3 and it has the width ∼ 3T which is smaller than in the BG method. We then calculate the estimator with f (x) = x and subtract the ultraviolet tail. The ultraviolet spectral function was taken in the form (13) with A and ω 0 from the BG method. In Fig. 4 we plot our results.

Discussion and conclusion
In this paper the temperature dependence of the bulk viscosity in gluodynamics is studied within lattice simulation. To carry out this study we measured the correlation functions of the energy-momentum tensor trace for a set of temperatures in the range T /T c ∈ (0.9, 1.5). To extract the bulk viscosity from the correlation function we applied the Backus-Gilbert method and the Tikhonov regularization method. The results obtained within both approaches are shown in Fig. 4. We also studied the finite volume effects and found that they are small for all the temperatures except the T = T c . Finite volume effects for T = T c might be important.
Let us now consider the results obtained in this paper. From Fig. 4 we see that the ratio ζ/s is small in the region T /T c 1.1 − 1.2 and in the vicinity of the transition T /T c 1.1 − 1.2 it quickly rises. This behavior is in agreement with a lot of phenomenological studies (see, for instance, [20,23,27]).
Below the critical temperature ζ/s continues to rise. It is in disagreement with some phenomenological studies of QCD [25,28]. This discrepancy might be explained as follows. In our study we convolute the spectral function with the resolution function which has the width ∼ (3 − 4)T c . Below the critical temperature there is the scalar glueball contribution to the spectral function which might alter our result. One might also find another explanation of this discrepancy. In is known that the confinement/deconfinement phase transition in the SU (3)-gluodynamics is of the first order, while this transition in QCD is a crossover. The rise of the ζ/s below the critical point might be assigned to the rapid decrease of the entropy density s below the transition in the gluodynamics. Better understanding of this discrepancy requires additional study of the bulk viscosity and the spectral function of the energy-momentum tensor both from the lattice side and in these phenomenological models.  Figure 5: The ratio ζ/s calculated in this paper and other studies: the results obtained in this paper within the Backus-Gilbert method and the Tikhonov regularization, the lattice results obtained in [6] and [18], perturbative results obtained in [19] and the results of [20].
In order to compare our results with the results of other approaches, in Fig. 5 we plot the ratios ζ/s for T /T c 1 calculated in our paper and in other studies. In particular, the blue circles and the red triangles represent the results obtained in this paper within the Backus-Gilbert method and the Tikhonov regularization correspondingly.
The black circles and the yellow squares represent lattice results obtained in [6] and [18] correspondingly. It is seen that our results are in agreement with the previous lattice studies of the bulk viscosity.
The blue band represents the perturbative results obtained in [19]. The uncertainty in this band is due to the variation of the scale in the region µ ∈ (2πT, 4πT ). It is seen that out results dramatically disagree with the perturbative results, what once again confirms that QGP is a strongly correlated system.
The results of [20] are represented with the violet diamonds. It is interesting to notice that the rise of ζ/s in [20] starts at T /T c ∼ 1.1 which agrees with our results. Finally let us consider the following question. The perturbative calculations of the ζ revealed the following relation between the bulk viscosity, the shear viscosity η and the speed of sound v s : ζ/η ∝ (1 − 3v 2 s ) 2 . On the other hand similar ratio in the AdS/CFT is predicted to be: ζ/η ∝ (1 − 3v 2 s ) [40,41]. In addition it was argued that there is an inequality ζ/η 2/3(1 − 3v 2 s ) which is valid for QGP [41]. To check these assumptions in Fig. 6 we plot our results obtained within the Tikhonov regularization 3 in the region T /T c ∈ (1.05, 1.425) as a function of the speed of sound in gluodynamics calculated in [42]. The viscosity η was taken from our paper [11]. If the temperature at which ζ is calculated is not present in [11] we take the average of the closest points. We fit our data with the linear and the quadratic fits. We also plot the line ζ/η = 2/3(1 − 3v 2 s ). Unfortunately the uncertainty of the calculation is rather large and one can not distinguish between these two hypotheses.