Coupled Boltzmann Transport Equations of Heavy Quarks and Quarkonia in Quark-Gluon Plasma

We develop a framework of coupled transport equations for open heavy flavor and 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 flavor observables in heavy-ion collision experiments and can account for both, uncorrelated and correlated recombination. Our recombination implementation depends on real-time open heavy quark and antiquark distributions. We carry out consistency tests to show how the interplay among open heavy flavor transport, quarkonium dissociation and recombination drives the system to equilibrium. We then apply our framework to study bottomonium production in heavy-ion collisions. We include $\Upsilon(1S)$, $\Upsilon(2S)$, $\Upsilon(3S)$, $\chi_b(1P)$ and $\chi_b(2P)$ in the framework and take feed-down contributions during the hadronic gas stage into account. Cold nuclear matter effects are included by using nuclear parton distribution functions for the initial primordial heavy flavor production. A calibrated $2+1$ dimensional viscous hydrodynamics is used to describe the bulk QCD medium. We calculate both the nuclear modification factor $R_{\mathrm{AA}}$ of all bottomonia states and the azimuthal angular anisotropy coefficient $v_2$ of the $\Upsilon(1S)$ state and find that our results agree reasonably with experimental measurements. Our calculations indicate that correlated cross-talk recombination is an important production mechanism of bottomonium in current heavy-ion experiments. The importance of correlated recombination can be tested experimentally by measuring the ratio of $R_{\mathrm{AA}}(\chi_b(1P))$ and $R_{\mathrm{AA}}(\Upsilon(2S))$.


Introduction
Heavy quarkonia are bound states of heavy quark-antiquark pairs QQ. The mass spectra of the ground and lower excited quarkonium states can be reasonably well described by the nonrelativistic Schrödinger equation with a potential model [1]. Inside a hot nuclear medium, i.e., the quark-gluon plasma (QGP), the attractive potential can be significantly suppressed due to the static screening effect in the plasma [2]. As a result, a bound QQ can "melt" at sufficiently high temperature [3][4][5]. Therefore, quarkonium suppression (compared to a non QGP baseline) can be used as a signal of the formation of a QGP in heavy-ion collisions. In general, shallower bound states melt at lower temperatures and one would expect a "sequential" suppression pattern.
However, this simple picture is complicated by other in-medium processes such as the dissociation of quarkonium states via dynamical scattering (i.e. the dynamical screening effect, which is related to the imaginary part of the QQ potential [6,7]) and heavy quark (re)combination [8,9]. To account for all these effects, semi-classical transport equations have been widely applied . These calculations typically contain three components: The first part is a temperature-dependent potential that is parametrized from lattice calculations of the free energy of a QQ singlet [32,33] or direct lattice calculations of the real part of the QQ potential [34]. The second input is the dissociation rate of each quarkonium state. Perturbative calculations of dissociation rates include both the gluo-dissociation process [35,36] and the inelastic scattering with the medium (Landau damping). An effective field theory of QCD, potential nonrelativistic QCD (pNRQCD) [37][38][39] has been applied to study both the static screening of the potential and the dissociation rate [40][41][42]. Similar effective theory has also been used to study dark matter bound state formation [43][44][45]. Some studies also included viscous effects [46][47][48] and modifications due to the quarkonium moving with respect to a static medium [49,50]. The final component is a recombination model, for example, a statistical hadronization model [9], a coalescence model based on Wigner functions [22] or a model that is based on detailed balance and a finite relaxation rate [24].
Unlike the first two effects, recombination contains significant model-dependencies in most studies. Here, we shall distinguish between two kinds of recombination: uncorrelated and correlated recombination. In uncorrelated recombination, the Q andQ originate from differential initial hard vertices. In proton-proton collisions, these heavy quarks and antiquarks would almost never (re)combine to form a quarkonium state due to their separation in phase space. However, in heavy-ion collisions, multiple QQ pairs are produced from the initial hard scatterings in one collision event. The momenta of these heavy quarks and antiquarks change continuously inside the QGP due to diffusion and energy loss. Therefore, the chance for an uncorrelated pair to come close to each other in phase space is higher. When they are close, they may combine to form a quarkonium state. The uncorrelated recombination rate rises with the number of open heavy quarks produced in the collision, which is crucial to explain the small amount of suppression observed for the J/ψ with rising collision energy. Naively, one would expect J/ψ to be more suppressed at higher collision energies due to the hotter medium and stronger plasma screening effects. Therefore, uncorrelated recombination is important for the phenomenology of charmonium. However, its effect on bottomonium production is expected to be negligible, since only a few bottom-antibottom quark pairs are produced in one collision.
In correlated recombination, the Q andQ originate from the same initial hard vertex. For example, a QQ pair in a pre-quarkonium resonance or emerging from a previous dissociation of a quarkonium state is considered a correlated pair. The phenomenological effect of correlated recombination for both charmonium and bottomonium has not been systematically explored yet.
Most recombination models depend on the uncorrelated open heavy quark distributions, they can thus only account for uncorrelated recombination. To explore the physical effects of correlated recombination, we develop a set of coupled transport equations for Q andQ's as well as for quarkonia, which allows us to study both, uncorrelated and correlated recombination. The transport equation of quarkonium provides consistent treatment of dissociation and recombination, in the sense that both of them are derived from QCD under systematic nonrelativistic expansions [51,52]. The derivation is based on the combination of the open quantum system framework (in which correlated recombination is taken into account) and the pNRQCD effective field theory. The application of the open quantum system framework to studying quarkonium in-medium dynamics [53][54][55][56][57][58][59][60][61][62][63][64][65] and its combination with pNRQCD [66,67] has recently drawn significant theoretical interest. A quarkonium transport coefficient has been defined in Ref. [68]. In the open quantum system approach, quarkonium dissociation is caused by the wavefunction decoherence, which also leads to correlated recombination at the same time. For example, if we start with the 1S state: |ψ(t = 0) = |1S , the wavefunction decoherence leads to | 1S|ψ(t) | 2 < 1 for t > 0. But at the same time, if the 2S state exists as a well-defined bound state, we will also have | 2S|ψ(t) | 2 > 0, i.e., some 2S state is regenerated. Since our transport equation for quarkonium is derived from the open quantum system approach, it can handle correlated recombination.
The coupled transport equations allow us to study the in-medium transport of both open and hidden heavy flavor states. In this paper, we will use this framework to study bottomonium production in heavy-ion collisions. We will include Υ(1S), Υ(2S), Υ(3S), χ b (1P ) and χ b (2P ) states in the transport network. By solving the coupled transport equations via test particle Monte Carlo simulations, we will explore the importance of correlated recombination in bottomonium phenomenology. This paper is organized as follows: in Sect. 2, we will introduce the set of coupled transport equations. Then in Sect. 3, some simulation tests inside a QGP box will be shown and compared to the system properties in thermal equilibrium. Details on applying the transport equations to the study of heavy-ion collisions will be explained in Sect. 4. Results on the nuclear modification factor (R AA ) and the azimuthal angular anisotropy coefficient v 2 of bottomonia will be discussed later in Sect. 5. Finally, we will draw conclusions in Sect. 6.

