Bottomonium Continuous Production from Unequilibrium Bottom Quarks in Ultrarelativistic Heavy Ion Collisions

We employ the Langevin equation and Wigner function to describe the bottom equark dynamical evolutions and their formation into a bound state in the expanding Quark Gluon Plasma (QGP). The additional suppressions from parton inelastic scatterings are supplemented in the regenerated bottomonium. Hot medium modifications on $\Upsilon(1S)$ properties are studied consistently by taking the bottomonium potential to be the color-screened potential from Lattice results, which affects both $\Upsilon(1S)$ regeneration and dissociation rates. Finally, we calculated the $\Upsilon(1S)$ nuclear modification factor $R_{AA}^{\rm rege}$ from bottom quark combination with different diffusion coefficients in Langevin equation, representing different thermalization of bottom quarks. In the central Pb-Pb collisions (b=0) at $\sqrt{s_{NN}}=5.02$ TeV, we find a non-negligible $\Upsilon(1S)$ regeneration, and it is small in the minimum bias centrality. The connections between bottomonium regeneration and bottom quark energy loss in the heavy ion collisions are also discussed.


I. INTRODUCTION
Since charmonium was proposed as a probe for the existence of the deconfined matter called "Quark-Gluon Plasma" (QGP) [1], its production mechanisms has been widely studied based on coalescence model [2][3][4] and transport models [5][6][7] in nucleus-nucleus collisions. The nuclear modification factor R AA is a measurement of cold and hot medium suppressions on quarkonium yields. The cold nuclear matter effects include the nuclear absorption [8], Cronin and shadowing effects [9][10][11]. The first one is negligible at LHC colliding energies due to strong Lorentz dilation, where "spectator" nucleons already move out of the colliding region before the formation of a quarkonium eigenstate. Cronin effect will shift the momentum distribution of primordially produced hidden-and open-charm(or bottom) states [12,13]. This can be included by a proper modification of their transverse momentum distributions in pp collisions [14,15]. Shadowing effect is weak at RHIC colliding energies, but important at the LHC colliding energies. All the cold nuclear matter effects can be included in the heavy flavor initial distributions prior to the hot medium effects. In nucleus-nucleus collisions at √ s N N = 2.76 TeV and 5.02 TeV, regeneration from charm and anti-charm quarks is widely believed to dominate the prompt charmonium yields [16,17]. This is supported by the enhancement of J/ψ R AA and the suppression of J/ψ mean transverse momentum square p 2 T J/ψ observed at 2.76 TeV and 5.02 TeV: the regenerated J/ψs from thermalized charm quarks carry small momentum compared with the primordially produced ones, this will pull down the p 2 T J/ψ of final prompt J/ψ in nucleus-nucleus collisions [18].
However for bottomonium, the situation seems not so clear. Transport model calculations suggest a nonnegligible bottomonium regeneration in √ s N N = 5.02 TeV Pb-Pb collisions [19]. Also, experimental data hinted a stronger bottomonium regeneration at 5.02 TeV compared with 2.76 TeV, but within its large uncertainty which prevents solid conclusions [20,21]. Considering that heavy quark mass is very large, it takes some time to reach kinetic thermalization [22,23] in the fast cooling QGP with an initial temperature ∼ 500 MeV in AA collisions. Non-thermalization of bottom quark momentum distribution will suppress the combination probability of b andb quarks in QGP. However, the ratio of hiddento open-bottom states is at the order of 0.1% which is smaller than the charm flavor. This may make the yield from (b +b → Υ(1S) + g, b +b + ζ → Υ(1S) + ζ with ζ = g, q,q) not negligible compared with primordially produced Υ(1S). It is very interesting to develop a realistic model to consider the dynamical evolutions of bottom quarks, bottomonium regeneration process and the following hot medium suppression after they are regenerated. We employ the "Langevin equation + Wigner function + (gluon, quasi-free) dissociations" for bottom quark and bottomonium evolutions in heavy ion collisions. Furthermore, we consider the hot medium modifications on bottomonium properties at finite temperature by taking the color-screened heavy quark potential (extracted from Lattice free energy F (r, T ) [24]). With color-screened heavy quark potential in time independent Schrödinger equation, we obtain the mean radius and binding energy of Υ(1S) at different temperatures. Each of them will be used in Wigner function (regeneration rate) and (gluon, quasi-free) dissociation rates.
Our paper is organized as follows. In Section II, we introduce the Langevin equation and Wigner function for heavy quark dynamical evolutions and their combination. The hot medium modifications on Υ(1S) properties (such as the mean radius and binding energy) are also studied based on potential model. In Section III, we introduce the hydrodynamic equations for QGP expansion in nucleus-nucleus collisions. The relevant inputs of heavy quarks are presented in Section IV. In Section V, we give the Υ(1S) nuclear modification factor R rege AA from the combination of b andb quarks in QGP. Different coupling strength between heavy quarks and QGP (controlling the heavy quark thermalization) are studied in the Υ(1S) regeneration. We summarize the work in Section VI.

