Quarkonium spectral function in medium at next-to-leading order for any quark mass

The vector channel spectral function at zero spatial momentum is calculated at next-to-leading order in thermal QCD for any quark mass. It corresponds to the imaginary part of the massive quark contribution to the photon polarization tensor. The spectrum shows a well defined transport peak in contrast to both the heavy quark limit studied previously, where the low frequency domain is exponentially suppressed at this order and the naive massless case where it vanishes at leading order and diverges at next-to-leading order. From our general expressions, the massless limit can be taken and we show that no divergences occur if done carefully. Finally, we compare the massless limit to results from lattice simulations.


Introduction
In heavy ion collisions, fireballs of quark-gluon plasma are formed. The QCD matter being strongly interacting, the quarks and gluons only escape as mesons or hadrons when the plasma cooled down sufficiently. To reconstruct what happened at early stages of the collision, we have to resort to probes that can be traced in experiment and whose properties are modified in the plasma. For several available probes, one quantity describes their fate in the plasma: the vector channel spectral function in medium.
For instance the way bound states such as charmonium or bottomonium decay 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 describes how the massive quarks will be diffused or slowed down by the plasma. This is another observable which is of interest for experiment, in fact the heavy quarks can be tagged and their transverse momentum distribution studied (see for instance [5,6]. For light quarks, the zero frequency limit of the spectral function describes the electric conductivity [7] and its momentum depencance the photon or dilepton emission rate [8,9,10]. In the vacuum, this spectral function is known at five loops [11] for massless fermions 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 massless result is known since a long time [16] and the case of large masses with respect to the temperature M ≫ T was calculated rather recently [17]. Here we extend the previous calculations to any mass and discuss the transport part of the spectrum, which is suppressed in the heavy quark limit and was not obtained in 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 corrections are of higher order [17] and will not be addressed here either.
Of course in QCD, the convergence of perturbation theory is slow due to the largeness of the coupling α s and moreover finite temperature gauge theories suffer from infrared problems so that the full infrared physics requires nonperturbative methods. Lattice computations contain the full physics but are performed in Euclidean time and do not have a direct access to the Minkowskian spectral function. After measuring the corresponding Euclidean correlator, an analytic continuation is needed to obtain the desired spectral function. In the case of discrete numerical data of finite precision the reconstruction of the spectrum is very hard to perform [19,20,21,22]. The challenge is even bigger here since the Euclidean correlator is not even continuous at zero Euclidean times and hence the full analytical continuation is ill defined [23]. That's where perturbation theory could again be of use, since the zero Euclidean time limit (or the corresponding large frequency limit of the spectral function) is weakly coupled and accessible to perturbation theory. This 'large' perturbative part (containing zero and possibly finite temperature contributions) could be subtracted from the 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 Euclidean corralator was computed recently together with its mass dependence in [26]. In this paper we complete our program and calculate the related spectral function i.e. perform the analytic continuation. This is not a trivial task even though the Euclidean correlator G(τ, p) of ref. [26] can be cross checked by convoluting the spectral function ρ(ω, p) with the finite temperature kernel K(τ, ω): After defining the observables in section 2, we discuss the calculation in sec 3. and refer to appendices for details. In section 4 we present our results for the spectrum, discuss the transport coefficients and derive the massless limit, which we compare to lattice results. Conclusions are given in section 5.

Correlators and spectral functions 2.1 Basic setup
We consider the vector current of a massive quark described by the operator 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 Following [17] we will start from the Euclidean correlator in frequency space (ω n = 2πnT denote the Matsubara frequencies), which can be calculated using conventional finite temperature Feynman rules [7,27,28]. The spectral function can be determined from the discontinuity of the Euclidean correlator along the imaginary axis: Note that for ω ∼ gT usual perturbation theory breaks down and would require to use Hard Thermal Loop (HTL) resummation, which will not be considered here. As was shown in [17], no infrared divergences occur so that HTL corrections are subleading for ω ≫ gT .

Possible applications
One observable defined by this spectral function is the production rate of muon pairs from the decay of massive quark pairs [1,2]. If we suppose that the quark pair is at rest we have in particular: 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 the threshold ω ∼ 2M , where the spectrum can be obtained more precisely with dedicated resummations [17,29]. When speaking of heavy flavor diffusion or electric conductivity, the diffusion coefficient D is obtained from the low energy behavior of the spectrum. In fact one expects the spectral function to look like a Lorentzian in the low energy limit 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 spectrum deviates from a Lorentzian 1 . The diffusion coefficient can then be extracted as For a thorough discussion of this formula and the zero frequency limit see for instance ref. [7]. If the onset of non-transport physics ω UV is well separated from the transport peak, another strategy can be used [30]. The idea is to calculate another observable, the momentum diffusion coefficient κ, which is proportional to the drag coefficient η D . It can be extracted from the falloff of the Lorentzian peak [31]: where M kin is the in medium kinetic mass defined by the low momentum limit of the dispersion relation. 2 The fluctuation-dissipation theorem finally relates this coefficient to the diffusion coefficient D = 2T 2 /κ and hence to the drag The momentum diffusion coefficient can be calculated in perturbation theory with dedicated resummations [32]. However even if the resummations seem to catch the relevant physics, the convergence of the perturbative series for κ is at best very slow. In the case of the heavy quarks for instance, the first nonvanishing contribution arises at O(α 2 ) and the correction O(α 2 g) is an order of magnitude larger for typical heavy ion plasmas [33].

Outline of the calculation
We now turn to the calculation of the spectral function, following refs. [17,26,29,34].

Leading order and notations
Performing the Wick contractions and the trace algebra, we get at leading order (LO): where Q = (−iω±ǫ, 0) will be set in the process of taking the discontinuity. Note that 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-leading order (NLO), we introduce the following notations: with ∆(P ) = P 2 − M 2 , so that at leading order, All the relevant I's are calculated in appendix B and in particular: where we introduced the notation E 2

Next to leading order
After performing the Wick contractions and the trace algebra, we get the NLO in terms of master sum-integrals defined in equations (11): Note that the first four terms are independent of ω and do not contribute to the spectral function.

Renormalization
The previous NLO expression (14) is UV divergent but is finite after redefinition of the mass. The counterterms for the currents read: and using the pole mass scheme as in ref. [26], we set where we denoted E pk = E p+k , k = |k| and ∆ ±± = k ± E p ± E pk .

Thermal correction to the mass
After subtraction of the counterterms, the spectral function is finite everywhere but at the threshold ω = 2M . The divergence there comes from thermal corrections, in fact one can rewrite the whole spectral function as and the term δM 2 is actually responsible for the divergence. In fact we can resum this contribution by redefining the mass M 2 → M 2 + δM 2 T : The explicit shift δM 2 T , is the thermal contribution to the dispersion relation, which, for a massless fermion is the well-known and for a massive fermion, the less well known [35] expression which actually depends on the integration variable p. As a summary, to avoid divergences at the threshold, we resum the mass correction. To calculate the spectral function of a fermion of mass squared M 2 , we calculate the spectral function at the mass M 2 + δM 2 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]. The divergence is however integrable and the shift was left out in ref. [26], where the Euclidean correlator was calculated, but the changes are easily tractable (see appendix D).

