Fireball tomography from bottomonia elliptic flow in relativistic heavy-ion collisions

We calculate the elliptic flow of bottomonia produced in Pb$\,+\,$Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. We consider temperature-dependent decay widths for the anisotropic escape of various bottomonium states and observe that the transverse momentum dependence of bottomonia elliptic flow provides a tomographic information about the QGP fireball at different stages of its evolution. For the space-time evolution of the fireball, we employ simulation results from the 3+1D quasiparticle anisotropic hydrodynamic model. We find that our results for transverse momentum dependence of bottomonia elliptic flow are in reasonable agreement with experimental results from the ALICE and CMS collaborations.

We calculate the elliptic flow of bottomonia produced in Pb + Pb collisions at √ sNN = 5.02 TeV. We consider temperature-dependent decay widths for the anisotropic escape of various bottomonium states and observe that the transverse momentum dependence of bottomonia elliptic flow provides a tomographic information about the QGP fireball at different stages of its evolution. For the space-time evolution of the fireball, we employ simulation results from the 3+1D quasiparticle anisotropic hydrodynamic model. We find that our results for transverse momentum dependence of bottomonia elliptic flow are in reasonable agreement with experimental results from the ALICE and CMS collaborations.

I. INTRODUCTION
Relativistic heavy-ion collisions create a strongly interacting hot and dense medium of deconfined quarks and gluons called the quark-gluon plasma (QGP). The thermodynamic and transport properties of QGP are governed by quantum chromodynamics (QCD). Therefore, experimental results from high-energy nucleus-nucleus collisions at the CERN Large Hadron Collider (LHC) and the BNL Relativistic Heavy Ion Collider (RHIC) are important to test the predictions of QCD. There are several experimental observables for studying the QGP, one of which is the study of the in-medium propagation of heavy quarks (charm and bottom quarks) and their quarkonium bound states. Heavy quark observables provide very useful internal probes of the QGP [1]. In particular, heavy quarkonia created in these collisions are significantly affected by the medium, leading to distinctive features in their spectra and yields. Therefore, the study of heavy quarks and quarkonia is important to understanding QGP evolution [1][2][3][4][5][6][7][8].
The study of the suppression of heavy quarkonia, as a signature for the creation of an equilibrated QGP, was first proposed by Matsui and Satz based on the in-medium dissociation of quarkonium bound states due to Debye screening of the quark-antiquark (qq) potential [2]. In this simple classical picture, the quarkonia bound states having the largest binding energy, i.e., the J/ψ and the Υ, possess the highest dissociation temperatures [2,3]. However, in recent years, it has been shown from first-principles finite-temperature QCD calculations that the in-medium qq potential also contains an imaginary part which is generally associated with the inmedium breakup rate of quarkonium states [9][10][11][12][13]. This in-medium break-up rate leads to large thermal widths for quarkonia states, resulting in their dissociation even at temperatures at which they would be expected to survive in the traditional scenario [14][15][16][17][18][19].
Another important aspect in the study of heavy quarkonia is the role of the QGP dynamics. In an expanding QGP, the temperature drop leads to a decrease of the Debye screening of the potential, thereby resulting in recombination of uncorrelated qq pairs into a stable bound state [20][21][22][23]. While this process is quite significant for charmonia yields [24] at RHIC and LHC, recombination seems to be marginal for bottomonia owing to the fact that bb pairs are produced much less abundantly in initial hard scatterings than charmonia [25,26]. Moreover, since dissociation is a non-instantaneous process, the survival probability of bottomonium states is in-medium path length dependent and is naturally described within an "escape mechanism" scenario [27][28][29][30][31].
For the spatially anisotropic QGP produced in a generic high energy heavy-ion collision, this path-length dependence leads to an anisotropic emission pattern of quarkonia, which in turn results in anisotropic momentum distribution and therefore flow-like signatures, as was first predicted for J/ψ [32].
Since bottomonia are very massive, it is expected that they are produced very early in the initial hard scattering and their propagation is not affected significantly by the medium collectivity. Therefore, the anisotropic escape mechanism would constitute the major contribution to the observed momentum anisotropy, which is usually expressed in terms of Fourier harmonics v n . Moreover, since they are produced early with a spatial distribution given by hard binary collisions, their propagation through the medium and escape provide a tomographic probe of the fireball. In this work, we calculate the elliptic flow of bottomonia produced in Pb + Pb collisions at √ s NN = 5.02 TeV at the LHC. We consider temperaturedependent decay widths for anisotropic escape of various bottomonium states and demonstrate that the transverse momentum dependence of bottomonia elliptic flow reflects properties of the fireball at different stages of its evolution.

