Quarkonium inside Quark-gluon Plasma: Diffusion, Dissociation, Recombination and Energy Loss

We consider the quarkonium diffusion, dissociation and recombination inside a quark-gluon plasma. We compute scattering amplitudes in potential nonrelativistic QCD for relevant processes. These processes include the gluon absorption/emission at the order $gr$, inelastic scattering at the order $g^2r$ and elastic scattering with medium constituents at the order $g^2r^2$. We show these amplitudes satisfy the Ward identity. We also consider one-loop corrections. The dipole interaction between the color singlet and octet is not running at the one-loop level. Interference between the tree-level gluon absorption/emission and its thermal loop corrections cancels the collinear divergence in the $t$-channel inelastic scattering. The inelastic scattering has no soft divergence because of the finite binding energy of quarkonium. We write out the diffusion, dissociation and recombination terms explicitly for a Boltzmann transport equation and define the dissociation and recombination rates. Furthermore, we calculate the diffusion coefficient of quarkonium. We find our result of diffusion coefficient differs from a previous calculation by two to three orders of magnitude. We explain this and can reproduce the previous result in a certain limit. Finally we discuss two mechanisms of quarkonium energy loss inside quark-gluon plasma.


I. INTRODUCTION
Since the early study of static screening on quarkonium [1], the bound state of a heavy quark antiquark pair, quarkonium has been used as a probe of quark-gluon plasma (QGP) in heavy ion collisions. In a high temperature QGP, the attractive potential between the heavy quark antiquark pair QQ is significantly suppressed and thus the bound state cannot exist. The screening effects have been widely investigated by computing the free energy of QQ [2,3] or spectral functions [4,5] on a lattice at finite temperature.
In addition to the static screening effect, the dynamical screening effect also exists inside QGP. It is the dissociation of quarkonium caused by collisions with medium constituents.
The dissociation process at leading order (LO) in the coupling constant is the gluon absorption process g + H → Q +Q where H indicates a quarkonium state. It was first investigated by using large-N c expansions [6,7]. At next-leading order (NLO), inelastic scattering between quarkonium and medium constituents contributes to the dissociation l + H → l + Q +Q where l can be a gluon, light quark or light antiquark. The inelastic scattering was first studied in the quasi-free limit where the Q andQ are treated as free particles and each of them scatters independently with medium constituents [8]. Later, the interference effect was taken into account. This leads to a dependence of the dissociation rate on the relative position of the heavy quark antiquark pair [9,10]. This maps into a dependence of the inelastic scattering on the bound state wave function [11], as in the case of gluon absorption. The dissociation rate can be interpreted as an imaginary part of the QQ potential. One can also interpret the dissociation as a decoherence of the quarkonium wave function in the language of open quantum systems [12]. More recently, these dissociation rates were studied from the thermal loop corrections of the singlet propagator in potential nonrelativsitic QCD (pNRQCD) [13][14][15] by systematic weak coupling and nonrelativistic expansions. Anisotropic corrections to dissociation rates have also been considered [16][17][18].
To describe the transport of quarkonium inside QGP, one also needs to consider the in-medium recombination from unbound QQ pairs [19]. This can be modeled by detailed balance and a phenomenological factor controlling how much open heavy quarks are thermalized [20]. Recombination from parametrized nonequilibrium heavy quark distributions has also been investigated [21]. For practice, one needs heavy quark distributions from real in-medium dynamics. To this end, one can couple the transport equations of quarkonium with those of heavy quarks [22]. The inverse of the gluon absorption, the gluon emission, Q +Q → g + H has been simulated in Ref. [22] with dynamically evolving heavy quark distributions and the approach to detailed balance and equilibrium was demonstrated.
Besides dissociation and recombination, quarkonium can also diffuse inside QGP because it is approximately a color dipole, and thus not exactly color neutral. It may elastically scatter with medium constituents. The diffusion coefficient of quarkonium has been estimated in both weak and strong coupling limits [23]. However, the diffusion has not been included in phenomenological studies using transport equations.
In this paper, we consider the dissociation, recombination and diffusion of quarkonium in the same theoretical framework. We apply pNRQCD to calculate the relevant terms in a Boltzmann transport equation. The transport equation can be derived from first principles by using the open quantum system formalism and effective field theory pNRQCD [24]. We compute directly the scattering amplitudes of gluon absorption/emission, inelastic scattering and elastic scattering in pNRQCD, in contract to previous studies of dissociation rates that use the optical theorem and loop corrections to the forward amplitudes [13][14][15]. By writing out the amplitudes explicitly, we can show these amplitudes satisfy the Ward identities, which has not been explicitly shown before. We also consider the loop correction to the gluon absorption/emission. We demonstrate that the t-channel inelastic scattering is infrared safe. Furthermore, we compute the diffusion coefficient of quarkonium and reproduced the previous result in a certain limit. Finally, we discuss two mechanisms of quarkonium energy loss inside QGP: one through diffusion and the other via dissociation first and then recombination later.
The paper is organized as follows. In Section II, we briefly introduce the Boltzmann transport equation of quarkonium. Then in Section III, we explain the effective field theory pNRQCD used in the calculations. The calculation results of dissociation, recombination and diffusion are shown in the following Sections IV and V. In Section V, two mechanisms of quarkonium energy loss are also discussed. Finally the conclusions are drawn in Section VI.