Explicit calculation of the NLO result
While the leading order is given in equations (12,13) the next-to-leading order requires significantly more work. The full expression is obtained adding to the NLO (14), the mass counterterm (15) and the contribution from the thermal mass shift (19). The first step is to carry out the sums (see Appendix A) and take the discontinuity (5). We are then left with the integrals and a delta function remaining from the discontinuity. Parts of the integrals can be performed analytically and the remaining ones have to be done numerically. The explicit expressions for the master integrals are given in appendix B and the details on their integration in appendix C.

Results
The final result can be split into three parts First the vacuum part [13,14,15] that can be integrated explicitly Secondly the first thermal part, that we will denote 'bosonic' thermal correction, calculated in ref. [17] for which one integral is left for numerical evaluation (for mass shift contribution see Appendix D). It is proportional to the Bose-Einstein distribution function n B (k) and does not contain any Fermi-Dirac distribution: The remaining contribution ρ f erm N LO 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 is the only piece containing a non-vanishing transport peak near ω → 0 so that it dominates the spectrum at low frequency. In this part, two integrals have to be performed numerically. The full expression is rather lengthy and can be read from appendix C. For ω < 2M only a few structures actually contribute and we can quote the result here: where This last term ρ f erm,1 N LO contribute to the transport peak and comes form gluon emission or absorption by the massive fermion.
The full result containing the LO and NLO contributions is plotted in fig. 1 (left) for T = 1.5T c and different masses ranging from M = 0.5T to M = 11.9T corresponding to a bottom quark of mass 4.65GeV. Note that the spectral function is in fact negative and for clarity we show −ρ V (ω) in the plots. In fig. 1 (right) we show the zero frequency limit of the spectral function lim ω→0 ρ V (ω)/ω, which enters the determination of the transport coefficient D. The transport coefficient itself requires division by the suscptibility χ, for which we use the NLO result of ref. [26], see fig. 2. Separate results for the three contributions to the next to leading order are shown in fig. 2 for the case M = 3.3T , representing the charm quark at T = 1.5T c studied on the lattice in [36] and for the generic case M = T in fig. 3. In these figures we scaled the result to ω 2 so that the vacuum contribution goes to a constant at large frequency.

