Heavy Quarkonium moving in a Quark-Gluon Plasma

By means of effective field theory techniques, we study the modifications of some properties of weakly coupled heavy quarkonium states propagating through a quark-gluon plasma at temperatures much smaller than the heavy quark mass, m_Q. Two different cases are considered, corresponding to two different hierarchies between the typical size of the bound state, r, the binding energy, E, the temperature, T, and the screening mass, m_D. The first case corresponds to the hierarchy m_Q>>1/r>>T>>E>>m_D, relevant for moderate temperatures, and the second one to the hierarchy m_Q>>T>>1/r, m_D>>E, relevant for studying the dissociation mechanism. In the first case we determine the perturbative correction to the binding energy and to the decay width of states with arbitrary angular momentum, finding that the width is a decreasing function of the velocity. A different behavior characterizes the second kinematical case, being the width of s-wave states a non-monotonic function of the velocity, increasing at moderate velocities and decreasing in the ultra-relativistic limit. We obtain a simple analytical expression of the decay width for T>>1/r>>m_D>>E at moderate velocities, and we derive the s-wave spectral function for the more general case T>>1/r, m_D>>E. A brief discussion of the possible experimental signatures as well as a comparison with the relevant lattice data are also presented.


I. INTRODUCTION
Heavy quarks produced in the early stage of relativistic heavy-ion collisions are valuable probes of the medium that develops at later stages. They can be used to resolve its energy density, and eventually to understand which are the fundamental degrees of freedom of the system. Indeed, at sufficiently high energy densities matter should form a quark-gluon plasma (QGP) and the propagating heavy quarks should be capable to convey this information to us. In relativistic heavy-ion collisions (HIC) it is expected that the variation of the interaction between heavy quarks due to the creation of a hot medium should be observable. In particular, the Debye screening of the Coulomb-like potential between a heavy quark and a heavy antiquark was proposed in Ref. [1,2] as a dissociation mechanism, resulting in the suppression of the yields of heavy quarkonium, QQ, states in HIC. The low-lying heavy quarkonium states are considered as the most powerful probes because they are the only hadronic states that are able to survive above the deconfinement temperature (see [3,4] for reviews). This is due to the fact that even at weak coupling, namely ignoring confinement, these states still exist. In addition, the vector states enjoy a rather clean dilepton decay channel, which makes them easy to identify experimentally, although in HIC the corresponding background is not yet completely understood, see e.g. [5].
Suppression of charmonium states was first observed in Pb-Pb collisions by the NA50 Collaboration [6] at relatively low center of mass energy (per nucleon), √ S N N = 158 GeV. However, in contrast with the naive Debye screening scenario, a complicated pattern emerged because of the various processes involving charm quarks, see e.g. [7] for a recent experimental analysis. Charm quarks can indeed be produced not only by hard scattering at the early stage of the collision (prompt production), but also later on by collisions inside the QGP (non-prompt production) and their density is sufficiently high that they can recombine in charmonium states. Bottomonium states give a much clearer signal, see e.g. the discussion in [8], because these states are very massive and can only be promptly produced by hard scatterings: the probability of generating these states within the QGP is extremely low. Moreover, since bottomonia are heavy and compact objects, they do not equilibrate with the expanding medium, and can be really considered as external probes.
Recent experimental results by the CMS Collaboration [9][10][11] and by the STAR Collaboration [12] indicate a clear suppression of Υ states, meaning that the nuclear modification factor R AA , expressing the ratio of the yields of a state in HIC with respect to p -p collisions, is less than 1. The R AA decreases with increasing centrality and/or with increasing p T and is higher for the fundamental state, meaning that it is less suppressed. Indeed, the experimental results of Pb -Pb collisions at √ S N N = 2.76 TeV [11] seem to indicate sequential suppression of bottomonium states [10][11][12], in particular, integrating over centrality it has been found that R AA (Υ(1s)) ≃ 0.6, R AA (Υ(2s)) ≃ 0.1 and R AA (Υ(3s)) < 0.1.
The effective field theory (EFT) techniques are very useful for the description of heavy quarkonia, because they are suited for handling systems with well separated energy scales. In the case of QQ in a thermal medium [13][14][15][16] two distinct kinds of scales appear, namely the non-relativistic scales and the thermal scales. The non-relativistic scales are given by the mass of the heavy quark, m Q , the typical momentum transfer, 1/r ∝ m Q α s (α s = g 2 /4π is the QCD coupling constant), and the binding energy, E ∝ m Q α 2 s [17]. We have assumed the weak coupling regime and identified the relative velocity between Q andQ with α s (see [18,19] for reviews). The relevant thermal scales to our analysis are the temperature T and the Debye mass m D ∝ gT . We shall discuss two possible hierarchies m Q ≫ 1/r ≫ T ≫ E ≫ m D and m Q ≫ T ≫ 1/r , m D ≫ E, which we shall refer to as Case I and Case II, respectively. If the bound state moves with respect to the medium, the EFT analysis becomes more complicated, because additional energy scales may appear [20]. We shall restrict our analysis here to the case of moderate velocities (v ≁ 1) for which no further scales are induced, so that Case I and Case II above can be safely addressed. However, at some instances we will push our results to the ultra-relativistic limit (v → 1). Although this gives the correct results in the QED case [20], and hence, we expect them to be sensible for QCD as well, one must keep in mind that they are on a less firm ground.
Using EFT techniques it has been shown that, at least in perturbation theory, the dissociation of heavy quarkonia is not due to the Debye screening but to the appearance of an imaginary part in the potential [13,14,21,22]. In other words, at high temperature heavy quarkonia disappear not because the binding energy vanishes, but because the thermal width becomes so large that the QQ state melts in the continuum. In QCD two different processes contribute to the thermal width: inelastic parton scattering, which is the dominant process for m D ≫ E, and the gluo-dissociation process that corresponds to the decay of a color singlet state into a color octet induced by a thermal gluon; this process is dominant for m D ≪ E [23,24] (see [25] for an early discussion). The inelastic parton scattering is often referred in the literature as Landau damping, the reason is that this scattering is always mediated by a spacelike gluon and can be related to the absorptive part of the gluon propagator. We shall also use this nomenclature from now on. In the strong coupling regime, the effect of an imaginary part in usual potential models has been addresed in [26], and in the so called T-matrix approach imaginary parts are incorporated in heavy quark self-energies through a set of Schwinger-Dyson equations, see for instance [27,28]. Recently, the imaginary part of the potential has also been calculated on the lattice [29,30] (see also [31,32], for a description in terms of open quantum systems).
The study of bound states propagating in the QGP at finite velocity is relevant for Υ states that are promptly produced in HIC and will cross the hot medium with a relative velocity v. In principle it might happen that heavy flavors are drifted by the expanding plasma. Indeed the PHENIX Collaboration [33,34] has observed a large v 2 of heavy-flavor electrons, suggesting that there is significant damping of heavy quarks while they travel across the medium. This picture has also received support from microscopic calculations of heavy quark diffusion in the quarkgluon plasma [35]. However, the elliptic flow of the Υ induced by the expanding medium should be negligible if the Debye length is larger than the typical distance between quarks, because heavy quarkonium at distances larger than its radius is colorless. Therefore, in both Case I and Case II the drift should be small and certainly less important for bottomonia than for lighter quarkonia, like the J/ψ, which can be non-promptly produced and are expected to roughly comove with the thermal bath. This is because before recombining both charm quarks have been drifted by the expanding QGP.
In the first study of moving QQs performed in [36], the hierarchy of scale of Case II was assumed, but only the real part of the potential was considered. The imaginary part of the potential was studied in QED in [20], where the velocity dependence of the cylindrically symmetric real and imaginary parts of the potential were determined. In the present paper we extend the analysis of [20] in two directions.
Regarding Case I, we consider QCD instead of QED; the main difference is that while in QED a proton and an electron will always form an electrically neutral state, in QCD a heavy quark and a heavy antiquark can be found in a singlet and an octet state, and this induces new terms in the computation. We determine the velocity dependence of the thermal width and of the energy shifts at the leading order. In particular, we find that at the leading order the energy shifts of the s-wave states do not depend on the velocity (like in QED), but the energy shifts of all the other states depend on the velocity (unlike in QED).
Regarding Case II, we extend the analysis of [20] by deriving an approximate analytical expression for the s-wave width as a function of the temperature and of v, valid for the particular hierarchy of scales T ≫ 1/r ≫ m D ≫ E. Moreover, considering the more general case, where T ≫ 1/r and m D ≫ E but the product rm D is arbitrary, we solve the corresponding Schrödinger equation numerically and determine the spectral representation of the two-point function.
This paper is organized as follows. In Sec. II we discuss Case I, corresponding to the hierarchy m Q ≫ 1/r ≫ T ≫ E ≫ m D . We derive the expression of the width and of the energy shifts as a function of the temperature and of the velocity of the bound state. In Sec. III we discuss Case II, corresponding to the hierarchy m Q ≫ T ≫ 1/r ∼ m D ≫ E, this section is divided in two subsections, in the first one we do an analytical analysis of the case 1/r ≫ m D while in Sec. III B we solve the Schrödinger equation numerically for 1/r ∼ m D and determine the spectral function for various values of T and v. In Sec. IV we present a brief discussion of the observable consequences of the velocity dependent thermal width, we compare our results with existing lattice simulations and we draw our conclusions. In the Appendix A we discuss the framework used to take into account the effect of a moving thermal medium. In the Appendix B we present some details and numerical checks of the procedure used in Sec. III B to derive the spectral amplitudes.

