Collisional and thermal dissociation of $J/\psi$ and $\Upsilon$ states at the LHC

We present new results for the suppression of high transverse momentum charmonium [$J/\psi, \psi(2S)$] and bottomonium [$\Upsilon(1S),\Upsilon(2S),\Upsilon(3S)$] states in Pb+Pb collisions at the Large Hadron Collider. Our theoretical formalism combines the collisional dissociation of quarkonia, as they propagate in the quark-gluon plasma, with the thermal wavefunction effects due to the screening of the $Q\bar{Q}$ attractive potential in the medium. We find that a good description of the relative suppression of the ground and higher excited quarkonium states, transverse momentum and centrality distributions is achieved, when comparison to measurements at a center-of-mass energy of 2.76 TeV is performed. Theoretical predictions for the highest Pb+Pb center-of-mass energy of 5.02 TeV at the LHC, where new experimental results are being finalized, are also presented.


I. INTRODUCTION
The fate of quarkonia -for example the J/ψ and the Υ meson families -in a thermal medium, such as the quark-gluon plasma (QGP) created in heavy ion collisions (HIC), can help us characterize its properties. In particular, quarkonia are sensitive to the space-time temperature profile and transport coefficients of the QGP, see [1][2][3] for recent reviews. Experimentally, a key observable that carries such information is the nuclear modification factor of the yields of quarkonia in nucleusnucleus (AA) collisions, when compared to their yields in nucleon-nucleon (N N ) collisions scaled with the number of binary interactions In HICs one expects that the short distance formation dynamics of a QQ pair is not affected since m Q ≫ T where T is the typical temperature of the QGP. To simplify the calculations, it is often also assumed that the matrix elements for transition of QQ to quarkonia is not modified, and for every binary collision the formation of a specific quarkonium state happens with the same probabilities as in N N collisions. This takes a time scale close to its inverse of its binding energy. However, due to the screening of the color interaction between Q and Q in a deconfined QGP [4], as well as processes leading to the dissociation [5] of quarkonium states, we expect the yields of quarkonia to be suppressed in heavy ion collisions (R AA < 1).
Several methods have been used to estimate screening and dissociation effects encountered by a quarkonium in a thermal medium. A widely used approach is based on the intuitive idea that the real part of the "finite temperature potential" between two (nearly) static heavy quarks captures the screened QQ interaction while the imaginary part of the potential captures dissociation. For T = 0, the real part can be quantitatively obtained by calculating the Polyakov loop correlation functions [6][7][8]. For T > 0, the connection between various correlators calculated on the lattice and the potential between Q and Q is subtle [6][7][8][9]. While the singlet free energy F 1 (r) of the QQ state as a function of the separation r, and internal energy U (r) can be extracted from the lattice data, the connection of either of the two with the QQ potential is indirect. An important step in clarifying this connection was taken in Ref. [10], which showed that the "potential" between two heavy quarks is complex, with the imaginary part connected to thermal processes that can lead to dissociation of quarkonia. Significant progress has been made in the perturbative calculation of the real and imaginary part of the QQ potential for QQ at rest in the medium [10][11][12] or moving slowly in the medium [13]. Analytic calculations can be performed in certain regimes by considering various hierarchies of energy scales. For example, one has to assume that T and m Q are large enough so that perturbation theory is valid all the way down to energy scales πT and the binding energy E b .
While the real part of the potential for T ≫ Λ QCD can be obtained using perturbation theory, non-perturbative effects are substantial near the crossover temperature and it is better to estimate this quantity from lattice calculations by matching the Euclidean correlators measured on the lattice [8,9,14,15] with those evaluated using the potential. The extraction of the imaginary part of the potential using this technique is challenging and often perturbative estimates are used. Extensive phenomenological study of quarkonium suppression by using this approach has been performed [16][17][18][19]. Recently, approaches treating the QQ as an open quantum system were developed, where a stochastic equation is written for the evolution of the QQ wavefunction [20][21][22][23]. We also note that the connection between the heavy quark correlators measured on the lattice and in a non-perturbative framework using the in medium T-matrix has been made in Refs. [24,25]. Quarkonium suppression [26][27][28], as well as low p T observables like heavy quark flow have been studied using this approach. Thermal properties of quarkonia have also been investigated in the strong coupling regime using AdS/CFT techniques, see Ref. [29].
In this paper we calculate the differential R AA as a function of p T by solving rate equations [30][31][32] describing the change in the yields as a function of time in HICs. Conceptually, our approach resembles treating the QQ as an open quantum system. It was, however, introduced earlier to describe the attenuation of open heavy flavor [30]. Gluon exchanges with the medium lead to the modification of the QQ state and, hence, a reduction in the overlap with the quarkonium wavefunction. However we do not connect the imaginary part of the potential calculated in approaches cited above to the decay rate. This is because we focus on R AA at high p T . For high p T partons traversing the medium, a very successful picture is that the interactions with the medium lead to transverse momentum broadening. The decay rate in our calculation is related to the accumulation of relative momenta between Q andQ.
The form of the rate equations is the same as used in our previous work [32] and they involve the dissociation time and the formation time as inputs. The formation time in our formalism is a measure of the time scale on which the proto-quarkonium QQ state develops interactions with the medium and its overlap with the quarkonia becomes substantial. We vary it in the neighborhood of O(1 fm). The dissociation time scale is computed as in [32] by calculating the survival probability of the quarkonia.
One conceptual change in our framework, when compared to [32], is that we use the real part of the lattice motivated thermal potentials [7,9,14] to solve for the quarkonium wavefunction and square the overlap with the thermal wavefunction to get the survival probability. This is justified if the time scales on which the medium screens the QQ interaction is smaller than the formation time as well as the dissociation time. A rough estimate of the Debye screening time is the inverse of the Debye screening mass (µ D ), which is ∼ gT . For g = 1.85, (1/µ D ) is numerically smaller than 1/2 fm for most of the evolution of the medium at the LHC. On the other hand, the formation is taken to be ∼ 1 fm. Finally, this puts a constraint that the minimum dissociation time consistent within our formalism is roughly 1 fm. We see that while there is no substantial hierarchy between 1/µ D and the formation and dissociation time scales in our calculation, numerically 1/µ D is smaller and, hence, we work in this approximation.
In addition, we make two technical improvements. First, we consider a 2 + 1 dimensional viscous hydrodynamic medium, rather than a simplified 1 + 1 dimensional Bjorken expansion, as a model for the QGP [33]. An important element of our approach is the use of nonrelativistic QCD (NRQCD) [34] to obtain the baseline nucleon-nucleon cross sections for quarkonia and understand the p T -dependent feed down. The second technical improvement is a refit of the long distance matrix elements for χ c and ψ(2S) to obtain a better description of the high p T part of their differential spectra. We are further motivated by the fact that new experimental results that extend the measurements of quarkonia to high p T have appeared. Improved measurements of the higher excited states, such as the ψ(2S) and the Υ(2S) and Υ(3S), that are differential in transverse momentum are particularly useful in constraining the mechanisms of quarkonium suppression.
In this paper we focus on the nuclear modification of the prompt J/ψ and Υ states around mid-rapidity at the LHC and at high transverse momentum [35][36][37][38][39][40]. In section II we describe the theoretical model that is employed in the calculation of quarkonium production in heavy ion collisions. Phenomenological results for the nuclear modification of the J/ψ and Υ states are presented in section III. Our conclusions and outlook are given in section IV.