Charm transport
In the insert of the same figures 2, 3, we scale the spectrum to the frequency and zoom on the low frequency region to see the transport peak. We see that ρ(ω) ω goes smoothly to a constant with vanishing slope at ω → 0 and typical Lorentzian curvature, hence defining a diffusion coefficient D according to In the case of the charm quark at T = 1.5T c , we see that it is in principle small 2πT D ≈ 0.1 in comparison to the lattice results suggesting 2πT D ≈ 2 [36] and to the perturbative heavy quark limit 2πT D ≈ 10 [33].
That said, the transport peak is not well separated from the UV physics -at least in this low order calculation -and there is no region where η D ≪ ω ≪ ω UV so that the momentum diffusion coefficient cannot be defined straight from formula (9). Of course for the only purpose of defining κ or η D one could just fit the low frequency part with equation (7) and get η D as a fit parameter adjusted to the curvature of ρ(ω) ω around ω = 0. In the previous case we get η D ∼ 2.5 which translates, if one would still trust the fluctuation-dissipation theorem and supposing that M kin = M , into 2πT D ∼ 0.8. This differ from the direct estimate made above (2πT D ≈ 0.1) showing that the fluctuationdissipation theorem does not seem to apply here. Note again that the present computation is not systematic [37] for ω < gT so that no strong conclusion should be made with this remark.

Massless limit and electric conductivity
Let's we consider a fermion that is massless in the vacuum. Even if in this case HTL correction would be of the same order as our NLO result, it is interesting to see how our result behaves. As we resummed the thermal mass correction according to equation (20), we in fact have to calculate the spectral function of a fermion of mass δM 2 T . In a typical quark-gluon plasma produced in heavy ion collision we have numerically (21) δM 2 The spectrum of such a fermion is shown in fig. 3. We see that the spectrum has no divergence neither at the threshold nor at zero frequency and shows a transport peak. This result differs from the old result of ref. [16]. There, the mass shift was noticed and performed in the LO result but was not made in the NLO calculation where the mass was set to zero. Their NLO contribution hence diverges at zero frequency, whereas ours has a threshold structure at ω = 2δM T and a  [16] where the mass shift was not fully taken into account (dotted red line). Note that the threshold structure appears negative as we show only the thermal part of the NLO. For the complete spectrum see fig. 1 transport peak, see fig. 3. At high frequency both results agree and merge to the asymptotic behavior derived in ref. [38], which reads 3 The transport coefficient obtained here is of order 2πT D ∼ 0.3, which is again on the low side, the perturbative resummation from [39,40] gets in this case 2πT D ∼ 25 and lattice results ranges from 2πT D ∼ 1 − 6 depending on the analytic continuation method used [41,22].