II. CASE I
The low-lying bottomonium states, Υ(1s) and η b , produced at early times in relativistic heavy-ion collisions are likely to have a typical size, r, smaller than the inverse temperature during most of their evolution in the QGP. At intermediate times, the temperature is also likely to be larger than the binding energy. Having in mind this possibility we shall study in detail the particular case This hierarchy of energy scales was considered in [15] for a thermal bath at rest. For a moving thermal bath, it was studied in full detail for the hydrogen atom in [20]; here we generalize those results to QCD. The general formalism to deal with a moving thermal medium is reviewed in the Appendix A. An important outcome is that in the bound state reference frame two additional energy scales should be considered: where v is the velocity of the medium with respect to the bound state. When v → 1 these scales are widely separated, a fact that must be taken into account in order to build the appropriate effective field theory (EFT). For instance, in [20] the appropriate EFT was constructed for the case T + ≫ 1/r ≫ T , which is different from the EFT obtained for the case T + ∼ T ∼ T − , valid for v ≁ 1. We shall mainly restrict ourselves to the latter case, and only comment on the limit v → 1. Note, indeed, that the QED analysis of [20] shows that the results obtained with the EFT theory valid for T + ∼ T ∼ T − and then naively extrapolated to the v → 1 case, coincide with the ones obtained with the proper EFT with T + ≫ 1/r ≫ T (if no large log resummations are performed). Hence, our results may hold for the v → 1 case as well.

