The Bulk Channel in Thermal Gauge Theories

We investigate the thermal correlator of the trace of the energy-momentum tensor in the SU(3) Yang-Mills theory. Our goal is to constrain the spectral function in that channel, whose low-frequency part determines the bulk viscosity. We focus on the thermal modification of the spectral function, $\rho(\omega,T)-\rho(\omega,0)$. Using the operator-product expansion we give the high-frequency behavior of this difference in terms of thermodynamic potentials. We take into account the presence of an exact delta function located at the origin, which had been missed in previous analyses. We then combine the bulk sum rule and a Monte-Carlo evaluation of the Euclidean correlator to determine the intervals of frequency where the spectral density is enhanced or depleted by thermal effects. We find evidence that the thermal spectral density is non-zero for frequencies below the scalar glueball mass $m$ and is significantly depleted for $m\lesssim\omega\lesssim 3m$.


Introduction
Understanding the emergence of bulk properties and collective behavior from the strong interaction of particle physics is a goal that has been pursued vigorously in recent years. Experimentally, heavy-ion collisions have provided a way to heat up nuclear matter to temperatures at which quarks and gluons are expected to dominate its macroscopic properties [1,2,3,4]. Theoretically, both analytic [5] and numerical methods [6] are being used to unravel the intricate dynamics of QCD at finite temperature and baryon density.
At temperatures of a few hundred MeV and small baryon density, lattice QCD remains the tool of choice to determine the equilibrium properties of the system. As we shall discuss in detail, access to the near-equilibrium properties, such as the transport properties, is limited in this approach, because lattice QCD employs the Euclidean formulation of thermal field theory. Real-time properties can thus only be determined by analytic continuation, which for numerical purposes represents an ill-posed problem. However, it is not excluded that we can learn about the gross features of the spectral density by suitably combining numerical and analytic methods. The purpose of this paper is to illustrate how such a combination of methods can be put to practice, taking as an example the bulk channel in the SU(3) Yang-Mills theory.
One of the lasting legacies of the Relativistic Heavy Ion Collider (RHIC) experiments is the success of the hydrodynamic description of the collisions, see for instance the review [7]. The early agreement between ideal hydrodynamics and experiment was gradually refined to incorporate the dissipative effects of shear viscosity η, the three-dimensional geometry [8,9,10] and to estimate the sensitivity to initial conditions [11,12]. Most recently, the effects of bulk viscosity ζ have been investigated [13,14]. At temperatures sufficiently high for weak coupling methods to be applicable, the bulk viscosity is three orders of magnitude smaller than the shear viscosity [15]. However, by analogy with many condensed matter systems (see for instance [16]), where a maxixum in the bulk viscosity is observed near liquidgas phase transitions, one generally expects the bulk viscosity to become appreciable, if anywhere, near the rapid crossover between the confined and the deconfined phases. It has been suggested that this sudden onset of dissipation under expansion could trigger the freeze-out in heavy-ion collisions due to cavitation [14]; the breakup of a fluid occurs when a component of the stress-energy tensor, say T zz , turns negative. The behavior of ζ near the QCD crossover was studied about two-and-a-half years ago, first using a QCD sum rule and a crude model for the bulk spectral function [17,18], and subsequently using Euclidean correlators obtained by Monte-Carlo methods [19]. Both types of analyses missed the fact that a homogeneous fluctuation in T µµ has a component which takes arbitrarily long to relax to equilibrium. It is therefore necessary to revisit these analyses, and this constitutes the second motivation for this work. A study based on AdS/CFT methods [20] also found an increase in the bulk viscosity near T c , however less pronounced than suggested by [17,19].
In the near future, the low-energy runs at RHIC, and further ahead, the heavy-ion experiments at FAIR will explore the QCD diagram at larger net baryon density. It is widely thought, although not beyond doubt, that a chiral critical point exists in the T vs. µ B plane [21]. The associated critical exponents [22] indicate that the bulk viscosity should exhibit a strong divergence near the critical point. If this increase in ζ is not too closely localized around T c , it could be one of the signatures of the critical point.
The outline of this paper is as follows. We start in section 2 with a review of retarded and Euclidean correlators, establishing the precise relations between them. Details are collected in the appendix. In section (3) we turn to the specific case of the bulk channel in the SU(N ) gauge theory. We describe how to use the bulk sum rule as a constraint on the vacuum-subtracted spectral density ∆ρ, and present the operator product expansion (OPE) prediction for the large-frequency behavior of the spectral density. In section 4 we give a numerical application of our results for N = 3. We first present data deep in the confined phase, as well as in the deconfined phase. The sum rule is combined with lattice data to yield information on the changes in the spectral density that occur between the confined and the deconfined phases. We end with the conclusions of this study.

