Quarkonium spectral function in medium at next-to-leading order for any quark mass
aa r X i v : . [ h e p - ph ] O c t Quarkonium spectral function in medium atnext-to-leading order for any quark mass
Yannis Burnier
Institute of Theoretical Physics, EPFL, CH-1015 Lausanne, Switzerland
October 18, 2018
Abstract
The vector channel spectral function at zero spatial momentum is cal-culated at next-to-leading order in thermal QCD for any quark mass. Itcorresponds to the imaginary part of the massive quark contribution tothe photon polarization tensor. The spectrum shows a well defined trans-port peak in contrast to both the heavy quark limit studied previously,where the low frequency domain is exponentially suppressed at this orderand the naive massless case where it vanishes at leading order and divergesat next-to-leading order. From our general expressions, the massless limitcan be taken and we show that no divergences occur if done carefully.Finally, we compare the massless limit to results from lattice simulations.
In heavy ion collisions, fireballs of quark-gluon plasma are formed. The QCDmatter being strongly interacting, the quarks and gluons only escape as mesonsor hadrons when the plasma cooled down sufficiently. To reconstruct whathappened at early stages of the collision, we have to resort to probes that canbe traced in experiment and whose properties are modified in the plasma. Forseveral available probes, one quantity describes their fate in the plasma: thevector channel spectral function in medium.For instance the way bound states such as charmonium or bottomoniumdecay into muon pairs [1, 2] could tell us how they are affected by the plasma[3, 4]. Another example is the heavy flavor diffusion coefficient, which describeshow the massive quarks will be diffused or slowed down by the plasma. Thisis another observable which is of interest for experiment, in fact the heavyquarks can be tagged and their transverse momentum distribution studied (seefor instance [5, 6]. For light quarks, the zero frequency limit of the spectralfunction describes the electric conductivity [7] and its momentum depencancethe photon or dilepton emission rate [8, 9, 10].In the vacuum, this spectral function is known at five loops [11] for masslessfermions and Taylor expansions in the mass are known to four loops [12, 13,14, 15]. In the presence of a finite temperature medium, the two loop masslessresult is known since a long time [16] and the case of large masses with respectto the temperature M ≫ T was calculated rather recently [17]. Here we extendthe previous calculations to any mass and discuss the transport part of the1pectrum, which is suppressed in the heavy quark limit and was not obtainedin the previous calculations. We still consider implicitly that the frequency ω ≫ gT is sufficiently large so that we can neglect hard thermal loop (HTL)corrections [18] in the fermion propagators and vertices. Other HTL correctionsare of higher order [17] and will not be addressed here either.Of course in QCD, the convergence of perturbation theory is slow due to thelargeness of the coupling α s and moreover finite temperature gauge theories suf-fer from infrared problems so that the full infrared physics requires nonperturba-tive methods. Lattice computations contain the full physics but are performedin Euclidean time and do not have a direct access to the Minkowskian spectralfunction. After measuring the corresponding Euclidean correlator, an analyticcontinuation is needed to obtain the desired spectral function. In the case ofdiscrete numerical data of finite precision the reconstruction of the spectrum isvery hard to perform [19, 20, 21, 22]. The challenge is even bigger here since theEuclidean correlator is not even continuous at zero Euclidean times and hencethe full analytical continuation is ill defined [23]. That’s where perturbationtheory could again be of use, since the zero Euclidean time limit (or the corre-sponding large frequency limit of the spectral function) is weakly coupled andaccessible to perturbation theory. This ’large’ perturbative part (containingzero and possibly finite temperature contributions) could be subtracted fromthe lattice data [24] or used as a prior to define the analytical continuation [25].In the case of the vector current spectral function considered here, the Eu-clidean corralator was computed recently together with its mass dependence in[26]. In this paper we complete our program and calculate the related spec-tral function i.e. perform the analytic continuation. This is not a trivial taskeven though the Euclidean correlator G ( τ, p ) of ref. [26] can be cross checkedby convoluting the spectral function ρ ( ω, p ) with the finite temperature kernel K ( τ, ω ): G ( τ, p ) = Z ∞ dωπ ρ ( ω, p ) K ( τ, ω ) , K ( τ, ω ) = cosh( ω ( τ − β/ ωβ/ . (1)After defining the observables in section 2, we discuss the calculation in sec3. and refer to appendices for details. In section 4 we present our results forthe spectrum, discuss the transport coefficients and derive the massless limit,which we compare to lattice results. Conclusions are given in section 5. We consider the vector current of a massive quark described by the operatorΨ( τ, x ) J µ ( τ, x ) = ¯Ψ( τ, x ) γ µ Ψ( τ, x ) , (2)with µ = 0 , .., d . The object we compute here is the spectral function in medium,which is given by the thermal average of the current commutator ρ V ( ω ) = Z dt e iωt Z d d x h
12 [ J µ ( t, x ) , J µ (0 , i T . (3)2ollowing [17] we will start from the Euclidean correlator in frequency space( ω n = 2 πnT denote the Matsubara frequencies), G E ( ω n ) = Z β dτ e iω n τ Z d d x h J µ ( τ, x ) J µ (0 , i , (4)which can be calculated using conventional finite temperature Feynman rules[7, 27, 28]. The spectral function can be determined from the discontinuity ofthe Euclidean correlator along the imaginary axis: ρ V ( ω ) = Disc [ G E ( − iω )] = 12 i lim ε → + [ G E ( − iω + ε ) − G E ( − iω − ε )] . (5)Note that for ω ∼ gT usual perturbation theory breaks down and would requireto use Hard Thermal Loop (HTL) resummation, which will not be consideredhere. As was shown in [17], no infrared divergences occur so that HTL correc-tions are subleading for ω ≫ gT . One observable defined by this spectral function is the production rate of muonpairs from the decay of massive quark pairs [1, 2]. If we suppose that the quarkpair is at rest we have in particular: dN µ ¯ µ d xd q = − e Z π ) ω m µ ω ! − m µ ω ! n B ( ω ) ρ V ( ω ) , (6)where Ze is the charge of the quark and n B the Bose-Einstein Distribution and ω = E µ + + E µ − . The main contribution to this observable comes form thethreshold ω ∼ M , where the spectrum can be obtained more precisely withdedicated resummations [17, 29].When speaking of heavy flavor diffusion or electric conductivity, the diffusioncoefficient D is obtained from the low energy behavior of the spectrum. In factone expects the spectral function to look like a Lorentzian in the low energylimit − ρ V ( ω ) ω <ω<ω UV ≈ χD η D η D + ω , (7)where χ is the susceptibility, η D another number called the drag coefficient and ω UV the scale at which other kind of physics enter and where the spectrumdeviates from a Lorentzian . The diffusion coefficient can then be extracted as D = − χ lim ω → ρ V ( ω ) ω . (8)For a thorough discussion of this formula and the zero frequency limit see forinstance ref. [7].If the onset of non-transport physics ω UV is well separated from the trans-port peak, another strategy can be used [30]. The idea is to calculate anotherobservable, the momentum diffusion coefficient κ , which is proportional to the Note also that ρ V ( ω > < η D . It can be extracted from the falloff of the Lorentzian peak[31]: κ = 2 M kin T η D ≈ − M kin ω χ ρ V ( ω ) ω (cid:12)(cid:12)(cid:12)(cid:12) η D ≪ ω ≪ ω UV , (9)where M kin is the in medium kinetic mass defined by the low momentum limitof the dispersion relation. The fluctuation-dissipation theorem finally relatesthis coefficient to the diffusion coefficient D = 2 T /κ and hence to the dragcoefficient η D = κ/ (2 M kin T ).The momentum diffusion coefficient can be calculated in perturbation theorywith dedicated resummations [32]. However even if the resummations seem tocatch the relevant physics, the convergence of the perturbative series for κ isat best very slow. In the case of the heavy quarks for instance, the first non-vanishing contribution arises at O ( α ) and the correction O ( α g ) is an order ofmagnitude larger for typical heavy ion plasmas [33]. We now turn to the calculation of the spectral function, following refs. [17, 26,29, 34].
Performing the Wick contractions and the trace algebra, we get at leading order(LO): G VE = 2 N c PZ { P } (cid:18) − − ǫ )∆( P ) + − M + 2 Q (1 − ǫ )∆( P )∆( P − Q ) (cid:19) so that ρ VLO ( ω ) = Disc [ G VE ( − iω )] (10)= 2 N c Disc [ PZ { P } (cid:18) − − ǫ )∆( P ) + − M + 2 Q (1 − ǫ )∆( P )∆( P − Q ) (cid:19) ] , where Q = ( − iω ± ǫ, ) will be set in the process of taking the discontinuity. Notethat the first term is independent of Q and will not contribute to the spectrum.To simplify the following expressions in both the LO and the next-to-leadingorder (NLO), we introduce the following notations: I ij = Disc PZ { P } P ) i ∆( P − Q ) j , (11) I nijklm = Disc PZ K, { P } ( K · Q ) n ( K ) i ∆( P ) j ∆( P − K ) k ∆( P − Q ) l ∆( P − K − Q ) m with ∆( P ) = P − M , so that at leading order, ρ VLO ( ω ) = − N c (4 M + 2 ω ) I ( ω ) . (12) Namely the velocity dependence of the free energy is expanded as F ( v ) = M rest + M kin v / O ( v ). I ( ω ) = θ ( ω − M ) √ ω − M πω tanh (cid:16) ω T (cid:17) (cid:20) ε (cid:18) µ ω − M (cid:19)(cid:21) + πωδ ( ω ) Z p n ′ F ( E p )2 E p , (13)where we introduced the notation E p = p + M , R p = R d d p (2 π ) d . After performing the Wick contractions and the trace algebra, we get the NLOin terms of master sum-integrals defined in equations (11): ρ VNLO N c C F g = 4(1 − ǫ ) I − − ǫ ) I + 8(1 − ǫ ) M I − − ǫ ) I − − ǫ ) I + 8(1 − ǫ ) I − M + ω (1 − ǫ )) I + 4(1 − ǫ ) I +8 M (cid:0) M + ω (1 − ǫ ) (cid:1) I − − ǫ ) (cid:0) M + ω (1 − ǫ ) (cid:1) I − − ǫ ) I − + 2 (cid:0) M ǫ + ω (cid:0) − ǫ − ǫ (cid:1)(cid:1) I +2 (cid:0) M − M ω ǫ − ω (1 − ǫ ) (cid:1) I − − ǫ ) I + 4(1 − ǫ ) (cid:0) M + ω (1 − ǫ ) (cid:1) I . (14)Note that the first four terms are independent of ω and do not contribute tothe spectral function. The previous NLO expression (14) is UV divergent but is finite after redefinitionof the mass. The counterterms for the currents read: ρ V,CTNLO N c C F g = δM g C F N c ∂ρ VLO ∂M (15)with 14 N c ∂ρ VLO ∂M = (4 M + 2 ω (1 − ǫ )) I − I − − ǫ ) I (16)and using the pole mass scheme as in ref. [26], we set δM = − g C F M (4 π ) (cid:18) ǫ + ln ¯ µ M + 43 (cid:19) = − g C F Z k (2 − ǫ )( E pk − k ) + M ∆ − + + M ∆ ++ kE pk , (17)where we denoted E pk = E p + k , k = | k | and ∆ ±± = k ± E p ± E pk .5 .4 Thermal correction to the mass After subtraction of the counterterms, the spectral function is finite everywherebut at the threshold ω = 2 M . The divergence there comes from thermal cor-rections, in fact one can rewrite the whole spectral function as ρ V ( ω, M ) = ρ VLO ( ω, M ) + δM T ∂ρ VLO ( ω, M ) ∂M + ¯ ρ VNLO ( ω, M ) + O ( g ) , (18)where ¯ ρ VNLO = ρ VNLO − δM T ∂ρ VLO ∂M (19)and the term δM T ∂ρ
VLO ( ω,M ) ∂M is actually responsible for the divergence. In factwe can resum this contribution by redefining the mass M → M + δM T : ρ V ( ω, M ) = ρ VLO ( ω, M + δM T ) + ¯ ρ VNLO ( ω, M + δM T ) + O ( g ) . (20)The explicit shift δM T , is the thermal contribution to the dispersion relation,which, for a massless fermion is the well-known δM T = g C F Z ∞ dkπ k ( n B ( k ) + n F ( k )) = g C F T δM T g C F = 2 Z k n B ( k ) k + n F ( E pk ) E pk (cid:18) − M ∆ ++ ∆ −− − M ∆ + − ∆ − + (cid:19) (22)= Z ∞ dkπ (cid:18) k n B ( k ) + k E k (cid:18) M pk ln (cid:12)(cid:12)(cid:12)(cid:12) k − pk + p (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) n F ( E k ) (cid:19) , which actually depends on the integration variable p .As a summary, to avoid divergences at the threshold, we resum the masscorrection. To calculate the spectral function of a fermion of mass squared M ,we calculate the spectral function at the mass M + δM as written in equation(20) and modify the NLO contribution as explained in equation (19).Note that the relevant bosonic part of this shift was performed in [17]. Thedivergence is however integrable and the shift was left out in ref. [26], where theEuclidean correlator was calculated, but the changes are easily tractable (seeappendix D). While the leading order is given in equations (12,13) the next-to-leading orderrequires significantly more work. The full expression is obtained adding to theNLO (14), the mass counterterm (15) and the contribution from the thermalmass shift (19). The first step is to carry out the sums (see Appendix A) andtake the discontinuity (5). We are then left with the integrals and a delta func-tion remaining from the discontinuity. Parts of the integrals can be performedanalytically and the remaining ones have to be done numerically. The explicitexpressions for the master integrals are given in appendix B and the details ontheir integration in appendix C. 6
Results
The final result can be split into three parts ρ VNLO ( ω ) = ρ vacNLO ( ω ) + ρ bosNLO ( ω ) + ρ fermNLO ( ω ) . (23)First the vacuum part [13, 14, 15] that can be integrated explicitly ρ vacNLO ( ω )4 N c C F g = 2 θ ( ω − M )(4 π ) ω (cid:20) (4 M − ω ) L ω − √ ω − M ω + √ ω − M ! (24)+(7 M + 2 M ω − ω )acosh (cid:16) ω M (cid:17) + ω p ω − M × (cid:18) ( ω + 2 M ) ln ω ( ω − M ) M −
38 ( ω + 6 M ) (cid:19)(cid:21) , where L = 4Li ( x ) + 2Li ( − x ) + ln( x ) [2 ln(1 − x ) + ln(1 + x )]. Secondly thefirst thermal part, that we will denote ’bosonic’ thermal correction, calculated inref. [17] for which one integral is left for numerical evaluation (for mass shift con-tribution see Appendix D). It is proportional to the Bose-Einstein distributionfunction n B ( k ) and does not contain any Fermi-Dirac distribution: ρ bosNLO N c C F g = 2(4 π ) ω Z ∞ dk n B ( k ) k (cid:26) θ ( ω ) θ (cid:18) k − M − ω ω (cid:19) (cid:20) (25)2 ω k s − M ω ( ω + 2 k ) + ( ω + 2 M ) p ω ( ω + 2 k ) p ω ( ω + 2 k ) − M − (cid:16) ω − M + 2 ωk ( ω + 2 M ) + 2 ω k (cid:17) acosh r ω ( ω + 2 k )4 M (cid:21) + θ ( ω − M ) θ (cid:16) ω − M ω − k (cid:17)(cid:20) ω k s − M ω ( ω − k )+( ω + 2 M ) p ω ( ω − k ) p ω ( ω − k ) − M − (cid:16) ω − M − ωk ( ω + 2 M ) + 2 ω k (cid:17) acosh r ω ( ω − k )4 M (cid:21) + θ ( ω − M ) (cid:20) − ω + 2 M ) ω p ω − M +4 (cid:16) ω − M + 2 ω k (cid:17) acosh (cid:18) ω M (cid:19)(cid:21)(cid:27) . The remaining contribution ρ fermNLO contains at least one Fermi-Dirac distribution n F ( E p ) or n F ( E pk ) hence it is suppressed in the limit M ≪ T . However it isthe only piece containing a non-vanishing transport peak near ω → ω < M only a few structures actually contribute and wecan quote the result here: ρ fermNLO ω< M = − n F (cid:16) ω (cid:17) ρ bosNLO + ρ ferm, NLO , (26)7 Ω (cid:144) T - Ρ (cid:144) Ω ² Spectrum, for different masses at T = T c M = TM = TM = T M = T M = T Ρ NLO Ρ LO - - (cid:144) T li m Ω ® Ρ H Ω L g ² T Ω Zero energy limit of the spectrum
Figure 1: (Left): LO (green dashed) and NLO (blue) spectral function( − ρ V ( ω ) /ω ) for M = (0 . , , . , , T , T = 1 . T c . The LO vanishes atfrequencies below the threshold ω = 2 M where the NLO shows a discontinuity.(Right) Value of the low energy limit of ρ V ( ω ) /ω as function of the quark mass.where ρ ferm, NLO N c C F g = 18 π Z ∞ ω/ Z ∞ E p ( (27)1 + M (cid:0) M + ω (cid:1) ω (cid:18) E p + 2 k − ω ) + 1(2 E p − ω ) (cid:19) + − k ω + (cid:0) kω + 2 M (cid:1) − ω ( ω − k ) kω (cid:18) E p + 2 k − ω − E p − ω (cid:19)) × [ n B ( k )( n F ( E p ) − n F ( E p + k − ω )) + ( n F ( E p ) − n F ( E p + k − ω )] . This last term ρ ferm, NLO contribute to the transport peak and comes form gluonemission or absorption by the massive fermion.The full result containing the LO and NLO contributions is plotted infig. 1 (left) for T = 1 . T c and different masses ranging from M = 0 . T to M = 11 . T corresponding to a bottom quark of mass 4.65GeV. Note that thespectral function is in fact negative and for clarity we show − ρ V ( ω ) in theplots. In fig. 1 (right) we show the zero frequency limit of the spectral functionlim ω → ρ V ( ω ) /ω , which enters the determination of the transport coefficient D .The transport coefficient itself requires division by the suscptibility χ , for whichwe use the NLO result of ref. [26], see fig. 2. Separate results for the three con-tributions to the next to leading order are shown in fig. 2 for the case M = 3 . T ,representing the charm quark at T = 1 . T c studied on the lattice in [36] and forthe generic case M = T in fig. 3. In these figures we scaled the result to ω sothat the vacuum contribution goes to a constant at large frequency. In the insert of the same figures 2, 3, we scale the spectrum to the frequencyand zoom on the low frequency region to see the transport peak. We see that ρ ( ω ) ω goes smoothly to a constant with vanishing slope at ω → D according to8 .5 1.0 5.0 10.0 50.0 - Ω (cid:144) T - Ρ (cid:144) g Ω M = T c Ω (cid:144) T - Ρ (cid:144) g Ω T Ρ NLO Ρ NLOvac Ρ NLObos Ρ NLOferm (cid:144) T Π T D (cid:144) g Transport coefficient
Figure 2: (Left): Different parts of the NLO contribution to the spectral functionfor M=3 . T and Lorentzian fit of the transport peak (black dot-dashed line)(Right): The generated transport coefficient D .equation (8). The value of D is plotted as a function of the quark mass in fig.2(right). We see that it is suppressed for large masses T D ∝ ( M/T ) − andbehaves as a power law T D ∝ ( M/T ) − at small M . In the case of the charmquark at T = 1 . T c , we see that it is in principle small 2 πT D ≈ . πT D ≈ πT D ≈
10 [33].That said, the transport peak is not well separated from the UV physics - atleast in this low order calculation - and there is no region where η D ≪ ω ≪ ω UV so that the momentum diffusion coefficient cannot be defined straight fromformula (9). Of course for the only purpose of defining κ or η D one could justfit the low frequency part with equation (7) and get η D as a fit parameteradjusted to the curvature of ρ ( ω ) ω around ω = 0. In the previous case we get η D ∼ . M kin = M , into 2 πT D ∼ .
8. This differ fromthe direct estimate made above (2 πT D ≈ .
1) showing that the fluctuation-dissipation theorem does not seem to apply here. Note again that the presentcomputation is not systematic [37] for ω < gT so that no strong conclusionshould be made with this remark.
Let’s we consider a fermion that is massless in the vacuum. Even if in this caseHTL correction would be of the same order as our NLO result, it is interestingto see how our result behaves. As we resummed the thermal mass correctionaccording to equation (20), we in fact have to calculate the spectral function ofa fermion of mass δM T . In a typical quark-gluon plasma produced in heavy ioncollision we have numerically (21) δM T = g C F T / ≈ T . The spectrum ofsuch a fermion is shown in fig. 3. We see that the spectrum has no divergenceneither at the threshold nor at zero frequency and shows a transport peak. Thisresult differs from the old result of ref. [16]. There, the mass shift was noticedand performed in the LO result but was not made in the NLO calculationwhere the mass was set to zero. Their NLO contribution hence diverges atzero frequency, whereas ours has a threshold structure at ω = 2 δM T and a9 .5 1.0 5.0 10.0 50.0 - - Ω (cid:144) T - Ρ (cid:144) g Ω M = T0.0 0.5 1.0 1.50.000.010.020.030.040.050.06 Ω (cid:144) T - Ρ (cid:144) g Ω T Ρ NLO Ρ NLOvac Ρ NLObos Ρ NLOferm Ω (cid:144) T - Ω Ρ N L O T (cid:144) g T M = T Figure 3: (Left): The different contributions to the NLO spectrum for M = T . (Right): The result of this paper for the total thermal NLO contribution(continuous blue line) compared to the result of ref. [16] where the mass shiftwas not fully taken into account (dotted red line). Note that the thresholdstructure appears negative as we show only the thermal part of the NLO. Forthe complete spectrum see fig. 1transport peak, see fig. 3. At high frequency both results agree and merge tothe asymptotic behavior derived in ref. [38], which reads − ρ TNLO = − ρ bosNLO − ρ fermNLO ω ≫ M,T = 4 N c C F g (cid:18) πT ω + T M πω (cid:19) . (28)The transport coefficient obtained here is of order 2 πT D ∼ .
3, which isagain on the low side, the perturbative resummation from [39, 40] gets in thiscase 2 πT D ∼
25 and lattice results ranges from 2 πT D ∼ − The full spectral function for M = δM T = g C F T / T = 1 . T c is shownin fig. 4 together with the Euclidean correlator scaled to the free correlator [8].Keeping in mind that our approximations are not consistent at low frequencies,we still compare the Euclidean correlator to the lattice data of ref. [41] fora massless fermion. Apart from an overall normalization we see a differentslope at large τ signaling a lack of power in the small ω region. This couldbe compensated for by an additional transport peak. Keeping for instance thesame width η D = 0 .
7, if we add an additional − ρ t /ω = Dχη D ω + η D for ω < δM T with 2 πT ¯ D = 0 .
05 we get the dashed curve on fig. 4b). We see that it now goesnicely parallel to the lattice data. The total diffusion coefficient would then be2 πT D ∼ .
4, still on the low side. The asymptotic behavior of ref. [38] was derived without the mass shift hence containsonly the first term, see also footnote 6 there Note that the normalization of the perturbative curve could be improved by adding higherorder in the vacuum spectral function. - Ω (cid:144) T Ρ (cid:144) Ω ² Spectrum, M = Ω (cid:144) T Ρ (cid:144) Ω T Ρ NLO Ρ LO ΡΡ+Ρ t Lattice Τ G ii H Τ L (cid:144) G fr ee ii H Τ L Euclidean Correlator G ii H Τ L Figure 4: (Left): Spectrum at leading and next-to-leading order for M = 0 . T .In the insert, we fit the low energy spectrum with a Lorentzian. (Right): Next-to-leading order Euclidean correlator for M = T (blue continuous line) togetherwith lattice data (dots). We note that the agreement can be greatly improvedadding more power to the transport peak, even if some normalization factorremains. (magenta dashed line) We calculated the thermal correction to the massive quark vector spectral func-tion at NLO in thermal QCD. The thermal corrections are small in comparisonto the vacuum corrections for large fermion masses and comparable to them for M ∼ T .The result shows some typical features: First, the threshold for pair pro-duction is smoothed by thermal corrections, the discontinuity in the vacuumspectrum being partly compensated by thermal effects. Secondly a transportpeak appears at this order, it is however small in comparison to other expecta-tions. We also see that the transport peak is broad and is not well separatedfrom other kind of physics, rendering its determination via the fluctuation-dissipation theorem difficult. It should however be stressed that higher looporders give corrections of the same order for ω ≤ gT .Setting the vacuum fermion mass to zero in our formulas do not lead todivergences. This contradict the calculation of ref. [16] where the spectral func-tion was calculated for massless fermions and the result for the NLO diverges atzero frequency hindering a definition of the corresponding Euclidean correlator.The difference can be traced back to how the thermal mass shift is introducedin the NLO. In the previous reference, the thermal mass shift was performed inthe leading order part as done here but not in the NLO contribution where thefermion remained purely massless, leading to divergences at zero frequency. Acknowledgements
The author would like to thank M. Laine for helpful discussions and imporantcomments on the manuscript. This work was supported by the SNF grantPZ00P2-142524. 11
Calculation of the master integrals
The calculation of the master sum-integrals is mostly standard (for details seeref. [28]), only one difficulty arises in double poles at zero frequency, which willbe explained below. Their calculation is subtle and is shown explicitly for thecase of I , which is the simplest from the amount of algebra but containsmost of the technical difficulties.We start by rewriting the master integral in a form where the elementarysummation formula (valid for 0 < τ < β ) T X p n e ip n τ p n + E p = n B ( E p )2 E p (cid:16) e τE p + e ( β − τ ) E p (cid:17) T X { p n } e ip n τ p n + E p = n F ( E p )2 E p (cid:16) e τE p − e ( β − τ ) E p (cid:17) (29)can be applied. To do that we shift P − K → K and introduced the newintegration variables S, R and the corresponding delta functions: I = Disc PZ K, { P } P )∆( P − K )∆( P − Q )∆( P − K − Q ) = Disc PZ { P,K,S,R } δ ( R − P + Q ) δ ( S − K + Q )∆( P )∆( K )∆( R )∆( S ) . (30)The temporal part of the delta function can be rewritten as an integral, where,we keep track of time ordering: δ ( r n − p n + q n ) δ ( s n − k n + q n ) (31)= e iβ ( p n + s n ) Z β dτ e i ( r n − p n + q n ) τ Z τ dτ e i ( − s n + k n − q n ) τ (32)+ e iβ ( p n + s n + q n ) Z β dτ e i ( − s n + k n − q n ) τ Z τ dτ e i ( r n − p n + q n ) τ . (33)The factors e iβ ( p n + k n ) = e iβ ( p n + s n + q n ) = 1 were added so that the total phaseof all the Matsubara frequencies are between 0 and β . Thanks to the timeordering and this last prescription the sums can be performed with formula(29) and then the integrals over τ , τ calculated. Remembering that q n is abosonic Matusubara frequency we can set e iq n β → . The result is a ratherlong expression containing many terms that can be simplified in rewriting allthe exponents in terms of Fermi-Dirac (or Bose-Einstein) distributions. Eachterm contains a product of Fermi-Dirac (or Bose-Einstein) distributions in thenumerator and products of all kinds of sums of the energies ( a p E p + a k E k + a r E r + a s E s + ia q q n ) with integer a i ’s in the denominator. Note that thereare no divergences (up to the poles for the q n variable on the complex axis)since when one sum of the energies ( a p E p + a k E k + a r E r + a s E s ) vanishes in This simplification has to be done, it is part of the prescription to get the correct analyticcontinuation. f ( E p , E r , E k , E s , ω )( ω − E p + E r )( ω + E s − E k ) . (34)This term contains a simple pole at ω = 0 when we enforce one of the deltafunctions contained in (30). Note that its contribution is finite even after re-placing the second delta function as the numerator f vanishes when energiesare set equal. However this term also contains a double pole when both deltasare set to zero. The simplest way to deal with this issue is to rewrite the deltafunctions in (30) as δ ( r − p ) δ ( k − s ) ∝ δ ( E r − E p ) δ (Ω r − Ω p ) δ ( E s − E k ) δ (Ω s − Ω k ) ∝ δ ( E r − E p ) δ (Ω r − Ω p ) δ ( E s − E k − E r + E p ) δ (Ω s − Ω k ) (35)and use the δ ( E s − E k − E r + E p ) to replace E s − E k by E r − E p (note that thisreplacement has to be done as a limit). After this, all the poles which contributeat ω ∼ ω ± ( E r − E p )). The discontinuity can be taken ineach term separately usingDisc (cid:20) ω − A (cid:21) = − πδ ( ω − A ) (36)for simple poles and Disc (cid:20) ω − A ) (cid:21) = πδ ′ ( ω − A ) (37)for double poles.The spatial integrals over r , s can be performed using the remaining deltafunctions δ ( r − p ) δ ( s − k ). They force us to perform the limit E r → E p andget as final result: I ω ∼ = Z p,k ( πωδ ( ω ) n ′ F1 E p E pk (1 − n F2 ) + n ′ F1 n ′ F2 E p E pk !) . (38) B Master integrals
Following the notations of ref. [17], we detail here the different master integrals(I denote E p = p + M , E k = k + M , E p − k = ( p + k ) + M , k being thegluon momentum). First the leading order sum-integrals: I = π Z p (cid:20) ωδ ( ω ) n ′ F1 E p + ( δ ( ω − E p ) − δ ( ω + 2 E p )) 1 − n F1 E p (cid:21) , (39) I = π Z p (1 − ǫ ) (cid:20) ωδ ( ω ) n ′ F1 p E p + ( δ ( ω − E p ) − δ ( ω + 2 E p )) 1 − n F1 p E p (cid:21) . I = Z k n B ( k )2 k I , (40) I = Z k n B ( k )2 k I , (41) I = Z k − n F ( E k )2 E k I , (42) I = Z k − n F ( E k )2 E k I . (43)The truly NLO sum-integrals are given below using the notation ∆ ± = E p ± E pk and ∆ στ = k + σE p + τ E pk . Here I only give the terms proportional to ωδ ( ω )when the other terms contributing at ω > I ω ∼ = I ω ∼ = 0 , (44) I ω ∼ = πωδ ( ω ) Z k,p n ′ F1 E p E pk k (cid:26) [1 + n B0 − n F2 ] (cid:18) − + + 1∆ ++ (cid:19) − [ n B0 + n F2 ] (cid:18) −− + 1∆ + − (cid:19)(cid:27) , (45) I ω ∼ = 0 , (46) I ω ∼ = πωδ ( ω ) Z k,p (cid:20) ( k − E p − E pk ) n ′ F1 n ′ F2 E p E pk ∆ ++ ∆ − + ∆ + − ∆ −− + n ′ F1 E pk E p k (cid:26) E pk (cid:18) ∆ ++ + E pk ∆ + ∆ − + + E pk ∆ − + (cid:19) (47)+ E pk n B0 (cid:18) − ∆ + 1∆ −− ∆ − + (cid:19) + n F2 k E pk − E p + E pk − k ∆ −− ∆ + 2∆ − − E p + E pk − k ∆ − + ∆ − !(cid:27)(cid:21) . The last sum-integrals are given for all ω for completeness (as they were not allwritten in a suitable form for our present purpose): I = π Z p,k " ωδ ( ω ) ( n ′ F1 E p E pk (1 − n F2 ) + n ′ F1 n ′ F2 E p E pk ) (48) − ( δ ( ω − E p ) − δ ( ω + 2 E p )) [1 − n F1 ][1 − n F2 ]8 E p E pk ∆ + ∆ − ,I − = 2 I + 12 (cid:0) ω − M (cid:1) I (49)14nd I = π Z k,p " ωδ ( ω ) E p E pk k [ n ′ F1 − E p n ′′ F1 ] × (cid:26) [1 + n B0 − n F2 ] (cid:18) − + + 1∆ ++ (cid:19) − [ n B0 + n F2 ] (cid:18) −− + 1∆ + − (cid:19)(cid:27) − n ′ F1 E p E pk k (cid:26) n B0 − n F2 ∆ ++ ∆ − + (cid:18) − + + 1∆ ++ (cid:19) − n B0 + n F2 ∆ + − ∆ −− (cid:18) −− + 1∆ + − (cid:19)(cid:27)! (50)+ [ δ ( ω − E p ) − δ ( ω + 2 E p )]32 k E p ( k + (1 − ǫ ) E pk − E p + 2 M ) − E pk p E p E pk p × (1 − n F1 ) (cid:26) [1 + n B0 − n F2 ] (cid:18) − + + 1∆ ++ (cid:19) − [ n B0 + n F2 ] (cid:18) −− + 1∆ + − (cid:19)(cid:27) − E p + E pk − k − M E p E pk p (1 − n F1 ) n ′ F2 (cid:20) ++ + 1∆ + − + 1∆ − + + 1∆ −− (cid:21) + 1 − n F1 E p E pk ((cid:20) ++ + 1∆ − + (cid:21) (cid:20) ++ − − + + 1 E p (cid:21) [1 + n B0 − n F2 ]+ (cid:20) −− + 1∆ + − (cid:21) (cid:20) −− − + − − E p (cid:21) [ n B0 + n F2 ] ) − − n F1 E p E pk ( ∆ + E p + E pk − E p − k p ! (cid:20) n B0 − n F2 ∆ + n B0 + n F2 ∆ −− (cid:21) + ∆ − E p + E pk − E p − k p ! (cid:20) n B0 − n F2 ∆ − + + n B0 + n F2 ∆ − (cid:21))! + [ δ ( ω − ∆ ++ ) − δ ( ω + ∆ ++ )] (1 + n B0 )(1 − n F1 − n F2 ) + n F1 n F2 E p E pk k ∆ ∆ − + + [ δ ( ω − ∆ −− ) − δ ( ω + ∆ −− )] − n B0 (1 − n F1 − n F2 ) + n F1 n F2 E p E pk k ∆ −− ∆ − + [ δ ( ω − ∆ + − ) − δ ( ω + ∆ + − )] n B0 n F1 − (1 + n B0 ) n F2 + n F1 n F2 E p E pk k ∆ −− ∆ − + [ δ ( ω − ∆ − + ) − δ ( ω + ∆ − + )] n B0 n F2 − (1 + n B0 ) n F1 + n F1 n F2 E p E pk k ∆ ∆ − + . C Full NLO result
To obtain the full NLO result form the previous formulas, the integrals over p , k still have to be performed. From the antisymmetry of the spectral function itis enough to consider the case ω ≥
0. Using rotation symmetry, it is obviousthat all the angular integrals are trivial up to the one containing the angle θ between p an k so that we have three non-trivial integrals to perform over | p | , | k | , cos θ . The full result can be split into three parts, a part proportional15o ωδ ( ω ) which we call ρ tNLO for transport, a part proportional to δ ( ω − E p ),called ρ fNLO for factorized and a last part proportional to one of the δ ( ω ± ∆ ±± )denoted by ρ pNLO for phase space, following the notation of [17]. For the termsproportional to δ ( ω ) the three integrals have to be performed but the angularintegral happens to be analytically doable [26]. For the ones proportional to δ ( ω − E p ) one can constrain the | p | integral and only two integrals are leftbut most of them have to be performed numerically. The terms proportional to δ ( ω ± ∆ ±± ) require more work as the domains where the ω = ± ∆ ±± constraintscan be satisfied are non-trivial. C.1 δ ( ω ) terms All the terms proportional to δ ( ω ) in the master integrals can be combinedaccording to equation (14). After performing an analogous work as in ref. [26],i.e. integrating by parts and performing the angular integrals we get: ρ tNLO N c C F g = ωδ ( ω )4 π Z ∞ dp n ′ F ( E p ) Z ∞ dk " ak n B ( k ) (cid:18) − p E p (cid:19) (51)+ n F ( E k ) E k ak − M − a k p E p − M k p E p E k − M kp E k E p (cid:26) M E k ( a − p − a + 1) E k + M (cid:27) ln (cid:12)(cid:12)(cid:12)(cid:12) p + kp − k (cid:12)(cid:12)(cid:12)(cid:12)! , where the parameter a keeps track of the thermal mass shift. Setting a = 0 isequivalent to performing the thermal mass shift (21), otherwise a = 1. C.2 δ ( ω − E p ) terms Summing all terms proportional to δ ( ω − E p ) occurring in equation (14) weget a formula of the form ρ fNLO ( ω ) = Z p Z k g ( k, E p , E pk , ω ) kE p E pk δ ( ω − E p ) . (52)While the function g can be easily constructed from the formulas of the abovesection, it is very long and will not be given here. Note that this integral containsdivergent terms and requires a careful renormalization. This proceeds as at zerotemperature [17] and will not be explained here. The thermal part we calculatehere is finite. We can rewrite the previous integral as: ρ fNLO ( ω ) = 18 π Z ∞ dk Z ∞ M dE p Z E + pk E − pk dE pk g ( k, E p , E pk , ω )2 δ ( ω − E p )= 18 π Z ∞ dk Z E + pk E − pk dE pk g (cid:16) k, ω , E pk , ω (cid:17) θ ( ω − M ) , (53)where E ± pk = q E p + k ± pk → q ω / k ± k √ ω − M after applyingthe delta function. The last two integrals are performed numerically and ρ fNLO ( ω ) = θ ( ω − M ) (cid:16) − n F (cid:16) ω (cid:17)(cid:17) (cid:16) ρ vac,fNLO + ρ bos,fNLO (cid:17) + ρ ferm, NLO ( ω ) , (54)16here ρ ferm, NLO ( ω )4 N c C F g = (55)+ θ ( ω − M )8 π (cid:16) − n F (cid:16) ω (cid:17)(cid:17) Z ∞ dk Z E + pk E − pk dE pk n ′ F2 M (cid:0) M + ω (cid:1) E pk ( ω − M ) ω × " k + (cid:0) M + E pk ω (cid:1) (cid:20) δ ++ + 1 δ −− (cid:21) + (cid:0) M − E pk ω (cid:1) (cid:20) δ + − + 1 δ − + (cid:21) + n F2 " M + ω ωE pk ( ω − M ) ( ωM E pk (cid:18) δ − δ − − δ − δ − + − δ + δ + δ + δ −− (cid:19) − M ( M − E pk (1 − a )) (cid:18) δ ++ + 1 δ + − + 1 δ − + + 1 δ −− (cid:19) + k ( aE pk − M ) ) + 12 ω ( − kωa + M (2 M + ω ) E pk (cid:18) δ + δ − + + δ − δ − δ − δ −− − δ + δ − (cid:19) ++ (cid:18) P kδ − + 1 kδ + (cid:19) (cid:0) k (cid:0) M + 3 ω (cid:1) + 8 M − ω (cid:1) + (cid:18) kδ −− + 1 kδ − + (cid:19) (cid:16) − k ω + (cid:0) kω + 2 M (cid:1) − ω ( k − ω ) (cid:17) + (cid:18) kδ + − + 1 kδ ++ (cid:19) (cid:16) k ω − (cid:0) M − kω (cid:1) + ω ( k + ω ) (cid:17) − ω (1 − a ) (cid:18) δ ++ + 1 δ + − + 1 δ − + + 1 δ −− (cid:19)!) with δ ±± = 2 k ± ω ± E pk , δ ± = ω ± E pk and again, setting a = 0 is equivalent toperforming the thermal mass shift (21). The symbol P means that we treat thepole at δ − = 0 in the principal value sense. Numerically this can be implementedas follows. We split the integral as R ω/ E − pk + R E + pk ω/ , perform a change of integrationvariable E pk → ω − E pk in the second integral and add it to the first one. C.3 δ ( ω ± ∆ ±± ) terms Summing all terms proportional to δ ( ω ± ∆ ±± ) occurring in equation (14) weget a formula of the form X ± , ± Z p Z k f ±± ( k, E p , E pk , ω ) kE p E pk ( δ ( ω − ∆ ±± ) − δ ( ω + ∆ ±± )) . (56)17f we restrict ourselves to ω >
0, only half of the δ ’s can actually be realized insome domain of the integrals so that the previous sum becomes θ ( ω − M )8 π Z k dk Z E p E − p dE p f ++ ( k, E p , ω − E p − k, ω ) − π Z ∞ k dk Z E p E − p dE p f −− ( k, E p , ω − E p + k, ω )+ 18 π Z ∞ k dk Z ∞ E p dE p f + − ( k, E p , − ω + E p + k, ω )+ 18 π Z ∞ k dk Z ∞− E − p dE p f − + ( k, E p , ω + E p − k, ω ) , (57)where the boundaries of the integrals are given by k = ω − M ω , k = max (cid:0) , − k (cid:1) , k = ω ,E ± p = ω − k ± k s − M ω ( ω − k ) ,E ± p = ω + k ± k s − M ω ( ω + 2 k ) (58)and the functions f ± , ± through f ++ ( k, E p , ω − E p − k, ω ) = π kM ω ( ω − E p ) (2( E p + k ) − ω ) (59)+ M (3 ω − E p ) ω ( ω − E p ) + 2( E p + k − ω ) (2 E p − ω )(2( E p + k ) − ω ) ! × [( n B0 + 1)(1 − n F ( ω − E p − k ) − n F1 ) + n F1 n F ( ω − E p − k )] ,f −− ( k, E p , ω − E p + k, ω ) = π − kM ω ( ω − E p ) (2( E p − k ) − ω ) (60)+ M (3 ω − E p ) ω ( ω − E p ) + 2( k − E p + ω ) (2 E p − ω )(2( E p − k ) − ω ) ! × [ n F1 n F ( ω − E p + k ) − n B0 (1 − n F ( ω − E p + k ) − n F1 )] ,f + − ( k, E p , E p + k − ω, ω ) = π kM ω ( ω − E p ) (2( E p + k ) − ω ) (61)+ M (3 ω − E p ) ω ( ω − E p ) + 2( E p + k − ω ) (2 E p − ω )(2( E p + k ) − ω ) ! × [ − ( n B0 + 1) n F ( E p + k − ω ) + n F1 n B0 + n F1 n F ( E p + k − ω )]18nd f − + ( k, E p , ω + E p − k, ω ) = π − kM ω (2 E p + ω ) (2( E p − k ) + ω ) (62)+ M (4 E p + 3 ω ) ω (2 E p + ω ) + 2( E p − k + ω ) (2 E p + ω )(2( E p − k ) + ω ) ! × [ n B0 n F ( E p − k + ω ) − n F1 ( n B0 + 1) + n F1 n F ( E p − k + ω )] . Note that at T = 0 only the first integral in (57) contributes. This term,setting T = 0, added to the expression ρ vac,fNLO in equation (54) gives the fullzero temperature result ρ vacNLO given in formula (24). For M ≫ T the firsttwo integrals contribute. If we take in these terms the part proportional to n B ( k ) and not containing any Fermi-Dirac distribution function and add themto the expression ρ bos,fNLO in equation (54), we get the ’bosonic’ thermal correction ρ bosNLO given in formula (25). The remaining terms are exponentially suppressedif M ≫ T but dominate the spectrum at small ω .Even if the full result is infrared finite, the different terms are not. Toavoid such problems one can first add the two last integrals in (57) after havingperformed the shift E p → E p − ω + k in the last integral. As a result of thatwe get the three last lines of formula (26) and then in all terms add an ǫ to thelower bound of the k integration and take the limit ǫ → C.4 Fermionic contribution
The full result is the sum of the vacuum part (24), the ’bosonic’ thermal cor-rections (25) and the fermionic contribution, which we can write as ρ fermNLO = ρ tNLO − n F (cid:16) ω (cid:17) ( ρ vacNLO + ρ bosNLO ) + ρ ferm, NLO + ρ ferm, NLO + ρ ferm, NLO , (63)where ρ tNLO is given in formula (51), ρ vacNLO in (24), ρ bosNLO in (25), ρ ferm, NLO in(27), ρ ferm, NLO in (55) and the remaining ρ frem, NLO below: ρ frem, NLO N c C F g = θ ( ω − M )8 π Z k dk Z E p E − p dE p kM ω ( ω − E p ) (2( E p + k ) − ω )+ M (3 ω − E p ) ω ( ω − E p ) + 2( E p + k − ω ) (2 E p − ω )(2( E p + k ) − ω ) ! (64) × [( n B0 + 1)(2 n F ( ω/ − n F ( ω − E p − k ) − n F1 ) + n F1 n F ( ω − E p − k )] − π Z ∞ k dk Z E p E − p dE p − kM ω ( ω − E p ) (2( E p − k ) − ω )+ M (3 ω − E p ) ω ( ω − E p ) + 2( k − E p + ω ) (2 E p − ω )(2( E p − k ) − ω ) ! × [ n F1 n F ( ω − E p + k ) − n B0 (2 n F ( ω/ − n F ( ω − E p + k ) − n F1 )] . Mass shift
The results of ref. [26] for the Euclidean correlator can be modified to accountfor the mass shift (19). We have to add to G V ( τ )4 N c C F g in equations (4.4-4.5) of [26]the following integral G MSV ( τ )4 N c C F g = Z p δM T " p E p (cid:0) D E p ( τ ) + 2 T n ′ F ( E p ) (cid:1) (65)+ (cid:18) M E p (cid:19) ∂ E p D E p ( τ )2 E p + M E p n ′′ F ( E p ) . In formula (25) the thermal mass shift has been performed. Without mass shift, ρ bosNLO N c C F g would contain an additional θ ( ω − M )(4 π ) ω Z ∞ dk n B ( k ) k (cid:18) ωk p ω − M − k ω (2 M + ω ) √ ω − M (cid:19) . (66) References [1] L. D. McLerran and T. Toimela, “Photon and Dilepton Emission from theQuark - Gluon Plasma: Some General Considerations,” Phys. Rev. D (1985) 545.[2] H. A. Weldon, “Reformulation of finite temperature dilepton production,”Phys. Rev. D (1990) 2384.[3] E. V. Shuryak, “Quark-Gluon Plasma and Hadronic Production of Leptons,Photons and Psions,” Phys. Lett. B , 150 (1978) [Sov. J. Nucl. Phys. ,408 (1978)] [Yad. Fiz. , 796 (1978)].[4] T. Matsui and H. Satz, “ J/ψ
Suppression by Quark-Gluon Plasma Forma-tion,” Phys. Lett. B (1986) 416.[5] B. B. Abelev et al. [ALICE Collaboration], “Azimuthal anisotropy of D me-son production in Pb-Pb collisions at √ s NN = 2 .
76 TeV,” arXiv:1405.2001[nucl-ex].[6] L. Adamczyk et al. [STAR Collaboration], “Observation of D me-son nuclear modifications in Au+Au collisions at √ s NN = 200 GeV,”arXiv:1404.6185 [nucl-ex].[7] J. I. Kapusta and C. Gale, Finite-Temperature Field Theory: Principlesand Applications ,Cambridge University Press, Cambridge, 2006[8] G. Aarts, J. M. Martinez Resco, “Continuum and lattice meson spectralfunctions at nonzero momentum and high temperature,” Nucl. Phys.
B726 (2005) 93-108. [hep-lat/0507004].[9] M. Laine, “NLO thermal dilepton rate at non-zero momentum,” JHEP (2013) 120 [arXiv:1310.0164 [hep-ph]].2010] I. Ghisoiu and M. Laine, “Interpolation of hard and soft dilepton rates,”arXiv:1407.7955 [hep-ph].[11] P. A. Baikov, K. G. Chetyrkin and J. H. Kuhn, “R(s) and hadronic tau-Decays in Order alpha**4(s): Technical aspects,” Nucl. Phys. Proc. Suppl. (2009) 49 [arXiv:0906.2987 [hep-ph]].[12] A. H. Hoang, V. Mateu and S. Mohammad Zebarjad, “Heavy Quark Vac-uum Polarization Function at O(alpha**2(s)) O(alpha**3(s)),” Nucl. Phys.B (2009) 349 [arXiv:0807.4173 [hep-ph]].[13] A. O. G. Kallen and A. Sabry, “Fourth order vacuum polarization,” Kong.Dan. Vid. Sel. Mat. Fys. Med. (1955) 17, 1.[14] R. Barbieri and E. Remiddi, “Infra-red divergences and adiabatic switching.fourth order vacuum polarization,” Nuovo Cim. A (1973) 99.[15] D. J. Broadhurst, J. Fleischer and O. V. Tarasov, “Two loop two pointfunctions with masses: Asymptotic expansions and Taylor series, in anydimension,” Z. Phys. C (1993) 287 [hep-ph/9304303].[16] T. Altherr and P. Aurenche, “Finite Temperature QCD Corrections toLepton Pair Formation in a Quark - Gluon Plasma,” Z. Phys. C (1989)99.[17] Y. Burnier, M. Laine, M. Vepsalainen, “Heavy quark medium polarizationat next-to-leading order,” JHEP (2009) 008. [arXiv:0812.2105 [hep-ph]].[18] E. Braaten, R. D. Pisarski and T. -C. Yuan, “Production of Soft Dileptonsin the Quark - Gluon Plasma,” Phys. Rev. Lett. (1990) 2242.[19] M. Asakawa, T. Hatsuda and Y. Nakahara, “Maximum entropy analysis ofthe spectral functions in lattice QCD,” Prog. Part. Nucl. Phys. (2001)459 [hep-lat/0011040].[20] A. Rothkopf, “Improved Maximum Entropy Analysis with an ExtendedSearch Space,” J. Comput. Phys. (2013) 106 [arXiv:1110.6285[physics.comp-ph]].[21] Y. Burnier and A. Rothkopf, “Bayesian Approach to Spectral FunctionReconstruction for Euclidean Quantum Field Theories,” Phys. Rev. Lett. (2013) 18, 182003 [arXiv:1307.6106 [hep-lat]].[22] Y. Burnier, M. Laine, L. Mether, “A Test on analytic continuationof thermal imaginary-time data,” Eur. Phys. J. C71 (2011) 1619.[arXiv:1101.5534 [hep-lat]].[23] G. Cuniberti, E. De Micheli and G.A. Viano, ” Reconstructing the thermalGreen functions at real times from those at imaginary times,” Commun.Math. Phys. 216 (2001) 59 [cond-mat/0109175].[24] Y. Burnier and M. Laine, “Towards flavour diffusion coefficient and elec-trical conductivity without ultraviolet contamination,” Eur. Phys. J. C (2012) 1902 [arXiv:1201.1994 [hep-lat]].2125] H. T. Ding, A. Francis, O. Kaczmarek, F. Karsch, H. Satz and W. Soeldner,“Charmonium properties in hot quenched lattice QCD,” Phys. Rev. D (2012) 014509 [arXiv:1204.4945 [hep-lat]].[26] Y. Burnier and M. Laine, “Massive vector current correlator in thermalQCD,” JHEP (2012) 086 [arXiv:1210.1064 [hep-ph]].[27] M. Le Bellac, ”Thermal Field Theory,” Cambridge University Press, Cam-bridge, 2000[28] R. E. Norton and J. M. Cornwall, “On the Formalism of Relativistic ManyBody Theory,” Annals Phys. (1975) 106.[29] Y. Burnier, M. Laine, M. Vepsalainen, “Heavy quarkonium in any channelin resummed hot QCD,” JHEP (2008) 043. [arXiv:0711.1743 [hep-ph]].[30] J. Casalderrey-Solana and D. Teaney, Phys. Rev. D (2006) 085012[hep-ph/0605199].[31] S. Caron-Huot, M. Laine and G. D. Moore, JHEP (2009) 053[arXiv:0901.1195 [hep-lat]].[32] G. D. Moore and D. Teaney, “How much do heavy quarks thermalize in aheavy ion collision?,” Phys. Rev. C (2005) 064904 [hep-ph/0412346].[33] S. Caron-Huot and G. D. Moore, “Heavy quark diffusion in perturba-tive QCD at next-to-leading order,” Phys. Rev. Lett. (2008) 052301[arXiv:0708.4232 [hep-ph]].[34] Y. Burnier and M. Laine, “Charm mass effects in bulk channel correla-tions,” JHEP (2013) 012 [arXiv:1309.1573 [hep-ph]].[35] D. Seibert, “The high frequency finite temperature quark dispersion rela-tion,” [nucl-th/9310008].[36] H. -T. Ding, A. Francis, O. Kaczmarek, H. Satz, F. Karsch, W. Soldner,“Charmonium correlation and spectral functions at finite temperature,”PoS LATTICE2010 (2010) 180. [arXiv:1011.0695 [hep-lat]].[37] G. D. Moore and J. M. Robert, hep-ph/0607172.[38] S. Caron-Huot, “Asymptotics of thermal spectral functions,” Phys. Rev. D (2009) 125009 [arXiv:0903.3958 [hep-ph]].[39] P. B. Arnold, G. D. Moore and L. G. Yaffe, “Transport coefficients in hightemperature gauge theories. 1. Leading log results,” JHEP (2000)001 [hep-ph/0010177].[40] P. B. Arnold, G. D. Moore and L. G. Yaffe, “Transport coefficients in hightemperature gauge theories. 2. Beyond leading log,” JHEP (2003)051 [hep-ph/0302165].[41] H.-T. Ding, A. Francis, O. Kaczmarek, F. Karsch, E. Laermann andW. Soeldner, “Thermal dilepton rate and electrical conductivity: An anal-ysis of vector current correlation functions in quenched lattice QCD,” Phys.Rev. D83