Coupled Boltzmann Transport Equations of Heavy Quarks and Quarkonia in Quark-Gluon Plasma
Xiaojun Yao, Weiyao Ke, Yingru Xu, Steffen A. Bass, Berndt Müller
PPrepared for submission to JHEP
MIT-CTP/5192
Coupled Boltzmann Transport Equations of HeavyQuarks and Quarkonia in Quark-Gluon Plasma
Xiaojun Yao, a,b, Weiyao Ke, c,b
Yingru Xu, b Steffen A. Bass, b Berndt M¨uller b a Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA, 02139,USA b Department of Physics, Duke University, Durham, NC 27708, USA c Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
E-mail: [email protected] , [email protected] , [email protected] , [email protected] , [email protected] Abstract:
We develop a framework of coupled transport equations for open heavy flavorand quarkonium states, in order to describe their transport inside the quark-gluon plasma.Our framework is capable of studying simultaneously both open and hidden heavy flavorobservables in heavy-ion collision experiments and can account for both, uncorrelated andcorrelated recombination. Our recombination implementation depends on real-time openheavy quark and antiquark distributions. We carry out consistency tests to show how theinterplay among open heavy flavor transport, quarkonium dissociation and recombinationdrives the system to equilibrium. We then apply our framework to study bottomoniumproduction in heavy-ion collisions. We include Υ(1 S ), Υ(2 S ), Υ(3 S ), χ b (1 P ) and χ b (2 P )in the framework and take feed-down contributions during the hadronic gas stage intoaccount. Cold nuclear matter effects are included by using nuclear parton distributionfunctions for the initial primordial heavy flavor production. A calibrated 2 + 1 dimensionalviscous hydrodynamics is used to describe the bulk QCD medium. We calculate both thenuclear modification factor R AA of all Υ( nS ) states and the azimuthal angular anisotropycoefficient v of the Υ(1 S ) state and find that our results agree reasonably with experi-mental measurements. Our calculations indicate that correlated cross-talk recombinationis an important production mechanism of bottomonium in current heavy-ion experiments. Corresponding author. a r X i v : . [ h e p - ph ] A p r ontents A.1 g + H ↔ Q + ¯ Q q + H ↔ q + Q + ¯ Q and g + H ↔ g + Q + ¯ Q |(cid:104) Ψ p rel | r | ψ nl (cid:105)| B Details on Feed-down Contributions 25
Heavy quarkonia are bound states of heavy quark-antiquark pairs Q ¯ Q . The mass spectraof the ground and lower excited quarkonium states can be reasonably well described bythe nonrelativistic Schr¨odinger equation with a potential model [1]. Inside a hot nuclearmedium, i.e., the quark-gluon plasma (QGP), the attractive potential can be significantlysuppressed due to the static screening effect in the plasma [2]. As a result, a bound Q ¯ Q can “melt” at sufficiently high temperature [3–5]. Therefore, quarkonium suppression(compared to a non QGP baseline) can be used as a signal of the formation of a QGP in– 1 –eavy-ion collisions. In general, shallower bound states melt at lower temperatures andone would expect a “sequential” suppression pattern.However, this simple picture is complicated by other in-medium processes such as thedissociation of quarkonium states via dynamical scattering (i.e. the dynamical screeningeffect, which is related to the imaginary part of the Q ¯ Q potential [6, 7]) and heavy quark(re)combination [8, 9]. To account for all these effects, semi-classical transport equationshave been widely applied [10–26]. These calculations typically contain three components:The first part is a temperature-dependent potential that is parametrized from lattice cal-culations of the free energy of a Q ¯ Q singlet [27, 28] or direct lattice calculations of the realpart of the Q ¯ Q potential [29]. The second input is the dissociation rate of each quarkoniumstate. Perturbative calculations of dissociation rates include both the gluo-dissociation pro-cess [30, 31] and the inelastic scattering with the medium (Landau damping). An effectivefield theory of QCD, potential nonrelativistic QCD (pNRQCD) [32–34] has been appliedto study both the static screening of the potential and the dissociation rate [35–37]. Sim-ilar effective theory has also been used to study dark matter bound state formation [38].Some studies also included viscous effects [39–41] and contributions due to the quarkoniummoving with respect to a static medium [42, 43]. The final component is a recombinationmodel, for example, a statistical hadronization model [9], a coalescence model based onWigner functions [19] or a model that is based on detailed balance and a finite relaxationrate [21].Unlike the first two effects, recombination contains significant model-dependencies inmost studies. Here, we shall distinguish between two kinds of recombination: uncorrelatedand correlated recombination. In uncorrelated recombination, the Q and ¯ Q originate fromdifferential initial hard vertices. In proton-proton collisions, these heavy quarks and anti-quarks would almost never (re)combine to form a quarkonium state due to their separationin phase space. However, in heavy-ion collisions, multiple Q ¯ Q pairs are produced from theinitial hard scatterings in one collision event. The momenta of these heavy quarks and anti-quarks change continuously inside the QGP due to diffusion and energy loss. Therefore, thechance for an uncorrelated pair to come close to each other in phase space is higher. Whenthey are close, they may combine to form a quarkonium state. The uncorrelated recombi-nation rate rises with the number of open heavy quarks produced in the collision, which iscrucial to explain the small amount of suppression observed for the J/ψ with rising collisionenergy. Naively, one would expect
J/ψ to be more suppressed at higher collision energiesdue to the hotter medium and stronger plasma screening effects. Therefore, uncorrelatedrecombination is important for the phenomenology of charmonium. However, its effect onbottomonium production is expected to be negligible, since only a few bottom-antibottomquark pairs are produced in one collision.In correlated recombination, the Q and ¯ Q originate from the same initial hard ver-tex. For example, a Q ¯ Q pair in a pre-quarkonium resonance or emerging from a previousdissociation of a quarkonium state is considered a correlated pair. The phenomenologicaleffect of correlated recombination for both charmonium and bottomonium has not beensystematically explored yet.Most recombination models depend on the uncorrelated open heavy quark distribu-– 2 –ions, they can thus only account for uncorrelated recombination. To explore the physicaleffects of correlated recombination, we develop a set of coupled transport equations for Q and ¯ Q ’s as well as for quarkonia, which allows us to study both, uncorrelated and corre-lated recombination. The transport equation of quarkonium provides consistent treatmentof dissociation and recombination, in the sense that both of them are derived from QCDunder systematic nonrelativistic expansions [44, 45]. The derivation is based on the com-bination of the open quantum system framework (in which correlated recombination istaken into account) and the pNRQCD effective field theory. The application of the openquantum system framework to studying quarkonium in-medium dynamics [46–58] and itscombination with pNRQCD [59, 60] has recently drawn significant theoretical interest. Aquarkonium transport coefficient has been defined in Ref. [61].The coupled transport equations allow us to study the in-medium transport of bothopen and hidden heavy flavor states. In this paper, we will use this framework to studybottomonium production in heavy-ion collisions. We will include Υ(1 S ), Υ(2 S ), Υ(3 S ), χ b (1 P ) and χ b (2 P ) states in the transport network. By solving the coupled transportequations via test particle Monte Carlo simulations, we will explore the importance ofcorrelated recombination in bottomonium phenomenology. This paper is organized asfollows: in Sect. 2, we will introduce the set of coupled transport equations. Then inSect. 3, some simulation tests inside a QGP box will be shown and compared to the systemproperties in thermal equilibrium. Details on applying the transport equations to the studyof heavy-ion collisions will be explained in Sect. 4. Results on the nuclear modificationfactor ( R AA ) and the azimuthal angular anisotropy coefficient v of bottomonia will bediscussed later in Sect. 5. Finally, we will draw conclusions in Sect. 6. The set of coupled Boltzmann transport equations for the distribution functions of unboundheavy quark-antiquark pairs Q ¯ Q and each quarkonium state with the quantum number nls ( n is for the radial excitation, l the orbital angular momentum and s the spin) is given by( ∂∂t + ˙ x Q · ∇ x Q + ˙ x ¯ Q · ∇ x ¯ Q ) f Q ¯ Q ( x Q , p Q , x ¯ Q , p ¯ Q , t ) = C Q ¯ Q − C + Q ¯ Q + C − Q ¯ Q (2.1)( ∂∂t + ˙ x · ∇ x ) f nls ( x , p , t ) = C + nls − C − nls , (2.2)where ˙ x = ∂ x ∂t . The left-hand sides of these equations describe the free streaming ofdistribution functions in phase space while the right-hand sides contain collision terms ofthe heavy particles interacting with the plasma. The collision terms with ± superscriptsrepresent the quarkonium dissociation ( − ) and recombination (+) while the term withoutany superscript describes the energy and momentum changes of the open heavy quark-antiquark pairs. In the following, we will explain these collision terms in detail. If we neglect the interaction between the heavy quark-antiquark pair, we can write C Q ¯ Q = C Q + C ¯ Q , (2.3)– 3 –.e., the heavy quark and antiquark interact independently with the medium. In potentialmodels, the interaction between the Q ¯ Q pair is attractive for the color singlet and repulsivefor the color octet. For a Coulomb potential, the color-averaged potential vanishes ifwe assume the color is in thermal equilibrium. Furthermore, in-medium Q ¯ Q potentialsare significantly suppressed, so we expect the potential interaction between the Q ¯ Q pairto be weak. Therefore, throughout this paper, we will neglect the interaction betweenthe Q ¯ Q pair in the open heavy flavor transport equations. As a remark, we note thata factorized Q ¯ Q distribution f Q ¯ Q ( x Q , p Q , x ¯ Q , p ¯ Q , t ) = f Q ( x Q , p Q , t ) f ¯ Q ( x ¯ Q , p ¯ Q , t ) indeedleads to Eq. (2.3). But the opposite is not true in general. Most recombination modelsimplicitly assume the factorization of the Q ¯ Q distributions and thus cannot study thecorrelated recombination. But here we only assume Eq. (2.3) in this work.We will use a weak coupling picture for the transport of open heavy quarks [62–65].The interaction of open heavy quarks (and antiquarks) with the medium is describedby scattering between open heavy quarks and medium partons (which include both light(anti)quarks and gluons, abbreviated as q and g respectively). The collision term C Q includes three types of scattering processes: the elastic 2 → q + Q → q + Q and g + Q → g + Q , the inelastic 2 → q + Q → q + Q + g and g + Q → g + Q + g and the inelastic 3 → q + Q + g → q + Q and g + Q + g → g + Q , and similarlyfor C ¯ Q . In this work, we will use the Monte Carlo simulations in the Lido package [66]to solve the transport equations of open heavy quark-antiquark pairs. The Lido packagecontains both a linearized Boltzmann transport description and a model that is based onthe Langevin equation with radiation corrections. The latter description has been reportedin Ref. [67]. We will only use the linearized Boltzmann description in this work. For the dissociation and recombination terms in the transport equations, we will use theexpressions in Ref. [68]. Detailed expressions of the relevant collision terms can be foundin Appendix A. The calculations therein are based on a version of pNRQCD under thehierarchy of scales M (cid:29) M v (cid:29)
M v (cid:38) T (cid:38) m D . Here M is the heavy quark mass, v thetypical relative velocity between the Q ¯ Q pair inside the bound state, T the temperatureof the QGP and m D the Debye mass. The typical size of quarkonium is roughly give by r ∼ Mv and the typical binding energy is about M v . The last inequality T (cid:38) m D meansthe QGP is weakly-coupled. For charmonium, we have v ∼ . v ∼ . M v ∼
500 MeV.One may worry that the hierarchy is not always true in real heavy-ion collisions. Inthe early time of the QGP expansion, the temperature can be (cid:38)
450 MeV. Due to thestatic screening effect, the in-medium binding energies of quarkonium states can be muchsmaller than their vacuum binding energies, especially for excited states. Thus, our assumedhierarchy indeed breaks down in the early stage and in principle one has to use a differentversion of pNRQCD if there still exists a hierarchy of scales. However, even before thebreakdown of the assumed hierarchy, as the temperature increases, the dissociation ratesof excited quarkonium states such as Υ(2S) and χ b (1P) blow up rapidly. What happensin our calculations is that after a very short time period in the early stage, the excited– 4 –uarkonium states have dissociated and evolve then as an unbound, correlated Q ¯ Q pair.This is the right physics: When our hierarchy of scales breaks down, we can either have T (cid:38) M v (cid:29)
M v or M v (cid:29) T (cid:29) M v . We do not consider M v (cid:29) T (cid:29) M v herebecause the real values of v do not allow this hierarchy to happen for both charmoniumand bottomonium. (For (cid:29) to be valid, one at least needs a factor of three in the ratio.) Sowe only need to consider T (cid:38) M v (cid:29)
M v here. Whenever this occurs, one expects both theplasma screening effects to be extremely strong ( Mv gives the rough size of the quarkoniumstate) and no bound states can be well-defined. In this case, quarkonium, even if it binds,binds very weakly and it is reasonable to assume that it behaves more like an unbound Q ¯ Q pair. Then the correct description is the transport of an unbound Q ¯ Q pair. Therefore, inour calculations, whenever the temperature is high enough to break our assumed hierarchy,the quarkonium state dissociates after a tiny time step and then evolves as an unbound,correlated Q ¯ Q pair. This is a good approximation of the correct physics. The transportof excited quarkonia as well-defined bound states is valid again in the later stage of theevolution, when the temperature drops and our hierarchy is resumed. In this sense, mostexcited quarkonium states are probably generated via recombination in the later stage ofthe QGP evolution. Their suppression mechanism is mainly the decorrelation of the Q ¯ Q pair in coordinate and momentum space (or decoherence of their wavefunction). Thoseexcited states that are observed may probably come from recombination of Q ¯ Q pairs thatare still correlated. We will discuss this in more detail in Sect. 5.For the current calculations, we work to the leading order in the nonrelativistic expan-sion (the multipole expansion under this hierarchy of scales is equivalent to a nonrelativisticexpansion). At this order, the higher Fock state | Q ¯ Qg (cid:105) of quarkonium (in which the Q ¯ Q is acolor octet) is suppressed at least by two powers of v with respect to the leading Fock state | Q ¯ Q (cid:105) (in which the Q ¯ Q is a color singlet) [69]. Therefore, in our calculations, quarkoniumis always a Q ¯ Q pair in the color singlet. Unbound Q ¯ Q pairs can be in the color singletor octet states. At the order of the nonrelativistic expansion that we are studying rightnow, the transition between a quarkonium and an unbound pair only occurs via a dipoleinteraction between the color singlet and octet states. In other words, the dissociationproduct is a Q ¯ Q in the color octet state and only color octet Q ¯ Q pairs can recombine asquarkonia via the dipole interaction. The dipole term works here because of the hierar-chy M v (cid:29) T . When the quarkonium or the Q ¯ Q pair is at rest with respect to the localmedium, the typical energy of a medium parton is E g ∼ T . The typical size of quarkoniumis given by r ∼ Mv . So the quarkonium size is small in the sense that rE g (cid:28) Q ¯ Q is moving with respect to the local medium (for example, when the quarkonium state hasa finite transverse momentum), the typical energy of medium partons in the quarkonium(or Q ¯ Q ) rest frame is boosted, E g ∼ γT where γ is the boost factor that depends on therelative velocity between the quarkonium (or the center-of-mass of the Q ¯ Q ) and the localmedium. The condition rE g (cid:28) V s = − C F α pot s r for color singletand V o = N c α pot s r where C F = N c − N c and N c = 3. We will take α pot s in the potentials tobe a parameter and choose its value as α pot s = 0 .
36 throughout this paper. The couplingconstant in the scattering vertices will be taken to be constant α s = 0 .
3. The effects ofrunning coupling and nonperturbative potentials will be left to future studies . Since theoctet potential is non-zero here, the wavefunction of the unbound octet pair is a Coulombscattering wave rather than a plane wave. In this way, we re-sum an infinite number ofCoulomb exchanges between the octet pair in the initial (for recombination) or final (fordissociation) state. At the leading order in the nonrelativistic expansion, the potential isindependent of the orbital angular momentum and spin. Thus the dissociation rates arethe same for states separated by fine and hyperfine splittings. The recombination rates arealso the same up to a spin degeneracy factor g s = or .We include the following scattering channels for the dissociation and recombination: g + H ↔ Q + ¯ Q , q + H ↔ q + Q + ¯ Q and g + H ↔ g + Q + ¯ Q where H can indicate anyquarkonium state. The first process is induced by real gluon absorption (for dissociation)and emission (for recombination). The last two processes are mediated by virtual gluons(inelastic scattering). The elastic scattering between quarkonia and medium gluons isneglected here because it occurs at higher orders in the multipole expansion (or equivalentlyhere, nonrelativistic expansion) [68]. The direct transitions between different quarkoniaspecies are also omitted due to the same reason. Expressions of the reaction rates of the1 S , 2 S and 1 P states can be found in Appendix A. For 3 S and 2 P states, we assume theycannot exist inside the QGP. In other words, they dissociate immediately after enteringthe QGP and cannot be (re)generated inside the hot medium.One important feature of our framework is the inclusion of correlated recombination(because we keep track of the evolution of the two-particle distribution of the Q ¯ Q pairrather than the distribution of a singlet heavy quark). The correlated recombination leadsto cross-talk between different quarkonium states. When an excited quarkonium state suchas Υ(2 P ) or Υ(3 S ) dissociates, the produced Q ¯ Q can form a lower excited state or theground state Υ(1 S ). As will be discussed later, this is an important production mechanismfor the ground state. When an excited state dissociates at high temperature, the dissociated Q ¯ Q pair may form the ground state and then survive the subsequent evolution. When thetemperature drops and a ground state dissociates occasionally (the dissociation rate is stillnon-vanishing even if the static screening effect is small), the dissociated Q ¯ Q pair may forman excited state. The time evolution of the whole system is a network of reactions amongunbound heavy quarks, antiquarks and all quarkonia states. Our framework can handlethis reaction network and study the physical impacts of the cross-talk recombination. The different values of α pot s and α s chosen here already hint at the importance of nonperturbativepotentials. – 6 – .3 Monte Carlo Simulations We solve the coupled transport equations by test particle Monte Carlo simulations. Wesample a certain number of Q ¯ Q pairs and quarkonia according to the distribution functionsat the initial time. Mathematically, the distribution function for each particle species isrepresented by f ( x , p , t ) = (cid:88) i δ ( x − x i ) δ ( p − p i ) , (2.4)and similarly for the two-particle distribution functions. Here x i and p i are the positionand momentum of the i -th sampled test particle. The integral of a distribution functionover the whole phase space gives the total number of the particle associated with thatdistribution. The positions and momenta of the sampled particles obey their initial distri-bution functions. Then we evolve the positions and momenta of all particle species stepby step. The time step size is chosen as ∆ t = 0 .
01 fm/c in the laboratory frame. At eachtime step, we consider each of the following processes.
The position of the particle changes according to x ( t + ∆ t ) = x ( t ) + ∆ t p ( t ) E ( t ) , (2.5)where E ( t ) and p ( t ) are the energy and momentum of the particle at the current time step. This process is implemented via the Lido package. For given heavy quark (antiquark)momentum and local temperature, the package calculates the scattering rate of the heavyquark (antiquark) with the medium in each scattering channel. If a certain scatteringprocess occurs, the package generates a light scattering partner utilizing local mediumproperties (temperature and flow field). Its outgoing momentum is sampled from thedifferential scattering rate. The outgoing momentum of the heavy quark (antiquark) canbe obtained from energy and momentum conservation. Then we can update the momentumof the heavy quark (antiquark).
For a quarkonium state nls with a certain momentum (velocity) and a position at somelocal temperature, we calculate its dissociation rate in the laboratory frame. The methodto obtain the dissociation rate from the collision term C − can be found in Ref. [68]. The ratetimes the time step size leads to the dissociation probability in this step. If it is determined(by Monte Carlo sampling) that the quarkonium state dissociates, we sample the incomingand/or outgoing momenta of the relevant light particles and obtain the momenta of theoutgoing Q ¯ Q pair from energy and momentum conservation. The positions of the unbound Q and ¯ Q are given by the position of the quarkonium before dissociation. Then we removethis quarkonium state from the relevant particle list and add the produced Q ¯ Q pair to thelist of heavy quarks and antiquarks. – 7 – .3.4 Recombination For each unbound Q ¯ Q pair, we need to calculate their recombination rate. In the MonteCarlo simulation, the two-particle distribution function of the unbound Q ¯ Q pair is repre-sented by f ( x Q , p Q , x ¯ Q , p ¯ Q , t ) = (cid:88) i,j δ ( x Q − x i ) δ ( p Q − p i ) δ ( x ¯ Q − ˜ x j ) δ ( p ¯ Q − ˜ p j ) . (2.6)The recombination rate of a specific pair can be obtained by the following replacement δ ( x Q − x i ) δ ( p Q − p i ) δ ( x ¯ Q − ˜ x j ) δ ( p ¯ Q − ˜ p j ) → δ (cid:16) x cm − x i + ˜ x j (cid:17) δ (cid:0) p cm − ( p i + ˜ p j ) (cid:1) δ (cid:16) p rel − p i − ˜ p j (cid:17) e − ( x i − ˜ x j )22 σ , (2.7)where x cm , p cm , x rel and p rel are the center-of-mass (cm) and relative positions andmomenta. This Gaussian ansatz is motivated by the recombination formula derived inRef. [44], in which the recombination term for a Q ¯ Q far away from each other is suppressedexponentially by the bound state wavefunction. The width of the Gaussian is chosen to bethe typical size of the bound state. More specifically, for 1 S , σ = a B where a B = α s C F M is the Bohr radius. For 2 S and 1 P state, we have σ = 2 a B and similarly for higher excitedstates. The recombination rate of this specific pair can be obtained by doing the followingreplacement in the recombination term C + in the Boltzmann equation (The expression of C + can be found in Appendix A.) f ( x Q , p Q , x ¯ Q , p ¯ Q , t ) → δ (cid:16) p rel − p i − ˜ p j (cid:17) e − ( x i − ˜ x j )22 σ . (2.8)In practice, for each unbound Q ¯ Q pair with positions x i , ˜ x j and momenta p i , ˜ p j , wefirst boost the pair into their cm frame. Then we compute their relative momentum andcalculate their recombination rate in the cm frame and then boost the rate back into thelaboratory. If the recombination into a specific quarkonium state occurs, we sample themomentum of the outgoing quarkonium state based on the differential rate and energy-momentum conservation. The position of the quarkonium state is given by the cm positionof the unbound pair before recombination. Finally we remove the unbound pair from thelist of heavy quarks and antiquarks and add a new quarkonium state to the relevant list ofquarkonium states. Before we study quarkonium production in heavy-ion collisions, we validate our test particleMonte Carlo simulations for the coupled transport equations. We solve the equations insidea cubic volume of QGP matter at fixed temperature with a side length L = 10 fm. Thebox has periodic boundary conditions, i.e., when a particle reaches the boundary of thebox, it appears on the opposite side. In other words, the medium behaves as an infinitelylarge QGP with a finite heavy flavor density (which include both open and hidden heavyflavor states). – 8 –e focus on the bottom system, since the nonrelativistic expansion works better forbottom than for charm. We will sample a fixed number of unbound b ¯ b pairs initially. Theirpositions are randomly distributed in the volume while their momenta obey thermal oruniform distributions (other distributions can also be specified). In the mode of uniformmomentum distribution, each component of the bottom quark’s momentum, p i , i = x, y, z ,is sampled from a uniform distribution between 0 and 3 GeV. Using the uniform momentumdistribution mode, if we turn off the transport equation of unbound pairs and only simulatethe transport equations of quarkonia, we find that the system does not properly thermalize.It is the transport of open heavy flavors that drives the thermalization of all heavy quarkstates. These findings have been reported in Refs. [24, 70]. Here we extend the previousstudies to the case of excited states. For simplicity, we will simulate all the cases with theopen heavy flavor transport equations turned on. The lessons we learn by comparing thecase with open heavy flavor transport equations turned on and that off have been discussedbefore. We focus on demonstrating the consistency of the numerical implementation ofdissociation and recombination here.We initialize N b, tot = 50 bottom quarks and N ¯ b, tot = 50 antibottom quarks in theQGP volume. Their positions are random and their momenta are sampled in the uniformdistribution mode, as described above. We consider the following three cases:1. The temperature of the QGP box is fixed to be 300 MeV throughout. We only simu-late the Υ(1 S ) channel for quarkonium, i.e., only the dissociation and recombinationof Υ(1 S ) are allowed.2. Only the Υ(2 S ) channel is turned on. Here, the temperature is fixed at 180 MeV.
3. Same as case 2, but only χ b (1 P ) is studied.We will compute the hidden bottom fraction as a function of time, which is defined by N b, hidden N b, tot = N b, hidden N b, open + N b, hidden = N H , (3.1)where N H is the number of all bottomonium states. For the three cases listed above,only one bottomonium state contributes to N H . We will compare simulation results ofthe hidden bottom fraction with that in thermal equilibrium, which can be calculatedas follows. In thermal equilibrium, the numbers of open bottom quarks and a specificquarkonium state H are given by N eq i = g i Vol (cid:90) d p (2 π ) λ i e − E i ( p ) /T , (3.2)with i = b , ¯ b or H , E i ( p ) = (cid:113) M i + p relativistically and M i + p M i nonrelativistically.Here g i is the spin and color degeneracy factor and λ i is the fugacity. We have g b = 6 In principle, since we use a Coulomb potential here, excited bottomonia states exist at high temperature.In the simulation tests, one can choose a higher temperature. But at high temperature, the hidden bottomfraction from excited bottomonia is tiny in thermal equilibrium and one requires large statistics to obtainreasonable results. Furthermore, reaction rates become bigger as temperature increases, so a smaller timestep is required. – 9 –or the bottom quark, g H = 3 for Υ( nS ) and averaged χ b ( nP ). Since the total number ofbottom quark is equal to that of antibottom, λ b = λ ¯ b . Using detailed balance, we have λ H = λ b λ ¯ b = λ b . With this relation, one can solve the fugacities from the balance equationwith a given volume: N b, tot = N eq b + N eq H . (3.3)Once we obtain the fugacities, we can compute the hidden bottom fraction in thermalequilibrium: N eq H N b, tot . (3.4)The comparison between our simulation results and the system properties in thermal equi-librium is depicted in Fig. 1 for the three cases listed above. The results are obtained fromaveraging 10000 simulation events. As can be seen from the plots, the interplay betweendissociation and recombination can drive the system to equilibrium. The relaxation rate ofthe system is about 7 − p T spectrum and density of the heavy particles, as well as the temperatureof the medium. Our simulations can reproduce the correct limit of the hidden bottomfraction in thermal equilibrium. These tests serve as consistency checks in our studies andwe are now ready to move on and simulate real collision systems. In this section we will describe our numerical setup to study bottomonium production inheavy-ion collisions. We include Υ(1 S ), Υ(2 S ), Υ(3 S ), χ b (1 P ) and χ b (2 P ) states in thequarkonium transport equations. To solve the coupled transport equations, we need aninitial condition and a bulk QGP medium. For the calculations of R AA , we also need toinclude all feed-down contributions from excited states to the ground and lower excitedstates in the hadronic gas stage. We will explain the calculation of the initial phase spacedistribution, the bulk QGP evolution and the feed-down network in a sequence below. The initial phase space distributions are determined as follows. The momenta of heavyquark-antiquark pairs and each quarkonium state are calculated and sampled utilizing
Pythia [71]. The quarkonium production calculation in
Pythia is based on the NRQCDfactorization [69]. We use the nuclear parton distribution function (nPDF) parametrizedby EPPS16 [72]. The nPDF is the only cold nuclear matter (CNM) effect we include. Thesuppression factors caused solely by the CNM effects R CNM for the Pb-Pb collisions at5 .
02 TeV are shown in Fig. 2 for different bottomonia states. In the plots, we omit theuncertainty bands for simplicity. We use the same bins in the transverse momentum p T andrapidity y as the CMS measurements, which will be shown later. The production of χ b (2 P )is not available yet in Pythia , so we will assume the CNM effect on χ b (2 P ) is the same– 10 – t (fm/c) − − − N b , h i dd e n / N b , t o t simulation w/ HQ energy lossnon-relativistic equilibriumrelativistic equilibrium (a) Case 1: Υ(1 S ) in a 300 MeV QGP box. t (fm/c) − − − N b , h i dd e n / N b , t o t simulation w/ HQ energy lossnon-relativistic equilibriumrelativistic equilibrium (b) Case 2: Υ(2 S ) in a 180 MeV QGP box. t (fm/c) − − N b , h i dd e n / N b , t o t simulation w/ HQ energy lossnon-relativistic equilibriumrelativistic equilibrium (c) Case 3: Υ(1 P ) in a 180 MeV QGP box. Figure 1 : Comparison of the hidden bottom fraction between simulations and properties in thermalequilibrium. as that on χ b (1 P ). Furthermore, Pythia with nPDF EPPS16 cannot describe the nuclearmodification factor R pAu = 0 .
82 measured by the STAR collaboration [73]. So we use
Pythia to calculate the p T -dependence of the CNM effect but use R CNM = ( R pAu ) = 0 . y -dependence of R CNM is not needed here because the STAR measurements are carried out in the mid-rapidityregion ( | y | < . s ( τ , x ⊥ ) ∝ (cid:18) ( T A ) p + ( T B ) p (cid:19) /p , (4.1)where T A = T A ( x ⊥ ) and T B = T B ( b ⊥ − x ⊥ ) are the nuclear thickness functions of thetwo projectiles separated by the impact parameter b ⊥ , p is a parameter that has beencalibrated on bulk observables of the QGP and τ is the thermalization time of the systemafter the initial collision. We choose τ = 0 . p T (GeV) . . . . . . R C N M (a) CNM effect as a function of transversemomentum. y . . . . . . . . . R C N M (b) CNM effect as a function of rapidity. Figure 2 : CNM effects on bottomonia at 5 .
02 TeV Pb-Pb collisions as functions of transverse momentumand rapidity. streaming. The parametrized entropy density will be used as the initial condition of thehydrodynamic equation, which will be explained in the next subsection. The productionof heavy quark-antiquark pairs and quarkonia are thought to be hard processes becauseof the large mass scale. Therefore, their initial production probability is assumed to beproportional to the binary collision density in the transverse plane, which is proportionalto T A T B . Since the nuclear thickness function is used in both the initial entropy densityand the initial heavy quark-antiquark production, the corona effect is taken into accountin our calculations.All heavy particles are assumed to be produced at τ = 0 and they propagate via freestreaming until τ = τ = 0 . τ = 0 is considered a valid assumptionfor heavy quark-antiquark pairs because their production time is about M in their restframe. But this may no longer be true for quarkonia, whose formation time is estimatedas Mv in their rest frame, which is the time for the heavy quark-antiquark pair to developthe quarkonium wavefunction of the relative motion . For bottom quarks and bottomonia, M ∼ .
05 fm/c while Mv ∼ . τ . Since we assume particles are free streaming before τ , itdoes not matter whether the ground state is free streaming as a fully formed quarkoniumstate or an unbound pre-resonant heavy quark-antiquark pair. But it matters for excitedquarkonium states, because they may be formed inside the thermal QGP. In other words, at τ , we may only have pre-resonant excited quarkonium states. But this is just a tiny effect inour calculations. The excited states have very large dissociation rates when the temperatureis high. They will dissociate in one time step after entering the QGP and becoming For heavy quark-antiquark pairs and quarkonia with finite transverse momenta, their formation timesin the laboratory frame will be boosted by a γ -factor. So even the ground quarkonium state may forminside the thermal medium. But it should be noted that the total dissociation rate in the laboratory frameincreases with the transverse momentum. The effect of the finite formation time is small. See the followingmain text for the arguments. – 12 –orrelated unbound pair to evolve further. It really does not matter if we have pre-resonantor completely formed excited quarkonium states, because they evolve as unbound Q ¯ Q pairwhen entering the thermal QGP. Improvements can be done by including the relativemomentum broadening of the Q ¯ Q pair in the pre-thermalization stage, which will be leftto future studies.For our calculation on bottomonia, we will only initialize quarkonia states but notunbound b ¯ b pairs, because the number of such pairs produced in current heavy-ion collisionexperiments is very small (for central Pb-Pb collisions at 5 .
02 TeV, the average numberof unbound b ¯ b pairs is less than one per rapidity). So unlike charmonium production,uncorrelated recombination is negligible for bottomonium production. Since most of theunbound heavy quark-antiquark pairs are produced back-to-back in the transverse planeinitially, correlated recombination of these unbound pairs is also negligible. Due to theeven smaller production cross section of bottomonium, we will assume in our calculationsat most one bottomonium state is produced initially in one heavy-ion collision event. We use a 2 + 1 dimensional viscous hydrodynamic model VISHNew [75, 76] to describe theevolution of bulk QGP matter. The package numerically solves ∂ µ T µν = 0 (4.2)with the energy-momentum tensor T µν = eu µ u ν − ( P + Π)( g µν − u µ u ν ) + π µν , (4.3)Π = − ζ ∇ · u, (4.4) π µν = 2 η ∇ (cid:104) µ u ν (cid:105) (4.5)for given initial conditions. Here e and P are the local energy density and pressure, g µν isthe metric and u µ is the local four-velocity of the medium. Π is the bulk stress with thebulk viscosity ζ , and π µν is the shear stress tensor with the shear viscosity η . The anglebracket means traceless symmetrization.Both the bulk and shear viscosities are parameters here and can be temperature-dependent. We will use the parametrizations and calibrations in Ref. [77]. Ref. [77] usesthe TRENTo model to calculate the initial entropy density s and then obtain the energydensity and pressure with an equation of state calculated by lattice QCD. To obtain theinitial stress-energy tensor at τ , Ref. [77] further assumes the initial flow velocity and theviscous terms vanish at τ . All parameters in the TRENTo model and the hydrodynamicequations are calibrated to experimental observables on light hadrons. For 2 .
76 TeV Pb-Pbcollisions, we will use the parameters calibrated to charged particles shown in Table IIIof Ref. [77] with the exception of the parameter p in the TRENTo model. We choose p = 0, which is consistent with the calibration in Ref. [77] (the calibrated value of p is0 . +0 . − . ). For the parameter calibrations in other collision energies and systems, we willfollow Ref. [67]. – 13 – .3 Feed-down Network We terminate the transport evolution when the local medium temperature drops to T c =154 MeV, i.e., we neglect the dissociation of quarkonium states in the hadronic gas. Thefinal nuclear modification factor R AA ( i ) for i = Υ(1 S ), Υ(2 S ), Υ(3 S ), χ b (1 P ) and χ b (2 P )includes both the CNM effect R CNM and the hot medium effect. The evaluation of theCNM effect has been discussed in Sect. 4.1. We will now discuss how we compute the hotmedium effect.Our calculation includes, unlike many others, correlated recombination and we findthat this effect plays a crucial role. We need to take into account the following situation:Some quarkonium state produced initially may end up as a different quarkonium state.For example, a 2 P state produced initially, melts inside the QGP and recombines as a 1 S state which survives the following in-medium evolution. To this end, for each centralitybin, we simulate N init i events in each of which one i -quarkonium state is initialized. Afterthe in-medium evolution, among these N init i events, there are N i → j events that have a j -quarkonium state in the end. We note in general (cid:80) j N i → j ≤ N init i . For example, N S → S isthe number of Υ(2 S ) states generated from those N init1 S events where the initial quarkoniumstate is a Υ(1 S ). N i → i is the “surviving” number whose physical meaning is well-known (italso has contribution from the correlated recombination), and has been calculated in manystudies. N i → j with j (cid:54) = i is the contribution from correlated, cross-talk recombination.As explained in Sect. 4.1, we only initialize one bottomonium state in each relevantcollision event. For simplicity of the expressions, we will assume all N init i ’s are the sameand equal to N init = 30000 in our calculations of R AA . The total final number of the i -quarkonium state produced from N init ’s collision events ( N init events for each state weconsider) is given by (we include CNM effect here) N final i = (cid:88) j R CNM ( j ) σ j σ i N j → i , (4.6)where R CNM ( j ) is the CNM effect on the primordial production of the j -quarkonium stateand σ i is the initial primordial production cross section of the i -quarkonium state (withoutany feed-down contributions). The ratios of the initial production cross sections are listedin Appendix B. Now we include the feed-down contributions for Υ(1 S ) and Υ(2 S ): N init , fd2 S = N init + (cid:88) j =3 S, P σ j σ S N init Br[ j → S ] (4.7) N init , fd1 S = N init + σ S σ S N init , fd2 S Br[2 S → S ] + (cid:88) j =1 P, S, P σ j σ S N init Br[ j → S ] (4.8) N final , fd2 S = N final2 S + (cid:88) j =3 S, P N final j Br[ j → S ] (4.9) N final , fd1 S = N final1 S + N final , fd2 S Br[2 S → S ] + (cid:88) j =1 P, S, P N final j Br[ j → S ] , (4.10)where we have used N init i = N init for all i ’s. The branching ratios in vacuum are listed in– 14 –ppendix B. Finally the nuclear modification factors for i = 1 S , 2 S are given by R AA ( i ) = N final , fd i N init , fd i . (4.11)For j = 3 S , 1 P and 2 P , we do not consider any feed-down contributions. Then theirnuclear modification factors are given by R AA ( j ) = N final j N init j = N final j N init . (4.12) We first show the results of R AA as a function of centrality at 5 .
02 Pb-Pb collisions. Thecomparison between our calculation and the measurements by the CMS collaboration isshown in Fig. 3. To explicitly demonstrate the importance of correlated recombination, weshow the full calculation results in Fig. 3a and the results without any contribution from thecross-talk recombination in Fig. 3b, i.e., in the latter case, the contribution from an initial i -quarkonium state ending up as a final j -quarkonium state ( j (cid:54) = i ) is excluded. Remainingsame-species ( i → i ) correlated recombination is still included in the latter case, thoughwe expect its contribution to be much smaller than the cross-talk recombination. As canbe seen from the comparison of the two figures, the cross-talk recombination is crucial todescribe the data. For the 1 S state production, most excited states dissociate quickly whenentering the QGP because of the initial high temperature. These melted correlated b ¯ b pairsmay form 1 S state after the dissociation because the 1 S state can exist at high temperature.For the 2 S production, most of the primordially produced 2 S states cannot survive thein-medium evolution in most centrality bins. They are regenerated when the temperaturecools down and allows their existence in the medium. At low temperature, a 2 S statecan be formed from a correlated unbound b ¯ b pair (which may come from a dissociatedquarkonium state such as a 1 S state). For the 3 S production, since we assume they cannotexist inside QGP (melting temperature is below 154 MeV), correlated recombination doesnot affect its production. The 3 S state is only produced in very peripheral collisions dueto the corona effect.Next we discuss the results of R AA as a function of the transverse momentum, shownin Fig. 4a. First we notice that the experimental results of R AA are almost flat as afunction of p T . This is a highly non-trivial result because the CNM effect has a dramaticdependence on p T , as shown in Fig. 2. This means for bottomonium, the p T -dependentCNM effect is washed out by the hot medium effect, which leads to a final flat R AA . Ourcalculations can describe the R AA for p T <
10 GeV and start to deviate from the dataas p T becomes bigger. This may indicate the breakdown of the nonrelativistic expansionthat our calculations are heavily replied on. For a finite p T , the typical energy of mediumlight particles in the rest frame of quarkonium is given by T √ − v T rather than T , where v T is the transverse velocity of the quarkonium with respect to the local medium. Ourcalculations are valid if rT √ − v T (cid:28) r ∼ Mv is the typical size of quarkonium.– 15 –
100 200 300 400 N part . . . . . . R AA (a) With cross-talk recombination. N part . . . . . . R AA (b) Without cross-talk recombination. Figure 3 : Bottomonia R AA as functions of centrality at 5 .
02 TeV Pb-Pb collisions. Experimental dataare taken from Ref. [78]. p T (GeV) . . . . . . R AA (a) R AA ( p T ). . . . . y . . . . . . R AA (b) R AA ( y ). Figure 4 : Bottomonia R AA as functions of transverse momentum and rapidity at 5 .
02 TeV Pb-Pbcollisions. Experimental data are taken from Ref. [78].
Under our assumed hierarchy rT (cid:28)
1, the condition rT √ − v T (cid:28) v T is small.Higher order contributions in the nonrelativistic expansion become important and haveto be included consistently as p T increases . Our calculations probably indicate that thishappens when p T >
10 GeV. One should note that for charmonium this can happen at amuch smaller p T because what matters in the argument is the transverse velocity ratherthan the transverse momentum. Since most bottomonia are produced at low transversemomentum, the deviation at mid and high p T does not affect the centrality and rapiditydependence of R AA in our calculations. Under our assumed hierarchy, the multipole expansion is equivalent to the nonrelativistic expansion.One should keep in mind that the dipole interaction vertex grows linearly with the quarkonium size andthus the reaction rate grows with the size. But physically we know there is an upper limit of the rate whenthe size is large: The rate cannot be bigger than twice the reaction rate of one heavy quark interacting withthe medium. When the quarkonium size is large enough, the heavy quark and antiquark interact with themedium almost independently. Therefore, we expect that including higher order terms in the nonrelativisticexpansion will reduce the rate and thus the R AA at large p T will go up. – 16 –hen we show the results of R AA as a function of the rapidity in Fig. 4b. Our calcula-tions are consistent with the experimental measurements. Since we use a 2 + 1 dimensionalviscous hydrodynamics, the hot medium effect is independent of y . Furthermore, we seein Fig. 2 that the y dependence of the CNM effect is also mild. Therefore, the final R AA is almost flat in the rapidity. N part . . . . R AA (a) R AA ( N part ) at 2 .
76 TeV Pb-Pb collisions. N part . . . . . . . R AA N coll uncertainty 2S+3S, theory2S+3S, syst2S+3S, stat (b) R AA ( N part ) at 200 GeV Au-Au collisions. p T (GeV) . . . . . . R AA (c) R AA ( p T ) at 2 .
76 TeV Pb-Pb collisions. p T (GeV) . . . . . . R AA (d) R AA ( p T ) at 200 GeV Au-Au collisions. . . . . y . . . . . . R AA (e) R AA ( y ) at 2 .
76 TeV Pb-Pb collisions.
Figure 5 : Bottomonia R AA as functions of centrality, transverse momentum and rapidity at 2 .
76 TeVPb-Pb and 200 GeV Au-Au collisions. Data of the CMS and STAR measurements are taken fromRefs. [79] and [73] respectively. The p T dependence of R AA for the 200 GeV Au-Au collisions is calculatedand compared for the centrality range 0 − – 17 –he comparison between our calculations and experimental measurements for 2 . R AA (2 S + 3 S ) in peripheral Au-Aucollisions at 200 GeV. In Fig. 5b, our calculations of R AA (2 S + 3 S ) are consistent withthe measurements at central collisions, though are slightly lower than the central values.But our calculation result has a large discrepancy with the data point at the peripheralcollision. The uncertainty of the experimental measurement there is quite large due tothe large uncertainty in the determination of N coll . The uncertainty associated with N coll is small for central collisions. This discrepancy with the single data point in peripheralcollisions leads to the discrepancy in the p T dependence of R AA (2 S + 3 S ), as depicted inFig. 5d. We expect that the future sPHENIX collaboration will provide data with highprecision and then the origin of the current discrepancy issue will be more clear: Either thecalculations are consistent with improved measurements or the calculations need improving.Improvements of the calculations can be carried out in understanding the CNM effect at 200GeV Au-Au collisions. Since we assume nPDF is the only CNM effect, we use a constant R CNM = 0 .
67 for all centrality bins. It is possible that the real CNM effect is centrality-dependent and R CNM goes up as the collision becomes more peripheral. An increase of R CNM in peripheral collisions will lessen the discrepancy we see here. p T (GeV) − . . . . . v theory, 10 − , | y | < . − , . < | y | < − , | y | < . Figure 6 : Azimuthal angular anisotropy coefficient v of Υ(1 S ) at 5 .
02 TeV Pb-Pb collisions. Theexperimental results of the ALICE and CMS collaborations are taken from Refs. [80] and [81].
Finally, we study the azimuthal angular anisotropy of Υ(1 S ). In particular, we computethe v coefficient which is defined by E d N S d p = 12 π d N S p T d p T d y (cid:0) v cos[2( φ − Ψ RP )] (cid:1) , (5.1)in which we neglect higher order harmonics. Here φ − Ψ RP is the azimuthal angle of the 1 S state with respect to the reaction plane, which is defined event-by-event. Our calculation– 18 –esult for the 5 .
02 Pb-Pb collision in the centrality range −
90% is shown in Fig. 6.We also plot the experimental results measured by the ALICE collaboration which arein a different rapidity and centrality regions, just to show the state-of-art of the currentbottomonium v measurements. We stop our calculations at p T = 24 GeV because as wehave seen in the R AA comparison, higher order corrections in the nonrelativistic expansionstart to become important as p T increases, which are neglected in our current setup. We cancalculate v at higher p T values in our current setup but their physical meaning is less robustand we cannot learn much from doing that. Our calculation result is consistent with thecurrent experimental data, though current measurements have large statistical errorbars.The last CMS data point has a quite large p T range: 10 −
50 GeV. Our nonrelativisticexpansion calculation would definitely break down at 50 GeV. Furthermore, at such a high p T , the fragmentation production would start to dominate and the suppression mechanismwould mainly be jet energy loss. Due to the steep p T spectrum in the primordial productionof bottomonium, most of the contribution to the last data point comes from p T ∼ − v . Thefirst contribution is the path dependence. In peripheral collisions, the QGP has a ellipticshape. Quarkonia moving along the longer axis will be more suppressed. Also, the reactionrates of quarkonium in the medium depend on the relative velocity between the quarkoniumand the local medium, which has a flow velocity. This also influences the shape of v asa function of p T . Finally, after the dissociation of quarkonium, the unbound Q ¯ Q pair candevelop some v by interacting with the medium. Later if they recombine, the v willbe partly or fully inherited by the regenerated quarkonium. Uncorrelated recombinationcan also contribute to v and is crucially important for charmonium production in heavy-ion collision. Open charm quarks that develop v during the in-medium evolution willcontribute to the charmonium v if they recombine. But this v generation mechanism isnegligible for bottomonium since the number of unbound b ¯ b pairs produced in one collisionevent is smaller than one per rapidity in mid-rapidity. Future precise measurements onthe azimuthal angular anisotropy will greatly help us understand the in-medium dynamicsof quarkonium, especially how quarkonia with finite transverse momenta interact with anexpanding medium. These non-equilibrium transport properties of quarkonium are noteasy to study via finite temperature lattice QCD calculations. The interplay betweentheory and phenomenology will help deepen our understanding on these, in particular onceexperimental data with high precision are available. In this paper, we develop a set of coupled transport equations for open heavy quark-antiquark pairs and quarkonia. Our framework is capable of calculating observables forboth open and hidden heavy flavors. Dissociation and recombination rates are calculated For the v observable, the CMS collaboration defines the centrality using hard probe triggers. So inthe v calculation, we define the centrality bin in TRENTo by using the number of binary collisions N coll rather than the multiplicity. – 19 –rom pNRQCD via a systematic weak-coupling and nonrelativistic expansion, with the as-sumption of a weakly-coupled QGP. Recombination rates depend on the real-time heavyquark-antiquark distribution function, which is solved by the transport equation for openheavy flavors. Our framework can handle both uncorrelated and correlated recombina-tion. By numerically solving the coupled transport equations via Monte Carlo techniques,we demonstrate how the ground and excited bottomonia states approach thermal equilib-rium inside a QGP box with a constant temperature. Furthermore, we study bottomoniaproduction in heavy-ion collisions using our framework. The initial phase space distribu-tions are calculated from Pythia with the nPDF EPPS16 and the TRENTo model. Themedium description utilizes a 2 + 1 dimensional viscous hydrodynamic model. The pa-rameters of the TRENTo model and the hydrodynamic equations have been calibrated toexperimental observables on light hadrons. We include Υ(1 S ), Υ(2 S ), Υ(3 S ), χ b (1 P ) and χ b (2 P ) states in our reaction network. Our calculations demonstrate the importance ofcorrelated cross-talk recombination in bottomonium phenomenology and can describe mostof the experimental data. Discrepancies are seen at mid and large transverse momenta,where omitted higher order terms in the nonrelativistic expansion become gradually moreimportant.The current formalism can be improved in several ways. First, one can include theeffect of the running coupling constant and higher order corrections in the nonrelativisticexpansion. One may also consider extending the calculation of the reaction rates to thecase of a strongly-coupled QGP, which will be more relevant at low temperature. Second,the simple Coulomb potential can be replaced by a more realistic nonperturbative potentialmodel. The nonperturbative potential may be parametrized and the parameters can becalibrated to experimental observables such as R AA and v , for example, by using therecently developed Bayesian analysis techniques [82]. A simultaneous description of bothopen and hidden heavy flavor observables is also possible in our framework and is worthexploring. With a nonperturbative potential, the technical challenge would be to develop afast numerical algorithm to sample the momenta of the outgoing Q ¯ Q pair (in dissociation)or quarkonium (in recombination) from the differential reaction rates. No previous studieshave done this, but it is crucial for a consistent microscopic treatment like our calculationhere. With a Coulomb potential, reaction rates have analytic expressions which are easierto handle in inverse transform and importance samplings (For the case of how to sampleefficiently with a Coulomb potential, see Chapter 4.2 of Ref. [83]). For a nonperturbativeparametrization of the screened potential, one has to include its dependence on the relativevelocity between the quarkonium state and the hydro-cell. We use Coulomb potentialhere so we do not need to consider the velocity dependence in the potential. Also, theCNM effect at 200 Au-Au collisions should be better understood and constrained by themeasurements in p-Au collisions. Furthermore, the transport in the pre-thermalizationstage should be investigated. For open heavy flavors, this has been done in Ref. [84]. If oneincludes the transport of heavy flavors in the pre-thermalization stage, for consistency, oneshould also take into account the thermalization process of the medium between the initialhard collision and the hydrodynamics starting time. Parameter calibrations should be re-analyzed with free streaming replaced by a more realistic thermalization process. Finally,– 20 –e will extend the current study of bottomonium production to charmonium production inheavy-ion collisions. There, we must also initialize unbound c ¯ c pairs, since many of themare produced primordially and uncorrelated recombination is enhanced. Acknowledgments
XY would like to thank the organizers of the workshop EMMI Rapid Reaction Task Force:Suppression and (re)generation of quarkonium in heavy-ion collisions at the LHC and theorganizers of the workshop Quarkonia as Tools 2020, at both of which great discussions onquarkonium in-medium transport came out. This work is supported by U.S. Departmentof Energy (DOE) research grants DE-FG02-05ER41367. KW acknowledges support fromthe DOE grant DE-AC02-05CH11231 and National Science Foundation under the grantACI-1550228 within the JETSCAPE Collaboration. XY acknowledges support from theDOE grant DE-SC0011090, Brookhaven National Laboratory and Department of Physics,Massachusetts Institute of Technology.
A Dissociation and Recombination Terms in Transport Equations
We will follow most of the notations used in Ref. [68]. We will use k to label the momen-tum of the quarkonium involved in a scattering process. q will be used to indicate themomentum of the gluon absorbed or emitted by the quarkonium. This gluon can be real inthe process g + H ↔ Q + ¯ Q or virtual in inelastic scattering. If the gluon is on-shell, q = | q | will be used to represent its energy. For the process q + H ↔ q + Q + ¯ Q , p and p areused indicate the momenta of the light quarks on the left and right respectively. Similarlyfor the process g + H ↔ g + Q + ¯ Q , we will use q and q to represent the momenta ofthe gluons on the left and right respectively. In the quarkonium rest frame, the energy ofthe quarkonium is given by −| E nl | where E nl is the binding energy. In the rest frame ofthe unbound Q ¯ Q pair, its energy is given by p M where p rel is their relative momentum.The expressions shown below are valid in the rest frame of quarkonium for dissociation orthe center-of-mass frame of a Q ¯ Q pair for recombination. These two frames are not equiv-alent, but their difference is negligible (suppressed by TM ). In real simulations, the gluondistribution is boosted from the local rest frame of the hydro-cell (where temperature isdefined) to these rest frames of the heavy particles. After calculating the reaction rates inthe rest frames of the heavy particles, one has to boost them back to the laboratory framewhere we keep track of the phase space distributions. In the following, for simplicity, wewill only show expressions for quarkonia or Q ¯ Q pairs whose center-of-mass motions are atrest with respect to the medium local rest frame. In practice, we account for the Lorentzboost properly, as described above. – 21 –or later convenience, we define a “ δ − derivative” symbol, first introduced in Ref. [85] δδ p i (cid:90) n (cid:89) j =1 d p j (2 π ) h ( p , p , · · · , p n ) (cid:12)(cid:12)(cid:12) p i = p ≡ δδw ( p ) (cid:90) n (cid:89) j =1 d p j (2 π ) h ( p , p , · · · , p n ) w ( p i )= (cid:90) n (cid:89) j =1 ,j (cid:54) = i d p j (2 π ) h ( p , p , · · · , p i − , p , p i +1 , · · · , p n ) , (A.1)where the second δ denotes the standard functional variation and h ( p , p , · · · , p n ) and w ( p i ) are arbitrary independent smooth functions. A.1 g + H ↔ Q + ¯ Q The dissociation and recombination terms in the transport equation of the quarkoniumstate nls can be written as C ± nls ( x , p , t ) = δ F ± nls δ k (cid:12)(cid:12)(cid:12) k = p , (A.2)where we used the “ δ − derivative” symbol defined above. The functions F ± nls are definedby F + nls ≡ g + (cid:90) d k (2 π ) d p Q (2 π ) d p ¯ Q (2 π ) d q q (2 π ) (1 + n B ( q )) f Q ¯ Q ( x Q , p Q , x ¯ Q , p ¯ Q , t )(2 π ) δ ( k + q − p cm ) δ ( −| E nl | + q − p M ) (cid:88) |M| (A.3) F − nls ≡ (cid:90) d k (2 π ) d p cm (2 π ) d p rel (2 π ) d q q (2 π ) n B ( q ) f nls ( x , k , t )(2 π ) δ ( k + q − p cm ) δ ( −| E nl | + q − p M ) (cid:88) |M| , (A.4)where n B is the Bose-Einstein distribution function, g + = N c g s and g s is the degeneracyfactor for spin: g s = for a quarkonium state with spin s = 1 and for spin s = 0. Foran arbitrary Q ¯ Q pair, its probability of being a spin-1 state is . Its probability of being acolor octet is N c − N c . When computing the scattering amplitude squared, we have summedover the relevant quantum numbers (color of the Q ¯ Q , gluon polarization and color), sowe only use N c to avoid double counting. In the definition of F − nls , p cm = p Q + p ¯ Q , and p rel = p Q − p ¯ Q are the center-of-mass and relative momenta of a Q ¯ Q pair with momenta p Q and p ¯ Q .After summing over the relevant quantum numbers (see Ref. [68] for details), thescattering amplitude squared is given by (cid:88) |M| ≡ g C F q |(cid:104) Ψ p rel | r | ψ nl (cid:105)| , (A.5)where C F = N c − N c , | ψ nl (cid:105) is the wavefunction of the quarkonium state nls (states withdifferent spins are degenerate) and | Ψ p rel (cid:105) is the Coulomb scattering wave for the unbound– 22 – ¯ Q pair. For non-S wave states, we average over the polarizations:12 l + 1 l (cid:88) m l = − l (cid:90) d p rel (cid:104) ψ nlm l | r i | Ψ p rel (cid:105)(cid:104) Ψ p rel | r j | ψ nlm l (cid:105) = 13 δ ij l + 1 l (cid:88) m l = − l (cid:90) d p rel |(cid:104) Ψ p rel | r | ψ nlm l (cid:105)| ≡ δ ij (cid:90) d p rel |(cid:104) Ψ p rel | r | ψ nl (cid:105)| . (A.6) A.2 q + H ↔ q + Q + ¯ Q and g + H ↔ g + Q + ¯ Q For the inelastic scattering channels, we have C ± nls, inel ( x , p , t ) = δ F ± nls, ineq δ k (cid:12)(cid:12)(cid:12) k = p + δ F ± nls, ineg δ k (cid:12)(cid:12)(cid:12) k = p . (A.7)The inelastic F ± nls functions are defined by F + nls, ineq ≡ g + (cid:90) d k (2 π ) d p Q (2 π ) d p ¯ Q (2 π ) d p p (2 π ) d p p (2 π ) n F ( p )(1 − n F ( p )) f Q ¯ Q ( x Q , p Q , x ¯ Q , p ¯ Q , t )(2 π ) δ ( k + p − p cm − p ) δ ( −| E nl | + p − p M − p ) (cid:88) |M ineq | (A.8) F − nls, ineq ≡ (cid:90) d k (2 π ) d p cm (2 π ) d p rel (2 π ) d p p (2 π ) d p p (2 π ) n F ( p )(1 − n F ( p )) f nls ( x , k , t )(2 π ) δ ( k + p − p cm − p ) δ ( −| E nl | + p − p M − p ) (cid:88) |M ineq | (A.9) F + nls, ineg ≡ g + (cid:90) d k (2 π ) d p Q (2 π ) d p ¯ Q (2 π ) d q q (2 π ) d q q (2 π ) n B ( q )(1 + n B ( q )) f Q ¯ Q ( x Q , p Q , x ¯ Q , p ¯ Q , t )(2 π ) δ ( k + q − p cm − q ) δ ( −| E nl | + q − p M − q ) (cid:88) |M ineg | (A.10) F − nls, ineg ≡ (cid:90) d k (2 π ) d p cm (2 π ) d p rel (2 π ) d q q (2 π ) d q q (2 π ) n B ( q )(1 + n B ( q )) f nls ( x , k , t )(2 π ) δ ( k + q − p cm − q ) δ ( −| E nl | + q − p M − q ) (cid:88) |M ineg | , (A.11)where n F is the Fermi-Dirac distribution for fermions. The scattering amplitudes squared,with summation over relevant quantum numbers are given by (see Ref. [68] for details) (cid:88) |M ineq | = 163 g T F C F |(cid:104) Ψ p rel | r | ψ nl (cid:105)| p p + p · p q (A.12) (cid:88) |M ineg | = 13 g C F |(cid:104) Ψ p rel | r | ψ nl (cid:105)| q · ˆ q ) q ( q + q ) , (A.13)in which T F = , ˆ q i ≡ q i q i for i = 1 , |(cid:104) Ψ p rel | r | ψ nl (cid:105)| . – 23 – .3 |(cid:104) Ψ p rel | r | ψ nl (cid:105)| For non-S wave, we average over the third component of the orbital angular momentumand |(cid:104) Ψ p rel | r | ψ nl (cid:105)| is defined in Eq. (A.6). The partial wave expansion of the Coulombscattering wave | Ψ p rel (cid:105) is given byΨ p rel ( r ) = (cid:104) r | Ψ p rel (cid:105) = 4 π (cid:88) (cid:96),m i (cid:96) e iδ (cid:96) F (cid:96) ( ρ ) ρ Y (cid:96)m (ˆ r ) Y ∗ (cid:96)m (ˆ p rel ) , (A.14)in which Y (cid:96)m denotes the spherical harmonics and ρ = p rel r (A.15) δ (cid:96) = arg Γ(1 + (cid:96) + iη ) (A.16) η = α s M N c p rel (A.17) F (cid:96) ( ρ ) = 2 (cid:96) e − πη/ | Γ(1 + (cid:96) + iη ) | (2 (cid:96) + 1)! ρ (cid:96) +1 e iρ F ( (cid:96) + 1 + iη ; 2 (cid:96) + 2; − iρ ) , (A.18)where F ( a ; b ; z ) is the confluent hypergeometric function. We use the Wigner-Eckarttheorem to simplify the calculations.For 1 S state, we have |(cid:104) Ψ p rel | r | ψ S (cid:105)| = 2 π ηa B p (1 + η )(2 + ηa B p rel ) (1 + a B p ) ( e πη − e η arctan ( a B p rel ) , (A.19)where a B = α s C F M is the Bohr radius and η is defined in Eq. (A.17).For 2 S state, we have |(cid:104) Ψ p rel | r | ψ S (cid:105)| = 2 π ηa B p (1 + η )(1 + 4 a B p ) ( e πη − e η arctan (2 a B p rel ) ( − − ηa B p rel + 8 a B p − η a B p + 4 ηa B p ) . (A.20)This expression can be further simplified by using ηa B p rel = N c − . These matrix elementsof 1 S and 2 S have been reported in Ref. [36].For 1 P state, we find |(cid:104) Ψ p rel | r | ψ P (cid:105)| = 2 π ηa B a B p ) ( e πη − e η arctan (2 a B p rel ) (cid:104)(cid:0) η ( η − a B p + 12(2 η − a B p + 18 ηa B p rel + 3 (cid:1) +16 a B p (3 + 2 ηa B p rel ) ( η + 5 η + 4) (cid:105) . (A.21)We do not need the dipole transition matrix elements for 3 S and 2 P states, since weassume they cannot be formed inside the QGP. In other words, if a 3 S or 2 P state entersthe QGP, it melts immediately. In practice, we force them to dissociate in the code when T >
154 MeV and do not allow them to be (re)generated inside the QGP.– 24 –
Details on Feed-down Contributions
We will use the notations σ tot nl and σ nl to denote the total and primordial cross sections ina nucleon-nucleon collision. The former contains the latter and feed-down contributions.To work out the feed-down contributions from excited states to the ground and lower-excited states, the first thing we need is the branching ratio in vacuum. The most recentresults reported by the Particle Data Group [86] are summarized in Table 1. For the P-wavestates, we need the averaged branching ratio since our transport equations are degeneratein spin. To this end, we follow Ref. [21] and assume for both 1 P and 2 P states σ χ b σ χ b = 0 .
85 (B.1) σ χ b σ χ b = 1 . . (B.2)Then we can work out the averaged branching ratiosBr[ χ b (1 P ) → Υ(1 S )] ≡ Br[ χ b (1 P ) → Υ(1 S )] σ χ b + Br[ χ b (1 P ) → Υ(1 S )] σ χ b + Br[ χ b (1 P ) → Υ(1 S )] σ χ b σ χ b + σ χ b + σ χ b = 0 .
159 (B.3)Br[ χ b (2 P ) → Υ(1 S )] ≡ Br[ χ b (2 P ) → Υ(1 S )] σ χ b + Br[ χ b (2 P ) → Υ(1 S )] σ χ b + Br[ χ b (2 P ) → Υ(1 S )] σ χ b σ χ b + σ χ b + σ χ b = 0 .
056 (B.4)Br[ χ b (2 P ) → Υ(2 S )] ≡ Br[ χ b (2 P ) → Υ(2 S )] σ χ b + Br[ χ b (2 P ) → Υ(2 S )] σ χ b + Br[ χ b (2 P ) → Υ(2 S )] σ χ b σ χ b + σ χ b + σ χ b = 0 .
083 (B.5)The next thing we need is the primordial cross section ratio. Following Ref. [21] andreferences therein, we have σ tot2 S = 0 . σ tot1 S (B.6) σ S = 0 . σ tot1 S (B.7) σ P = 1 . σ tot1 S (B.8) σ P = 0 . σ tot1 S . (B.9)The total cross sections of 1 S and 2 S can be written as σ tot1 S = σ S + σ tot2 S Br[Υ(2 S ) → Υ(1 S )] + σ S Br[Υ(3 S ) → Υ(1 S )]+ σ P Br[ χ b (1 P ) → Υ(1 S )] + σ P Br[ χ b (2 P ) → Υ(1 S )] (B.10) σ tot2 S = σ S + σ S Br[Υ(3 S ) → Υ(2 S )] + σ P Br[ χ b (2 P ) → Υ(2 S )] . (B.11)Using Eqs. (B.6, B.7, B.8, B.9) and the branching ratios, we find σ S = 0 . σ tot1 S (B.12) σ S = 0 . σ tot1 S . (B.13)– 25 –hannel Branching ratioΥ(2 S ) → Υ(1 S ) 0.265Υ(3 S ) → Υ(1 S ) 0.066 χ b (1 P ) → Υ(1 S ) 0.019 χ b (1 P ) → Υ(1 S ) 0.352 χ b (1 P ) → Υ(1 S ) 0.180 χ b (2 P ) → Υ(1 S ) 0.004 χ b (2 P ) → Υ(1 S ) 0.115 χ b (2 P ) → Υ(1 S ) 0.077Υ(3 S ) → Υ(2 S ) 0.106 χ b (2 P ) → Υ(2 S ) 0.014 χ b (2 P ) → Υ(2 S ) 0.181 χ b (2 P ) → Υ(2 S ) 0.089 Table 1 : Branching ratios of different bottomonium states.Primordial cross section ratio Value σ S /σ S σ S /σ S σ P /σ S σ P /σ S σ S /σ S σ P /σ S σ P /σ S σ S /σ P σ P /σ P Table 2 : Ratios of primordial cross sections.Experimentally it is known that the feed-down contributions to σ tot1 S is about 67% (averagedover the transverse momentum), consistent with our prescription here. Now we are readyto compute the ratios of the primordial cross sections. The results are listed in Table. 2. References [1] C. Quigg and J. L. Rosner, Phys. Lett. , 153 (1977).[2] T. Matsui and H. Satz, Phys. Lett. B , 416 (1986).[3] F. Karsch, M. T. Mehr and H. Satz, Z. Phys. C , 617 (1988).[4] L. D. McLerran and B. Svetitsky, Phys. Rev. D , 450 (1981).[5] A. Mocsy and P. Petreczky, Phys. Rev. Lett. , 211602 (2007) [arXiv:0706.2183 [hep-ph]].[6] M. Laine, O. Philipsen, P. Romatschke and M. Tassler, JHEP , 054 (2007)[hep-ph/0611300]. – 26 –
7] A. Beraudo, J.-P. Blaizot and C. Ratti, Nucl. Phys. A , 312 (2008) [arXiv:0712.4394[nucl-th]].[8] R. L. Thews, M. Schroedter and J. Rafelski, Phys. Rev. C , 054905 (2001)[hep-ph/0007323].[9] A. Andronic, P. Braun-Munzinger, K. Redlich and J. Stachel, Phys. Lett. B , 259 (2007)[nucl-th/0701079 [NUCL-TH]].[10] L. Grandchamp, R. Rapp and G. E. Brown, Phys. Rev. Lett. , 212301 (2004)[hep-ph/0306077].[11] L. Grandchamp, S. Lumpkins, D. Sun, H. van Hees and R. Rapp, Phys. Rev. C , 064906(2006) [hep-ph/0507314].[12] L. Yan, P. Zhuang and N. Xu, Phys. Rev. Lett. , 232301 (2006) [nucl-th/0608010].[13] Y. Liu, Z. Qu, N. Xu and P. Zhuang, Phys. Lett. B , 72 (2009) [arXiv:0901.2757[nucl-th]].[14] T. Song, K. C. Han and C. M. Ko, Phys. Rev. C , 034907 (2011) [arXiv:1103.6197[nucl-th]].[15] T. Song, K. C. Han and C. M. Ko, Phys. Rev. C , 014902 (2012) [arXiv:1109.6691[nucl-th]].[16] R. Sharma and I. Vitev, Phys. Rev. C , no. 4, 044905 (2013) [arXiv:1203.0329 [hep-ph]].[17] F. Nendzig and G. Wolschin, J. Phys. G , 095003 (2014) [arXiv:1406.5103 [hep-ph]].[18] B. Krouppa, R. Ryblewski and M. Strickland, Phys. Rev. C , no. 6, 061901 (2015)[arXiv:1507.03951 [hep-ph]].[19] B. Chen and J. Zhao, Phys. Lett. B , 819 (2017) [arXiv:1704.05622 [nucl-th]].[20] J. Zhao and B. Chen, Phys. Lett. B , 17 (2018) [arXiv:1705.04558 [nucl-th]].[21] X. Du, R. Rapp and M. He, Phys. Rev. C , no. 5, 054901 (2017) [arXiv:1706.08670[hep-ph]].[22] S. Aronson, E. Borras, B. Odegard, R. Sharma and I. Vitev, Phys. Lett. B , 384 (2018)[arXiv:1709.02372 [hep-ph]].[23] E. G. Ferreiro and J. P. Lansberg, JHEP , 094 (2018) [arXiv:1804.04474 [hep-ph]].[24] X. Yao, W. Ke, Y. Xu, S. Bass and B. M¨uller, Nucl. Phys. A , 755 (2019)[arXiv:1807.06199 [nucl-th]].[25] J. Hong and S. H. Lee, Phys. Lett. B , 135147 (2020) [arXiv:1909.07696 [nucl-th]].[26] B. Chen, M. Hu, H. Zhang and J. Zhao, Phys. Lett. B , 135271 (2020) [arXiv:1910.08275[nucl-th]].[27] O. Kaczmarek, F. Karsch, P. Petreczky and F. Zantow, Phys. Lett. B , 41 (2002)[hep-lat/0207002].[28] A. Bazavov et al. [TUMQCD Collaboration], Phys. Rev. D , no. 5, 054511 (2018)[arXiv:1804.10600 [hep-lat]].[29] Y. Burnier, O. Kaczmarek and A. Rothkopf, Phys. Rev. Lett. , no. 8, 082001 (2015)[arXiv:1410.2546 [hep-lat]].[30] M. E. Peskin, Nucl. Phys. B , 365 (1979). – 27 –
31] G. Bhanot and M. E. Peskin, Nucl. Phys. B , 391 (1979).[32] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Nucl. Phys. B , 275 (2000)[hep-ph/9907240].[33] N. Brambilla, A. Pineda, J. Soto and A. Vairo, Rev. Mod. Phys. , 1423 (2005)[hep-ph/0410047].[34] S. Fleming and T. Mehen, Phys. Rev. D , 034502 (2006) [hep-ph/0509313].[35] N. Brambilla, J. Ghiglieri, A. Vairo and P. Petreczky, Phys. Rev. D , 014017 (2008)[arXiv:0804.0993 [hep-ph]].[36] N. Brambilla, M. A. Escobedo, J. Ghiglieri and A. Vairo, JHEP , 116 (2011)[arXiv:1109.5826 [hep-ph]].[37] N. Brambilla, M. A. Escobedo, J. Ghiglieri and A. Vairo, JHEP , 130 (2013)[arXiv:1303.6097 [hep-ph]].[38] T. Binder, B. Blobel, J. Harz and K. Mukaida, [arXiv:2002.07145 [hep-ph]].[39] A. Dumitru, Y. Guo and M. Strickland, Phys. Lett. B , 37 (2008) [arXiv:0711.4722[hep-ph]].[40] A. Dumitru, Y. Guo and M. Strickland, Phys. Rev. D , 114003 (2009) [arXiv:0903.4703[hep-ph]].[41] Q. Du, A. Dumitru, Y. Guo and M. Strickland, JHEP , 123 (2017) [arXiv:1611.08379[hep-ph]].[42] H. Liu, K. Rajagopal and U. A. Wiedemann, Phys. Rev. Lett. , 182301 (2007)[hep-ph/0607062].[43] M. A. Escobedo, J. Soto and M. Mannarelli, Phys. Rev. D , 016008 (2011)[arXiv:1105.1249 [hep-ph]].[44] X. Yao and T. Mehen, Phys. Rev. D , no. 9, 096028 (2019) [arXiv:1811.07027 [hep-ph]].[45] X. Yao, W. Ke, Y. Xu, S. A. Bass, T. Mehen and B. M¨uller, [arXiv:2002.04079 [hep-ph]].[46] C. Young and K. Dusling, Phys. Rev. C , 065206 (2013) [arXiv:1001.0935 [nucl-th]].[47] N. Borghini and C. Gombeaud, Eur. Phys. J. C , 2000 (2012) [arXiv:1109.4271 [nucl-th]].[48] Y. Akamatsu and A. Rothkopf, Phys. Rev. D , 105011 (2012) [arXiv:1110.1203 [hep-ph]][49] Y. Akamatsu, Phys. Rev. D , 056002 (2015) [arXiv:1403.5783 [hep-ph]].[50] J. P. Blaizot, D. De Boni, P. Faccioli and G. Garberoglio, Nucl. Phys. A , 49 (2016)[arXiv:1503.03857 [nucl-th]].[51] R. Katz and P. B. Gossiaux, Annals Phys. , 267 (2016) [arXiv:1504.08087 [quant-ph]].[52] S. Kajimoto, Y. Akamatsu, M. Asakawa and A. Rothkopf, Phys. Rev. D , no. 1, 014003(2018) [arXiv:1705.03365 [nucl-th]].[53] D. De Boni, JHEP , 064 (2017) [arXiv:1705.03567 [hep-ph]].[54] J. P. Blaizot and M. A. Escobedo, JHEP , 034 (2018) [arXiv:1711.10812 [hep-ph]].[55] J. P. Blaizot and M. A. Escobedo, Phys. Rev. D , no. 7, 074007 (2018) [arXiv:1803.07996[hep-ph]]. – 28 –
56] Y. Akamatsu, M. Asakawa, S. Kajimoto and A. Rothkopf, JHEP , 029 (2018)[arXiv:1805.00167 [nucl-th]].[57] T. Miura, Y. Akamatsu, M. Asakawa and A. Rothkopf, Phys. Rev. D , no. 3, 034011(2020) [arXiv:1908.06293 [nucl-th]].[58] R. Sharma and A. Tiwari, Phys. Rev. D , no.7, 074004 (2020) [arXiv:1912.07036[hep-ph]].[59] N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, Phys. Rev. D , no. 3, 034021 (2017)[arXiv:1612.07248 [hep-ph]].[60] N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, Phys. Rev. D , no. 7, 074009 (2018)[arXiv:1711.04515 [hep-ph]].[61] N. Brambilla, M. A. Escobedo, A. Vairo and P. Vander Griend, Phys. Rev. D , no. 5,054025 (2019) [arXiv:1903.08063 [hep-ph]].[62] P. B. Gossiaux and J. Aichelin, Phys. Rev. C , 014904 (2008) [arXiv:0802.2525 [hep-ph]].[63] P. B. Gossiaux, R. Bierkandt and J. Aichelin, Phys. Rev. C , 044906 (2009)[arXiv:0901.0946 [hep-ph]].[64] J. Uphoff, O. Fochler, Z. Xu and C. Greiner, J. Phys. G , no. 11, 115106 (2015)[arXiv:1408.2964 [hep-ph]].[65] S. Cao, T. Luo, G. Y. Qin and X. N. Wang, Phys. Rev. C , no. 1, 014909 (2016)[arXiv:1605.06447 [nucl-th]].[66] W. Ke, Y. Xu and S. A. Bass, Phys. Rev. C , no. 6, 064901 (2018) [arXiv:1806.08848[nucl-th]].[67] Y. Xu, J. E. Bernhard, S. A. Bass, M. Nahrgang and S. Cao, Phys. Rev. C , no. 1, 014907(2018) [arXiv:1710.00807 [nucl-th]].[68] X. Yao and B. M¨uller, Phys. Rev. D , no. 1, 014008 (2019) [arXiv:1811.09644 [hep-ph]].[69] G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D , 1125 (1995) Erratum: [Phys.Rev. D , 5853 (1997)] [hep-ph/9407339].[70] X. Yao and B. M¨uller, Phys. Rev. C , no. 1, 014908 (2018) Erratum: [Phys. Rev. C ,no. 4, 049903 (2018)] [arXiv:1709.03529 [hep-ph]].[71] T. Sj¨ostrand et al. , Comput. Phys. Commun. , 159 (2015).[72] K. J. Eskola, P. Paakkinen, H. Paukkunen and C. A. Salgado, Eur. Phys. J. C , no. 3, 163(2017) [arXiv:1612.05741 [hep-ph]].[73] P. Wang [STAR Collaboration], Nucl. Phys. A , 723 (2019).[74] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu and U. Heinz, Phys. Rev. C , no. 2,024907 (2016) [arXiv:1605.03954 [nucl-th]].[75] H. Song and U. W. Heinz, Phys. Rev. C , 064901 (2008) [arXiv:0712.3715 [nucl-th]].[76] C. Shen, Z. Qiu, H. Song, J. Bernhard, S. Bass and U. Heinz, Comput. Phys. Commun. ,61 (2016) [arXiv:1409.8164 [nucl-th]].[77] J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu and U. Heinz, Phys. Rev. C , no. 2,024907 (2016) [arXiv:1605.03954 [nucl-th]]. – 29 –
78] A. M. Sirunyan et al. [CMS Collaboration], Phys. Lett. B , 270 (2019). [arXiv:1805.09215[hep-ex]].[79] V. Khachatryan et al. [CMS Collaboration], Phys. Lett. B , 357 (2017) [arXiv:1611.01510[nucl-ex]].[80] S. Acharya et al. [ALICE Collaboration], Phys. Rev. Lett. , no. 19, 192301 (2019).[81] CMS Collaboration, CMS-PAS-HIN-19-002; J. Park, Quark Matter 2019.[82] J. E. Bernhard, J. S. Moreland and S. A. Bass, Nature Phys. , no. 11, 1113 (2019).[83] X. Yao, ProQuest Number: 13882375 (2019) [arXiv:1911.08500 [nucl-th]].[84] S. Mrowczynski, Eur. Phys. J. A , no. 3, 43 (2018) [arXiv:1706.03127 [nucl-th]].[85] X. Yao and B. M¨uller, Phys. Rev. D , no. 7, 074003 (2018) [arXiv:1801.02652 [hep-ph]].[86] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D , no. 3, 030001 (2018)., no. 3, 030001 (2018).