Coupled Transport Equations
The set of coupled Boltzmann transport equations for the distribution functions of unbound heavy quark-antiquark pairs QQ 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 whereẋ = ∂x ∂t . The left-hand sides of these equations describe the free streaming of distribution functions in phase space while the right-hand sides contain collision terms of the heavy particles interacting with the plasma. The collision terms with ± superscripts represent the quarkonium dissociation (−) and recombination (+) while the term without any superscript describes the energy and momentum changes of the open heavy quarkantiquark pairs. In the following, we will explain these collision terms in detail.

Transport of Open Heavy Quark-Antiquark Pairs
If we neglect the interaction between the heavy quark-antiquark pair, we can write i.e., the heavy quark and antiquark interact independently with the medium. In potential models, the interaction between the QQ pair is attractive for the color singlet and repulsive for the color octet. For a Coulomb potential, the color-averaged potential vanishes if we assume the color is in thermal equilibrium. Furthermore, in-medium QQ potentials are significantly suppressed, so we expect the potential interaction between the QQ pair to be weak. Therefore, throughout this paper, we will neglect the interaction between the QQ pair in the open heavy flavor transport equations. As a remark, we note that . But the opposite is not true in general. Most recombination models implicitly assume the factorization of the QQ distributions and thus cannot study the correlated 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 [69][70][71][72]. The interaction of open heavy quarks (and antiquarks) with the medium is described by 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 → 2 scattering q + Q → q + Q and g + Q → g + Q, the inelastic 2 → 3 scattering q + Q → q + Q + g and g + Q → g + Q + g and the inelastic 3 → 2 scattering q + Q + g → q + Q and g + Q + g → g + Q, and similarly for CQ. In this work, we will use the Monte Carlo simulations in the Lido package [73] to solve the transport equations of open heavy quark-antiquark pairs. The Lido package contains both a linearized Boltzmann transport description and a model that is based on the Langevin equation with radiation corrections. The latter description has been reported in Ref. [74]. We will only use the linearized Boltzmann description in this work.

Transport of Quarkonia
For the dissociation and recombination terms in the transport equations, we will use the expressions in Ref. [75]. Detailed expressions of the relevant collision terms can be found in Appendix A. The calculations therein are based on a version of pNRQCD under the hierarchy of scales Here M is the heavy quark mass, v the typical relative velocity between the QQ pair inside the bound state, T the temperature of the QGP and m D the Debye mass. The typical size of quarkonium is roughly give by r ∼ 1 M v and the typical binding energy is about M v 2 . The last inequality T m D means the QGP is weakly-coupled. For charmonium, we have v 2 ∼ 0.3 while for bottomonium, v 2 ∼ 0.1 [76]. They both give M v 2 ∼ 500 MeV.
One may worry that the hierarchy is not always true in real heavy-ion collisions. In the early time of the QGP expansion, the temperature can be 450 MeV. Due to the static screening effect, the in-medium binding energies of quarkonium states can be much smaller than their vacuum binding energies, especially for excited states. Thus, our assumed hierarchy indeed breaks down in the early stage and in principle one has to use a different version of pNRQCD if there still exists a hierarchy of scales. However, even before the breakdown of the assumed hierarchy, as the temperature increases, the dissociation rates of excited quarkonium states such as Υ(2S) and χ b (1P) blow up rapidly. What happens in our calculations is that after a very short time period in the early stage, the excited quarkonium states have dissociated and evolve then as an unbound, correlated QQ pair. This is the right physics: When our hierarchy of scales breaks down, we can either have We do not consider M v T M v 2 here because the real values of v do not allow this hierarchy to happen for both charmonium and bottomonium. (For to be valid, one at least needs a factor of three in the ratio.) So we only need to consider T M v M v 2 here. Whenever this occurs, one expects both the plasma screening effects to be extremely strong ( 1 M v gives the rough size of the quarkonium state) 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 QQ pair. Then the correct description is the transport of an unbound QQ pair. Therefore, in our 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 QQ pair. This is a good approximation of the correct physics. The transport of excited quarkonia as well-defined bound states is valid again in the later stage of the evolution, when the temperature drops and our hierarchy is resumed. In this sense, most excited quarkonium states are probably generated via recombination in the later stage of the QGP evolution. Their suppression mechanism is mainly the decorrelation of the QQ pair in coordinate and momentum space (or decoherence of their wavefunction). Those excited states that are observed may probably come from recombination of QQ pairs that are 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 expansion (the multipole expansion under this hierarchy of scales is equivalent to a nonrelativistic expansion). At this order, the higher Fock state |QQg of quarkonium (in which the QQ is a color octet) is suppressed at least by two powers of v with respect to the leading Fock state |QQ (in which the QQ is a color singlet) [76]. Therefore, in our calculations, quarkonium is always a QQ pair in the color singlet. Unbound QQ pairs can be in the color singlet or octet states. At the order of the nonrelativistic expansion that we are studying right now, the transition between a quarkonium and an unbound pair only occurs via a dipole interaction between the color singlet and octet states. In other words, the dissociation product is a QQ in the color octet state and only color octet QQ pairs can recombine as quarkonia via the dipole interaction. Keeping only the dipole interaction in the calculation works here because of the hierarchy M v T . When the quarkonium or the QQ pair is at rest with respect to the local medium, the typical energy of a medium parton is E g ∼ T . The typical size of quarkonium is given by r ∼ 1 M v . So the quarkonium size is small in the sense that rE g 1 and thus the dipole vertex is a weak-coupling interaction. However, if the quarkonium or the QQ is moving with respect to the local medium (for example, when the quarkonium state has a finite transverse momentum), the typical energy of medium partons in the quarkonium (or QQ) rest frame is boosted, E g ∼ γT where γ is the boost factor that depends on the relative velocity between the quarkonium (or the center-of-mass of the QQ) and the local medium. The condition rE g 1 may no longer be true and our calculations would break down then. Rigorously speaking, our calculations only apply to quarkonium states at low transverse momentum. In practice, we should be careful when interpreting our results in the mid and high transverse momentum regions. We will come back to this issue later in Sect. 5.
The potentials used in the calculations are Coulombic: 2Nc and N c = 3. We will take α pot s in the potentials to be a parameter and choose its value as α pot s = 0.36. The coupling constant in the scattering vertices will be taken to be constant α s = 0.3. We will vary these two coupling constants and discuss the calculation uncertainties in Section 5. The effects of running coupling and nonperturbative potentials will be left to future studies 1 . Since the octet potential is non-zero here, the wavefunction of the unbound octet pair is a Coulomb scattering wave rather than a plane wave. In this way, we re-sum an infinite number of Coulomb exchanges between the octet pair in the initial (for recombination) or final (for dissociation) state. At the leading order in the nonrelativistic expansion, the potential is independent of the orbital angular momentum and spin. Thus the dissociation rates are the same for states separated by fine and hyperfine splittings. The recombination rates are also the same up to a spin degeneracy factor g s = 1 4 or 3 4 . 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 any quarkonium 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 is neglected here because it occurs at higher orders in the multipole expansion (or equivalently here, nonrelativistic expansion) [75]. The direct transitions between different quarkonia species are also omitted due to the same reason. Expressions of the reaction rates of the 1S, 2S and 1P states can be found in Appendix A. For 3S and 2P states, we assume they cannot exist inside the QGP. In other words, they dissociate immediately after entering the 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 QQ pair rather than the distribution of a singlet heavy quark). The correlated recombination leads 1 The different values of α pot s and αs chosen here already hint at the importance of nonperturbative potentials.
to cross-talk between different quarkonium states. When an excited quarkonium state such as Υ(2P ) or Υ(3S) dissociates, the produced QQ can form a lower excited state or the ground state Υ(1S). As will be discussed later, this is an important production mechanism for the ground state. When an excited state dissociates at high temperature, the dissociated QQ pair may form the ground state and then survive the subsequent evolution. When the temperature drops and a ground state dissociates occasionally (the dissociation rate is still non-vanishing even if the static screening effect is small at low temperature), the dissociated QQ pair may form an excited state. The time evolution of the whole system is a network of reactions among unbound heavy quarks, antiquarks and all quarkonia states. Our framework can handle this reaction network and study the physical impacts of the cross-talk recombination.

