Heavy Quarkonium in a weakly-coupled quark-gluon plasma below the melting temperature

We calculate the heavy quarkonium energy levels and decay widths in a quark-gluon plasma, whose temperature T and screening mass m_D satisfy the hierarchy m alpha_s>>T>>m alpha_s^2>>m_D (m being the heavy-quark mass), at order m alpha_s^5. We first sequentially integrate out the scales m, m alpha_s and T, and, next, we carry out the calculations in the resulting effective theory using techniques of integration by regions. A collinear region is identified, which contributes at this order. We also discuss the implications of our results concerning heavy quarkonium suppression in heavy ion collisions.


Introduction
Heavy quarkonium has been suggested since long time as a thermometer for the medium that forms at the core of heavy-ion collision experiments [1]. The early arguments were based on the naïve expectation that above the deconfinement temperature the confining part of the quark-antiquark potential vanishes and the Coulomb part turns into a Yukawa potential due to screening. Since the Yukawa potential supports a finite number of bound states depending on the screening (Debye) mass, and the latter is linear in the temperature, it is then clear that the relative fraction of the different heavy quarkonium states observed will depend on the temperature of the medium. In addition, the electromagnetic decays of these states provide a clean experimental signature. The gross picture above appears to be supported by experiments [2].
In the last few years, significant progress has been made in deriving the quark-antiquark potential on a rigorous basis. A model independent study of the real-time static potential was initiated for large temperatures (T ≫ 1/r > ∼ m D ) in [3,4,5,6] and its implications for a QED and QCD plasma studied. For a wider range of temperatures, including lower temperatures, an effective field theory (EFT) study of non-relativistic bound states in a plasma was initiated in [7] for QED and in [8] for QCD in the static limit. The potential obtained in this way differs in many respects from the most commonly used phenomenological potentials (for some reviews see [9,10]). Most remarkably, it develops an imaginary part. At least two mechanisms have been identified that are responsible for the appearance of a thermal width: the Landau-damping phenomenon [3] and the quark-antiquark coloursinglet break up [8]. In particular, it has been pointed out that quarkonium dissociation due to the former rather than screening may be the dominant mechanism at the origin of heavy quarkonium dissociation in a medium [7,11]. These developments motivate us to revisit the physics of heavy quarkonium states in a thermal bath in a more systematic way. We shall focus here on temperatures for which πT is smaller than the typical momentum transfer in the bound states: such temperatures are those reachable at present days colliders [10].
Heavy quarkonium in a medium is characterized both by the scales typical of a nonrelativistic bound state and by the thermal scales. The non-relativistic scales are the inverse of the typical radius of the system 1/a 0 and its typical binding energy E. The thermal scales are the temperature (or multiple of πT ) and the electric screening mass m D , among other lower energy scales, which are not relevant to our discussion. In the weak-coupling regime, which we will assume throughout this work, these scales may be expressed in terms of the strong coupling constant g ≪ 1, the heavy quark mass m, and the temperature T : m D ∼ gT , 1/a 0 ∼ mα s and E ∼ mα 2 s , where α s = g 2 /(4π). Nonrelativistic scales and thermal scales are hierarchically ordered. This allows to investigate the quarkonium properties in a medium using the same systematic framework provided by non-relativistic effective field theories at zero temperature [12].
In this work, we aim at studying heavy quarkonium at finite temperature including the contribution induced by a large but finite quark mass, in this way merging and completing the findings of Refs. [7,8]. We will adopt the same real-time EFT framework of [7,8] and assume for definitiveness the following hierarchy between the thermodynamical and the non-relativistic scales: With this choice, the thermal bath affects the Coulombic bound state as a small perturbation, yet modifying the Coulomb potential. We remark that this temperature is below the melting temperature, which is of order mα 2/3 s [7]. Moreover, this may indeed correspond to the situation of interest in present day colliders. As a consequence of (1.1), in the weakcoupling regime, we have that mg 3 ≫ T ≫ mg 4 , corresponding to mg 4 ≫ m D ≫ mg 5 . We furthermore assume that Λ QCD , the QCD scale, is smaller than m D (although results that do not involve a weak-coupling expansion at the scale m D , which are all the results of the paper before Sec. 6, are valid also for m D ∼ Λ QCD ). A number of different inequalities has been addressed in the Abelian case in [13].
We will concentrate on the energy levels and decay widths and we will determine how they get modified in a thermal bath whose temperature is such that it satisfies the conditions (1.1). In order to be definite, we will further assume (m D /E) 4 ≪ g, in this way keeping small the number of required corrections suppressed by powers of m D /E, and we will evaluate the spectrum with an accuracy of order mα 5 s . The strongest limitation for the practical application of our final results to actual bottomonium and charmonium systems comes from the fact that we use perturbation theory at the ultrasoft scale mα 2 s . Still, we expect them to be relevant for the ground states of bottomonium and, to a lesser extent, charmonium. Some intermediate expressions, for which perturbation theory is only used at the scale T ≫ mα 2 s may have a wider range of applicability. We also assume a vanishing charm quark mass in the bottomonium case (effects of a non-vanishing mass are discussed in [13]).
The paper is organized in the following way. In Sec. 2, we briefly review the Feynman rules of QCD at finite temperature in the real-time formalism. In Sec. 3, we set up the effective field theory that follows from QCD by integrating out the scales m and mα s in the heavy quark-antiquark sector. In Sec. 4, we calculate the contributions to the spectrum coming from the scale T , in Sec. 5, those coming from the scale E and, finally, in Sec. 6, those coming from the scale m D . In Sec. 7, we summarize our results giving the thermal energy shifts and widths up to order mα 5 s .

