Heavy Quarkonium at finite temperature and chemical potential

We generalize known results for heavy quarkonium in a thermal bath to the case of a finite baryonic density, and provide a number of formulas for the energy shift and decay width that hold at weak coupling for sufficiently large temperature and/or chemical potential. We find that a non-vanishing decay width requires a temperature larger than the typical binding energy, no matter how large the chemical potential is. This implies that at zero temperature the dissociation mechanism of heavy quarkonium is due entirely to screening, unlike in the finite temperature case. We use several effective theories in order to sort out the contributions of the relevant energy and momentum scales. In particular, we consider contributions of the so called quasi-static magnetic modes. The generalization to the case of a finite isospin/strangeness chemical potential is trivial. We discuss possible applications to the SIS and NICA conditions, and compare with available lattice results.

We generalize known results for heavy quarkonium in a thermal bath to the case of a finite baryonic density, and provide a number of formulas for the energy shift and decay width that hold at weak coupling for sufficiently large temperature and/or chemical potential. We find that a non-vanishing decay width requires a temperature larger than the typical binding energy, no matter how large the chemical potential is. This implies that at zero temperature the dissociation mechanism of heavy quarkonium is due entirely to screening, unlike in the finite temperature case. We use several effective theories in order to sort out the contributions of the relevant energy and momentum scales. In particular, we consider contributions of the so called quasi-static magnetic modes. The generalization to the case of a finite isospin/strangeness chemical potential is trivial. We discuss possible applications to the SIS and NICA conditions, and compare with available lattice results.