Matching the massless limit with lattice results
The full spectral function for M 2 = δM 2 T = g 2 C F T /4 with T = 1.46T c is shown in 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] for a massless fermion. Apart from an overall normalization 4 we see a different slope at large τ signaling a lack of power in the small ω region. This could be compensated for by an additional transport peak. Keeping for instance the same width η D = 0.7, if we add an additional −ρ t /ω =  fig. 4b). We see that it now goes nicely parallel to the lattice data. The total diffusion coefficient would then be 2πT D ∼ 0.4, still on the low side. In the insert, we fit the low energy spectrum with a Lorentzian. (Right): Nextto-leading order Euclidean correlator for M = T (blue continuous line) together with lattice data (dots). We note that the agreement can be greatly improved adding more power to the transport peak, even if some normalization factor remains. (magenta dashed line)

Conclusion
We calculated the thermal correction to the massive quark vector spectral function at NLO in thermal QCD. The thermal corrections are small in comparison to 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 production is smoothed by thermal corrections, the discontinuity in the vacuum spectrum being partly compensated by thermal effects. Secondly a transport peak appears at this order, it is however small in comparison to other expectations. We also see that the transport peak is broad and is not well separated from other kind of physics, rendering its determination via the fluctuationdissipation theorem difficult. It should however be stressed that higher loop orders give corrections of the same order for ω ≤ gT .
Setting the vacuum fermion mass to zero in our formulas do not lead to divergences. This contradict the calculation of ref. [16] where the spectral function was calculated for massless fermions and the result for the NLO diverges at zero frequency hindering a definition of the corresponding Euclidean correlator. The difference can be traced back to how the thermal mass shift is introduced in the NLO. In the previous reference, the thermal mass shift was performed in the leading order part as done here but not in the NLO contribution where the fermion remained purely massless, leading to divergences at zero frequency.

A Calculation of the master integrals
The calculation of the master sum-integrals is mostly standard (for details see ref. [28]), only one difficulty arises in double poles at zero frequency, which will be explained below. Their calculation is subtle and is shown explicitly for the case of I 0 01111 , which is the simplest from the amount of algebra but contains most of the technical difficulties.
We start by rewriting the master integral in a form where the elementary summation formula (valid for 0 < τ < β) can be applied. To do that we shift P − K → K and introduced the new integration variables S, R and the corresponding delta functions: 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 ) The factors e iβ(pn+kn) = e iβ(pn+sn+qn) = 1 were added so that the total phase of all the Matsubara frequencies are between 0 and β. Thanks to the time ordering and this last prescription the sums can be performed with formula (29) and then the integrals over τ 1 , τ 2 calculated. Remembering that q n is a bosonic Matusubara frequency we can set e iqnβ → 1 5 . The result is a rather long expression containing many terms that can be simplified in rewriting all the exponents in terms of Fermi-Dirac (or Bose-Einstein) distributions. Each term contains a product of Fermi-Dirac (or Bose-Einstein) distributions in the numerator 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 there are 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 the denominator, the numerator vanishes as well. This cancellation obviously occurs since the initial integral (33) has no poles. Now, it is possible to extract the discontinuity in each of these terms. The difficulty arises in terms proportional to This term contains a simple pole at ω = 0 when we enforce one of the delta functions contained in (30). Note that its contribution is finite even after replacing the second delta function as the numerator f vanishes when energies are set equal. However this term also contains a double pole when both deltas are set to zero. The simplest way to deal with this issue is to rewrite the delta functions in (30) as and use the δ(E s − E k − E r + E p ) to replace E s − E k by E r − E p (note that this replacement has to be done as a limit). After this, all the poles which contribute at ω ∼ 0 are of the type (ω ± (E r − E p )). The discontinuity can be taken in each term separately using for simple poles and for double poles. The spatial integrals over r, s can be performed using the remaining delta functions δ 3 (r − p)δ 3 (s − k). They force us to perform the limit E r → E p and get as final result:

B Master integrals
Following the notations of ref. [17], we detail here the different master integrals k being the gluon momentum). First the leading order sum-integrals: From them a few NLO sum-integrals are derived: 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 ω > 0 can be read from [17].
The last sum-integrals are given for all ω for completeness (as they were not all written in a suitable form for our present purpose): and

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 it is enough to consider the case ω ≥ 0. Using rotation symmetry, it is obvious that 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 proportional to ωδ(ω) which we call ρ t N LO for transport, a part proportional to δ(ω − 2E p ), called ρ f N LO for factorized and a last part proportional to one of the δ(ω ± ∆ ±± ) denoted by ρ p N LO for phase space, following the notation of [17]. For the terms proportional to δ(ω) the three integrals have to be performed but the angular integral happens to be analytically doable [26]. For the ones proportional to δ(ω − 2E p ) one can constrain the |p| integral and only two integrals are left but most of them have to be performed numerically. The terms proportional to δ(ω ±∆ ±± ) require more work as the domains where the ω = ±∆ ±± constraints can be satisfied are non-trivial.

C.1 δ(ω) terms
All the terms proportional to δ(ω) in the master integrals can be combined according to equation (14). After performing an analogous work as in ref. [26], i.e. integrating by parts and performing the angular integrals we get: where the parameter a keeps track of the thermal mass shift. Setting a = 0 is equivalent to performing the thermal mass shift (21), otherwise a = 1.

C.2 δ(ω − 2E p ) terms
Summing all terms proportional to δ(ω − 2E p ) occurring in equation (14) we get a formula of the form While the function g can be easily constructed from the formulas of the above section, it is very long and will not be given here. Note that this integral contains divergent terms and requires a careful renormalization. This proceeds as at zero temperature [17] and will not be explained here. The thermal part we calculate here is finite. We can rewrite the previous integral as: where E ± pk = E 2 p + k 2 ± 2pk → ω 2 /4 + k 2 ± k √ ω 2 − 4M 2 after applying the delta function. The last two integrals are performed numerically and where with δ ±± = 2k±ω ±E pk , δ ± = ω ±2E pk and again, setting a = 0 is equivalent to performing the thermal mass shift (21). The symbol P means that we treat the pole at δ − = 0 in the principal value sense. Numerically this can be implemented as follows. We split the integral as ω/2 E − pk + E + pk ω/2 , perform a change of integration variable E pk → ω − E pk in the second integral and add it to the first one.
If we restrict ourselves to ω > 0, only half of the δ's can actually be realized in some domain of the integrals so that the previous sum becomes Note that at T = 0 only the first integral in (57) contributes. This term, setting T = 0, added to the expression ρ vac,f N LO in equation (54) gives the full zero temperature result ρ vac N LO given in formula (24). For M ≫ T the first two 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 them to the expression ρ bos,f N LO in equation (54), we get the 'bosonic' thermal correction ρ bos N LO given in formula (25). The remaining terms are exponentially suppressed if M ≫ T but dominate the spectrum at small ω.
Even if the full result is infrared finite, the different terms are not. To avoid such problems one can first add the two last integrals in (57) after having performed the shift E p → E p − ω + k in the last integral. As a result of that we get the three last lines of formula (26) and then in all terms add an ǫ to the lower bound of the k integration and take the limit ǫ → 0 after having added all terms.

C.4 Fermionic contribution
The full result is the sum of the vacuum part (24), the 'bosonic' thermal corrections (25) and the fermionic contribution, which we can write as below:

D Mass shift
The results of ref. [26] for the Euclidean correlator can be modified to account for the mass shift (19). We have to add to GV (τ ) 4NcCFg 2 in equations (4.4-4.5) of [26] the following integral In formula (25)