A. Matching between pNRQCD and pNRQCDHTL
Since 1/r ≫ T we can take as the starting point the pNRQCD Lagrangian at T = 0 [37,38], which is obtained from QCD by sequentially integrating out energy scales of order m Q and of order 1/r, where L light is the QCD Lagrangian for light quarks, g is the coupling constant, E is the chromo-electric field, S and O are the quark-antiquark singlet and octet fields respectively, and (r = |r|) correspond to the singlet and octet Hamiltonians (the dots stand for 1/m Q corrections); hereafter C A = 3 and C F = 4/3. Thermal corrections to this Lagrangian are exponentially suppressed because the energy scales integrated out (m Q and 1/r) are much larger than T . Note also that no dependence on the velocity appears at this stage because the velocity enters in the calculation through the scales T + and T − in (2) only. Because we are assuming that the binding energy, E, is much smaller than the temperature, we may integrate out energy scales of the order of T as well. If we do so, we obtain an EFT which is temperature and velocity dependent. This EFT was called pNRQCD HTL in [15], where it was used in the case of vanishing velocity. We consider here the general case of non-vanishing velocity, following the analogous QED calculation developed in [20]. At the order we are considering, the pNRQCD HTL Lagrangian is obtained from L pNRQCD in Eq. (3) by replacing L light → L vHTL , and h s → h s + δV s , where L vHTL is the Hard Thermal Loop Lagrangian for a plasma moving with a velocity v [39], and the potential δV s encodes thermal contributions to the singlet potential, which depend on the velocity as well. The expression of δV s can be obtained by a standard matching procedure, using dimensional regularization (DR) to regulate the IR divergences arising from the expansion T ≫ E. In this case we have to consider the pNRQCD diagram in Fig. 1, where the dipole vertices and the octet propagator can be read off from (3), with D µν (k) the gluon propagator. Since T ≫ E and we use DR, the following expansion can be performed which corresponds to a temperature expansion, meaning that upon substituting this expression in Eq. (5) we can expand the thermal contribution of the singlet potential as follows where δV s,T n ∝ T n . In the Coulomb gauge where f B (x) is the Bose-Einstein distribution function, see the Appendix A. Notice that δV s,T 3 vanishes because the integrand is an odd function of k 0 ; also δV s,T vanishes because the integrand differs from an odd function by a δ(k 0 ), which forces the integral to be zero in DR. Hence, the only non-vanishing contribution is given by where we use the following notation, already introduced in [20], We can manipulate r i (E − h 0 )r j in the same way as it was done in the v = 0 case in [15], obtaining a more compact expression where N c is the number of colors. The correction to the singlet potential in the pNRQCD HTL Lagrangian is where the O(α 2 s r 2 T 3 ) contributions above arise from α s corrections to the diagram in Fig. 1. They have been calculated in [15] for the v = 0 case. Differently from the hydrogen atom case [20], the correction to the potential depends explicitly on the velocity. This could be expected from the fact that the Gromes relation (which is deduced by assuming Poincaré invariance) is violated at finite temperature [40]. However, as we shall detail in the next section, for the s-wave states the corresponding velocity dependence in the energy shifts cancels out at first order in perturbation theory.

B. Computation in pNRQCDHTL and final results
With the obtained pNRQCD HTL Lagrangian we can evaluate the thermal corrections to the binding energy and to the decay width of the various hydrogen-like states. Since we are using perturbation theory, we shall assume that the wave-functions are given by the unperturbed hydrogen atom solutions, which can be identified by the principal quantum number, n, the angular momentum, l, and the magnetic quantum number, m. For a given heavy quarkonium state the binding energy at the leading order in the perturbative expansion is given by where E c n is the binding energy taking into account only the Coulombic part of the potential and Σ s is the self-energy of the singlet component of the heavy quarkonium. Clearly, the Coulombic part of the potential does not contribute to the decay width, which is nonzero only because of the thermal corrections, and at the leading order in the perturbative expansion it is given by The singlet self-energy can be determined computing the diagram in Fig. 1 but this time in pNRQCD HTL . In order to properly take into account the moving thermal bath, the boosted Bose-Einstein distribution function has to be used, and since T ≫ E, we expand hence the self-energy can be written as follows This integral is very similar to the one evaluated in the QED case [20] and can be written as follows where ℜJ ij is given in Eq. (48) of [20], By manipulating r i (E − h o ) 2 r j we obtain the following result and at the considered order it is pure imaginary. On the other hand, the correction to the singlet potential in Eq. (15) is real, with no imaginary part. Therefore, at the leading order The correction to the binding energy in the regime 1/r ≫ T ≫ E ≫ m D is given by where l ′′ l ′ m ′′ m ′ |lm are the Clebsch-Gordan coefficients and is the expectation value of the radial position operator in the hydrogen atom, with a 0 = 1/m Q C F α s the Bohr radius. In the expression of the binding energy we can distinguish three different contributions, corresponding to the three terms in the square bracket of Eq. (24). The first one is an overall energy shift, independent of the quantum state. The second term is a shift of the binding energy that removes the degeneracy in l associated to the so-called "accidental" symmetry of the hydrogen atom. The third term in the square bracket is an energy shift that depends not only on n and l, but also on m and it is thereby related to the breaking of rotational invariance. There exists a privileged direction corresponding to v, thus the binding energy depends on the relative orientation between the angular momentum and the velocity. For s-wave states, having zero angular momentum, there is no such dependence and this contribution to the binding energy vanishes, and as already anticipated it does not depend on the velocity of the plasma. The latter result is surprising, because one would have naively expected that a v dependence should arise because for the moving bound state the effective temperature depends on v, see the Appendix A. However, our calculation shows that this is not the case, at least at the first order in perturbation theory. Regarding the width, from the expression of the self energy in Eq. (22), we obtain where is a negative and decreasing function of v. It can be easily shown, using the expression above, that for any state the width is a decreasing function of the velocity, vanishing for v → 1, meaning that an ultra-relativistic velocity has the effect of stabilizing the system. This behavior is due to the √ 1 − v 2 prefactor in Eq. (27), which can be traced back to the expansion of the boosted Bose-Einstein distribution function in Eq. (18) and is therefore due to the "Doppler shift" of the temperature, see the Appendix A. As we shall see in the next section, an analogous behavior is obtained for the hierarchy of energy scales of Case II in the ultra-relativistic limit, although the microscopic description appears to be different.
In the expression of the width we can further distinguish two contributions, corresponding to the two terms in the square bracket in Eq. (27). Both are velocity dependent, but the first one does only depend on n, meaning that it originates from terms that do not break the rotational and the accidental symmetries. The second term depends on all the quantum numbers and vanishes for s-wave states. Thus, for s-wave states the width simplifies to which, as observed above for the general case, is a decreasing function of the velocity, vanishing for v → 1.