I. INTRODUCTION
High energy heavy-ion collision experiments (HIC) have shown the existence of collective behavior in the strong interactions, namely a new state of matter that is usually refer to as quark-gluon plasma (QGP) (see [1] for a review). In order to study the properties of the QGP, the so called hard probes have been extremely useful. Among them, the suppression of heavy quarkonium states in the products of the HIC were proposed as a signal of QGP formation long ago [2] (see [3][4][5] for reviews). Nowadays the sequential suppression of Υ(1, 2, 3) has been clearly observed [6,7]. However, the QCD dynamics which is actually responsible for this suppression is not so easy to identify. The original proposal of Matsui and Satz [2] that screening would be the mechanism behind sequential melting of quarkonium bound states is not entirely correct. In [8], it was shown that for a weakly interacting QGP the same dynamics that produces screening, also produces an imaginary part to the potential, as a consequence of the so called Landau damping. In [9], it was emphasized that this imaginary part is parametrically larger than the real part and hence Landau damping rather than screening should be regarded as the key mechanism for heavy quarkonium dissociation. Imaginary potentials were also obtained in strongly coupled QGP settings [10][11][12] and included in models of in-medium heavy quarkonium [13][14][15]. Later on, within the weakly coupled QGP, detailed analysis were made taking into account the interplay between temperature, screening mass and the different scales in the quarkonium bound state dynamics, the main lesson being that finite temperature effects cannot always be incorporated in a phenomenological potential [9,[16][17][18]. The effects of a relative velocity of the heavy quarkonium with respect to the medium have also been analysed [19,20] (see also [21]). Recently, part of these findings have been embedded in a more realistic framework of an expanding QGP, no matter whether this is weakly or strongly coupled [22][23][24][25][26][27][28][29].
High energy HIC experiments are essentially gluon colliders, and the resulting medium has a negligible baryonic chemical potential with respect to the temperature scales attained. In the near future, however, there are planned HIC experiments at lower energies that aim at attaining large values of the chemical potential at SIS (CBM) [30] and Nica (MPD) [31] (see [32] for a recent overview). These colliders will have energy enough to produce charmonium bound states [33]. It is then worth exploring in a solid theoretical framework, namely using QCD at weak coupling and the well-known effective field theories for heavy quarkonium, the fate of these states at non-zero baryon chemical potential. This is so even if the charm quark mass may not be high enough to apply weak coupling techniques beyond the ground state, or if the values of chemical potential actually attained in the experiments may not be large enough to justify a weak coupling analysis. Indeed, the weak coupling analysis may unravel qualitative new features that may then be incorporated in more realistic models or settings. This was the case, for instance, when an imaginary part of the potential at finite temperature was uncovered in [8].
We shall then restrict ourselves to the study of a heavy quarkonium propagating in a QGP in thermodynamical equilibrium such that the temperatures T and baryon chemical potentials µ fulfill m T g T Λ QCD and m µ gµ Λ QCD , where g is the QCD coupling constant and m is the heavy quark mass. We shall also assume that the heavy quarkonium is weakly coupled, namely m mα s mα s 2 Λ QCD , where α s = g 2 /4π, and is at rest in the QGP rest frame. Recall that p ∼ mα s is the typical heavy quark momentum in the bound state (and r ∼ 1/mα s its typical radius) and E ∼ mα s 2 the typical size of the binding energy. We aim at the calculation of the leading order effects both in µ and in T on the mass (binding energy) and decay width. These calculations are non-trivial because on the one hand at energy scales of the order of or below the typical quarkonium binding energy Coulomb ressummations must be carried out, and on the other hand at momentum scales of the order of or below the Debye mass HTL resummations must be carried out. In that respect, the use of suitable effective field theories [34][35][36] is very convenient.
We shall distribute the paper as follows. In Sec. II we set the basic formalism, in Sec. III and IV we address the two most significant cases, p max(T , µ) E and max(T , µ) p E, respectively. Each section contains subsections where the cases T µ and µ T are separately addressed. Sec. V discusses our results as well as other cases that are not addressed in full detail. The more technical developments are relegated to the appendices.

II. BASIC FORMALISM
Throughout this work we will use the real-time formalism for thermal field theory (for reviews see eg. [37][38][39]), including both the effects of temperature and chemical potential, following the lines of [16,18]. We shall consider the heavy quarks and heavy quarkonium as probe particles, and hence absent in the medium. This means in practice that the real-time non-relativistic propagator reduces for them to the 11 components only, and those take the form of the usual non-relativistic retarded propagators at zero temperature and chemical potential.
In the real-time formalism, the longitudinal and transverse gluon propagators are four-by-four matrices diagonal in color space which at tree level in the Coulomb gauge read (color indices omitted), respectively [40] where K = (k 0 , k) and k = |k|. Note that the longitudinal part of the gluon propagator in Coulomb gauge does not depend on the temperature. We can write them in terms of advanced/retarded propagators, n B being the (temperature-dependent) bosonic occupation number. This formula holds at any order. At tree level, where "R" stands for retarded and "A" for advanced. Note that the property is inherited by the full propagators, and will be often used in the following. When integrating out the hard scale, namely for K max(T, µ), the gauge sector reduces to that of the wellestablished Hard Thermal Loop (HTL) effective theory, whose longitudinal and transverse propagators read respectively, where For a QCD medium at finite temperature and density with N f light quark flavors the expression for the Debye mass m D reads [41] where N c is the number of colors, T F = 1/2 and we have separated the bosonic (m 2 D(B) ∼ N c ) and fermionic (m 2 D(F ) ∼ N f ) contribution to it for later use. At energy and momentum scales smaller than the Debye mass m D , the longitudinal gluons are screened and can be integrated out. However, a particular class of transverse gluons, the so called quasi-static magnetic modes, survive at those lower scales. They fulfill m D k k 0 . Hence, in this case the transverse propagator can be approximated by Having introduced the formalism we will use throughout the paper, we now move on to investigate the effects of a dense medium on quarkonium states.
We start by considering, and extending to finite chemical potential, the case discussed in [17,18,23], namely m p max(T, µ) E. In this case, the T and µ will affect the binding energy and decay width of the quarkonium, but not its size. The heavy quarkonium essentially remains a Coulombic bound state, and the medium effects are perturbations to it. Since scales much larger than max(T, µ) lead to exponentially suppressed Boltzmann factors, we can start our considerations directly from the (vacuum) pNRQCD Lagrangian [34,42], namely with E i = F i0 chromo-electric field. The singlet/octet Hamiltonians are where P and p are the center-of-mass and relative momentum respectively, and the various V (n) are potentials known up to a certain order. In our calculations, P, V B and the subleading potentials (n > 0) can be neglected and we may approximate α s /r, and V A 1. Now we may integrate out the largest scale, T or µ, and get to another EFT which is valid at the lower scales E, m D . The outcome will be a new contribution to the singlet potential: where D µν (k 0 , k) stands for the full gluon propagator in the Coulomb gauge and we are using dimensional regularization (DR) with D = d + 1 = 4 + 2 , ν the DR subtraction scale, and Ω d denoting the solid angle in d spatial dimensions. As mentioned in the previous section, quarkonium and heavy quarks are considered in this work as test particles outside of the medium. As a consequence, in the real-time formalism only the 11 components of the gluon propagators, [D µν ] 11 , will couple to them. In the following, we will omit for brevity the 11 indices and only label explicitly the retarded and advanced components when they appear.

A. Integrating out the hard scale
Integrating out the hard scale (T or µ) will give us a new EFT which we will refer to as pNRQCD HT L , following [18]. In addition, if the hard scale is much larger than E ∼ h o , we can expand the octet propagator 1 , From the general expression Eq. (13), we will consider two contributions: 1. One-loop hard contribution (Fig. 1).
In this case, D µν (k 0 , k) = D  µν (−k 0 , k), the first term in the expansion Eq. (15) leads to a vanishing contribution. The contribution from the second term has been computed in [16,18]. It does not contain any fermionic occupation number -the only medium dependence enters in the n B from the 11 gluon propagator prescription. So there will not be any µ dependence, and we can just take those results. One obtains Since the contribution from the first term of Eq. (15) vanishes for symmetry reasons and Eq. (16) comes from the second-order term in the expansion, there could be higher-loop diagrams that give a contribution comparable to it. We analyze them in the next section.
At two loops, the longitudinal gluons may now contribute because the one-loop self-energy provides them with a T and µ dependence. The transverse gluons give subleading contributions because the would-be leading term in Eq. (15) vanishes for the same symmetry reasons as in the previous section. Hence, Eq. (13) reduces to, where D 00 (k 0 , k) must be calculated at one loop. Moreover, we can write where P denotes the principal value integral. Using again the symmetry properties of D, Eq. (5), we see that only the delta function piece survives. Hence, as long as we work at the lowest order of the expansion Eq. (15), the symmetry of the problem forces k 0 → 0 and we only need to calculate Π 00 (k 0 → 0, k). For the one-loop longitudinal gluon self-energy Π 00 (k 0 , k) we have to sum the gluonic and fermionic loop contributions shown in Fig. 3. With our definitions we can write where "F" labels the contribution coming from the loops of N f massless quarks (first diagram of Fig. 3) and "G" labels the contribution from the second, third and fourth diagram of Fig. 3. The gluon contribution to the self-energy is unchanged by the presence of a chemical potential and can be taken directly from [16]. Our focus will then be on the fermionic contribution, which at µ = 0 and finite T is given by [16], Note that Π 00, F (−k 0 , k) = Π 00, F (k 0 , k). We need to generalize this fermionic contribution to finite µ. Recall that in our real-time Feynman rules (see eg. [44]) the occupation number enters as a sgn( . For µ = 0 this reduces to 1 ± 2n B/F (|q 0 |). In a dense medium the fermionic occupation number instead is given by n F (q 0 − µ) and we face expressions like and if f (−q 0 ) = f (q 0 ), as in Eq. (20), we can just replace in our expressions Furthermore, from the discussion after Eq. (18), we know that we only need the small k 0 limit of the longitudinal gluon self-energy. If we expand Π R 00 (k) and Π A 00 (k) for k 0 k and keep terms up to order k 0 , the result for its real and imaginary parts is Taking the additional k → 0 limit here would lead to the familiar HTL self-energy. However, we have to keep k arbitrary here since we are calculating the hard contribution. We can then write the contribution of our self-energy correction to the 11 propagators as with δD We then end up with At this point we can integrate out the next larger scale. The leading contribution is given by an expression analogous to Eq. (13) in which the gluon propagators correspond to those of the HTL effective theory Depending on whether the next larger scale is the Debye mass m D ∼ g max(T, µ) or the binding energy E ∼ mα s 2 the approximations to be carried out differ. If m D E then E − h o can be expanded in the above expression. If instead m D E, then one can expand the self-energies in the HTL gluon propagators. If m D ∼ E, it is not possible to proceed analytically beyond extracting the UV divergences that cancel the IR ones in Eq. (33), see [9,17]. We shall not further consider this last case.
Qualitatively we can single out two cases: if T is large enough we can extend the formulas obtained in [16,18] for vanishing µ to the case of nonzero chemical potential, be it large (µ ∼ T ) or small (µ T ). The small T case requires some extra care, as we will see in Sec. III C.

B. Large T
We start by computing the hard contribution. For large T (T µ) , as long as we work at the lowest order of the expansion Eq. (15) we can use the k 0 → 0 limit. Then The hard contribution at finite temperature and vanishing chemical potential has been calculated in [16,18]: where ζ is the Riemann Zeta function and we can easily isolate the real and imaginary contributions δV R G , δV I G coming from the gluon loops (∼ N c ) from the the contributions δV R F , δV I F (∼ N f ) coming from the fermionic one. The former will be unchanged by the presence of a chemical potential, so we will focus on the latter.
Eq. (29) is obtained by working out first the k integral using DR in Eq. (26), which helps putting all scale-less integrals to zero, then performing the q integral in Eq. (23) and Eq. (24). Note that these contributions are suppressed by a factor rT with respect to the purely real ones obtained in (16). The real part is finite, while the imaginary one has an IR log divergence. Now let us compute these fermionic contributions at finite µ. We have, where Li denotes the polylogarithm function. Note also that the real part is finite. This will not the case for the imaginary part, which reads where we reconstructed the contribution proportional to the Debye mass, see Eq. (9), which goes together with an additional ∼ T 3 factor multiplying a more involved piece containing the derivative of the polylogarithm functions with respect to their first argument. Note that δV I above contributes to the decay width at leading order whereas δV R is subleading with respect to Eq. (16). Putting together the results above with the gluonic contributions in Eq. (29) and including the leading contribution Eq. (16), we can work out the corrections to the energy levels δE nl = δV R nl as well as the decay rate Γ nl = −2 δV I nl for a given n, l state. We have where we have introduced r nl = a 0 [3n 2 − l(l + 1)]/2 and r 2 nl = a 2 0 n 2 [5n 2 + 1 − 3l(l + 1)]/2, a 0 = 2/(mC F α s ) being the Bohr radius.
The expressions above so far hold for arbitrary µ, as long as T µ. We consider next the expansion for small µ, µ T . Note that the µ dependence in n F (k 0 ) is analytic, and the expansion in µ does not modify its UV and IR behaviour. Hence, the scale µ will not introduce extra singularities in our loop calculations. As a consequence, the results in this case can be obtained by just expanding in µ Eq. (32) and Eq. (33) above.
For the real part, since the leading term Eq. (16) does not depend on µ so it remains the same. The µ dependence arises from the next-to-leading term Eq. (30), The first term in the expansion indeed corresponds to the µ = 0 result (cfr. Eq. (29)). For the imaginary part, we have from Eq. (31) and we recover the fermionic part (∼ N f ) of Eq. (29) from the µ → 0 limit of Eq. (34) and Eq. (35). Keeping terms up to O(µ 2 /T 2 ) only, we have for the energy and the decay rate contributions The expressions Eq. (32) and Eq. (33), which reduce to Eq. (36) and Eq. (37) above in the T µ limit, are the outcome of integrating out the hard scale in the heavy quarkonium sector. In the gluonic sector the outcome is the celebrated HTL effective theory. This effective theory has exactly the same form for µ = 0 as for µ = 0, the only difference being that the Debye mass acquires a µ dependence, as displayed in Eq. (9). In the case T µ, m D ∼ gT and m D can also be expanded in a series of (µ/T ) 2 .
Beyond the contributions at the hard scale, there will be additional contributions to the energy shifts and decay widths from lower scales. The form of these contributions will depend on the relative size between m D and mα s 2 , but not on the size of µ because of its analytic dependence. Below the hard scale T we can use HTL for the light degrees of freedom. Hence all the µ dependence will be in m D except for the case µ m D in which there will be additional analytic dependences arising from the fermionic distribution function in HTL fermion loops. The latter however will be suppressed by g 2 factors. Let us next discuss the two most extreme cases.

mD E
In this case, we can further integrate out the scale m D to get additional modifications to the potentials. These quantities basically depend on m D (except for the T factor in the imaginary part coming from the n B , which is unchanged), so the inclusion of a chemical potential simply amounts to considering the appropriate expression for the Debye mass in a dense medium. This has been worked out in [17] (see also [23]) for QED. We simply take the results from there, Eqs. (10)- (11), and correct for QCD color factors. We display directly the corrections to the energy shift and decay width below, Putting together the results above with Eq. (32) and Eq. (33), we get the final result for this case for the energy and the decay rate, Note that the 1/ pole in the imaginary part Eq. (39) cancels with the one of Eq. (33).
For T µ, we can just add the soft (∼ m D ) scale contributions, Eqs. (38) and (39), to the hard contribution Eq. (36) and Eq. (37) to obtain our final result, again only up to order µ 2 /T 2 : If we consider even lower scales, we find that the contribution at the scale E ∼ mα s 2 m D may only be due to quasi-static magnetic photons Eq. (10) and is of order α s r 2 T E(Em D ) 1/3 , and hence suppressed with respect to the contributions calculated so far. Therefore, our final results in this case are Eq. (40) and Eq. (41), which reduce to Eq. (42) and Eq. (43) for T µ. In order to get a feeling on the contributions computed in this section, we plot in Fig. 4 the results for the energy shift and the decay rate as function of the chemical potential for different values of α s .

E mD
If E ∼ mα s 2 > m D we should be integrating out this scale first rather than the Debye mass. This has been worked out at µ = 0 in [17] (Eqs. (6)- (7)) for QED, and in [18] for QCD. In this case the denominator Eq. (15) cannot be expanded, so that we are no longer fixed to the k 0 → 0 limit, but we can still make use of the expansion Eq. (28) for the bosonic occupation number. Furthermore, we can employ the HTL gluon self-energies expanded in powers of m D /E. The leading energy shift is given by the longitudinal gluon contribution, Eq (5.18) in [18], whereas both longitudinal and transverse gluons contribute to the decay width, which is given by (5.25) in [18], where E n = −mC 2 F α s 2 /(4n 2 ), n = 1, 2, . . . , is the energy of the state and I n,l a numerical constant dependent on the n, l state given in [18]. We introduced the shorthand notation the latter being a subleading T -independent contribution which we will neglect in the following. Again the generalization to finite µ is straightforward: we can clearly distinguish single factors of T , which come from the expansion of the bosonic distribution function in the gluon propagator and thus are unchanged by the introduction of density, whereas all the remaining medium dependence is expressed in terms of the Debye mass. We can just replace m D by the appropriate µ−dependent value for the Debye mass, given by Eq. (9) .
We then get out final result for this case by adding the above expressions to Eq. (32) and Eq. (33), If we consider lower scales, in this case the scale m D , we find that it gives contributions of the order α s r 2 T m 3 D /E which are suppressed with respect to the ones calculated so far. Hence, our final results in this case are given by Eq. (48)   Quarkonium energy shift and decay rate for the ground state as function of the ratio µ/T including E scale corrections, for N f = 2. We plot in the left panel the result (δE10 − δE10|µ=0)/(4αs 2 CF T 3 r 2 10/3), and in the right one (Γ10 − Γ10|µ=0)/(−4παs 2 T 3 CF r 2 10/9). This time our ratio for the energy shift turns out to be αs-independent, while for the decay rate we need to specify the value of the ratio T /mαs 2 : here we chose for illustration 2 (black), 5 (blue) and 10 (red).

C. Small T
Let us next consider the case µ T . This case is technically more involved as it cannot be obtained by just taking the small T limit of the general case (recall that the distribution functions are not analytic in T ). Below the scale µ, HTL must be used but the approximation N B (k 0 ) ∼ 2T /k 0 for the Bose distribution function does not hold in general. This is because the constraint k 0 µ still allows k 0 T . Let us study the two extreme cases, T E and T E below.

T E
The contributions at the hard (µ) scale are now restricted to the quark loop in Fig. 2, see sec. III A 2. The leading contribution can be worked out by just replacing in Eq. (30) and Eq. (31) n F (q − µ) + n F (q + µ) → θ(µ − q). The leading T dependence requires more effort, one can nevertheless work it out. The real part becomes up to exponentially suppressed terms, and the imaginary part reads We see that the overall factor of T coming from the bosonic occupation number makes this imaginary contribution parametrically smaller than its real counterpart. At scales below µ the HTL effective theory must be used for gluons and light quarks. If T m D ∼ gµ , E, we can next integrate T out. The selfenergies can then be expanded in powers of m D in the HTL propagators, which at leading order become the QCD ones. Hence, we obtain an extra contribution, which coincides with Eq. (16). Since the last result is finite, we still need to integrate out the next larger scale in order to cancel the 1/ pole in Eq. (53). This is done in Sec. III C 1 a and III C 1 b below. We address the case m D T E in Sec. III C 1 c. Note that for µ T the gluon distribution functions are exponentially suppressed at the hard scale, and hence they do not contribute to the HTL selfenergies. Therefore, in the following subsections, all the Debye masses will only have the fermionic contributions.

a. The case T m D ∼ gµ E
We can take the result of integrating out m D from Eq. (38) and Eq. (39). Putting everything together we finally obtain Note that, parametrically, in the energy shift, the two first terms and the third term in the first line compete to be the leading contribution, whereas the remaining ones are suppressed. The last term in the decay width is also suppressed.

b. The case T E m D ∼ gµ
We can take the result of integrating out E from Eq. (44) and Eq. (45). Putting everything together we finally obtain Parametrically, in the energy shift, all terms may compete for the leading order, except for term linear in µ, which is always smaller than the term cubic in µ. In the decay width, the terms proportional to C (1) nl and m 2 D(F ) are the leading and next-to-leading ones respectively, while the one proportional to T 3 is suppressed. c. The case µ m D ∼ gµ T E In this case, the next scale to be integrated out is m D . This produces the same contribution as in the general case, namely Eq. (38) and Eq. (39) with m D(F ) instead of the full m D , as in the previous subsection. Putting everything together we obtain Parametrically, the first term both in the energy shift and in the decay width above is the leading one. Concerning the T -dependent terms, since we have not considered contributions at the scale T so far, we may wonder whether the terms above are the most important ones, or contributions at lower scales may provide larger T -dependent terms. In order to resolve this question, let us next consider the contributions at lower scales. Only quasi-static magnetic modes survive below m D . These are obtained by approximating the transverse HTL self-energy to the case k 0 k m D , see Eq. (10). When m D T , these modes contribute at lower scales from the diagram in Fig. 6. At the scale T , their contribution is of the order α s r 2 T 2 (T m 2 D ) 1/3 , and hence it becomes the most important T -dependent contribution to the real part of the potential 2 . It reads which produces a further energy shift to be added to Eq. (58), There are also contributions at the scale E from the quasi static magnetic modes, which are of order α s r 2 T E(Em 2 D ) 1/3 , and hence suppressed with respect to the ones considered so far.

E T
The case T E ∼ mα s 2 deserves a special treatment. Let us focus first on the the longitudinal contribution of Eq. (13), which becomes where we introduced for brevity a ≡ E − h 0 . In order to single out real and imaginary parts we consider the combinations δV R = (δV + δV * )/2 and δV I = −i(δV − δV * )/2.
After writing the denominator as i a−k0+iη = −iP 1 k0−a + πδ(a − k 0 ), we thus get to The k 0 integral of δV I can be carried out and leads to Recall now that a = E − h o ∼ E T . We can then approximate N B (a) ∼ sgn(a) = −1, for any a < 0, as it is the case for a bound state, since E < 0 and h o is positive definite. Hence the imaginary part is zero, irrespectively of the form of the longitudinal gluon propagator.
A similar reasoning can be done with the transverse contribution, which also leads to a vanishing contribution for the imaginary part. The real part can be obtained from Eq. (63) by replacing k 2 → k 2 0 and D R,A 00 (k 0 , k) → D R,A ii (k 0 , k). Therefore no imaginary part to the potential is generated when T is smaller than the binding energy scale. Note that the same argument is valid both for the two-loop hard contribution as well as for the HTL one, as it does not depend on the the details of the gluon propagator.
Let us first consider the hard contribution k ∼ µ to the real part. The first term in Eq. (63) gives Eq. (52), as expected, since a µ. The second term is subleading. Indeed, for k 0 ∼ k ∼ µ, the denominator can be expanded in a, the would-be-leading order vanishes (odd in k 0 ), and hence the leading contribution is a/µ suppressed. For k 0 k ∼ µ, the difference between retarded and advanced propagators is proportional to k 0 /k, and hence also suppressed. Nevertheless, the region k 0 ∼ T provides the leading T -dependent contribution, where ψ(z) is the Digamma function. Note that δV R hard | T −dependent is not really a potential, as it depends on the external energy E. The soft regions, k 0 , k µ, give subleading contributions, but they do contribute to the leading T -dependence. Let us display the two extreme cases, a. The case µ m D , E T Let us next consider the contributions at the scale k ∼ m D . The first term in Eq. (63) gives Eq. (38), which together with Eq. (53) leads to Eq. (58). However, the leading T -dependence is not given by this expression but by the region k ∼ m D , k 0 ∼ T in the second term of Eq. (63), which together with Eq. (66) leads to the following T -dependent energy shift, The matrix element above has been calculated in [45,46] (see also [47]). The contribution of the transverse photons is T 2 /m 2 D suppressed with respect the one of the longitudinal gluons that we have just displayed.
b. The case µ E T m D In this case the region k ∼ k 0 ∼ T gives the leading T -dependence. In fact, it is due to the transverse gluons because their propagator at tree level already contributes. It gives the same result as for the µ = 0 case, namely, The abelian limit of the expression above agrees with the one of [9]. The size of this term is parametrically larger than the hard contribution in Eq. (66). Nevertheless, it is important to carry out the calculation at the soft scale that cancels the 1/ pole in Eq. (66). This is realized by the longitudinal gluons at the scale k 0 ∼ k ∼ T , which give (70) Putting together Eq. (69), Eq. (66) and Eq. (70), we obtain for the T -dependent energy shift in this case, So far the energy scales associated with the thermal medium were assumed to be smaller than the typical momentum exchanges p between the constituents of the bound state. This implies that thermal effects can be treated as perturbations to the bound state dynamics. The melting of the bound state may still occur, because it can develop a medium decay width comparable to the binding energy. One may wonder however, in which conditions the medium effects will be so strong that they will affect the leading-order bound state dynamics, namely the leading order potential. When p ∼ max(T, µ), this is not the case yet. This is because the longitudinal gluon propagator is not sensitive to the medium at tree level, and hence the Coulomb-like potential remains as the LO potential. The one loop correction is suppressed by a g 2 factor, and hence medium effects are still a perturbation.
For µ = 0, this case is analyzed in Sec. IV of [9] and in Sec. IIb/Appendix D of [17] for QED. We shall not develop it further, since it does not bring in any qualitative difference with respect to the previous section. In contrast, the case max(T, µ) p ∼ m D introduces modifications in the LO potential, and hence in the full bound state dynamics. For µ = 0, in the static limit of QCD (m → ∞, p ≡ 1/r), this case was addressed in [8] and [16], and in the full dynamical case of QED in [9,17]. In the following we extend these results to finite chemical potential.
The suitable starting point now is Non-Relativistic QCD (NRQCD) [35,43], since the heavy quark mass is still larger than the remaining scales in the problem, and hence it can be integrated out. We will only need the leading order Lagrangian, where ψ is a non-relativistic field that annihilates heavy quarks, and c.c. stands for the charge conjugated term, namely the analogous terms for the heavy antiquarks, see [48,49].

A. Integrating out the hard scale
In the gluon and light quark sector, the integration of the largest scale max(T, µ) produces the HTL effective theory. In the heavy quark sector, it produces a shift of the heavy quark mass δm. In the static limit, the leading contribution corresponds to the two-loop diagram in Fig. 7, which is O(α s 2 max(T, µ)) and turns out to be suppressed by a factor of g with respect to lower energy contributions. However, when 1/m corrections are considered, there is a leading order contribution from the diagram of Fig. 8 We now proceed to integrating out the lower scales. As before, it is useful to treat separately the cases in which T is large and small, respectively.

B. Large T (T µ)
This corresponds to calculating the mass shift and potentials (O(α s m D ∼ mα s 2 ∼ E) ) using HTL. The result can then be just read from [8,16]. For the mass shift we get, while for the potential shift, where now m D may depend on both µ and T . Recall that the potential develops an imaginary part, first uncovered in [8]. If T ∼ µ, our final result for the potential plus mass shift is just the addition of twice the (complex) mass shifts in Eq. (73) and in Eq. (74), and the potential in Eq. (75). This is also the case if T µ. Then the leading µ dependence is obtained by expanding m D in µ 2 /T 2 . Note that in these cases the imaginary part of the potential is parametrically larger than the real part if m D ∼ 1/r ∼ p, hence bound states can only exist if 1/r m D , and cease to exist when this imaginary part takes over the real part, that is before the screening mechanism r ∼ m D is sizable [9]. Even if heavy quarkonium cannot be considered a bound state anymore, its spectral function can be calculated from the evolution by its non-hermitian Hamiltonian [50].
The case µ T , however, deserves a separate discussion. First of all, the hard contribution from Fig. 8 should be dropped since T is not hard anymore. When T → 0 the imaginary part of the potential Eq. (75) and mass shift Eq. (74) vanish, and one may naively think that the bound state is stable. But this need not be so. On the one hand, there could be subleading contributions that do not vanish in this limit, and on the other hand, before T reaches zero, there are additional scales that play a role, in particular the binding energy E. Let us analyze the following two cases separately, p ∼ m D T E and p ∼ m D E T . The case T p ∼ m D E reduces to expanding the results for the T ∼ µ case in T /µ.

T E
At the scale m D we still get the same result as in Eq. (75). The T factor in the imaginary part comes from the Bose enhancement in the gluon distribution function n B (k 0 ) ∼ T /k 0 for k 0 T , which still holds since k 0 ∼ E, the typical energy transfer, and E T . However, now the imaginary part of the potential is parametrically smaller than the real part, and one may wonder whether T -independent contributions to the imaginary part exist that compete in size with Eq. (75). The leading T -independent contributions to the imaginary part of the mass shift come from Fig.  7 (hard scale) and Fig. 9 (m D scale), when the internal heavy quark line is on-shell. They are O(g 2 m 2 D /m). Since The T dependence from the hard scale is encoded in m D . The leading T dependence in the real part of Eq. (75) is ∼ ET 2 /µ 2 . Since µ is the largest scale in the problem, one may wonder whether other contributions from lower scales are larger. In order to address this question we must take into account that below the m D scale the only low energy degrees of freedom in the light sector are the quasi-static magnetic gluons Eq. (10). Furthermore, below the scale p ∼ 1/r we can use pNRQCD with the mass shifts and singlet potential given in Eq. (74) and Eq. (75) respectively, and similar modifications to the octet potential, At the scale T , there is a contribution from Fig. 6, in which singlet and octet propagators must be understood with the potentials described above, This contribution is O(α s r 2 ET (T m 2 D ) 1/3 ) and hence parametrically larger than ∼ ET 2 /µ 2 (recall that m D ∼ p ∼ 1/r ∼ mα s implies that µ ∼ gm). Then the leading T -dependence to the energy shift is given by the expectation value of the expression above and the decay width by minus twice the expectation value of the imaginary part of Eq. (75), The expectation values above are calculated with the real part of Eq. (75) in the Hamiltonian. We have used that r i (h o − E)r i = Ncαs 2 re −m D r + 3 m on physical states in Eq. (77). Let us finally mention that there are parametrically larger T -dependent contributions to the mass shift, O(g 2 m 2 D /m) from the one-loop self-energy diagram in the region k ∼ m D and k 0 ∼ T . However these contributions are logarithmic in T and hence very smooth. In addition, they are difficult to calculate. We have displayed in Appendix B the logarithmically enhanced contributions. There are similar subleading contributions (∼ g 4 E) from the region k 0 ∼ T , k ∼ µ that in some particular cases may compete with Eq. (77) as well.

E T
In this case, the imaginary part of the tree-level potential turns out to be zero. This is because all relevant scales are bigger than T and hence N B (k 0 ) ∼ sgn(k 0 ). Then the imaginary part of the tree-level potential becomes proportional to absolute value of the transfer energy, which is zero for on-shell heavy quarks in the center of mass frame. Regarding the contribution from the heavy quark selfenergy, one can then work along the same lines as in Sec. III C 2 in order to prove that the imaginary part vanishes at one loop, as mass-shift contributions would be of the same form as Eq. (65), now with a = E − k 2 /2m (and similarly for the transverse contribution). The leading corrections to the imaginary part may arise from the vertex correction ( fig. 10), the two gluon exchange diagrams ( fig. 11) and the two-loop heavy quark selfenergy (fig 12). We prove in Appendix A that they also vanish. Therefore the imaginary part of the potential and of the mass shift vanish at leading order and including O(α s ) corrections.
Let us next focus on the temperature dependence of the energy shift. The potential depends on temperature through the Debye mass, which gives a T -dependent contribution to the energy shift of O(g 2 m D T 2 /µ 2 ). Since µ is the largest scale in the problem after the heavy quark mass, we may expect more important contributions from lower scales. We find that the leading T -dependent contribution comes from the one-loop self-energy diagram in which the longitudinal gluon propagator has k ∼ m D and k 0 ∼ T , which is O(g 2 T 2 /E). This contribution is difficult to calculate. On the one hand the energy scale k 0 ∼ T E, and hence bound state effects cannot be ignored. On the other hand, pNRQCD cannot be straightforwardly used since k ∼ m D ∼ p ∼ 1/r, and hence the multipole expansion does not hold. In order to avoid the last problem we shall restrict ourselves to the particular case m µ p m D E T . The Tdependent part of the energy shift can be obtained from the second term in Eq. (63), where now a = E − h o − k 2 /4m. Notice that we have included the quarkonium center of mass recoil energy, which was negligible in Sec. III. For k ∼ m D and k 0 ∼ T the HTL propagator must be used. We obtain, The 1/ arises from an UV divergence in k. It should be compensated by an IR divergence of a contribution at a higher k scale. In Eq. (63), the scale k ∼ −4m(E − h o ) ∼ p m D is also relevant. It allows to make an expansion in m D in the HTL propagators, which induces the IR divergence we are looking for. We obtain, There is a problem with the result above: for k ∼ p the multipole expansion, on which pNRQCD is based, does not hold. Nevertheless, if we are only interested in the IR behaviour k → 0, namely k p, then it can be used. That means that our calculation above gets the correct IR behaviour, and hence the correct log, but the finite pieces are not reliable. Putting together Eq. (80) and Eq. (81), we then obtain, where the O(1) means there is an unknown number that adds to the logarithmically enhanced contribution.

V. DISCUSSION
We have worked out the modifications in the binding energy and decay width that a QGP at high temperature and/or chemical potential induces in a heavy quarkonium state, generalizing earlier work done in the limit of a vanishing chemical potential. This was done from QCD at weak coupling in the real-time formalism with approximations that are well under control, relying on the hierarchy of scales in the problem. This is in contrast with earlier work on heavy quarkonium at finite chemical potential, in which some modeling is introduced [51]. In particular, we have shown that the rather usual assumption that the medium effects can be encoded in a modified potential, as assumed in the early days [52,53], is not always true.
We have restricted ourselves to heavy quarkonia at rest. The effects of a relative velocity with respect to the thermal bath may eventually be addressed along the lines of refs. [19,20]. In fact, when the effects of the medium can entirely be encoded in a potential, they have already been addressed in [54]. We have focused on a number of cases in which analytic results can be produced. However, it should be clear from our general formulas that numerical results can be also obtained for the remaining cases.
When the temperature and chemical potential are smaller than the typical momentum exchange between the heavy quarks, the medium effects are a perturbation that, in general, cannot be encoded in a potential. This has been already emphasized for zero chemical potential in [9,[16][17][18]. In this case, Coulomb resummations must be always carried out, and the medium effects enter through gluons emitted by chromoelectric dipole transitions which turn a color-singlet quarkonium into a color-octet one or viceversa. In that respect it is very helpful to use pNRQCD. Depending on the energy and momentum of the emitted gluon, HTL resummations may also be necessary. If the chemical potential and the temperature have the same size, the results we obtain are similar to the ones of the zero chemical potential case, but include non-trivial functions of µ/T . If T µ, we can just expand our results in µ/T , as the distribution functions are analytic in µ/T . However, if µ T , the distribution functions are not analytic in T /µ, and this requires extra care. In this limit, the Debye mass m D ∼ gµ may be comparable to T and hence accounting properly for the leading temperature effects requires HTL resummations. We find that if the temperature is larger than the binding energy, the decay width is proportional to T , but it vanishes otherwise.
When the temperature or the chemical potential are larger than the typical momentum exchange between the heavy quarks, the medium effects modify the leading order potential. This is the case addressed in the pioneering works [2], in which the screening was proposed as the mechanism leading to J/ψ suppression. Later on, an important imaginary part due to Landau damping was uncovered for this potential which changed the picture [8]. When T µ the imaginary part of the potential is proportional to g 2 T and parametrically larger than the real part (∼ g 2 m D ). Due to this imaginary part, the heavy quarkonium melts before noticing the screening effects, as in the case of zero chemical potential [9]. When T m D ∼ gµ, screening and Landau damping compete for being the leading effect. The imaginary part of the potential exists as long as the temperature is larger than the binding energy, but it vanishes otherwise. We have been able to prove it at next-to-leading order in α s .
Our analysis turned out to be technically challenging, as Coulomb and/or HTL resummations have been necessary in several instances. The use of effective field theories has been invaluable to keep track of the important terms in a systematic manner. Dimensional regularization has been used to regulate both the IR and UV divergencies that arise in the intermediate steps of the calculations when we factorize the contributions of the different scales. We have obtained contributions from energy and momentum regions that had been ignored so far. In that respect the method of integration by regions developed in [55] (see [56] for a review) has also been very useful. For instance, in order to get the leading temperature effects in the binding energy when µ p ∼ m D E T we needed gluons of energy ∼ T and momentum ∼ m D . These gluons are on the one hand sensitive to the binding energy, and hence Coulomb resummations are required, and on the other hand have a momentum large enough so that the multipole expansion cannot be applied, and hence the calculation cannot be carried out entirely in pNRQCD. We circumvented these difficulties by introducing the extra hypothesis p m D . Another non-trivial example is the contribution of quasistatic magnetic modes [57,58] when µ p ∼ m D T E that give an important T -dependent piece of the binding energy 3 . Finally, let us mention the logarithmic T -dependence in the same case, which is log enhanced and requires the introduction of an extra regularization to factorize the energy scale from the momentum scale. We have chosen an analytic regularization similar to ref. [59], see Appendix B.
Our results are obtained entirely in the weak coupling regime of QCD and thus may not be straightforwardly applied to realistic experimental situations, especially for charmonium, as some of the scales in the problem may not be large enough. Nevertheless, we believe they provide important constraints to models, as they fix quite a number of asymptotic behaviours for large µ of more realistic models. In the absence of a definitive approach to address real-time phenomena in general and large chemical potentials in particular in lattice QCD, complementary approaches based on weak coupling QCD should be helpful.
Let us then consider J/ψ, which will be observed in most of the planned experiments [32]. If we take m c ∼ 1.6 GeV 4 , then the experimental value of the J/ψ mass delivers E ∼ 0.1 GeV. If we associate this value with a Coulombic state, we obtain α s (p) ∼ 0.4 and p ∼ 1/r ∼ 0.4 GeV. We see that the value of p is very low even if α s is relatively small 5 . For the maximum expected values of the baryon chemical potential µ B quoted in ref. [32], we have µ = µ B /3 0.3 GeV 6 . It means that most of the times we would be in the case of Sec. III, and only when µ ∼ 0.3 GeV, the case of Sec. IV may be relevant. This is of course provided that T µ.
Although analyzing bottomonium does not seem to be in the future experimental plans, some of the colliders feeding the relevant experiments (e.g. NICA, RHIC, SPS) are energetic enough to produce it. If we take m b ∼ 4.9 GeV 7 , then the experimental value of the Υ(1S) mass delivers E ∼ 0.34 GeV. If we associate this value with a Coulombic state, we obtain α s (p) ∼ 0.4 and p ∼ 1/r ∼ 1.2 GeV. Then, for the expected values of the chemical potential, Υ(1S) would always be in the case of Sec. III.
If we stick to qualitative features of our results, the most relevant one is that a temperature larger than the size of the binding energy T > E appears to be necessary for heavy quarkonium to develop a decay width. No decay width is developed if T < E , no matter how large is the chemical potential (provided it is smaller than the heavy quark mass). This may be understood in terms of the Fermi sea: In order to dissociate quarkonium, a light quark of the Fermi sea must provide an energy larger than the binding energy to the bound state. But then it becomes less energetic in the final state, and since all the states with less energy are occupied in the Fermi sea, the process cannot take place. Hence at large chemical potential and small temperature, we only expect modifications in the heavy quarkonium mass (through the binding energy). The dissociation mechanism would be screening, namely the one originally proposed in [2].
For sufficiently heavy quark mass, chemical potential and/or temperature, our results are reliable. In the case of small (zero) temperature and large chemical potential, one should observe in the quarkonium spectral function a shift in the location of each bound state peak with no modifications in the width when we increase µ. This is in contrast with what happens at large temperature and small (zero) chemical potential, in which case, apart from the shift in the location of the bound state peaks, a widening of the peaks is observed when the temperature is increased. In fact, the melting of the bound states occurs because the peaks corresponding to different bound states overlap and lose their identity. This can be understood at weak coupling in terms of the Landau damping [8]. In the case of large chemical potential one would just observe bound states peaks disappearing when we increase the chemical potential. It would be interesting to cross-check our results in lattice QCD simulations, but this would require having overcome the difficulties of dealing with a large chemical potential (see [3,[61][62][63][64][65] for reviews).
However, we can compare with the results of ref. [66], a NRQCD lattice simulation for N c = 2 and N f = 2 and heavy quark mass ma = 5, 4, 3, where a is the lattice spacing. They consider 0 ≤ µa ≤ 1.1 and 1/24 ≤ T a ≤ 1/12, hence we can probe the µ T regime. If we assume that the binding energies are Coulombic, from the values for different masses of ∆Ea at µ = 0 in their Fig. 1, we obtain that α s ∼ 0.65 − 0.7 at the scale of the typical relative momentum p. This implies E ∼ 0.3a and p ∼ 1/r ∼ 1.2a. Hence, most of the data displayed in their Fig. 1 is in the region m p µ E T , and we should compare it with the results in our Eq. (58) and Eq. (68). The left panel of their Fig. 1 shows the binding energy as a function of µa for three values of the heavy quark mass. For these plots to be compatible with Eq. (58), we need the total (i.e. including the one hidden in the Debye mass) coefficient of the µ 3 term to be positive (the temperature can be neglected). This is achieved if α s (m D ) 0.86. If so, our expression qualitatively describes the rising observed from µa ∼ 0.6 to µa ∼ 1. We can also understand the bending downwards around µa ∼ 1: in this region µ ∼ p and with our values of α s (m D ), µ ∼ m D , hence Eq. (74) and Eq. (75) should better be used for the energy shift. If we expand Eq. (75) for m D r 1, the first correction to the Coulomb potential is negative, which may explain the above mentioned downward trend. However, we cannot explain the mild decreasing from µa ∼ 0.3 to µa ∼ 0.6. We would probably need expressions for µ ∼ E that we have not worked out, or it may simply happen that α s becomes too large at those low scales so that our weak coupling description is not appropriate even qualitatively. In any case, the behavior of the curves with the mass is easy to understand as the dependence on the chemical potential goes as µ 3 /m 2 . Hence the smaller the mass is, the more noticeable the effects are, as clearly shown in the left panel. The temperature effects are displayed in the right panel of their Fig. 1. Those should be encoded in Eq. (68), and we indeed see in the plot that rising the temperature increases the binding energy, although we do not observe the quadratic increase of Eq. (68).
Finally, our results can also be applied to the case of a non-vanishing isospin chemical potential µ I rather than a baryon chemical potential. We only have to replace N f µ by 2|µ I | in our equations. This is because our expressions are symmetric under µ ↔ −µ and each light quark contributes the same amount of µ at finite baryon chemical potential. At finite isospin chemical potential, the u-quark contributes by µ I , the d-quark by −µ I . We could then try to compare with the two-flavor lattice results of ref. [67]. However, the results displayed in that reference correspond to µ I 0.3 GeV, a too low scale to apply our weak coupling calculation 8 .

HTL correction to the vertex
Beyond leading order, a possible source of imaginary contributions to the potential is the vertex correction of Fig. 10. Writing the full vertex function as W a = T a (1 + δW ), δW reads whereẼ = E − (k + p) 2 /(2m) andẼ = E − (k + p ) 2 /(2m), (E, p) and (E , p ) being the incoming and outgoing heavy-quark energy and three-momentum respectively. Eventually, we will use thatẼ,Ẽ are much smaller than k, p and p , and that in a bound state E, E < 0. However this limit must be taken once the k 0 integral has been carried out, otherwise we are left with ill-defined expressions. Since in the small-temperature limit (here and in the following we will use the shorthand notation D R/A 00 (k 0 , k) ≡ R/A(k 0 , k) andṘ(k 0 , k) = dR(k 0 , k)/dk 0 ) D 00 (k 0 , k) = 1 2 R(k 0 , k) + A(k 0 , k) + sgn(k 0 ) (R(k 0 , k) − A(k 0 , k)) , we write δW = δW 1 + δW 2 , with δW 2 containing the terms proportional to sgn(k 0 ) and δW 1 all the rest. δW 1 can be evaluated by contour integration, and the smallẼ,Ẽ limit gives, which is purely imaginary. The imaginary part of δW 2 can also be evaluated using the formula, It leads to which cancels exactly (A3). Hence no imaginary part arises from the vertex correction at one loop.
2. HTL two-gluon exchange contributions to the potential The two gluon exchange contributions of Fig. 11 may also provide imaginary parts. Consider first the diagram on the left projected on color singlet states. This diagram contains the iteration of the leading order potential, which must be subtracted, where we recall that D 00 stands for the 11 component of the real-time temporal gluon propagator. Upon writing it in terms of the retarded and advanced propagators we obtain δV = δV 1 + δV 2 , with For δV 1 the integral over k 0 can be carried out, which turns the two heavy quark propagators into a i/(E + iη) quarkonium propagator and replaces the k 0 in the retarded and advanced propagators by E/2 and −E/2 respectively. Using E k , q, we finally get Im(δV 1 ) = − iC 2 F 2 d d k (2π) d R(0, k)Ṙ(0, k + q) +Ṙ(0, k)R(0, k + q) , where we have also used R(0, k) = A(0, k) andṘ(0, k) = −Ȧ(0, k). For δV 2 , it is tempting to take E → 0 in the integrand, and then formally show that Im(δV 2 ) = 0. However, it turns out that the integral is ill-defined in that limit, and this naive result is wrong. Instead, one can show that, for E < 0, Im(δV 2 ) = − iC 2