A. Langevin Equation for Heavy Quark Diffusion
The heavy quark diffusion in Quark Gluon Plasma can be treated as a Brownian motion, which is widely studied with Langevin equation [25,26]. During the evolution of heavy quarks, they can combine into a quarkonium [2,3,27,28]. The probability of the combination of heavy quark Q andQ depends on their distribuitons in phase space and also the properties of the produced quarkonium at finite temperature [29,30]. Instead of dealing with heavy quark distributions, we employ the Langevin equation to simulate Q andQ evolutions in the hot medium. Combination process (Q +Q → (QQ) bound +g, Q+Q+ζ → (QQ) bound +ζ with ζ = g, q,q) can be included through the Wigner function W QQ→Υ [3]. The Langevin equation is written as where η D (p) and ξ are the drag force and the noise of the hot medium on heavy quarks. ξ satisfies the correlation relation κ is the diffusion coefficient of heavy quarks in momentum space, which is connected with spatial diffusion coefficient D by κ = 2T 2 /D. The drag force in Langevin equation can then be determined by the fluctuationdissipation relation [31] η D (p) = κ 2T E T is the temperature of the bulk medium and E = m 2 Q + | p| 2 is the heavy quark energy. In order to numerically solve the Langevin equation for heavy quark diffusions in QGP, it is discretized as below As long as the time step ∆t for numerical evolutions is small enough, one can assume free motions for heavy quarks with a constant velocity v Q = p/E during this time step, and update the heavy quark momentum by Eq.(4) at the end of each ∆t due to hot medium effects. The medium-induced radiative energy loss [32,33] and parton elastic collisions [34] of heavy quarks can be included in the terms of the drag force η D (p) and the noise ξ. ξ i (t) i=1,2,3 in Eq.(6) is sampled randomly based on a Gaussian function with the width κ/∆t. The initial transverse momentum distribution as an input of Eq.(4) is obtained from PYTHIA simulations. As the initial energy density of QGP changes with coordinates, the heavy quark initial distribution in QGP affects their evolutions, which will finally affect the heavy quark thermalization degree and quarkonioum regeneration. Heavy quark pairs are produced from parton hard scatterings, their density (within rapidity region ∆y) is proportional to the number of binary collisions, where σ QQ pp (∆y) is the heavy quark production cross section in proton-proton collisions within rapidity ∆y.
is the thickness function of lead. Nucleus density ρ Pb ( x T , z) is taken to be the Woods-Saxon distribution. When colliding energy is at the order of ∼TeV, theoretical and experimental studies indicate a strong nucleus (anti-)shadowing effect on heavy quark (quarkonium) production [11]. Furthermore, this effect depends on the nucleon density, which gives different modification on heavy quark production at different positions of the nucleus. We employ a theoretical model (EPS09 NLO) to obtain a shadowing factor r S ( x T , p T , y). The initial distribution of heavy quarks (quarkonium) in Pb-Pb collisions is then taken