II. BOLTZMANN TRANSPORT EQUATIONS
The dynamical evolution of quarkonium inside a QGP can be described by Boltzmann transport equations for the distribution function of each quarkonium state with the quantum number nl and spin s, In the following when we compute the square of the scattering amplitudes, we will average over the polarizations of non-S wave quarkonium states. So we omit the quantum number m throughout the paper. The three collision terms C + nls , C − nls and C nls represent the recombination, dissociation and diffusion of quarkonium in the medium respectively. The Boltzmann transport equation can be derived from QCD by using the open quantum system formalism and effective field theory under the assumption of Markovian process and weak coupling between the quarkonium and the medium [24]. When the quarkonium size is smaller than the screening length of the medium, the weak coupling assumption is valid.
Both the nonrelativisitic and weak coupling expansion parameters are smaller for bottomonium than charmonium. This implies that the leading order results are more reliable for bottomonium.
In the following sections, we will compute these collision terms in pNRQCD.

III. POTENTIAL NRQCD
The effective field theory pNRQCD can be systematically derived from QCD under the separation of scales: M M v M v 2 [25]. Here M denotes the mass of heavy quark, assumed to be large and v is the typical relative velocity between the heavy quark antiquark pair inside a quarkonium. For charmonium, v 2 ∼ 0.3 while for bottomonium, v 2 ∼ 0.1 [26].
To derive pNRQCD, one integrates out the scales M and M v in sequence, and applies a double expansion in the heavy quark velocity v (nonrelativistic expansion) and the inter- In this hierarchy, a bound quarkonium state can be well-defined. If T and m D are on the order of M v, the Debye screening will be so strong that the potential cannot support any bound state. In cases with rT ∼ rm D ∼ 1, the multipole expansion also breaks down.
The Lagrangian density of pNRQCD is given by where E represents the chromoelectric field and We will work to the lowest order in the expansion of v. By the virial theorem, p 2 rel /M ∼ V The potentials and Wilson coefficients V A,B in the chromoelectric dipole vertices can be obtained by matching pNRQCD with NRQCD [25,27,28]. Up to order g 2 r, The chromomagnetic vertices are suppressed by powers of v. The potential is Coulomb, which is approximately valid inside QGP since the confining part is flattened. One can improve the potentials by doing matching calculations at finite temperature or using potential models motivated from lattice calculations. 1 If the medium is moving with respect to the QQ pair at a velocity v med , the c.m. kinetic energy is still suppressed at least by one power of v if v med √ 1 − v. We assume the medium is static with respect to the QQ pair in this paper. Generalization to the case of moving medium with v med √ 1 − v can be easily worked out.
Under a gauge transformation U (R, t), therefore the Lagrangian density is invariant under a gauge transformation associated with the c.m. motion. It is worth noting that the relative motion is not gauged due to the multipole expansion. The potentials may be gauge-dependent at the matching calculation.
To make the wavefunction associated with the relative motion explicit, we do a change of basis in the relative motion by defining where N c = 3 and T F = 1 2 . We define the quadratic Casimir of the fundamental representa- Then the Lagrangian density of the singlet and octet part can be written as [28] L pNRQCD (R, t) = L kin,s + L kin,o + L int,so + L int,oo + · · · (12) The bra-ket notation saves us from writing the integral over the relative position explicitly.
The singlet and octet composite fields are quantized by where E is the eigenenergy of the state under the Hamiltonians, Eq. (5). The whole Hilbert space factorizes into two parts: one for the c.m. motion and the other for the relative motion. The wavefunctions of the relative motion can be obtained by solving Schrödinger equations, which are part of the equations of motion of the free composite fields. They can be hydrogen-like wavefunctions |ψ nl for bound singlets with the eigenenergy −|E nl |, or Coulomb scattering waves |ψ p rel and |Ψ p rel for unbound singlets and octets with the eigenenergy p 2 rel /M . No bound state exists in the octet channel due to the repulsive potential. The operators a p rel (p cm ) act on the Fock space to annihilate (create) a composite particle with c.m. momentum p cm and corresponding quantum numbers in the relative motion. These annihilation and creation operators satisfy the following commutation rules: All other commutators are zero. The Feynman rules are summarized in Fig. 1. We use the Euclidean summation and r µ = 0 when µ = 0.

