Approach to equilibrium of quarkonium in quark-gluon plasma

We calculate the dissociation and recombination rates of $\Upsilon(1S)$ in quark-gluon plasma by using potential non-relativistic QCD. We then study the dynamical in-medium evolution of the $b\bar{b}-\Upsilon$ system in a periodic box via the Boltzmann equation and explore how the system reaches equilibrium. We find that interactions between the free heavy quarks and the medium are necessary for the system to reach equilibrium. We find that the angular distribution of $\Upsilon(1S)$ probes the stages at which recombination occurs. Finally, we study the system under a longitudinal expansion and show that different initial conditions evolve to distinct final ratios of hidden and open $b$ flavors. We argue that experimental measurements of the ratio could address open questions in the quarkonium production in heavy ion collisions.

Heavy quarkonium is used as a probe of quark-gluon plasma (QGP) produced in relativistic heavy ion collisions. In their pioneering work, Matsui and Satz [1] argued that Debye screening of the color attraction between heavy quarks (Q) and anti-quarks (Q) leads to the dissociation of the bound state at sufficiently high temperature. As a result, the quarkonium yield was predicted to be suppressed with respect to the scaled yield measured in proton-proton (pp) collision. This ratio is defined as the nuclear modification factor R AA .
The static screening mechanism proposed in the original work is obscured by several factors: nuclear modification of initial production, dynamical screening (dissociation caused by scattering), recombination and feed-down contribution after hadronization. Recent progress in understanding the full complexity of quarkonium suppression includes the calculation of the imaginary potential of QQ, which is related to the dissociation rate [2,3]; the calculation of the viscous (anisotropic) corrections to the real and imaginary parts of the potential [4][5][6]; the calculation of the dissociation rate in potential nonrelativistic QCD (pNRQCD) [7][8][9] and holographic gravity models of QCD [10,11]; the description of the time evolution in the open quantum system approach [12][13][14][15][16][17][18][19]; the use of Langevin equations to study the impact of QQ diffusion [20] on quarkonia recombinations [21][22][23][24][25]; and phenomenological studies that can reproduce the experimentally measured R AA of bottomonia with [26] and without recombination [27].
The process of recombination is less well-understood than dissociation. Different studies incorporate recombination in varying ways. In the open quantum system formalism, an initially unbound QQ pair will evolve to have non-vanishing overlap with bound state wave functions due to the action of stochastic potentials. Onedimensional numerical studies in the abelian case have been carried out [14,16,18,19] and the non-abelian case has been developed formally [15]. Other approaches use coalescence models based on Wigner functions [21,25] * xiaojun.yao@duke.edu or invoke detailed balance to model the recombination rate as the dissociation rate times equilibrium quarkonium fraction [26], which is only true near chemical and thermal equilibrium. Here we aim to present a calculation that incorporates dissociation and recombination in a consistent way and without assuming the quarkonium distribution is close to equilibrium.
Recombination is believed to be less important for bottomonium than charmonium because fewer b-quarks are produced in the collision. This would be consistent with experimental measurements if the dissociation is the dominant in-medium process, which is based on the assumption of small suppression of initial quarkonium production. However, we only know the nuclear modification on parton distribution functions is small [28], which is not sufficient for the assumption. Quarkonium production in pp collisions factorizes into short-distance production of heavy quarks and long-distance coalescence into quarkonium [29]. It is unlikely that the long-distance physics completes before the formation of QGP in heavy ion collisions. Therefore the fraction of quarkonia formed in the initial stage is not well-determined.
Furthermore, if the initial QGP temperature is higher than the melting temperature above which certain quarkonium state cannot exist due to static screening (which has been studied from the temperature dependence of the binding energy [30] or spectral functions [31] in potential models), correlated QQ pairs will remain unbound when entering QGP and may (re)combine later. This type of process is also defined as recombination throughout this paper. In short, recombination of correlated and uncorrelated heavy quark anti-quark pairs may be more important than originally thought.
To illustrate these points, we propose a dynamical inmedium transport model based on Boltzmann evolution with dissociation and recombination. This approach allows us to explore how experiments can help address various unknown aspects of quarkonium production mechanisms in heavy ion collisions. Here we consider the Υ(1S) as an example. The physical picture is as follows: if the local temperature is higher than the melting temperature, Υ(1S) dissociates and locally only unbound bb exist. However, if the local temperature is lower, on one hand, Υ(1S) can exist and propagate in the medium and may be dissociated by scattering with medium gluons and light quarks; on the other hand, unbound b andb propagate and diffuse, and at any time, they may recombine into Υ(1S) by scattering with medium constitutes, if they are sufficiently close to each other and their relative momentum favors recombination.
We calculate the dissociation and recombination rates to lowest order in pNRQCD [32,33]. The effective theory can be derived from QCD under the hierarchy of scales M M v M v 2 , T, m D where M = 4.65 GeV is the b-quark mass, v ∼ 0.3 is the relative velocity of bb inside Υ(1S), T is the temperature, and m D is the Debye screening mass. The pNRQCD Lagrangian is given by where E represents the color electric gauge field. The Lagrangian of gluon and light quark is just QCD with momenta M v. The pNRQCD is a systematic expansion in v or 1/M (NR expansion) and r, the relative distance between bb (multipole expansion). The degrees of freedom are the color singlet S(R, r, t) and color octet O(R, r, t) states where R denotes the center-of-mass (c.m.) position and r the relative coordinate. The color singlet and octet Hamiltonians are expanded in powers of 1/M : By the virial theorem, Higherorder terms of potentials including relativistic corrections, spin-orbital and spin-spin interactions are further suppressed by extra powers of v. The c.m. kinetic energy is also suppressed because momenta ∼ M v have been integrated out in the construction so P cm M v. We only work to order M v 2 because v is small and also heavy ion experiments do not resolve hyperfine structures currently. In the following, we only consider temperatures at which the Υ(1S) exists (T C < T < 2.5T C , T C = 155 MeV). In this domain the confining potential is flattened and potentials can be approximated by Coulomb interactions The singlet-octet and octet-octet vertices are color dipole The Feynman diagram of the transition between the singlet Υ(1S) and the unbound bb octet via absorption or emission of a gluon is shown in Fig. 1. For simplicity, we here consider only the interaction with on-shell gluons in the QGP. Transitions mediated by virtual gluons (inelastic scattering with medium constitutes) are at next order in α s and can be easily included within our formalism. But we are also aware that when m D E 1S , the inelastic scattering dominates [9]. Our future full calculations will include both. The scattering amplitude is given by where T F = 1/2, |ψ 1S is the hydrogen-like 1S wave function for Υ(1S) and |Ψ p rel is the Coulomb wave function for unbound bb octet. The 1S binding energy is given by E 1S = −α 2 s C 2 F M/4 and the gluon energy is q = |q|. Throughout our paper we set α s = 0.3. Here k 1,2 are c.m. momenta and their associated kinetic energies will be neglected according to the power counting.
The set of Boltzmann equations for the b,b and Υ(1S) distribution functions f i (x, p, t) is given by For b andb quarks, the collision term C i describes their scattering with thermal constituents of QGP. This process has been described either as diffusion in the framework of the Langevin equation [34,35] or as two-body scattering in the framework of the linearized Boltzmann equation [36][37][38]. Here we use, for simplicity, the relaxation-time approximation with The relaxation rate is assumed to be Γ r = T 2 /M [34] and f eq i is the relativistic Boltzmann distribution. For Υ(1S), the gain term C + is from recombination by gluon emission and the loss term C − is from dissociation by gluon absorption 1 : where the second equation defines the gluo-dissociation rate Γ d . The scattering amplitude is calculated in the rest frame of Υ for dissociation and that of bb for recombination, where the pNRQCD is valid. The Bose distribution of medium gluons n (q) B is boosted into the rest frame of Υ or bb accordingly. The two frames are not equivalent but since the gluon energy is small compared to M (T M ), the difference is suppressed by T /M . The overline indicates an average over initial-state and sum over final-state quantum numbers (color and spin). The phase space measure is relativistic for gluons and non-relativistic for b-quarks, which is consistent with our field definitions.
For dissociation, the rest-frame rate is then boosted back into the medium frame by the factor γ −1 = √ 1 − v 2 where v is the Υ velocity. For recombination, the gamma factor cancels out, as explained further below. The quark momenta p 1 and p 2 are related to the relative momentum p rel used in the amplitude calculation via 1 2 (p 1 − p 2 ) where the primed momenta are in the bb rest frame. The factor 8/9 ensures that only a color octet bb pair can form a singlet bound state by emitting a gluon.
We solve the Boltzmann equations (5) by stochastic simulation. Here we study the evolution inside a periodic box of QGP with side length L = 10 fm. The QGP tem-perature is constant throughout the box but can change with time. A certain number of b,b and Υ (N b , Nb and N Υ ) are initialized by random sampling of their positions and momenta, assuming a given initial momentum distribution. At each time step ∆t, we consider three types of processes: First, for each Υ with a given velocity, we determine whether it dissociates according to the probability Γ d ∆t.
If it dissociates, we sample the incoming gluon momentum in the rest frame of Υ according to the integrand of the rate integral and then calculate the outgoing relative momentum of bb by energy-momentum conservation. Finally we boost the momenta of b andb back into the medium frame; their positions are set to be the position of Υ before dissociation.
For each b-quark with position y i and momentum k i we need to determine the total recombination rate with neighboringb quarks with position z j and momenta k j . However, the quark and anti-quark distributions in the expression (6) should be evaluated at the same position, but the product of two delta functions is ill-defined. We introduce a position-dependence of the recombination probability by means of a Gaussian function with a width σ chosen to be the Υ(1S) Bohr radius. This ensures that the recombination probability of a widely separated bb pair vanishes. The product of local distributions in (6) is thus replaced with where the sum over i, j runs over all b,b contained in the box. The position-dependence (including the choice of σ) disappears if one averages over many spatial configurations. The theta function assures that only approaching quark pairs can recombine. Thus, a bb that has just dissociated cannot recombine until at least one of them scatters once. For a given b-quark, theb density is Lorentz boosted into the rest frame of bb, and when the rate is transformed back into the medium frame, the two γ factors cancel. This explains the absence of gamma factor in (6). If the b-quark is found to recombine, ab-quark is chosen based on the value of recombination probability. We then sample the outgoing gluon and replace the bb pair with a Υ whose momentum is determined from energy-momentum conservation. Its position is given by the c.m. position of the quark pair as indicated in (8).
Third, the diffusion of unbound b orb quark in the QGP is implemented by re-sampling its momentum from thermal distribution at a probability of Γ r ∆t in each time step. We exclude elastic scattering between medium particles and bound Υ because it cannot happen at the order we are working.
As a first application, we study how the bb-Υ system reaches chemical equilibrium. We set up the system at a constant temperature T = 350 MeV with N b,tot = N b + N Υ = 50 in three different initial conditions: 3. as case 2, but including HQ diffusion.
The results of N Υ /N b,tot are plotted in Fig. 2. At equilibrium  and solved from N eq b + N eq Υ = N b,tot . The simulations converge to the NR lines because the rates are calculated in a pNRQCD. If excited states are included, the Υ(1S) equilibrium fraction will decrease but only insignificantly. As the lower part of Fig. 2 shows, HQ diffusion is necessary for the system to reach equilibrium starting from a non-thermal initial distribution.
We next study the azimuthal angular anisotropy of Υ produced from recombinations of b andb with certain azimuthal momentum distributions simulating elliptic flow of the QGP, which is gradually transmitted to the unbound heavy quarks by diffusion during the QGP phase. Since quarkonia can form at any time below the melting temperature and not necessarily have to wait until the QGP hadronizes [40], measurements of the quarkonium elliptic flow can, in principle, tell us at what time quarkonia are formed by recombination. Therefore it is important to understand how the elliptic flow transmits from heavy quarks to quarkonia. In our study, the momentum distributions of b andb are chosen as: where φ is the angle around the z-axis. The initial p T distribution is taken from the FONLL calculation for 2.76 TeV Pb-Pb collision at rapidity y = 0 [39]. Pairs of b and b are sampled and recombined by gluon emission according to the rate at T = 250 MeV assuming they are at the same position. The v 2 of produced Υ is computed by averaging cos(2φ) in each p T bin with size 1 GeV. The results are plotted in Fig. 3. At low p T , the distribution becomes isotropic as expected. As p T increases, the curves are flatten out. We note that at high p T fragmentation becomes the dominant mechanism, which will be studied in future work. In the plotted p T range where recombination dominates, the quarkonium v 2 is sensitive  to that of heavy quarks. Finally we study the dynamics of the system under a boost invariant longitudinal expansion [41]. The temperature dependence is the Bjorken model given by Here we assume t 0 = 1 fm/c, T 0 = 350 MeV and a speed of sound c 2 s = 1/3. The various rates are plotted as a function of temperature in Fig. 4: the expansion rate defined as |dT /dt|/T , the dissociation rate of a static Υ, the thermally averaged recombination rate, and the HQ relaxation (thermalization) rate. We simulate the system starting either at N b = 5, N Υ = 0 or at N b = 0, N Υ = 5 with HQ diffusion. The initial momenta of b andb are randomly sampled angularly with the magnitude distributed according to the p T spectrum in the same FONLL calculation as above. The momentum distribution of Υ is given by the convolution of those of b andb. The evolution of the Υ(1S) fraction is shown in Fig. 5. The fraction in pp collision is roughly 1.76 × 10 −3 [26] and is indicated by the dotted horizontal line. For comparison, we take the weighted average of the two simulations with initial fraction 0 or 1 so that the initial fraction starts at the pp value.
For the curve starting at N Υ /N b,tot = 1 (all b,b initially bound) dissociation is the dominant process. Because the dissociation and expansion rates are on the same order, as shown in Fig. 4, the survival probability of Υ is large and the curve stays far away from equilibrium at the end of expansion.
On the other hand, the curve starting at N Υ /N b,tot = 0 (all b,b initially free) always fall below the equilibrium. The reason is two-fold: The recombination is significantly slower than the expansion as shown in Fig. 4 and the thermalization of HQ is not fast enough. We also studied simulations without HQ diffusion and find that the influence of HQ diffusion is small in this scenario. However, if all the rates except the expansion rate were larger, the curve including HQ diffusion would be closer to the equilibrium curve, though it would still lag behind.
Lastly, the curves starting at the pp fraction happen to approach the equilibrium line in the end, but this does not indicate the system reaches equilibrium. The recombination contribution here is negligible. However, we note that both the equilibrium fraction and recombination contribution depend on the value of N b,tot . For much larger values of N b,tot the equilibrium and recombination curves would move up and the recombination contribution would be significant, similar to what is observed in the charm sector.
It can be seen that different initial conditions lead to largely distinct final ratios. Since we do not fully understand the initial production of quarkonia in heavy ion collisions, we can hope to learn this together with the in-medium evolution from experiments. The change during the hadronic phase should be small due to the small cross sections [42]. Therefore, it is important to measure the final ratios of hidden and open heavy flavors in various p T and rapidity ranges as a function of centrality and collision energy. We will gain more information from these measurements on the quarkonium production mechanism. Thus, it is essential to do the measurements at both the CERN Large Hadron Collider and BNL Relativistic Heavy Ion Collider.
In summary, we have used pNRQCD to calculate the dissociation and recombination rates of Υ(1S) in a thermal QGP. We studied the dynamics of the bb − Υ system in QGP via the Boltzmann transport equation, which we solved by Monte-Carlo simulation. We showed how the system reaches equilibrium starting from different initial conditions. We demonstrated the importance of HQ diffusion in the medium: It is necessary for the system to reach equilibrium. We then calculated the elliptic flow of Υ produced from recombinations. We argued that measurements of v 2 (Υ) probe stages of quarkonia production by recombination. Finally we studied the system under a Bjorken expansion. We showed that different initial Υ fractions evolve to widely different final results and argued that measurements on the hidden-to-open heavy flavor ratio could address open questions in quarkonium production in heavy ion collisions.
We acknowledge helpful discussions with Steffen Bass, Weiyao Ke, Michael Strickland and Yingru Xu. XY thanks Nora Brambilla, Miguel Escobedo, Jacopo Ghiglieri, Péter Petreczky and Antonio Vairo for discussions on pNRQCD and acknowledges the hospitality of the nuclear theory group at Brookhaven National Laboratory where part of this work was completed. XY acknowledges support from U.S. Department of Energy (Research Grant No. DE-FG02-05ER41367) and Brookhaven National Laboratory.