Monte Carlo Simulations
We solve the coupled transport equations by test particle Monte Carlo simulations. We sample a certain number of QQ pairs and quarkonia according to the distribution functions at the initial time. Mathematically, the distribution function for each particle species is represented by and similarly for the two-particle distribution functions. Here x i and p i are the position and momentum of the i-th sampled test particle. The integral of a distribution function over the whole phase space gives the total number of the particle associated with that distribution. The positions and momenta of the sampled particles obey their initial distribution functions. Then we evolve the positions and momenta of all particle species step by step. The time step size is chosen as ∆t = 0.01 fm/c in the laboratory frame. At each time step, we consider each of the following processes.

Free Streaming
The position of the particle changes according to where E(t) and p(t) are the energy and momentum of the particle at the current time step.

Momentum Change
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 heavy quark (antiquark) with the medium in each scattering channel. If a certain scattering process occurs, the package generates a light scattering partner utilizing local medium properties (temperature and flow field). Its outgoing momentum is sampled from the differential scattering rate. The outgoing momentum of the heavy quark (antiquark) can be obtained from energy and momentum conservation. Then we can update the momentum of the heavy quark (antiquark).

Dissociation
For a quarkonium state nls with a certain momentum (velocity) and a position at some local temperature, we calculate its dissociation rate in the laboratory frame. The method to obtain the dissociation rate from the collision term C − can be found in Ref. [75]. The rate times 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 incoming and/or outgoing momenta of the relevant light particles and obtain the momenta of the outgoing QQ pair from energy and momentum conservation. The positions of the unbound Q andQ are given by the position of the quarkonium before dissociation. Then we remove this quarkonium state from the relevant particle list and add the produced QQ pair to the list of heavy quarks and antiquarks.

Recombination
For each unbound QQ pair, we need to calculate their recombination rate. In the Monte Carlo simulation, the two-particle distribution function of the unbound QQ pair is represented by The recombination rate of a specific pair can be obtained by the following replacement where x cm , p cm , x rel and p rel are the center-of-mass (cm) and relative positions and momenta. This Gaussian ansatz is motivated by the recombination formula derived in Ref. [51], in which the recombination term for a QQ far away from each other is suppressed exponentially by the bound state wavefunction. The width of the Gaussian is chosen to be the typical size of the bound state. More specifically, for 1S, σ = a B where a B = 2 αsC F M is the Bohr radius. For 2S and 1P state, we have σ = 2a B . The recombination rate of this specific pair can be obtained by doing the following replacement in the recombination term C + in the Boltzmann equation (The expression of C + can be found in Appendix A.) In practice, for each unbound QQ pair with positions x i ,x j and momenta p i ,p j , we first boost the pair into their cm frame. Then we compute their relative momentum and calculate their recombination rate in the cm frame and then boost the rate back into the laboratory. If the recombination into a specific quarkonium state occurs, we sample the momentum of the outgoing quarkonium state based on the differential rate and energymomentum conservation. The position of the quarkonium state is given by the cm position of the unbound pair before recombination. Finally we remove the unbound pair from the list of heavy quarks and antiquarks and add a new quarkonium state to the relevant list of quarkonium states.

