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 polarisation 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][5][6]. Another example is the heavy flavour diffusion coefficient, which describes how the massive quarks will be diffused or slowed down by the plasma. This observable is of interest for experiment as heavy quarks can be tagged and their transverse momentum distribution studied; see for instance [7,8]. For light quarks, the zero frequency limit of the spectral function describes the a e-mail: yannis.burnier@epfl.ch electric conductivity [9] and its momentum dependence the photon or dilepton emission rate [10][11][12].
In the vacuum, this spectral function is known at five loops [13] for massless fermions and Taylor expansions in the mass are known to four loops [14][15][16][17]. In the presence of a finite temperature medium, the two loop massless result has been known since a long time [18] and the case of large masses with respect to the temperature M T was calculated rather recently [19]. 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 [20] in the fermion propagators and vertices. Other HTL corrections are of higher order [19] 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 [21][22][23][24]. 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 [25]. 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 [26] or used as a prior to define the analytical continuation [27].
In the case of the vector current spectral function considered here, the Euclidean correlator was computed recently together with its mass dependence in [28]. 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. [28] can be cross checked by convoluting the spectral function ρ(ω, p) with the finite temperature kernel K (τ, ω): After defining the observables in Sect. 2, we discuss the calculation in Sect. 3. and refer to appendices for details. In Sect. 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 Sect. 5.

Basic setup
We consider the vector current of a massive quark described by the operator (τ, x) 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 [19] 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 [9,29,30]. 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 us to use hard thermal loop (HTL) resummation, which will not be considered here. As was shown in [19], 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 [19,31]. When speaking of heavy flavour diffusion or electric conductivity, the diffusion coefficient D is obtained from the low energy behaviour 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 ω U V 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: For a thorough discussion of this formula and the zero frequency limit see for instance Ref. [9]. If the onset of non-transport physics ω U V is well separated from the transport peak, another strategy can be used [32]. 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 [33]: where M kin is the in-medium kinetic mass defined by the low momentum limit of the dispersion relation. 2 The fluctuationdissipation theorem finally relates this coefficient to the diffusion coefficient D = 2T 2 /κ and hence to the drag coefficient The momentum diffusion coefficient can be calculated in perturbation theory with dedicated resummations [36]. 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 non-vanishing contribution arises at O(α 2 ) and the correction O(α 2 g) is an order of magnitude larger for typical heavy ion plasmas [37].
Here we head for another approach, namely a lattice determination of the low frequency spectral features and we aim at calculating the perturbative large frequency part, which is needed to properly analyse the Euclidean lattice data and extract the low frequency part [26].

Outline of the calculation
We now turn to the calculation of the spectral function, following Refs. [19,28,31,38].

Leading order and notations
At leading order the only diagram is where the squares denotes the current insertion and the lines the quark propagator. Performing the Wick contractions and the trace algebra, we get at leading order (LO) (10) where and denotes the sum integrals over bosonic and fermionic matsubara frequencies. The spectrum is then given by where Q = (−iω ± , 0) will be set in the process of taking the discontinuity. Note that the first term in Eq. (10) is 2 Namely the velocity dependence of the free energy is expanded as independent of Q and does 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 are calculated in Appendix B and in particular: where we introduced the notation E 2

Next-to-leading order
At leading order the diagrams are: where the curly line denotes the gluon propagator. After performing the Wick contractions and the trace algebra, we get the NLO in terms of the master sum integrals defined in Eqs. : Note that the first four terms are independent of ω and do not contribute to the spectral function.

