Bottomonium suppression in PbPb collision at energies available at the CERN large hadron collider

We had been gradually working towards building a comprehensive quarkonia suppression formalism to explain all 3 dependencies of quarkonium suppression obtained from heavy-ion collision experiments. We present here the improved version of quarkonia suppression framework. It assumes bottomonia produced in the early stage which dissociates due to color screening, gluonic dissociation, and collisional damping in addition to the shadowing as an initial state effect. The QGP medium formed in the collisions is assumed to evolve under ($3+1$)-dimensional relativistic viscous hydrodynamics which is modeled using ECHO-QGP. This replaces the Bjorken's hydrodynamics which we had used in our earlier work where we determined the centrality and transverse momentum dependent suppression. The correlated bottom quark and bottom anti-quark could recombine in the plasma. A rate equation is employed, whose solution gives the final number of bottomonium after dissociation and recombination under $(3+1)$-dimensional expansion of the QGP medium. The Shadowing effect, which is the dominant Cold Nuclear Matter effect at LHC energies, has now been modified by employing the most recent parton distribution functions obtained from CT14 global analysis and shadowing factors from EPPS16. Using this improved formalism we determine the centrality, transverse momentum, and rapidity, dependencies of bottomonium suppression for $\Upsilon(1S)$, and $\Upsilon(2S)$ at the LHC's energies of $2.76$ TeV and $5.02$ TeV. We find a fairly good agreement between theoretically calculated survival probability and the measured nuclear modification factor($R_{AA}$) at the two energies.