II. HYDRODYNAMICAL MODEL
For the initial spatial distribution of the bottomonium states, we assume that, due to their large masses, they are produced in initial hard scattering when the heavy nuclei collide. Therefore, we distribute the bottomonium production points in the transverse plane according to the number of binary collisions, N coll (x, y), obtained from the optical Glauber model. Along the longitudinal direction (z), we assume a flat distribution. For the initial transverse momentum (p T ) distribution of the bottomonia, we adopt a power law dependence obtained from numerical PYTHIA simulations for proton-proton (p + p) collisions, scaled by the corresponding mass number of colliding nuclei. 1 The initial Υ distribution in p + p collision from PYTHIA is given by [35] where Y is the longitudinal rapidity in momentum space, ) the maximum value of forward rapidity of the bottomonia. Finally, the momentum rapidity follows a Gaussian distribution given by: In the present calculations, we integrate over Y and consider the resulting p T distribution.
Owing to the quantum mechanical uncertainty principle, each bottomonium bound state has a finite formation time τ form . In the bottomonium rest frame, the value of intrinsic formation time τ 0 form is assumed to be inversely proportional to the binding energy in vacuum of each state. We use τ 0 form = 0.2, 0.4, 0.6, 0.4, 0.6 fm/c for Υ(1S), Υ(2S), Υ(3S), χ b (1P ) and χ b (2P ) states, respectively [17]. Due to time dilation, the formation time in the plasma frame is given by τ form = E T τ 0 form /M , where M is the mass of the bottomonium states and E T = p 2 T + M 2 is the transverse energy with p T being the transverse momentum. After their formation, the color-neutral bb bound states suffer few elastic scatterings in QGP, and due to their large mass, they propagate quasi freely with nearly straight-line trajectories.
In order to accurately evaluate the survival of the bottomonium states propagating through the QGP medium, we need a realistic simulation of the space-time evolution of the expanding QGP. In the present work, 1 We do not try to account for the difference between the parton distribution functions in a Pb nucleus and a proton. The effect of including nuclear PDFs (and their uncertainties) is expected to be small for v 2 , while the reasonable agreement of our calculated R AA with the experimental results validate our approach [33,34]. we employ the simulation results from the recently developed 3+1D quasiparticle anisotropic hydrodynamics (aHydroQP) framework [36]. This framework, which includes both shear and bulk viscosities as well as an infinite number of transport coefficients, has been quite successful in describing multitude of experimental observables for both Pb+Pb collisions at √ s NN = 2.76 TeV [37,38] and Au+Au collisions at √ s NN = 200 GeV collisions [39]. In this version of the aHydroQP code, the shear viscosity to entropy density ratio (η/s) is assumed to be constant and all other transport coefficients are determined self-consistently within a quasiparticle frameowrk. The temperature dependence of the quasiparticle mass is extracted from lattice QCD results for the entropy density [40]. For details of this framework we refer the reader to Ref. [41] which presents a comprehensive review of the approach. The dissipative hydrodynamical evolution of QGP is started at τ i = 0.25 fm/c. In order to construct the initial energy density profile in the transverse plane, we use optical Glauber model with a linear combination of binary collisions and wounded-nucleon. The fraction of binary collisions is set to κ binary = 0.15. The nucleonnucleon inelastic collision cross section is taken to be σ NN = 67 mb. For the density distribution of the incoming nuclei, we use a standard Woods-Saxon profile, with saturation density n 0 = 0.17 fm −3 , nuclear radius R A = (1.12A 1/3 − 0.86A −1/3 ) fm and skin depth d = 0.54 fm. We take the initial central temperature to be T 0 = 630 MeV and the specific shear viscosity to be 4πη/s = 2.0. These values were determined by fits to identified particle spectra and successfully reproduce many bulk observables in 5.02 TeV Pb-Pb collisions [42,43]. The space-time evolution of the effective temperature in each centrality class was obtained by using 4D interpolating functions constructed from the aHydroQP model evolution for a given centrality. For details concerning the hydrodynamical framework and fits to 2.76 TeV experimental data, please see Ref. [42].
Once we know the temperature of the medium along the bottomonium trajectory from the aHydroQP model, one can calculate the breakup rates of the bottomonium states by using temperature-dependent thermal decay widths. We adopt the in-medium dissociation of different bound bb states from recent state-of-the-art numerical solution of 3D Schrödinger equation with a temperaturedependent complex heavy-quark potential [14,19]. This potential is based on a generalization of Karsch-Mehr-Satz potential [3] obtained from the internal energy. It also includes an imaginary part emerging from Landau damping of the exchanged gluons within the framework of hard-thermal-loop approximation [12][13][14]. This imaginary part of the potential results in a temperaturedependent in-medium breakup disassociation rate for each state. The change in sign of the real part of the binding energy of a given state, on the other hand, allows one to determine when the state becomes completely unbound. For completeness, we illustrate in Fig. 1 the temperature dependence of the decay width (Γ(T )) of directly produced Υ(1S) states, as used in our present calculations.
Note that, we set the temperature scale to T c = 160 MeV, below which we assume that the screening effects due to plasma rapidly decreases because of transition to the hadronic phase. This temperature corresponds to the temperature of the assumed QGP phase transition. This is a reasonable approach since, in the hadronic phase, the widths of the bottomonium states become approximately equal to their vacuum widths, which are negligible in comparison to the in-medium widths. Taking all these aspects into consideration, the thermal decay width of formed bottomonia is modeled as Γ(T (x, y, τ )) = 2 for > 0, 10 GeV for ≤ 0, where we use the notations ≡ Re[E bind (T (x, y, τ ))] and ≡ Im[E bind (T (x, y, τ ))] for the real and imaginary parts of the in-medium binding energy, respectively. The temperature depends on the instantaneous transverse position (x, y) of the bottomonium state at time τ along its trajectory. The value of 10 GeV in the second case in Eq. (3) is chosen to ensure that there is rapid dissociation of states which are fully unbound. However, in practice, the results are insensitive to this value provided that it is large enough to quickly dissociate the state under consideration.
In order to account for the finite time for bound states to develop from the primordial bb wave packet, we include a formation-time effect. This effect tends to decrease the suppression rate which one can intuitively associate with a geometric expansion of the bb wave packet from almost point-like production to the bound-state size [25]. Therefore we correct the decay width accordingly as where, τ form is the bottomonium formation time in the plasma rest frame. We further assume that bottomonia only dissociate when the QGP is formed, i.e. Γ eff (T ) = 0 for τ < τ i . In order to study the propagation of bottomonium states, we consider that the QGP is formed at τ i and the bottomonia are created with transverse coordinates (x, y) and transverse momentum p T along the azimuthal φ p direction. The position of this state at any future time τ is then given by 2 where v T ≡ p T /E T is the bottomonium transverse velocity and τ ≡ τ − τ i . The instantaneous temperature along a given bottomonium trajectory are obtained from the aHydroQP model predictions for the full 3+1D evolution of the QGP effective temperature. Subsequently, the effective in-medium decay width is determined for these temperatures using Eqs. (3) and (4). The in-medium decay width, as obtained above, determines the survival probability of a given state as it propagates through the medium. The final transmittance probability for a bottomonium state, labeled by j, is given by [14][15][16][17] where the final time τ f is determined in the aHydroQP model as the proper time at which the local temperature of the expanding medium drops below the freeze-out temperature T f = 130 MeV. Using the above equation, we calculate the transmitted spectra as a function of the azimuthal angle and the transverse momentum of all of the produced bottomonium states, dN j p T dp T dφ p = dx dy N coll (x, y) (7) Using the above equation, one can calculate the nuclear modification factor and elliptic flow of the j-th resonance state as a function of transverse momentum using the standard definitions.