B. Heavy Quark Recombination Process with Wigner Function
In the nucleus collisions, Heavy quark pairs are produced and evolve inside the QGP as a Brownian motion. During QGP evolutions, Q andQ can meet each other and combine into a bound state, which may survive from the hot medium due to its large binding energy. If Q and Q can reach kinetic equilibrium, their relative momentum ( p Q − pQ) and relative distance ( X Q − XQ) will be small, which enhances the combination probability of Q andQ [35]. The discussion about connections between heavy quark thermalization and the quarkonium regeneration is left to the next secton. Wigner function is widely used for hadron production in the coalescence model. It gives the probability of Q andQ combining into a bound state with relative distance r = X Q − XQ, and momentum q = p Q − pQ, see the function below [3] f A 0 is the normalization factor. Here we neglect the contributions of momentum carried by light partons in the heavy quark formation formation process, and use the "bb → Υ" to represent both b +b → Υ(1S) + g and b +b + ζ → Υ(1S) + ζ. In realistic simulations, Only 0.01% ∼ 0.1% of bottom quarks can form a bound state in the expanding QGP with a lifetime of ∼ 10 f m/c and an initial temperature of ∼ 500 MeV [36]. Most of them become open heavy-flavor hadrons with a light (anti-)quark. Considering the large ratio of σ bb pp /σ Υ(1S) pp ∼ 1000 (much larger than the charm flavor ∼ 200), this combination process may be important for the bottomonium nuclear modification factor R AA . Bottomonium properties such as binding energy and shape of the wave function can be modified by the hot medium. This can affect the probability of heavy quark combination. We include hot medium effects on bottomonium properties by introducing the temperature dependence of Gaussian width σ(T ). It is connected with the bottomonium mean radius square at finite temperature by σ(T ) = 8 r 2 Υ (T )/3 [3].
In order to study the bottomonium regeneration in QGP, we put one b andb in QGP and evolve each of them with Langevin equation to obtain the probability W bb→Υ of their combination to regenerate a Υ(1S). The total yield of regenerated bottomonium within rapidity ∆y in nucleus-nucleus collisions is scaled by the number of bottom pairs as The average number of bb pairs produced in central rapidity region is only ∼ 1, the regenerated bottomonium yield is proportional to the number of bb pairs. The initial momentum and position of heavy quarks are randomly generated based on the distributions from PYTHIA simulations and Eq.(7) (both with modifications of the shadowing effect). At each time step, one can obtain their relative distance r and relative momentum q, and their combination probability P (r, q) = r 2 q 2 f (r, q). If the probability is larger than a random number between 0 and 1, then the formation process of bottomonium happens. Otherwise, they continue the evolutions with Langevin equation independently untill moving out of QGP (hadronization as open bottom hadrons). After Υ(1S) is regenerated inside QGP, it will decay due to parton inelastic scatterings and color screening from QGP. We supplement this part by the rate equation is the decay rate from gluon and quasi-free dissociations. It is connected with Wigner function (see Eq. (8)).
After one Υ(1S) is regenerated at the time t 0 , the initial condition of Eq.(10) becomes N Υ (t 0 ) = 1, and N Υ (t) decreases with time based on the rate equations, Note that N Υ (t ≥ t 0 ) ≤ 1 where t 0 is the time of Υ(1S) regeneration. P Υ ≈ p b + pb and R Υ are the momentum and coordinate of the center of the regenerated Υ(1S). From hydrodynamic equations, the local temperature of QGP is different at different coordinate R and time. Therefore, we update Υ(1S) position at each time step, and take a new local temperature T of QGP at R Υ (t 1 + ∆t) to recalculate the decay rate Γ diss for Υ(1S) evolution at the next time step. We continue Eq.(10-11) from the time t 0 of Υ(1S) regeneration untill it moves out of QGP (where local temperature is smaller than the critical temperature T c ). Doing sufficient events, N events bb , of putting one b andb in the expanding QGP, and sum Υ(1S) final yields to be N rege Υ (tot), one can obtain the Υ(1S) regeneration probability from b andb evolutions in QGP, W bb→Υ = N rege Υ (tot)/N events bb .