I. INTRODUCTION
We have made considerable progress in our understanding of the strong interaction since the advent of heavy-ion collision experiments especially at Relativistic Heavy-ion Collider(RHIC) and at Large Hadron Collider (LHC). Quantum Chromodynamics (QCD), the theory of strong interaction predicts that a thermalized medium of Quark-Gluon Plasma (QGP) should form as a transient stage of heavy-ion collisions when the temperature and/or matter density of the fireball exceeds a certain threshold value. The largest portion of the speculative phase diagram of the QCD is occupied by this QGP [1,2]. It is pertinent for us to investigate the nature of this medium in order to make advancement in our understanding of QCD and its limits. And with consensus that QGP does exist in ultra-relativistic energy nucleus-nucleus collisions, the focus is now shifting towards characterizing this medium with the help of the probes which acted as its signatures. A good starting point is to use the wide range of available experimental data from different systems of collisions and energies to constrain the input parameters of the dynamics at work.
So far, the experimental findings indicate the formation of an inviscid medium [3,4]. A theoretical lower limit for the shear viscosity to entropy density ratio(η/s) of strongly interacting matter has been know for a while [5]. But the estimation of such low η/s for quark matter in heavy ion collision has become possible only recently using lattice QCD [6] and bayesian parameter estimation methods [7]. Hydrodynamics has been widely applied to simulated heavy ion collisions of various systems and energies. Phenomenological models based on hydrodynamics better explain the observables from high energy collisions with thermalization started early, τ<1fm/c [8,9]. What makes this problem of characterizing QGP so difficult is the sheer complexity of the system. In high energy collisions, the only variables in our control are particle species and the energy of collision. All the rest has to be inferred indirectly from the observables. Predicting the transport and thermodynamic properties of a medium formed in the collision by analysing patterns in the produced particle yields is a scrupulous task. Hence, modelling of such a complex system can only be dealt phenomenologically [10]. Right after the collision, we get a fireball where the medium consisting of quarks and gluons expands against the surrounding vacuum and cools down rapidly. We assume that this system thermalizes quickly after the collision and we mark the time required for the system to thermalize as QGP formation time. Hydrodynamics is switched off when the system temperature falls below the pseudo critical temperature(T c ) of the QGP. The recently agreed upon value of T c for QGP computed by Lattice QCD collaborations is 156 MeV [11][12][13]. One way of testing a given model is by calculating a physical quantity which could act as a theoretical counterpart of an observable measured by detectors. If QGP does exist as a transient stage of heavyion collision, then we should be able to notice an agreement between these two quantities [14]. Among many of such signatures of QGP, one is quarkonia suppression, on which we are concentrating in the current work. Quarkonia are mesonic bound states of heavy quark and heavy anti quark, which are produced in the early stage of collision. They could dissociate due to various types of interactions with the partons in the medium and would be detected comparatively lesser in number than in collision systems where we do not expect QGP, like at low energies and in p-p collisions [15]. In order to quantify this suppressed production of quarkonium, experimentalists measure a physical quantity called nuclear modification factor(R AA ). It is the quarkonia yield in heavy-ion collision divided by the yield of the quarkonium in p-p collision scaled by N coll . Its value less than one, greater than one and equal to one indicates suppression, enhancement and no medium effect, respectively. Surprisingly, there are cases which contribute to suppression in collision systems where we do not expect a thermalized medium. These non-QGP effects arise due to situations before collisions (initial state effects where the system is said to be cold), even though small needs to be modelled into the total suppression scheme [16]. These non-QGP effects are called cold nuclear matter (CNM) effects. The R AA is measured over a wide range of collision energies as a function of centrality, N part , transverse momentum, p T and rapidity, y. Various phenomenological models have tried to explain the measured values of the suppression consistently over a wide range of beam center of mass energies and collision systems [17]. Few of them attempted explaining centrality and p T dependencies of [18,19]. And even fewer predict all the three dependencies of the suppression over a wide range of available center of mass energies [20]. In our earlier work published in Eur. Phys. J C 79, 147 (2019), we had explained the p T and N part dependence of R AA over a range of LHC energies. It was based on the suppression due to colour screening, gluonic dissociation and collisional damping under (1 + 1)-dimensional Bjorken's expansion of the thermalized medium. The net quarkonium yield was determined using a rate equation which combines suppression and recombination due to correlated quark anti quark pairs.
In the current study, we start with the initially produced bottomonia yields which evolves in (3 + 1)dimensional hydrodynamic medium.
This initial yield is influenced by a CNM effect used here called "Shadowing", which has been updated for the newly available parton distribution functions and shadowing factors. The bottomonium bound state while drifting through QGP could dissociate due to gluonic dissociation, collisional damping and color screening which are adapted to the (3 + 1)-dimensional hydrodynamic expansion. Color screening has been streamlined by eliminating the need of assuming a pressure profile of collision. A lattice QCD based equation of state(EOS) from Wuppertal-Budapest Collaboration has been utilized. The input parameters for the hydrodynamics are constrained using the transverse momenta and rapidity spectra for poins from the ALICE experiment. And lastly after considering the possibility of recombination of correlated bottom quark and anti-quark pair, we find the final number of bottomonia for ground state and excited states [21]. We then find a quantity called survival probability(S p ) which is theoretically equivalent to the experimentally measured R AA . We determine this as a function of transverse momentum, centrality and rapidity at 2.76 TeV and 5.02 TeV energies and then compare with the corresponding R AA values for Υ(1S) and Υ(2S) states. On comparison, we find a reasonably well agreement between S p and R AA at two LHC energies.
The arrangement of the topics in the current paper is as follows. Sec. I describes the general introduction of the proposed work. In Sec. II, we briefly describe the (3 + 1)-dimensional hydrodynamics used to model the bulk of medium using ECHO-QGP. It is followed by the quarkonia suppression formalism and the various effects incorporated into the scheme. Sec. III presents the results and discussions on the suppression of two bottomonium states, (1S) and (2S) versus N part , p T , and rapidity y at the two LHC energies. Finally in Sec. IV, we summarize our results obtained and conclude the work.

II. FORMALISM
Here we describe the formalism in brief. The suppression formalism which has been developed in our previous work [22], has been adapted for the (3 + 1)-dimensional viscous hydrodynamics. More details about individual medium effects used can be found in [23][24][25].