III. RESULTS AND DISCUSSIONS
In order to calculate the nuclear modification factor R AA and elliptic flow v 2 of bottomonia due to escape probability through an anisotropic medium, we numerically evaluate Eq. (7). We consider Pb + Pb collisions at √ s NN = 5.02 TeV in various centrality classes. We obtain the "raw" spectra for each bottomonium state by integrating over the transverse temperature profile in Eq. (7). The feed down contribution from excited states are calculated using a p T -averaged feed down fraction taken from a recent compilation of p + p data at LHC [17]: .1% and f 3P →1S = 1.5%. Subsequently, the inclusive spectrum for Υ(1S) is obtained by linear superposition of the raw spectra for each state and given by dN all The inclusive spectrum thus obtained is then used to calculate the nuclear modification factor and elliptic flow of Υ(1S). In Fig. 2, we show our model predictions for the transverse momentum dependence of nuclear modification factor R AA (p T ) of Υ(1S) in Pb + Pb collisions at √ s NN = 5.02 TeV, together with results from ALICE (0-90% centrality class, forward rapidity 2.5 < |y| < 4) [44] and CMS (0-100% centrality class, mid-rapidity |y| < 2.4) [45]. The model calculations are performed in the same |y| < 2.4 rapidity interval as CMS in ten centrality bins: 0-10%, 10-20%, · · · , 90-100% centrality. The results are averaged with a weight proportional to e −C/20 (C being the centre point of the respective centrality bin) and accounting for the inclusive Υ(1S) yield in each bin, yielding our reported results for the 0-100% centrality interval [17]. Our calculated R AA is in reasonable agreement with the available measurements from both CMS and ALICE. Although our calculations were performed at central rapidity, the weak rapidity dependence of charged particle multiplicity results in similar dynamics at mid- and forward rapidities [46]. We note that, comparing to the previous calculations of R AA of Υ(1S) within the anisotropic hydrodynamics framework [17], these prior results did not include pre-resonance absorption. The effect of the pre-resonance absorption is larger for excited states, which have a longer formation time, and at higher p T , due to time dilation.
In Fig. 3, we show our results for the transverse momentum dependence of elliptic flow parameter for Υ(1S) state, together with results from ALICE [47] and CMS [48]. Both model calculations and CMS results are for the 5-60% centrality class and at mid-rapidity (|y| < 2.4), while the ALICE measurements are for the same centrality interval but at forward rapidity (2.5 < y < 4.0). Our calculations predict a negative v 2 for low p T , which could be explained as follows. The bulk of bottomonia which survive the QGP are those produced near the "edge" of the fireball, see Fig. 4. This stems from the fact that the binary overlap profile produces pairs with elliptic initial geometry depending on the impact parameter. Those with momentum pointing into the QGP are strongly suppressed during their traversal. In the transverse expansion, the eccentricity of the fireball tends to decrease due to faster expansion along xaxis and can even overtake states which were previously "escaped". As a result, the low-momentum bottomonia undergo more suppression along the x-axis as they have to travel larger effective distance in the medium. 3 On the other hand, bottomonia with high transverse momentum quickly traverse the fireball and feel the initial, predominantly elliptic, geometry. This leads to positive elliptic flow at large p T as shown in Fig. 3. In order to explicitly verify our argument of negative v 2 at low p T , in presence of transverse expansion, we calculated the v 2 of directly produced Υ(1S) states by considering only longitudinal expansion of the fireball. In absence of transverse expansion the negative component of v 2 at low p T vanishes and v 2 → 0 as p T → 0, a feature common for all bottomonium states. Therefore, the transversemomentum dependence of elliptic flow can provide valuable information about the time evolution of the fireball's shape in the transverse direction. This is analogous to taking snapshots of fireball anisotropy at different proper time which is sampled differently by states with different p T . With this interpretation, one sees that bottomonium elliptic flow can provide tomographic information about the evolution of the transverse geometry of the fireball.
Before we conclude, let us briefly discuss the sensitivity of the predicted observables to systematic uncertainties associated with our model. The intrinsic formation times of the different bottomonium states are assumed to be inversely proportional to their vacuum binding energies. But the choice is not unique and highly model dependent.
To quantify the resulting uncertainties on R AA (p T ) and v 2 (p T ) of the inclusive Υ(1S) states, we vary the intrinsic formation time of each resonance state by ±50%, from their default values. We find that the resulting change of either observable is small. Generically, decreasing the formation time results in a larger suppression, i.e. a smaller R AA , and a slightly higher v 2 . This affects bottomonia with a large p T , as already observed in our study at √ s NN = 2.76 TeV [31]. Accounting for bottomonia from the recombination of independent b andb will also have a very small impact, since the mechanism is responsible for only a small fraction of the total bottomonium yield. The variation of initial conditions like τ i , spatial profile, viscosity etc. of the hydrodynamic evolution are beyond the scope of the present work. As mentioned earlier, these parameters are adopted from [42,43], where they are determined by fitting the identified hadron spectra and can explain a multitude of bulk observables in 5.02 TeV Pb+Pb collisions. The relative motion between the bottomonia and the expanding QGP has not taken into consideration while calculating the in-medium decay widths. As discussed in [31], the effect on the estimated v 2 of the correction to the dissociation rates due to Doppler shifting [53,54] is expected to be small, especially at small p T . Finally we examine the role of pre-resonance absorption, that remains operative during the spatial expansion of the initially compact bb pair to a full grown resonance state. The choice of Γ(T ) eff is motivated from the linear increase in time of the size of the QQ dipole. To test the influence of the ansatz, we show in Fig. 5 how turning off the dissociation of the pre-resonant bb pairs affects the shapes of R AA and v 2 of the inclusive Υ(1S) states. The suppression of the nascent pairs becomes important for the "fast" (p T > m) bottomonium states, which in the absence of pre-resonance absorption can escape the plasma before they are suppressed.