IV. DISSOCIATION AND RECOMBINATION
In this section, we consider the contributions to dissociation and recombination terms in the Boltzmann equation at the order gr and g 2 r. All contributing Feynman diagrams are shown in Fig. 2.

A. Contributions at the order gr
At the order gr, only the diagram in Fig. 2a contributes. It is a gluon absorption process for dissociation and a gluon emission process for recombination. The scattering amplitude is given by The gluon is on-shell so q 0 = |q| ≡ q. The Ward identity can be easily verified: So we can compute the amplitude in any gauge we want. In Coulomb gauge, we define the square of the amplitude magnitude, summed over the gluon and octet colors a, b and the gluon polarizations The grey blob indicates the dipole interaction.
Since we average over the polarizations of quarkonium for l > 0, we expect ψ nl |r i |Ψ p rel ∝ p i rel . When integrating over the relative momentum of the heavy quark antiquark pair from the dissociation, we obtain allowing us to write To simplify the notation of the dissociation and recombination terms for a quarkonium state nls in the Boltzmann equation, we define where n B is the Bose-Einstein distribution function, g + = 1 N 2 c g s and g s is the multiplicity factor in spin: g s = 3 4 for a quarkonium with spin s = 1 and 1 4 for spin s = 0. In the definition of F + nls(a) , p cm = p Q + pQ, and p rel = are the c.m. and relative momenta of a pair of heavy quark and antiquark with momenta p Q and pQ.
We further define a "δ−derivative" symbol, first introduced in Ref. [29] δ δp i n j=1 where the second δ denotes the standard functional variation and h(p 1 , p 2 , · · · , p n ) and w(p i ) are arbitrary independent smooth functions. Then the C ± nls terms in the Boltzmann equation (1) can be written as For C + nls(a) , we further require 1 2 (x Q +xQ) = x, i.e., the position of a recombined quarkonium is given by the c.m. position of the heavy quark antiquark pair. Recombinations from the heavy quark antiquark pairs with a distance much larger than the Bohr radius |x Q − xQ| a B can be shown to be negligible [24]. In practical numerical simulations, a Gaussian dependence on |x Q −xQ| a B may be applied [22].
The dissociation rate of the quarkonium state nls from the diagram in Fig. 2a is given The recombination rate of a heavy quark into the quarkonium state nls surrounded by heavy antiquarks with the distribution fQ(xQ, pQ, t) is given by B. Contributions at the order g 2 r The light quarks are assumed massless. First we check the amplitude is independent of the gauge choice. The gauge invariance reflects in the invariance of the amplitude when the µ in the gluon propagator is replaced with µ +q µ . It has been shown in the last section IV A that the dipole interaction between the singlet and octet is invariant under such a replacement, by virtue of the Ward identity Eq. (24). What remains to show is the invariance of the vertex of the light quark. This is guaranteed by the Dirac equation:

