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

Abstract

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.

1 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–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 electric conductivity [9] and its momentum dependence the photon or dilepton emission rate [10–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–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\gg 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 $\omega \gg 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 $\alpha _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–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(\tau ,\mathbf{p})$ of Ref. [28] can be cross checked by convoluting the spectral function $\rho (\omega ,\mathbf{p})$ with the finite temperature kernel $K(\tau ,\omega )$:

$$\begin{aligned}&G(\tau ,\mathbf{p})=\int _{0}^{\infty } \frac{\mathrm{d}\omega }{ \pi }\;\rho (\omega ,\mathbf{p})K(\tau ,\omega ),\nonumber \\&K(\tau ,\omega )=\frac{\cosh (\omega (\tau -\beta /2))}{\sinh (\omega \beta /2)}.\nonumber \end{aligned}$$

(1)

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.

2 Correlators and spectral functions

2.1 Basic setup

We consider the vector current of a massive quark described by the operator $\Psi (\tau ,\mathbf{x})$

$$\begin{aligned} J_\mu (\tau ,\mathbf{x})=\bar{\Psi }(\tau ,\mathbf{x})\gamma _\mu \Psi (\tau ,\mathbf{x}), \end{aligned}$$

(2)

with $\mu =0,\ldots ,d$. The object we compute here is the spectral function in medium, which is given by the thermal average of the current commutator

$$\begin{aligned} \rho ^V(\omega )=\int \mathrm{d}t\; e^{i\omega t}\int \mathrm{d}^{d}x \left\langle \frac{1}{2} [J^\mu (t,x),J_\mu (0,0)]\right\rangle _T. \end{aligned}$$

(3)

Following [19] we will start from the Euclidean correlator in frequency space ($\omega _n=2\pi n T$ denote the Matsubara frequencies),

$$\begin{aligned} G_E(\omega _n)=\int _0^\beta \mathrm{d}\tau e^{i\omega _n \tau }\int \mathrm{d}^{d}x\langle J^\mu (\tau ,x)J_\mu (0,0)\rangle , \end{aligned}$$

(4)

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:

$$\begin{aligned} \rho ^V(\omega )= & {} \mathop {\text{ Disc } }[G_E(-i\omega )]\nonumber \\= & {} \frac{1}{2i}\lim _{{\varepsilon }\rightarrow 0^+}[G_E(-i\omega +{\varepsilon })-G_E(-i\omega -{\varepsilon })]. \end{aligned}$$

(5)

Note that for $\omega \sim 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 $\omega \gg gT$.

2.2 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:

$$\begin{aligned} \frac{\mathrm{d}N_{\mu \bar{\mu }}}{\mathrm{d}^4x \mathrm{d}^4q}= & {} \frac{-2e^4Z^2}{3(2\pi )^5 \omega ^2}\left( 1+\frac{2m^2_\mu }{\omega ^2}\right) \left( 1-\frac{4m^2_\mu }{\omega ^2}\right) ^{\frac{1}{2}}\nonumber \\&\times n_\mathrm{B}(\omega )\rho ^V(\omega ), \end{aligned}$$

(6)

where Ze is the charge of the quark and $n_\mathrm{B}$ the Bose–Einstein distribution and $\omega =E_{\mu ^+}+E_{\mu ^-}$. The main contribution to this observable comes form the threshold $\omega \sim 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

$$\begin{aligned} -\frac{\rho ^V(\omega )}{\omega }\overset{0<\omega <\omega _{UV}}{\approx }3 \chi D\frac{\eta _D^2}{\eta _{D}^2+\omega ^2}, \end{aligned}$$

(7)

where $\chi$ is the susceptibility, $\eta _D$ another number called the drag coefficient and $\omega _{UV}$ the scale at which other kind of physics enter and where the spectrum deviates from a Lorentzian.¹ The diffusion coefficient can then be extracted:

$$\begin{aligned} D=-\frac{1}{3\chi }\lim _{\omega \rightarrow 0}\frac{\rho ^V(\omega )}{\omega }. \end{aligned}$$

(8)

For a thorough discussion of this formula and the zero frequency limit see for instance Ref. [9].

If the onset of non-transport physics $\omega _{UV}$ is well separated from the transport peak, another strategy can be used [32]. The idea is to calculate another observable, the momentum diffusion coefficient $\kappa$, which is proportional to the drag coefficient $\eta _D$. It can be extracted from the falloff of the Lorentzian peak [33]:

$$\begin{aligned} \kappa =2 M_\mathrm{kin} T \eta _D\approx \left. -\frac{2M_\mathrm{kin}^2\omega ^2}{3\chi }\frac{\rho ^V(\omega )}{\omega } \right| _{\eta _D\ll \omega \ll \omega _{UV}}, \end{aligned}$$

(9)

where $M_\mathrm{kin}$ is the in-medium kinetic mass defined by the low momentum limit of the dispersion relation.² The fluctuation–dissipation theorem finally relates this coefficient to the diffusion coefficient $D=2 T^2/\kappa$ and hence to the drag coefficient $\eta _D=\kappa /(2M_\mathrm{kin} T)$.

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 $\kappa$ is at best very slow. In the case of the heavy quarks for instance, the first non-vanishing contribution arises at $\mathcal {O}(\alpha ^2)$ and the correction $\mathcal {O}(\alpha ^2g)$ 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].

3 Outline of the calculation

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

3.1 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\omega \pm \epsilon ,\mathbf 0)$ will be set in the process of taking the discontinuity. Note that the first term in Eq. (10) is 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:

(11)

(12)

with $\Delta (P)=P^2-M^2$, so that at leading order,

$$\begin{aligned} \rho ^V_{LO}(\omega )=-2N_\mathrm{c}(4 M^2+2 \omega ^2)I_{11}(\omega ). \end{aligned}$$

(13)

All the relevant I are calculated in Appendix B and in particular:

$$\begin{aligned} I_{11}(\omega )= & {} \theta (\omega -2M)\frac{\sqrt{\omega ^2-4M^2}}{16\pi \omega }\tanh \left( \frac{\omega }{4T}\right) \end{aligned}$$

(14)