A review of finite-temperature correlators
In this section, we give the definitions of the retarded correlator and the Euclidean correlator, and show their relation in detail. This allows us to fix our conventions, but also to highlight a subtlety that arises in the zero-frequency case. The retarded correlator, evaluated at a purely imaginary frequency, is the analytic continuation of the Euclidean correlator, which is only available at the Matsubara frequencies. We show that the odd derivatives of the retarded correlator along the imaginary axis, evaluated at the origin, coincide with the derivatives of the spectral density at the origin. The even derivatives, on the other hand, correspond to integrals over all frequencies of the spectral density. We illustrate the relations that we derive on the example of the bulk channel in SU(N ) gauge theory. Working in frequency space rather than in time-coordinate space, we show the relation between the bulk viscosity and the frequency-space Euclidean correlator, Fig. (1).
We work in the canonical ensemble, i.e. expectation values of operators are defined by where the partition function Z(β) is such that 1 = 1 and β = 1/T is the inverse temperature. In the following we consider a generic Hermitian operator O. We use the Heisenberg representation, withÔ We start by introducing the real-time correlators, Its Fourier transformG plays an important role in the following. With these conventions we note thatG R is analytic in the half complex plane Im(ω) > 0. The Euclidean correlator is defined as It obeys the Kubo-Martin-Schwinger relation and can be expressed as a Fourier series on the interval 0 ≤ t < β, where The spectral representation of the correlators is useful to prove simple relations between the Euclidean and real-time correlators. One finds that (see the appendix) The relation of the Matsubara zero-mode ℓ = 0 with the zero-frequency retarded correlator requires special care. We shall see that their difference can be expressed in terms of the spectral function.
The spectral function is defined as For real ω and in finite volume, ρ is a distribution, namely a discrete sum of delta functions. However, in the infinite-volume limit of an interacting theory, it is typically a smooth function, except possibly at a finite number of points. The most useful relations among the correlators are those that survive the infinite-volume limit. In the appendix, we show that In the limit ω I → 0, only the integration region around the origin contributes on the RHS. If the spectral function is finite at the origin, the RHS is O(ω I ), and therefore lim ǫ→0GR (iǫ) = G E . On the other hand, if the spectral function contains a delta-function at the origin, ρ(ω)/ω = Aδ(ω), lim ǫ→0GR (iǫ) = G (0) E − A. We have written Eq. (2.11) with an explicit UV cutoff to make it clear that one should take the limit ω I → 0 first.
The connection between the retarded and Euclidean correlators is thus precisely established in frequency space. In appendix we also rederive the relation between them in coordinate space (well-known in particular to lattice practitioners), which we shall use in section (4).

The retarded correlator along the imaginary axis
If we assume for a moment that converges in the UV, a simple calculation based on the explicit expression for the spectral density (A.13), or alternatively a contour integration, then leads to the Kramers-Kronig relation ∀ω I > 0. (2.14) In general, it is necessary to perform a subtraction from ρ to make the dispersion integral converge in quantum field theory. Starting from (2.14), it is easy to derive by recursion the following equalities, Thus at the origin, the even derivatives of the retarded correlator along the imaginary axis are given by integrals over the spectral density, while the odd derivatives are equal to the corresponding derivatives of the spectral function.

Remarks on the Euclidean correlator
Combining Eq. (2.14) with Eq. (2.11), one obtains Equations (2.11) and (2.17) show that while G R (iω I ) is insensitive to the presence of a ωδ(ω) term in the spectral function for any ω I > 0, G E definitely receives a contribution from it.
Secondly we remark that the positivity of ρ(ω)/ω implies the positivity of certain linear combinations of the Euclidean correlator. For instance, . (2.19)