A. (3+1)-dimensional Hydrodynamical Expansion of the Medium
Hydrodynamics has been quite successful in explaining bulk observables from ultra-relativistic heavy-ion collisions for a wide range of system and energies [26][27][28][29]. The agreement of hydrodynamical predictions with experimental results has been taken as an indirect evidence for the correctness of dynamics. The Bjorken's hydrodynamics assumes that the fireball expands only along the longitudinal direction and is restricted to the plateau region of rapidity spectra (dN/dy vs y). This leads to the key variables of dynamics, e.g. temperature, pressure, energy density, entropy density to become an explicit function of proper time [30]. Thus, although models based on Bjorken's evolution are adequate to estimate the observables at midrapidity, yet they are ineffective in providing the complete spacetime evolution of the system. In order to simulate the true dynamics of a collision that holds up expansion along the transverse directions and larger rapidities, one has to switch to the complete (3 + 1)-dimensional hydrodynamics. ECHO-QGP is a FORTRAN based code to find the solutions of the conservation equations, d µ N µ = 0 and d µ T µν = 0, where, d µ is the covariant derivative, N µ is the four current and T µν is the energy momentum tensor. To solve these equations numerically in (3+1)dimensions with relativistic speeds and viscous conditions, Israel-Stewart's second order formalism has been used in ECHO-QGP [31]. A Cooper-Frye prescription handles the freeze-out stage, where the produced particles are assigned their momenta at the constant temperature hypersurface. We vary the input parameters of the ECHO-QGP so that the particle momentum spectra calculated here matches with the measured spectra from experiments as explained below. We mark the end of the QGP phase at the value of proper time when the maximum temperature of the system drops below pseudo critical temperature, T c .
The equation of state (EOS) from Wuppertal-Budapest(WB) Collaboration [32] replaces the earlier quasiparticle EOS [33]. WB EOS is spline interpolation with the hadron Resonance Gas (HRG) EOS [34] for the hot and dense hadronic matter after hadronization. WB EOS computed from the lattice QCD is a better choice as the pseudo critical temperature range of QGP predicted by their analysis lie close to the presently agreed upon value [35]. We ran the hydrodynamics code for 11 values of the impact parameter covering 0 − 100% centrality range. We chose the geometric Glauber initialization in ECHO-QGP [36]. A rapidity profile of p−p collisions is also employed as an input. The two parameters characterizing this profile are ∆ s ; which is the extension of the rapidity plateau and σ η ; which is the width of the Gaussian falloff of the profile. Values for both of these parameters are varied until the shape of pions (π + ) rapidity spectra matches with that from the experimental poins rapidity spectra for the two mentioned LHC energies as shown in Fig(1) [37,38]. The values of the relaxation time coefficient for viscosity of second order, τ π and the shear viscosity to entropy density ratio, η/s are taken from [39,40]. The thermalization time in the code is set for both the energies viz. 2.76 TeV and 5.02 TeV at 0.20 fm/c [26,39,41]. The inelastic nucleon-nucleon crosssection is taken to be 61.8 mb and 70 mb for 2.76 TeV and 5.02 TeV, respectively [36]. Lattice QCD predicts formation of a thermalized medium at energy density above 1.0 GeV/fm 3 [42]. Initial energy density which goes as an input in ECHO-QGP was at first calculated roughly using an approximate relation ǫ 0 = 1 A T τ 0 J(y, η) dE T dy involving overlap area, initial thermalization time and the differential transverse energy [30,43]. But the peak values of the pion p T -spectra for these values fall shorter than the experimental values. Hence, we varied the initial energy density at each centrality such that the pions p T -spectra from ECHO-QGP matches with that from experiment values [44,45]. The comparison of the pion spectras at few of these centralities is shown in Fig(2).
The key parameters used in the ECHO-QGP hydrodynamics are summarized in Table I.
The suppression formalism, described in the next section, requires temperature of the medium at different centralities and rapidities as a function of proper time  in the transverse plane, which are obtained from the ECHO-QGP. Calculation of suppression at all transverse (x,y) points is computationally infeasible. Hence, temperatures are integrated over the transverse plane at all centrality, rapidity and proper time with a Gaussian weight factor and are taken as an input in the suppression formalism. The standard deviation of this Gaussian profile is varied within a specific range which is explained in the results and discussion section.

B. Suppression formalism
Bottom quark anti-quark pair produced in the hard scattering of colliding nuclei in the early stage of collision combines to form the bottomonia mesons. These heavy flavor mesons will drift in the medium and their decay products are eventually detected. During their time in the medium, these meson bound states are affected by various medium dependant dissociation mechanisms like color screening, collisional damping and gluonic dissociation, which are individually explained below along with the possible recombination due to correlated b −b pairs and the non-medium effect of shadowing.

