Predictions for bottomonia suppression in 5.023 TeV Pb-Pb collisions

We compute the suppression of the bottomonia states Upsilon(1S), Upsilon(2S), Upsilon(3S), chi_b(1P), chi_b(2P), and chi_b(3P) states in Large Hadron Collider (LHC) sqrt(s_NN)) = 5.023 TeV Pb-Pb collisions. For the background evolution we use 3+1d anisotropic hydrodynamics with conditions extrapolated from sqrt(s_NN) = 2.76 TeV and we self-consistently compute bottomonia decay rates including non-equilibrium corrections to the interaction potential. For our final results, we make predictions for R_AA as function of centrality, rapidity, and p_T for the Upsilon(1S) and Upsilon(2S) states, including feed down effects. In order to assess the dependence on some of the model assumptions, we vary the shear viscosity-to-entropy density ratio, 4 pi eta/s = (1, 2, 3), and the initial momentum-space anisotropy parameter, xi_0 = (0, 10, 50), while holding the total light hadron multiplicity fixed.


Introduction
Ultra-relativistic heavy ion collisions (URHICs) carried out at the Large Hadron Collider (LHC) at CERN and the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory create super-hot and dense matter. The observed thermal particle spectra and anisotropic collective flow seen in such collisions has been taken as evidence of the creation of a deconfined high-temperature quark-gluon plasma (QGP). The existence of the QGP itself is a prediction of finite-temperature quantum chromodynamics (QCD), and the QGP is believed to have been the state of the early universe.
Comparisons between theory and experiment suggest that LHC URHICs produce a QGP with an initial temperature on the order of T 0 = 500−600 MeV for √ s NN = 2.76 TeV Pb-Pb collisions [1,2]. The LHC has recently performed URHICs with nucleon-nucleon center-of-mass energies of √ s NN = 5.023 TeV, which are expected to produce initial temperatures on the order of T 0 = 600−700 MeV. Analysis of collective flow data indicate that the QGP behaves like a nearly-perfect fluid with shear viscosity-to-entropy density ratio η/s approaching the conjectured lower bound of 1/(4π). Despite being nearly ideal, one complication that must be faced is that, although the QGP is close to thermal equilibrium at late times, URHICs produce a QGP which is momentum-space anisotropic in the local rest frame with the degree of momentum-space anisotropy decreasing slowly as a function of proper time. This can have an impact on the various signatures for QGP formation, such as heavy quarkonium suppression [3,4]. Heavy quark-bound states are of particular interest because they can survive into the deconfined phase due to their large binding energy. They are expected to survive up to a few times the pseudo-critical temperature, T c ; however, at sufficiently high temperatures, even heavy quark bound states should disassociate. As a result, in the late 1980's Karsch, Matsui, Mehr, and Satz predicted that heavy quark-antiquark bound states, known as quarkonia, would be suppressed in heavy ion collisions compared to proton-proton (p-p) collisions [5,6]. The resulting suppression of heavy quarkonia can be used to probe the entire history of the QGP, since heavy quarkonium states are produced early in the QGP evolution. There has been a significant amount of work dedicated to the study of heavy quarkonia suppression. For recent related work, see References [7][8][9][10]. In this paper, we focus on bottomonia states which survive up to plasma temperatures of T d ∼ 600, 230, and 170 MeV for Υ(1S), Υ(2S), and Υ(3S) states, respectively [11]. At these temperatures, the in-medium width of the state becomes on the order of the real part of its binding energy, and the state quickly dissociates; however, prior to this disassociation point, the states have large widths (∼50-100 MeV), with the widths for the excited states being larger than those for the ground state. As a consequence, one expects sequential melting of heavy bottomonia states in QGP, with excited states melting before the ground state, and lighter states dissociating before heavier states [12].
For the √ s NN = 5.023 TeV initial conditions, we extrapolate the initial conditions used at lower collision energies and make predictions for the inclusive R AA for the Υ(1S) and Υ(2S) states, including feed down contributions. To assess the dependence on some of the model assumptions, we vary the shear viscosity-to-entropy density ratio, 4πη/s ∈ {1, 2, 3}, and the initial momentum-space anisotropy parameter, ξ 0 ∈ {0, 10, 50}, while holding the total light hadron multiplicity fixed.