An example
Here we anticipate the results of section (3) to illustrate the general formulas of this section. In Figure (1), we sketch a possible form of the retarded correlator for the operator T µµ projected onto zero spatial momentum. In order to deal with a UV finite function, we subtract the corresponding zero-temperature correlator, and leave the temperature-dependence implicit. This subtraction is sufficient to make the dispersion integral (2.14) converge. As we shall see in section (3.1), the spectral function of the operator d 3 x T µµ (x) admits a delta function located at the origin, with a weight given by a thermodynamic quantity. Applying Eq. (2.11), we obtain ∆G R (iǫ) = ∆G Furthermore an exact sum rule fixes ∆G (0) E in terms of thermodynamic functions, and ∆G R (iǫ) turns out to be negative. Using the operator-product expansion (section 3.2), one can show that ∆G (ℓ) E is negative at large Matsubara frequencies, implying the same for ∆G R (iω) at large ω in view of Eq. (2.9). Finally, Eq. (2.16) and (3.3) show that (d/dω)∆G R (iω) at the origin is proportional to the bulk viscosity, which has to be positive by the second law of thermodynamics. All this information is summarized in Fig. (1). The behavior at intermediate frequencies ω is not known to us, in principle ∆G R (iω) could even change sign in view of the subtraction made. It is also not known whether the plotted function ever exceeds unity for ω a multiple of 2πT .
In a numerical Euclidean approach, the available information isG R (iω ℓ ), ℓ = 0, 1, . . . The low-frequency real-time properties of the system, such as the bulk viscosity and the associated relaxation time, are determined by the odd derivatives ofG R (iω) at the origin, see Eq. (2.16). The latter are difficult to determine, because ω = 0 is a non-analytic point ofG R . In this respect, Eq. (2.16) shows that if the spectral density has a gap m, as it does at zero temperature in a confining theory, thenG R (iω) can be expanded in integer powers of ω 2 . Therefore subtracting the zero-temperature correlator only affects the even derivatives of ∆G R (iω) at the origin.
Working in frequency space has the advantage that the high-frequencies are decoupled from the discussion at this stage. There are, however, technical difficulties in determining Euclidean correlators in frequency space. Unlike the t-dependent correlators, they are off-shell correlators, and therefore require additive counterterms in order to approach the continuum limit with corrections that are suppressed with a power of the lattice spacing. An example of this complication arises in the derivation of sum rules [23]. In this special case the sum rule precisely dictates what the counterterms are. In general, renormalization conditions have to be imposed in order to determine them.
Finally we make the observation that a polynomial of degree n in ω going through the n + 1 points , has a slope at the origin given by When expressed in terms of dispersion integrals as in Eq. (2. 18-2.19), the degree of convergence in the UV of this expression improves with the order of the polynomial: with some combinatorics one can show that it is asymptotically ρ(ω)/ω n+3 for n even and ρ(ω)/ω n+2 for n odd. In particular, the slope at the origin of the n = 2 polynomial is given by where the linear combination (2.19) appears. In the specific case of the bulk channel, the spectral integral (2.22) is UV-convergent without any subtraction to ρ(ω). Thus the slope of the quadratic polynomial is negative-definite, just as the actual slope of the retarded correlator d dωG R (iω). One might be tempted to use it as a Euclidean estimator for the bulk viscosity ζ. However, if the spectral function has a transport peak at the origin of width ω 0 ω M , expression (2.22) amounts parametrically to ζω 0 /ω M , while d dωG R (iω) ∼ ζ. It is clear why that is: a polynomial interpolating between the Euclidean points is only a good approximation to the retarded correlator if the latter does not contain a scale ω 0 ≪ ω M over which it varies significantly.

The reconstructed Euclidean correlator
When analyzing Euclidean correlators obtained from Monte-Carlo simulations, it is instructive to study in what way the Euclidean correlators differ from what they would be if the spectral density was unaffected by thermal effects [24]: allows us to derive in particular the exact relations These expressions are useful in practice when one is interested in a difference of spectral functions. We will exploit Eq. (2.25). Even if the zero-temperature Euclidean correlator can usually not be determined for arbitrarily large separations, in that regime it is completely dominated by the lightest state in the channel, so that its functional form is known and completely determined by one matrix element and an energy. The contribution to G rec E of a single state of energy E is easily obtained, Secondly, the correlator falls off very rapidly, so that the terms with large values of m in Eq. (2.25) make a very small contribution.