Contributions from diagram 2b
In Coulomb gauge, We define |M (b) | 2 summed over the octet color a, the spins s 1 , s 2 , colors i, j and flavors of the incoming and outgoing light quarks and include also the contribution from antiquarks as Next we check the infrared sensitivity of the term inside the square bracket. The energymomentum conservation gives q 0 = p 1 − p 2 = |p 1 | − |p 2 | = |E nl | + p 2 rel /M , q = p 1 − p 2 = p cm − k. Assume the angle between p 1 and p 2 is θ. We have In both terms, there is no soft divergence because the binding energy |E nl | serves as a soft regulator: The first term has no collinear divergence either. The collinear divergence happens in the second term when cos θ → 1.
Physically this occurs when the momenta of both the incoming and outgoing light quarks are in the same direction. The transferred gluon is on-shell. In this case, the inelastic scattering cannot be distinguished from the real gluon process shown in Fig. 2a. As we show below in Section IV B 4, the interference between the diagram in Fig. 2a and its thermal loop correction Fig. 2l cancels this collinear divergence.

Contributions from diagrams 2c, 2d, 2e and 2f
The processes of inelastic scattering with gluons in the medium are depicted in Fig. 2c, 2d, 2e and 2f. All the four diagrams are needed for the gauge invariance. First we consider the gauge transformation of the internal gluon line in Fig. 2c. If we cut the diagram into two halves by cutting the internal gluon line, we need to show both the upper and lower parts vanish when contracted with q ρ . For the dipole interaction in the lower part, this has been shown by the Ward identity Eq. (24). For the three-gluon vertex in the upper part, it can be shown that by using q 1 · 1 = q 2 · 2 = 0.
Next we consider the gauge transformation of the external gluon line. We fix the internal gluon line to be in the Lorentz gauge.
where in the last line we have used E p = p 2 rel M + p 2 cm 4M and kept only the relative kinetic energy of the octet by neglecting the c.m. kinetic energy. This is consistent with our power counting.
Therefore we can compute these diagrams in any gauge we want. In Coulomb gauge, the zero component of the gauge field is not dynamical. So we only need to compute the diagram in Fig. 2c, We define the square of the amplitude magnitude, summed over all color indexes and polarizations: where the transverse polarization tensor is defined as P T ij (q) = δ ij −q iqj and P T 00 = P T 0i = P T i0 = 0. As in the process of inelastic scattering with light quarks, the first term in the square bracket is infrared safe because of the finite binding energy. The second term is collinear divergent when the momenta of the incoming and outgoing gluons are in the same direction. In that case, the transferred gluon is on-shell. As will be shown in Section IV B 4, the interference between the diagrams in Figs. 2a and 2m will cancel this divergence. 3. Contributions from diagrams 2g, 2h, 2i, 2j and 2k The diagrams in Figs. 2g, 2h, 2i, 2j and 2k  The phase spaces constraint by p 1 + p 2 = q 0 in Fig. 2g and q 1 + q 2 = q 0 in Fig. 2h are much smaller than those of p 1 − p 2 = q 0 in Fig. 2b and q 1 − q 2 = q 0 in Fig. 2c. The suppression of processes with two incoming light quarks or gluons has been noted before [21].

