Energy hierarchies governing quarkonium dynamics in Heavy Ion Collisions

In this paper, we critically examine hierarchies between energy scales that determine quarkonium dynamics in the quark gluon plasma. A particularly important role is played by the ratio of the binding energy of species ($E_b$) and the medium scales; temperature ($T$) and Debye mass ($m_D$). It is well known that if these ratios are much larger than one then the dominant process governing quarkonium evolution is dissociation by thermal gluons (gluo-dissociation). On the other hand, if this ratio is much smaller than one then quarkonium dynamics is dominated by screening and Landau damping of the exchanged gluons. Here we show that over most of the evolution, the scale hierarchies do not fall in either limit and one needs to use the full structure of the gluonic spectral function to follow the dynamics of the $Q\bar{Q}$ pair. This has a significant bearing when we follow the quantum dynamics of quarkonia in the medium. The inverse medium relaxation time is also $\sim T$ and if $E_b$ is comparable (or larger) in magnitude to $T$, the quantum evolution of $Q\bar{Q}$ is non-local in time within the Brownian approximation.


I. INTRODUCTION
In the vacuum, the bound states of heavy quarks (Q which can be b or c) and anti-quarks ( Q) feature three prominent momentum scales: the heavy quark masses M , the inverse relative separation 1  r , and the binding energies E b (see Ref. [1] for a comprehensive review).These scales satisfy the hierarchies M 1 r E b , which justifies non-relativistic treatments of these states.An additional relevant scale for their description is the scale of quantum Chromodynamics (Λ QCD ) which may or may not be somewhat smaller than E b [1] but can be assumed to be significantly smaller than 1/r (especially for bb pair).The hierarchy of these scales allows one to integrate out modes at the scale M , and 1/r systematically, and derive a low energy effective field theory (EFT) valid at the scale E b .At the lowest order in rE b , the EFT consists of nonrelativistic quarks bound by a potential [2].At higher order, the theory features interactions mediated by gluons of wavelength 1/E b .Effects of higher order terms are suppressed by positive powers of rE b , where factors of r can be seen as arising from a long wavelength expansion of the fields.This framework is called pNRQCD [1].
In a thermal medium at temperature T , new scales appear which govern the dynamic properties of quarkonia in the medium.It was pointed out in a classic paper [3] that the screening of the Q Q interaction on an inverse length scale m D could lead to the "melting" of the bound states.Moreover, later on, it was realized that scattering between bound states and the thermal constituents of the medium plays a major role in the dissociation of the quarkonium states.This leads to the generation of the imaginary part of the quarkonium potential in a thermal medium [4].Additionally, absorption of thermal gluons in the medium could lead to gluo-dissociation [5].It was shown [6] that pNRQCD naturally incorporates processes leading to gluo-dissociation [5,7] as its dynamical degree of freedom that includes low energy gluonic degrees of freedom (and other light degrees of freedom if any) in addition to the wavefunctions of Q Q pair.The corresponding emergent scale Γ ∼ 1/τ R (where τ R is the relaxation time of quarkonia) is related to the dynamics of inelastic interactions of QQ with the medium.Furthermore, E b and 1 r might themselves be modified from their vacuum values by these thermal effects.
Two medium scales that play a role in quarkonium dynamics in the quark gluon plasma (QGP) are T and m D .In the weak coupling limit, there is a hierarchy between m D and T [8].In this regime the coupling g is small and m D ∼ gT T .However, for the temperatures of interest (150 − 500 MeV), lattice results suggest that m D /T ∼ 2 [9].Using 2πT as the relevant energy scale [10] at which α s is computed also gives similar values of g.
This implies that in this regime leading order weak coupling expressions in g are not quantitatively reliable.For example, it is known that the higher order corrections to the momentum diffusion coefficient are larger than the leading order value [11].However, non-perturbative calculations of some relevant dynamical processes is still challenging and weak-coupling calculations are still useful.An important result in weak-coupling was obtained in Ref. [4] which showed that the potential between quark-antiquark pair is complex at finite T .
Such weak-coupling calculations have given insight into the problem and results from these calculations can be used to obtain estimates for experimental observables of interest: for example R AA in heavy ion collisions (HIC).
In this paper, we will focus on bottomonia and use leading order expressions for the gluon polarization tensor.But we will not assume m D /T is small.While this is not a formal expansion in g but might better capture some important qualitative dynamical properties of the QGP.This has been used in other papers for open heavy flavor [40].
The next question is how the thermal scales compare to the scales associated with bound states.Clearly, M T and a non-relativistic treatment is applicable for quarkonia slowly moving in the medium.For bottomonia, the values of 1/r are comparable to ∼ 1GeV [41] and we will assume that 1/r T and hence use pNRQCD to describe the system.However, we do not assume that the hierarchy between 1/r and the screening mass m D is so strong that we can ignore the screening of the Q Q potential when calculating quarkonium properties in the medium.In practice, we see that the effect of screening on the Υ(1S) wavefunction is small, but the Υ(2S) and Υ(3S) states are affected by screening.
On the other hand, there is no clear separation between E b and m D , T (FIG. 4 below).Moreover, their relative order depends on the species and can change as the medium cools down as it evolves.In this paper, we will take all three to be of the same order.Therefore, a non-relativistic treatment is still applicable.However, the integration of modes from 1/r to E b includes thermal effects.
Here we would like to point out that further assuming scale separations between E b , m D , and T can allow us to write simpler EFTs assuming specific choices of these hierarchies.These have been investigated in detail in a series of papers [6,[42][43][44].Our goal in this paper is to avoid assuming a clear separation between the three scales.In specific regimes where the separations exist, our results will clearly reduce to results from [42,44].However, we shall see that in a wide range of parameters, physics lies in an intermediate regime where clear separations do not exist.
To do this we use the full perturbative form of the gluon spectral function.We include contributions from the transverse and longitudinal gluons both in the Landau damping (LD) regime and in the space-like regime where gluo-dissociation occurs.This gives a clear framework to include both processes in a unified language and allows us to compare the contributions to decay from the various process.This is the first time the full gluonic spectral function applicable in both kinematic regimes has been used to compute the total decay rates.In our calculation, the singlet wavefunction is approximated to be the instantaneous eigenstate of a lattice inspired thermal potential and hence incorporates screening.For the octet state, the spectrum is fixed by the constraint that at large r the real part of the octet potential approaches the real part of the singlet potential.For the wavefunction, we systematically compare two limiting cases.One where the screening is strong that the potential is flat in r and the other where the screening is very weak.These can be seen as limiting cases of the physical situation where the screening length is comparable to that in the singlet channel [45].
The comparison between E b and T is shown in FIG. 4 which clarifies that these two scales are close to each other.The consequence of this is shown in FIG.7 where we show for the Υ(1S) state gluo-dissociation dominates in a wide temperature region of interest.FIG. 8 shows that for the Υ(2S) state also both contributions are comparable.
Finally, we find that (FIG.6) the imaginary potential over-predicts the contribution from LD substantially.
The plan of the paper is as follows.In Sec.II we will review the formalism and highlight the assumptions and approximations involved in our method.In Sec.III, we discuss the connection between EE correlator and the momentum diffusion coefficient of heavy quark, particularly, in the static limit.In Sec.IV, we discuss the implementation of the real part of the singlet potential to obtain the singlet wave function at a given T. We also discuss the two extreme cases of complete screening and no screening for octet interactions.Finally, in Sec.V, we discuss our results followed by the conclusion and future directions in Sec.VI.

II. FORMALISM
pNRQCD [1] is an EFT for bound states of quarkonia.In vacuum, it relies on the hierarchy of scales M One can think of it as a two-step process where relativistic dynamics of Q and Q are integrated out first, to obtain NRQCD at scales 1/r [46].If 1 r T, m D , the energies corresponding to the thermal scales is much smaller than the relative momentum between Q Q (∼ 1/r) then NRQCD at this scale is unaffected by T, m D , and hence this theory is the same as the theory in vacuum [46].
The pNRQCD lagrangian is obtained by integrating out modes from 1/r to E b .The structure of the EFT is governed only by the symmetries and the particle content of the theory.In the rest frame of Q Q in vacuum, this theory has been extensively studied [1].
While the medium introduces a new vector u µ associated with the medium rest frame which can lead to additional operators in the lagrangian [43,44], we only consider the case where the quarkonium is (nearly) at rest in the medium, and hence the form of the lagrangian is unchanged from that in vacuum.The lagrangian is of the form Here S(O) is the singlet (octet) wavefunction in the relative coordinate between Q and Q.
The lagrangian (Eq. 1) is obtained by systematically performing a multipole expansion which encodes the factorization of wavelengths of the order of 1/E b compared to the short distance r.The dots represent higher order terms in this expansion.
The low energy coefficients (LEC's) V A (r), V B (r) are 1 at leading order in perturbation theory and are expected to be close to 1 at a short distance.In our paper, we will take them to be 1.The other input to the theory are the potentials, V s (r) and V o (r).If E b ∼ T, m D , then the integration of modes from 1/r to E b is affected by the medium and hence the functional forms of V s (r) and V o (r) is different from their forms in vacuum.If the Q and Q can be treated as static (for example if their mass is so high that their kinetic energy can be ignored), then one can run the integration of modes all the way to zero energy.It is well known that in this limit V s and V o are complex [4,6].The real and imaginary parts of the static potentials have been calculated in weak coupling limit [4,6,47].Moreover, we also expect that non-perturbative contributions to the potential are substantial especially for the excited states of bottomonia, because while 1/r is large compared to Λ QCD the hierarchy is not very strong.Additionally, neither E b nor T are much larger than Λ QCD and hence the medium itself at this scale is strongly coupled.
Let us note that the static calculation gives the Q Q potential under the assumption that E b is the smallest scale in the problem and the kinetic energies of the Q and Q are negligible.However, if E b and T are comparable, one needs to go beyond the static approximation.In this case the thermal losses can not be captured by an imaginary potential.
In this work we assume that the real part of the potential at scale E b is captured by the static value.Thus, we take the real part of the potentials to be in the range considered by lattice inspired potentials [38].On the other hand to estimate losses due to thermal process, we compute the imaginary part of the singlet self energy diagram in the multipole expansion.The key assumption here is that in the pNRQCD lagrangian (Eq.1) at scale ∼ E b , the imaginary part of the potential is small and losses predominantly arise from dynamics at scales ∼ E b .The advantage of this approach is that it captures finite frequency effects of thermal excitation due to the medium.However, this approach misses finite frequency effects in the real parts of the self-energy.These corrections will change the energy of the singlet states and similarly for the octet state and hence change the binding energy of the singlet states.However, the effect of this contribution on the decay rate of the singlet state is higher order in the multipole expansion (r 4 instead of r 2 ) and can be safely ignored in our calculation.
Finally, we assume that the octet state, once formed, decoheres rapidly and can no longer lead to a reformation of the singlet state.In quantum calculations [36,39,53] these processes can be taken into account but this is beyond the scope of our paper.
FIG. 1: Cut diagrams contributing to the decay width.Single solid line is for singlet and double lines for octet.
The gluon line corresponds to a dressed gluon.
With this setup, let us start from a singlet state.The dissociation is given by the imaginary part of the singlet self-energy correction.The corresponding diagram is shown in FIG. 1.Here, the gluon line is resummed and gets contributions both from LD which arises from the imaginary part of the gluon self-energy and pole of the gluon propagator.In FIG. 1, gluon momentum (k 0 ,k) is directed inward at the first vertex and octet momentum (q 0 ,q) is directed outward from the same vertex.The incoming momentum of the singlet is p µ = q µ − k µ .
In order to calculate the imaginary part of the singlet self-energy, we follow the cutting rules at the finite temperature are given in Refs.[54,55].There are two cut diagrams as shown in FIG. 1 .Following the cut rules and implementing appropriate propagator for each cut diagram, the imaginary part of the self-energy reads as To avoid this lengthy expression, from here onwards, we use the following shorthand notation for the above equation In Eq. 2, C F (= 4/3) is color factor.ρ 00 (k 0 , k) and ρ jj (k 0 , k) are gluon spectral functions that are discussed in the next section.p µ = (p 0 , 0) is the four-momentum of the incoming singlet state and ρ o is the tree level Q Q spectral function in the octet channel which can be obtained from the octet propagator (see Eq. 1) Hence, To proceed further we need information about the temporal and spatial gluonic spectral functions, ρ 00 and ρ ii , which we discuss below.

A. Gluon-polarization tensor
In this section, we review the well known expressions for the gluon polarization tensor [56] that are essential inputs to the evaluation of the imaginary part of quarkonium self-energy.
The general form of the gluon self-energy is given as where P L µν (P T µν ) are longitudinal (transverse) projection operators and Π L (Π T ) are component of the gluon selfenergy along these directions.In order to evaluate quarkonium decay width using Eq. 3, we need the imaginary part of the gluon propagator.This may come from the pole of the propagator in the region of phase space where the imaginary part of the gluon self-energy is zero, and from the region where the imaginary part of the gluon self-energy is finite.The pole contribution is nonvanishing in the limit k 0 > k, and the latter contribution that requires the real and imaginary parts of the gluon self-energy is finite when k 0 < k [56].Below we discuss various components of the gluon self-energy.
Let us first consider the regime k 0 < k (space-like).This we call the Landau damping regime.
The gluon loop contribution to the imaginary part of the longitudinal component of the gluon self-energy is given as where N = 3 and f (q) is Bose-Einstein distribution function.Let us note that with an expansion in k 0 /T in the distribution function and by taking the limit k q, Eq.7 goes to its hard thermal loop (HTL) counterpart.
Similarly, The quark loop contribution with N f (light) quark flavors to the imaginary part of the longitudinal component of the gluon self-energy is given as where f (q) is Fermi-Dirac distribution function.The total imaginary part of the longitudinal gluon self-energy can be obtained by summing Eqs.7 and 8.
The real part of the longitudinal component of the selfenergy is where In obtaining Eq. 9, we have dropped terms of the order of k 0 /T, k/T .These terms are important when k 0 , k T but we drop these terms because of the following reason.The exact forms of these higher order terms depend on the gauge (for eg.see [56]) while Eq. 9 is the HTL form and is gauge invariant [8].This expression is valid for both k 0 > k and k 0 < k.Moreover, from Eq. 3 it is clear that the contribution to Σ from k 0 , k T is exponentially suppressed and hence making this approximation will not cause a significant error in our result.
Similarly, for the transverse gluon, the imaginary contributions to the gluon self-energy are, and, The real part of the transverse component of the gluon self-energy is Let us now consider the regime k 0 > k (time-like).This we call the pole regime.In this regime, the imaginary part of Π L and Π T are zero and the real parts are as above.At order g 3 the widths of these modes are finite [56,57] but we ignore this in our calculation.Now we can calculate the gluon spectral function which goes in Eq. 3. Below we discuss it for both time-like and space-like gluons.

B. Gluon spectral functions
The general form of the gluon spectral function in a medium reads as Here Below we discuss the form of these spectral functions for both k 0 > k as well as k 0 < k.
For k 0 < k (LD), we use the gluon self energies (shown in the previous section) to write the resummed gluon propagator and obtain where Π L ( Π L ) is sum of both gluon and quark contributions.
Similarly, the transverse component of the spectral function can be written as (15) It is worth mentioning here that the above form of the spectral functions reproduces the momentum diffusion coefficients obtained within the kinetic theory framework in Ref. [40].
For k 0 > k the gluon propagator is simply a pole.The quarkonium dissociation in this regime is due to the absorption of a gluon from thermal medium.This process is known as gluo-dissociation in the literature.In the limit k 0 T , the spectral function is given by the imaginary part of free gluon retarded propagator and gluodissociation in this case has been studied in Refs.[39,43] However, for realistic situations one needs to take full resummed propagator.Thus, similar to Eq. 13 the general form of the spectral function reads as where p stands for pole.The longitudinal spectral function in this regime is given by The transverse spectral function is given by The imaginary part of Σ 11 gets contribution from both Eqs. 13 and 16.While in the low frequency limit, LD gives the dominant contribution, pole contributions are significantly large in the intermediate and high frequency limit.The overall pole contribution merges with their free spectral function counterpart at an asymptotically large frequency.This we show in FIG. 3.

III. CONNECTION WITH THE MOMENTUM DIFFUSION COEFFICIENT
In this section, we relate Eq. 3 with the standard definition of the momentum diffusion coefficient in terms of the electric field correlator which is given as [58,59] where U is the Wilson line in the fundamental representation, E i = ∂ i A 0 − ∂ 0 A i is the color electric field and trace over color degrees of freedom.In Eq. 19, the infinite integration limit represents the zero frequency limit of the correlator.For the leading order results, one needs to replace the Wilson lines by identity (i.e., U = 1) to obtain for more details see Refs.[40,59,60].
It is useful to compare this quantity (Eq.20) to the expression Eq. 3. Imposing the condition q 0 > 0 in Eq. 3 and performing energy integration in Eq. 2 using the energy delta function we rewrite the imaginary part of the singlet self-energy as where k 0 = p 0 − q 0 .For future use, we define a quantity A single heavy quark, traversing through the thermal medium, gets uncorrelated random kicks from the medium constituents that give rise to κ.However, for quarkonium bound state, not only scattering but also absorption of thermal gluons contribute to the dissociation.The latter process is kinematically forbidden for a single heavy quark.Therefore, in the frequency regime where gluo-dissociation dominates κ is not the same as κ.However, in the static limit where dissociation via scattering (i.e., LD) is dominant, the two coefficients defined in Eqs.19 and 22 seem identical, at the leading order.We have checked that the LD part of κ agrees with the one obtained in Ref. [40].
We note that at higher order there is no reason for these two coefficients to be identical.The reason is that for quarkonium, chromo-electric field correlator is defined with Wilson lines in the adjoint representation [61].On the other hand for a single heavy quark, Wilson lines are in fundamental representation.
It is useful to note here that in Eq. 22 if one makes k 0 , k T approximation before integrating over k, the integrand is of the HTL form and is ultraviolet (UV) divergent.This is due to the fact that the applicability of the HTL resummation is restricted to the low frequency limit.This divergence can be cured by either using a cutoff k ∼ m D or by adding the UV contribution to the integral carefully [44].However, the contribution from Eqs.14 and 15 vanishes in the high frequency limit.Therefore, if we use the imaginary part of the longitudinal self-energy given by Eqs. 7 and 8 and evaluate κ, the integral is convergent and can be computed numerically which we do next.In FIG. 2, we plot the electric field correlator arising from various contributions that appear in the evaluation of singlet self-energy diagram in FIG. 1 as a function of frequency.Here we take T = 0.25 GeV and Debye mass m D = 0.5 GeV.The black curves are for the longitudinal gluon with a solid line for LD and a dashed one for pole contributions.The red curves are for transverse gluon where solid and dashed lines are for LD and pole contributions, respectively.As anticipated, in the small frequency limit the dominant contribution comes from longitudinal gluon Landau damping.In this limit, other contributions are either zero or very small.Pole contributions switch on at a somewhat larger frequency, i.e., k 0 ∼ 0.2 GeV.Moreover, the transverse gluon pole contribution is larger (in magnitude) compared to the longitudinal one.Finally, at high frequency, transverse pole contribution dominates and eventually approaches the corresponding free limit.
In FIG. 3, we have plotted κ/T 3 as a function of k 0 .The red (dashed) curve here gets contribution from both k 0 < k as well as k 0 > k phase space regions.In the static limit, i.e., k 0 ≈ 0, we have checked that κ/T 3 agrees with that in Ref. [40].The black (solid) line is in the free limit which as expected is zero at zero frequency.
It is well known that the perturbative result for κ is too low by roughly a factor of 5 − 10 than the nonperturbative value [62][63][64][65].For example, recent lattice results for κ/T 3 are estimated as [66] 1.99 For finite k 0 there are no lattice results available in the literature.Naively, we expect them to be different from the perturbative estimates.Therefore, we expect that our predictions for R AA are underestimated.However, our results capture the qualitative features of relative contributions of pole and LD in the range of temperatures available in HIC.Motivated by lattice QCD calculations of m D , we choose g = 2 [38,67].
In the intermediate frequency regime, the peak structure in κ/T 3 is from the transverse pole contribution.Finally, in the high frequency limit, it merges with its free limit counterpart.

IV. DECAY WIDTH
At any given temperature T , the decay width of a singlet state |φ at leading order is given by where Σ 11 is given in Eq. 3. Inserting a complete set of octet states |o o| at the right bracket of Eq. 3, one obtains where summation is over all octet states allowed by the selection rule.In operator form O(p 0 , p) is same as defined in Eq. 3 and now p 0 = q 0 − k 0 .Let us note that the octet state lies in the continuum.q 0 is the energy of the octet state, which is given by, One point to note is that momentum conservation implies that the center of mass momentum of the octet state is −k.This implies that the energy of the |o state has an additional contribution k 2 /(4M ) which should be added to the right hand side of Eq. 25.The value of k is governed by T since it is the region in k space where the gluon spectral function multiplied by the Bose-Einstein distribution function is not exponentially suppressed.In the hierarchy we are working, E b and T are both small scales compared to M and hence quantities of the order of T 2 /(4M ) are suppressed by an extra power of M compared to the right hand side of Eq. 25 and hence can be safely dropped.

A. Modelling the singlet state
To complete the evaluation of Eq. 24, we need the functional form of the singlet state |φ as well as octet state |o .Below we discuss the prescriptions to obtain these wave functions.
For a state created in vacuum and "dropped" into the QGP, a natural choice for |φ is the wavefunction in vacuum.Further, if the thermal effects are weak then |o can be taken to as octet states in vacuum.If the initial formation of quarkonia is not affected by the medium (for example if the formation of the quarkonium states occurs on a time scale much shorter than the formation of the QGP) then this is a well motivated model for |φ and |o .This picture has been previously used for phenomenology [68,69].
At the LHC and the RHIC, the formation time of the QGP is a fraction of a fm/c and is not substantially larger than the formation time of quarkonia, of the order of 46   I: Binding energies, mass and r 2 of various bounds state at T = 0 using Eq.27.All dimensions are in GeV.
1/E b .One can expect the formation dynamics of quarkonia to be substantially affected by the medium.One natural way to include these effects is to start the evolution from a narrow initial Q Q state of width ∼ M and follow its quantum evolution from very early time [36,38].In this paper, we do not study the quantum dynamics and this analysis is beyond the scope of the paper.If dissociation can be modelled by the imaginary part of the potential (i.e. in the E b T regime) another possible approach is to assume that the evolution dynamics is slow (adiabatic approximation) and the quarkonium state is initially formed in the eigenstate of the complex potential and at each instant the quarkonium state is in the eigenstate of the complex potential [24,70].In this paper, dissociation is calculated using Eq. 3 which can not be captured by a complex potential and hence the adiabatic method is not applicable.We model the effect of the medium on the formation of quarkonia by making the maximal approximation that the initial state and subsequent to formation is determined by the real part of the instantaneous thermal potential [17].
More concretely, for the singlet states wavefunction, we use the eigenstates where V s (r, T ) is the real part of the thermal potential.
Here we have subtracted the rest energy from all the Q Q states.Similarly, |o is given by Eq. 25 with V o given by the real part of the octet potential in the thermal medium.In summary, to calculate |φ and |o we need the real parts of the potentials V s , V o .
For the singlet potential we use the lattice inspired potential which is given by Refs.[38,71] The effective coupling a = 0.409 and the string tension σ = 0.21 GeV 2 are fixed from the vacuum masses and binding energies (see TABLE I) with bottom mass M = 4.7 GeV.Here we take m D = 0 for obtaining the vacuum spectrum.
For finite T we keep a and σ the same as in T = 0 and m D = (1 + N f /6)gT .For g = 2 Eq.27 gives potentials consistent with those used for bottomonium phenomenology with lattice based potentials [27,72].The Q Q potential in the medium is screened, as a result of which E b becomes smaller with increasing temperature.At sufficiently high temperature the bound state is dissolved [3].It is worth mentioning that for the 1S state, the wavefunction does not depend on the temperature of the medium up to T ∼ 480 MeV and it remains approximately the same as that of vacuum Coulombic state while the excited states dissolve earlier [41].
In FIG. 4, we plot the binding energy of Υ(1S), Υ(2S) and Υ(3S) states as a function of medium temperature.For a given potential, E b is given by E b = 2M −M M +V ∞ , where M is bottom current mass, M M is bound state mass and V ∞ is the asymptotic value of real part of singlet potential.The key point we want to highlight in FIG. 4 is that for the temperature range relevant for HICs, the hierarchies E b T or E b T may not be satisfied, at least for 1S and 2S states.It is worth mentioning that E b obtained here agrees with the one in Ref. [73].On the other hand for higher states, binding energy approaches zero around this temperature and E b T .A consequence of our model is that we can not address the observed phenomenology of the 3S state [74] as for the central bins R AA for 3S state (in our model) is zero.
The key dynamics missing from the classical model are (1) quantum formation dynamics and (2) processes that allow for reformation of bound states which are important for capturing 3S dynamics.FIG.4: Binding energies for Υ(1S), Υ(2S) and Υ(3S) states as a function of temperature.Here we take constant coupling g = 2.

B. Modelling the octet state
The octet states are also affected by the thermal medium.At short distances we know that the potential is repulsive Coulombic.In perturbation theory both the real and the imaginary parts of the medium modified octet potential have been computed [47].Recently, both the real and imaginary parts of the octet potential have also been computed in pure gluonic theory on the lattice [45].An important outcome from these pa-pers is that in both perturbative and non-perturbative calculations one finds that at large r the singlet potential approaches the octet potential.This is still an active area of research but the form of the potential for 2 + 1 flavor QCD is not yet known.
Based on these considerations, we take two limiting cases for the octet potential.The first is when the screening is strong and one can ignore the octet repulsion at the distance scale of interest and hence where V ∞ = 2σ m D is the asymptotic value of singlet potential.In this case, the |o wavefunctions are the same as the wavefunctions of the free particle but the energy levels start from V ∞ .
The other limit is where the screening is weak and We expect the true physics to be between these two limiting cases.
Finally, the octet potential also has an imaginary piece, which corresponds to the change of the octet state to the singlet state.In this paper, we assume that these processes do not regenerate bound states because the octet states are much broader than the singlet state.
For the repulsive Coulombic potential, the general form of the radial wave function is given as [75] where 1 F 1 is confluent hypergeometric function, ρ = rp, ν = 1 8a0p with a 0 = 2 αM as Bohr radius.In this work, we take a 0 = 1/1.334[36].The normalization factor C l in Eq.30 reads as Let us note that for the above form of C l , the wavefunction obeys the following form for normalisation With this choice of normalisation, the decay width has the form given in Eqs.36,37 and 39.With the above form of the radial wave function, the general form of the octet wave function can be written as Here Y l m is spherical harmonics.For the 1P state, replacing l = 1 in the obove equation and summing over quantum number m, the octet wavefunction |o reads as For s(d) states, the wave function can be obtained by replacing l = 0(2) in Eq. 33.
In the case of no final state interaction, we take free wave function which in terms of Bessel function is given as Below we discuss the contribution from the Landau damping and gluon absorption in the decay width.

C. Landau damping contribution
One way to organize the kinematic regimes which contribute to the decay width (Eq.23) is space-like and time-like.Dissociation (of quarkonia) via scattering of the bound state with the thermal partons occurs when the momentum of the exchanged gluon is space-like.As mentioned, this mechanism is known as LD.It is estimated by taking the resummed propagator for the gluon line in FIG. 1.The LD contribution is dominant when . The overall contribution is both from the longitudinal as well as transverse gluon.Moreover, in the small k 0 limit, the longitudinal gluon contribution turns out to be dominant.This is easily understood by noting that while both ρ L and ρ T go as k 0 at small k 0 , ρ T is multiplied by (k 0 ) 2 in Eq. 3 while ρ L is multiplied by (k i ) 2 .For E b m D , T , k 0 k, and hence the longitudinal gluon contribution turns out to be dominant.However, we evaluate both of them numerically and do not drop the transverse contribution.
The explicit expression of Γ L for gluonic system is given (see Eq. 24) by, For finite N f , Π 00 and Π 00 gets contribution from the gluon as well as quark loop.Here p is octet momentum and k 0 = p 2 M − E where E is binding energy of the bound state.
Assuming a strong hierarchy between E b and T , this can be simplified to the one obtained in Ref. [6].
The transverse gluon contribution arises from the transverse part of the spectral function given in Eq.15.
As may be noticed in Eq. 3, for finite k 0 this term can not be ignored.Moreover, at significantly large k 0 , this term can dominate over the longitudinal gluon contribution.This situation may arise when the hierarchy between E b and medium temperature is not very strong.This can be observed from FIG. 2 that for k 0 sufficiently large, the transverse contribution to Σ 11 can be larger than the longitudinal.It is also clear that in the kinematic regime where the transverse LD contribution is comparable to the longitudinal LD contribution, the binding energy and the final state interactions can not be ignored in either.
Following the similar prescription as for the case of longitudinal gluon, the contribution to the decay width from the transverse gluon reads as Total contribution from Landau damping is sum of Eqs.36 and 37.

D. Pole contribution
In the kinematic regime E b T , the dominant contribution to quarkonium dissociation from singlet to unbound octet state occurs by absorbing a time-like gluon from the thermal medium.This process is known as gluodissociation [43].For the decay width evaluation, the contribution to the quarkonium self-energy arises from the pole of the gluon propagator.In the free limit, i.e., T E b , the medium contribution to the singlet to octet thermal breakup appears in the thermal weight only.However, in the intermediate temperature range, HTL effects also become important and one needs to use resummed propagator.In the HTL resummed propagator, both longitudinal and transverse gluon contribute to the imaginary part of the gluon propagator.After adding both of these contributions and performing k 0 integration using energy delta function, Σ 11 reads as where In the limit T E b , Eq. 38 can be solved analytically by making an expansion in the bosonic distribution and replacing the spectral function by the free spectral function.Using Eq. 38, decay width is given as The contribution of transverse gluon to the decay width of the bound states is dominant over the longitudinal one.This can be observed from the frequency behavior of Σ 11 shown in FIG. 2. Below we discuss the results obtained for the decay width and R AA .

V. RESULTS FOR CLASSICAL DYNAMICS
In this section, we discuss the results for Υ(1S) and Υ(2S) states using the decay width obtained in the previous section.We mainly focus on the relative contributions of pole and LD within the perturbative limit and show the overall effect on R AA .
At any given time t, the survival probability of a given singlet state can be obtained by using the rate equation where N (t) is the number of bound states at time t and summation is over the two contributions arising from LD and pole.The total number of states produced after some time t f may be obtained from Eq. 40 and is given as Here t 0 is the initial time, N 0 is number of bound states at time t 0 , and Γ(t) = i Γ i (t).

A. Medium model
To calculate Γ as a function of time we need a model for the background evolution of the thermal medium.In this paper we use a simple model and take the medium by a Bjorken expanding medium [76] with a temperature We are interested in quarkonia with zero rapidity and hence have replaced the proper time with the local time.T (t 0 ) is the temperature at a reference proper time t 0 .We take t 0 = 0.6 fm which is a little later than the start of hydrodynamics for LHC energies [77].This is starting time for the quarkonium evolution and is comparable to the formation time for quarkonia.Similar numbers were taken in Ref. [37].The temperature at time t 0 depends on the impact parameter or equivalently the centrality.We use centrality bins analogous to the one used in ALICE [78].For this purpose, we use the Glauber model (we do not use the Monte-Carlo Glauber model here though we see that the difference between our centrality bins and the bins obtained from the Monte-Carlo Glauber model [79] is small) to relate the impact parameter to the number of participants (N part ) and the number of binary collisions(N bin ).Following ALICE, bins in the observed dN ch /dη are related to bins in N part using the relation where λ = 2.75 and f = 0.212 [78].With these values of the parameters λ and f , we quantitatively agree with  With dN ch /dη in hand for each centrality bin, the initial temperature at t 0 is obtained by using following prescription [81].The value of dN /dy is related to the experimentally measured charged particle multiplicity by the relation, where J = 1.12 is the jacobian for y to η transformation [82].N is the multiplicity of the particles produced at the end, y is pseudo-rapidity.The factor 3/2 in Eq. 44 incorporates the contribution from neutral particles.In order to be consistent within our model, for each centrality bin, we take dN ch /dη obtained by using Eq.43 with the parameters mentioned above.
In the nuclear collision experiments, the number density of partons or entropy is decided by the rapidity distribution of the produced particles.Assuming that initially the system is at chemical equilibrium, and assuming that the entropy does not change during the evolution, we can write the initial parton density as where A ⊥ is transverse size of the system.We estimate it from the Glauber model for each centrality via [83] A where x 2 , y 2 is size along the transverse directions.
A rough estimate of T 0 can be made from the initial number density assuming a non-interacting QGP.Then, n 0 = (β 1 +2β 2 )T 3 0 with β 1 = 8π 2 /15 and β 2 = 7π 2 N f /40 represents equilibrium density for gluon and quark.Plugging it back in Eq.45 along with Eq.44, we obtain the initial temperature at time t 0 to be For √ s = 2.76 TeV, Eq.47 gives T 0 ∼ 250 MeV for the most central bin (0 − 2.5%).This is significantly lower than estimates of the temperatures obtained at t 0 ≈ 0.6 fm in hydrodynamic simulations [77].We note that Eq.45 ignores various effects such as particalization, viscous effects, interaction in the QGP, and expansion in the transverse direction.Moreover, losses in the work done by the system during the expansion are also not taken into account.These effects may lead to entropy generation and energy loss restricting the applicability of Bjorken expansion.Thus the initial parton density will be somewhat larger than estimated by Eq. 45.In order to take these losses into account we redefine the initial parton density by multiplying Eq. 45 with a fudge factor (f ) i.e., n 0 = f n 0 .This fudge factor is adjusted in such a way that we obtain the initial temperature in the most central bin to range from T 0 ≈ 450 MeV (f = 1.46) to T 0 ≈ 480 MeV (f = 1.78) at t 0 = 0.6 fm.For our final results of R AA we provide the band corresponding to these two values of the initial temperatures.While this is a highly simplified model for the background medium, we hope that this band of variation captures important features of its hydrodynamic evolution.
In FIG. 5, we show the variation of temperature as a function of time starting with t 0 = 0.6 fm till the temperature reaches the final value of 190 MeV.
Finally, for completeness, we also consider the final state feed-down effect, we follow the prescription given in Ref. [38].Therefore, for a given bound states, R AA is given by where N h is number of higher states at the end of evolution and α is feed-down parameter from higher state (N h ) to the state being considered.Feed-down matrix is given in Eq.(5.2) of Ref. [38].N 0 is number of states at the initial time t 0 .In the results presented here, for feed-down , we only consider higher states with contributions more than 10%.Thus for Υ(1S), we take the contribution from Υ(2S), χ b0 (1P),χ b1 (1P) and χ b1 (2P) states.We follow the same for Υ(2S) state as well.For more details on feed down see Ref. [38].Below we discuss the results in somewhat more detail.

B. Results
In FIG. 6, we compare the decay width (of Υ(1S) and Υ(2S) states) from Eq. 36 with the (imaginary) potential [84] formalism which is extensively used in the literature.In the latter case the decay width is obtained from where V s is the imaginary part of the singlet potential.In order to be consistent with the r 2 expansion in Eq. 1, we also make this approximation in the imaginary part of the potential.In this limit, the momentum integration in the potential becomes divergent, limiting its applicability in the small momentum ranges.We therefore put an upper cut off m D to get finite results.The corresponding decay width for Υ(1S) and Υ(2S) states are shown by the solid line in FIG. 6.In the same figure, we show the decay width obtained from Eq. 36 with E b = 0.6 GeV (for 1S) and 0.2 GeV for 2S.As may be noticed, V s (r, T ) gives a larger decay width which can be understood as follows.
Eq. 49 assumes the binding energy of the singlet state, and the octet state energy (q 2 /M ) of the species is negligible compared to the temperature.This can be seen by inserting a complete set of octet states in Eq. 36.The energy delta function gives k 0 = q 2 /M − E b .Taking the limit k 0 → 0, Eq.36 gives the result obtained from the imaginary potential expanded in r to order 2. We emphasize here that for realistic values, this approximation over-predicts the decay rate.We have checked that if we do not make an expansion in r in the imaginary potential, the decay width turns out to be even larger.
In FIG. 7, we plot the scattering (LD) and gluo dissociation (pole) contributions to Υ(1S) state.For Υ(2S) we plot the same quantity in FIG. 8.The bands correspond to no screening (Eq.33 In FIG. 9, we combine both gluon absorption and scattering processes and plot R AA for Υ(1S) state as a function of N part for each centrality bin.
For a given centrality, we first evaluate singlet state wave functions at each temperature by solving the Schrodinger equation given in Eq.26.With temperature dependent wave functions at hand, we estimate matrix element i.e., | φ|r|o | 2 (see Eq.36) and obtain corresponding decay width for all processes.Here we use temperature dependent binding energy shown in FIG. 4. Finally, we use Eq.48 to estimate R AA for a given state.As mentioned earlier, we consider two cases ; (a) when there is no screening in the final state octet interaction (show in blue color)and (b) when octet states are completely screened (shown in red color).Here, upper and lower edge of the band corresponds to the lower and higher initial temperature given in TABLE II.
Let us note that in the realistic situations, octet potential is neither completely screened nor has form of pure vacuum potential.Therefore, we expect the true value of suppression to lie somewhere between the blue and the red bands.It is evident that the suppression is underpredicted for these states.This may largely be due to the fact that the spectral functions are obtained at the leading order (LO) accuracy.Further improvements of our results require non-perturbative estimate of these spectral functions on lattice which so far have been computed only in the static limit.

VI. CONCLUSION
In this work, we perform a comprehensive analysis to quantify the relative contributions arising from LD and gluo-dissociation to the Υ(1S) and Υ(2S) states within leading order perturbation theory for the medium.We assume that the real part of the singlet potential is screened in a thermal medium, and we take lattice motivated real parts of the Q Q potential and estimate the binding energy and wave function as a function of T for Υ and χ states.By comparing the binding energy (shown in FIG. 4) and temperature it can be seen that for the 1S and 2S states both E b and T are comparable to each other (although for the 2S state for T > 300MeV one could take T E b ).We, therefore, expect that neither E b T nor E b T hierarchy is strictly satisfied in the interesting range of temperatures for the QGP, at least for 1S and 2S.However, for other higher excited states, E b T seems to be satisfied quite well.For Υ(3S) this can be seen in FIG. 4.
The consequence for the 1S state (and for the 2S in the later stages of the evolution) is that LD can not adequately describe the decay of the state.It is well known that in the static limit i.e., k 0 → 0, decay dominantly happens via longitudinal LD which can be captured by an imaginary potential between Q Q.However, this limit is valid if the binding energy of the quarkonium species is negligible compared to the medium temperature.As discussed in the text below FIG.6, this also requires the kinetic energy of the octet states to be small.For E b ∼ T , the full finite frequency region of the spectral function for the chromo-electric field correlator becomes important.
There is a caveat to the above conclusion.The perturbative value of the k 0 → 0 value of the spectral function is a factor of 5 − 10 times smaller than the non-perturbative value calculated on the lattice.If this enhancement persists till k 0 values of a few 100 MeV (see FIG. 3) then our conclusion about the relative contribution of the LD and the other contributions will be still be true.However, if the spectral function rapidly drops down towards the perturbative estimate (at k 0 ∼ GeV we expect the leading order perturbative result to be more reliable) then the LD contribution dominates for all the states and the static results can be used to capture the physics of interest.This provides a motivation to study the spectral function at finite frequency non-perturbatively, although it is a difficult problem.
The hierarchy between E b and T (T is related to the inverse of the environment time scale) plays an important role while performing a quantum calculations of quarkonia dynamics within the open quantum system framework.It has been shown [36,85,86] that if E b T quarkonium system is in the quantum Brownian motion regime and the evolution is local in time.In this case, one can obtain a Lindblad equation for the density matrix evolution for the Q Q system.On the other hand if E b T then this proof fails.We show here that this is the case for the Υ(1S) state throughout the evolution if one considers the state to be an eigenstate of the instantaneous (real) potential.E b is also T for the excited states in the early part of the evolution.However, If substantial regeneration of the excited states happens at the later time (low temperature) then the hierarchy E b T might be satisfied.In this regime if a quantum Brownian description of the dynamics is attempted, one might need dynamics that are correlated in time.Although a quantum optical description might become applicable in this case [87,88].We do not discuss quantum evolution in this work and leave this for the future.

1 r 1 r
E b .The scale separation M ensures that the Q and Q are non-relativistic, and 1 r E b means that the interactions between Q and Q (at leading order in 1/M ) can be written as potentials.
of the spectral function and D R(A) L is longitudinal component of resummed retarded (advanced) gluon propagator.Similarly, one can obtain the transverse component of the spectral function (ρ T ) by using ρ

3 FIG. 2 :
FIG. 2: Various contributions to scaled EE correlator as a function of frequency.Here we take T = 0.25 GeV and m D = 0.5 GeV.

3 FIG. 3
FIG.3: κ/T 3 as a function of frequency for constant coupling g = 2, T = 0.25 GeV and N f = 3. Black (solid) line is the free limit and the dashed (red) is resummed one.

FIG. 5 :
FIG. 5: Temperature as a function of time.Red band is for the most central bin and, central and green one are for 20-30% and 40-50% centrality, respectively.Here the initial time t 0 = 0.6 fm for all centrality bins.The band here covers upper and lower temperature correspond to two different values of the fudge factor.

FIG. 7 :
FIG. 7: Decay width of the Υ(1S) state as a function of medium temperature.The band corresponds to the free wave function (upper curve) (Eq.35) and Coloumbic repulsive wavefunction (lower curve) (Eq.33).Red color band is LD and blue color band is gluo-dissociation.

TABLE II :
[80] value of N part , impact parameter (b) and initial temperature (T 0 ) for various centrality bins.T 0 is temperature at time t 0 = 0.6 fm.The left column in T 0 is for f = 1.462 and the right one is for f = 1.782(see text below Eq. 47).the centrality dependence of dN ch /dη as a function of N part given in FIG.10of Ref.[80].For each centrality bin, we have mentioned the mean value of N part and impact parameter b in TABLEII.