The bulk channel in QCD
In this section the goal is to provide the theoretical background for the analysis of the bulk channel in SU(N ) gauge theory. Thus we consider the operator is the energymomentum tensor. Its spectral function ρ θ (ω, T ) and the difference play an important role in the following. The spectral function determines the bulk viscosity through the Kubo formula In what follows we review the behavior that hydrodynamics predicts for the lowfrequency part of the associated spectral function ρ θ . It turns out that it contains a term ωδ(ω). We then define an 'auxiliary' operator whose spectral function coincides with ρ θ except for not having the delta function at the origin. Next we present the asymptotic large-ω behavior of ρ θ based on the operator-product expansion. Finally, the spectral function of the auxiliary operator obeys a sum rule given in section (3.2).

Low-frequency prediction from hydrodynamics
In this section we consider more general connected correlation functions of the energymomentum tensor, as well as their associated spectral function ρ µνρσ (ω, k, T ) via (2.12). The spectral function for T 33 (assuming the spatial momentum k is pointing in the direction3) reads [25] (3.5) Now one shoud recall that for any smooth function such that ∞ −∞ f (x)dx = 1, 1 ǫ f (x/ǫ) provides a representation of the delta function. Here we will exploit this fact for the function and b = c s /(Γ s k) and ǫ = Γ s k 2 . In this way one finds that, in the sense of distributions, This is in contrast with what is obtained if one simply sets k = 0 in Eq. (3.5), in which case one misses the delta function. In the shear channel on the other hand, and the delta function is suppressed by a power of k 2 . Thus combining the shear and sound channels, Using the exact relations When one constrains the spectral function using sum rules and Euclidean correlators, one is always dealing with integrals of the spectral function over all frequencies. Therefore it is useful to subtract the contribution of the delta function in (3.12). One way to achieve this is to work with the operator Its spectral function, which we denote by ρ ⋆ , satisfies and, in view of Eq. (3.12), is free of the delta function at the origin. At equilibrium, recalling that c 2 and therefore its linear fluctuations are not correlated with those of the total energy 1 . This is the physical reason why a delta function at the origin is absent in its spectral function.

Large-frequency behavior from the operator product expansion
In section (3.1) we have analyzed the low-frequency behavior of the bulk spectral function using the general hydrodynamic predictions. We now turn to an analysis of its asymptotic high-frequency behavior. The leading ρ θ ∼ α 2 s (ω)ω 4 behavior was given for instance in [26]. In the OPE language, this contribution corresponds to the insertion of the unit operator, and it is therefore temperature independent. Because it makes a large contribution to the coordinate-space Euclidean correlator, it is best to remove it by subtracting the vacuum correlator from the thermal correlator. We therefore need to know what the leading behavior of the subtracted correlator and of the corresponding subtracted spectral function is.
The operator product expansion of local gluonic operators was calculated a long time ago in the context of QCD sum rules [27], and the results of more recent calculations are summarized in [28]. In many cases, the Wilson coefficients are known at next-to-leading order. However, only the Lorentz scalar operators were kept on the right-hand side of the OPE, since the correlator was evaluated in the vacuum. Since we want to apply the OPE to a system at finite-temperature, the traceless part of the energy-momentum tensor, for instance, can also contribute. The OPE results therefore need to be generalized to the finite-temperature context [29]. We have reproduced the result of [29] for the bulk channel using the background field method combined with the Fock-Schwinger gauge [30], We remark that while the two-point function ∆G R (iω) of the operator T µν − 1 4 δ µν θ contains a logarithmic divergence (which manifests itself in the Wilson coefficient of θ being O(1/g 2 0 ∼ log 1/a) in lattice regularization [23]), ∆G R (iω) is UV-finite for the operator θ. Therefore, assuming the validity of the operator-product expansion, the only relevant scale in its Wilson coefficients is ω, and the dispersion relation (2.14) then leads to Eq. (3.16) [31]. In section (4), we will use this information to constrain the form of the spectral function at large frequencies.

The Bulk Sum Rule
The idea of constraining the form of finite-temperature spectral functions with sum rules is not new [32,31,17]. Here we focus on the bulk channel. In Euclidean space, the following sum rule holds [33], where e and p are the energy density and pressure of the finite-temperature system. We now convert identity (3.17) into a sum rule for the bulk spectral function with the goal to constrain the latter. Combining equations (3.17) and (2.17), we obtain Next, we separate the contribution of the ωδ(ω) term in ρ θ from the rest of the spectral function. In other words, we want to turn Eq. (3.18) into a sum rule for ρ ⋆ , the spectral function of the operator O ⋆ (3.13), which is smooth at the origin. Using Eq. (3.14), one finds 2 A few remarks are in order. It would have been less straightforward to directly derive a sum rule for ρ ⋆ , because the Euclidean correlator of O ⋆ contains a short-distance singularity. We remark that Eq. (3.18) was already obtained in [17], but the presence of the ωδ(ω) term in ρ θ was missed. In [34], a finite spatial momentum k was used as an infrared regulator and the contribution of the zero-frequency delta function correctly identified for the first time. If one insists on expressing Eq. (3.19) in terms of ρ θ , one should strictly write lim ǫ→0