III. CASE II
The dissociation of heavy quarkonium is expected to occur for as it happens in a thermal bath at rest. In the color-screening model introduced in [1] the dissociation takes place because the number of bound states supported by a Yukawa potential decreases with the range of the potential, which is proportional to the screening length (1/m D ). The effect of screening becomes important when the screening length is of the order of the size of the system (1/m D ∼ r). However in the actual real-time potential computed in [21] (and confirmed by the EFT computations [13,14]) the dissociation takes place because the potential develops an imaginary part (Landau damping), and bound states turn into wide resonances as the temperature increases. This effect becomes important at a parametrically different scale, 1/(T m 2 D ) 1/3 ∼ r (up to logarithms) [13]. It is then particularly interesting to address the question whether the mechanism of dissociation (screening versus Landau damping) remains the same when the bound state moves with respect to the thermal bath.
An EFT study of this situation in the QED case for muonic hydrogen submerged in a bath of massless electrons was already performed in [20](see [16] for the v = 0 case). For heavy quarkonium the results are analogous and can be obtained by changing the value of the Debye mass from the QED to the QCD value and by correcting for trivial color factors, as it was already pointed out in [20]. We briefly review the basic steps of the derivation below.
• Since m Q ≫ T , the starting point of the calculation can be the NRQCD Lagrangian at zero temperature [17].
By integrating out the temperature scale we arrive at the NRQCD HTL , an EFT whose Lagrangian is the sum of NRQCD for the heavy quark sector (with thermal, velocity independent, corrections to the heavy quark mass) [13], and the HTL Lagrangian for gluons and light quarks, which now depends on the relative velocity v between the thermal bath and the bound state [39].
• We can also integrate out the scales 1/r and m D , which leads to pNRQCD HTL . Since gluons and light quarks in the HTL Lagrangian develop a mass gap of the order of m D , this effective theory does not contain them as explicit degrees of freedom, and hence it reduces to a singlet heavy quark-antiquark field interacting through a potential V s (which depends on v as well). The main conceptual difference with respect to the case discussed in the previous section is that now, in general, the thermal contributions cannot be considered as a perturbation in the potential.
The potential V s coincides with the one that was computed numerically in Sec.V of [20]. In that paper qualitative arguments were put forward on how the dissociation mechanism is modified when the velocity of the bound state with respect to the plasma increases. We will quantify those arguments here, by focusing on the effects of this potential on the physics of the 1s state. We shall discuss two different cases.
1. We consider the particular case 1/r ≫ m D . The thermal contributions can still be considered as a perturbation to the Coulomb potential. This allows us to compute the leading thermal corrections to the decay width (almost) analytically and to derive some explicit expressions for the velocity dependence. In particular we can parametrically estimate how the dissociation temperature depends on the velocity if 1 − v ≫ m D a 0 . Then, as in the v = 0 case, the dissociation mechanism is dominated by Landau damping effects. This has to be contrasted to Case I in Sec. II where the decay width is entirely due to gluo-dissociation.
2. We consider the general case in which the relative size between 1/r and m D is left arbitrary (1/r ∼ m D ) and compute the spectral function. Although the concept of dissociation temperature is useful for qualitative estimates, there is no universal definition for it, and hence it is of limited usefulness for a quantitative comparison of our results with other approaches. On the contrary, the spectral function is a well defined quantity so that our results can be straightforwardly compared with those obtained by different methods, in particular by lattice computations. Furthermore, it is related to a physical observable, the thermal dilepton production rate [41,42]. In current HIC experiments, however, the heavy quarkonium states are not expected to be thermalized, but rather to act as hard probes of the medium, and hence the connection of the spectral function to the dilepton spectrum in this case is not straightforward. In the spectral function a bound state with zero decay width appears as a delta function whereas scattering states produce a smooth curve. The spectral function allows us to observe all the intermediate situations which occur when changing the thermal bath temperature and velocity.
For the expressions of the coupling constant and the Debye mass, we will use the following parameterization (we set N c = 3, N f = 3), where the MS renormalization scheme has been adopted with Λ MS = 250 MeV; we also fix m Q = 4.881 GeV and the Bohr radius of Υ(1s) is given by a 0 ≃ 0.74 GeV −1 , both values are taken from [43]. This choice is motivated from the fact that computing higher order corrections to the potential would introduce a dependence on the renormalization scale of the type log n (rµ); on the other hand, computing higher order corrections to the Debye mass would introduce terms proportional to log n (2πT /µ).
A. The particular case T ≫ 1/r ≫ mD ≫ E In this case the potential can be considered as the Coulomb potential plus a perturbation, and hence the following formula provides a good approximation to the decay width of a s-wave state where ψ n (r) is the wave-function for a s-wave state in the Coulomb potential and ∆ S is the symmetric part of the 00 component of the gluon field propagator in the Coulomb gauge, which has been computed in [20] for QED. Its generalization to QCD can be straightforwardly obtained by introducing a color factor C F and substituting the value of m D by the corresponding QCD one, and depends on v, k = |k| and on θ, the angle between the vectors k and v. In the above equation and we have made the dependence on the Debye mass explicit by defining and Π R (z, v) is the retarded self-energy of the 00 component of the gluon field in the Coulomb gauge with k 0 = 0, which was first computed in [36] 1 . In the reference frame where the bound state is at rest with and In principle we can obtain the decay width for any s-wave state, but for illustrative purposes we shall focus on the ground state (n = 1). It is convenient to start the computation in Eq. (32) by first performing the integration over r, where we have switched to cylindrical coordinates in momentum space and performed as well the trivial integration over the azimuthal angle. The above expression can be numerically integrated, but we first obtain an approximate expression valid at moderate velocities. In this case, for any angle, g(z, v) is of order 1 and there are only two scales in the previous integral, 1/a 0 and m D . Moreover they fulfill the relation 1/a 0 ≫ m D , so that the technique of threshold expansion [44] can be used to work out the integral, thus obtaining This equation can be further simplified by taking into account that the log 2 mD a0 is logarithmically bigger than the rest of the terms in the parenthesis. With this approximation we arrive at the following result, or, equivalently, which holds up to O(1/ log(m D a 0 )) accuracy and is independent of the heavy quark mass and of the temperature. Then, in the regime T ≫ 1/r ≫ m D , the decay width increases with the velocity, as far as it remains moderate (v ≁ 1). Note that this behavior is opposite to the one observed in the regime 1/r ≫ T ≫ E ≫ m D in Sec. II B, see Eq. (29). If we take into account that m 2 D is proportional to T 2 , the decay width at temperature T and velocity v is the same as the one that we would observe at v = 0 but with Results beyond the logarithmic accuracy of (42) can be obtained by evaluating (40), or even better (39). The numerical values of Γ(v)/Γ(0) for the Υ(1s) state are reported in the left panel of Fig. 2, for three different temperatures, together with the approximate expression (42), solid black line, for 0 ≤ v 1. The approximate expression correctly reproduces the numerical values for v 0.5, but for larger values of v the ratio of the width decreases and becomes temperature dependent, departing from Eq.(42). This is due to the fact that for v → 1 further scales must be considered (T + ≫ T ≫ T − ) and (40) does not hold. This expression relies on the fact that 1/a 2 0 ≫ m 2 D |g(z, v)|, which not always holds for v → 1, even if 1/a 2 0 ≫ m 2 D does. In order to ascertain the reliability of the expression in (39) and the origin of the difference between (39) and (42), let us scrutinize the velocity and angular dependence of m 2 D g(z, v). We can interpret the square root of m 2 D g(z, v) with positive real part [20,36] as the velocity dependent Debye mass m D (v, θ) (m D (v, θ) should not be mistaken for the parameter m D that we used before, they coincide at v = 0 only and have the same size for moderate velocities v ≁ 1 only). In Fig. 3 we present the plots of ℜ[m D (v, θ)]/m D and of ℑ[m D (v, θ)]/m D as a function of θ for v = 0.1, 0.5, 0.9, 0.99. The real part is peaked at θ = π/2, corresponding to a vanishing value of the imaginary part, which is instead peaked at a value of θ that with increasing v approaches π/2. For θ ≁ π/2, m D (v, θ) is small and the imaginary part is of O(m D ) for any value of v, meaning that the bound state can be approximated with a Coulombic wave function, and therefore in this region (32) represents a good approximation. Moreover, for v 0.5 both the real and the imaginary part of m 2 D g(z, v) are of O(m D ), irrespective of the value of θ, and therefore the approximate expression (42) is reliable. This approximation is still qualitatively good up to v ≃ 0.9, although the increased value of m D (v, θ) for θ ∼ π/2, suggests that the quantitative agreement might be lost, as indeed can be observed in the left panel of Fig. 2.
An angular region that may jeopardize the perturbative expansion about the Coulombic wave function is only present for v > 0.9, around θ ∼ π/2. Indeed for θ = π/2 the real part of the Debye mass has a peak, and for θ ∼ π/2 the imaginary part is large. However, this angular region is small. In order to quantify this region we consider θ)], which has two maxima for and therefore for v > 0.9 the angular region around θ = π/2 where the real and the imaginary parts are large is given by which clearly shrinks to zero for v → 1. In order to clarify that the contribution of this angular region is small, we plot in the right panel of Fig. 2 the thermal correction to the Coulomb binding energy 100|ℜδV (r)|100 , δV being the thermal contribution to the singlet potential, normalized to the Coulomb binding energy. This quantity should be small for (39) to be reliable, as it turns out to be the case. Note, however, that the angular region ∆θ gives for v → 1 the largest contribution to (39), and in this case , the approximate expression in (40) does not hold anymore, in agreement with the results in the left panel of Fig. 2. We find in this case, using the same techniques of integration by regions, that for v → 1 the decay width goes to zero like α s T √ 1 − v 2 whereas the energy shift goes to a constant, consistent with the results reported in the right panel of Fig. 2.
The velocity dependence of the width in the ultra-relativistic limit is similar to the one discussed in the Case I, in Sec.II B, although the microscopic mechanism is different, Landau damping in the present case, and gluo-dissociation in the former. The reason for the decrease in the decay width is probably related to the fact that a moving bound state feels a plasma with a non-isotropic effective temperature see Appendix A for more details. Actually, the effective temperature is higher than T in the forward direction, and lower than T in the backward direction, thus it is not obvious that the width of the moving bound state should increase -or be modified at all -when the bound state moves with respect to the thermal medium. However, in the ultra-relativistic case the effective temperature is almost everywhere less than T , see Fig. 7, and it is higher than T only in a narrow region θ ∼ 0. According to the previous discussion this angular region does not give the leading  contribution to the width, which is instead dominated by the θ ∼ π/2 region, where the heat bath is effectively cold. Thus, a velocity close to 1 tends to stabilize the system. The presence of an imaginary Debye mass for any v > 0 can be related to the collisionless transfer of energy between the heavy quarks and the gauge fields. An accurate description of this phenomenon would require the discussion of the propagating modes, but an imaginary part of the Debye mass does in any case signal an instability. This phenomenon is akin to the plasma instabilities generated by a charged current in a plasma. A similar result was indeed obtained in [45,46] where the destabilizing effect of a single heavy quark propagating in a thermalized QGP was studied.