Test Simulation inside QGP Box
Before we study quarkonium production in heavy-ion collisions, we validate our test particle Monte Carlo simulations for the coupled transport equations. We solve the equations inside a cubic volume of QGP matter at fixed temperature with a side length L = 10 fm. The box has periodic boundary conditions, i.e., when a particle reaches the boundary of the box, it appears on the opposite side. In other words, the medium behaves as an infinitely large QGP with a finite heavy flavor density (which include both open and hidden heavy flavor states).
We focus on the bottom system, since the nonrelativistic expansion works better for bottom than for charm. We will sample a fixed number of unbound bb pairs initially. Their positions are randomly distributed in the volume while their momenta obey thermal or uniform distributions (other distributions can also be specified). In the mode of uniform momentum 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. In the uniform momentum distribution mode, if we turn off the transport equation of unbound pairs and only simulate the transport equations of quarkonia, we find that the system does not properly thermalize. It is the transport of open heavy flavors that drives the kinematic thermalization of all heavy quark states. These findings have been reported in Refs. [27,77]. Here we extend the previous studies to the case of excited states. For simplicity, we will simulate all the cases with the open heavy flavor transport equations turned on. The lessons we learn by comparing the case with open heavy flavor transport equations turned on and that off have been discussed before. We focus on demonstrating the consistency of the numerical implementation of dissociation and recombination here.
We initialize N b,tot = 50 bottom quarks and Nb ,tot = 50 antibottom quarks in the QGP volume. Their positions are random and their momenta are sampled in the uniform distribution 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 simulate the Υ(1S) channel for quarkonium, i.e., only the dissociation and recombination of Υ(1S) are allowed.
2. Only the Υ(2S) channel is turned on. Here, the temperature is fixed at 180 MeV. 2 3. Same as case 2, but only χ b (1P ) is studied.
We will compute the hidden bottom fraction as a function of time, which is defined by 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 bottom fraction from excited bottomonia is tiny in thermal equilibrium and one requires large statistics to obtain reasonable results. Furthermore, reaction rates become bigger as temperature increases, so a smaller time step is required.
where N H is the total number of all bottomonium states. For the three cases listed above, only one bottomonium state contributes to N H . We will compare simulation results of the hidden bottom fraction with that in thermal equilibrium, which can be calculated as follows. In thermal equilibrium, the numbers of open bottom quarks and a specific quarkonium state H are given by Here g i is the spin and color degeneracy factor and λ i is the fugacity. We have g b = 6 for the bottom quark, g H = 3 for Υ(nS) and averaged χ b (nP ). Since the total number of bottom quark is equal to that of antibottom, λ b = λb. In detailed balance, we have With this relation, one can solve the fugacities from the balance equation with a given volume: Once we obtain the fugacities, we can compute the hidden bottom fraction in thermal equilibrium: The comparison between our simulation results and the system properties in thermal equilibrium is depicted in Fig. 1 for the three cases listed above. The results are obtained from averaging 10000 simulation events. As can be seen from the plots, the interplay between dissociation and recombination can drive the system to equilibrium. The relaxation rate of the system is about 7 − 8 fm/c for the system conditions used here. This number is on the order of the lifetime of QGP in real collisions. But in general, the relaxation rate depends on the initial p T spectrum and density of the heavy particles, as well as the temperature of the medium. Our simulations can reproduce the correct limit of the hidden bottom fraction in thermal equilibrium. These tests serve as consistency checks in our studies and we are now ready to move on and simulate real collision systems.

Treatment of Relativistic Heavy-Ion Collisions
In this section we will describe our numerical setup to study bottomonium production in heavy-ion collisions. We include Υ(1S), Υ(2S), Υ(3S), χ b (1P ) and χ b (2P ) states in the quarkonium transport equations. To solve the coupled transport equations, we need an initial condition and a bulk QGP medium. For the calculations of R AA , we also need to include all feed-down contributions from excited states to the ground and lower excited states in the hadronic gas stage. We will explain the calculation of the initial phase space distribution, the bulk QGP evolution and the feed-down network in a sequence below.

Initial Conditions
The initial phase space distributions are determined as follows. The momenta of heavy quark-antiquark pairs and each quarkonium state are calculated and sampled utilizing Pythia [78]. The quarkonium production calculation in Pythia is based on the NRQCD factorization [76]. We use the nuclear parton distribution function (nPDF) parametrized by EPPS16 [79]. The nPDF is the only cold nuclear matter (CNM) effect we include. The suppression factors caused solely by the CNM effects R CNM for the Pb-Pb collisions at 5.02 TeV are shown in Fig. 2 for different bottomonia states. We use the same bins in the transverse momentum p T and rapidity y as the CMS measurements, which will be shown later. The uncertainty bands are estimated by the following method: We calculate the CNM effect for each set of the nPDF's in the EPPS16 parametrization, by using the LHAPDF interface [80]. The set 1 corresponds to the central value while the sets 2 − 41 are error sets. (In the LHAPDF interface, the set 0 corresponds to the central value while the sets 1 − 40 are error sets.) The uncertainty is given by (41) of Ref. [79]: where R CNM (i) denotes the CNM effect calculated from the i-th error set.
The production of χ b (2P ) is not available yet in Pythia, so we will assume the CNM effect on χ b (2P ) is the same as that on χ b (1P ). Furthermore, Pythia with nPDF EPPS16 cannot describe the nuclear modification factor R pAu = 0.82 ± 0.10(stat) ± 0.08(syst) measured by the STAR collaboration [81]. So we use R CNM = (R pAu ) 2 = 0.67 as the central value of the CNM effect for the 200 GeV Au-Au collision. We estimate the uncertainty in the R pAu as ∆R pAu = √ 0.11 2 + 0.08 2 ≈ 0.128, and the uncertainty of R CNM in Au-Au collisions as 2 × R pAu × ∆R pAu = 0.21. We assume the R CNM here is p T -independent since the p T range covered by the STAR measurement is limited. The y-dependence of R CNM is not needed here because the STAR measurements are carried out in the mid-rapidity region (|y| < 0.5).
The position distribution of the initial production vertices of heavy quark-antiquark pairs and quarkonia are calculated and sampled using the T R ENTo model [82]. The T R ENTo model assumes the entropy density deposited by the collision at mid rapidity is given by are the nuclear thickness functions of the two projectiles separated by the impact parameter b ⊥ , p is a parameter that has been calibrated on bulk observables of the QGP and τ 0 is the thermalization time of the system after the initial collision. We choose τ 0 = 0.6 fm/c, before which the system is just free streaming. The parametrized entropy density will be used as the initial condition of the hydrodynamic equation, which will be explained in the next subsection. The production of heavy quark-antiquark pairs and quarkonia are thought to be hard processes because of the large mass scale. Therefore, their initial production probability is assumed to be proportional to the binary collision density in the transverse plane, which is proportional to T A T B . Since the nuclear thickness function is used in both the initial entropy density and the initial heavy quark-antiquark production, the corona effect is taken into account in our calculations.
All heavy particles are assumed to be produced at τ = 0 and they propagate via free streaming until τ = τ 0 = 0.6 fm/c. Production at τ = 0 is considered a valid assumption for heavy quark-antiquark pairs because their production time is about 1 M in their rest frame. But this may no longer be true for quarkonia, whose formation time is estimated as 1 M v 2 in their rest frame, which is the time for the heavy quark-antiquark pair to develop          For the ground state, it is probably formed before τ 0 . Since we assume particles are free streaming before τ 0 , it does not matter whether the ground state is free streaming as a fully formed quarkonium state or an unbound pre-resonant heavy quark-antiquark pair. But it matters for excited quarkonium states, because they may be formed inside the thermal QGP. In other words, at τ 0 , we may only have pre-resonant excited quarkonium states. But this is just a tiny effect in our calculations. The excited states have very large dissociation rates when the temperature is high. They will dissociate in one time step after entering the QGP and becoming correlated unbound pair to evolve further. It really does not matter if we have pre-resonant or completely formed excited quarkonium states, because they evolve as unbound QQ pair when entering the thermal QGP. Improvements can be done by including the relative momentum broadening of the QQ pair in the pre-thermalization stage, which will be left to future studies.
For our calculation on bottomonia, we will only initialize quarkonia states but not unbound bb pairs, because the number of such pairs produced in current heavy-ion collision experiments is very small (for central Pb-Pb collisions at 5.02 TeV, the average number of unbound bb pairs is less than one per rapidity). So unlike charmonium production, uncorrelated recombination is negligible for bottomonium production. Since most of the unbound heavy quark-antiquark pairs are produced back-to-back in the transverse plane initially, correlated recombination of these unbound pairs is also negligible. Due to the even smaller production cross section of bottomonium, we will assume in our calculations at most one bottomonium state is produced initially in one heavy-ion collision event.