II. THEORETICAL FORMALISM
Quarkonia are bound states of heavy quarks (Q) and antiquarks (Q), and a concrete picture of the dynamics of the heavy quark pair (QQ) in N N collisions is given by non-relativistic quantum chromodynamics (NRQCD) [34]. In this effective theory, the initial hard collision produces a short distance (∼ 1/m Q ) QQ pair in a color-singlet or an octet state with a specific spin and orbital structure. The production cross-section for this short-distance state can be calculated using perturbative QCD. This QQ state evolves into a quarkonium state with probabilities that are given by long distance matrix elements (LDMEs). For color-octet states, this evolution process also involves the emission of soft partons to form a net color-singlet object which we assume occurs on a time scale which is shorter than O(1 fm). This framework has been successfully used to calculate the unpolarized differential yields of quarkonia versus the transverse momentum (p T ) in proton-proton (pp) or proton-antiproton (pp)] collisions [41][42][43][44]. An accurate description of both the cross sections and polarization of quarkonia in hadronic reactions still remains a challenge [45]. Recently, it has been suggested [46][47][48] that new experimental measurements of quarkonium production inside jets [49] may help better constrain the relevant LDMEs. The focus of this paper is, however, different. Our study concentrates on the production, propagation, and dissociation of quarkonium states in strongly interacting matter. We follow the NRQCD calculation outlined in [32] and use the LDMEs extracted there to give good description of the cross sections for bottomonia for pp and pp collisions for p T in the range of 5 to 30 GeV. For charmonia we improve the χ c fitting procedure by allowing the singlet matrix element as a free parameter. We also refit LDMEs for the ψ(2S) by fitting the LHC 7 TeV and CDF 1.8 TeV data. Both these changes make the spectra for χ c and ψ(2S) in the p T ∼ 10 − 20 GeV region softer and improve the description of data without spoiling the agreement at lower p T .
In order to address quarkonium attenuation in heavy ion reactions, we need to understand the J/ψ and Υ states' behavior and melting at finite temperature in the QGP, and the dissociation processes due to collisional interactions with the quasi-particles of the QCD medium. To accomplish this, a detailed knowledge of the wavefunctions at zero and finite temperatures is necessary. We start by solving the Schrödinger equation by separating the radial and angular parts of the wavefunction, ψ(r) = Y m l (r)R nl (r). The reduced equation for the radial part can be written as where µ red = mQ 2 is the reduced mass, and n and l are the principal and orbital quantum numbers, respectively. V (r) is the potential between the two heavy quarks, which can be estimated from the lattice [7]. The binding energy of the meson is then For T = 0 we take the form of the potential to be of the standard Cornell form The form of the potential has been validated in multiple lattice calculations by evaluating the Polyakov loop correlator as a function of the quark-antiquark separation. The Cornell long-distance part is cut off at r max = 1.1 GeV [9] to model string breaking on this length scale. This defines the value of V (∞) ≈ 1.2 GeV. For the mass of the heavy quark we take m c = 1.34 for the charm quark and m b = 4.5 for the bottom quark. For T > 0 the form of the real part of the potential we use is the internal energy found in [6]. We have checked that the potentials we use are quite close to the internal energies computed in the more recent work [7] and the potentials used in [14]. Internal energies provide a stronger binding potential compared to the single free energies [7]. The use of the internal energies can be justified if the dynamics of the quarkonia are too rapid for the QQ potential to fully thermalize, which might be especially likely for high p T quarkonia in the medium. With this setup, solutions are obtained for a variety of temperatures and for the Swave and P-wave states. Even though we are primarily interested in the J/ψ, ψ(2S), and Υ(nS), the P-wave χ c and χ b contribute via feed down.  l refers to the angular momentum of the QQ state, while n is the radial quantum number. E b nl is the binding energy, r 2 is the root mean square (RMS) radius of the quarkonium state, and k 2 is the mean square momentum. Pertinent results for the quarkonium and bottomonium wavefunctions are presented in Table I and Table II for zero temperature (T ) and T = 192 MeV, respectively. We show the binding energy E b nl and the root mean square (RMS) size r 2 of the quarkonium state. The RMS size changes only slowly with T , except near the dissociation temperature. We have chosen the finite temperature value ∼ 190 MeV to illustrate that several states, i.e. ψ(2S), Υ(3S), and χ b (3P ), cease to exist even with our selected strong binding potential. The χ c and χ b (2P ) states are very near dissociation. If we Fourier transform to momentum space, the mean squared momentum k 2 can also be evaluated and is given in Tables I and II. We note that while there is a correlation between the width of the wavefunctions of quarkonia and their binding energies, this correlation is highly non-linear. The widths change rapidly only near dissociation when E b nl → 0. In this work we are interested in large transverse momentum quarkonia and it is convenient to work in light cone momenta and with light-cone wavefunctions. The relation between the instant form and the light-cone form of the momentum space wavefunctions for mesons was discussed in detail in [31,32]. For the case of quarkonia, the color-singlet contribution can be understood as one matching to the lowest order (n = 2) Fock component of the state. A color-octet initial state must emit at least one gluon for a color neutral hadron to be produced. In either case the heavy meson or proto-quarkonium state of momentum P + = (P + , P) can be approximated as: where a † (b † ) represent an "effective" heavy quark (antiquark) in the 3 (3) state, c 1 , c 2 being the color indices [31,32]. The light cone wavefunction ψ(x, k) in Eq. (4), which describe the longitudinal momentum fraction x (1 − x) distribution and the transverse momentum k (−k) distribution of heavy quarks (antiquarks) is given by In Eq. (5) Λ(T ) is the transverse momentum width of the light-cone wavefunction which needs to be constrained at 0 and finite temperature to the mean transverse momentum squared from the solution to the Schrödinger equation, which remains invariant under boost. If we introduce the notation ∆k = k 1 − k 2 = 2k, the equation for Λ(T ) reads The factor 2/3 comes from the 2D projection of the mean squared transverse momentum κ 2 from the instant-form wavefunction. The temperature and/or density profiles of the medium, which play an important role in the dissociation of quarkonia, can be obtained from hydrodynamic simulations of the QGP [33,50,51] in 2.76 TeV and 5.02 TeV Pb+Pb collisions at the LHC. Specifically, we use the iEBE-VISHNU (2+1)-dimensional event-byevent viscous hydrodynamic package [33]. A sample tem- perature distribution at time t = 2 fm when the interplay between the formation and dissociation of quarkonia is important in setting the final observed J/ψ and Υ multiplicities is shown in Fig. 1. By comparing the temperature in the different points in the (x, y) plane perpendicular to the collision axis to the results in Table II one can get a sense of how the different quarkonium states will be attenuated in heavy ion collisions relative to proton collisions. The propagation of a QQ state in matter is accompanied by collisional interactions mediated at the partonic level, as long as the momentum exchanges between the medium quasi-particles and the heavy quarks can resolve the partonic structure of the meson. The related modification of the quarkonium wavefunction in Eq. (5) can lead to the dissociation of J/ψs and Υs in addition to the thermal effects. The cumulative one dimensional momentum transfer for a quarkonium state that starts at transverse position x 0 and propagates with velocity β, such that x(τ ) = x 0 + β(τ − t 0 ), reads Here µ 2 D = g 2 T 2 (1 + N f /6) is the Debye screening scale and we use 2 active light quark flavors N f = 2. The scattering inverse length of the quark is 1/λ q = σ qq ρ q +σ qg ρ g , where ρ q and ρ g are the partial densities of light quarks and gluons in the QGP. We label the cumulative one dimensional momentum transfer χµ 2 D ξ in analogy with a uniform static medium where the Debye scale is fixed and χ = L/λ q is the opacity. The elastic scattering cross sections are given by Last but not least, ξ is a parameter related to the heavy quark broadening from multiple scattering in the QGP [30,52]. In the limit of strictly soft interactions ξ = 1 and an enhancement of ξ ∼few may arise from the power law tails of the Moliere multiple scattering. We initialize the wavefunction ψ i (∆k, x) of the protoquarkonium QQ state with a width Λ 0 ≡ Λ(T = 0). This is a natural choice since in the absence of a medium it will evolve on the time-scale of O(1fm) or greater into the observed heavy meson. By propagating in the medium this initial wavefunction accumulates transverse momentum broadening χµ 2 D ξ. The probability that this QQ configuration will transition into a final-state heavy meson with thermal wavefunction ψ f (∆k, x) with Λ(T ) is given by In Eq. (9) Norm i is the normalization of the initial state, including the transverse momentum broadening from collisional interactions, and Norm f is the normalization of the final state. The dissociation rate for the specific quarkonium state can then be expressed as It will enter the time evolution of the J/ψs and Υs in the medium. To visualize the dissociation rate, we present 1/τ diss. as a function of time in Fig. 2 for 0-10% central Pb+Pb collisions at the LHC at √ s per nucleon pair of 2.76 TeV. We have weighed the dissociation rates with the binary collision density in the plane perpendicular to the collision axis and averaged over multiple hydrodynamic events. The top panel shows the dissociation rates of charmonium states and the bottom panel shows the corresponding rates for bottomonium state. We find that the modification of the ground states, such as J/ψ, Υ(1S) is dominated by dissociation at the early stages of the evolution of the QGP, whereas excited weakly bound states, such as are ψ(2S), χ b (3P ), can be strongly modified until they escape the medium. Finally, we note that effects that are suppressed at high transverse momentum, such as recombination [53] of unbound Q andQ because of the Boltzmann suppression factor [54], the Cronin effect [55], and power  corrections [56] do not play a role. Cold nuclear matter (CNM) energy loss might affect production cross sections [57], but its effects become significant at very high p T near the kinematic threshold [58]. Last but not least, it was also found that leading-twist shadowing effects near mid-rapidity at high transverse momentum at the LHC are small [59]. While experimental results in p+Pb collisions at the LHC cannot exclude CNM effects at the 5 − 10% level at p T > 5 GeV around mid-rapidity, they are also consistent with no nuclear modification [60][61][62]. For these reasons, for our study at |y| < 2.4, p T > 5 GeV we neglect these effects.