Gluonic Dissociation
Gluonic dissociation is referred to the process where a bottomonium color singlet state absorbs a soft gluon in the medium and gets excited to a color octet state. The cross-section for this process is calculated as [46]; where, E g is the soft gluon energy, E nl is the eigenvalues corresponding to the bottomonium wavefunction (g nl (r)), m q is bottom quark mass in GeV, N c is the number of color charges, α u Above equation gives the probability density, where g * nl and h qi ′ (r) are the singlet and octet wavefunction of bottomonium respectively, obtained after numerically solving the 3-dimensional Schrödinger's equation. We integrate the cross-section in equation (1) with Bose-Einstein distribution as weight factor over gluon momentum to calculate the dissociation width due to gluonic dissociation i.e., Γ gd,nl . The validity of the above cross-section assumes T≪ 1/r, where T is the medium temperature and r is the distance between quark and antiquark [47]. This regime is different than the one for which collisional damping holds, due to which our final results for T≪ 1/r may have less accuracy.

Collisional Damping
We expect bottomonium to dissociate while it traverses through the plasma, due to the momentum transfer arising out of collision. To account for this effect, we use a potential non-relativistic QCD (pNRQCD) formalism which depends on the imaginary part of the color potential between quark anti-quark pair. The complex potential between quark anti-quark pair located inside the QGP medium as determined by Laine et al., [48] using effective field theory, is given as; Where, α eff = 4α s s /3, σ is the string tension, whose value is 0.192 GeV 2 , m D is the Debye mass which is expressed as; Here, N c and N f are the number of color charges and number of flavors, respectively. This potential is valid in the regime T≫ 1/r m D [49]. We take the expectation value of the imaginary part of this potential to get the dissociation width corresponding to collisional damping [50] as; where, g nl (r) is the singlet wavefunction of the bottomonium.
In our previous work [22], the explicit τ dependence of gluonic dissociation widths and collisional damping widths(4) arose from the analytical expression for temperature, T(τ), obtained by solving Bjorken's dynamics.
We have now replaced it by the tabulated temperature values from ECHO-QGP at each centrality and rapidity integrated over the transverse plane for each of the bottomonium states.

Color Screening
Due to the presence of free color charges in QGP medium, the bottom quark and anti-quark (bandb) experiences a short range Yukawa-type color charge potential dependent on the medium temperature. As a consequence of this the formation of bound states in the medium is suppressed if the medium temperature goes beyond a certain temperature (dissociation temperature). This phenomenon is called the color screening [15] in an analogous way to the Debye charge screening in the quantum electrodynamics plasma. Different bottomonia species take different time, after collision, to form the corresponding bottom anti-bottom bound states termed as their respective formation times(τ f ). Consider a screening region in the fireball as a sphere with a screening radius(r s ) defined for a bottomonium state depending on its dissociation temperature(T D ). Suppose a bb quark pair forms at the position r Q , then it will likely form a bound state if it escapes the screening region in time equal to its formation time. The condition for which could be given by | r Q + v T τ F | ≥ r s , where v T is the transverse drift velocity of the heavy quark in the medium. Here τ F is the vacuum formation time of bottomonium.
Considering the transverse motion of this bottomonium state, we can simplify the above condition as; Here, m Q is the mass of the bottomonium state under consideration. Let us now consider a transverse radial distribution of bb produced after the hard scattering as: Where R T is the maximum fireball radius obtained for different centralities from Modified Glauber analysis [51].
The color screening survival probability is given by: Where, φ max (r) is the maximum positive azimuthal angle allowed by the condition in Eq(5).
In our previous work [22], we equated an assumed pressure profile in the transverse plane with the cooling law to obtain the screening time (time to drop initial pressure to the pressure corresponding to the dissociation temperature). Screening time is then equated to the bound state formation time (at the boundary of the screening region) to obtain the screening radii [25]. Here we eliminated the need of this assumption of pressure profile in the transverse plane by directly finding the screening radii from the transverse temperature evolution using ECHO-QGP. We take the minimum radius of the oval shaped contour shown in the Fig(3) as the screening radius at the evolution time equal to the formation time for a given bottomonium species. The temperature contours are marked for the dissociation temperature of all the bottomonium states at all centralities. The dissociation temperature (T D ) of different bottomonium states are borrowed from the analysis given in [52].