QCD at finite temperature in the real-time formalism
In this section, we review the Feynman rules of QCD with static quarks at finite temperature in the real-time formalism. Real-time expectation values depend on how the contour of the time integration in the partition function is deformed to include real times. In the paper, we adopt a contour that goes from an initial time t i to a real final time t f , from t f to t f − i0 + , from t f − i0 + to t i − i0 + and from t i − i0 + to t i − i/T . The propagators will be given with this conventional choice of contour. Since the contour has two lines moving along the real time axis, degrees of freedom double in real time and propagators are represented by 2 × 2 matrices. Furthermore we define which is the Bose-Einstein distribution. The non-relativistic propagator of an unthermalised quark-antiquark pair interacting through a potential V (r) reads (see [8], we have added here the kinetic energy) (2. 2) The expression for V (r) depends on whether the quark-antiquark pair is in a color singlet or in a color octet configuration. In the last case, an identity matrix in color space must be understood, and, in either case, an identity matrix in spin space is implicit. Since the [S(k 0 , k)] 12 component vanishes, the quark-antiquark fields of type "2" never enter in any amplitude of physical fields, i.e. fields of type "1". As a consequence, the fields "2" decouple and may be ignored when considering physical amplitudes. Throughout this paper we adopt the Coulomb gauge for our calculations. The free gluon propagator reads [14]: where k stands for the modulus of the three momentum k i . Due to the decoupling of the heavy quarks of type "2", in the following we will need only the "11" component of the heavy quark-antiquark propagators. At the order we are calculating, this is also the case for the gluon propagators, which will be needed up to one loop. All our equations will refer to this component, unless explicitly stated otherwise. In particular, we recall that at equilibrium the "11" component of the gluon propagator can be written in terms of the retarded (R) and advanced (A) propagators, as which holds for the tree level propagator as well as for the full one. The second term on the right-hand side, proportional to the difference between the retarded and advanced propagators, is often termed the symmetric propagator.