3 + 1d Anisotropic Hydrodynamics
In aHydro, the local rest frame (LRF) partonic distribution is assumed to be of the form [27,28,32] where −1 ≤ ξ(x) < ∞ is the local momentum-space anisotropy parameter, Λ(x) is the local transverse temperature of the plasma, and f eq is an arbitrary isotropic equilibrium distribution function which here we take to be a Boltzmann distribution. The local anisotropy ξ(x) encodes the viscous corrections to an ideal isotropic QGP. Note that, in the most recent applications of aHydro, one allows for the existence of three independent anisotropy parameters; see e.g., Reference [30]. This takes into account both bulk and shear viscous effects in the QGP. In this paper, we use the simpler "one anisotropy parameter" version. This is sufficient for our purposes, since the additional anisotropy parameters are small and, therefore, subleading. For details concerning the aHydro dynamical equations and the numerical methods used, see References [33,34]. We start the hydrodynamic evolution of the system at τ 0 = 0.3 fm/c. For the initial conditions, we use the Glauber model with a linear combination of wounded-nucleon and binary collisions. The fraction of the binary collision component is taken to be κ binary = 0.145. The inelastic cross-section for √ s NN = 2.76 TeV collisions is taken to be σ 2.76 TeV NN = 62 mb, while for √ s NN = 5.023 TeV collisions, we extrapolate from the results presented in Reference [35] to obtain σ 5.023 TeV NN = 67 mb. In both case,s we use a standard Woods-Saxon profile for the incoming nuclei with n A (r) = n 0 /(1 + e (r−R)/d ), where n 0 = 0.17 fm −3 is the central nucleon density and is determined via the normalization lim A→∞ d 3 r n A (r) = A, R = (1.12A 1/3 − 0.86A −1/3 ) fm is the nuclear radius, and d = 0.54 fm is the skin depth of the nucleus [36].
The initial longitudinal profile is approximated by distribution featuring a nearly boost-invariant plateau region for mid-rapidities and a tail characterized by a Gaussian which reflects limiting fragmentation [37] For √ s NN = 2.76 TeV collisions, fits were made to experimental particle pseudorapidity profiles which give ∆ς 2.76 TeV = 2.5 and σ 2.76 TeV ς = 1.4. For √ s NN = 5.023 TeV collisions, we use the predictions of Reference [38], which determined the rapidity profile in the high energy limit. Using the results of Reference [38], one finds that going from 2.76 TeV collisions to 5.023 TeV collisions, there is an approximately 12% increase in plateau halfwidth, which gives ∆ς 5.023 TeV = 2.8. Based on [38], we found no observable change of the Gaussian halfwidth in 5.023 TeV collisions, and therefore we take The combination of Glauber in the transverse plane and the longitudinal profile function above gives us the full 3d initial energy density profile for the QGP.
The initial plasma temperatures, T 0 , for 2.76 TeV collisions are taken from our previous work [24] and are consistent with the analysis of elliptic flow coefficients [39]. In Table 1, we show the initial conditions used at 2.76 TeV. The values in the center of Tables 1 and 2 are the initial transverse temperatures of the plasma, which enter into the one-particle distribution function, Equation (1). Note that unless ξ 0 = 0, the transverse temperature cannot be interpreted as temperature. The entry with ξ 0 = 0 and 4πη/s = 1 was obtained via fits to collective flow, and the other entries were obtained by varying either ξ 0 or η/s while holding the total particle multiplicity fixed. Note that, when ξ 0 = 0, we set the initial condition for the transverse temperature Λ 0 instead of the initial temperature due to the momentum-space anisotropy.  For 5.023 TeV collisions, we assume that the temperature of the plasma scales proportionally with the fourth root of the collision energy-i.e., T 0 ∝ s 1/8 NN . Using this assumption, we predict a 16% increase in the initial central temperature when going from 2.76 TeV to 5.023 TeV collisions. We apply this scaling to the 2.76 TeV initial conditions with ξ 0 = 0 and 4πη/s = 1 and then, as before, we fill out the rest of the initial conditions table by holding the total particle multiplicity fixed. The resulting table of initial conditions is shown in Table 2.