Shadowing
The shadowing correction to R AA applied in our formalism is a modified version of the similar work by R. Vogt [53]. We have replaced the shadowing factors used for gluons from EPS09 [54] to the more recent EPPS16 [55]. The central fit set is selected out of various available error sets in EPPS16. The parton distribution function of gluons have been updated to CT14 [56] from the earlier PDFs CTEQ6 [57]. The contribution of suppression arising due to the shadowing effect is expressed as [58]; Shadowing effect influences the number of initially produced bottomonia(N Q ). Hence, the shadowing corrected initial number of bottomonium is calculated as,N i Q (τ 0 , b) = N Q (τ 0 , b)S sh (p T , b).

C. Recombination mechanisms
We have incorporated the possibility of recombination of bb due to de-excitation from octet to singlet state with a gluon emission, even though it will be negligible for the case of bottomonium. We find the recombination cross-section in QGP using detailed balance from gluonic dissociation cross-section as [59] where, s is the Mandelstan variable, M nl , m q and mq are the masses of bottomonia states, bottom quark and bottom anti-quark, respectively. We then define a recombination factor as the thermal average of the product of the above cross-section and relative velocity between b andb as, Γ F,nl = σ f ,nl v rel k .

D. Final Number of Bottomonium
Due to all of the above effects, the bottomonia can dissociate or the correlated bb pair can recombine again into bound states. We assume that this interplay of dissociation and recombination is governed by a simple first order differential equation given as [60]; Here, N Q is the bottomonia yield at a given value of proper time (τ). First and second terms on the right hand side of this equation correspond to the recombination and dissociation terms, respectively. Γ F,nl and Γ D,nl are the corresponding recombination and dissociation rates. N q and Nq are the number of heavy quark and anti-quark produced in p−p collision. V(τ) is the instantaneous volume of the expanding fireball. The solution for the above first order differential equation under the approximation that, N Q < N q ,Nq is given by, where τ QGP is the QGP lifetime and τ 0 is the initial time at which we start hydrodynamics and which also marks the beginning of the QGP stage. The first term inside the bracket on the right hand side in Equation (11) is the bottomonia produced at the initial hard scattering stage. N Q (τ 0 ) is calculated as [61]; where, σ NN Q is the production cross-section of the bottomonium at a given collisional energy. T AA (b) is the nuclear overlap function. V(τ) in equation (11) is volume of fireball given by the formula; where, R T is the radius of colliding nuclei and R is the Reynold's number. ǫ(τ, p T ) in equation (11) is a suppression factor integrated cumulatively in τ values and is given by; Equation (11) also has multiplicative suppression factor which is integrated for the complete QGP lifetime and is calculated as; (15) In the above equation, τ ′ nl is the time required for the formation of a given bottomonium state within QGP.
We then calculate the ratio, S ′ p = N Q /N i Q referred to as the survival probability. The color screening has been considered as an independent effect and hence the total survival probability is calculated as S p = S ′ p × S cs . We find this survival probability for 5 bottomonia states which are Υ(1S), Υ(2S), χ b0 (1P), χ b0 (2P) and Υ(3S). A feed-down scheme dictates the total yield after all possible decays from higher excited states [23]. The Survival probability(S p ) obtained after feed down for Υ(1S) and Υ(2S) is plotted and compared with the respective R AA obtained from experiments.