III. PHENOMENOLOGICAL RESULTS
In this section we present the phenomenological results of our theoretical model for quarkonium dissociation due to thermal wavefunction effects and collisional breakup. For every centrality class and hydrodynamic event we distribute the production of the proto-quarkonium states according to the binary collision density in the 2D plane perpendicular to the collision axis. The azimuthal distribution in the directions of quarkonium propagation is uniform and we also average over multiple fluctuating hydrodynamic events that describe the QGP background.
For each quarkonium state the dynamics of production sus ALICE [35] and CMS measurements [38]. Bottom panel: RCP as a function of Npart versus ATLAS measurements [37].
and propagation through the QCD medium described above is given by a set of differential equations: d dt Here we denote by dσ QQ (t; p T )/dp T the cross section to produce the proto-quarkonium states that evolve into an interacting with the medium meson on the time scale of t form. . This time scale for the heavy quarks to interact with the QGP is taken in the relativistic limit.
The initial condition at t ≈ 0 includes the short distance perturbative QQ production cross sections and the long-distance matrix elements for the particular quarkonium state. dσ meson (t; p T )dp T is the cross section for that quarkonium state as a function of time t. Finally, dσ diss. (t; p T )/dp T is the cross section of the dissociated QQ pairs that will not produce quarkonia. In the absence of a medium t diss. → ∞ and the proto-quarkonium states evolve into the corresponding meson with unit probability. The system of equations Eqs. solved for each of the quarkonium states listed in Table I. Feed-down is then performed using the NRQCD cross sections and branching ratios.
We start by first discussing results at the lower energy of √ s = 2.76 TeV. In Fig. 3 we present the centrality dependence of prompt J/ψs in Pb+Pb collisions. The bands reflect the combined uncertainty of the interaction onset time t form. and the collisional dissociation of the quarkonium states. In the evaluation of the latter we keep the coupling between the heavy quarks and the medium fixed at g = 1.85 [32] but vary the broadening parameter ξ. The upper edge of the uncertainty band corresponds to t form. = 1.5 fm, ξ = 1. The lower edge of the uncertainty band corresponds to t form. = 1 fm, ξ = 2. The upper panel of Fig. 3 shows comparison to the AL-ICE [35] and CMS [38] prompt J/ψ measurements. We find improved description in the most central N part. bins relative to the case when thermal screening effects were not included [32]. Around N part. = 100 the data lies on the upper edge of the theoretical error band. The bottom panel of Fig. 3 shows the J/ψ R CP , where the 40%-80% peripheral collisions are used as a baseline. The ATLAS collaboration measured inclusive J/ψ [37]. However, in the p T < 10 GeV interval which dominates the cross section, the non-prompt B → J/ψ contribution is limited to 20-30% [32] and will not noticeably affect the theoretical results.
Recently, experimental results for the differential suppression of the Υ(nS) family have appeared at high p T [39]. Theoretical calculations for the Υ(1S) (red band) and Υ(2S) (blue band) in minimum bias √ s = 2.76 TeV Pb+Pb reactions are shown in Fig. 4. We have evaluated the cross sections for quarkonia in 10 centrality classes (labeled i) and The experimental data is described well, including its magnitude and p T dependence. We note that collisional dissociation mostly affects the ground Υ state, while thermal wavefunction effects dominate the attenuation pattern of the excited Υ states. The CMS collaboration also put an upper limit on the Υ(3S) cross section in Pb+Pb reactions, corresponding to an upper limit on its R AA [39]. Our calculated Υ(3S) cross section is consistent with this limit. While the theoretical approach presented in this Letter is applicable at large transverse momenta, we observe in Fig. 4 that the nuclear modification factor is approximately constant. This allows us to compare in Fig. 5 the centrality dependence in the lowest p T = 5 GeV bin, not very different form the mean p T of bottomonia at the LHC, to the experimentally measured Υ(1S) and Υ(2S) R AA dependence on the number of participants [39]. Very good agreement between data and theory is observed. In Figs. 6 and 7, we present theoretical predictions for the R AA of various quarkonium species as a function of p T in Pb+Pb collisions at √ s = 5.02 TeV. The top and bottom panels display results for charmonium and bottomonium states, respectively. We find a clear separation in suppression based on how tightly bound the quarkonium state is. We also find a flat or slightly increasing R AA with p T . By comparing Fig. 6 to Fig. 7, we observe that the attenuation of quarkonia in minimum bias collisions is only slightly smaller than in the most central collisions. The reason for that behavior is that minimum bias collisions are strongly dominated by the first 3 most central classes, as given by the weights W i in Eq. (14).
Very recently, at √ s = 5.02 TeV, measurements of relative suppression ratios of excited to ground quarkonium states have appeared [40,63]. The data for ψ(2S)/J/ψ is publicly available and shown in Fig. 8. Theoretically, the double suppression ratio can be obtained from the results in Fig. 7 and is compatible with the experimental data within the statistical and systematic error bars.