B.
The general case T ≫ 1/r , mD ≫ E: the spectral function The s-wave spectral function was computed in [47,48] for the case at rest (v = 0). The procedure developed in Ref. [47] can be easily generalized to a moving bound state. We shall use the expression of the potential determined in [20], which takes into account the relative velocity, v, between the bound state and the expanding plasma, and consider that the system has cylindrical symmetry, with its symmetry axis in the direction of v.
The formalism introduced in [47] can be generalized to cylindrical coordinates, resulting in the following expression of the spectral function where ψ(t, r, z) = u(t, r, z)/r, and u(t, r, z) is the solution of the Schrödinger equation with the initial condition u(0, r, z) = −6N c rδ 2 (r)δ(z). In order to numerically handle the two-dimensional Schrödinger equation, we use an operator-splitting method, namely we split the differential equation in two differential equations each containing derivatives with respect to only one variable (z or r). Therefore, the whole Hamiltonian is divided in two pieces, H = H 1 + H 2 , where and we then solve the corresponding Schrödinger equation recursively in a discrete space-time, see the Appendix B for more details and for a check of the numerical code. In Fig. 4 we report the spectral functions at vanishing velocity for certain values of the temperature as a function of ω/m Q , where ω ≪ m Q is the non-relativistic energy (ω = 0 corresponds to a relativistic energy of 2m Q ). For the sake of comparison, in the right panel of Fig. 4 we also report the spectral functions obtained in [47] with a different choice of α s and m D , see Refs. [47][48][49] for more details. At T = 250 MeV the spectral function is given by the superposition of a peak, corresponding to the Υ(1s) bound state and of a continuum. The bound state has a thermal width which is determined by the imaginary part of the potential and is dominated by the Landau damping. The width of the bound state increases with increasing temperature, and correspondingly, the contribution of the continuum increases. At T ≃ 500 MeV no peak of the spectral function is visible, meaning that the bound state has dissolved into the continuum.
With increasing temperature the position of the peak slightly drifts away to the left, meaning that the Υ(1s) mass decreases. The binding energy decreases as well, since the absolute value of the real part of the potential for r, z → ∞ increases. In fact, if one subtracts such an asymptotic value from ω, in both figures the peak drifts to the right as the temperature increases. However, the bound state does not disappear because the binding energy vanishes, but because for T ≥ 400 MeV Landau damping prevents the formation of a bound state. Then, we consider the effect of a non-vanishing velocity. In Fig. 5 we report the spectral functions for the Υ(1s) at T = 250 MeV (left panel) and at T = 400 MeV (right panel) for a few values of the velocity of the plasma. Comparing Fig. 4 and Fig. 5 it is apparent that the effect of an increasing velocity -at least qualitatively -is akin to the effect of an increasing temperature. Actually, in this case the position of the peak of the spectral function at T = 250 MeV does not seem to change at all. But the main effect is that with increasing velocity the height of the peak decreases and the corresponding width increases; a behavior that emulates an increase of the temperature of the medium. At T = 400 MeV and v = 0, the spectral function has a small peak, which almost disappears at v ≃ 0.9. The result at this temperature is qualitatively similar to the one observed in [50] at their highest temperature, although one has to take into account that we use different reference frames.
It is interesting to notice that the tendency of the peak to become smaller and wider as the velocity increases, changes when going from v = 0.9 to v = 0.99. At v = 0.99 the Υ(1s) peak is slightly higher and slightly narrower than at v = 0.9, as shown for two different temperatures in Fig. 5. This behavior is consistent with the result, already discussed in the previous sections, that the bound states become stable at ultra-relativistic velocities. As in the particular case considered in the previous section, this behavior can be related to the fact that the effective temperature of the plasma is the one given in Eq. (46), and therefore for large v the plasma is almost everywhere cold. The fact that m D (z, v) becomes purely imaginary, implies that the potential ceases to be Yukawa-like and becomes oscillatory, as already observed numerically in Ref. [20].
Except for this peculiar behavior at v → 1, both the spectral function analysis and the computation of the width through (32) show that the width increases as the velocity of the plasma increases, as far as v 0.9. This is just the opposite of the results of Eq. (29) in Case I. The reason is that the two results refer to different energy regions, which are dominated by different processes. In Case I the thermal width is dominated by gluo-dissociation processes. In the present case, the dominant contribution is determined by Landau damping, which is a collisionless process. We shall further comment on this issue in the Sec. IV.