III. RESULTS AND DISCUSSIONS
Results have been generated for centrality range of 0 − 100%, transverse momentum range of 1 − 30 GeV/c and rapidity range of |y| < 3.6 which covers the ranges of experimentally available data from CMS and ALICE at both 2.76 TeV and 5.02 TeV [62][63][64][65].
ALICE suppression data at forward rapidity complements well with the broader rapidity range of CMS data especially for p T and N part dependence of Υ(2S) at √ s NN = 5.02 TeV. The standard deviation (σ T ) of the Gaussian weight factor used for integrating temperatures from ECHO-QGP as mentioned in Section (II-A) has been varied to obtain an uncertainty patch in the theoretically calculated suppression results as shown in the figures. For 2.76 TeV the standard deviation value lies between 0.7 < σ T < 2.8 and for 5.02 TeV, it spans as 1.4 < σ T < 3.9. Selecting σ T greater than the upper limits for the two energies makes the QGP lifetime, in peripheral collisions at extreme rapidity values, smaller than the formation time of Υ(1S) and Υ(2S) states. Selecting σ T smaller than the lower limits undermines the temperatures at large x-y values rendering the purpose of integration futile. Fig.(4) depicts the centrality dependence of suppression for Υ(1S) and Υ(2S) states at 2.76 TeV as calculated by our present model. The corresponding experimental suppression data are shown for comparison. Our calculated values of the survival probability for Υ(1S) lie very close to the CMS data and are following the trend of AL-ICE data. Whereas for Υ(2S), our predicted values are slightly less suppressed but mostly following the CMS data. Fig.(5) shows the variation of our predicted values of suppression for Υ(1S) and Υ(2S) with respect to the transverse momentum at the 2.76 TeV center of mass energy. We find that the agreement among our calculated and measured values for Υ(1S) are reasonably well. However, for Υ(2S), calculated suppression is slightly lesser than the values from available data from CMS.
In Fig.(6), we have plotted our theoretical results of rapidity dependence of Υ(1S) and Υ(2S) suppression along with the corresponding experimental data. We find a quite reasonable agreement in the case of Υ(1S) than Υ(2S) within the uncertainty limit. Our predicted values show lesser suppression for Υ(2S) as compared with the corresponding experimental data. This rapidity dependence of suppression has become possible due to the interfacing of our earlier model with the ECHO-QGP's hydrodynamic expansion. Overall, for 2.76 TeV LHC center of mass energy, we find good agreement between our calculated bottomonium suppression values and the corresponding experimentally available data under the theoretical and measured uncertainty limit. Fig.(7) depicts the variation of survival probability for both the bottomonium states with respect to the centrality at 5.02 TeV center of mass energy. The measured values of suppression from ALICE and CMS data shown on the same plot for comparison for Υ(1S) associated with the centrality dependence are quite close to each other. We clearly see a quite good agreement from the Fig.(7) between our predicted values and the measured ones for both the bottomonium states over the whole range of centrality. Since at 5.02 TeV LHC energy.
Transverse momentum dependence of the survival probability values are plotted in Fig.(8) along with the measured values of R AA for both the bottomonium states at 5.02 TeV, center of mass energy. The agreement for Υ(1S) state is reasonably well especially for low p T values, whereas for Υ(2S), the predicted values are slightly less suppressed as compared to the CMS data.
Finally, rapidity dependent survival probability at 5.02 TeV energy is shown in Fig.(9) for both the bottomonium states. In comparison with the corresponding experimental data, our results are close to the measured values for Υ(2S). For Υ(1S), our calculated suppression is not following the trend of ALICE and CMS data. We saw similar disagreement for the corresponding rapidity results at 2.76 TeV energy. These results might improve with further refining of the input parameters taken in ECHO-QGP.
Thus a complete dependence of bottomonium suppression for both the states spanning over two LHC energies has shown a convincing agreement with the measured values.

IV. SUMMARY AND CONCLUSION
To summarize, we have used a quarkonia suppression formalism to explain bottomonium suppression data at 2.76 TeV and 5.02 TeV LHC energies. ECHO-QGP has allowed us to find (3 + 1)-dimensional evolution of the relevant physical quantities associated with the medium formed just after the collisions. The temperatures at different centralities and rapidities have been extracted from the ECHO-QGP and fed into the suppression formalism. This facilitated us to include the rapidity dependence of suppression in the formalism. The experimental data on transverse momentum(p T ) and rapidity(y) spectra has been utilized to fix the input parameters of the hydrodynamics. We have used here the EOS from the Wuppertal-Budapest collaboration which is computed from first principle Lattice QCD. Modified color screening formalism has been used with the newer dynamics eliminating the need to assume a pressure profile for collisions in the transverse plane. Shadowing effect has been updated with the recent gluon PDFs and shadowing factors. The uncertainty in temperature integration from ECHO-QGP has translated into an uncertainty patch in the final suppression results. This formalism has adequately explained the centrality and transverse momentum dependent suppression of Υ(1S) and Υ(2S) at both the LHC energies under consideration. Although, the calculated rapidity dependence of suppression does not follow the trend of data for Υ(1S), yet is quite satisfactory for Υ(2S) state at both the energies. Looking ahead, it will be interesting to do an open heavy flavor mesons evolution study in the QGP medium with a formalism suited well for such heavy-light mesons. A simultaneous study of open and hidden heavy flavor and will give a more reliable constraint on the model parameters since it will account for suppression for all heavy quarks bound states produced in the collisions.