Numerical results from the lattice
In this section we describe the lattice setup and the numerical results obtained by Monte-Carlo simulations. We employ the isotropic Wilson action [35], where the 'plaquette' P µν is the product of four link variables U µ (x) ∈ SU(3) around an elementary cell in the (µ, ν) plane. As a discretization of L 3/2 O θ , we use up to an irrelevant additive constant (its contribution vanishes in connected correlation functions). As a local update algorithm, we use the standard combination of heatbath and over-relaxation [36,37,38,39] sweeps in a ratio increasing from 3 to 5 as the lattice spacing is decreased. The multi-level algorithm is used to reduce the variance of the correlator [40,41]. To set the scale and to evaluate the lattice beta function dg −2 0 /d log a we use the parametrization of the Sommer scale r 0 ≈ 0.5fm given in [42]. We denote the mass gap by m, given by the lightest scalar glueball, whose mass is 3.96(5)/r 0 = m = 5.30(8)T c . Finally we relate the notation used in the rest of this section to the general definitions of sections (2) and (3), The quantities with a ⋆ subscript are associated with the operator O ⋆ defined in Eq. (3.13).
The Euclidean correlators C ... are dimensionless, renormalization-group invariant quantities. The goal of this section is to constrain the function ∆ρ ⋆ , the main results being Eq. (4.23) and Eq. (4.24). We start by describing a method we have used to evaluate the right-hand side of the bulk sum rule (3.19). We also evaluate the 'condensates' appearing in the OPE (3.16) of the bulk correlator. We obtain C ⋆ at finite temperature by subtracting the contribution of the δ function. We then describe the properties of the correlator at temperatures deep in the confined phase. This correlator is subsequently used to produce a 'reconstructed' correlator based on (2.25). From then on we are directly dealing with integrals of the subtracted spectral function, ∆ρ ⋆ . We reuse lattice data presented in [19], and combine it with the bulk sum rule to determine the gross features of the bulk spectral function.

Right-hand side of the bulk sum rule
When using the bulk sum rule (3.19), the question arises, how to best evaluate the combination of thermodynamic quantities appearing on the right-hand side of the equation. Note that the enthalpy e + p and the conformality measure e − 3p are relatively straightforward to obtain from the expectation value of the electric and magnetic field strengths. What is then needed is a reliable way to obtain the velocity of sound c 2 s , or alternatively, the specific heat c v .
The Euclidean sum rule (3.17) can be derived on the lattice. Here we will use that lattice sum rule [23] to determine the velocity of sound via where the discretization Θ of the trace anomaly operator is given by Eq. (4.2). We used the leading perturbative expression for the expression inside the large brackets, since based on that estimate it only departs from 4 by 0.008 or less. The right-hand side of the bulk sum rule (3.19), as calculated on N τ = 6 lattices, is displayed in units of e + p in Fig. (2) as a function of temperature. It is negative for a significant range of temperatures above T c , and its leading behavior in perturbation theory is [43] 3 It is therefore likely to be negative throughout the deconfined phase. We also note that for α s = 0.25, the right-hand side of (4.7) amounts to about -0.03, a very small value compared to those appearing in Fig. (2).