C. Bottomonium Dissociation and Regeneration
Rates at Finite Temperature QGP color screening reduces the binding energy of bottomonium, and increases its decay rate at finite temperature. With heavy quark potential to be V bb = F or V bb = U [37], One can obtain the Υ(1S) binding energy from solving time-independent radial Schrödinger equation (with = c = 1) Here m µ = m b /2 is the the reduced mass in the center of Υ mass frame. E Υ (T ) and ψ(r) are the binding energy and radial wave function of a Υ eigenstate. The mean radius and binding energy (see Fig.1) obtained consistently from Eq.(12) will be used in the Wigner function for bottomonium regeneration and the parton dissociations, respectively. The decay rates from gluon dissociation and quasi-free dissociation [38,39] with temperature dependent binding energy are plotted in Fig.2. When the Υ(1S) is strongly bound, gluon dissociation dominates the Υ(1S) decay rate, such as at the temperature region of T < 0.2 GeV with V bb = U (see red dashed and solid lines). However at a high temperature like T = 300 MeV, strong color screening effect reduces the Υ(1S) binding energy to be around 0.25 GeV with V=U and 0.04 GeV with V=F, which are far below the vacuum value ∼ 1.1 GeV [40]. This makes quasi-free dissociation dominates the decay rate. We include both contributions on the hot medium suppression of regenerated Υ(1S) at the entire temperature region T > T c .

III. HYDRODYNAMIC MODEL
We employ the (2+1) dimensional ideal hydrodynamics to simulate strong expansion of finite sized QGP produced in ultrarelativistic heavy ion collisions, Here T µν = (e + p)u µ u ν − g µν p is the energy-momentum tensor, and (e, p, u µ ) are the energy density, pressure and four-velocity of fluid cells. The equation of state of the deconfined medium is taken as an ideal gas of massless (u, d) quarks, 150 MeV massed s quarks and gluons [41]. Hadron phase is an ideal gas of all known hadrons and   [44,45]. Therefore, we still take the same value of τ 0 = 0.6 fm/c at 5.02 TeV Pb-Pb collisions due to its weak dependence on colliding energy. The transverse expansion of QGP controlled by Eq.(13) starts from τ 0 .