Integrating out the scales m and mα s
Our aim is to calculate the quarkonium spectrum in a thermal bath of temperature T . We take advantage of the hierarchy of scales (1.1) by constructing a hierarchy of effective field theories that follow from QCD by systematically integrating out the largest scales. The EFTs are constructed as series of operators whose matrix elements scale like the lower scales and that are suppressed by powers of the large scales, which have been integrated out.
The first scale to be integrated out from QCD is the heavy quark mass m. In the matching procedure, smaller scales are expanded. Thus, at this stage, the presence of the thermal scales does not affect the matching of the Lagrangian, which is the Lagrangian of non-relativistic QCD (NRQCD) [15].
The next scale to be integrated out is the scale of the inverse of the typical distance of the heavy quark and antiquark, which is of order mα s . According to (1.1), it is larger than the temperature. We are thus allowed to integrate out mα s from NRQCD setting to zero all thermodynamical scales. Furthermore, under the assumption that mα s ≫ Λ QCD , this integration can be carried out in perturbation theory order by order in α s . The EFT we obtain is potential non-relativistic QCD (pNRQCD) [16,17]. Its Lagrangian reads The fields S = S 1 c / √ N c and O = O a T a / √ T F , are the quark-antiquark singlet and octet fields respectively, n f is the number of light quarks q i , N c = 3 is the number of colours, The trace is intended over colour and spin indices. Gluon fields depend only on the center-ofmass coordinate R and on time; this is achieved by a multipole expansion in the relative distance r. The dots in the last line of Eq. (3.1) stand for higher orders in this expansion.
The dependence on the hard and soft scales m and 1/r is encoded in the Wilson coefficients; V A and V B are at leading order V A = V B = 1, whereas the singlet and octet Hamiltonians have the form where m is the heavy quark mass, p = −i∇ r . The dots stand for higher-order terms in the expansion in 1/m, both for the kinetic terms (relativistic corrections) and for the potentials, as well as for terms that depend on the center of mass three momentum.
The static potentials read where C F = (N 2 c − 1)/(2N c ) and α Vs and α Vo are series in α s and at the leading order α Vs = α Vo = α s . α Vs is known up to three loops [18,19], whereas α Vo to two loops [20]. Starting from order α 4 s , α Vs is infrared divergent. This divergence was first identified in [21] and analyzed in the framework of pNRQCD in [22], where it was shown to cancel against an ultraviolet (UV) divergence coming from the ultrasoft degrees of freedom (the scale E).
The non-static potentials V (1) s and V (2) s can be read from [23,12]. V consists of a sum of many terms, such as a p-dependent term, terms depending on the angular momentum, on the heavy quark-antiquark spins and a spin-orbit term. Some of these terms, as well as V (1) s , have an infrared divergence. The leading logarithmic dependence on ln(µr) accompanying these divergences can be read from [23,24]. In the octet sector, the non-static potentials are not known beyond tree level. Fortunately, for the present analysis only the leading order expression in α s for the static octet potential will be needed.
The power-counting of the pNRQCD Lagrangian (3.1) goes as following: the relative momentum p and the inverse distance 1/r have a size of O (mα s ), whereas the time derivative, the gluon fields and the center-of-mass momentum P scale like the lower energy scales. Therefore, the largest term in the singlet potential expanded in α s and 1/m is the Coulomb potential −C F α s /r. The spectrum of the corresponding singlet Hamiltonian is given by the Coulomb levels Subleading terms in the expansions in α s and 1/m are treated in quantum-mechanical perturbation theory. The corresponding shifts of the Coulomb levels have been computed in [23,24,25]. The infrared divergences mentioned above affect the spectrum at order mα 5 s . For what concerns the propagators of the singlet and octet fields in the real-time formalism, we have shown in Eq. (2.2) that the "12" component of a quark-antiquark propagator in a potential V (r) vanishes and the unphysical "2" component decouples. We are thus allowed to drop also here the real-time formalism indices and write only the "11" component of the propagator. For the rest of the paper, all amplitudes will be intended as the "11" components of the real-time matrices unless otherwise specified. In particular, for what concerns the singlet propagator, we thus have where E is the singlet energy. From Eq. (3.2), the singlet Hamiltonian h s reads where the dots stand for higher-order terms. In order to have a homogeneous power counting in the propagator, it is convenient to expand it around the leading-order Hamiltonian (3.6), which is of order mα 2 s . Similarly the octet propagator is where the octet Hamiltonian h o reads 8) and the dots stand for terms smaller than mα 2 s . The Wilson line in (3.7) can be expanded in powers of g. We will only need the leading order in such expansion,