Medium Evolution
We use a 2 + 1 dimensional viscous hydrodynamic model VISHNew [83,84] to describe the evolution of bulk QGP matter. The package numerically solves with the energy-momentum tensor π µν = 2η∇ µ u ν (4.6) for given initial conditions. Here e and P are the local energy density and pressure, g µν is the metric and u µ is the local four-velocity of the medium. Π is the bulk stress with the bulk viscosity ζ, and π µν is the shear stress tensor with the shear viscosity η. The angle bracket means traceless symmetrization. Both the bulk and shear viscosities are parameters here and can be temperaturedependent. We will use the parametrizations and calibrations in Ref. [85]. Ref. [85] uses the T R ENTo model to calculate the initial entropy density s and then obtain the energy density and pressure with an equation of state calculated by lattice QCD. To obtain the initial stress-energy tensor at τ 0 , Ref. [85] further assumes the initial flow velocity and the viscous terms vanish at τ 0 . All parameters in the T R ENTo model and the hydrodynamic equations are calibrated to experimental observables on light hadrons. For 2.76 TeV Pb-Pb collisions, we will use the parameters calibrated to charged particles shown in Table III of Ref. [85] with the exception of the parameter p in the T R ENTo model. We choose p = 0, which is consistent with the calibration in Ref. [85] (the calibrated value of p is 0.03 +0. 16 −0.17 ). For the parameter calibrations in other collision energies and systems, we will follow Ref. [74].

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. The final nuclear modification factor R AA (i) for i = Υ(1S), Υ(2S), Υ(3S), χ b (1P ) and χ b (2P ) includes both the CNM effect R CNM and the hot medium effect. The evaluation of the CNM effect has been discussed in Sect. 4.1. We will now discuss how we compute the hot medium effect.
Our calculation includes, unlike many others, correlated recombination and we find that 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 2P state produced initially, melts inside the QGP and recombines as a 1S state which survives the following in-medium evolution. To this end, for each centrality bin, we simulate N init i events in each of which one i-quarkonium state is initialized. After the in-medium evolution, among these N init i events, there are N i→j events that have a jquarkonium state in the end. We note in general j N i→j ≤ N init i . For example, N 1S→2S is the number of Υ(2S) states generated from those N init 1S events where the initial quarkonium state is a Υ(1S). N i→i is the "surviving" number whose physical meaning is well-known (it also has contribution from the correlated recombination), and has been calculated in many studies. N i→j with j = i is the contribution from correlated, cross-talk recombination.
As explained in Sect. 4.1, we only initialize one bottomonium state in each relevant collision event. For simplicity of the expressions, we will assume all N init i 's are the same and 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 we consider) is given by (we include CNM effect here) where R CNM (j) is the CNM effect on the primordial production of the j-quarkonium state and σ i is the initial primordial production cross section of the i-quarkonium state (without any feed-down contributions). The ratios of the initial production cross sections are listed in Appendix B. Now we include the feed-down contributions for Υ(1S) and Υ(2S):  For j = 3S, 1P and 2P , we do not consider any feed-down contributions. Then their nuclear modification factors are given by (4.13)

Uncertainty Estimates
We will discuss three uncertainty sources here. The first source is the uncertainty in the EPPS16 parametrizations of the nPDF. The uncertainty bands of the CNM effect caused by the nPDF have been estimated in Section 4.1 and will be included in the plots to be shown in the next subsection. The second uncertainty originates in the values of α s = 0.3 and α pot s = 0.36. We will vary these two coupling constants by ±10% to estimate the uncertainty caused by the parameter values. In the next subsection, we will show three curves for the R AA We expect the first two uncertainty sources will dominate over the last one, so we exclude the last uncertainty source in the following analysis.