Evaluation of the correlator C ⋆ (t, T )
Since we want to work with the correlation functions of the operator O ⋆ , we have to evaluate its Euclidean correlation function C ⋆ (t). Lattice perturbation theory [44] and numerical studies [19] show that it is advantageous to work with the operator O θ . To convert its correlation function to the correlation function C ⋆ , we use the equation (see 3.14 and 2.12) Here too a combination of thermodynamic quantities appears. We have obtained it by the same method as the right-hand side of the bulk sum rule, see the previous subsection. Close to the phase transition where c 2 s is very small, the result is rather sensitive to the value of c 2 s , since it appears in the denominator. That is the main reason we opted for a direct calculation of c 2 s from the lattice sum rule (4.6), rather than parametrizing the pressure and obtaining c 2 s by taking derivatives. In this way we avoid the dependence on a parametrization ansatz, which can be significant in a region where the derivatives are very   4.9) and the difference between C θ and C ⋆ is parametrically negligible, since both are O(α 2 s ). The resulting correlators are displayed in Fig. (3). The lattice data for C θ is the same as in [19]. In particular, the lattice spacing is given by 1/aT = 8. At high temperatures, there is very little numerical difference between them, confirming the parametric expectation (4.9). On the other hand we observe that the difference between C θ and C ⋆ becomes appreciable very close to T c . At that point, C ⋆ is significantly smaller than C θ .

Low-temperature correlators
The C θ correlator at about T c /2 is displayed in Fig. (4). The agreement between the data obtained at two different lattice spacings suggests that the cutoff effects on the correlators are small. The dotted curve at small separation, which represents the leading perturbative prediction for the choice α s = 0.30, shows that by comparison the non-perturbative correlator is rather large. The right panel shows that starting at about 0.4fm, the lightest glueball represents a large contribution of the full correlator. The data is not accurate enough to tell how well that state saturates the correlator beyon 0.6fm.
Although the finite size of the lattice implies that the temperature is about T c /2, rather than zero, it is very likely that the spectral function is already close to the zerotemperature spectral function, albeit in finite volume. This statement can be justified as follows. The confined phase can be thought of as a dilute gas of glueballs propagating freely, up to occasional collisions. Evidence in favor of this picture comes from a precision thermodynamic calculation in the confined phase [45]. Thus, to leading order at low temperatures, where m is the mass gap. Based on these approximate expressions, we find that the contribution of the exact delta function to the correlator C θ at t = β/2 is (4.12) Numerically, at T = T c /2 this evaluates to about 0.006, while C θ (T, β/2) = 50 (11). The area under the transport peak can be estimated using kinetic theory. The general expression is given in [25], where Λ is much larger than the inverse relaxation time, but lower than the thermal scale. At low temperatures this contribution is thus parametrically smaller than the contribution of the exact delta function, and it is indeed less than 15% of the latter at T c /2. In summary, there are strong reasons to believe that the bulk spectral function at T = T c /2 is close to the zero-temperature spectral function. Therefore it is legitimate to use these Euclidean correlators to obtain the reconstructed correlators at temperatures above T c , where they can be compared with the thermal correlators.

Evaluation of the reconstructed correlator
The 'reconstructed' Euclidean correlator is by definition what the Euclidean correlator would be if the spectral function remained unchanged with respect to the in vacuo spectral function. We want to evaluate this reconstructed Euclidean correlator (Eq. 2.23) in order to compare it with, and eventually subtract it from the thermal correlator. This evaluation is simplified by the fact that the θ correlator falls off very fast at zero temperature. A first approximation is thus provided by (see Eq. 2.25) 14) The state of smallest energy contributing to the zero-temperature correlator is the lightest glueball. Its mass and the matrix element F S ≡ vac|O θ |glueball have been calculated on the lattice (in this definition the norm of the states are set to unity). When βm is large, we can therefore estimate the leading correction to Eq. (4.14), As these examples show, (4.15) is a small correction to (4.14), in fact smaller than the latter's statistical error. In subsequent figures we will therefore display (4.14) as a simple estimate of the reconstructed correlator. We note that expression (4.14) is a strict lower bound on the reconstructed correlator.

Large-frequency contributions to spectral sums
For a large value of ω 2 , Eq. (3.16) implies to leading order, For instance, if we choose ω 2 = 12/r 0 ≈ 4.8GeV, then α s (ω 2 ) ≃ 0.14 [46], and the second term on the RHS of Eq. The contribution of the tail of the spectral function to the Euclidean correlator at t = β/2 is given in leading approximation by At T = 2T c , the exponential suppression is quite strong, e −ω 2 /2T ≈ 0.018 and ∆ρ ⋆ (ω 2 ) itself is of order (-0.02). This makes the large-frequency contribution to the Euclidean correlator negligible, and in any case much smaller than its current statistical error.