$$\begin{aligned}&\times \left[ 1+{\varepsilon }\left( 2+\ln \frac{\bar{\mu }^2}{\omega ^2-4 M^2}\right) \right] \nonumber \\&\quad +\pi \omega \delta (\omega )\int _p \frac{n_{\mathrm {F}}'(E_p)}{2E_p^2}, \end{aligned}$$

(15)

where we introduced the notation $E_p^2=\mathbf{p}^2+M^2$, $\int _p=\int \frac{\mathrm{d}^dp}{(2 \pi )^d}$ and set $d=3-2{\varepsilon }$.

3.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. (11):

$$\begin{aligned} \frac{\rho ^{V}_\mathrm{NLO}}{4 N_\mathrm{c}C_\mathrm{{F}}g^2}= & {} 4 (1-\epsilon )^2I_{12000}^0 -4 (1-\epsilon )I_{11100}^0\nonumber \\&+\,8 (1-\epsilon )M^2 I_{12100}^0-4 (1-\epsilon )^2 I_{02100}^0 \nonumber \\&+\,8(1-\epsilon ) I^0_{01110} \nonumber \\&-\,8(2 M^2+\omega ^2(1-\epsilon ))I^0_{11110}\nonumber \\&+\,4(1-\epsilon )I_{10110}^0 +8M^2 (2 M^2\nonumber \\&+\,\omega ^2 (1-\epsilon )) I_{12110}^0\nonumber \\&-\,4(1-\epsilon ) (2 M^2+\omega ^2 (1-\epsilon ))I_{02110}^0\nonumber \\&-\,8 (1-\epsilon )^2 I^1_{11110}\nonumber \\&-\,2(1-\epsilon )I_{-11111}^0+2 (2 M^2 \epsilon \nonumber \\&+\,\omega ^2 (2-\epsilon -\epsilon ^2))I_{01111}^0 \nonumber \\&+\,2 (4 M^4-2 M^2 \omega ^2 \epsilon -\omega ^4 (1-\epsilon ))I_{11111}^0\nonumber \\&-\,4(1-\epsilon )I_{11010}^0+ 4 (1-\epsilon ) (2 M^2\nonumber \\&+\,\omega ^2 (1-\epsilon ))I^0_{12010}. \end{aligned}$$

(16)

Note that the first four terms are independent of $\omega$ and do not contribute to the spectral function.

3.3 Renormalisation

The previous NLO expression (16) is UV divergent but is finite after redefinition of the mass. The counterterms for the currents read

$$\begin{aligned} \frac{\rho ^{V,CT}_\mathrm{NLO}}{4N_\mathrm{c}C_\mathrm{{F}}g^2}= \frac{\delta M^2}{g^2C_\mathrm{{F}}}\frac{ 1}{4N_\mathrm{c}}\frac{ \partial \rho ^V_{LO}}{\partial M^2} \end{aligned}$$

(17)

with

$$\begin{aligned} \frac{ 1}{4N_\mathrm{c}}\frac{ \partial \rho ^V_{LO}}{\partial M^2}=(4 M^2+2 \omega ^2 (1-\epsilon )) I_{21}-2I_{11}-2(1-\epsilon ) I_{20} \end{aligned}$$

(18)

and using the pole mass scheme as in Ref. [28], we set

$$\begin{aligned} \delta M^2= & {} -\frac{6g^2C_\mathrm{{F}}M^2}{(4\pi )^2}\left( \frac{1}{\epsilon }+\ln \frac{\bar{\mu }^2}{M^2}+\frac{4}{3} \right) \nonumber \\= & {} -g^2C_\mathrm{{F}}\int _{k} \frac{(2-2\epsilon ) (E_{pk}-k)+\frac{2 M^2}{\Delta _{-+}}+\frac{2 M^2}{\Delta _{++}}}{2 k E_{pk}}, \end{aligned}$$

(19)

where we denoted $E_{pk}=E_\mathbf{p+\mathbf k}$, $k=|\mathbf k|$ and $\Delta _{\pm \pm }=k\pm E_p\pm E_{pk}$ and of course the last integral is independent of p.

3.4 Thermal correction to the mass

After subtraction of the counterterms, the spectral function is finite everywhere but at the threshold $\omega =2M$. The divergence there comes from thermal corrections, in fact one can rewrite the whole spectral function as

$$\begin{aligned} \rho ^V(\omega ,M^2)= & {} \rho ^V_{LO}(\omega ,M^2)+\delta M_T^2 \frac{\partial \rho ^V_{LO}(\omega ,M^2)}{\partial M^2}\nonumber \\&+\bar{\rho }^V_\mathrm{NLO}(\omega ,M^2)+\mathcal {O}(g^4), \end{aligned}$$

(20)

where

$$\begin{aligned} \bar{\rho }^V_\mathrm{NLO}=\rho ^V_\mathrm{NLO}-\delta M_T^2 \frac{\partial \rho ^V_{LO}}{\partial M^2} \end{aligned}$$

(21)

and the term $\delta M_T^2 \frac{\partial \rho ^V_{LO}(\omega ,M^2)}{\partial M^2}$ is actually responsible for the divergence. In fact we can resum this contribution by redefining the mass $M^2\rightarrow M^2+\delta M^2_T$:

$$\begin{aligned} \rho ^V(\omega ,M^2)= & {} \rho ^V_{LO}(\omega ,M^2+\delta M_T^2)\nonumber \\&+\bar{\rho }^V_\mathrm{NLO}(\omega ,M^2+\delta M_T^2)+\mathcal {O}(g^4). \end{aligned}$$

(22)

The explicit shift $\delta M^2_T$, is the thermal contribution to the dispersion relation, which, for a massless fermion is the well-known expression

$$\begin{aligned} \delta M^2_T=g^2 C_\mathrm{{F}}\int _0^\infty \frac{\mathrm{d}k}{\pi ^2} k(n_{\mathrm {B}}(k)+n_{\mathrm {F}}(k))=\frac{g^2 C_\mathrm{{F}}T^2}{4} \end{aligned}$$

(23)

and for a massive fermion, the less well-known [39] expression

$$\begin{aligned} \frac{\delta M^2_T}{g^2C_\mathrm{{F}}}= & {} 2\int _k\left[ \frac{n_{\mathrm {B}}(k)}{k}\right. \nonumber \\&\left. +\frac{n_{\mathrm {F}}(E_{pk})}{ E_{pk}}\left( 1-\frac{M^2}{\Delta _{++}\Delta _{--}}-\frac{M^2}{\Delta _{+-}\Delta _{-+}}\right) \right] \nonumber \\ {}= & {} \int _0^\infty \frac{\mathrm{d}k}{\pi ^2} \bigg [k\,n_{\mathrm {B}}(k)\nonumber \\&+\frac{k^2n_{\mathrm {F}}(E_k)}{E_k} \left( 1+\frac{M^2}{2 p k} \ln \left| \frac{k-p}{k+p}\right| \right) \bigg ], \end{aligned}$$

(24)

which actually depends on the integration variable p. In the heavy quark limit $M\gg T$, the expression simplifies again [34, 35] and we get

$$\begin{aligned} \delta M^2_T=\frac{g^2 C_\mathrm{{F}}T^2}{6}. \end{aligned}$$

(25)

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+\delta 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).

3.5 Explicit calculation of the NLO result

While the leading order is given in Eqs. (13, 15) the next-to-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.

4 Results

The final result can be split into three parts,

$$\begin{aligned} \rho _\mathrm{NLO}^V(\omega )=\rho _\mathrm{NLO}^\mathrm{vac}(\omega )+\rho _\mathrm{NLO}^\mathrm{bos}(\omega )+\rho _\mathrm{NLO}^\mathrm{ferm}(\omega ). \end{aligned}$$

(26)

First the vacuum part [15–17], which can be integrated explicitly,

$$\begin{aligned} \frac{\rho ^\mathrm{vac}_\mathrm{NLO}(\omega )}{4N_\mathrm{c}C_\mathrm{{F}}g^2}= & {} \frac{ \theta (\omega -2M)}{(4\pi )^3 \omega } \nonumber \bigg [ \sqrt{\omega ^2-4 M^2} \left( -\frac{3(\omega ^2+6 M^2)}{4}\right. \nonumber \\&\left. +2(\omega ^2+2M^2) \ln \frac{\omega (\omega ^2-4 M^2)}{M^3} \right) \nonumber \\&+\frac{14 M^4+4 M^2\omega ^2-6\omega ^4}{\omega } a\mathrm {cosh}\left( \frac{\omega }{2 M}\right) \nonumber \\&+\frac{8 M^4-2\omega ^4}{\omega }L_2\left( \frac{\omega -\sqrt{\omega ^2-4 M^2}}{\omega +\sqrt{\omega ^2 - 4 M^2}}\right) \bigg ],\nonumber \\ \end{aligned}$$

(27)

where $L_2=4\mathrm {Li}_2(x) + 2\mathrm {Li}_2(-x) + \ln (x)[2\ln (1-x) +$ $\ln (1+x)].$

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_{\mathrm {B}}(k)$ and does not contain any Fermi–Dirac distribution:

$$\begin{aligned}&\frac{\rho ^\mathrm{bos}_\mathrm{NLO}}{4N_\mathrm{c}C_\mathrm{{F}}g^2 } = \frac{2}{(4\pi )^3 \omega ^2} \int _0^\infty \mathrm{d}k \frac{n_{\mathrm {B}}(k)}{k}\nonumber \\&\quad \times \biggl \{\theta (\omega ) \theta \left( k-\frac{4M^2-\omega ^2}{2\omega }\right) \nonumber \\&\quad \times \biggl [+(\omega ^2 + 2 M^2) \sqrt{\omega (\omega +2 k)}\nonumber \\&\quad \times \sqrt{\omega (\omega +2 k)- 4M^2} \nonumber \\&\quad - 2 \Omega _+ a\mathrm {cosh} \sqrt{\frac{\omega (\omega +2k)}{4M^2}}\nonumber \\&\quad +\, 2 \omega ^2 k^2 \sqrt{1-\frac{4 M^2}{\omega (\omega +2 k)}} \biggr ] +\theta (\omega - 2M)\theta \Bigl (\frac{\omega ^2-4M^2}{2\omega }-k\Bigr )\nonumber \\&\quad \times \biggl [+(\omega ^2 + 2 M^2)\sqrt{\omega (\omega -2 k)}\sqrt{\omega (\omega -2 k)- 4M^2} \nonumber \\&\quad - \,2 \Omega _- a\mathrm {cosh} \sqrt{\frac{\omega (\omega -2k)}{4M^2}} +2 \omega ^2 k^2 \sqrt{1-\frac{4 M^2}{\omega (\omega -2 k)}} \biggr ] \nonumber \\&\quad +\, \theta (\omega - 2M) \biggl [ - 2 (\omega ^2 + 2 M^2)\, \omega \sqrt{\omega ^2 - 4 M^2} \nonumber \\&\quad +\, 4 (\omega ^4 - 4 M^4 + 2 \omega ^2 k^2)a\mathrm {cosh} \biggl ( \frac{\omega }{2M} \biggr ) \biggr ] \biggr \}, \end{aligned}$$

(28)

where $\Omega _\pm =\omega ^4 - 4 M^4 \pm 2 \omega k (\omega ^2 + 2 M^2) + 2 \omega ^2 k^2$. The remaining contribution $\rho ^\mathrm{ferm}_\mathrm{NLO}$ contains at least one Fermi–Dirac distribution $n_{\mathrm {F}}(E_p)$ or $n_{\mathrm {F}}(E_{pk})$; hence it is suppressed in the limit $M\ll T$. However this piece contains a non-vanishing transport peak near $\omega \rightarrow 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 $\omega <2M$ only a few structures actually contribute and we can quote the result here:

$$\begin{aligned} \rho ^\mathrm{ferm}_\mathrm{NLO}\overset{\omega <2M}{=} -2n_{\mathrm {F}}\left( \frac{\omega }{2}\right) \rho ^\mathrm{bos}_\mathrm{NLO}+\rho ^{\mathrm{ferm},1}_\mathrm{NLO}, \end{aligned}$$

(29)

where

$$\begin{aligned} \frac{\rho ^{\mathrm{ferm},1}_\mathrm{NLO}}{4N_\mathrm{c}C_\mathrm{{F}}g^2}= & {} \frac{1}{8 \pi ^3}\int _{\omega /2}^\infty \int _{E_p^{1+}}^\infty \Biggl \{ 1+\frac{M^2 (2 M^2+\omega ^2)}{\omega ^2}\nonumber \\&\times \left( \frac{1}{(2 E_p+2 k-\omega )^2}+\frac{1}{(2 E_p-\omega )^2}\right) \nonumber \\&+\frac{-2 k^2 \omega ^2+(k \omega +2 M^2)^2-\omega ^2 (\omega -k)^2}{2 k \omega ^2}\nonumber \\&\times \left( \frac{1}{2 E_p+2 k-\omega }-\frac{1}{2 E_p-\omega }\right) \Biggr \} \nonumber \\&\times \{n_{\mathrm {B}}(k) (n_{\mathrm {F}}(E_p)-n_{\mathrm {F}}(E_p+k-\omega ))\nonumber \\&+(n_{\mathrm {F}}(E_p)-1) n_{\mathrm {F}}(E_p+k-\omega )\}. \end{aligned}$$

(30)