Integrating out the temperature
In this section, we proceed to integrate out modes of energy and momentum of the order of the temperature T . This amounts to modifying pNRQCD into a new EFT where only modes with energies and momenta lower than T are dynamical. We may denote the new EFT with pNRQCD HTL [26]. The EFT can be used for mα s ≫ T ≫ E, m D no matter what the relation between E and m D is. Its Lagrangian will get additional contributions with respect to pNRQCD. For our purposes, we are interested in the modifications to the singlet sector, corresponding to a thermal correction δV s to the singlet potential, and to the Yang-Mills sector, amounting to the Hard Thermal Loop (HTL) Lagrangian L HTL [27]: where we have set to one the matching coefficients of the dipole terms, whose quantum corrections are beyond the accuracy of the present paper. We calculate the correction δV s to the singlet potential. As in [8,7], the leading thermal correction is due to the dipole vertices O † r · gE S + S † r · gE O in the pNRQCD Lagrangian (3.1). These terms induce the diagram depicted in Fig. 1, where a colour-singlet state emits and reabsorbs a chromoelectric gluon through the dipole vertex and an intermediate colour-octet state. The amplitude reads (see [23,28,24] for the T = 0 case) where E is the energy of the singlet; we recall that this expression corresponds to the "11" component in the real-time formalism. Integrals over momenta are regularized in dimensional regularization, with D = 4 + ǫ and µ being the subtraction point. In Coulomb gauge, with the free propagators given in Eqs. (2.3) and (2.4), the contribution of the longitudinal gluon vanishes in dimensional regularization, whereas that of the transverse gluon can be divided into a vacuum and a thermal part: the first term in the square brackets is the vacuum part and the second term is the thermal part. The expression depends on the scales T and E. In order to single out the contribution from the scale T , which comes from the momentum regions k 0 ∼ T and k ∼ T , we recall that T ≫ (E − h o ) and expand the octet propagator as The contribution of the vacuum part of the propagator is scaleless for all the terms of the expansion and thus it vanishes. Conversely, in the thermal part, we have the Bose-Einstein distribution giving a scale to the integration. The zeroth-order term in the expansion (4.4) gives a vanishing integral [8], whereas the following terms contribute to the potential. The linear and the cubic terms in and can be shown to contribute to the real part of the potential. Since in our counting (4.5) behaves as mg 8 ≫ α s T 2 Er 2 ≫ mg 10 and (4.6) as α s E 3 r 2 ∼ mg 10 , further terms in the E/T expansion are not needed. Finally, the square term in the expansion, which would give an imaginary contribution to the potential, vanishes in dimensional regularization: We now evaluate the linear term defined in Eq. (4.5). The integration yields Matching the singlet propagator in pNRQCD with the singlet propagator in pNRQCD HTL we obtain where the left-hand part of the equality corresponds to the pNRQCD part of the matching and the right-hand side to the pNRQCD HTL part of the matching: δV s is the thermal correction to the singlet potential in pNRQCD HTL and δZ s the thermal correction to the singlet normalization in pNRQCD HTL . Our purpose is solely the evaluation of δV s , which is necessary for the spectrum.
is given by the difference between the octet and singlet potentials: ∆V is organized as an expansion in α s and 1/m. At the leading order, it is the difference between the tree-level static potentials: Higher-order terms are easily shown to contribute to the spectrum beyond our accuracy.
Similarly, for what concerns the singlet Hamiltonian, only the leading terms displayed in Eq. (3.6) are necessary. Hence ; the second term is easily identified as contributing to δV s , whereas plugging the first term back into Eq. (4.9) yields The term {r 2 , (E − h s )} contributes to the normalization of the wave function. We are thus left with 1 The commutator can be easily computed from the Hamiltonian displayed in Eq. (3.6).
Plugging the result back into Eq. (4.9), we obtain from the matching condition (4.14) The first term is the contribution of ∆V and was first obtained in [8]. The second term is the contribution of the kinetic term; a similar term appears in the Abelian case of Ref. [7]. Using first-order quantum-mechanical perturbation theory and the expectation values r n,l on the eigenstates of the Coulomb potential (n and l stand for the principal and angular momentum quantum numbers respectively, see, for instance, [29]) we obtain the following correction to the Coulomb energy levels We now move to the cubic term, as defined in Eq. (4.6). We have where I T comes from the evaluation of the integral. It reads [7] where γ E is the Euler's gamma. The divergence of this expression is of infrared (IR) origin: it arises when integrating over the Bose-Einstein distribution at momenta much smaller than the temperature. Since we are integrating out the temperature, i.e. getting the contribution for k ∼ T , this divergence is an artifact of our scale separation. We identify two possible schemes in which the cancellation of this divergence may be interpreted.
1. In the first scheme, the divergence is cancelled by an opposite ultraviolet divergence from a lower scale, in our case the binding energy. In the next section, we will indeed show that the thermal part of this very same diagram, when evaluated for loop momenta of the order of the binding energy, yields an ultraviolet divergence that exactly cancels the one here, whereas the vacuum part of that diagram gives an opposite UV divergence that cancels the IR divergence of the pNRQCD potentials, yielding a finite spectrum.
2. Alternatively one can observe that the pole of the divergence is exactly opposite to the infrared pole of the pNRQCD potentials, which can be read from [23] and the two therefore cancel. More precisely, the scaleless, and hence vanishing in dimensional regularization, integral of the vacuum part of Eq. (4.3), with the octet propagator expanded at the cubic order, can be rewritten as the sum of an infrared and an ultraviolet divergent integral. The infrared pole cancels with the one in Eq. (4.16) coming from the thermal part, whereas the ultraviolet one cancels the IR divergence of the pNRQCD potentials.
The two interpretation schemes are equivalent and produce at the end a finite spectrum, which is the relevant observable. The evaluation of r i (E − h o ) 3 r i in (4.16), in analogy to what has been performed previously in Eqs. (4.12) and (4.13), can be read from [24] 1 where the dots stand for wave function renormalizations. Matching to the right-hand side of Eq. (4.9), we obtain the corresponding contribution to the singlet potential δV s of pNRQCD HTL : Using first-order quantum-mechanical perturbation theory and the value of the Coulomb wave function at the origin, |ψ n,l (0)| 2 = δ l0 /(πn 3 a 3 0 ), we obtain the shift of the energy levels (4.20) Another possible contribution to the potential up to order mα 5 s is given by radiative corrections to the diagram shown in Fig. 1. At the next order in α s , corresponding to two loops, a sizable number of diagrams appears. In [8], it was shown that in Coulomb gauge only one diagram needs to be considered. It consists of a one-loop self-energy insertion in the longitudinal part of the chromoelectric correlator and it is shown in Fig. 2. It contributes at order α s T m 2 D r 2 , corresponding, in our scale hierarchy, to a magnitude in between mg 9 and mg 12 . Therefore, this term contributes to the spectrum up to order mα 5 s only if mg 3 ≫ T ≥ mg 10/3 . This makes clear that non-static contributions that were not considered in the analysis of Ref. [8], such as vertices originating from the spatial centerof-mass covariant derivative in the octet sector and higher-order singlet-octet vertices in the 1/m expansion (see [30]), contribute to terms smaller than mα 5 s only. At our accuracy, the diagram we are considering can again be written expanding the octet propagator for k 0 ∼ T ≫ (E − h o ) and retaining only the first term, independent of E − h o . Therefore, the result is the same as the one derived in [8]. It reads where ζ is the Riemann zeta function (ζ(2) = π 2 /6) and the Debye mass m D is defined as Equation (4.21) contains an imaginary part. It comes from the imaginary part of the gluon self-energy, which is related to the Landau-damping phenomenon, i.e. the scattering of particles carrying momenta of order T in the thermal bath with virtual, space-like longitudinal gluons. Furthermore, the imaginary part is infrared divergent. In the EFT framework, this divergence has to be cancelled by an opposite ultraviolet divergence coming from a lower scale. In the following section, we will indeed show that the same diagram, when integrated over momenta of the order of the binding energy, yields the desired UV divergence. Finally, we remark that the result in Eq. (4.21) comes from dimensionally regularizing only the integral over k while keeping the thermal part of the gluon self energy, which is finite, in exactly four spacetime dimensions. Using the same regularization when calculating the contribution coming from the binding-energy scale guarantees that the final result for the spectrum is finite and scheme independent. This is not the case for the potential, however, whose expression depends on the adopted scheme.