Contributions from diagrams 2l, 2m, 2n and 2o
These diagrams are the one-loop corrections of the gluon propagator. If resummed, they will give a thermal mass to the in-medium gluon. The loop correction is at the order g 2 so the whole diagram is at the order g 3 r. Their interference with the diagram in Fig. 2a will give contributions equivalent to amplitudes at the order g 2 r. We will show the interference cancels the collinear divergence in Eqs. (36) and (51). Thus there is no need to resum these diagrams here.
In Coulomb gauge, where for on-shell massless particle q 0 = |q| ≡ q and Π (l) µν is the gluon polarization tensor contributed from the fermion loop: (54) We define the interference term summed over colors and gluon polarizations In thermal field theory, the integral over l 0 is a summation in Matsubara frequency. After the summation where E 1 = |l + q|, E 2 = |l|. To see the cancellation of collinear divergence, consider the following integral that is used in the dissociation in the case of two flavors (up and down quarks) We focus on the second term in the square bracket in Eq. (57). Under a change of variables where · · · means the first term in the square bracket in Eq. (57). In the collinear limit, p 1 and p 2 are in the same direction, so p 1 −p 2 = |p 1 −p 2 |. With n B (p 1 −p 2 )(n F (p 1 )−n F (p 2 )) = −n F (p 1 )(1 − n F (p 2 )), one immediately sees I 1 cancels the collinear divergence in where |M (b) | 2 is given in Eq. (36). This is for the dissociation contribution from |M (b) | 2 .
The cancellation of the collinear divergence in the recombination can be shown similarly.
We next consider the interference between the diagrams in Figs. 2a and 2m. The amplitude iM (m) is exactly the same as iM (l) under the replacement of Π (l) with Π (m) . Π

(m)
µν is the gluon polarization tensor in Fig. 2m. After summing over the Matsubara frequencies where T A = N c , E 1 = |l+q| and E 2 = |l|. We will focus on the term with (n B (E 1 )−n B (E 2 )) and neglect the other term. Under a change of variables q 1 = l + q, q 2 = l, we can show that the collinear divergence of the integral cancels the collinear divergence in by the virtual of n B (q 1 − q 2 )(n B (q 1 ) − n B (q 2 )) = −n B (q 1 )(1 + n B (q 2 )). Cancellation of the divergence in the recombination process can be similarly shown.

Contributions from diagrams 2p, 2q, 2r and 2s
These diagrams are the one-loop corrections to the dipole interaction between the singlet and the octet. The correction is at the order g 3 r but its interference with the diagram in iM (r) = 0 + · · · (66) where T A = N c and the off-shell scheme has been used to extract the logarithmic divergence.
Only the terms with the poles are shown. Finite terms and corrections at higher orders in v 2 are omitted. Therefore which means the dipole interaction term between the singlet and octet is independent of scale at the one-loop level This has been already noted in Ref. [30]. From the matching condition, Eq. (6), we may set V A = 1 in the following, no matter of the scale involved.

Contributions from diagram 2t
The diagram in Fig. 2t describes the one-loop correction to the octet propagator at the order g 2 . The whole diagram is at the order g 3 r but its interference with the diagram 2a is equivalent to an amplitude at the order g 2 r. The one-loop correction to the octet propagator is given by 4M is the energy of the external octet field. The c.m. kinetic energy can be dropped because it is suppressed by v 2 compared with the relative kinetic energy. Then in dimensional regularization The power divergence is proportional to v 2 and neglected here, consistent with the power counting.

