A calculation of the shear viscosity in SU(3) gluodynamics
MIT-CTP 3830
A calculation of the shear viscosity in SU(3) gluodynamics
Harvey B. Meyer∗
Center for Theoretical Physics
Massachusetts Institute of Technology
Cambridge, MA 02139, U.S.A.
(Dated: October 22, 2018)
We perform a lattice Monte-Carlo calculation of the two-point functions of the energy-momentum
tensor at finite temperature in the SU(3) gauge theory. Unprecedented precision is obtained thanks
to a multi-level algorithm. The lattice operators are renormalized non-perturbatively and the clas-
sical discretization errors affecting the correlators are corrected for. A robust upper bound for the
shear viscosity to entropy density ratio is derived, η/s < 1.0, and our best estimate is η/s = 0.134(33)
at T = 1.65Tc under the assumption of smoothness of the spectral function in the low-frequency
region.
PACS numbers: 12.38.Gc, 12.38.Mh, 25.75.-q
Introduction.— Models treating the system produced in
heavy ion collisions at RHIC as an ideal fluid have had
significant success in describing the observed flow phe-
nomena [1, 2]. Subsequently the leading corrections due
to a finite shear viscosity were computed [3], in parti-
cular the flattening of the elliptic flow coefficient v2(pT)
above 1GeV. It is therefore important to compute the
QCD shear and bulk viscosities from first principles to
establish this description more firmly. Small transport
coefficients are a signature of strong interactions, which
lead to efficient transmission of momentum in the system.
Strong interactions in turn require non-perturbative com-
putational techniques. Several attempts have been made
to compute these observables on the lattice in the SU(3)
gauge theory [4, 5]. The underlying basis of these calcu-
lations are the Kubo formulas, which relate each trans-
port coefficient to a spectral function ρ(ω) at vanishing
frequency. Even on current computers, these calcula-
tions are highly non-trivial, due to the fall-off of the re-
levant correlators in Euclidean time (as x−50 at short dis-
tances), implying a poor signal-to-noise ratio in a stan-
dard Monte-Carlo calculation. The second difficulty is
to solve the ill-posed inverse problem for ρ(ω) given the
Euclidean correlator at a finite set of points. Mathemati-
cally speaking, the uncertainty on a transport coefficient
χ is infinite for any finite statistical accuracy, because
adding ǫωδ(ω) to ρ(ω) merely corresponds to adding a
constant to the Euclidean correlator of order ǫ, while ren-
dering χ infinite. Therefore smoothness assumptions on
ρ(ω) have to be made, which are reasonable far from the
one-particle energy eigenstates, and can be proved in the
hard-thermal-loop framework [6].
In this Letter we present a new calculation which dra-
matically improves on the statistical accuracy of the Eu-
clidean correlator relevant to the shear viscosity through
the use of a two-level algorithm [16]. This allows us to
derive a robust upper bound on the viscosity and a use-
∗Electronic address: meyerh@mit.edu
ful estimate of the ratio η/s, which has acquired a special
significance since its value 1/4π in a class of strongly cou-
pled supersymmetric gauge theories [7] was conjectured
to be an absolute lower bound for all substances [8].
Methodology.— In the continuum, the energy-
momentum tensor Tµν(x) = F
να − 14δµνF
ρσ , be-
ing a set of Noether currents associated with translations
in space and time, does not renormalize. With L0 = 1/T
the inverse temperature, we consider the Euclidean two-
point function (0 < x0 < L0)
C(x0) = L
d3x 〈T 12(0)T 12(x0,x)〉. (1)
The tree-level expression is Ct.l.(x0) =
f(τ)− π
with τ = 1 − 2x0
, dA = 8 the number of gluons and
f(z) =
ds s4 cosh2(zs)/ sinh2 s. The correlator C(x0)
is thus dimensionless and, in a conformal field theory,
would be a function of Tx0 only.
The spectral function is defined by
C(x0) = L
coshω(1
L0 − x0)
sinh ωL0
dω. (2)
The shear viscosity is given by [4, 13]
η(T ) = π
. (3)
Important properties of ρ are its positivity, ρ(ω)/ω ≥ 0
and parity, ρ(−ω) = −ρ(ω). The spectral function that
reproduces Ct.l.(x0) is
ρt.l.(ω) =
At.l. ω
tanh 1
+BL−40 ωδ(ω), (4)
At.l. =
(4π)2
, B =
dA. (5)
While the ω4 term is expected to survive in the inter-
acting theory with only logarithmic corrections, the δ-
function at the origin corresponds to the fact that gluons
http://arxiv.org/abs/0704.1801v1
mailto:meyerh@mit.edu
are asymptotic states in the free theory and implies an
infinite viscosity.
On the lattice, translations only form a discrete group,
so that a finite renormalization is necessary, Tµν(g0) =
Z(g0)T
(bare)
µν . We employ the Wilson action [9], Sg =
x,µ6=ν Tr {1 − Pµν(x)}, on an L0 · L3 hypertoroidal
lattice, and the following discretized expression of the
Euclidean energy:
(bare)
00 (x) ≡
a4g20
ReTrPkl(x)−
ReTrP0k(x)
One of the lattice sum rules [10] can be interpreted as a
non-perturbative renormalization condition for this par-
ticular discretization, from which we read off Z(g0) =
g20(cσ − cτ ). The definition of the anisotropy coeffi-
cients cσ,τ can be found in [12], where they are computed
non-perturbatively. With a precision of about 1%, a Padé
fit constrained by the one-loop result [11] yields
Z(g0) =
1− 1.0225g20 + 0.1305g40
1− 0.8557g20
, (6/g20 ≥ 5.7). (6)
Numerical results.— We report results obtained on a
β ≡ 6/g20 = 6.2, 8 · 203 lattice and on a β = 6.408, 8 · 283
lattice. The first is thus at a temperature of 1.24Tc, the
second at T = 1.65Tc. We use the results for aTc obtained
in [14] and the non-perturbative lattice β-function of [15]
to determine this. We employ the two-level algorithm
described in [16]. The computing time invested into the
1.65Tc simulation is about 860 PC days. Following [4], we
discretize 1
〈(T 11 −T 22)(T 11 −T 22)〉 instead of 〈T 12T 12〉
(the two are equal in the continuum) to write C(x0) =
〈Oη(0)Oη(x0)〉+O(a2), where
Oη(x0) ≡ 12a
{T 11 − T 22}(g0, x)
2Z(g0)
a g20
ReTr {P10 + P13 − P20 − P23}(x).
The three electric-electric, magnetic-magnetic and
electric-magnetic contributions to C(x0) are computed
separately and shown on Fig. 1. We apply the follow-
ing technique to remove the tree-level discretization er-
rors [17] separately to CBB, CEE and CEB . Firstly, x̄0
is defined such that Ct.l.cont(x̄0) = C
lat (x0). The improved
correlator is defined at a discrete set of points through
C(x̄0) = C(x0), and then augmented to a continuous
function via C(x̄
0 ) = α + βC
cont(x̄
0 ), i = 1, 2, where
0 and x̄
0 correspond to two adjacent measurements.
The resulting improved correlator, normalized by the
continuum tree-level result, is shown on Fig. 2. One ob-
serves that the deviations from the tree-level result are
surprisingly small, while deviations from conformality are
visible. The latter is not unexpected at these tempera-
tures, where p/T 4 is still strongly rising [18]. Finite-
volume effects on the T = 1.65Tc lattice are smaller than
 1000
 0  0.1  0.2  0.3  0.4  0.5