Summary
Summing up Eqs. (4.14), (4.19) and (4.21) we obtain the thermal correction to the potential in pNRQCD HTL up to terms whose contribution to the spectrum is smaller than mα 5 s : where the first two terms come from the linear part of Fig. 1, the terms in square brackets come from the cubic term and the last three lines originate from the diagram in Fig. 2.
This correction to the potential can be used for T ≫ E, m D no matter what the relative size between E and m D is. Analogously, the total contribution to the energy levels coming from the scale T is (4.26) The first and the second lines originate from the diagram in Fig. 1, and correspond to the linear and cubic terms in the expansion (4.4). The last line originates from the gluon self-energy diagram in Fig. 2, which also gives the full contribution of the scale T to the width:

Contribution to the spectrum from the scale E
After having integrated out the temperature in the previous section, many different scales (E, m D , Λ QCD , . . .) still remain dynamical in pNRQCD HTL . In our hierarchy, the binding energy is much larger than the Debye mass and Λ QCD is smaller than all other scales. Our purpose is to compute the correction to the spectrum and the width coming from the scales E and m D . This is achieved by computing loop corrections to the singlet propagator in pNRQCD HTL . We recall that the gauge sector of pNRQCD HTL coincides with the Hard Thermal Loop effective Lagrangian. The longitudinal and transverse gluon propagators in Coulomb gauge are given in the Hard Thermal Loop effective theory by [31] and and the upper sign refers to the retarded propagator and the lower sign to the advanced one. The "11" component can be obtained from the relation (2.7). We start by evaluating the diagram shown in Fig. 1, whose general expression is given in Eq. (4.2), but now the longitudinal and transverse gluon propagators are given by Eqs. (5.1) and (5.2). As we shall see, this is the only diagram we need to consider to get the spectrum at order mα 5 s . At the energy scale, we have k 0 ∼ (E − h o ) and therefore we have to keep the octet propagator unexpanded. However two expansions are still possible.
1. Since k ∼ E ≪ T , the Bose-Einstein distribution can be expanded in In the following, we will call δΣ s (E) the contribution of the diagram in Fig. 1 to the singlet self energy; the corresponding energy shift and width for the state |n, l are given by δE n,l = n, l|Re δΣ s (E n,l )|n, l and Γ n,l = −2 n, l|Im δΣ s (E n,l )|n, l .
We now proceed to the evaluation of Eq. (4.2) for loop momenta of the order of the binding energy, with the HTL propagators defined in Eqs. (5.1) and (5.2). We find convenient to compute separately the contributions coming from the transverse and longitudinal gluons.