4.6
The temperature dependence of the bulk spectral function Fig. (6) compares the thermal Euclidean correlators C θ (t, T ) and C ⋆ (t, T ) with the reconstructed correlator C rec ⋆ (t, T ; T 0 ). We see again that the difference between C θ and C ⋆ is small, except close to T c . However, the difference between the reconstructed correlator and C ⋆ is small too. Therefore the t-independent shift between C ⋆ and C θ is nonetheless crucial to correctly evaluate C ⋆ (t, T ) − C rec ⋆ (t, T ). From Fig. (6), we learn in particular that the central value of this difference is positive at t = β/2. While on a finer lattice one could reliably exploit more data points, here we will restrict ourselves to this middle-point, which carries the most information on the low-frequency part of ∆ρ ⋆ . Numerically we find, in terms of the subtracted spectral function (see Eq. (2.12) and Eq. (2.25)), (4.23) On the other hand, the bulk sum rule (3.19) tells us that (4.24) These two equations suggest an enhancement of the thermal spectral weight at low frequencies and a suppression at the higher frequencies. Indeed the higher frequencies contribute much more to the bulk sum rule than to Eq. (4.23).
As a technical remark we note that since a large cancellation takes place when subtracting the reconstructed correlator, the quantity X(T ) is susceptible to a large discretization error in addition to the quoted statistical error. In obtaining the results (4.23), we have used treelevel-improvement [44] separately for the finite-temperature and the zero-temperature correlators before calculating X(T ). In [19] it was shown that after treelevel-improvement the discretization errors affecting C θ (t = β/2, T ) are mild. Also, we have ignored possible finite-volume effects, which could potentially induce a large relative error in X(T ). New data on finer and physically larger lattices will be necessary to obtain X(T ) with completely controlled systematic errors.

Thermal rearrangement of the spectral weight
We have collected a certain amount of information on the thermal spectral function relative to its zero-temperature counterpart. Firstly, we know the asymptotic large-frequency behavior, thanks to the operator-product expansion. Secondly, two integrals over all frequencies of the subtracted spectral function are known. We also know that the difference ∆ρ ⋆ is positive for ω below the mass gap m, since the zero-temperature spectral function vanishes in that region, and we know the position and strength of the glueball peak.
When determining properties of the spectral function, it is technically advantageous to separate the exact delta functions it contains from the contributions which are smooth in the infinite volume limit. Identifying the delta function at the origin in the bulk spectral function was an important step in this direction. Since we are interested in low-energy properties, it is also advantageous to work with a subtracted spectral function, such that the difference goes to zero at high frequencies. This is achieved by subtracting the zerotemperature spectral function. However, doing so one introduces, with negative weight, delta functions that correspond to the stable one-particle states at zero-temperature [47]. In the pure SU(N ) gauge theory, these are the glueballs. The position and weight of these delta functions can in principle be determined. They have been determined for the lightest glueball [48], but only rough estimates are known for the second (and last) stable scalar glueball. Therefore we will not attempt to systematically isolate the stable-glueball contributions in the spectral integrals X(T ) and Y (T ), but it is instructive to estimate the negative contribution of the lightest scalar glueball, These contributions are quite large compared to X(T ) and Y (T ) themselves. We have seen that Y (T ) is negative. Based on the OPE, only a small fraction of its magnitude can be explained by a depletion of the spectral density for ω > ω 2 = 12/r 0 ≈ 4.8GeV. Therefore the bulk of the contribution to Y (T ) must come from the region m ≤ ω ≤ ω 2 . The depletion of the spectral density in the region of those glueballs that are stable or lie slightly above the two-particle threshold must therefore be the explanation for the sign and magnitude of Y (T ). Whether the glueballs are completely 'melted' in the deconfined phase, in the sense that the spectral function shows no narrow peaks in this region, cannot be answered on the basis of the available information. While no detailed properties of ∆ρ ⋆ can be determined based on the limited information at hand, we want to show somewhat more quantitatively how Eq. (4.23) and (4.24) dictate the rearrangement of the spectral weight. For instance, in an interval of frequencies, one may ask whether ω 2 ω 1 dω ω ∆ρ ⋆ is positive or negative. In order to address this question in a systematic way, we apply the following procedure (it is of course not unique). We select a binning in the frequency variable ω. We choose [0, ω 1 ] and [ω 1 , ω 2 ]. The upper end of the second bin, ω 2 , is chosen large enough for the OPE to be reliable. Since ρ ⋆ (ω, T = 0) vanishes for ω < m, it is of physical interest to choose ω 1 on the order of m.
Second, since we expect the (infinite-volume) spectral function to be smooth at least at high-frequencies, we parametrize ∆ρ ⋆ by a second-order polynomial in ω 2 in each bin. We require continuity and differentiability at the boundaries between the bins. For ω > ω 2 , we use the expression provided by the OPE. This then leaves 6 coefficients to be determined. Four equations are provided by the matching conditions at the edges of the bins, and two are provided by the known frequency integrals of ∆ρ ⋆ , Eq. (4.23) and (4.24). Determining this simple parametrization of the spectral function thus amounts, in the present case, to inverting a 6 × 6 matrix.
The result of this procedure is show in Fig. (5). The bands indicate the statistical error obtained by standard error propagation. These curves show qualitatively what was expected from the remarks below Eq. (4.24), namely the appearance of a positive spectral weight for ω < m, and a depletion relative to the vacuum spectral density for ω > m. Incidentally, this is qualitatively consistent with the prediction of perturbation theory, where the O(α 2 s ) low-frequency and O(α 2 s ) high-frequency contributions cancel eachother, resulting in a spectral sum dω∆ρ ⋆ /ω that is O(α 3 s ) [34,29]. We note that the values of the intercept in Fig. (5) at 1.24T c and 1.65T c , although strongly dependent on the fact that we have smoothened the spectral function, are in good agreement with our previous estimates for ζ/s [19], where the unsubtracted spectral function ρ θ was directly parametrized. This may well be because in both cases we made a smooth ansatz for the spectral function, but in view of the large cancellation that occurs in ∆ρ ⋆ , the consistency is not totally trivial.