C**(x0)
FIG. 1: The correlators that contribute to C(x0) =
(CBB +
CEE + 2CEB). Filled symbols correspond to T = 1.65Tc,
open symbols to 1.24Tc. Error bars are smaller than the data
symbols.
 0.2  0.3  0.4  0.5
C(x0) / C
t.l.(x0)
1.65Tc
1.24Tc
FIG. 2: The tree-level improved correlator C(x0) normalized
to the tree-level continuum infinite-volume prediction. The
four points in each sequence are strongly correlated, but their
covariance matrix is non-singular.
one part in 103 at tree-level. Non-perturbatively, at the
same temperature with resolution L0/a = 6, increasing
L/a from 20 to 30 reduces C(L0/2) by a factor 0.922(73).
While not statistically compelling as it stands, the effect
deserves further investigation.
The entropy density is obtained from the relation
s = (ǫ + p)/T and the standard method to compute
ǫ + p ([12], Eq.1.14). We find s/T 3 = 4.72(3)(5) and
5.70(2)(6) respectively at T/Tc = 1.24 and 1.65 (the first
error is statistical and the second is the uncertainty on
Z(g0)). The Stefan-Boltzmann value is 32π
2/45 in the
continuum and 1.0867 times that value [12] at L0/a = 8.
Unsatisfactory attempts to extract the viscosity.— In
order to compare with previous studies [4, 5], we fit C(x0)
with a Breit-Wigner ansatz
ρ(ω)/ω =
1 + b2(ω − ω0)2
1 + b2(ω + ω0)2
, (7)
although it clearly ignores asymptotic freedom, which im-
plies that ρ(ω) ∼ ω4 at ω ≫ T [6]. The result of a cor-
related fit at T = 1.65Tc using the points at Tx0 = 0.5,
0.35 and 0.275 is a3F = 0.78(4), (b/a)2 = 240(30) and
aω0 = 2.36(4), and hence η/s|T=1.65Tc = 0.33(3). A
comparison of this to the results of Ref. [5] illustrates the
progress made in statistical accuracy.
An ansatz motivated by the hard-thermal-loop frame-
work is [6]
ρ(ω)/ω =
1 + b2ω2
+ θ(ω − ω1)
tanhω/4T
. (8)
It is capable of reproducing the tree-level prediction,
Eq. 4, and it allows for a thermal broadening of the delta
function at the origin. Fitting the T = 1.65Tc points
shown on Fig. 2, the χ2 is minimized for b = 0 (effec-
tively eliminating a free parameter), A/At.l. = 0.996(8),
ω1/T = 7.5(2) and η/s = 0.25(3), with χ
min = 4.0. Thus
while the ansatz is hardly compatible with the data, it
shows that the data tightly constrains the coefficient A
to assume its tree-level value.
A bound on the viscosity.— The positivity property of
ρ(ω) allows us to derive an upper bound on the viscosity,
based on the following assumptions:
1. the contribution to the correlator from ω > Λ is
correctly predicted by the tree-level formula
2. the width of any potential peak in the region ω < T
is no less than O(T ).
The standard QCD sum rule practice is to use perturba-
tion theory from the energy lying midway between the
lightest state and the first excitation. With this in mind
we choose Λ = max(1
[M2 +M2∗ ] ≈ 2.6GeV, 5T ), where
M2(∗) are the masses of the two lightest tensor glueballs.
Perturbation theory predicts a Breit-Wigner centered at
the origin of width Γ = 2γ [6], where γ ≈ αsNT is the
gluon damping rate. To derive the upper bound we con-
servatively assume that for ω <
2T , ρ(ω)/ω is a Breit-
Wigner of width Γ = T centered at the origin. From
L0) ≥ L50
ρBW (ω) +
ρt.l.(ω)
sinhωL0/2
obtain (with 90% statistical confidence level)
η/s <
0.96 (T = 1.65Tc)
1.08 (T = 1.24Tc).
The spectral function.— As illustrated above, it is rather
difficult to find a functional form for ρ(ω) that is both
physically motivated and fits the data. In a more model-
independent approach, ρ(ω) is expanded in an orthogo-
nal set of functions, which grows as the lattice resolution
on the correlator increases, and becomes complete in the
 0  5  10  15  20  25