Transverse gluon contribution
The contribution of transverse gluons to Eq. (4.2) is in pNRQCD HTL δΣ (trans) We start by evaluating the contribution of the symmetric part, which turns out to be the leading one: where we have expanded the Bose-Einstein distribution. The expansion of the HTL propagators for m D ≪ k 0 , k needs to be performed with care in the region around the light cone, where the gluon propagator becomes singular. We refer to Appendix A for details on the expansion and the evaluation of the integral, whose final result reads The suppressed term of order α s r 2 E 4 /T comes from the k/(12 T ) term in the expansion of the thermal distribution, whereas the term of order α s T m 4 D r 2 /E 2 comes from subleading terms in the expansion of the propagator. 2 We now consider the first term in the square brackets in Eq. (5.5); it does not depend on the Bose-Einstein distribution and, when expanded for k 0 , k ∼ E ≫ m D , gives where P stands for the principal value prescription. Plugging Eq. (5.8) back into Eq. (5.5) yields Summing up Eqs. (5.7) and (5.9) we obtain the complete contribution of the transverse modes We remark that the contribution of the transverse modes at the energy scale is imaginary and finite, in contrast with what happens at zero temperature, where it is real and UV divergent, the divergence cancelling the infrared divergences appearing in the static, 1/m and 1/m 2 potentials at the scale 1/r. This is related to the discussion made in the previous section regarding the cancellation of the IR divergence in Eqs. (4.16) and (4.20) and can be understood in the following way. For E ≫ m D , the Hard Thermal Loop transverse propagator can be expanded for small m D , giving, at the zeroth order, . When plugged in Eq. (5.5) we obtain Eq. (4.3). Evaluated at the binding energy scale, the vacuum part is UV divergent and can be read from [28,23]: 11) where the logarithm of the energy gives rise to the so-called QCD Bethe logarithm in the spectrum [28,24]. On the other hand, the temperature-dependent part gives where the term proportional to r i (E − h o ) 2 r i comes from the first term in the expansion of the Bose-Einstein distribution and the one proportional to r i (E − h o ) 3 r i comes instead from the second term in that expansion, see (5.4). In the sum of Eqs. (5.11) and (5.12) the real parts, divergences included, cancel out and the imaginary parts combine to give the two m D -independent terms of Eq. (5.10). This shows that the binding energy scale contribution produces two opposite UV divergences. In terms of the two interpretation schemes discussed in the previous section, we may understand the cancellation of divergences in two possible ways. In the first way, the vacuum divergence in Eq. (5.11) cancels the IR divergences of the potentials, whereas the UV matter divergence in (5.12) cancels the IR matter divergence from the scale T in (4.16). In the second way, we consider the real part of the potential in pNRQCD HTL as finite, the IR divergences from the scales 1/r and T cancelling each other, and no UV divergences coming from the energy scale, which, as shown by Eq. (5.10), is indeed the case. We stress that the cancellation of the divergences between the vacuum and thermal parts in Eqs. (5.11) and (5.12) is due to the second term in the low-momentum expansion of the Bose-Einstein distribution, i.e. −1/2, which is known in thermal field theory to cause cancellations with the vacuum contribution. Finally, we observe that an analogous cancellation is also obtained in the Abelian case [7]. In order to obtain the contribution to the width from Eq. (5.10), we need to evaluate r i (E − h o ) 2 r i . We proceed as in the previous section and rewrite ( where the dots stand for contributions that vanish on the physical state. The width thus reads where the first line is the contribution from the term proportional to r i (E − h o ) 2 r i , the second line comes from the cubic term and has been obtained using Eqs. (4.18) and (4.20), and the third line is the contribution from the last term in the first line of Eq. (5.10).
The leading contribution to Eq. (5.14) is given by the first three terms, which are of the same size. The first term comes from the static potential and agrees with the one calculated in [8]. The second and third terms come from the kinetic energy; the second one agrees with the one calculated in [7]. This contribution to the thermal decay width originates from the possible break up of a quark-antiquark colour-singlet state into an unbound quark-antiquark colour-octet state: a process that is kinematically allowed only in a medium [8]. Clearly, the singlet to octet break up is a different phenomenon with respect to the Landau damping, which, in the previous section, provided another source for the in medium thermal width. In the situation E ≫ m D , which is the situation of interest for this work, the singlet to octet break up provides the dominant contribution to the thermal width. Indeed, comparing the Landau-damping width (4.24) with the singlet to octet break-up width (5.14), we see that the latter is larger than the former by a factor (mα 2 s /m D ) 2 .

Longitudinal gluon contribution
The contribution of the longitudinal gluons to Eq. (4.2) is where D R,A 00 (k) is the HTL propagator in (5.1). The first term in square brackets, i.e. (D R 00 + D A 00 )/2, does not depend on the Bose-Einstein distribution; therefore only the expansion in m D ≪ E, corresponding to m D ≪ k 0 , k, is possible. We then have (D R 00 + D A 00 )/2 = i/k 2 + O m 2 D /k 4 . The first term is the free propagator, which gives a scaleless integration, whereas the second one can be shown to contribute at order α s Em 2 D r 2 , which is smaller than mα 5 s . For what concerns the symmetric part of the propagator, i.e. (1/2 + n B (k 0 ))(D R 00 − D A 00 ), it should be noted that the retarded and advanced propagators depend on k 0 only through the HTL self-energy; therefore, imaginary parts in their denominators can enter only through the logarithm appearing in Eq. (5.1). Hence, the symmetric propagator is non-zero solely in the spacelike k 2 > k 2 0 region, which is related to the Landau-damping phenomenon. At leading order in the expansions of the Bose-Einstein distribution and of the propagator for m 2 D /k 2 ≪ 1, we thus have 16) The first term contributes to the spectrum at order α s T m 2 D r 2 , so further terms in Eq.
(5.16) are not needed (see footnote 2). We then have For what concerns the width, we observe that the divergence is of ultraviolet origin and cancels the one in Eq. (4.24), yielding a finite width; some care is, however, required in the handling of the logarithms of the energy, which give rise to an analogue of the Bethe logarithm. We have where E 1 = −mC 2 F α 2 s /4 is the energy of the ground state and n, l|r|k is the matrix element between a (bound) eigenstate |n, l of h s and a continuum eigenstate |k of h o . This expression can be reduced to a single integral using the techniques of [28,24]. We obtain for a singlet nS state and an octet P wave (the matrix element introduces a ∆l = 1 selection rule) The definitions of Y E n , X 2 n for n = 1, 2, 3 and ρ n can be found in [28] and [24], the latter reference correcting some misprints in the former. A numerical evaluation of these integrals for the three most tightly bound l = 0 states yields:

Summary
In summary, the contribution to the energy levels coming from the binding energy scale is entirely due to the longitudinal part of the chromoelectric correlator, where the first two lines come from the first two in Eq. (5.14) and the last two from Eq. (5.19) and from the last term in (5.14). I n,l is defined in Eq. (5.20).

Contributions to the spectrum from the scale m D
In our hierarchy of energy scales, the next scale after the binding energy is the Debye mass.
We thus have to evaluate Eqs. (5.5) and (5.15) for momenta of the order of m D . In detail, we have two regions to analyze: the first one is k 0 ∼ E − h o , k ∼ m D , corresponding to having the octet propagator unexpanded and conversely expanding the HTL propagators for k 0 ≫ k. It can be easily shown that both the transverse and the longitudinal parts result in a series of scaleless integrations over k, which vanish in dimensional regularization. The second region corresponds to having k 0 ∼ m D and k ∼ m D : the octet propagator then needs to be expanded, whereas the HTL propagators are to be kept in their resummed form. The resulting integrals are quite involved, however, by power counting arguments, it can be easily seen from Eqs. (5.5) and (5.15) that, once the octet propagator is expanded, the largest term comes again from the symmetric part of the gluon propagator, due to the T /k 0 enhancement factor. The size of this term turns out to be of order α s T m 3 D r 2 /E and, since we have assumed (m D /E) 4 ≪ g, it is beyond mα 5 s .