Results for Υ(nS)
We first show the results of R AA as a function of centrality at 5.02 TeV Pb-Pb collisions. The comparison between our calculation and the measurements by the CMS collaboration is shown in Fig. 3. The three curves in the plot correspond to the three sets of parameters, as discussed in the previous subsection. The uncertainty band is solely from the uncertainty of the nPDF and is centered at the middle curve, which corresponds to the central values of the parameters α s = 0.3 and α pot s = 0.36. As can be seen in the plot, the nPDF uncertainty dominates over the uncertainty of the parameters. To explicitly demonstrate the importance of correlated recombination, we show the full calculation results in Fig. 3a and the results without any contribution from the cross-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 = i) is excluded. Remaining same-species (i → i) correlated recombination is still included in the latter case, though we expect its contribution to be much smaller than the cross-talk recombination. As can be seen from the comparison of the two figures, the cross-talk recombination is crucial to describe the data, even if one  takes into account the uncertainties from the nPDF and parameter values. Furthermore, it seems unlikely to describe the data for 1S and 2S simultaneously without including the cross-talk recombination, by further increasing the coupling constant in the potential, since the change of R AA (2S) from the parameter variation is tiny. For the 1S state production, most excited states dissociate quickly when entering the QGP because of the initial high temperature. These melted correlated bb pairs may form 1S state after the dissociation because the 1S state can exist at high temperature. For the 2S production, most of the primordially produced 2S states cannot survive the in-medium evolution in most centrality bins. They are regenerated when the temperature cools down and allows their existence in the medium. At low temperature, a 2S state can be formed from a correlated unbound bb pair (which may come from a dissociated quarkonium state such as a 1S state). For the 3S production, since we assume they cannot exist inside QGP (melting temperature is below 154 MeV), correlated recombination does not affect its production. The 3S state is only produced in very peripheral collisions due to the corona effect.
Next we discuss the results of R AA as a function of the transverse momentum, shown in Fig. 4a. First we notice that the experimental results of R AA are almost flat as a function of p T . This is a highly non-trivial result because the CNM effect has a dramatic dependence on p T , as shown by the central values in Fig. 2, though the uncertainty band is quite large. The p T variation of the CNM effect can be understood as follows: At midrapidity, the energy fraction carried by the two gluons that fuse to produce quarkonium (in the collinear factorization, for a non-zero p T , a recoil particle is produced back-toback to the quarkonium state in the transverse plane), is given by Bottomonium mass is about M H ≈ 10 GeV. So at 5.02 TeV Pb-Pb collisions, x goes from 0.004 to 0.009 as p T increases from 0 to 20 GeV. As x increases, the nuclear modification factor R g (x) of the gluon PDF increases. The CNM effect on the production cross section grows as [R g (x)] 2 . This growing trend is true for every fitting set in the  EPPS16. Furthermore, the scale m T used in the perturbative calculation also increases as p T increases. At higher energy scales, the nPDF modification factor R g is expected to be closer to unity, which further increases its p T dependence.
The final flat R AA means for bottomonium, the p T -dependent CNM effect is washed out by the hot medium effect. To show this more explicitly, we plot the calculation results of R AA without the CNM effects in Fig. 4b. We see the hot medium effect leads to a raise of R AA at low transverse momentum. The reason behind is related with correlated recombination: Suppose a quarkonium state dissociates and then recombines later (via correlated recombination). After dissociation but before recombination, the correlated QQ pair interacts with the medium and loses energy. When the pair recombines, it has a smaller p T than that right after the dissociation. This leads to an increase of quarkonium states at low p T . Furthermore, low p T quarkonium states stay inside the medium for a longer time and probably go through the low temperature part of the QGP expansion. Since recombination at high temperature is ineffective, quarkonium states at low p T have a higher chance to recombine (via correlated recombination) if they dissociate than those at high p T . To demonstrate this point more clearly, we plot the calculation results without any CNM effects and without any cross-talk recombination in Fig. 4c, where the suppression has only a mild p T dependence.
Our calculations can describe the R AA for p T < 10 GeV and start to deviate from the data as p T becomes bigger. This may indicate the breakdown of the nonrelativistic expansion that our calculations are heavily replied on. For a finite p T , the typical energy of medium light particles in the rest frame of quarkonium is given by where v T is the transverse velocity of the quarkonium with respect to the local medium. Our calculations are valid if rT √ M v is the typical size of quarkonium. Under our assumed hierarchy rT 1, the condition rT √ Higher order contributions in the nonrelativistic expansion become important and have to be included consistently as p T increases 4 . Our calculations probably indicate that this happens when p T > 10 GeV. One should note that for charmonium this can happen at a much smaller p T because what matters in the argument is the transverse velocity rather than the transverse momentum. Since most bottomonia are produced at low transverse momentum, the deviation at mid and high p T does not affect the centrality and rapidity dependence of R AA in our calculations. Then we show the results of R AA as a function of the rapidity in Fig. 4d. Our calculations are consistent with the experimental measurements. Since we use a 2 + 1 dimensional viscous hydrodynamics, the hot medium effect is independent of y. Furthermore, we see in Fig. 2 that the y dependence of the CNM effect is also mild. Therefore, the final R AA is almost flat in the rapidity.
The comparison between our calculations and experimental measurements for 2.76 TeV Pb-Pb collisions and 200 GeV Au-Au collisions is shown in Fig. 5. Our calculations can describe most of the experimental data except R AA (2S + 3S) in peripheral Au-Au collisions at 200 GeV. In Fig. 5b, our calculations of R AA (2S + 3S) are consistent with the 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 peripheral collision. The uncertainty of the experimental measurement there is quite large due to the 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 peripheral collisions leads to the discrepancy in the p T dependence of R AA (2S + 3S), as depicted in Fig. 5d. We expect that the future sPHENIX collaboration will provide data with high precision and then the origin of the current discrepancy issue will be more clear: Either the calculations are consistent with improved measurements or the calculations need improving. 4 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 and thus the reaction rate grows with the size. But physically we know there is an upper limit of the rate when the size is large: The rate cannot be bigger than twice the reaction rate of one heavy quark interacting with the medium. When the quarkonium size is large enough, the heavy quark and antiquark interact with the medium almost independently. Therefore, we expect that including higher order terms in the nonrelativistic expansion will reduce the rate and thus the RAA at large pT will go up.   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.
Finally, we study the azimuthal angular anisotropy of Υ(1S). In particular, we compute the v 2 coefficient which is defined by in which we neglect higher order harmonics. Here φ − Ψ RP is the azimuthal angle of the 1S state with respect to the reaction plane, which is defined event-by-event. Our calculation result for the 5.02 Pb-Pb collision in the centrality range 5 10 − 90% is shown in Fig. 6. We only show the result from the central values of the parameters. We also plot the experimental results measured by the ALICE collaboration which are in a different rapidity and centrality regions, just to show the state-of-art of the current bottomonium v 2 measurements. We stop our calculations at p T = 24 GeV because as we have seen in the R AA comparison, higher order corrections in the nonrelativistic expansion start to become important as p T increases, which are neglected in our current setup. We can calculate v 2 at higher p T values in our current setup but their physical meaning is less robust and we cannot learn much from doing that. Our calculation result is consistent with the current 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 nonrelativistic expansion 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 mechanism would mainly be jet energy loss. Due to the steep p T spectrum in the primordial production of bottomonium, most of the contribution to the last data point comes from p T ∼ 10 − 20 GeV. So the comparison we show still has some physical meaning for the last data point. Several physical processes contribute to the development of the quarkonium v 2 . The first contribution is the path dependence. In peripheral collisions, the QGP has a elliptic shape. Quarkonia moving along the longer axis will be more suppressed. Also, the reaction rates of quarkonium in the medium depend on the relative velocity between the quarkonium and the local medium, which has a flow velocity. This also influences the shape of v 2 as a function of p T . Finally, after the dissociation of quarkonium, the unbound QQ pair can develop some v 2 by interacting with the medium. Later if they recombine, the v 2 will be partly or fully inherited by the regenerated quarkonium. Uncorrelated recombination can also contribute to v 2 and is crucially important for charmonium production in heavyion collision. Open charm quarks that develop v 2 during the in-medium evolution will contribute to the charmonium v 2 if they recombine. But this v 2 generation mechanism is negligible for bottomonium since the number of unbound bb pairs produced in one collision event is smaller than one per rapidity in mid-rapidity. Future precise measurements on the azimuthal angular anisotropy will greatly help us understand the in-medium dynamics of quarkonium, especially how quarkonia with finite transverse momenta interact with an expanding medium. These non-equilibrium transport properties of quarkonium are not easy to study via finite temperature lattice QCD calculations. The interplay between theory and phenomenology will help deepen our understanding on these, in particular once experimental data with high precision are available.

