Phenomenological study of quarkonium suppression and the impact of the energy gap between singlets and octets

The study of heavy quarkonium suppression in heavy-ion collisions represents an important source of information about the properties of the quark-gluon plasma produced in such collisions. In a previous paper, we have considered how to model the evolution of a quarkonium in such a way that the solution of the resulting equations evolves toward the correct thermal equilibrium distribution for an homogeneous and static medium. We found that it is crucial to take into account the energy gap between singlet and octet configurations when the temperature is not much greater than this energy gap. In this manuscript, we explore in more detail the phenomenological consequences of this observation in the more realistic situation of an expanding system. We consider two different scenarios, based on the same approximation scheme, but on different choices of parameters. In the first case, we rely on a Hard Thermal Loop approximation, while the second case is based on a recent determination of the static potential in lattice QCD. In both cases, we compute the decay width and the nuclear modification factor, both taking the energy gap into account and ignoring it. We find that the impact on the predictions is as large in the expanding medium as it is in the static case. Our conclusion is that this energy gap should be taken into account in phenomenological studies.


I. INTRODUCTION
Several physical phenomena are commonly invoked to explain the production of quarkonia in ultra-relativistic heavy ion collisions that are presently intensively studied at LHC [1][2][3]. Aside from the initial suggestion of the color screening of the binding potential by the quark-gluon plasma [4], collisions between the plasma constituents and the quarkonia could also lead to a suppression of the production rate. Various ways can be used to take these collisions into account. In the regime in which potential models are valid, it has been shown that collisions induce an imaginary part in the potential [5,6] (somewhat analogous to the imaginary potential of the nuclear optical potential model). In the effective theory picture where the interactions of quarkonia with the plasma are dominated by color dipolar interactions, collisions induce singlet to octet transitions [7]. In addition to these "suppression" mechanisms, there is also the possibility, when the number of heavy quarks created in a heavy-ion collision is large enough, that uncorrelated heavy quarks recombine to form a bound state inside the medium. This process, usually called recombination [8,9], appears to be an important contribution to charmonium production at present collider energies.
Recent reviews can be found in [23,24]. This formalism allows us to compute the evolution of bound states in a medium, taking into account quantum mechanical effects. In this way, it has been understood that quantum coherence [25] is essential at high temperatures when the binding energy is of the order of the decay width. At the same time, semi-classical equations can be derived rigorously from the open quantum system framework in special limits [13,18,[20][21][22]. This has improved the understanding of the validity region of these semi-classical approaches. The common assumption in most, if not all, approaches is that the medium is moderately affected by the presence of quarkonia. Various types of effective theories can then be obtained for the heavy quark dynamics, with further approximations becoming available, depending on the energy or time scales.
In the limit in which the temperature, T , is much bigger than the binding energy, E, the interaction between a quarkonium and the medium is Markovian. In this case, the dynamics is governed by a Gorini Kossakowski Sudarshan Lindblad (GKSL) equation [26,27]. Phenomenological predictions of quarkonium in heavy-ion collisions using the GKSL equation in the T ≫ E limit have been discussed in [14,15,28]. However, as was understood in [21], this limit leads to a maximization of the entropy without taking into account energy conservation constraints. In other words, this approximation does not lead to the correct thermalization. In [21], we showed that, when T ∼ E, the decay width is modified to include a dependence on the binding energy. With a slight abuse of language, we can say that the imaginary part of the potential depends on the energy. This dependence is a key element which guarantees that the evolution equation complies with the fluctuation-dissipation theorem. Phenomenologically, it is important to take this dependence into account when the temperature is not much larger than the binding energy, as was shown numerically in [29].
The aim of this manuscript is to emphasize the importance of the energy dependence of the decay width, not only in the case of a static medium, but also in a more realistic case of an expanding medium. We do this by computing the nuclear modification factor, R AA , in different scenarios. For each scenario, we compare the results obtained with and without the energy dependence of the imaginary part of the potential. We observe that, depending on the centrality of the collision, the productions of Υ(1S) and Υ(2S) are significantly less suppressed when we take into account the energy dependence. We emphasize that our purpose here is just to obtain an estimate of the effect, and show that it is indeed quite significant. We shall not go into a detailed analysis of experimental data.
The manuscript is organized as follows. In section II, we review the model that we developed in [21] and we discuss the imaginary part of the potential and its energy dependence.
In section III, we use this model to obtain phenomenological predictions of R AA using two different scenarios, one based on (resummed) perturbation theory, and another scenario in which we use lattice QCD inputs. Finally, we present our conclusions in section IV.

II. DISCUSSION OF THE MODEL
The model that we consider in the present paper is an extension of that introduced in the last section of [21]. Let us recall briefly its origin. One starts from the evolution equation for the reduced density matrix D Q of a pair of heavy quarks in a quark-gluon plasma, the interaction of the heavy quarks with the plasma, supposed to be weak, being limited to a one-gluon-exchange. The plasma is assumed to be in thermal equilibrium at a temperature T . This equation reads [21] Here, H Q is the Hamiltonian that governs the dynamics of the heavy quark pair in vacuum, with the corresponding evolution operator given by U Q (t) = e −iH Q t . n A x is the temporal component of the color current, whose explicit expression can be found in [20]. The interaction between the heavy quarks and the surrounding medium is captured by the correlators of thermal fluctuations of the gluon fields a A 0 , where D pl is the density matrix of the plasma, and A, B are color indices of the adjoint representation. We ignore here the color magnetic interactions, which are subleading for heavy quarks.
As was shown in detail in [21] this general equation can be simplified through a number of steps: • A part of the right-hand side of Eq. (1), which is hermitian, can be absorbed in a redefinition of H Q . This can be done by adding to both sides of Eq. (1) In the left-hand side, this is considered as a correction to the real part of the potential, while in the right-hand side, it remains as a compensating correction. From now on, H Q denotes the heavy quark hamiltonian after this redefinition.
• By projecting Eq. (1) on the eigenstates of H Q , with H Q |i = E i |i , one gets where D ij = i|D Q |j , E ij = E i − E j and L ij,kl D kl is a rewriting of the right-hand side of Eq. (1) from which the hermitian contribution to H Q has been subtracted.
• When the typical time between collisions of the quarkonium with plasma constituents is much bigger than the inverse of the binding energy, which we shall assume to be the case, L ij,kl can be treated as a perturbation. The multiple-scale analysis described in appendix B of [21] provides a way to handle the secular terms which appear in the perturbative expansion. Assuming for simplicity that there are no degenerate states, we obtain then an equation for the diagonal part of the density matrix, namely for the • A final simplification arises from the analysis of the color degrees of freedom. We are interested in the probability to find the quarkonium in a singlet state at the end of the evolution. By emitting or absorbing gluons, a color singlet can decay into an octet, and vice versa a singlet can emerge from an octet. We shall restrict ourselves to the dilute limit where there are only a few heavy quarks in the medium. This situation is appropriate for bottomonium production at the LHC. The study of charmonium would concern the dense limit with a large number of heavy quarks present in the medium. This is a much more difficult problem, which has only been addressed so far in special limits [20] (see also [13] for the analogous abelian case), or using the Boltzmann equation [18,19]. In the case of bottomonium, we have checked numerically that the contribution from octet decays is tiny. The reason is that octets are unbound and may occupy all the volume of the medium, while a bound color singlet occupies a small volume. Therefore, in the dilute limit, the rate of the octet to singlet transitions decreases as the volume of the medium increases. Besides, this transition is also suppressed in the large N c limit [30]. Finally, we arrive at a simple rate equation for the probability p s of finding the quarkonium in a singlet state. In the next subsection we analyse the properties of the decay rate Γ s .
A. Properties of Γ s .
The decay width Γ s is given by the following formula [21] where the operator sin q·r 2 , withr an operator acting on the heavy quark relative coordinates, describes the interaction of the quarkonium with the gluons of the plasma, after the center of mass of the quarkonium has been integrated out [21]. In the limit of small momentum transfer q, it reduces to a dipolar interaction, as we shall discuss shortly. In Eq. (7), |s is the wave function of the singlet and E s its binding energy. The ket |o, p is the state of a heavy quark pair in an octet state with energy E o p = p 2 M . In the large N c limit the octet potential is zero, therefore, the wave function of an octet can be approximated by a plane wave. As a consequence, the expectation value s| sin qr 2 |o, p is easily obtained as a Fourier transform of the wave function of the singlet bound state. The last ingredient needed to compute the decay width is the correlator ∆ > (w, q) (see Eq. (2)). We shall estimate it using the Hard Thermal Loop (HTL) approximation [31][32][33]. where and is the longitudinal component of the gluon polarization tensor. It is proportional to the Debye mass m D , given in leading order perturbation theory by where N c = 3, T F = 1 2 and N F is the number of flavours that can be considered light when computing m D (we take N F = 3). A plot of this Debye mass as a function of the temperature, as obtained with different scales for the running coupling constant α s , is given in Fig. 1. One can see that, for this range of temperatures, m D ≃ 2.5T , a result that is not truly in the perturbative regime. Thus we should regard our use of the HTL result as a convenient phenomenological guide, providing us with a reasonable approximation for the overall momentum and frequency dependence of the longitudinal tensor, while m D , to which it is proportional, could be viewed as a phenomenological parameter.
If one ignores the gap between singlets and octets, which may be legitimate when the temperature is high enough, one obtains the simpler expression which we refer to as the static limit [34]. This expression (12) can be read as twice the imaginary part of the potential, as originally computed in [5]. One may then interpret Eq. (7) in similar terms, as an energy dependent imaginary potential. This energy dependence, present in the exponential factor as well as in the correlator, plays a crucial role at low temperature since it ensures that the fluctuation-dissipation theorem is fulfilled. The impact of this energy dependence of the imaginary potential was studied numerically in [21], and is illustrated in Fig. 2. One can see that ignoring this energy dependence typically increases the decay width by a factor 2. The purpose of this note is to quantify this effect in the more realistic situation of an expanding plasma. We note that the importance of the energy gap was recognized in early studies of quarkonium interactions with matter [35,36].
Before closing this subsection, we note that the energy dependence in Eq. (7) originates from the energy gap between singlets and octets. We determine the singlet binding energy by solving the Schrödinger equation The potential V s is a screened potential which will be specified shortly. At this point we just comment on the interplay between screening and decay rate: • For simplicity, we set the origin of energies as that of the lowest energy octet state. Then the absolute value of the gap between a singlet state and the octet is equal to the binding energy of the singlet. Then, for the model to be valid, the condition |E s | ≫ Γ s must be fulfilled. This is equivalent to saying that the time scale for which thermal effects become important must be much larger than the inverse of typical energy difference between states. If this is not so, we cannot consider thermal effects as a perturbation using multiple-scale analysis (see discussion before Eq. (5)), which is crucial to get to the simple form of Eq. (6).
• The decay width increases rapidly as the gap decreases. Therefore, the phenomena Eq. (7) (with energy dependency) or from Eq. (12) (without energy dependency). Note that the effect of energy dependence is sizeable and amounts to typically a factor 2 reduction. Picture taken from [21].
of screening and decay due to the collisions influence each other. As the temperature increases, so does screening. Screening reduces the energy gap between singlets and octets, and this results in an increase of the decay width.
Because of this strong interplay, it is important to have a consistent and simultaneous treatment of both phenomena.

B. Relation with the dipole approach
At this point we wish to make a comment on the dipole approximation used in most EFT calculations based on the assumption of a clean separation of scales. In the (strict) limit This expression is plagued by an ultraviolet divergence whose origin can be found in the HTL approximation used for the gluon propagator. This approximation is only valid for q < ∼ m D , which is in principle compatible with our assumption. However, in expanding the sine function, assuming that qr ≪ 1, we have produced a factor q i q j whose magnitude is not controlled anymore by any factor. The HTL correlator results then integrated over momenta where it is not valid. This is the origin of the divergence. This divergence can be coped with by introducing a cutoff Λ in the q-integration, using the HTL approximation of the correlator for q < Λ and the full one-loop correlator for q > Λ. Since the divergence is logarithmic, the dependence on the cutoff eventually disappears, leaving contributions proportional to ln(T /m D ). We should note however that such a logarithm, while positive in the strict perturbative regime where T ≫ m D , may turn negative in the temperature range probed by present experiments (where the coupling constant may remain of order unity), leading to unphysical values for the decay rate.
An alternative to the procedure outlined above is to use a better approximation for the correlator ∆ > (0, q) in Eq. (14). We can rewrite this equation as [14,15] with κ a momentum diffusion coefficient [37]. This formula relates this diffusion coefficient to the correlator of electric field fluctuations, and this can be calculated beyond the HTL approximation, including possibly non perturbative effects [38]. A detailed discussion of these issues can be found in [23]. Here we note that in the important regime where the various scales are overlapping, a static diffusion constant is in any case not enough, since it ignores the frequency dependence which we claim is important.
We now present details of the implementations of the model and examine two scenarios to fix the parameters. For lack of a better terminology, we shall refer to these scenarios as the perturbative scenario and the lattice scenario. It should be understood however that the perturbative scenario departs from strict perturbation theory (as we have already alluded to) and the lattice scenario is not a lattice calculation.

C. The perturbative scenario
We obtain the binding energy and the wave-functions needed to compute Γ s and Γ 0 s by solving numerically the Schrödinger equation, using the algorithm described in [39]. We focus on Υ(1S) because it is only for this state that, within this scenario, the condition that the binding energy is much bigger than the decay width is fulfilled. This condition is needed for Eq. (6) to be valid. In the case of Υ(2S), the decay width is of the same magnitude as the binding energy at the lowest temperature that we consider. We need to specify first the hamiltonian H Q in Eq. (13). For M, we take as definition a naive version of the 1S of mass [40], which consists in setting M = We think that, given the accuracy we are aiming at, this is a reasonable option.
For the real part of the potential, we use a screened Yukawa potential where µ r and µ T are (independent) subtraction points. The scale µ r is related to the exchange of gluons between the heavy quark and the antiquark and enters the running coupling α s in Eqs. (7) and (16). A natural choice would be µ r ∼ 1 r , but this would lead to difficulties at large distance r in the numerical implementation of the running coupling.
Therefore we fix µ r ∼ 1 a 0 , where a 0 is the Bohr radius, which we define with the following self-consistent equation The scale µ T is the scale at which α s is evaluated when computing the Debye mass using Eq.  E s , that is when Γ s /E s becomes of order unity. We note that both < rm D > and Γ s /E s are increasing functions of the temperature. Fig. 3 also illustrates the influence of the gap by comparing the ratios Γ s /E and Γ 0 s /E s . We remind that the model is valid as long as this ratio is well below 1.
Consider now the uncertainty that results from the choices of µ r and µ T in the computation of the decay width. We only computed Γ s and Γ 0 s for a finite set of temperatures (increasing T from 200 to 450 MeV by steps of 50 MeV) to save on computational cost.  However, we have observed that the following function provides a good fit of our results, There is no theoretical reason for choosing this fit function apart from the fact that at very high energies one expects the decay to grow linearly with the temperature (neglecting the running of the coupling). This fit function will be used later when we calculate the R AA for the expanding system. Here this is used as an interpolation formula to estimate the error associated with the choices of the various subtraction scales. Figs. 4 and 5 illustrate the sensitivity of respectively Γ s and Γ 0 s to variations of µ r and µ T around nominal values, chosen to be respectively 1/a 0 and 2πT . Clearly the uncertainty associated with variations of µ T is bigger than that associated to changes in µ r . Overall, we observe that Γ 0 s > Γ s , as expected. As an alternative possibility to fix the basic ingredients of the model, we shall rely on the lattice data from a recent computation [41]. We refer to this second strategy as the "lattice scenario". As in the perturbative scenario, the binding energy and the wave function of the quarkonia are obtained by solving a Schrödinger equation. However, now we take as heavy quark mass the same as in [41], M = 4882 MeV, slightly larger than that of the perturbative scenario. As for the real part of the potential, it is given by where α, σ and m D are obtained from the fit performed in [41]. The values of m D differ from those in the perturbative scenario, and they are given at a different set of temperatures.
These values are listed in  [42]. This is why no result is given in The decay width is obtained from Eq. (7), with now the values of m D listed in Table I,  TABLE II. Parameters obtained by fitting the decay width in the three cases (from S1 to S3) and α s = α/C F with α taken from [41]. With these parameters, it becomes possible to study both Υ(1S) and Υ(2S) since for both states the condition E s > Γ s is fulfilled. The results of the calculation are displayed in Fig. 6 for the binding energy and in Figs. 7 and 8 for the decay widths Γ s and Γ 0 s respectively. In the case of the Υ(1S), the decay width is well fitted by the formula motivated in part by the fact that at high temperatures the decay width is linear with the temperature while at small temperatures we expect an exponential suppression due to the gap. For the case of Υ(2S), however, we found that Eq. (18) works better as a fit function.
These formulae are used in the analysis presented in Figs. 7 and 8. We note that the fact that we use the same fit function for Υ(2S) as we did for Υ(1S) in the perturbative scenario does not imply that our computation for Υ(2S) is perturbative. The biggest source of error in this lattice scenario comes from the value of m D . We consider three cases, referred to as S1, S2 and S3. The case S1 uses the central value of m D .
The case S2 and S3 correspond respectively to the the lowest (S2) or the largest (S3) values of m D that are compatible with the error given in [41]. The resulting fit parameters are given in table II. Using these parameters we obtain the decay widths of Υ(1S) and Υ(2S) shown in Fig. 7. Note that these figures contain an extrapolation to temperatures that are larger than those available in the lattice calculations (which are limited to T ≤ 260 MeV).
The results we have obtained for the binding energy and the decay width can be compared to recent results in the literature. The binding energies that we have obtained for Υ(1S) and Υ(2S) are shown in Figs. 6. For Υ(1S) these are qualitatively similar and of the same order of magnitude as those in [43][44][45][46][47]. For Υ(2S), our results are qualitatively similar to those in [43,44] but about five times smaller than those found in [46], in which the decay width is not based in the HTL model but rather on the computation in [48]. Concerning the decay width, we note that, in general, we obtain results for the Υ(1S) that are of the same order of magnitude as those in [43][44][45][46][47], however, the qualitative behavior is different since our results are more suppressed at small temperatures. The reason is that we take into account the suppression due to the energy gap between singlet and octet states, which was not taken into account in the other studies. We observe a similar behaviour in the case of Υ(2S), with the difference that the effect of the gap is much less pronounced here.
We have repeated the analysis for Γ 0 s (Fig. 8). The values that we obtained for Γ 0 s are substantially larger than those of Γ s , as it should, as long as we remain in the temperature region used in [41]. We observe, however, that at large temperatures Γ s in Fig. 7 is actually larger than Γ 0 s in Fig. 8. We view this as an artefact of extrapolating the results using Eq. (20) far beyond the region in which we performed the fit. This indicates that we have to be cautious when interpreting our results in the lattice scenario for collisions that probe high temperatures.

III. PHENOMENOLOGICAL IMPLICATIONS
We turn now to the main objective of this paper, which is to estimate the quantitative impact of the expansion on the previous results. To do so, we shall provide a crude estimate of the observable R AA . We assume that the quarkonium state, initially in a state n starts interacting with the plasma at some (small) finite time which we choose to be t 0 = 0.6 fm (value taken from [49]). The survival probability is then obtained from the rate equation with the initial condition p n (t 0 ) = 1. The survival probability of a given quarkonium state is , where t f is the time spent by the quarkonium in the plasma. Its value depends on the scenario considered and it will be specified shortly. The above calculation assumes that the temperature evolves slowly with time so that the interaction with the plasma can be treated in an adiabatic approximation. That is, one assumes that the quarkonium interacts with the medium as in the static case, with the medium characterized by its instantaneous temperature T (t).
For simplicity, we assume that the center of mass of the quarkonium does not move in the transverse plane. Thus the survival probability depends only on the time spent in the plasma, and on the local conditions in which it evolves [50]. These depend on the impact parameter, b, and on the point in the transverse plane, s, (with respect to the collision axis) in which the quarkonium is produced. We set the origin of coordinates in the transverse plane to coincide with the center of one of the nucleus. The initial temperature of the medium, denoted by T 0 (b, s), is related to the local energy density ε by the standard relation ε ∼ T 4 .
The energy density itself is taken to be proportional to the surface density of participants, which we estimate as a function of b and s using a Glauber model [51]. Thus we write where T A is the overlap function in the Glauber model, A = 207 and (in the type of collisions considered here) σ = 70 mb [52]. We take T 0 (0, 0) = 475 × 1.05 TeV for √ s = 5.02 TeV collisions at LHC. This number is obtained by using the standard value of [49] and taking into account that the temperature increases by about 5% when going from 2.76 TeV collisions to 5.02 TeV [53]. Regarding the time dependence of the temperature, we use Bjorken hydrodynamics, and ignore the transverse expansion, Now we have all the ingredients to compute R AA . According to the optical Glauber model, R AA for a given impact parameter is computed as where S(b, s) is the survival probability.

A. Perturbative scenario
The survival probability can be computed analytically if we use for the decay width the approximate expression (18). One gets then We can use this formula to compute R AA . In this case, we set T f = 200 MeV, arguing that a perturbative computation is only valid for temperatures above that of the phase transition.
We are aware that physics at lower temperatures might modify the survival probability.
However, we ignore those effects. At the moment, our aim is not to obtain a state-of-the-art phenomenological prediction, but to highlight the importance of the energy gap between singlets and octets.
In Fig. 9 we plot the results obtained by applying this formula, taking the dependence on µ T as a measure of our theoretical uncertainty. Although the uncertainty due to µ T is quite significant, we can clearly see the difference between Γ s and Γ 0 s . We see that taking into account the gap increases R AA by typically the same factor ∼ 2 as in the non expanding case. We note that in Fig. 9, and in the rest of the figures, we have plotted our results in the full range of centralities although our model is only valid when E ≫ Γ s .

B. The lattice scenario
The computation is completely analogous to what was already explained in the perturbative scenario. The only difference is that, in the case of Υ(1S) we used a fitting function of the type of Eq. (20) instead of Eq. (18). In the case of a decay width which follows Eq. (20) we can also obtain an analytic solution for S In the case of Υ(2S), we can also apply Eq. (25). Another difference with respect to the perturbative scenario is that now we use T f = 180 MeV. As we explained before, for temperatures smaller than this value, the fit obtained in [41] gives values for the energy gap between Υ(1S) and unbound states that are bigger than the experimental value of the dif- ference between the mass of Υ(1S) and two B mesons. Therefore, we assume that at such small temperatures there are no sizeable thermal effects.
Knowing this, we obtain the R AA results shown in Fig. 10. In the case of Υ(1S), the condition that the binding energy is much bigger than the decay width is fulfilled for a wide range of centralities. However, for Υ(2S) this is not the case. In fact, the validity range of our approximations for Υ(2S) in the lattice scenario is similar to the one of Υ(1S) in the perturbative case. In Fig. 11 we do the same but ignoring the gap. As in the perturbative case, we see that the difference is quite substantial (again a factor ∼ 2). Comparing Fig. 11 and Fig. 10 one can see that this factor can be essential to bring the R AA in agreement with its experimental value. This is approximately the case in Fig. 10, although no strong conclusion can be drawn from this remark, given the simplifications made in the present analysis. In fact, the results in Fig. 10 slightly underpredict the experimentally observed suppression, but this includes, among other things, a sizeable contribution from cold nuclear matter effects that are ignored in the present analysis.

IV. CONCLUSIONS
In this work, we have explored the phenomenological consequences of the observations made in [21], where we highlighted the importance of taking into account the energy gap between singlet and octet states when computing the decay width of a quarkonium bound state in a medium. The role of this energy gap can be understood as an energy dependence in the imaginary potential used to determine the bound state properties. The discussion in [21] was limited to the case of a static medium. In the present study, we have extended the analysis to the case of an expanding medium, and explored the predictions of a simple model in two different scenarios, one based on perturbative formulae, and a scenario based on the lattice QCD data of [41]. Our results have corroborated our initial findings, and indicate that the role of the energy dependence of the imaginary potential is as important in the expanding case as it is in the static medium: the value of R AA is significantly lower when the energy dependence is ignored. The effect is sizeable, the energy dependence reducing the expected suppression by typically a factor 2 for bottomonium at present LHC energies.

ACKNOWLEDGMENTS
This study was triggered by discussions initiated in the EMMI Rapid Reaction Task