Conclusions
We have computed the heavy quarkonium energy levels and widths in a quark-gluon plasma of temperature T such that mα s ≫ T ≫ mα 2 s ≫ m D . Assuming (m D /E) 4 ≪ g, the spectrum is accurate up to order mα 5 s . The thermal shift of the energy levels induced by the medium is obtained by summing the contribution from the scale T , given in Eq. (4.26), with the thermal part of the contribution from the energy scale. We remark that the contribution from the energy scale, given in Eq. (5.24), is the sum of both vacuum and thermal contributions, which, in the transverse sector, cancel. The thermal contribution of the transverse modes can be derived from Eq. (5.12). The complete thermal contribution to the spectrum up to order mα 5 s reads (we recall that E n = −mC 2 F α 2 s /(4n 2 ) and a 0 = 2/(mC F α s )) where L n,l is the QCD Bethe logarithm, defined as [28,24] L n,l = 1 We refer to [28,24] for details on the numerical evaluation of this integral. We furthermore remark that the thermal contribution to the spectrum is finite, the IR divergence in Eq. (4.26) having cancelled against the UV divergence coming from Eq. (5.12). The thermal width is obtained by summing the contribution from the scale T , given in Eq. (4.27), with the one coming from the energy scale as given in (5.25), the IR divergence in the former cancelling against the UV divergence in the latter. We then have × a 2 0 n 2 5n 2 + 1 − 3l(l + 1) where I n,l is defined in Eq. (5.20). We remark that, up to the order considered here, the thermal contribution to the spectrum and to the width is independent of the spin. Our results are expected to be relevant for the ground states of bottomonium (Υ(1S) and η b ), and to a lesser extent to those of charmonium (J/ψ and η c ), for a certain range of temperatures in the quark-gluon plasma for which (1.1) is fulfilled. Let us now try to figure out what our results imply for the electromagnetic decays to lepton pairs or to two photons. First of all, the masses of the heavy quarkonium states increase quadratically with the temperature at leading order (first line of (7.1)), which would translate into the same functional increase in the energy of the outgoing leptons and photons if produced by the quarkonium in the plasma. Second, since electromagnetic decays occur at short distances (∼ 1/m ≪ 1/T ), the standard NRQCD factorization formulas hold, and, at leading order, all the temperature dependence is encoded in the wave function at the origin. The leading temperature correction to it comes from first-order quantum-mechanical perturbation theory of the first term of (4.25). The size of this correction is ∼ n 4 T 2 /(m 2 α s ). Hence, a quadratic dependence on the temperature should also be observed in the frequency in which leptons or photons are produced by the quarkonium in the plasma. Finally, at leading order, a decay width linear with temperature is developed (first line of (7.3)), which implies a tendency to decay to the continuum of colour-octet states. Hence, a smaller number of vector and pseudoscalar ground states is expected to be in the sample with respect to the zero temperature case.
i.e. the region where the gluon is far from being on shell. The second region is called the collinear region. In this region, we have We observe that the collinear scale m 2 D /(E − h o ) has, in our energy scale hierarchy, a magnitude in between mg 4 and mg 6 . It is, therefore, smaller that the Debye mass by a factor of m D /E ≪ 1 and still larger than the non-perturbative magnetic mass, which is of order g 2 T , by a factor T /E ≫ 1. For simplicity, we separate the two regions by a cut-off Λ, such that We start by analyzing the off-shell region. Here k 2 0 − k 2 = λ(2k + λ) ≫ m 2 D and we can thus expand the retarded propagator propagator in Eq. (5.2) as (A.5) Terms contributing to the real part of this propagator and hence to ∆ R − ∆ A can come either from the poles of the denominators, yielding a δ(k 2 0 − k 2 ), or from the imaginary part of the logarithm. However, δ(k 2 0 − k 2 ) = 0 over the whole off-shell region. We can safely discard these terms and obtain (A.6) Note that the principal value prescription is irrelevant since our integration region excludes the poles. From Eq. (A.1), we get δΣ (trans, symm) s,off shell This integral does not need to be dimensionally regularized, so we can set D = 4 at this point and obtain δΣ (trans, symm) s,off shell (A.8) We consider, now, the collinear region. We start again from the retarded propagator introduced in Eq. (5.2). We perform the change of variables k 0 − k = λ and we expand for λ ∼ m 2 D /k ≪ k, thereby implementing the collinear hierarchy. We then have (∆ R − ∆ A )(k 0 > 0) = ∆ 1 + ∆ 2 + ∆ 3 + ∆ 4 + ∆ 5 + O λ k 3 , (A.9) where the ∆ i are defined as (A.14) We start by plugging ∆ 1 in Eq. (A.1). We then have δΣ (trans, symm) s,1 The contribution of ∆ 2 is δΣ (trans, symm) s,2 (E) = g 2 C F T r i 3 the leading order term in the expansion of ((E − h o + iη) 2 − (k + λ) 2 ) −1 , which would contribute at order α s T m 2 D r 2 , vanishes because the integral over λ is zero. The contribution of ∆ 3 is δΣ (trans, symm) s,3 (E) = − ig 2 C F T m 2 D r i 12 where the dots stand for higher orders. We remark that the dependence on the cut-off scale Λ has disappeared. The contribution of ∆ 4 is δΣ (trans, symm) s,4 The needed λ integrals are (E) has only contributions that are suppressed by powers of 1/Λ. Finally, the contribution of ∆ 5 is δΣ (trans, symm) s,5