Model Potential
Modeling bottomonia states in a QGP which is anisotropic in momentum-space requires going beyond the Karsch, Matsui, Mehr, and Satz model which describes the free energy of a static heavy quark-antiquark pair in an isotropic plasma. It has been shown in References [7,25] that free-energy-based potential models of the quarkonium pair do not agree with experimental R AA data.
For this reason, we use the internal energy of the bottomonia pair U = F + TS to provide the model potential. The masses of heavy quarks allow bottomonia pairs to be treated using pNRQCD methods which allow us to systematically include relativistic corrections. We use a model which represents the short-and medium-range gluonic screening of the heavy-quark potential in a momentum-space anisotropic plasma [19,21,22,25,26]. For the long range real part of the potential, we use the original Karsh-Mehr-Satz (KMS) form for the free energy. The resulting internal-energy-based model has the following form for the real part of the heavy quark potential where µ = G(ξ, θ)m D [19,22,25], with θ being the angle between the line connecting the quark-antiquark pair and the beam-line direction, m D is the isotropic leading-order Debye mass, a = 0.385, σ = 0.223 GeV 2 [40], and m b = 4.7 GeV is the mass of the bottom quark. The last term taken from Reference [41] accounts for the temperature-and spin-independent finite-quark-mass correction, to obtain our final complex-valued potential (This term is not particularly important for bottomonia states. We include it for historical continuity with our previous works, where we also considered charmonium suppression). For the imaginary part of the potential, we use a small-ξ expansion of the heavy-quark potential [18,20,21] [ and ψ 1 and ψ 2 are defined as follows The algorithm from Reference [23,42] is used to solve the resulting 3d Schrödinger equation on a regular lattice by transforming to imaginary time and using the finite difference time domain (FDTD) method. Using this method, we compute the real and imaginary parts of the binding energy over a range of Λ from 144 MeV to 1037 MeV. For each Λ, we compute the real and imaginary binding energies for a range of anisotropies, ξ, from −0.3 to 200, with an irregular spacing which accounts for the fact that the majority of the time the system probes small values of ξ. We use a N 3 = 256 3 lattice with lattice spacing a = 0.1 GeV −1 ≈ 0.02 fm and a = 0.15 GeV −1 ≈ 0.03 fm for a total box length of L = Na ≈ 5.04 fm and L = Na ≈ 7.56 fm, for Υ and χ b states, respectively. The imaginary-time step size for the algorithm is taken to be ∆τ = a 2 /8.
The real and imaginary binding energies are extracted using [25] where Negative values of [E bind ] only occur for large values of ξ, which is a consequence of the small-ξ expansion. Large values of ξ correspond to a nearly free streaming quark-gluon plasma, so it is expected that the widths of bottomonia states return to vacuum values, which are on the order of ∼ keV, which effectively allows us to set [E bind ] = 0 for this specific case.

Bottomonium Suppression
We use the following rate equation to account for in-medium bottomonia state decay, where all variables depend on longitudinal proper time τ = √ t 2 − z 2 , the transverse coordinate x ⊥ , and the spatial rapidity ς = arctanh(z/t). The rate of decay is computed by [25] where γ dis is a large value which is chosen such that a completely unbound state decays quickly. We emphasize that E bind and hence Γ are local quantities of τ, x ⊥ , and ς through the 3 + 1d background evolution of the transverse temperature Λ(x) and the local momentum-space anisotropy ξ(x).

Survival Probability
The survival probability is determined by integrating Equation (10) over proper time with the integration limits set dynamically. The lower integration limit for the proper-time integration is max(τ form , τ 0 ), where τ 0 is the initial proper time for plasma evolution and τ form is the time-dilated formation time of the state in question, computed via τ form (p T ) = γτ 0 form = E T τ 0 form /M, where M is the mass of the state. The rest-frame formation time τ 0 form for each of the states is taken to be inversely proportional to the vacuum binding energy of each state [43], and thus, τ 0 form = 0.2, 0.4, 0.6, 0.4, 0.6, and 0.6 fm/c, for Υ(1S), Υ(2S), Υ(3S), χ b (1P), χ b (2P), and χ b (3P), respectively. The upper limit for the proper-time integration is determined to be the proper time at which the local energy density becomes less than the energy density of an N c = 3 and N f = 2 ideal gas of quarks and gluons with a temperature of T f = 192 MeV. This is the temperature at which screening effects are assumed to turn off rapidly due to the transition to the hadronic phase.
Transverse momentum cuts were implemented by assuming that the transverse momentum distribution of bottomonia states is proportional to E −4 T , Once the transverse momentum cut is applied, we average R AA over the transverse plane where T A (x, y) = ∞ −∞ dz n A x 2 + y 2 + z 2 is the nuclear thickness function and n AB = T A (x + b/2, y)T B (x − b/2, y) is the overlap density function. In the above relations, it is assumed that n AA sets the probability for bottomonia production at a given point in the transverse plane. After the p T cuts and spatial averaging are performed, centrality averaging is performed by converting the impact parameter, b, to centrality, C, using the Glauber model, and then integrating over the centrality cut with a probability distribution proportional to e −C/20 with 0 < C < 100 [44].

Excited State Feed Down
A certain fraction of Υ(1S) and Υ(2S) states produced in URHICs are formed via the decay of excited states. To compute the post feed down R AA for the Υ(1S) and Υ(2S) states, R AA is computed taking into account the suppression of excited states which undergo late time feed down. For this purpose, we use p T -averaged feed down fractions obtained recently from a compilation of p-p data available from ATLAS, CMS, and LHCb [45]. The fractions f Υ(1s) i and f Υ(2s) i are given in Tables 3 and 4, which specify the fraction of a particular state which contributes to the Υ(1S) and Υ(2S) states. The inclusive feed down R AA of each of the states is calculated using a linear combination of the primordial R AA using the respective feed down fractions R Υ(nS) AA is the primordial suppression for the ith state.  Table 4. Feed down fractions to the Υ(2S) state.

Results
We now turn to our results and predictions. In Figure 1, we show the "primordial" R AA for the six states as a function of N part for the case of 4πη/s = 1. The primordial suppression is the result obtained prior to taking into account feed down effects. The left panel shows the result for √ s NN = 2.76 TeV, and the right panel shows our prediction for √ s NN = 5.023 TeV. From this figure, one can immediately see the sequential suppression of states, with excited states showing more suppression, which is a consequence of the lower dissociation temperatures compared to, for example, the Υ(1S) state. Going from 2.76 TeV to 5.023 TeV, we see increased primordial suppression for all states considered, with the "primordial" R AA decrease by approximately 31%, 48%, 44%, 51%, 45%, and 44% for the Υ(1S), Υ(2S), Υ(3S), χ b (1P), χ b (2P), and χ b (3P) states, respectively, for central collisions. Finally, we note that at both energies, the model predicts significant primordial suppression of the Υ(1S), even though the temperatures probed in the plasma are below the state's disassociation temperature for most of the plasma evolution. The primordial suppression seen is a result of the fact that the decay rate remains large even below the naive disassociation temperature for a given state. In Figures 3 and 4, we compare the inclusive suppression of the Υ(1S) state at 2.76 TeV and 5.023 TeV as a function of p T and y, respectively. The labelling and line types are the same as in Figure 2. Due to the large size of the experimental p T bins, for Figure 3 we have binned our predictions into the experimental bins. From Figures 3 and 4, we find that the model once again predicts enhanced suppression of the Υ(1S) when going from 2.76 to 5.023 TeV as a function of p T and y. In the lowest p T bin, we predict a decrease in R AA of approximately 25% for the case 4πη/s = 1. For |y| = 0, we predict a decrease in R AA of approximately 26% for the same case. We also note that there is a slight increase in suppression for forward rapidities, which is due to the increased plateau halfwidth used in the initial conditions.  In Figures 5-7, we collect our predictions for the Υ(1S) and Υ(2S) inclusive R AA for √ s NN = 5.023 TeV collision energy. The line styles are the same as in the previous figures. As can be seen from these figures, we also predict further suppression of the Υ(2S) at √ s NN = 5.023 TeV as a function of N part , p T , and y. Finally, in Figure 8, we plot the inclusive Υ(1S) suppression as a function of N part for the case 4πη/s = 1, but now varying the initial momentum-space anisotropy ξ 0 of the QGP. The solid (black) line shows the case ξ 0 = 0, which corresponds to a QGP that is perfectly isotropic at τ 0 . The short-dashed (red) line and long-dashed (blue) lines correspond to ξ 0 = 10 and ξ 0 = 50, respectively. The finite values of ξ 0 map to initial pressure anisotropies in the local rest frame of P L /P T = 0.13 and 0.03, respectively. The presence of an initial momentum-space anisotropy is predicted by both weak and strong coupling approaches, with weak coupling approaches predicting larger initial momentum-space anisotropies than the strong coupling approaches [3,4].  As Figure 8 demonstrates, the effect of initial momentum-space anisotropy is to decrease the amount of Υ(1S) suppression. Going from ξ 0 = 0 to ξ 0 = 50, we see an approximately 14% increase in the inclusive Υ(1S) R AA for a central collision. The effect of initial momentum-space anisotropy on the Υ(2S) is qualitatively similar, however, the effect is larger. We find that the inclusive Υ(2S) R AA increases by approximately 27% for ξ 0 in the same range.

Conclusions
The primary goal of this work was to make predictions for bottomonia suppression at the latest LHC URHIC energy of √ s NN = 5.023 TeV. In order to make our predictions, we had to make some extrapolations of the initial conditions for the background. The key changes in going from 2.76 TeV to 5.023 TeV were a 16% increase in the initial central temperature, an increase in the width of the central rapidity plateau [38], and a small change in the nucleon-nucleon scattering cross section [35].
With these assumptions, we then made predictions with the same model employed in our prior works on this subject [24][25][26]. The model includes the effect of in-medium dissociation of bottomonia, uses a full 3 + 1d anisotropic viscous hydrodynamics background, and takes into account non-equilibrium modifications of the heavy quark potential associated with large momentum-space anisotropy. We did not include any cold nuclear matter effects or regeneration, since these are expected to be small for the bottomonia states. Comparing higher energy 5.023 TeV collisions to 2.76 TeV collisions, our model predicts that one should see enhanced suppression of bottomonia states as a function of centrality, transverse momentum, and rapidity. We have provided quantitative calculations of R AA for both the Υ(1s) and Υ(2s), including feed down effects which can be compared with forthcoming data from the LHC.
Finally, in addition to extrapolating our results to the higher energy, we looked at the effect of initial momentum-space anisotropy on bottomonia suppression. For levels of momentum-space anisotropy predicted by theoretical models, we find an approximately 14% increase in the Υ(1s) inclusive R AA . We plan to study this effect in more detail in the future, since it can affect the extracted value of η/s when comparing with experimental data.