IV. CONCLUSIONS
In summary, we presented theoretical results for the p T -differential suppression of charmonia and bottomonia in Pb+Pb collisions at the LHC. The dynamics of QQ pairs, which evolve into the observed quarkonium states, is governed in HICs by the formation and dissociation time scales. A key element of our formalism that addresses this dynamics is that the formation time of protoquarkonia is ∼ 1 fm. We assume that this time scale is long enough that the QQ interact via a color screened thermal potential [6,9]. Therefore we employ the wavefunctions obtained by solving the Schrödinger equation for QQ interacting via a screened potential to calculate the dissociation time scale, using the theoretical setup described in [32]. The technical advances that are further incorporated in the calculation are better constraints on the NRQCD matrix elements that are relevant for the production of high-p T ψ(2S) and χ c states, and their feed-down to J/ψ, and 2+1 dimensional event-by-event hydrodynamic modeling of the QGP background [33].
We explored the phenomenological implications of this theoretical model for quarkonium production and propagation in the QGP background created in heavy ion collisions, first around mid-rapidity in Pb+Pb collisions at the LHC. We found good separation in the magnitude of the suppression between the ground and excited charmonium and bottomonium states, compatible with recent experimental measurements. Our results indicate that effects of thermal screening of the QQ attractive potential fully dominate the attenuation of ψ(2S), Υ(2S) and Υ(3S). On the other hand, J/ψ and Υ(1S) are also sensitive to the dissociation processes due to collisional interactions. The approximately constant or slightly decreasing R AA with p T predicted by this model arises from the early O(1 fm) formation of the interacting quarkonium state. The uncertainty of the phenomenological results was estimated by varying the formation time and the strength of the collisional broadening of the QQ pair. We found that the charmonium suppression measurements are better described by the upper edge of the R AA uncertainty band, whereas bottomonium suppression measurements are better described by its lower edge. While an illustrative subset of results was presented here, detailed predictions are available that will allow to test this model versus upcoming experimental measurements of quarkonium suppression at √ s = 5.02 TeV.
In the future, we plan to address data at finite rapidity with the same parameters and test the model further. This will likely require inclusion of CNM effects since p−Pb data at finite rapidity seems to show non-trivial nuclear modification patterns. It will also be interesting and instructive to investigate non-prompt J/ψ production, which probes the complementary but different physics of in-medium modification of heavy-quark parton showers [64,65].