This last term, $\rho ^{\mathrm{ferm},1}_\mathrm{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.5 T_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 $-\rho ^V(\omega )$ in the plots. In Fig. 1 (right) we show the zero frequency limit of the spectral function $\lim _{\omega \rightarrow 0}\rho ^V(\omega )/\omega$, which enters the determination of the transport coefficient D. The transport coefficient itself requires division by the susceptibility $\chi$, for which we use the NLO result of Ref. [28]; see Fig. 2.

Fig. 1

Left LO (green dashed) and NLO (blue) spectral function ($-\rho ^V(\omega )/\omega ^2$) for $M=(0.5,1,3.3,5,10)T$, $T=1.5 T_c$. The LO vanishes at frequencies below the threshold $\omega =2M$ where the NLO shows a discontinuity. Right Value of the low energy limit of $-\rho ^V(\omega )/\omega$ as a function of the quark mass. In this plot we used $g=1.6$ as in [28]

Fig. 2

Left Different parts of the NLO contribution to the spectral function for $M=3.3T$ and Lorentzian fit of the transport peak (black dot-dashed line). Right The generated transport coefficient D

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.5 T_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 $\omega ^2$ so that the vacuum contribution goes to a constant at large frequency.

Fig. 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. [18] 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

4.1 Charm transport

Keeping in mind that the present computation is not systematic³ [41] for $\omega <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 $\frac{\rho (\omega )}{\omega }$ goes smoothly to a constant with vanishing slope at $\omega \rightarrow 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\propto (M/T)^{-2}$ and behaves as a power law $T D\propto (M/T)^{-\frac{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\pi T D\approx 0.1$, in comparison to the lattice results, suggesting $2\pi T D\approx 2$ [40] and in comparison to the perturbative heavy quark limit, $2\pi T D\approx 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 $\eta _D\ll \omega \ll \omega _{UV}$ so that the momentum diffusion coefficient cannot be defined straight from Eq. (9).

4.2 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 $\delta M_T^2$. In a typical quark–gluon plasma produced in heavy ion collision we have numerically (23) $\delta M_T^2=g^2 C_\mathrm{{F}}T^2/4 \approx T^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. [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 $\omega =2\delta 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⁴

$$\begin{aligned} -\rho ^T_\mathrm{NLO}= & {} -\rho ^{bos}_\mathrm{NLO}-\rho ^{\mathrm{ferm}}\nonumber \\ \overset{\omega \gg M,T}{}= & {} 4 N_\mathrm{c}C_\mathrm{{F}}g^2\left( \frac{\pi T^4}{36\omega ^2}+\frac{ T^2 M^2}{8 \pi \omega ^2}\right) . \end{aligned}$$

(31)

The transport coefficient obtained here is of order $2 \pi T D\sim 0.3$, which is again on the low side, the perturbative resummation from [43, 44] gets in this case $2 \pi T D\sim 25$ and lattice results ranges from $2 \pi T D\sim 1$–6 depending on the analytic continuation method used [24, 45–47].

4.3 Matching the massless limit with lattice results

The full spectral function for $M^2=\delta M_T^2=g^2C_\mathrm{{F}}T/4$ with $T=1.46T_c$ is shown in 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.

Fig. 4

Left Spectrum at leading and next-to-leading order for $M=0.9T$ (respectively green dashed and blue continuous line, errorband given by renormalisation scale variation as in [28]). Inset Fit of the low energy domain with a Lorentzian. Right NLO Euclidean correlator for $M=0.9T$ (green continuous line) together with continuum extrapolated quenched lattice data [45] (blue dots). In this particular plot, to improve the comparison to lattice data we use a four loop running to define the frequency dependence of the coupling ($\Lambda _{QCD}=0.216$, renormalisation scale $\mu =\mathrm{max}(\omega ,\pi 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

In Fig. 4 right, we see first a quite good match at small $\tau$ but then a different slope at large $\tau$, signalling a lack of power in the small $\omega$ region. To picture how far we are from the lattice data, we multiply the $\omega <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 $\omega <2M$ part has been removed). We see that the difference is sizeable but still much smaller than the given lattice error bars.

5 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\sim 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 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 $\omega \le 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.

Appendices

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_{01111}^0$, 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<\tau <\beta$)

$$\begin{aligned} T\sum _{p_n}\frac{e^{i p_n \tau }}{p_n^2+E_p^2}= & {} \frac{n_{\mathrm {B}}(E_p)}{2 E_p}(e^{\tau E_p}+e^{(\beta -\tau )E_p})\nonumber \\ T\sum _{\{p_n\}}\frac{e^{i p_n \tau }}{p_n^2+E_p^2}= & {} \frac{n_{\mathrm {F}}(E_p)}{2 E_p}(e^{\tau E_p}-e^{(\beta -\tau )E_p}) \end{aligned}$$

(A.1)

can be applied. To do that we shift $P-K\rightarrow K$ and introduce the new integration variables S, R and the corresponding delta functions:

(A.2)

The temporal part of the delta function can be rewritten as an integral, where we keep track of time ordering:

$$\begin{aligned}&\delta (r_n-p_n+q_n)\delta (s_n-k_n+q_n)=e^{i \beta (p_n+s_n)}\nonumber \\&\quad \times \int _0^\beta \mathrm{d}\tau _1 e^{i(r_n-p_n+q_n)\tau _1}\int _0^{\tau _1} \mathrm{d}\tau _2 e^{i(-s_n+k_n-q_n)\tau _2}\nonumber \\&\quad +e^{i \beta (p_n+s_n+q_n)}\nonumber \\&\quad \times \int _0^\beta \mathrm{d}\tau _1 e^{i(-s_n+k_n-q_n)\tau _1}\int _0^{\tau _1} \mathrm{d}\tau _2 e^{i(r_n-p_n+q_n)\tau _2}. \end{aligned}$$

(A.3)

The factors $e^{i \beta (p_n+k_n)}=e^{i \beta (p_n+s_n+q_n)}=1$ were added so that the total phase of all the Matsubara frequencies are between 0 and $\beta$. Thanks to the time ordering and this last prescription the sums can be performed with Eq. (A.1) and then the integrals over $\tau _1,\tau _2$ calculated. Remembering that $q_n$ is a bosonic Matusubara frequency we can set $e^{i q_n \beta }\rightarrow 1$.⁵ 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+i a_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

$$\begin{aligned} \frac{f(E_p,E_r,E_k,E_s,\omega )}{(\omega -E_p+E_r)(\omega +E_s-E_k)}. \end{aligned}$$

(A.4)

This term contains a simple pole at $\omega =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

$$\begin{aligned}&\delta ^3(\mathbf{r}-\mathbf{p})\delta ^3(\mathbf{k}-\mathbf{s})\end{aligned}$$

(A.5)

$$\begin{aligned}&\quad \propto \delta (E_r-E_p)\delta (E_s-E_k)\delta ^2(\Omega _r-\Omega _p)\delta (\Omega _s-\Omega _k)\nonumber \\&\quad \propto \delta (E_r-E_p)\delta (E_s-E_k-E_r+E_p)\\&\qquad \times \delta ^2(\Omega _r-\Omega _p)\delta (\Omega _s-\Omega _k)\nonumber \end{aligned}$$

(A.6)

and use the $\delta (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 $\omega \sim 0$ are of the type $(\omega \pm (E_r-E_p))$. The discontinuity can be taken in each term separately using

$$\begin{aligned} \mathop {\text{ Disc } }\left[ \frac{1}{\omega -A}\right] =-\pi \delta (\omega -A) \end{aligned}$$

(A.7)

for simple poles and

$$\begin{aligned} \mathop {\text{ Disc } }\left[ \frac{1}{(\omega -A)^2}\right] =\pi \delta '(\omega -A) \end{aligned}$$

(A.8)

for double poles.

The spatial integrals over $\mathbf r,\mathbf s$ can be performed using the remaining delta functions $\delta ^3(\mathbf{r-\mathbf p})\delta ^3(\mathbf{s - \mathbf k})$. They force us to perform the limit $E_r\rightarrow E_p$ and get as the final result

$$\begin{aligned} I^{0~ (\omega \sim 0)}_{01111}= & {} \int _{p,k} \pi \omega \delta (\omega )\\&\times \left( \frac{n_\mathrm{{F{1}}}'}{4E_p^2 E_{pk}^3}(1-2n_\mathrm{{F{2}}})+\frac{n_\mathrm{{F{1}}}'n_\mathrm{{F{2}}}'}{4 E_p^2 E_{pk}^2}\right) .\nonumber \end{aligned}$$

(A.9)

Appendix B: Master integrals

Following the notations of Ref. [19], we detail here the different master integrals (denoting $E_p^2=p^2+M^2,E_k^2=k^2+M^2,~E_{p-k}^2=(p+k)^2+M^2$, k being the gluon momentum). First the leading order sum integrals:

$$\begin{aligned} I_{11}= & {} \pi \int _p\left[ \omega \delta (\omega )\frac{n_\mathrm{{F{1}}}'}{2E_p^2}\right. \nonumber \\&\left. + (\delta (\omega -2E_p)-\delta (\omega +2E_p))\frac{1-2n_\mathrm{{F{1}}}}{4 E_p^2}\right] , \end{aligned}$$

(B.10)

$$\begin{aligned} I_{21}= & {} \pi \int _p(1-2\epsilon )\left[ \omega \delta (\omega )\frac{n_\mathrm{{F{1}}}'}{8p^2 E_p^2}\right. \nonumber \\&\left. +(\delta (\omega -2E_p)-\delta (\omega +2E_p))\frac{1-2n_\mathrm{{F{1}}}}{16 p^2 E_p^2 }\right] . \end{aligned}$$

(B.11)

From them a few NLO sum integrals are derived:

$$\begin{aligned} I_{11010}^0= & {} \int _k\frac{1+2n_{\mathrm {B}}(k)}{2k} I_{11},\end{aligned}$$

(B.12)

$$\begin{aligned} I_{12010}^0= & {} \int _k\frac{1+2n_{\mathrm {B}}(k)}{2 k} I_{21},\end{aligned}$$

(B.13)

$$\begin{aligned} I_{01110}^0= & {} \int _k\frac{1-2n_{\mathrm {F}}(E_k)}{2E_k} I_{11},\end{aligned}$$

(B.14)

$$\begin{aligned} I_{02110}^0= & {} \int _k\frac{1-2n_{\mathrm {F}}(E_k)}{2E_k} I_{21}. \end{aligned}$$

(B.15)

The truly NLO sum integrals are given below using the notation $\Delta _\pm =E_p\pm E_{pk}$ and $\Delta _{\sigma \tau }= k+\sigma E_p+\tau E_{pk}$. Here I only give the terms proportional to $\omega \delta (\omega )$ when the other terms contributing at $\omega >0$ can be read from [19],

$$\begin{aligned} I_{10110}^{0~(\omega \sim 0)}= & {} I_{11001}^{0~(\omega \sim 0)}=0,\end{aligned}$$

(B.16)

$$\begin{aligned} I_{11110}^{0~(\omega \sim 0)}= & {} \pi \omega \delta (\omega )\int _{k,p} \frac{n_\mathrm{{F{1}}}'}{8 E_p^2 E_{pk} k}\nonumber \\&\times \left\{ [1+n_\mathrm{{B{0}}}-n_\mathrm{{F{2}}}]\left( \frac{1}{\Delta _{-+}}+\frac{1}{\Delta _{++}}\right) \right. \nonumber \\&\left. -[n_\mathrm{{B{0}}}+n_\mathrm{{F{2}}}]\left( \frac{1}{\Delta _{--}}+\frac{1}{\Delta _{+-}}\right) \right\} , \end{aligned}$$

(B.17)

$$\begin{aligned} I_{11110}^{1~(\omega \sim 0)}= & {} 0,\end{aligned}$$

(B.18)

$$\begin{aligned} I_{11111}^{0~(\omega \sim 0)}= & {} \pi \omega \delta (\omega )\int _{k,p}\biggl [ \frac{(k^2-E_p^2-E_{pk}^2)n_\mathrm{{F{1}}}'n_\mathrm{{F{2}}}'}{4 E_p^2 E_{pk}^2\Delta _{++}\Delta _{-+}\Delta _{+-}\Delta _{--}}\nonumber \\&+\frac{n_\mathrm{{F{1}}}'}{2 E_{pk} E_p^2 k}\nonumber \\&\times \biggl \{ \frac{1}{4E_{pk}^2}\left( \frac{\Delta _{++}+E_{pk}}{\Delta _{++}^2}+\frac{\Delta _{-+}+E_{pk}}{\Delta _{-+}^2}\right) \nonumber \\&+E_{pk} n_\mathrm{{B{0}}}\left( \frac{1}{\Delta _{+-}^2\Delta _{++}^2}+\frac{1}{\Delta _{--}^2\Delta _{-+}^2}\right) \nonumber \\&+\frac{n_\mathrm{{F{2}}} k}{2E_{pk}^2}\biggl (\frac{2\Delta _+^2-E_p^2+E_{pk}^2- k^2}{\Delta _{--}^2\Delta _{++}^2}\nonumber \\&+\frac{2\Delta _-^2-E_p^2+E_{pk}^2- k^2}{\Delta _{-+}^2\Delta _{+-}^2}\biggr )\biggr \}\biggr ]. \end{aligned}$$

(B.19)

The last sum integrals are given for all $\omega$ for completeness (as they were not all written in a suitable form for our present purpose):

$$\begin{aligned} I_{01111}^0= & {} \pi \int _{p,k}\Biggl [\omega \delta (\omega )\nonumber \end{aligned}$$

(B.20)

$$\begin{aligned}&\times \left\{ \frac{n_\mathrm{{F{1}}}'}{4E_p^2 E_{pk}^3}(1-2n_\mathrm{{F{2}}})+\frac{n_\mathrm{{F{1}}}'n_\mathrm{{F{2}}}'}{4 E_p^2 E_{pk}^2}\right\} \nonumber \\&\quad -\left( \delta (\omega -2E_p)-\delta (\omega +2E_p)\right) \nonumber \\&\times \frac{[1-2n_\mathrm{{F{1}}}][1-2n_\mathrm{{F{2}}}]}{8 E_p^2 E_{pk}\Delta _+\Delta _-} \Biggr ], \\ I_{-11111}^0= & {} 2 I^0_{01110}+\frac{1}{2}(\omega ^2-4M^2)I^0_{01111} \end{aligned}$$

(B.21)

and

$$\begin{aligned} I_{12110}^{0}= & {} \pi \int _{k,p}\left[ \omega \delta (\omega )\Biggl (\frac{1}{16 E_p^4 E_{pk} k} \left[ n_\mathrm{{F{1}}}'-\frac{1}{2}E_pn_\mathrm{{F{1}}}''\right] \right. \nonumber \\&\times \left\{ [1+n_\mathrm{{B{0}}}-n_\mathrm{{F{2}}}]\left( \frac{1}{\Delta _{-+}}+\frac{1}{\Delta _{++}}\right) \right. \nonumber \\&\left. -[n_\mathrm{{B{0}}}+n_\mathrm{{F{2}}}]\left( \frac{1}{\Delta _{--}}+\frac{1}{\Delta _{+-}}\right) \right\} \nonumber \\&-\frac{n_\mathrm{{F{1}}}'}{8 E_p^2 E_{pk} k} \left\{ \frac{1+n_\mathrm{{B{0}}}-n_\mathrm{{F{2}}}}{\Delta _{++}\Delta _{-+}}\left( \frac{1}{\Delta _{-+}}+\frac{1}{\Delta _{++}}\right) \right. \nonumber \\&\left. -\frac{n_\mathrm{{B{0}}}+n_\mathrm{{F{2}}}}{\Delta _{+-}\Delta _{--}}\left( \frac{1}{\Delta _{--}}+\frac{1}{\Delta _{+-}}\right) \right\} \Biggr ) \nonumber \\&+\frac{[\delta (\omega -2E_p)-\delta (\omega +2E_p)]}{32k}\nonumber \\&\times \Biggl (\frac{E_p^2(k^2{+}(1{-}4\epsilon ) E_{pk}^2{-}E_p^2{+}2M^2){-}4 E_{pk}^2 p^2}{4 E_p^4 E_{pk}^3 p^2}\nonumber \\&\times (1-2n_\mathrm{{F{1}}})\left\{ [1+n_\mathrm{{B{0}}}-n_\mathrm{{F{2}}}]\left( \frac{1}{\Delta _{-+}}+\frac{1}{\Delta _{++}}\right) \right. \nonumber \\&\left. -[n_\mathrm{{B{0}}}+n_\mathrm{{F{2}}}]\left( \frac{1}{\Delta _{--}}+\frac{1}{\Delta _{+-}}\right) \right\} \nonumber \\&-\frac{E_p^2+E_{pk}^2-k^2-2M^2}{4 E_p^2 E_{pk}^2 p^2}\left( 1-2n_\mathrm{{F{1}}}\right) n_\mathrm{{F{2}}}'\nonumber \\&\times \left[ \frac{1}{\Delta _{++}}+\frac{1}{\Delta _{+-}}+\frac{1}{\Delta _{-+}}+\frac{1}{\Delta _{--}}\right] \nonumber \\&+\frac{1-2n_\mathrm{{F{1}}}}{ E_p^3 E_{pk}}\Biggl \{\left[ \frac{1}{\Delta _{++}}+\frac{1}{\Delta _{-+}}\right] \left[ \frac{1}{\Delta _{++}}\right. \nonumber \\&\left. -\frac{1}{\Delta _{-+}}+\frac{1}{E_p}\right] \left[ 1+n_\mathrm{{B{0}}}-n_\mathrm{{F{2}}}\right] \nonumber \\&+\left[ \frac{1}{\Delta _{--}}+\frac{1}{\Delta _{+-}}\right] \left[ \frac{1}{\Delta _{--}}-\frac{1}{\Delta _{+-}}\right. \nonumber \\&\left. -\frac{1}{E_p}\right] \left[ n_\mathrm{{B{0}}}+n_\mathrm{{F{2}}}\right] \Biggr \} \nonumber \\&-\frac{1-2n_\mathrm{{F{1}}}}{2E_p^2E_{pk}^2}\left\{ \left( \frac{\Delta _+}{E_p}+\frac{E_{pk}^2-E_p^2-k^2}{2 p^2} \right) \right. \nonumber \\&\times \left[ \frac{1+n_\mathrm{{B{0}}}-n_\mathrm{{F{2}}}}{\Delta _{++}^2}+\frac{n_\mathrm{{B{0}}}+n_\mathrm{{F{2}}}}{\Delta _{--}^2} \right] \nonumber \\&+ \left( \frac{\Delta _-}{E_p}+\frac{E_{pk}^2-E_p^2-k^2}{2 p^2} \right) \nonumber \\&\left. \times \left[ \frac{1+n_\mathrm{{B{0}}}-n_\mathrm{{F{2}}}}{\Delta _{-+}^2}+\frac{n_\mathrm{{B{0}}}+n_\mathrm{{F{2}}}}{\Delta _{+-}^2} \right] \right\} \Biggr ) \nonumber \\&+[\delta (\omega -\Delta _{++})-\delta (\omega +\Delta _{++})]\nonumber \\&\times \frac{(1+n_\mathrm{{B{0}}})(1-n_\mathrm{{F{1}}}-n_\mathrm{{F{2}}})+n_\mathrm{{F{1}}}n_\mathrm{{F{2}}}}{8 E_p E_{pk} k\Delta _{++}^2\Delta _{-+}^2} \nonumber \\&+[\delta (\omega -\Delta _{--})-\delta (\omega +\Delta _{--})]\nonumber \\&\times \frac{-n_\mathrm{{B{0}}}(1-n_\mathrm{{F{1}}}-n_\mathrm{{F{2}}})+n_\mathrm{{F{1}}}n_\mathrm{{F{2}}}}{8 E_p E_{pk} k\Delta _{--}^2\Delta _{+-}^2} \nonumber \\&+[\delta (\omega -\Delta _{+-})-\delta (\omega +\Delta _{+-})]\nonumber \\&\times \frac{n_\mathrm{{B{0}}}n_\mathrm{{F{1}}}-(1+n_\mathrm{{B{0}}})n_\mathrm{{F{2}}}+n_\mathrm{{F{1}}}n_\mathrm{{F{2}}}}{8 E_p E_{pk} k\Delta _{--}^2\Delta _{+-}^2} \nonumber \\&+[\delta (\omega -\Delta _{-+})-\delta (\omega +\Delta _{-+})]\nonumber \\&\left. \times \frac{n_\mathrm{{B{0}}}n_\mathrm{{F{2}}}-(1+n_\mathrm{{B{0}}})n_\mathrm{{F{1}}}+n_\mathrm{{F{1}}}n_\mathrm{{F{2}}}}{8 E_p E_{pk} k\Delta _{++}^2\Delta _{-+}^2}\right] . \end{aligned}$$

(B.22)

Appendix C: Full NLO result

To obtain the full NLO result form the previous formulae, the integrals over $\mathbf{p},\mathbf{k}$ still have to be performed. From the antisymmetry of the spectral function it is enough to consider the case $\omega \ge 0$. Using rotation symmetry, it is obvious that all the angular integrals are trivial up to the one containing the angle $\theta$ between $\mathbf{p}$ and $\mathbf{k}$ so that we have three non-trivial integrals to perform over $|\mathbf{p}|, |\mathbf{k}|, \cos \theta$. The full result can be split into three parts, a part proportional to $\omega \delta (\omega )$, which we call $\rho ^t_\mathrm{NLO}$ for transport, a part proportional to $\delta (\omega -2E_p)$, called $\rho ^f_\mathrm{NLO}$ for factorised and a last part proportional to one of the $\delta (\omega \pm \Delta _{\pm \pm })$ denoted by $\rho ^p_\mathrm{NLO}$ for phase space, following the notation of [19]. For the terms proportional to $\delta (\omega )$ the three integrals have to be performed but the angular integral happens to be analytically doable [28]. For the ones proportional to $\delta (\omega -2E_p)$ one can constrain the $|\mathbf{p}|$ integral and only two integrals are left but most of them have to be performed numerically. The terms proportional to $\delta (\omega \pm \Delta _{\pm \pm })$ require more work as the domains where the $\omega =\pm \Delta _{\pm \pm }$ constraints can be satisfied are non-trivial.

1.1 Appendix C.1: $\delta (\omega )$ terms

All the terms proportional to $\delta (\omega )$ in the master integrals can be combined according to Eq. (16). After performing an analogous work as in Ref. [28], i.e. integrating by parts and performing the angular integrals we get

$$\begin{aligned} \frac{\rho ^t_\mathrm{NLO}}{4N_\mathrm{c}C_\mathrm{{F}}g^2}= & {} \frac{ \omega \delta (\omega )}{4\pi ^3}\int _0^\infty \mathrm{d}p\, n_{\mathrm {F}}'(E_p)\int _0^\infty \mathrm{d}k\,\Biggr [ \nonumber \\&+ a k\,n_{\mathrm {B}}(k)\left( 1-\frac{3p^2}{E_p^2}\right) \nonumber \\&+\frac{n_{\mathrm {F}}(E_k)}{E_k}\Biggl (a k^2-M^2-a\frac{3 k^2 p^2}{E_p^2}-\frac{M^2 k^2 p^2}{E_p^2E_k^2} \nonumber \\&-\frac{M^2 k p}{2 E_k^2E_p^2}\biggl \{\frac{M^2 E_k^2(a-1)}{p^2} \nonumber \\&-2(a+1)E_k^2+M^2\biggr \}\ln \left| \frac{p+k}{p-k}\right| \Biggr )\Biggr ], \end{aligned}$$

(C.23)

where the parameter a keeps track of the thermal mass shift. Setting $a=0$ is equivalent to performing the thermal mass shift (23), otherwise $a=1$.

1.2 Appendix C.2: $\delta (\omega -2E_p)$ terms

Summing all terms proportional to $\delta (\omega -2E_p)$ occurring in Eq. (16) we get a formula of the form

$$\begin{aligned} \rho ^f_\mathrm{NLO}(\omega )= & {} \int _p\int _k \frac{ g(k,E_p,E_{pk},\omega )}{k E_p E_{pk}} 2 \delta (\omega -2E_p). \end{aligned}$$

(C.24)

While the function g can easily be constructed from the formulae 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 renormalisation. This proceeds as at zero temperature [19] and will not be explained here. The thermal part we calculate here is finite. We can rewrite the previous integral as

$$\begin{aligned} \rho ^f_\mathrm{NLO}(\omega )= & {} \frac{1}{8\pi ^4}\int _0^\infty \mathrm{d}k \int _M^\infty dE_p \int _{E_{pk}^-}^{E_{pk}^+}\mathrm{d}E_{pk} \nonumber \\&\times g(k,E_p,E_{pk},\omega ) 2\delta (\omega -2E_p)\nonumber \\= & {} \frac{1}{8\pi ^4}\int _0^\infty \mathrm{d}k \int _{E_{pk}^-}^{E_{pk}^+}dE_{pk} \nonumber \\&\times g\left( k,\frac{\omega }{2},E_{pk},\omega \right) \theta (\omega -2M), \end{aligned}$$

(C.25)

where $E_{pk}^\pm =\sqrt{E_p^2+k^2\pm 2 p k}$ is replaced by $\sqrt{\omega ^2/4+k^2\pm k \sqrt{\omega ^2-4 M^2}}$ after applying the delta function. The last two integrals are performed numerically and

$$\begin{aligned} \rho ^f_\mathrm{NLO}(\omega )= & {} \theta (\omega -2 M)\left( 1-2n_{\mathrm {F}}\left( \frac{\omega }{2}\right) \right) (\rho ^{\mathrm{vac},f}_\mathrm{NLO}+\rho ^{\mathrm{bos},f}_\mathrm{NLO})\nonumber \\&+\rho ^{\mathrm{ferm},2}_\mathrm{NLO}(\omega ), \end{aligned}$$

(C.26)

where

$$\begin{aligned} \frac{\rho ^{\mathrm{ferm},2}_\mathrm{NLO}(\omega )}{4 N_\mathrm{c}C_\mathrm{{F}}g^2}= & {} \frac{\theta (\omega -2 M)}{8 \pi ^3}\left( 1-2 n_{\mathrm {F}}\left( \frac{\omega }{2}\right) \right) \int _0^\infty \mathrm{d}k\nonumber \\&\times \int _{E_{pk}^-}^{E_{pk}^+}dE_{pk}\Biggl (n_\mathrm{{F{2}}}'\frac{ M^2 \left( 2 M^2+\omega ^2\right) }{E_{pk} (\omega ^2-4M^2) \omega } \nonumber \\&\times \Biggl [k+ \frac{\left( 2M^2+E_{pk} \omega \right) }{2}\left[ \frac{1}{\delta _{++}}+\frac{1}{\delta _{--}}\right] \nonumber \\&+ \frac{\left( 2M^2-E_{pk} \omega \right) }{2}\left[ \frac{1}{\delta _{+-}}+\frac{1}{\delta _{-+}}\right] \Biggr ] \nonumber \\&+n_\mathrm{{F{2}}}\Biggl [\frac{2M^2+\omega ^2}{\omega E_{pk}^2(\omega ^2-4 M^2)}\Biggl \{\frac{\omega M^2 E_{pk}}{2}\nonumber \\&\times \left( \frac{\delta _-}{\delta _{+-}^2}-\frac{\delta _-}{\delta _{-+}^2}-\frac{\delta _+}{\delta _{++}^2}+\frac{\delta _+}{\delta _{--}^2}\right) \nonumber \\&-M^2(M^2-E_{pk}^2(1-a))\left( \frac{1}{\delta _{++}}+\frac{1}{\delta _{+-}}\right. \nonumber \\&\left. +\frac{1}{\delta _{-+}}+\frac{1}{\delta _{--}}\right) +k(a E_{pk}^2-M^2)\Biggr \} \nonumber \\&+\frac{1}{2\omega ^2}\Biggl \{-2 k \omega a+\frac{M^2(2M^2+\omega ^2)}{E_{pk}}\nonumber \\&\times \left( \frac{\delta _+}{\delta _{-+}^2}+\frac{\delta _-}{\delta _{++}^2}-\frac{\delta _-}{\delta _{--}^2}-\frac{\delta _+}{\delta _{+-}^2}\right) + \nonumber \\&+\left( \mathcal {P}\frac{1}{k \delta _-}+\frac{1}{k \delta _+}\right) \left( k^2 \left( 4 M^2+3 \omega ^2\right) \right. \nonumber \\&\left. +8 M^4-2 \omega ^4\right) \nonumber \\&+\left( \frac{1}{k\delta _{--}}+\frac{1}{k\delta _{-+}}\right) \left( -2 k^2 \omega ^2\right. \nonumber \\&\left. +\left( k \omega +2 M^2\right) ^2-\omega ^2 (k-\omega )^2\right) \nonumber \\&+\left( \frac{1}{k\delta _{+-}}+\frac{1}{k\delta _{++}}\right) \left( 2 k^2 \omega ^2\right. \nonumber \\&\left. -\left( 2 M^2-k \omega \right) ^2+\omega ^2 (k+\omega )^2\right) \nonumber \\&-2M^2\omega (1-a)\left( \frac{1}{\delta _{++}}+\frac{1}{\delta _{+-}}+\frac{1}{\delta _{-+}}\right. \nonumber \\&\left. +\frac{1}{\delta _{--}}\right) \Biggr \} \Biggr ]\Biggr ) \end{aligned}$$

(C.27)

with $\delta _{\pm \pm }=2k\pm \omega \pm E_{pk}$, $\delta _\pm =\omega \pm 2E_{pk}$ and, again, setting $a=0$ is equivalent to performing the thermal mass shift (23). The symbol $\mathcal {P}$ means that we treat the pole at $\delta _-=0$ in the principal value sense. Numerically this can be implemented as follows. We split the integral as $\int _{E_{pk}^-}^{\omega /2}+\int ^{E_{pk}^+}_{\omega /2}$, perform a change of integration variable $E_{pk}\rightarrow \omega -E_{pk}$ in the second integral and add it to the first one.

1.3 Appendix C.3: $\delta (\omega \pm \Delta _{\pm \pm })$ terms

Summing all terms proportional to $\delta (\omega \pm \Delta _{\pm \pm })$ occurring in Eq. (16) we get a formula of the form

$$\begin{aligned}&\sum _{\pm ,\pm }\int _p\int _k \frac{f_{\pm \pm }(k,E_p,E_{pk},\omega )}{k E_p E_{pk}} (\delta (\omega -\Delta _{\pm \pm }) \nonumber \\&\quad -\delta (\omega +\Delta _{\pm \pm })). \end{aligned}$$

(C.28)

If we restrict ourselves to $\omega >0$, only half of the $\delta$’s can actually be realised in some domain of the integrals so that the previous sum becomes

$$\begin{aligned}&\frac{\theta (\omega -2M)}{8\pi ^4}\int _0^{k^{1}}\mathrm{d}k \int _{E_p^{1-}}^{E_{p}^{1+}}\mathrm{d}E_p\,f_{++}(k,E_p,\omega -E_p-k,\omega ) \nonumber \\&\quad -\frac{1}{8\pi ^4}\int _{k^{2}}^{\infty }\mathrm{d}k \int _{E_p^{2-}}^{E_{p}^{2+}}\mathrm{d}E_p\,f_{--}(k,E_p,\omega -E_p+k,\omega )\nonumber \\&\quad +\frac{1}{8\pi ^4}\int _{k^3}^{\infty }\mathrm{d}k \int _{E_p^{1+}}^{\infty }\mathrm{d}E_p\,f_{+-}(k,E_p,-\omega +E_p+k,\omega )\nonumber \\&\quad +\frac{1}{8\pi ^4}\int _{k^3}^{\infty }\mathrm{d}k \int _{-E_p^{1-}}^{\infty }\mathrm{d}E_p\,f_{-+}(k,E_p,\omega +E_p-k,\omega ),\nonumber \\ \end{aligned}$$

(C.29)

where the boundaries of the integrals are given by

$$\begin{aligned}&k^{1}=\frac{\omega ^2-4M^2}{2\omega }, \quad k^{2}=\max (0,-k^1) ,\quad k^{3}=\frac{\omega }{2}, \nonumber \\&E_p^{1\pm }=\frac{\omega -k}{2}\pm \frac{k}{2}\sqrt{1-\frac{4M^2}{\omega (\omega -2k)}},\nonumber \\&E_p^{2\pm }=\frac{\omega +k}{2}\pm \frac{k}{2}\sqrt{1-\frac{4M^2}{\omega (\omega +2k)}} \end{aligned}$$

(C.30)

and the functions $f_{\pm ,\pm }$ through

$$\begin{aligned}&f_{++}(k,E_p,\omega -E_p-k,\omega )\nonumber \\&\quad = \pi \Biggl (\frac{4 k M^4}{\omega ^2 (\omega -2 E_p)^2 (2 (E_p+k)-\omega )}\nonumber \end{aligned}$$

(C.31)

$$\begin{aligned}&\qquad +\frac{M^2 (3 \omega -4 E_p)}{\omega (\omega -2 E_p)^2}+\frac{2 (E_p+k-\omega )^2}{(2 E_p-\omega ) (2 (E_p+ k)-\omega )}\Biggr ) \nonumber \\&\qquad \times \left[ (n_\mathrm{{B{0}}}+1) (1-n_{\mathrm {F}}(\omega -E_p-k)-n_\mathrm{{F{1}}})\right. \nonumber \\&\qquad \left. +n_\mathrm{{F{1}}} n_{\mathrm {F}}(\omega -E_p-k )\right] ,\\&f_{--}(k,E_p,\omega -E_p+k,\omega )\nonumber \\&\quad = \pi \Biggl (-\frac{4 k M^4}{\omega ^2 (\omega -2 E_p)^2 (2(E_p-k)-\omega )}\nonumber \end{aligned}$$

(C.32)

$$\begin{aligned}&\qquad +\frac{ M^2 (3 \omega -4 E_p)}{\omega (\omega -2 E_p)^2}+\frac{2 (k-E_p+\omega )^2}{(2 E_p-\omega ) (2 (E_p- k)-\omega )}\Biggr ) \nonumber \\&\qquad \times \left[ n_\mathrm{{F{1}}} n_{\mathrm {F}}(\omega -E_p+k)\right. \nonumber \\&\qquad \left. -n_\mathrm{{B{0}}} (1-n_{\mathrm {F}}(\omega -E_p+k )-n_\mathrm{{F{1}}})\right] ,\\&f_{+-}(k,E_p,E_p+k-\omega ,\omega )\nonumber \\&\quad = \pi \Biggl (\frac{4 k M^4}{\omega ^2 (\omega -2 E_p)^2 (2 (E_p+k)-\omega )}\nonumber \end{aligned}$$

(C.33)

$$\begin{aligned}&\qquad +\frac{ M^2 (3 \omega -4 E_p)}{\omega (\omega -2 E_p)^2}+\frac{2 (E_p+k-\omega )^2}{(2 E_p-\omega ) (2 (E_p+ k)-\omega )}\Biggr ) \nonumber \\&\qquad \times [-(n_\mathrm{{B{0}}}+1) n_{\mathrm {F}}(E_p+k-\omega )\nonumber \\&\qquad +n_\mathrm{{F{1}}} n_\mathrm{{B{0}}}+n_\mathrm{{F{1}}} n_{\mathrm {F}}(E_p+k-\omega )], \\&f_{-+}(k,E_p,\omega +E_p-k,\omega )\nonumber \\&\quad = \pi \Biggl (-\frac{4 k M^4}{\omega ^2 (2 E_p+\omega )^2 (2 (E_p- k)+\omega )}\nonumber \\&\qquad +\frac{ M^2 (4 E_p+3 \omega )}{\omega (2 E_p+\omega )^2}+\frac{2 (E_p-k+\omega )^2}{(2 E_p+\omega ) (2 (E_p- k)+\omega )}\Biggr ) \nonumber \\&\qquad \times \left[ n_\mathrm{{B{0}}}n_{\mathrm {F}}(E_p-k+\omega )\right. \nonumber \\&\qquad \left. -n_\mathrm{{F{1}}} (n_\mathrm{{B{0}}}+1)+n_\mathrm{{F{1}}} n_{\mathrm {F}}(E_p-k+\omega )\right] . \end{aligned}$$

(C.34)

Note that at $T=0$ only the first integral in (C.29) contributes. This term, setting $T=0$, added to the expression $\rho ^{vac,f}_\mathrm{NLO}$ in Eq. (C.26) gives the full zero temperature result $\rho ^{vac}_\mathrm{NLO}$ given in Eq. (27). For $M\gg T$ the first two integrals contribute. If we take in these terms the part proportional to $n_{\mathrm {B}}(k)$ and not containing any Fermi–Dirac distribution function and add them to the expression $\rho ^{\mathrm{bos},f}_\mathrm{NLO}$ in Eq. (C.26), we get the ’bosonic’ thermal correction $\rho ^\mathrm{bos}_\mathrm{NLO}$ given in Eq. (28). The remaining terms are exponentially suppressed if $M\gg T$ but dominate the spectrum at small $\omega$.

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\rightarrow E_p-\omega +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 $\epsilon$ to the lower bound of the k integration and take the limit $\epsilon \rightarrow 0$ after having added all terms.

1.4 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

$$\begin{aligned} \rho _\mathrm{NLO}^{\mathrm{ferm}}= & {} \rho _\mathrm{NLO}^{t}-2n_{\mathrm {F}}\left( \frac{\omega }{2}\right) (\rho _\mathrm{NLO}^\mathrm{vac}+\rho _\mathrm{NLO}^\mathrm{bos}) \nonumber \\&+\rho _\mathrm{NLO}^{\mathrm{ferm},1}+\rho _\mathrm{NLO}^{\mathrm{ferm},2}+\rho _\mathrm{NLO}^{\mathrm{ferm},3}, \end{aligned}$$

(C.35)

where $\rho _\mathrm{NLO}^{t}$ is given in Eq. (C.23), $\rho _\mathrm{NLO}^{vac}$ in (27), $\rho _\mathrm{NLO}^{bos}$ in (28), $\rho _\mathrm{NLO}^{\mathrm{ferm},1}$ in (30), $\rho _\mathrm{NLO}^{\mathrm{ferm},2}$ in (C.27) and the remaining $\rho _\mathrm{NLO}^{frem,3}$ by

$$\begin{aligned} \frac{\rho _\mathrm{NLO}^{\mathrm{frem},3}}{4N_\mathrm{c}C_\mathrm{{F}}g^2}= & {} \frac{\theta (\omega -2M)}{8\pi ^3}\int _0^{k^{1}}\mathrm{d}k \int _{E_p^{1-}}^{E_{p}^{1+}}\mathrm{d}E_p\, \nonumber \\&\times \Biggl (\frac{4 k M^4}{\omega ^2 (\omega -2 E_p)^2 (2 (E_p+k)-\omega )} \nonumber \\&+\frac{M^2 (3 \omega -4 E_p)}{\omega (\omega -2 E_p)^2}\nonumber \\&+\frac{2 (E_p+k-\omega )^2}{(2 E_p-\omega ) (2 (E_p+ k)-\omega )}\Biggr ) \nonumber \\&\times \bigl [(n_\mathrm{{B{0}}}+1) (2n_{\mathrm {F}}(\omega /2)\nonumber \\&-n_{\mathrm {F}}(\omega -E_p-k)-n_\mathrm{{F{1}}}) \nonumber \\&+n_\mathrm{{F{1}}} n_{\mathrm {F}}(\omega -E_p-k )\bigr ] \nonumber \\&-\frac{1}{8\pi ^3}\int _{k^{2}}^{\infty }\mathrm{d}k \int _{E_p^{2-}}^{E_{p}^{2+}}\mathrm{d}E_p\, \nonumber \\&\times \Biggl (-\frac{4 k M^4}{\omega ^2 (\omega -2 E_p)^2 (2(E_p-k)-\omega )} \nonumber \\&+\frac{ M^2 (3 \omega -4 E_p)}{\omega (\omega -2 E_p)^2}\nonumber \\&+\frac{2 (k-E_p+\omega )^2}{(2 E_p-\omega ) (2 (E_p- k)-\omega )}\Biggr ) \nonumber \\&\times [n_\mathrm{{F{1}}} n_{\mathrm {F}}(\omega -E_p+k) \nonumber \\&-n_\mathrm{{B{0}}} (2n_{\mathrm {F}}(\omega /2)-n_{\mathrm {F}}(\omega -E_p+k )-n_\mathrm{{F{1}}})].\nonumber \\ \end{aligned}$$

(C.36)

Appendix D: Mass shift

The results of Ref. [28] for the Euclidean correlator can be modified to account for the mass shift (21). We have to add to $\frac{G_V(\tau )}{4 N_\mathrm{c}C_\mathrm{{F}}g^2}$ in Eqs. (4.4)–(4.5) of [28] the following integral:

$$\begin{aligned} \frac{G^{MS}_V(\tau )}{4 N_\mathrm{c}C_\mathrm{{F}}g^2}= & {} \int _p \delta M_T^2\left[ \frac{p^2}{2E_p^4}\left( D_{2E_p}(\tau )+ 2 T n_{\mathrm {F}}'(E_p)\right) \right. \nonumber \\&\left. +\left( 1+\frac{M^2}{2E_p^2}\right) \frac{\partial _{E_p}D_{2E_p}(\tau )}{2E_p}+\frac{M^2}{2E_p^3}n_{\mathrm {F}}''(E_p)\right] .\nonumber \\ \end{aligned}$$

(D.37)

In Eq. (28) the thermal mass shift has been performed. Without mass shift, $\frac{\rho ^{bos}_\mathrm{NLO}}{4N_\mathrm{c}C_\mathrm{{F}}g^2}$ would contain additionally

$$\begin{aligned}&\frac{\theta (\omega -2M)}{(4\pi )^3\omega ^2}\int _0^\infty \mathrm{d}k\frac{2 n_{\mathrm {B}}(k)}{k} \nonumber \\&\quad \times \left( 4 \omega k^2\sqrt{\omega ^2-4 M^2}-4 k^2\omega \frac{(2M^2+\omega ^2)}{\sqrt{\omega ^2-4 M^2}}\right) . \end{aligned}$$

(D.38)