Renormalisation
The previous NLO expression (16) is UV divergent but is finite after redefinition of the mass. The counterterms for the currents read (18) and using the pole mass scheme as in Ref. [28], we set where we denoted E pk = E p+k , k = |k| and ±± = k ± E p ± E pk and of course the last integral is independent of p.

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 wherē is actually responsible for the divergence. In fact we can resum this contribution by redefining the mass The explicit shift δ M 2 T , is the thermal contribution to the dispersion relation, which, for a massless fermion is the wellknown expression and for a massive fermion, the less well-known [39] expression which actually depends on the integration variable p. In the heavy quark limit M T , the expression simplifies again [34,35] and we get 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 Eq. (22) and modify the NLO contribution as explained in Eq. (21). Note that the relevant part of this shift (25) was performed in [19]. The divergence is, however, integrable and the shift was left out in Ref. [28], 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 Eqs. (13,15) the nextto-leading order requires significantly more work. The full expression is obtained adding to the NLO (16), the mass counterterm (17) and the contribution from the thermal mass shift (21). 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 [15][16][17], which can be integrated explicitly, where Second, the first thermal part, which we will denote the 'bosonic' thermal correction, calculated in Ref. [19] for which one integral is left for numerical evaluation (for the 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: where ± = ω 4 − 4M 4 ± 2ωk(ω 2 + 2M 2 ) + 2ω 2 k 2 . The remaining contribution ρ ferm NLO 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 this piece contains a nonvanishing transport peak near ω → 0 and 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, ρ ferm,1 NLO , contributes to the transport peak and comes from 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 (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 susceptibility χ , for which we use the NLO result of Ref. [28]; see  Separate results for the three contributions to the next-toleading 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 [40] 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
Keeping in mind that the present computation is not systematic 3 [41] for ω < gT , we shall still briefly discuss its low frequency limit. In the insert of the same figures, Figs. 2 and 3, we scale the spectrum to the frequency and zoom in 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 Eq. (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 ) −2 and behaves as a power law T D ∝ (M/T ) − 1 4 at small M.
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 [40] and in comparison to the perturbative heavy quark limit, 2π T D ≈ 10 [37]. 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 ω ω U V so that the momentum diffusion coefficient cannot be defined straight from Eq. (9).

Massless limit and electric conductivity
Let us consider a fermion that is massless in the vacuum. Even if at low frequency, the 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 Eq. (22), 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 (23) 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. [18]. 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 transport peak; see Fig. 3. At high frequency both results agree and merge to the asymptotic behaviour derived in Ref. [42], which reads 4 The transport coefficient obtained here is of order 2π T D ∼ 0.3, which is again on the low side, the perturbative resummation from [43,44] 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 [24,[45][46][47].

Matching the massless limit with lattice results
The full spectral function for Fig. 4 together with the Euclidean correlator scaled to the free correlator. Keeping in mind that our approximations are not consistent at low frequencies, we still compare the Euclidean correlator to the continuum extrapolated quenched lattice data of Ref. [45] for a massless fermion.
In Fig. 4 right, we see first a quite good match at small τ but then a different slope at large τ , signalling a lack of power in the small ω region. To picture how far we are from the lattice data, we multiply the ω < 2M region by a constant so that to fit the lattice data. We get an excellent match with a factor 2.26. To see how the width of the transport peak shows up in the correlator we compare this first limit to an infinitely thin transport peak fitted to reproduce the lattice data (here the existing ω < 2M part has been removed). We see that the difference is sizeable but still much smaller than the given lattice error bars.

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. Second, a transport peak appears at this order, it is, however, small in comparison to other expectations. We also see that the transport peak is  , π T )). Note that the agreement can be improved adding more power to the transport peak. The red dashed curve is obtained increase the height of the existing transport peak by a factor 2.26 and the black dotted one contains instead of a transport peak an infinitely thin delta peak. Note that the precision of the lattice data has to be very high to distinguish between the two extreme cases broad and is not well separated from other kind of physics. 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 formulae do not lead to divergences. This contradict the calculation of Ref. [18] 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.
Acknowledgments The author would like to thank M. Laine for helpful discussions and important comments on the manuscript. This work was supported by the SNF Grant PZ00P2-142524.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Appendix A: Calculation of the master integrals
The calculation of the master sum integrals is mostly standard (for details, see Ref. [30]); 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 case as regards 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 introduce 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 ) = e iβ( p n +s n ) The factors e iβ( p n +k n ) = e iβ( p n +s n +q n ) = 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 Eq. (A.1) and then the integrals over τ 1 , τ 2 calculated. Remembering that q n is a bosonic Matusubara frequency we can set e iq n β → 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 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 (A.3) 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 (A.2). 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 (A.2) 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 the final result

Appendix B: Master integrals
Following the notations of Ref. [19], we detail here the different master integrals (denoting E 2 k being the gluon momentum). First the leading order sum integrals: , (B.10) From them a few NLO sum integrals are derived: (B.14) 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 [19], 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 8E p E pk k 2 × n B0 n F2 − (1 + n B0 )n F1 + n F1 n F2 8E p E pk k 2 where E ± pk = E 2 p + k 2 ± 2 pk is replaced by ω 2 /4 + k 2 ± k √ ω 2 − 4M 2 after applying the delta function. The last two integrals are performed numerically and with δ ±± = 2k ±ω± E pk , δ ± = ω±2E pk and, again, setting a = 0 is equivalent to performing the thermal mass shift (23). 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. get the 'bosonic' thermal correction ρ bos NLO given in Eq. (28). 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 (C.29) 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 Eq. (29) 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.

Fermionic contribution
The full result is the sum of the vacuum part (27), the 'bosonic' thermal corrections (28) and the fermionic contribution, which we can write as