Prediction: χ b (1P ) More Suppressed than Υ(2S)
One consequence of the important contribution from correlated recombination in bottomonium production would be that the χ b (1P ) state is more suppressed than the Υ(2S) state. If there were no correlated recombination, one would expect R AA (χ b (1P )) to be similar to R AA (Υ(2S)), since the two states have similar binding energies and sizes. However, correlated recombination can alter dramatically this naive expectation. Since χ b (1P ) and Υ(2S) have similar binding energies and sizes, their recombination rates from a correlated bb pair are also close. In particular, the probability of an initial χ b (1P ) state ending up as a Υ(2S) state (via first dissociation and then correlated recombination) is similar to that of an initial Υ(2S) state ending up as a χ b (1P ) state. But the primordial production cross section of χ b (1P ) is 4 − 5 times that of Υ(2S). Much more χ b (1P ) states are produced initially and thus, the number of Υ(2S) states regenerated from initial χ b (1P ) states is much larger than the number of χ b (1P ) states regenerated from initial Υ(2S) states. Therefore, Υ(2S) is less suppressed than χ b (1P ).
To demonstrate this idea, we plot R AA (χ b (1P )) and R AA (Υ(2S)) in Fig. 7 6 for two cases: with cross-talk recombination and without it. It can be clearly seen that in the former case, Υ(2S) is less suppressed while in the latter case, Υ(2S) and χ b (1P ) have  similar R AA 's. So we propose a new measurement to test experimentally the importance of correlated recombination, which is the ratio of R AA (χ b (1P )) and R AA (Υ(2S)). We plot this ratio as a function of transverse momentum in Fig. 7c for both cases. The difference is dramatic: in the case without cross-talk recombination, the double ratio is approximately unity, consistent with the argument based on the quarkonium binding energy and size given above, while in the case with cross-talk recombination, the double ratio can be as small as 0.2 at low p T and gradually increases at higher p T . The uncertainties due to the choice of parameters are relatively small in both cases, as can be seen from the variation of the three curves. In the case without cross-talk recombination, the nPDF uncertainty does not cancel out in the ratio, and is very large at low p T , which reflects the difference of the nPDF effects on 2S and 1P states as shown in Fig. 2. In the case with cross-talk recombination, the nPDF uncertainty band becomes much narrower. The reason of the narrowness is twofold. First, the CNM effect and the hot medium effect are multiplicative (see Eq. (4.7)). When the central value of the ratio is smaller, as in the case with cross-talk recombination, the nPDF uncertainty band would look narrower in a plot with a linear scale in the vertical axis. Furthermore, because of cross-talk recombination, different bottomonium states can turn to each other during the evolution and thus the difference of the CNM effects in the initial production is gradually washed out. In a nutshell, after accounting for the calculation uncertainties, the suppression of χ b (1P ) relative to Υ(2S) is still manifest, especially at low and intermediate p T . As can be seen from Fig. 7c, this double ratio observable is of huge model discriminatory power. However, measurements on the double ratio may be challenging at the moment. Reconstructing the χ b (1P ) state requires detecting the low energy photon emitted in its decay to the Υ(1S) state. But a large number of photons at low p T are produced in heavy-ion collisions due to the neutral pion decay, which leads to a large combinatorial background in the χ b (1P ) reconstruction. Furthermore, the photon resolution in the calorimeter is limited. To overcome the resolution problem, one can reconstruct the photon from its conversion to an electron positron pair due to the interaction with detector materials, but the conversion reduces the signal by a substantial factor of order 10 [90]. Moving to intermediate and high p T would help detecting the emitted photon, but at the same time the interesting signal in the double ratio dies away as p T increases, needless to say the available experimental statistics at high p T . A tradeoff must be made between the interesting physics, detector resolution and experimental statistics. Considering these factors and the large nPDF uncertainty at low p T in the case without cross-talk recombination, an optimal p T window for such a double ratio measurement might be roughly between 5 and 15 GeV. Detailed experimental analysis on the p T searching window is needed. With higher luminosity and improved detector efficiencies, Run 3 of LHC and the sPHENIX program at RHIC may be able to measure this double ratio observable.