IV. INPUTS OF BOTTOM FLAVOR
The bottomonium regeneration requires the number of bottom quarks in nucleus-nucleus collisions. We determine this by using the production cross section in pp collisions and binary collision scaling in Pb-Pb collisions, N bb PbPb = σ bb pp N coll (b). Lack of experimental data about σ bb pp at 5.02 TeV pp collisions, we extract its value by the linear interpolation between the cross sections at central rapidity at 1.96 TeV and 2.76 TeV collisions. At 1.96 TeV pp collisions, CDF collaboration published the cross section of b-hadrons integrated over all transverse momenta in the rapidity |y| < 0.6 to be 17.6 ± 0.4(stat) +2.5 −2.3 (syst) µb [46]. With this we extract the central value of the differential cross section to be dσ bb pp /dy = 14.7 µb. Combined with the dσ bb pp /dy = 23.28 ± 2.70(stat) +8.92 −8.70 (syst) µb in the central rapidity at 2.76 TeV [47] , we obtain the differential cross section dσ bb pp /dy = 47.5 µb in the central rapidity at 5.02 TeV pp collisions. Our purpose is to study the contribution of regeneration component in the experimentally measured inclusive Υ(1S) yields, presented as the nuclear modification factor [15],  [50].
The coupling strength between bottom quarks and QGP is indicated by the drag coefficient in Langevin equation. The spatial diffusion coefficient is taken to be D(2πT ) = C [51]. Different values of C will be employed to study the effects of bottom quark thermalization on bottomonium regeneration.
In the previous work [36], we employ the Langevin equation plus Wigner function to calculate J/ψ and ψ(2S) regeneration from dynamical evolutions of (anti-)charm quarks, with an assumption that they are produced at each certain temperature, without the following suppression from hot medium after their regeneration. In this work, we improve our approach of "Langevin equation +Wigner function" by considering hot medium modifications on quarkonium properties, which change quarkonium regeneration and dissociation rates through the mean radius r Υ (T ) and binding energy E Υ (T ) of quarkonium.
In the realistic simulations, we generate one b andb randomly in the coordinate and momentum space based on the probability distributions given in previous sections. Then we evolve them separately with two indi- But it is only ∼ 4% at minimum bias centrality (with N p ≈ 200). Note that transport model calculations gave the regeneration R AA to be ∼ 8% in semi-central (N p ≈ 200) and most cenral Pb-Pb collisions in the rapidity |y| < 2.4 at 2.76 TeV [19,50].
One way to testify the contribution of Υ(1S) regeneration in heavy ion collisions is to study the rapidity dependence of the p T −integrated R Υ(1S) AA (y). The bottom quark differential cross section decreases with rapid- The heavy quark potential is taken as internal energy V=U (used to determine the mean radius and binding energy of Υ(1S)). Lines with triangle, square and circle markers are for diffusion coefficient D(2πT ) = 2, 4, 7 respectively. Note that these D values satisfy the relation of 1 < ∼ D(2πT ) < ∼ 7 from pQCD and Lattice calculations [52][53][54].
ity, which will suppress the Υ(1S) regeneration at forward rapidity. For charmonium, the decreasing tendency of R J/ψ AA (y) with rapidity is very strong and explained well by the regeneration mechanism [15]. As charmonium regeneration mainly dominates at the low p T bins and drops to zero at middle and high p T bins. Therefore, R J/ψ AA shows strong decreasing tendency with rapidity at p T > 0 (where regeneration dominates) and almost no rapidity dependence at p T > ∼ 4 GeV/c. Considering the fraction of Υ(1S) regeneration is only 10 ∼ 15% in its inclusive R AA at the impact parameter b=0, we expect a weaker rapidity dependence of R Υ(1S) AA (y). In the minimum bias centrality, This tendency should be much weaker and R
We also did the calculations in Υ(1S) weak binding scenario (V=F), see Fig.4. With weak binding (V=F), Υ(1S) can only be regenerated at T < 400 MeV where its binding energy is non-zero. Also, the dissociation rate of regenerated Υ(1S) is larger for V=F compared with V=U, which makes regenerated Υ(1S) easier to be dissociated. Υ(1S) regeneration with V=F is smaller (see Fig.4). This is also consistent with transport model calculations. The Υ(1S) R AA from regeneration contribution is around 6% (with D(2πT ) = 4) and 3% (with D(2πT ) = 7) in the most central collisions. In all centralities, its contribution is smaller than the value of V = U , see Fig.4.

VI. SUMMARY
We employ the Langevin equation to describe the diffusion of bottom quarks, and Wigner function for the b andb quarks to regenerate Υ(1S)s during the QGP evolutions. After the regeneration of a Υ(1S), it also suffers the color screening and parton inelastic scatterings from QGP. The dissociation and regeneration rates of Υ(1S) are connected with each other by the temperature dependent binding energy E Υ (T ) and mean radius r Υ (T ), which can be obtained simultaneously from Schrödinger equation with the color screened heavy quark potential extracted from Lattice calculations. Supplement the cold nuclear matter effects, we give the full calculations of Υ(1S) regeneration in Pb-Pb collisions at √ s N N = 5.02 TeV. With different drag coefficient in Langevin equation, the energy loss of b andb quarks is different. Strong coupling between bottom quarks and QGP can make b andb lose more energy, and increases the probability of their formation into a bound state. In the scenario of V=U (V=F), we obtain the Υ(1S) regeneration R rege AA to be 0.1 ∼ 0.2 (0.06 ∼ 0.1) in the most central collisions, but negligible in minimum bias centrality. With realistic evolutions of bottom quarks and their hadronization process, we study the regeneration contribution to the Υ(1S) nuclear modification factor R AA measured in experiments, and also build the connection between bottom quark energy loss and bottomonium regeneration in this work.