Conclusion
We have collected information, obtained by analytic methods, on the spectral density of the bulk channel in SU(N ) gauge theory. We have reviewed the precise connection between the Euclidean and the retarded correlator in frequency space, illustrated in Fig. (1), recalculated the OPE of the bulk channel correlator, Eq. (3.16), and presented a new way, based on (Eq. 2.25), to compute the reconstructed correlator. The latter technique allows one to directly constrain the difference of thermal to vacuum spectral densities, ∆ρ.
Numerically, we have found that in the deconfined phase of SU(3) gauge theory, a significant depletion of the spectral density must take place in the interval m ω 3m, while we have presented evidence that a non-vanishing spectral density appears at frequencies below the mass gap, obviously a finite-temperature effect. This qualitative conclusion comes from combining the bulk sum rule with numerical lattice correlators. With numerical data of higher accuracy, it would be possible to quantify more precisely how much spectral weight is clustered around the origin with frequencies ω T .
Unfortunately, the behavior of the bulk spectral function at low frequencies near the deconfining temperature remains largely undetermined. The Euclidean correlator C(t, T ≈ T c ) clearly exhibits an enhancement compared to higher temperatures, even after subtraction of the t-independent contribution, see Eq. (4.8), which is quite large just above the transition. However this enhancement can be explained to a large extent by the properties of the vacuum, as the comparison with the reconstructed correlator shows. At the moment we can state that a depletion of the spectral density has to take place for frequencies above the mass gap, in view of the bulk sum rule, but are unable at present to give a useful estimate of the spectral weight below the mass gap. This more careful reanalysis of the bulk channel has taught us that the effects of a potentially large bulk viscosity near T c are more subtle to detect in spectral integrals than previous work suggested [17,19]. More accurate lattice data, in particular for the speed of sound near the phase transition, would allow one to determine the cumulated spectral weight below the mass gap. We believe that this goal can be achieved in the foreseeable future.

A. Spectral representation of finite-temperature correlators
In this appendix we give the spectral representation of various correlators. We use this representation to obtain Eq. (2.11).

A.1 Real-time correlators
Inserting two complete sets of energy eigenstates in the expression (2.3) for G R (t), one obtains G R (t) = iθ(t) Z m,n |O nm | 2 e −βEn (e iEnmt − e −iEnmt ) , (A.1) where we employ the notation In the same representation, the Fourier transform (2.4) reads, for Im ω > 0,

A.2 Euclidean correlators
The spectral representation of the Euclidean correlator (2.5) reads G E (t) = 1 Z n,m |O nm | 2 e −βEn e Enmt , (A.7) and after some algebraic manipulations, one also finds for the Fourier coefficients (2.8) In the next section, we shall see that the right-hand side can be expressed in terms of the spectral function.