IV. DISCUSSION AND CONCLUSIONS
In this section we first discuss how the relative velocity v used throughout is related to measurable quantities in HIC experiments, like the momentum of the heavy quarkonium state in the lab frame, P µ , and the local velocity of the QGP, w, in that frame, and make a rough estimate of the importance of the relative motion in the yields. Next we compare our result with lattice computations, earlier weak-coupling analysis, and AdS/CFT calculations. We close it with the conclusions.
The clearest experimental signal of the velocity dependence in the in-medium heavy quarkonium properties should be in the dilepton yields at fixed rapidity and transverse momentum. In order to have an estimate of the effect, we assume that in a central collision the produced medium expands at a constant velocity, w, with respect to the lab frame. Typical values for w quoted in the literature are w ∼ 1 and w ⊥ ∼ 0.6 for RHIC and w ⊥ ∼ 0.66 for LHC. We further assume that the system has had enough time to thermalize and that it is isotropic. A heavy quarkonium produced with a certain P µ in that frame, moves with respect to the plasma with a velocity v = −P 0 w + P ·w which is the velocity appearing in the formulas of the previous sections. Notice that for a given longitudinal momentum this velocity is not totally fixed, it still depends on the transverse momentum and on the modulus of the parallel and perpendicular velocities of the plasma in the lab frame, and on the angle ϕ between w ⊥ and P ⊥ in the transverse plain. The modulus of the velocity can be written as where M is the heavy quarkonium mass. Assuming a uniform distribution for the angle ϕ, the modification of the dilepton yields can be estimated by the following formula where Γ(v) is the velocity dependent decay width calculated in the previous sections evaluated for v ≡ v(ω, P ) given in Eq. (52), τ is the lifetime of the QGP, about 10 fm/c for RHIC or about 15 fm/c for LHC, and Y (v) stands for the yields. In order to estimate the size of the effect, we display in Fig. 6 the results obtained for P = 0 and w ∼ w ⊥ ∼ 0.66 for two different temperatures. We plot the ratio between the velocity dependent yield and the yield at v = 0 as a function of P ⊥ . The yield has a non-trivial dependence on the transverse momentum, which modifies with increasing temperature from a monotonic decreasing behavior at T 250 MeV to a non-monotonic behavior at higher temperature. A number of oversimplifications have been employed in Eq. (53). We have assumed that the heavy quarkonium decays in the medium and that the medium temperature and the expanding velocity are constant. These approximations should be reasonable if the decay is much shorter than τ , and since 1/Γ(v) ≃ 3 fm for Υ(1s), this seems the case. In principle these aspects can be corrected for along the approaches of [51,52] or [53] (see also [54][55][56]). We have as well assumed a constant P for the whole evolution, that is we have neglected the damping of the heavy quarkonium, which should be a reasonable approximation because the drift by the expanding medium is expected to be small. We have also ignored the velocity dependence in the production mechanism and the contribution of the continuum to the yield. Eq. (53) is a reasonable approximation if all the above-mentioned corrections factor out in the yield; which might not be the case. However, we postpone to future work a more reliable estimate of this quantity. In any case, we believe that Eq. (53), together with Fig. 6, is enough to pinpoint the importance of the relative velocity between the heavy quarkonium states and the thermal medium.
The spectral function of the bottomonium states in a moving thermal bath have been studied with different lattice methods in [50,[57][58][59]. It is important to take into account that while our computation is performed in the heavy quarkonium rest frame, lattice computations are done in the thermal bath rest frame. The spectral function in the heavy quarkonium rest frame ρ HQ (k 0 ) where k 0 = p 0 − M and p 2 ∼ M 2 , is related to the spectral function in the plasma rest frame ρ plasma by the following equation where γ = 1/ √ 1 − v 2 and η(v) is a function of the velocity that is not important for the discussion below. If the thermal modifications are a perturbation ρ HQ (k 0 ) in the vicinity of a peak is well approximated by a Breit-Wigner distribution hence   figure 4 of [59] for the ηc state with the prediction of our Eq. (42). The first column p is the momentum in the units used in [59] (∼0.5 GeV). The second column is the velocity of the plasma deduced by looking at figure 5 in the same reference. The third column is the width obtained by assuming that the spectral function can be approximated by a Breit-Wigner distribution and comparing the highest point of the peak with the points where the value is half the maximum. Finally the fourth column is the value of the width predicted by using Eq. (42), where, for Γ s−wave

1
(v = 0), we have used the corresponding value in the third column (106 MeV).
This means that even if heavy quarkonium is not modified by the velocity of the plasma in the frame where it is at rest, we would still see a modification of the spectral function in the plasma rest frame. This modification will lead to an increase of the energy where the peak is located and a broadening of the peak.
In the case of Sec. II, we can compare with the lattice results of [58]. Unfortunately, the velocity range explored in [58] was at most of the order of v = 0.2, for which no velocity dependent width change has been observed. This null result is compatible with our results taking into account the error in the lattice computations and the fact that their temperatures are not very high. Note that the ansätze made in [58] for the binding energy and the decay width as a function of the velocity, based on the hydrogen atom computation in [20], holds for s-wave states only according to our results for QCD. At zero velocity the lattice results of the same group [60,61] turned out to be compatible with our results and with those of [15]. The results of this section for the decay width also appear to be compatible with the weak coupling calculation of Ref. [62] at leading order (LO), which also shows a decreasing behavior of the decay width with the velocity for small screening masses. At that order, only gluo-dissociation diagrams contribute, like in our case.
In the case of Sec. III A, the results only hold if the thermal corrections can be considered as a perturbation. This implies that we can only compare to spectral functions that have approximately a Breit-Wigner form. Because of this we cannot compare with the lattice results in [50] but we can compare with those in [59]. By analyzing Figs. 4 and 5 in [59] we can obtain approximate values of the decay width of the η c for several momenta. The decay widths we obtain from those figures and the corresponding prediction from Eq. (42), which is flavor independent and so holds also for charmonium states, are shown in Table I. We observe a similar qualitative behavior, since in both cases the width is a non-monotonic function of the velocity, i.e. it increases for low values of v and decreases for larger v, but the value of v at which it starts to decrease is lower in [59] than in our case. However, one has to consider that a more detailed statistical analysis would be necessary to disentangle possible MEM artifacts from the actual width in the plots of Ref. [59], see [63] 2 . Moreover, the temperatures at which this comparison is done may be too close to the deconfinement phase transition for the weak coupling expansion used in our computations to be reliable. The equivalent temperature regime for bottomonium (i.e. higher temperature) would be much safer. The results of this section agree with the weak coupling estimate of the dissociation temperature in Ref. [64] and also appear to be compatible with the contributions to the decay width at next-to-leading order (NLO) displayed in Ref. [62]. In the last reference, it was found that the NLO contribution was much larger than LO one 3 . In our EFT approach this can be easily understood if the system is in the kinematical regime of Sec. III, in which the gluo-dissociation processes contributing to their LO are parametrically suppressed. This is also consistent with the arguments and results presented in Ref. [65].
A number of analysis on the velocity dependence of the screening length have been carried out for strong coupling using the AdS/CFT approach [66][67][68][69][70][71][72] (see [73] for a review). It is not straightforward to compare these results to ours, as they do not obtain an imaginary part in the potential. This would be a first important difference. Furthermore, in momentum space, what plays the role of the screening mass for us is the complex, angle and velocity dependent, Debye mass m D (v, θ), see Fig.3 , which translates into a non-trivial potential in coordinate space for which no simple analytical form has been found. Hence we cannot make further statements on this respect than those already made in Ref. [20] 4 . However, we can certainly compare with the two AdS/CFT calculations of the heavy quarkonium spectral function at non-vanishing velocity we are aware of [74,75]. These spectral functions qualitatively agree with ours in Case II at moderate velocities, in the sense that the bound state peaks become smaller and wider as the velocity increases. Let us finally remark that we observe in the ultrarelativistic limit an oscillatory behavior of the potential rather than an exponential damping, that would lead to the stabilization of the bound states, which is not observed in the AdS/CFT approach.
In summary, we have analyzed heavy quarkonium states moving in a weakly coupled QCD plasma. In the Case I, corresponding to the hierarchy m Q ≫ 1/r ≫ T ≫ E ≫ m D , we have found that the thermal decay width decreases as the velocity increases, like in QED [20]. The decay width is in this case dominated by gluo-dissociation processes [23]. However, unlike in QED, the thermal energy shift becomes velocity dependent, except for the s-wave states. In the Case II, corresponding to the hierarchy m Q ≫ T ≫ 1/r , m D ≫ E, we have found a different behavior for the decay width, namely it increases as the velocity increases, except for ultra-relativistic velocities for which it starts decreasing again. This non-trivial behavior was overlooked in Ref. [20]. The decay width is in this case dominated by the Landau damping. Putting all together, we conclude that the decay width depends in a nontrivial way on the temperature and on the velocity, which complicates the interpretation of HIC experimental data, as we tried to illustrate by Fig. 6. Our results are consistent and in qualitative and semi-quantitative agreement with the few available lattice data [58,59], and also appear to be compatible with the weak coupling analysis of refs. [62,64].
where we have defined the effective temperature which is plotted in Fig. 7 for few values of v. Fig. 7 helps to clear away the misconception that a bound state moving with non-vanishing velocity in a thermal bath feels a higher temperature. Indeed, the effective temperature is in most of the directions smaller than T ; for v ∼ 1 we find that T eff (θ, v) > T only for 0 < θ < √ 2(1 − v 2 ) 1/4 . Intuitively, the dependence of the effective temperature on v and θ can be understood as a Doppler effect. While at v = 0 it is clear that the thermal medium introduces a new scale T in the problem, it is not clear a priori how many scales a moving thermal medium introduces. This can be understood by using light-cone coordinates. We choose v in the z direction and define k + = k 0 + k 3 and Then, we have that β µ k µ = 1 2 where Therefore, in light-cone coordinates, it becomes explicit that the distribution function actually depends on two scales, T + and T − . Obviously, for any value of v, one has that T + ≥ T ≥ T − , and moreover T + corresponds to the highest temperature measurable by the observer, while T − corresponds to the lowest temperature measurable by the observer. In this work we consider always that T + and T − are of the same order of magnitude. Even though this is not so for very large velocities, in all the cases considered in [20] we found the results obtained assuming T + ∼ T − were indeed correct even for v → 1. With increasing values of the Debye mass the peak of the 1s state moves to higher values of energy, meaning that the corresponding binding energy decreases. Although the peak height decreases, for the same reason explained above, note that no appreciable broadening of the spectral function appears, meaning that the numerical procedure does not produce a fictitious increase of the width. Indeed, in this simple model (and in any model with a real potential) the dissociation happens when the peak of the spectral function approaches zero. Our numerical results indicate that the 1s state of the Yukawa potential dissociates at λ ≃ 1.2, in good agreement with the numerical results of [80,81].
Regarding the 2s state, at λ = 0.1 (blue line) it is still visible, with binding energy ω ≃ −0.05m e α 2 (in good agreement with the results of [80,81]), but for larger values of λ the 2s state is no more visible, although it is known that it only dissociates at λ ≃ 0.31, see [80,81]. The reason, as explained above, is that the corresponding peak is very small and cannot be identified with the used numerical accuracy.
In summary, from the analysis of the Yukawa potential, we conclude that the algorithm correctly reproduces the binding energy of the 1s state at any value of the Debye mass, but the decrease of the peak height observed in the left panel of Fig. 8 for increasing values of m D is an artifact due to the combined effect of the numerical discretization and of the reduction of the strength of the interaction channel. Remarkably, the algorithm does not produce a fictitious width. The analysis of the dissociation of the excited states for this model with the present method is problematic, because of the reduction of the peak height.
On the right panel of Fig. 8, we show the spectral functions obtained with the potential reported in [20] considering vanishing velocity. This potential has an imaginary component for any non-vanishing value of the temperature. At T = 0, i.e. for λ = 0, the potential is real and Coulombic, and the standard peaks of the hydrogen atom for the 1s state at ω ≃ −0.5m e α 2 , the 2s state at ω ≃ −0.125m e α 2 and the 3s state at ω ≃ −0.05m e α 2 are reproduced with a good accuracy. As before, the peaks of states with high principal quantum number are suppressed. Increasing the temperature, the binding energy of the 1s state decreases, but the corresponding spectral function not only moves to higher energies, it also becomes wider. From the insight gained in the analysis of the Yukawa potential, we conclude that the broadening of the peak is due to the imaginary part of the potential and not to the numerical procedure. Moreover, improving the discretization procedure, we checked that the reduction of the peak height of the 1s state (at nonvanishing temperature) is not an artifact, but it is instead a genuine effect related to the imaginary part of the potential. The reason is that with a finite imaginary potential, the spectral functions are not delta-functions, but smoother functions which can be resolved with the used discretization method.