Charm mass effects in bulk channel correlations
aa r X i v : . [ h e p - ph ] N ov Charm mass effects in bulk channel correlations
Y. Burnier a and M. Laine ba Institute of Theoretical Physics, EPFL, CH-1015 Lausanne, Switzerland b Institute for Theoretical Physics, Albert Einstein Center, University of Bern,Sidlerstrasse 5, CH-3012 Bern, Switzerland
Abstract
The bulk viscosity of thermalized QCD matter at temperatures above a few hundred MeVcould be significantly influenced by charm quarks because their contribution arises four per-turbative orders before purely gluonic effects. In an attempt to clarify the challenges of alattice study, we determine the relevant imaginary-time correlator (of massive scalar densities)up to NLO in perturbation theory, and compare with existing data. We find discrepanciesmuch larger than in the vector channel; this may hint, apart from the importance of takinga continuum limit, to larger non-perturbative effects in the scalar channel. We also recallhow a transport peak related to the scalar density spectral function encodes non-perturbativeinformation concerning the charm quark chemical equilibration rate close to equilibrium.October 2013 . Introduction
Viscosities play an important role in the hydrodynamics of finite-size systems, such as thosegenerated in heavy ion collision experiments. In contrast to thermodynamic functions like thepressure or energy density, the dominant contributions to them arise from the slowest (mostweakly interacting) processes relevant for equilibrating energy and momentum flows. Thismay lead to counter-intuitive results; for instance, in cosmology, neutrinos or dark matterparticles could play a dominant role for determining viscosities of the cosmic fluid [1, 2].In this paper we are concerned with the bulk viscosity of a QCD plasma similar to thatgenerated in heavy ion collision experiments [3]. Although less prominent than shear viscosity,it also affects the hydrodynamics of the system in an interesting way: indeed the bulk viscositycould grow rapidly as the temperature decreases below the QCD crossover [4, 5] and hasthen been speculated to contribute to “clusterization” [6]–[9] which might be viewed as athermodynamical precursor to the chemical freezeout process referred to as hadronization.The physical processes relevant for the bulk viscosity are those associated with the breakingof scale invariance. We suspect that at temperatures
T > ∼
300 MeV a significant contributionmay be given by massive quarks, in particular charm quarks. Despite enhancement fac-tors [10] the charm quarks are unlikely to reach chemical equilibrium within the lifetime ofheavy ion collisions; however it is expected that they depart from equilibrium on the side ofbeing too many (cf. e.g. ref. [11]). As suggested by experiments [12, 13] and theoretical de-terminations [14]–[18] (similar conclusions have also been reached through phenomenologicalstudies [19]–[23] as well as investigations in other gauge theories [24]–[26]), they also rapidlyequilibrate kinetically. Nevertheless, to the best of our knowledge, their contribution to bulkviscosity has not been studied in detail, even though it has been stressed that charm quarksmay affect a related quantity, the speed of sound, in a substantial way [27].In this paper, we do not address the bulk viscosity per se; rather, we consider the imaginary-time correlator from which it can be extracted non-perturbatively. Our goal is to computea mesonic part of this correlator (cf. eq. (3.1)) up to next-to-leading order (NLO) in pertur-bation theory. Given that thermal perturbation theory in general works better for mesoniccorrelators than for gluonic ones, and that the existence of a non-zero mass scale brings us fur-ther towards the asymptotically free regime, we should expect to find quantitatively accurateresults at T > ∼
300 MeV, as has previously been demonstrated in the vector channel [28].The plan of this paper is the following. After elaborating on the general physics of thebulk channel (sec. 2), the setup of the computation is outlined (sec. 3) and the main analyticresults are presented (sec. 4). Numerical illustrations and comparisons with non-perturbativedata comprise sec. 5, whereas sec. 6 collects together our findings. Two appendices containvarious details related to the NLO computation. In ref. [3] massive quarks were considered but results were only worked out for M ≪ πT . . Physics background Making use of dimensional regularization (the spacetime dimension is denoted by D = 4 − ǫ )and expressing the QCD action as S = Z d t Z d − ǫ x (cid:26) − F aµν F aµν + ¯ ψ (cid:16) i ←→ /D − M B (cid:17) ψ (cid:27) , (2.1)the trace of the energy-momentum tensor, ˆ θ µν , is [29, 30]ˆ θ µµ = c θ g B F aµν F aµν | {z } ≡ θ + ¯ ψM B ψ , (2.2)where equations of motion were used for the quark fields. Here c θ is a numerical factor ; g B is the bare gauge coupling; and g is a dimensionless renormalized gauge coupling, evaluatedin the MS scheme at the renormalization scale ¯ µ . The bare mass parameter M B is assumedto be a diagonal N f × N f matrix.The trace ˆ θ µµ has a non-zero thermal expectation value, (cid:10) ˆ θ µµ (cid:11) T = e − p = T dd T (cid:16) pT (cid:17) , (2.3)where e is the energy density, p is the pressure, T is the temperature, and chemical potentialsare assumed zero. In a scale-invariant theory, in which p ∝ T , this expectation valuevanishes. In QCD, in contrast, it is non-zero, both because of dimensional transmutationand because of non-zero mass parameters. The former effect originates from loop corrections,and indeed e − p = O ( c θ g T ) in massless QCD at T ≫
150 MeV. In contrast, in thepresence of masses, the expectation value is non-zero even in a free theory: (cid:10) ˆ θ µµ (cid:11) T = 4 N c N f X i =1 M i Z p n F ( E p,i ) E p,i + O ( g ) , Z p ≡ Z d p (2 π ) , E p,i ≡ q p + M i , (2.4)where n F is the Fermi distribution and M i is a renormalized quark mass of flavour i . Theexpectation value in eq. (2.4) vanishes in the chiral limit M i ≪ πT and, because of Boltzmannsuppression, also for large masses, M i ≫ πT . Yet it can give a rather substantial contributionfor M i ∼ πT ; for temperatures relevant for heavy ion collision experiments, this could be thecase with charm quarks [31]–[34], which have a mass M c < M D = 1 .
86 GeV.In the following, we consider 2-point correlation functions of ˆ θ µµ . Of particular interest isthe real-time correlator ( X ≡ ( t, x )) ζ = 19 lim ω → + (cid:26) ω Z X e iωt (cid:28) h ˆ θ µµ ( X ) , ˆ θ µµ (0) i(cid:29) T (cid:27) , (2.5) More precisely, c θ = lim ǫ → D − g = − b − b g + . . . , b = N c − T F π ) , b = N − N c T F − C F T F π ) . Theusual group theory factors are N c = 3, C F ≡ ( N − / (2 N c ) and T F ≡ N f / ζ ∼ T g (cid:16) − c s (cid:17) ∼ c θ g T , (2.6)where c s is the speed of sound. In a system with vanishing chemical potentials,13 − c s = T p ′′ dd T (cid:16) p ′ T (cid:17) . (2.7)Even though different from h ˆ θ µµ i T in eq. (2.3), eq. (2.7) also measures the breaking of scaleinvariance, and shows a similar parametric behaviour as h ˆ θ µµ i T . Therefore, we may expectthat ζ is significantly influenced by quark masses.The argument can be made more precise by relating ζ to the heavy quark chemical equili-bration rate. For M i ≫ πT , ˆ θ µµ of eq. (2.2) is dominated by the same term ¯ ψM i ψ that alsodominates the heavy quark Hamiltonian. The shape of the corresponding spectral functionwas discussed in ref. [35]; if we take the limit in eq. (2.5) all the way down to frequencies ω < ∼ Γ chem , where Γ chem is the heavy quark chemical equilibration rate, then the heavy quarkcontribution to the bulk viscosity may be estimated as δζ = 118 T lim ω → (cid:26) M i χ f Γ chem ω + Γ chem (cid:27) = M i χ f T Γ chem . (2.8)Here χ f denotes the heavy flavour susceptibility. Recalling that at weak coupling [10, 28]Γ chem = g C F πM i (cid:16) N f + 2 C F − N c (cid:17)(cid:16) T M i π (cid:17) e − M i /T (cid:18) O (cid:16) TM i , r g M i T (cid:17)(cid:19) , (2.9) χ f = 4 N c (cid:16) M i T π (cid:17) e − M i /T (cid:18) O (cid:16) TM i , g TM i (cid:17)(cid:19) , (2.10)it is seen that δζ ∼ M i /g T . Therefore δζ exceeds the gluonic contribution in eq. (2.6) byfour perturbative orders, O (1 /g ), as is the case also for 0 < M i ≪ πT [3].The purpose of the present study is to consider the imaginary-time correlator correspondingto eq. (2.5): (cid:28)Z x ˆ θ µµ ( X )ˆ θ µµ (0) (cid:29) T , (2.11)where X ≡ ( τ, x ), 0 < τ < /T , and heavy quarks are assumed to be in full equilibrium. Recalling eq. (2.2), this correlator contains three terms. The 2-point correlator of c θ θ hasbeen computed up to NLO in the weak-coupling expansion [36, 37] and compared with Given that in real-world heavy ion collisions heavy quarks appear in overabundance, this can in somesense be considered a conservative assumption. ω ∼ fm /c that can play a role in the hydrodynamical evolution of the system. If the width of the charmquark transport peak is narrower than this (Γ chem < ω ), then only gluons and light quarkscontribute to the bulk viscosity. If Γ chem ∼ ω , then heavy quarks should be included but to dothis consistently requires going beyond a hydrodynamical description, perhaps by employingkinetic theory. If Γ chem > ω , a hydrodynamical description applies and the heavy quarkcontribution should be added to the bulk viscosity. In each case the chemical equilibrationrate Γ chem is seen to be a fundamental quantity, whose non-perturbative determination as afunction of the heavy quark mass would be more than welcome.
3. Setup of the computation
For simplicity, we consider a situation in the following in which there is one quenched heavyquark, of bare mass M B and renormalized mass M (different schemes are specified presently),and N f massless dynamical quarks. (This means that from now on the quarks in sec. 2 shouldbe thought of as having N f + 1 flavours.) The imaginary-time correlator considered is G S ( τ ) ≡ M B Z x D ( ¯ ψψ )( τ, x ) ( ¯ ψψ )(0 , ) E T , < τ < β , β ≡ T , (3.1)where ψ denotes a single-flavour heavy quark Dirac spinor. Defined this way, the NLO scalardensity correlator is finite after mass and gauge coupling renormalization.We compute the correlator by first determining the corresponding correlator in momentumspace, with an external four-momentum Q ≡ ( ω n , ) , (3.2)where ω n is a bosonic Matsubara frequency. Denoting that result by ˜ G S ( ω n ), the correlatorof eq. (3.1) is obtained from G S ( τ ) = T X ω n e − iω n τ ˜ G S ( ω n ) . (3.3)4iven that the definition in eq. (3.1) involves a bare parameter, the issue of renormalizationand quark mass definitions plays an important role. The bare mass parameter, M B , can beexpressed as M B = M + δM , where the choice of δM defines a scheme. We write δM = − g C F M (4 π ) (cid:18) ǫ + ln ¯ µ M + 43 + δ (cid:19) + O ( g ) , (3.4)where ¯ µ is the scale parameter of the MS scheme, and terms of O ( ǫ ) were omitted. For δ = 0, M corresponds to a pole mass, which tends to compactify analytic expressions butis ambiguous on the non-perturbative level and also leads to problems of convergence (seebelow). If we choose δ = − ln ¯ µ M − , (3.5)then M stands for the MS mass, which we denote by m (¯ µ ); its asymptotic running reads m (¯ µ ) = m (¯ µ ref ) (cid:20) ln(¯ µ ref / Λ MS )ln(¯ µ/ Λ MS ) (cid:21) C F11 N c − T F , (3.6)where typically ¯ µ ref ≡ δ open for the moment.
4. Analytic results
Denoting propagators by ∆ P ≡ P + M (4.1)and Matsubara sum-integrals by Σ R { P } ≡ T P { p n } R p , where { P } stands for fermionic Mat-subara momenta, the tree-level correlator reads= − N c M PZ { P } (cid:26) − P + Q + 4 M ∆ P ∆ P − Q (cid:27) . (4.2)The counterterm contribution is= − N c δM PZ { P } (cid:26) − P + 2 M ∆ P + Q + 8 M ∆ P ∆ P − Q − M ( Q + 4 M )∆ P ∆ P − Q (cid:27) , (4.3)5hereas the “genuine” 2-loop graphs amount to+ = 4 g N c C F M PZ K { P } (cid:26) − D − K ∆ P + D − P ∆ P − K + 2 K ∆ P ∆ P − K − M K ∆ P ∆ P − K + ( D − Q + 4 M ) K ∆ P ∆ P − Q − ( D − Q + 4 M )∆ P ∆ P − K ∆ P − Q − M − D − K · Q + 4 Q K ∆ P ∆ P − K ∆ P − Q + 4 M ( Q + 4 M ) K ∆ P ∆ P − K ∆ P − Q − D − M + ( D − Q ∆ P ∆ P − K ∆ P − Q ∆ P − K − Q + 8 M + 6 M Q + Q K ∆ P ∆ P − K ∆ P − Q ∆ P − K − Q (cid:27) . (4.4) All the Matsubara sums (cf. appendix A) as well as some of the angular integrals appearingin eqs. (4.2)–(4.4) can be carried out analytically; in addition partial integrations permit forsimplifications (cf. ref. [40] for the massless case). To display the results, we employ thefunctions D E k +1 ··· E n E ··· E k ( τ ) ≡ e ( E + ··· + E k )( β − τ )+( E k +1 + ··· + E n ) τ + e ( E + ··· + E k ) τ +( E k +1 + ··· + E n )( β − τ ) [ e βE ± · · · [ e βE n ± , (4.5)where the sign in the denominator is chosen according to whether the particle is a boson ora fermion. The energy variables ǫ k ≡ | k | , E p ≡ p p + M , E pk ≡ p ( p − k ) + M , E ± pk ≡ p ( p ± k ) + M (4.6)appear frequently, and we denote D E p ≡ D E p E p . Then the leading-order (LO) result reads G LOS ( τ ) | τ -dep. = 2 N c M Z p p D E p ( τ ) E p , (4.7) G LOS ( τ ) | const. = − N c M Z p M T n ′ F ( E p ) E p , (4.8)where a τ -independent part stemming from an approximate transport peak in the corre-sponding spectral function has been separated. Scheme dependence can be expressed as∆ δ G NLOS ( τ )4 g N c C F M (cid:12)(cid:12)(cid:12)(cid:12) τ -dep. = 3 δ π Z p D E p ( τ ) (cid:18) M E p − (cid:19) , (4.9)∆ δ G NLOS ( τ )4 g N c C F M (cid:12)(cid:12)(cid:12)(cid:12) const. = − δ π Z p M T n ′ F ( E p ) p (cid:18) M E p − (cid:19) . (4.10)6ntroducing the shorthand notations∆ στ ≡ ǫ k + σE p + τ E pk , ∆ σ = E p + σE pk , (4.11)the τ -dependent part of the NLO correction reads (after renormalization, cf. appendix B) G NLOS | τ -dep. g N c C F M = Z p D E p ( τ )8 π (cid:20) M E p (cid:18) pE p ln E p + pE p − p (cid:19) − − (cid:18) − M E p (cid:19) Z ∞ d k θ ( k ) k (cid:21) (4.12)+ Z p,k P (Z z D ǫ k E p E pk ( τ ) M ǫ k E p E pk ∆ + − ∆ − + (cid:20) ǫ k + ( E p + E pk ) M + ∆ −− ∆ ++ − ǫ k ∆ + − ∆ − + + 4 ǫ k M ∆ ∆ + − ∆ − + (cid:21) + Z z D ǫ k E p E pk ( τ ) M ǫ k E p E pk ∆ + − ∆ − + (cid:20) ǫ k + ( E p + E pk ) M + ∆ ++ ∆ −− − ǫ k ∆ + − ∆ − + + 4 ǫ k M ∆ −− ∆ + − ∆ − + (cid:21) − Z z D E p ǫ k E pk ( τ ) M ǫ k E p E pk ∆ ++ ∆ −− (cid:20) ǫ k + ( E p − E pk ) M + ∆ + − ∆ − + − ǫ k ∆ ++ ∆ −− + 4 ǫ k M ∆ − + ∆ ++ ∆ −− (cid:21) + D E p ( τ )2 ǫ k (cid:18) − M E p (cid:19) (cid:20) − E p ( E + pk − E − pk ) − pǫ k ( E + pk + E − pk )2 p ( E p − ǫ k ) − ǫ k M ( E + pk − E − pk ) p ( E p − ǫ k ) E + pk E − pk − E p − M pE p (cid:18) ln (cid:12)(cid:12)(cid:12)(cid:12) ( E p + p )(2 p − ǫ k )( E p − p )(2 p + ǫ k ) (cid:12)(cid:12)(cid:12)(cid:12) + ln (cid:12)(cid:12)(cid:12)(cid:12) − ǫ k / ( E p + E − pk ) − ǫ k / ( E p + E + pk ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) + θ ( k ) (cid:21) + D E p ( τ ) n B ( ǫ k ) ǫ k (cid:20) − pE p ln E p + pE p − p + (cid:18) − M E p (cid:19)(cid:18) ǫ k − p − E p − M pE p ǫ k ln E p + pE p − p (cid:19) (cid:21) + D E p ( τ ) n F ( E k ) E k (cid:20) − E p − E p + M p − k ) E p + M pkE p ln (cid:12)(cid:12)(cid:12)(cid:12) p + kp − k (cid:12)(cid:12)(cid:12)(cid:12) + E p ( E p + E k − M ) + M pk ( E p − E k ) E p (cid:18) ln (cid:12)(cid:12)(cid:12)(cid:12) p + kp − k (cid:12)(cid:12)(cid:12)(cid:12) + E k E p + E k ln (cid:12)(cid:12)(cid:12)(cid:12) M + E p E k + pkM + E p E k − pk (cid:12)(cid:12)(cid:12)(cid:12)(cid:19)(cid:21)) , where P denotes a principal value, R z an integral over the angles between p and k (with R z n B the Bose distribution, and the function θ ( k ) is specified in eq. (4.16). Theconstant contribution reads G NLOS | const. g N c C F M = Z p T n ′ F ( E p ) M E p (cid:26) π + Z k n B ( ǫ k ) ǫ k (cid:20) p (cid:21) + Z k n F ( E k ) E k (cid:20) p (cid:18) − M k (cid:19) + 1 k + 1 E k + 4 E k − M pkE k ln (cid:12)(cid:12)(cid:12)(cid:12) p + kp − k (cid:12)(cid:12)(cid:12)(cid:12)(cid:21)(cid:27) . (4.13)Note that for M ≫ πT the leading thermal correction in eq. (4.13), from R k n B ( ǫ k ) /ǫ k , agreeswith the leading NLO thermal correction to χ f [28], as is to be expected from integrationover the transport peak in eq. (2.8) with the standard relation in eq. (5.16).7 .3. Infrared and ultraviolet regimes For the numerical evaluation of eq. (4.12) it must be kept in mind that individual parts of theexpression contain divergences; only the sum is well-defined. In particular, at small k there isa divergence originating from the terms integrated over z in eq. (4.12) and from terms wherethe integral had already been carried out. For the latter type the small- k part reads D E p ( τ ) ǫ k (cid:18) − M E p (cid:19) (cid:18) − E p − M pE p ln E p + pE p − p (cid:19)(cid:20)
12 + n B ( ǫ k ) (cid:21) , (4.14)containing both vacuum and Bose-enhanced structures, leading to logarithmic and powerlikedivergences. When terms of both origins are summed together, the small- k divergences cancel.At large k the leading asymptotic behaviour is D E p ( τ ) ǫ k (cid:18) − M E p (cid:19) (cid:18) E p ǫ k (cid:19) . (4.15)Although integrable this expression is only power-suppressed, and for numerical handling itmay be advantageous to accelerate the decrease. The leading tail can be subtracted withthe help of the auxiliary function θ ( k ) appearing on the 1st and 6th rows of eq. (4.12), forinstance θ ( k ) ≡ − E p Θ( k − k min ) k + λ , Z ∞ d k θ ( k ) k = − E p λ ln (cid:16) λ k min (cid:17) . (4.16)The resolution obtained can be increased with a suitable tuning of λ and k min ; as typicalvalues we have used λ ∼ M/ k min ∼ T .
5. Numerical evaluations
For a comparison with lattice simulations certain parameter values need to be fixed. Boththe gauge coupling and the masses are running parameters; in accordance with ref. [39] weset Λ MS ≃
216 MeV to fix the gauge coupling (note that the simulations in ref. [39] arefor N f = 0). As a rule the gauge coupling is evolved with 3-loop running (cf. ref. [42]and references therein), with some exceptions as specified below. For the running massesa standard choice according to ref. [41] is to evaluate them at ¯ µ ref = 2 GeV, and then thecharm mass is m c (¯ µ ref ) = 1 . m c ( m c ) = 1 . m c (¯ µ ref ) = 0 . µ ,which can be used for estimating the uncertainties of a fixed-order computation. We always8valuate the result at some “optimal” scale ¯ µ opt as well as at 0 . µ opt and 2 . µ opt ; the maximaand minima among these three numbers are then used for constructing “error bands” for theperturbative predictions. Obviously such bands only serve as lower bounds for the systematicuncertainties related to the perturbative computation.In order to display the results in a useful way, they should be normalized to an appropriateexpression. We have considered three options for this purpose: a free result with a light quarkmass ( M ≪ πT ); the scalar correlator in pure Yang-Mills theory; as well as a “reconstructed”correlator following from a zero-temperature spectral function. The LO result for the scalar channel correlator is given in eqs. (4.7), (4.8). It vanishes in thechiral limit, and cannot be evaluated in a closed form for a general mass. Nevertheless, if weassume that M ≪ πT , we can set the mass to zero within the integrand, yielding [44] G naive S ( τ ) ≡ N c T M (cid:20) π (1 − τ T ) 1 + cos (2 πτ T )sin (2 πτ T ) + 2 cos(2 πτ T )sin (2 πτ T ) (cid:21) . (5.1)It turns out, however, that this expression does not compare well with lattice data even for τ ≪ /T ; the situation is illustrated in fig. 1. The reason can be understood as being relatedto running effects, as we now explain.Let us start by considering eq. (5.1) at short distances, G naive S ( τ ) τ ≪ T ≈ N c M π τ . (5.2)Then compute the NLO correction to this result; making use of the spectral function ineq. (5.13), the LO + NLO expression reads G LO+NLOS ( τ ) τ ≪ M , T ≈ N c M π τ (cid:20) g C F (4 π ) (cid:18) ln τ M − δ −
32 + 2 γ E (cid:19)(cid:21) . (5.3)We note that if M is kept fixed, G LO+NLOS turns negative at small τ , but that simultaneouslythe loop expansion breaks down because the NLO correction overtakes the LO term. Theproblem can be rectified if we use running parameters: going to the MS scheme; choosing δ according to eq. (3.5) and subsequently replacing M → m (¯ µ ) according to eq. (3.6); andtaking ¯ µ to scale with τ according to a “fastest apparent convergence” criterion,¯ µ → βe − γ E τ ( β − τ ) , (5.4)9 .0 0.1 0.2 0.3 0.4 0.5 τ T G Slattice [N τ = 48] / G Sfree G Slattice [N τ = 48] / G Snaive
T = 1.45 T c , T c = 1.25 Λ MS _ Figure 1:
The ratio G lattice S /G free S (with m (¯ µ ref ) = 967 MeV; open symbols) as well as G lattice S /G naive S (with M/T = 3 . G free S (cf. eq. (5.5)) yields results closer to unity atshort distances; the very smallest distances are affected by lattice artifacts. such that the NLO correction in eq. (5.3) always remains small, we are led to define G free S ( τ ) ≡ N c T m (¯ µ ref ) ( ln (cid:2) ¯ µ refΛ MS (cid:3) ln (cid:2) βe − γ E τ ( β − τ )Λ MS (cid:3) ) C F11 N c − T F × (cid:20) π (1 − τ T ) 1 + cos (2 πτ T )sin (2 πτ T ) + 2 cos(2 πτ T )sin (2 πτ T ) (cid:21) . (5.5)Let us stress that this is a definition rather than an exact result; indeed only the asymptoticsat τ ≪ β can be fixed unambiguously thanks to asymptotic freedom. Nevertheless, as seenin fig. 1, eq. (5.5) agrees much better with lattice data at small distances than eq. (5.1).The full NLO results normalized to eq. (5.5) are shown in fig. 2 as a function of m (¯ µ ref ).As illustrated with the example of m (¯ µ ref ) = 0 .
25 GeV, the renormalization scale dependencegets significantly reduced when going from the LO to the NLO level. The NLO results for m (¯ µ ref ) = 1 GeV agree with the lattice data at small τ , but at large τ a discrepancy sets in.We return to the discrepancy in connection with fig. 3.10 .0 0.1 0.2 0.3 0.4 0.5 τ T G S / G S fr ee LO [0.25 GeV]NLOlattice [N τ = 48]T = 1.45 T c , T c = 1.25 Λ MS _ Figure 2:
The ratio G S /G free S (cf. eq. (5.5)) as a function of τ T and the quark mass m (¯ µ ref ), indicatednext to the coloured bands. The bands reflect the uncertainty related to the choice of the renormal-ization scale as specified in the text. The lattice simulations are from ref. [39] and correspond to m (¯ µ ref ) ≃
967 MeV. The data at τ T > .
15 are suppressed with respect to the continuum prediction.
The spectral function related to the gluonic part of the trace anomaly, i.e. the operator c θ θ defined in eq. (2.2), was computed up to NLO in ref. [37]. At LO it reads (for ω > ∼ πT ) ρ LO θ ( ω ) = N c C F c θ g ω π (cid:20) n B (cid:16) ω (cid:17)(cid:21) , (5.6)which yields the correlator G naive θ ( τ ) c θ g N c C F ≡ π T (cid:20) π (1 − τ T ) 2 cos(2 πτ T ) + cos (2 πτ T )sin (2 πτ T ) + 1 + 2 cos (2 πτ T )sin (2 πτ T ) (cid:21) . (5.7)However, there are again loop corrections which imply that the running of the coupling needsto be taken into account; the asymptotics reads [45] ρ vac θ ( ω ) ω ≫ πT ≈ N c C F c θ ω π (cid:26) g (¯ µ ) + g (¯ µ ) N c (4 π ) (cid:18)
223 ln ¯ µ ω + 733 (cid:19) + O ( g ) (cid:27) , (5.8)with thermal corrections strongly suppressed at ω ≫ πT [46, 36]. Inserting into eq. (5.16) itis seen that at τ ≪ β the effects from the running can be captured through a choice analogous11 .0 0.1 0.2 0.3 0.4 0.5 τ T G S / G θ fr ee LO [0.25 GeV]NLOlattice [N τ = 48]T = 1.45 T c , T c = 1.25 Λ MS _ Figure 3:
The ratio G S /G free θ (cf. eq. (5.10)) as a function of τ T and the quark mass m (¯ µ ref ),indicated next to the coloured bands. The bands reflect the uncertainty related to the choice of therenormalization scale. Recalling m c (¯ µ ref ) = 1 . − m (¯ µ ref ) ≃
967 MeV are suppressed with respect to the perturbative prediction. to eq. (5.4), viz. ¯ µ → βe − γ E τ ( β − τ ) , (5.9)and this leads us to define G free θ ( τ ) ≡ c θ N c C F T ( π (11 N c − T F ) ln (cid:2) βe − γ E τ ( β − τ )Λ MS (cid:3) ) × (cid:20) π (1 − τ T ) 2 cos(2 πτ T ) + cos (2 πτ T )sin (2 πτ T ) + 1 + 2 cos (2 πτ T )sin (2 πτ T ) (cid:21) . (5.10)Only the asymptotics at τ ≪ β is unambiguously fixed, otherwise eq. (5.10) represents achoice. (To be concrete, for this definition we have evaluated the coupling appearing in c θ ,cf. footnote 2, at 1-loop level at the scale indicated by eq. (5.9).)The NLO result for G S , normalized to G free θ , is shown in fig. 3. With this normalizationthe agreement at small τ and the discrepancy at large τ of the NLO expression and latticedata becomes clearly visible. Independently of lattice data, it is also seen that the relativemagnitude of the scalar density correlator in the infrared domain τ ∼ β/ (¯ µ ref ) ≈ . − .
25 GeV, quite close to the experimental value m c (¯ µ ref ) = 1 . ζ remains to be seen but we nevertheless considerthe pattern seen in fig. 3 to be intriguing. On the lattice side it has become fashionable to normalize thermal imaginary-time correla-tors to a “reconstructed” correlator. Even though we suspect that this increases systematicuncertainties (cf. below), we have worked out this case as well for completeness.Within perturbation theory, the massive scalar channel spectral function at zero tempera-ture reads (the structure is similar to the vector channel, cf. ref. [47]) ρ vac S ( ω ) = θ ( ω − M ) N c M ( ω − M ) πω + θ ( ω − M ) 4 g N c C F M (4 π ) ω (cid:26) (5.11)( ω − M )( ω − M ) L (cid:18) ω − √ ω − M ω + √ ω − M (cid:19) + (cid:18) ω − ω M − M (cid:19) acosh (cid:18) ω M (cid:19) − ω ( ω − M ) (cid:20) ( ω − M ) (cid:18) ln ω ( ω − M ) M + 34 δ (cid:19) −
38 (3 ω − M ) (cid:21)(cid:27) + O ( g ) , where the function L is defined as L ( x ) ≡ ( x ) + 2 Li ( − x ) + [2 ln(1 − x ) + ln(1 + x )] ln x . (5.12)The ultraviolet asymptotics from here is ρ vac S ( ω ) ω ≫ M ≈ N c ω M π (cid:26) g C F (4 π ) (cid:18) ln M ω − δ + 32 (cid:19) + O ( g ) (cid:27) . (5.13)We note again that in the pole mass scheme, where M is constant and δ = 0, eq. (5.13)turns negative at very large ω and the NLO correction overtakes the LO term, implyinga breakdown of the perturbative series. If, however, we go to the MS scheme, choosing δ according to eq. (3.5) and setting M → m (¯ µ ) according to eq. (3.6), we obtain anotherrepresentation: ρ vac S ( ω ) ω ≫ M ≈ N c ω m (¯ µ ref )8 π (cid:20) ln(¯ µ ref / Λ MS )ln(¯ µ/ Λ MS ) (cid:21) C F11 N c − T F × (cid:26) g (¯ µ ) C F (4 π ) (cid:18) ln ¯ µ ω + 176 (cid:19) + O ( g ) (cid:27) . (5.14)Taking ¯ µ to scale with ω , a series is obtained which becomes increasingly convergent as ω grows. The result is independent of the precise choice of ¯ µ up to higher-order corrections; inpractice we set ¯ µ → max( πT, ω e − / ) . (5.15) The infrared cutoff πT only affects the smallest masses, m (¯ µ ref ) < ∼ πT . .0 0.1 0.2 0.3 0.4 0.5 τ T G S / G S r ec LO [0.25 GeV]NLOlattice [N τ = 48]T = 1.45 T c , T c = 1.25 Λ MS _ Figure 4:
The ratio G S /G rec S (cf. eq. (5.16)) as a function of τ T and the quark mass m (¯ µ ref ),indicated next to the coloured bands. The bands reflect the uncertainty related to the choice of therenormalization scale. The lattice simulations are from ref. [39]. The comparison contains uncontrolleduncertainties from both sides as discussed in the text. Subsequently the “reconstructed” correlator is defined by taking the vacuum spectral functionand using it to compute a thermal Euclidean correlator from G rec S ( τ ) ≡ Z ∞ d ωπ ρ vac S ( ω ) cosh (cid:16) β − τ (cid:17) ω sinh βω . (5.16)In fig. 4 results normalized to the reconstructed correlator are shown. As expected, theresults approach unity at short distances for all masses (the small- τ result is determined bylarge ω and then thermal corrections are strongly suppressed with respect to the vacuumterm [46, 48]). However the match to lattice data is not as good as in figs. 2, 3; amongpossible explanations we can envisage the following: • On the continuum side it is questionable whether perturbation theory, which misses allresonance contributions, can reflect the vacuum scalar spectral function even qualita-tively in the ω -range dominating eq. (5.16) at moderate τ , ω > ∼ πT ∼ • On the lattice side the theoretical formulae (e.g. eq. (5.16)) used for relating measureddata to the reconstructed correlator assume a continuous τ -variable; this necessitatestaking a continuum limit which has not been reached.14 The “reconstructed” correlator of ref. [39] was not measured at T = 0 but at T = 0 . T c .Because of these issues we refrain from speculating on the origins of the discrepancy in fig. 4.
6. Conclusions
The purpose of this paper has been to investigate the influence of a finite charm quark masson a 2-point imaginary-time correlator in the so-called bulk (or scalar) channel, correspondingto the trace of the energy-momentum tensor in continuum QCD. We could have anticipated arather substantial influence, given that in massless QCD the corresponding correlator is non-zero only because of an “anomalous” breaking of scale invariance, and is therefore proportionalto a high power of running coupling constants (cf. eq. (5.6)).Remarkably, the relative influence of a massive quark on the scalar channel correlator atlarge imaginary-time separations appears to peak for m (¯ µ ref ) = 1 . − .
25 GeV (cf. fig. 3), notfar from the physical charm mass m c (¯ µ ref ) = 1 . ζ for all quark masses as wellas refined lattice investigations of scalar correlators in the future, including a systematicapproach to the continuum limit.An interesting principal application of the scalar correlator is that it can be used fordetermining the charm quark chemical equilibration rate close to equilibrium [35], a quantitythat experiences a non-trivial (Sommerfeld) enhancement even deep in the non-relativisticregime [10]. The determination poses however a numerical challenge because physics resides ina very narrow transport peak which is difficult to resolve from Euclidean data. (For M ≫ πT the peak is exponentially narrow, cf. eqs. (2.8), (2.9), rather than only 1 /M -suppressed likein the vector channel [49]; however for m (¯ µ ref ) ∼ T > ∼
400 MeV thedifference might not be dramatic. Note also that there are significant challenges in resolvingthe transport peak even in pure Yang-Mills theory [4], however with connected mesoniccorrelators a higher numerical precision can typically be reached.)Taking the current lattice data at face value, we find discrepancies much larger than in thevector channel where an identical data set was used (cf. figs. 2, 3 vs. ref. [28]). Part of thereason is surely that the quark mass, with associated uncertainties in its determination, nowplays a more prominent role. In order to get the uncertainties under control a continuumlimit needs to be taken. Nevertheless it is also tempting to speculate that the scalar channelmight experience larger non-perturbative effects than the vector one. Of course, these maypartly originate from quarkonium physics rather than the transport regime, cf. e.g. ref. [39]and references therein.For the full trace of eq. (2.2), an open challenge both on the lattice and in continuum isthe study of the “mixed” channel, with correlators involving crossterms between mesonic and15urely gluonic operators. On the lattice this correlator might be rather noisy numerically,but hopefully suitable methods will eventually be developed. A similar problem concerns thedisconnected contraction of mesonic scalar densities.It may finally be wondered why we have concentrated on the bulk rather than the moreprominent shear channel. The reason is that for bulk viscosity the effects from massless QCDare “anomalously” suppressed, so that the contribution from the explicit breaking of scaleinvariance through quark masses is relatively speaking more important. In the shear channelthe charm mass also has an influence, however we do not expect it to be larger than ∼ T > ∼
400 MeV before [31]–[34].
Note added
We have recently explored possible reasons for the discrepancy between NLO expressions andlattice data in the scalar channel, by varying the corresponding spectral function, and foundthat the discrepancy is likely to originate from quarkonium-related physics, more preciselyfrom the proper location of the quark-antiquark threshold [50].
Acknowledgements
We thank H.-T. Ding for helpful discussions and providing us with lattice data from ref. [39],and G.D. Moore for helpful discussions. This work was partly supported by the Swiss NationalScience Foundation (SNF) under the grant 200021-140234 as well as under the Ambizionegrant PZ00P2-142524.
Appendix A. Master sum-integrals
The way that the computation is organized is by first using substitutions of sum-integrationvariables in order to express the result in terms of a small number of “master” sum-integrals,or basis functions, cf. eq. (4.4). After transformation to coordinate space the basis functionscan be defined as I m m m n n n n n ( τ ) ≡ T X ω n e − iω n τ PZ K { P } ( M ) m ( Q ) m (2 K · Q ) m ( K ) n ∆ n P ∆ n P − K ∆ n P − Q ∆ n P − K − Q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Q =( ω n , ) . (A.1)In ref. [28] expressions obtained after carrying out the Matsubara sums were given for all thebasis functions appearing in the present computation, and these have been employed in orderto arrive at eqs. (B.1), (B.2) below. Here we just remark, for completeness, that one of the16asis functions discussed in ref. [28] is actually not independent of the others: I − ( τ ) = 2 I ( τ ) − I ( τ ) − I ( τ ) . (A.2) Appendix B. Renormalization of the scalar channel correlator
Expressing eq. (4.4) in terms of the basis functions listed in appendix A of ref. [28], we obtainresults for the scalar correlator in which the Matsubara sums have been carried out: G NLOS | τ -dep. g N c C F M = Z p , k D ǫ k E p E pk ( τ ) M ǫ k E p E pk ∆ + − ∆ − + (cid:20) ǫ k + ( E p + E pk ) M + ∆ −− ∆ ++ − ǫ k ∆ + − ∆ − + + 4 ǫ k M ∆ ∆ + − ∆ − + (cid:21) + Z p , k D ǫ k E p E pk ( τ ) M ǫ k E p E pk ∆ + − ∆ − + (cid:20) ǫ k + ( E p + E pk ) M + ∆ ++ ∆ −− − ǫ k ∆ + − ∆ − + + 4 ǫ k M ∆ −− ∆ + − ∆ − + (cid:21) − Z p , k D E p ǫ k E pk ( τ ) M ǫ k E p E pk ∆ ++ ∆ −− (cid:20) ǫ k + ( E p − E pk ) M + ∆ + − ∆ − + − ǫ k ∆ ++ ∆ −− + 4 ǫ k M ∆ − + ∆ ++ ∆ −− (cid:21) + Z p D E p ( τ ) (cid:26) “eq. (B.3)”+ Z k n B ( ǫ k ) ǫ k (cid:20) M E p − M ( E p − M )2 E p E pk (cid:18) + 1∆ −− − − − − + (cid:19) + E p + E pk − M E p E pk (cid:18) ++ ∆ −− − + − ∆ − + (cid:19) + M (2 E p − M ) E p (cid:18) ++ ∆ −− + 1∆ + − ∆ − + (cid:19)(cid:21) + Z k n F ( E pk ) E pk P (cid:20) M E p − M ( E p − M )2 ǫ k E p (cid:18) − + + 1∆ −− − − − (cid:19) + E p + E pk − M E p (cid:18) + ∆ ++ ∆ −− + 1∆ − ∆ + − ∆ − + (cid:19) − M E p ∆ + ∆ − + E pk M (2 E p − M ) E p (cid:18) + ∆ ++ ∆ −− − − ∆ + − ∆ − + (cid:19)(cid:21)(cid:27) + Z p E p ∂ E p D E p ( τ ) (cid:26) Z k n B ( ǫ k ) ǫ k (cid:20) E p − M E p − M ( E p − M )2 E p (cid:18) ++ ∆ + − + 1∆ −− ∆ − + (cid:19)(cid:21) + Z k n F ( E pk ) E pk (cid:20) E p − M E p − M ( E p − M )2 E p (cid:18) ++ ∆ −− + 1∆ + − ∆ − + (cid:19)(cid:21)(cid:27) . (B.1)17he constant contribution, in turn, can be expressed as G NLOS | const. g N c C F M = Z p , k T n ′ F ( E p ) n ′ F ( E pk )2 E p E pk (cid:26) − M + 2 M (cid:18) ++ ∆ −− + 1∆ + − ∆ − + (cid:19)(cid:27) + Z p T n ′ F ( E p ) (cid:26) “eq. (B.4)”+ Z k n B ( ǫ k ) ǫ k (cid:20) M E p − ǫ k M E p (cid:18) ∆ − − −− ∆ − + (cid:19) + M (2 E p − M ) E p (cid:18) ++ ∆ + − + 1∆ −− ∆ − + (cid:19)(cid:21) + Z k n F ( E pk ) E pk (cid:20) M E p + M E p E pk + 2 M E p E pk (cid:18) ∆ ∆ ∆ −− − ∆ − ∆ − ∆ − + (cid:19) + (cid:18) M (2 E p − M ) E p − M E p E pk (cid:19)(cid:18) ++ ∆ −− + 1∆ + − ∆ − + (cid:19)(cid:21)(cid:27) + Z p T E p n ′′ F ( E p ) (cid:26) Z k n B ( ǫ k ) ǫ k (cid:20) − M E p + M E p (cid:18) ++ ∆ + − + 1∆ −− ∆ − + (cid:19)(cid:21) + Z k n F ( E pk ) E pk (cid:20) − M E p + M E p (cid:18) ++ ∆ −− + 1∆ + − ∆ − + (cid:19)(cid:21)(cid:27) . (B.2)The “0”s in eqs. (B.1), (B.2) represent vacuum contributions that vanish after renormal-ization. The coefficients of D E p ( τ ) and T n ′ F ( E p ) are also related to renormalization, but donot vanish:“eq. (B.3)” = Z k P (cid:26) (1 − ǫ )(2 M − E p )2 E p M (cid:18) ǫ k − E pk (cid:19) − ǫE p + (1 − ǫ ) M E p E pk ( E pk − E p )+ ǫ ǫ k + E pk ǫ k E pk [( ǫ k + E pk ) − E p ] + 2 M ( ǫ k + E pk )( E p − M ) ǫ k E pk E p [( ǫ k + E pk ) − E p ] + ( ǫ k + 2 E pk )( E p − M )(2 E p − M ) ǫ k E pk E p [( ǫ k + E pk ) − E p ]( E pk − E p ) (cid:27) , (B.3)“eq. (B.4)” = Z k P (cid:26) − ǫE p (cid:18) ǫ k − E pk (cid:19) − (1 − ǫ ) M E p E pk + [ − E p + ( ǫ k + E pk )( ǫ k + 3 E pk )] M E p E pk [( ǫ k + E pk ) − E p ] (cid:27) . (B.4)Reducing to a basis of independent structures as indicated in eqs. (B.17)–(B.20) of ref. [28],18nd carrying out the integrals in the structures with one or two propagators, we obtain“eq. (B.3)” = (cid:26)Z K P (cid:20) (cid:18) − M E p (cid:19)(cid:18) − ǫM ∆ K − − ǫM K + 1∆ K − ǫ ∆ K ∆ K − Q (cid:19) + 2 M E p (cid:18) K ∆ P − K − K ∆ K − Q (cid:19) + (cid:18) − M E p (cid:19) M K ∆ P − K + (cid:18) − M E p (cid:19) E p − M ) K ∆ P − K ∆ P − K − Q (cid:21)(cid:27) p = iE p ,Q =(2 iE p , ) + O ( ǫ )= 18 π (cid:20) M E p (cid:18) pE p ln E p + pE p − p (cid:19) − (cid:21) + (cid:18) − M E p (cid:19) Z k P (cid:26) − M E pk (cid:20) ǫ k ( ǫ k + E p )∆ + 12 ǫ k ( ǫ k − E p )∆ − + − ǫ k − E p ) E pk (cid:21) − E p − M ǫ k E p E pk (cid:18) + + 1∆ − − ++ + 1∆ − + (cid:19)(cid:27) + O ( ǫ ) , (B.5)“eq. (B.4)” = (cid:26)Z K P (cid:20) − ǫ ) E p (cid:18) K − K (cid:19) + 2 M E p K ∆ P − K − − ǫ ) M E p ∆ K (cid:21)(cid:27) p = iE p ,Q =(2 iE p , ) + O ( ǫ )= 38 π M E p + O ( ǫ ) . (B.6)In the k -integral of eq. (B.5) the angular integration can be carried out, which leads to thevacuum part of eq. (4.12); eq. (B.6) in turn yields the vacuum part of eq. (4.13). Note thatboth eq. (B.5) and (B.6) contain a counterterm contribution from eq. (4.3) and are thereforefinite in the ultraviolet regime (large k ). References [1] S. Weinberg,
Entropy generation and the survival of protogalaxies in an expanding uni-verse,
Astrophys. J. 168 (1971) 175.[2] S. Weinberg,
Gravitation and Cosmology (Wiley, New York, 1972).[3] P.B. Arnold, C. Dogan and G.D. Moore,
The bulk viscosity of high-temperature QCD,
Phys. Rev. D 74 (2006) 085021 [hep-ph/0608012].[4] G.D. Moore and O. Saremi,
Bulk viscosity and spectral functions in QCD,
JHEP 09(2008) 015 [0805.4201]. 195] E. Lu and G.D. Moore,
The Bulk Viscosity of a Pion Gas,
Phys. Rev. C 83 (2011)044901 [1102.0017].[6] G. Torrieri, B. Tom´aˇsik and I. Mishustin,
Bulk viscosity driven clusterization of quark-gluon plasma and early freeze-out in relativistic heavy-ion collisions,
Phys. Rev. C 77(2008) 034903 [0707.4405].[7] A. Monnai and T. Hirano,
Effects of Bulk Viscosity at Freezeout,
Phys. Rev. C 80 (2009)054906 [0903.4436].[8] K. Rajagopal and N. Tripuraneni,
Bulk Viscosity and Cavitation in Boost-InvariantHydrodynamic Expansion,
JHEP 03 (2010) 018 [0908.1785].[9] H. Song and U.W. Heinz,
Interplay of shear and bulk viscosity in generating flow inheavy-ion collisions,
Phys. Rev. C 81 (2010) 024905 [0909.1549].[10] D. B¨odeker and M. Laine,
Sommerfeld effect in heavy quark chemical equilibration,
JHEP01 (2013) 037 [1210.6153].[11] A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel,
Statistical hadronizationof heavy quarks in ultra-relativistic nucleus-nucleus collisions,
Nucl. Phys. A 789 (2007)334 [nucl-th/0611023].[12] M. Mustafa [STAR Collaboration],
Measurements of Non-photonic Electron Productionand Azimuthal Anisotropy in √ s NN = 39, 62.4 and 200 GeV Au+Au Collisions fromSTAR at RHIC, Nucl. Phys. A904-905 (2013) 665c [1210.5199].[13] G. Ortona [ALICE Collaboration],
Open-Charm Meson Elliptic Flow Measurement inPb-Pb Collisions at √ s NN = 2 . TeV with ALICE at the LHC,
Heavy quark diffusion in perturbative QCD at next-to-leading order,
Phys. Rev. Lett. 100 (2008) 052301 [0708.4232].[15] S. Caron-Huot, M. Laine and G.D. Moore,
A Way to estimate the heavy quark thermal-ization rate from the lattice,
JHEP 04 (2009) 053 [0901.1195].[16] H.B. Meyer,
The errant life of a heavy quark in the quark-gluon plasma,
New J. Phys.13 (2011) 035008 [1012.0234].[17] A. Francis, O. Kaczmarek, M. Laine and J. Langelage,
Towards a non-perturbative mea-surement of the heavy quark momentum diffusion coefficient,
PoS LATTICE 2011 (2011)202 [1109.3941]. 2018] D. Banerjee, S. Datta, R. Gavai and P. Majumdar,
Heavy Quark Momentum DiffusionCoefficient from Lattice QCD,
Phys. Rev. D 85 (2012) 014510 [1109.5738].[19] G.D. Moore and D. Teaney,
How much do heavy quarks thermalize in a heavy ion colli-sion?,
Phys. Rev. C 71 (2005) 064904 [hep-ph/0412346].[20] H. van Hees, M. Mannarelli, V. Greco and R. Rapp,
Nonperturbative heavy-quark diffu-sion in the quark-gluon plasma,
Phys. Rev. Lett. 100 (2008) 192301 [0709.2884].[21] M. Laine, G.D. Moore, O. Philipsen and M. Tassler,
Heavy Quark Thermalization inClassical Lattice Gauge Theory: Lessons for Strongly-Coupled QCD,
JHEP 05 (2009)014 [0902.2856].[22] F. Riek and R. Rapp,
Quarkonia and Heavy-Quark Relaxation Times in the Quark-GluonPlasma,
Phys. Rev. C 82 (2010) 035201 [1005.0769].[23] W.M. Alberico, A. Beraudo, A. De Pace, A. Molinari, M. Monteno, M. Nardi, F. Prinoand M. Sitta,
Heavy flavours in AA collisions: production, transport and final spectra,
Eur. Phys. J. C 73 (2013) 2481 [1305.7421].[24] C.P. Herzog, A. Karch, P. Kovtun, C. Kozcaz and L.G. Yaffe,
Energy loss of a heavyquark moving through N = 4 supersymmetric Yang-Mills plasma, JHEP 07 (2006) 013[hep-th/0605158].[25] S.S. Gubser,
Drag force in AdS/CFT,
Phys. Rev. D 74 (2006) 126005 [hep-th/0605182].[26] J. Casalderrey-Solana and D. Teaney,
Heavy quark diffusion in strongly coupled N = 4 Yang Mills,
Phys. Rev. D 74 (2006) 085012 [hep-ph/0605199].[27] G. Torrieri and J. Noronha,
Flavoring the Quark-Gluon Plasma with Charm,
Phys. Lett.B 690 (2010) 477 [1004.0237].[28] Y. Burnier and M. Laine,
Massive vector current correlator in thermal QCD,
JHEP 11(2012) 086 [1210.1064].[29] S.L. Adler, J.C. Collins and A. Duncan,
Energy-Momentum-Tensor Trace Anomaly inSpin 1/2 Quantum Electrodynamics,
Phys. Rev. D 15 (1977) 1712.[30] J.C. Collins, A. Duncan and S.D. Joglekar,
Trace and Dilatation Anomalies in GaugeTheories,
Phys. Rev. D 16 (1977) 438.[31] M. Laine and Y. Schr¨oder,
Quark mass thresholds in QCD thermodynamics,
Phys. Rev.D 73 (2006) 085009 [hep-ph/0603048]. 2132] M. Cheng [RBC-Bielefeld Collaboration],
Charm Quarks and the QCD Equation of State,
PoS LAT2007 (2007) 173 [0710.4357].[33] C. DeTar, L. Levkova, S. Gottlieb, U.M. Heller, J.E. Hetrick, R. Sugar and D. Toussaint,
QCD thermodynamics with nonzero chemical potential at N t = 6 and effects from heavyquarks, Phys. Rev. D 81 (2010) 114504 [1003.5682].[34] S. Borsanyi, G. Endrodi, Z. Fodor, S.D. Katz, S. Krieg, C. Ratti, C. Schroeder andK.K. Szabo,
The QCD equation of state and the effects of the charm,
PoS LATTICE2011 (2011) 201 [1204.0995].[35] D. B¨odeker and M. Laine,
Heavy quark chemical equilibration rate as a transport coeffi-cient,
JHEP 07 (2012) 130 [1205.4987].[36] M. Laine, M. Veps¨al¨ainen and A. Vuorinen,
Ultraviolet asymptotics of scalar and pseu-doscalar correlators in hot Yang-Mills theory,
JHEP 10 (2010) 010 [1008.3263].[37] M. Laine, A. Vuorinen and Y. Zhu,
Next-to-leading order thermal spectral functions inthe perturbative domain,
JHEP 09 (2011) 084 [1108.1259].[38] H.B. Meyer,
The Bulk Channel in Thermal Gauge Theories,
JHEP 04 (2010) 099[1002.3343].[39] H.-T. Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz and W. Soeldner,
Charmoniumproperties in hot quenched lattice QCD,
Phys. Rev. D 86 (2012) 014509 [1204.4945].[40] M. Nishimura and Y. Schr¨oder,
IBP methods at finite temperature,
JHEP 09 (2012) 051[1207.4042].[41] J. Beringer et al. [Particle Data Group Collaboration],
Review of Particle Physics (RPP),
Phys. Rev. D 86 (2012) 010001.[42] T. van Ritbergen, J.A.M. Vermaseren and S.A. Larin,
The four-loop β -function in Quan-tum Chromodynamics, Phys. Lett. B 400 (1997) 379 [hep-ph/9701390].[43] J.A.M. Vermaseren, S.A. Larin and T. van Ritbergen,
The 4-loop quark mass anomalousdimension and the invariant quark mass,
Phys. Lett. B 405 (1997) 327 [hep-ph/9703284].[44] W. Florkowski and B.L. Friman,
Spatial dependence of the finite temperature mesoncorrelation function,
Z. Phys. A 347 (1994) 271.[45] H.B. Meyer,
Energy-momentum tensor correlators and spectral functions,
JHEP 08(2008) 031 [0806.3914]. 2246] S. Caron-Huot,
Asymptotics of thermal spectral functions,
Phys. Rev. D 79 (2009) 125009[0903.3958].[47] D.J. Broadhurst, J. Fleischer and O.V. Tarasov,
Two-loop two-point functions withmasses: Asymptotic expansions and Taylor series, in any dimension,
Z. Phys. C 60(1993) 287 [hep-ph/9304303].[48] Y. Burnier, M. Laine and M. Veps¨al¨ainen,
Heavy quark medium polarization at next-to-leading order,
JHEP 02 (2009) 008 [0812.2105].[49] P. Petreczky and D. Teaney,
Heavy quark diffusion from the lattice,
Phys. Rev. D 73(2006) 014508 [hep-ph/0507318].[50] Y. Burnier and M. Laine,