Conclusions
In this paper, we develop a set of coupled transport equations for open heavy quarkantiquark pairs and quarkonia. Our framework is capable of calculating observables for both open and hidden heavy flavors. Dissociation and recombination rates are calculated from pNRQCD via a systematic weak-coupling and nonrelativistic expansion, with the assumption of a weakly-coupled QGP. Recombination rates depend on the real-time heavy quark-antiquark distribution function, which is solved by the transport equation for open heavy flavors. Our framework can handle both uncorrelated and correlated recombination. By numerically solving the coupled transport equations via Monte Carlo techniques, we demonstrate how the ground and excited bottomonia states approach thermal equilibrium inside a QGP box with a constant temperature. Furthermore, we study bottomonia production in heavy-ion collisions using our framework. The initial phase space distributions are calculated from Pythia with the nPDF EPPS16 and the T R ENTo model. The medium description utilizes a 2 + 1 dimensional viscous hydrodynamic model. The parameters of the T R ENTo model and the hydrodynamic equations have been calibrated to experimental observables on light hadrons. We include Υ(1S), Υ(2S), Υ(3S), χ b (1P ) and χ b (2P ) states in our reaction network. Our calculations demonstrate the importance of correlated cross-talk recombination in bottomonium phenomenology and can describe most of the experimental data. Discrepancies are seen at mid and large transverse momenta, where omitted higher order terms in the nonrelativistic expansion become gradually more important. We propose a new measurement on the ratio of nuclear modification factors of χ b (1P ) and Υ(2S) to test the importance of correlated recombination in experiments.
The current formalism can be improved in several ways. First, one can include the effect of the running coupling constant and higher order corrections in the nonrelativistic expansion. One may also consider extending the calculation of the reaction rates to the case 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 potential model. The nonperturbative potential may be parametrized and the parameters can be calibrated to experimental observables such as R AA and v 2 , for example, by using the recently developed Bayesian analysis techniques [91]. A simultaneous description of both open and hidden heavy flavor observables is also possible in our framework and is worth exploring. With a nonperturbative potential, the technical challenge would be to develop a fast numerical algorithm to sample the momenta of the outgoing QQ pair (in dissociation) or quarkonium (in recombination) from the differential reaction rates. No previous studies have done this, but it is crucial for a consistent microscopic treatment like our calculation here. With a Coulomb potential, reaction rates have analytic expressions which are easier to handle in inverse transform and importance samplings (For the case of how to sample efficiently with a Coulomb potential, see Chapter 4.2 of Ref. [92]). For a nonperturbative parametrization of the screened potential, one has to include its dependence on the relative velocity between the quarkonium state and the hydro-cell. We use Coulomb potential here so we do not need to consider the velocity dependence in the potential. Also, the CNM effect at 200 Au-Au collisions should be better understood and constrained by the measurements in p-Au collisions. Furthermore, the transport in the pre-thermalization stage should be investigated. For open heavy flavors, this has been done in Ref. [93]. If one includes the transport of heavy flavors in the pre-thermalization stage, for consistency, one should also take into account the thermalization process of the medium between the initial hard collision and the hydrodynamics starting time. Parameter calibrations should be reanalyzed with free streaming replaced by a more realistic thermalization process. Finally, we will extend the current study of bottomonium production to charmonium production in heavy-ion collisions. There, we must also initialize unbound cc pairs, since many of them are produced primordially and uncorrelated recombination is enhanced. DOE 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. [75]. We will use k to label the momentum of the quarkonium involved in a scattering process. q will be used to indicate the momentum of the gluon absorbed or emitted by the quarkonium. This gluon can be real in the 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 1 and p 2 are used indicate the momenta of the light quarks on the left and right respectively. Similarly for the process g + H ↔ g + Q +Q, we will use q 1 and q 2 to represent the momenta of the gluons on the left and right respectively. In the quarkonium rest frame, the energy of the quarkonium is given by −|E nl | where E nl is the binding energy. In the rest frame of the unbound QQ pair, its energy is given by p 2 rel M where p rel is their relative momentum. The expressions shown below are valid in the rest frame of quarkonium for dissociation or the center-of-mass frame of a QQ pair for recombination. These two frames are not equivalent, but their difference is negligible (suppressed by T M ). In real simulations, the gluon distribution is boosted from the local rest frame of the hydro-cell (where temperature is defined) to these rest frames of the heavy particles. After calculating the reaction rates in the rest frames of the heavy particles, one has to boost them back to the laboratory frame where we keep track of the phase space distributions. In the following, for simplicity, we will only show expressions for quarkonia or QQ pairs whose center-of-mass motions are at rest with respect to the medium local rest frame. In practice, we account for the Lorentz boost properly, as described above.
The dissociation and recombination terms in the transport equation of the quarkonium state nls can be written as where we used the "δ−derivative" symbol defined above. The functions F ± nls are defined by where n B is the Bose-Einstein distribution function, g + = 1 N 2 c g s and g s is the degeneracy factor for spin: g s = 3 4 for a quarkonium state with spin s = 1 and 1 4 for spin s = 0. For an arbitrary QQ pair, its probability of being a spin-1 state is 3 4 . Its probability of being a color octet is N 2 . When computing the scattering amplitude squared, we have summed over the relevant quantum numbers (color of the QQ, gluon polarization and color), so we only use 1 N 2 c to avoid double counting. In the definition of F − nls , p cm = p Q + pQ, and p rel = p Q −pQ 2 are the center-of-mass and relative momenta of a QQ pair with momenta p Q and pQ.
After summing over the relevant quantum numbers (see Ref. [75] for details), the scattering amplitude squared is given by 2Nc , |ψ nl is the wavefunction of the quarkonium state nls (states with different spins are degenerate) and |Ψ p rel is the Coulomb scattering wave for the unbound QQ pair. For non-S wave states, we average over the polarizations: 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) = We do not need the dipole transition matrix elements for 3S and 2P states, since we assume they cannot be formed inside the QGP. In other words, if a 3S or 2P state enters the 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.

B Details on Feed-Down Contributions
We will use the notations σ tot nl and σ nl to denote the total and primordial cross sections in a 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 lowerexcited states, the first thing we need is the branching ratio in vacuum. The most recent results reported by the Particle Data Group [95] are summarized in Table 1. For the P-wave states, we need the averaged branching ratio since our transport equations are degenerate  σχ b1 (1P ) in Ref. [96] and follow the assumptions made in Ref. [24] to write The next thing we need is the primordial cross section ratio. Using the experimental Cross sections in proton-proton collisions Experimental results (nb) from Ref. [86] B × σ(Υ(1S)), |y| < 2.4, 5.02 TeV 3.353±0.081(stat)±0.167(syst) B × σ(Υ(2S)), |y| < 2.4, 5.02 TeV 0.873±0.031(stat)±0.046(syst) B × σ(Υ(3S)), |y| < 2.4, 5.02 TeV 0.404±0.017(stat)±0.022(syst)  inputs listed in Table 2 and following Ref. [24] and references therein, we assume Experimentally it is known that the feed-down contributions to σ tot 1S is about 67% (averaged over the transverse momentum), consistent with our prescription here. Now we are ready to compute the ratios of the primordial cross sections. The results are listed in Table 3.