IV. SUMMARY AND CONCLUSION
In this paper, we calculated the elliptic flow of bottomonia produced in Pb + Pb collisions at √ s NN = 5.02 TeV at the Large Hadron Collider. We employed PYTHIA and the Glauber model to generate initial distribution of momentum and position of produced bottomonia. For the space-time evolution of the fireball, we employed a recently developed 3+1D quasiparticle anisotropic hydrodynamic model which has been tuned to reproduce experimental observables such has identified hadron spectra and elliptic flow at √ s NN = 5.02 TeV.
We considered the effect of temperature-dependent decay widths on the traversal of various bottomonium states through a hot QGP medium with an spatially anisotropic geometry. We have also accounted for the feed down contribution from higher excited states to the bottomonium ground state. We propose that the transverse momentum dependence of bottomonia elliptic flow provides tomographic probe of the QGP fireball transverse expansion. We found that our results for transverse-momentum dependence of bottomonia elliptic flow are in good agreement with experimental results from the ALICE and CMS collaborations. They also agree with recent calculations [55,56] simulating the real time quantum evolution of bottomonium states in the QGP, relaxing the assumption of an adiabatic evolution of the bb system [57]. Looking forward, one can anticipate greatly increased statis-tics in future runs from the experimental collaborations at LHC. Hopefully, with these increased statistics, we can more quantitatively assess the efficacy of our model.