ρ(ω) K(x0=1/2T,ω)/T
T=1.24Tc
T=1.65Tc
FIG. 3: The result for ρ(ω). The meaning of the error bands
and the curves is described in the text. The area under them
equals C(L0/2) = 8.05(31) and 9.35(42) for 1.24Tc and 1.65Tc
respectively.
limit of L0/a → ∞. We proceed to determine the func-
tion ρ̄(ω) ≡ ρ(ω)/ tanh(1
ωL0) by making the ansatz
ρ̄(ω) = m(ω) [1 + a(ω)], (10)
where m(ω) > 0 has the high-frequency behavior of
Eq. 4, and correspondingly define K̄(x0, ω) = coshω(x0−
L0)/ cosh
ωL0. Suppose that m(ω) already is a
smooth approximate solution to ρ̄(ω); inserting (10) into
(Eq. 2), one requires that a(ω) =
ℓ cℓaℓ(ω), with
{aℓ} a basis of functions which is as sensitive as pos-
sible to the discrepancy between the lattice correla-
tor and the correlator generated by m(ω). These are
the eigenfunctions of largest eigenvalue of the symmet-
ric kernel G(ω, ω′) ≡
M(x0, ω)M(x0, ω
′), where
M(x0, ω) ≡ K̄(x0, ω)m(ω). These functions satisfy
dωuℓ(ω)uℓ′(ω) = δℓℓ′ and have an increasing number
of nodes as their eigenvalue decreases. Thus the more
data points available, the larger the basis and the finer
details of the spectral function one is able to determine.
To determine the spectral function from N points of
the correlator, we proceed by first discretizing the ω vari-
able into an Nω-vector. The final spectral function is
given by the last member ρ(N) of a sequence whose first
member is ρ(0) = m and whose general member ρ(n) re-
produces n points (or linear combinations) of the lattice
correlator. For n ≥ 1, ρ(n) = ρ(n−1)[1 +
ℓ=1 c
and the functions a
ℓ (ω) are found by the SVD decom-
position [19] of the Nω × n matrix M (n)t, where M (n)ij ≡
0 , ωj)ρ̄
(n−1)(ωj). The ‘model’ m(ω) is thus up-
dated and agrees with ρ(ω) at the end of the procedure.
We first performed this procedure on coarser lattices
with L0/a = 6 at the same temperatures, starting from
m(ω) = At.l.ω
4/(tanh(1
ωL0) tanh(
ωL0) tanh
2(cωL0))
with 1
≤ c ≤ 1
, and then recycled the output as seed for
the L0/a = 8 lattices. On the latter we used the N = 4
points shown on Fig. 2.
The next question to address is the uncertainty on
ρ(ω). It is important to realize that even in the absence
of statistical errors, a systematic uncertainty subsists due
to the finite number of basis functions we can afford to
describe ρ(ω) with. A reasonable measure of this uncer-
tainty is by how much ρ(ω) varies if one doubles the re-
solution on C(x0). This can be estimated by ‘generating’
new points by using the computed ρ(N)(ω). On the other
hand we perform a two-point interpolation in x0-space
(we chose the form (α + β(x0 − 12L0)
2)/ sin5(πx0/L0)),
and take the difference between these and the generated
ones as their systematic uncertainty. In practice this dif-
ference is added in quadrature with the statistical uncer-
tainty. Next we repeat the procedure to find ρ described
above with N → 2N : if we use as seed ρ(N), then by
construction it is left invariant by the iterative proce-
dure, but the derivatives of ρ(2N) with respect to the
2N points of the correlator can be evaluated. The er-
ror on ρ(ω) is then obtained from a formula of the type
(δρ)2 =
)2(δCi)
2 which however keeps track of
correlations in x0 and Monte-Carlo time. This is the er-
ror band shown on Fig. 3 and the corresponding shear
viscosity values are
η/s =
0.134(33) (T = 1.65Tc)
0.102(56) (T = 1.24Tc).
It is also interesting to check for the stability of the so-
lution under the use of a larger basis of functions. If in-
stead of starting from ρ(N)(ω) we restart from ρ(0) (the
output of the L0/a = 6 lattice) and fit the 2N (depen-
dent) points using 2N basis functions {aℓ}, we obtain the
curves drawn on Fig. 3. As one would hope, the oscilla-
tions of ρ(2N)(ω) are covered by the error band.
Conclusion.— Using state-of-the-art lattice techniques,
we have computed the correlation functions of the energy-
momentum tensor to high accuracy in the SU(3) pure
gauge theory. We have calculated the leading high-
temperature cutoff effects and removed them from the
correlator relevant to the shear viscosity, and we nor-
malized it non-perturbatively, exploiting existing results.
We obtained the entropy density with an accuracy of 1%.
The most robust result obtained on the shear viscosity
is the upper bound Eq. (9), which comes from lumping
the area under the curve on Fig. 3 in the interval [0, 6T ]
into a peak of width Γ = T centered at the origin. Sec-
ondly, our best estimate of the shear viscosity is given by
Eq. (11), using a new method of extraction of the spectral
function. The errors contain an estimate of the systema-
tic uncertainty associated with the limited resolution in
Euclidean time. We are extending the calculation to finer
lattice spacings and larger volumes to further consolidate
our findings.
The values (11) are intriguingly close to saturating the
KSS bound [8] η/s ≥ 1/4π. We note that in perturba-
tion theory the ratio η/s does not depend strongly on the
number of quark flavors [20]. Our results thus corrobo-
rate the picture of a near-perfect fluid that has emerged
from the RHIC experiments, with the magnitude of the
anisotropic flow incompatible with η/s & 0.2 [3].
Acknowledgments.— I thank Krishna Rajagopal and
Philippe de Forcrand for their encouragement and many
useful discussions. This work was supported in part by
funds provided by the U.S. Department of Energy under
cooperative research agreement DE-FC02-94ER40818.
[1] P. F. Kolb, P. Huovinen, U. W. Heinz and H. Heiselberg,
Phys. Lett. B 500, 232 (2001); P. Huovinen, P. F. Kolb,
U. W. Heinz, P. V. Ruuskanen and S. A. Voloshin, Phys.
Lett. B 503, 58 (2001).
[2] D. Teaney, J. Lauret and E. V. Shuryak, Phys. Rev. Lett.
86, 4783 (2001).
[3] D. Teaney, Phys. Rev. C 68, 034913 (2003).
[4] F. Karsch and H.W.Wyld, Phys. Rev. D 35, 2518 (1987).
[5] A. Nakamura and S. Sakai, Phys. Rev. Lett. 94, 072305
(2005).
[6] G. Aarts and J.M. Martinez Resco, JHEP 0204, 053
(2002).
[7] G. Policastro, D.T. Son and A.O. Starinets, Phys. Rev.
Lett. 87, 081601 (2001).
[8] P. Kovtun, D.T. Son and A.O. Starinets, Phys. Rev. Lett.
94, 111601 (2005).
[9] K.G. Wilson, Phys. Rev. D 10, 2445 (1974).
[10] C. Michael, Phys. Rev. D 53 (1996) 4102.
[11] F. Karsch, Nucl. Phys. B 205, 285 (1982).
[12] J. Engels, F. Karsch and T. Scheideler, Nucl. Phys. B
564 (2000) 303 [arXiv:hep-lat/9905002].
[13] A. Hosoya, M.A. Sakagami and M. Takao, Annals Phys.
154, 229 (1984).
[14] B. Lucini, M. Teper and U. Wenger, JHEP 0401, 061
(2004).
[15] S. Necco and R. Sommer, Nucl. Phys. B 622 (2002) 328.
[16] H.B. Meyer, JHEP 0401, 030 (2004).
[17] R. Sommer, Nucl. Phys. B 411 (1994) 839.
[18] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Leg-
eland, M. Lutgemeier and B. Petersson, Nucl. Phys. B
469, 419 (1996).
[19] R.K. Bryan, Eur. Biophys. J., 18 (1990) 165.
[20] P. Arnold, G.D. Moore and L.G. Yaffe, JHEP 0305, 051
(2003).
http://arxiv.org/abs/hep-lat/9905002