Summary
To write the dissociation and recombination terms in the Boltzmann equations explicitly, we define F ± nls(b) and F ± nls(c) as The g + factor and the relation between p cm , p rel and p Q , pQ are defined in Section IV A.
The collinear divergent parts in the square of amplitudes have been shown to be cancelled by the interference between the tree-level process of gluon absorption/emission and its oneloop corrections. After regularization, we can drop the terms that are originally collinear divergent if they are small, so we can write (by setting V A = 1) The dissociation and recombination terms in the Boltzmann equation from the inelastic scattering with light quarks and gluons are given by The dissociation rate of the quarkonium state nls from the inelastic scattering is given by about α s T . The c.m. momentum of the octet is at least q 1 ∼ T α s T . So the effect from the virtual octet diffusion is small and there is no need to resum gA 0 into the virtual octet.
We define the square of the total amplitude, summed over colors and polarizations of gluons To write the diffusion term in the Boltzmann equation explicitly, we define The diffusion term in the Boltzmann equation (1) can be written as We can also define the diffusion coefficient as the square of momentum transferred per unit time 3κ ≡ d 3 k 2 (2π) 3 After some simplications κ = 32 729π 5 α 2 s dq q 8 n B (q)(1 + n B (q)) P dp rel p 2 rel | Ψ p rel |r|ψ nl | 2 (|E nl + For 1S state, if the q 2 in the denominator inside the integral over p rel is neglected and the octet relative wavefunction is a plane wave, one can show where the Bohr radius a B = 2 αsC F M . If one takes the large-N c approximation, C F = 3/2, and multiplies the expression (87) by a factor of 9/8 (because when we sum over colors, there is a factor of 8, in large-N c , the factor is 9), then Eq. (87) agrees with the previous estimate using perturbative calculations in another effective field theory [23]. the denominator dominates over M v 2 and we expect κ ∝ T 5 . Furthermore, if we assume the octet relative wavefunction is a plane wave, we expect κ(1S) ∝ (M a 2 B ) 2 ∝ M −2 at high temperatures. The mass dependence of κ is also plotted for three difference heavy quark masses. At high temperature the mass scaling is approximately valid.
In principle, a quarkonium has two ways to lose energy inside a QGP. One way is the elastic scattering or diffusion. The other way is to dissociate first, lose energy as an unbound heavy quark antiquark pair and then recombine later. The former mechanism only works when the quarkonium is a well-defined bound state inside QGP. So it only makes sense when the temperature is below the quarkonium melting temperature. As shown in Fig. 4, the rate of momentum transfer due to the diffusion is very slow for Υ(1S), compared with that of open heavy quarks (κ/T 3 of heavy quarks is on the order of 1 or 10 [32]). This is also true for J/ψ. Because for J/ψ, we expect its diffusion coefficient is 10 times larger than that of Υ(1S) from the mass scaling. But J/ψ has lower melting temperature, so the diffusion coefficient would probably only make sense below 300 MeV. Therefore the latter mechanism (dissociation followed by energy loss and recombination) probably dominates the quarkonium energy loss, even though not every quarkonium finally observed has to go through this sequence of processes. Some of the primordially produced quarkonia may survive the in-medium evolution and lose almost no energy.

VI. CONCLUSION
In this paper we calculated the dissociation, recombination and diffusion terms in the Boltzmann transport equation of quarkonium. We considered the processes of gluon absorption/emission, inelastic scattering and elastic scattering with medium constituents. We computed scattering amplitudes directly in pNRQCD and showed they satisfy the Ward identities. Loop corrections were also considered. The dipole interaction is not running at the one-loop level. The interference between the gluon absorption/emission and its thermal loop corrections cancels the collinear divergence in the inelastic scattering. The inelastic scattering amplitude is infrared safe.
By choosing the Coulomb gauge, we explicitly wrote down expressions for the dissociation rate of quarkonium, the recombination rate of a heavy quark with an arbitrary heavy antiquark distribution, and the diffusion coefficient of quarkonium. We found that the diffusion coefficient of quarkonium is much smaller than that of the heavy quark. This implies that the dominant energy loss mechanism of quarkonium inside QGP is not diffusion, but rather a sequence of processes: first dissociation, then energy loss as unbound heavy quarks and later recombination.
The calculations presented here can be generalized to study the effect of a turbulent plasma on quarkonium in the early stage of heavy ion collisions, as is done for heavy quarks [33]. For a complete description of quarkonium production in heavy ion collisions, the quarkonium transport equation needs to be coupled with transport equations of heavy quarks. The Boltzmann equations of heavy quarks have been constructed and used in phenomenology [34][35][36]. By coupling these transport equations, the recombination of quarkonium will be calculated from the real-time dynamical heavy quark distributions rather than phenomenological models. The coupled Boltzmann transport equations have been used to study Upsilon production at both RHIC and LHC and can describe the experimental data [37]. In future work, we will solve the coupled Boltzmann equations and study charmonium production in heavy ion collisions.