Quarkonium Semiclassical Transport in Quark-Gluon Plasma: Factorization and Quantum Correction

We study quarkonium transport in the quark-gluon plasma by using the potential nonrelativistic QCD (pNRQCD) effective field theory and the framework of open quantum systems. We argue that the coupling between quarkonium and the thermal bath is weak using separation of scales, so the initial density matrix of the total system factorizes and the time evolution of the subsystem is Markovian. We derive the semiclassical Boltzmann equation for quarkonium by applying a Wigner transform to the Lindblad equation and carrying out a semiclassical expansion. We resum relevant interactions to all orders in the coupling constant at leading power of the nonrelativistic and multipole expansions. The derivation is valid for both weakly coupled and strongly coupled quark-gluon plasmas. We find reaction rates in the transport equation factorize into a quarkonium dipole transition function and a chromoelectric gluon distribution function. For the differential reaction rate, the definition of the momentum dependent chromoelectric gluon distribution function involves staple-shaped Wilson lines. For the inclusive reaction rate, the Wilson lines collapse into a straight line along the real time axis and the distribution becomes momentum independent. The relation between the two Wilson lines is analogous to the relation between the Wilson lines appearing in the gluon parton distribution function (PDF) and the gluon transverse momentum dependent parton distribution function (TMDPDF). The centrality dependence of the quarkonium nuclear modification factor measured by experiments probes the momentum independent distribution while the transverse momentum dependence and measurements of the azimuthal angular anisotropy may be able to probe the momentum dependent one. We discuss one way to indirectly constrain the quarkonium in-medium real potential by using the factorization formula and lattice calculations. The leading quantum correction to the semiclassical transport equation of quarkonium is also worked out. The study can be easily generalized to quarkonium transport in cold nuclear matter, which is relevant for quarkonium production in eA collisions in the future Electron-Ion Collider.


Introduction
Heavy quarkonium has been used as a probe of the quark-gluon plasma (QGP) in heavy ion collisions for many years. The basic idea is the static screening effect in the hot plasma [1,2]: the real part of the attractive potential between the heavy quark-antiquark pair (QQ) is significantly suppressed at high temperature and quarkonium "melts". Therefore, suppression of quarkonium production in heavy ion collisions can be used as a signature of the QGP formation. The melting temperature can be defined for each quarkonium state as the minimum temperature when the state becomes unbound in the plasma. Since all quarkonium states have varying sizes, they are influenced by the static screening differently and thus have distinct melting temperatures, ordered by their sizes (or binding energies). A sequential melting pattern is expected where shallower quarkonium states melt at lower temperatures. However, the simple static screening picture is complicated by several other factors: cold nuclear matter effects, as well as quarkonium dissociation and recombination inside the plasma. One cold nuclear matter effect is the nuclear modification of parton distribution functions (PDF) inside heavy nuclei. Dissociation and recombination are hot medium effects. It has been shown that static screening and dissociation can be generated simultaneously from thermal loop correction to the quarkonium in-medium propagator [3,4]. Thus if the static screening is included in the study, dissociation should also be taken into account for consistency. Furthermore, if a quarkonium state can still exist as a well-defined bound state above T c ∼ 150 MeV, which is a rough estimate of the transition temperature from the QGP phase to the hadronic phase 1 , then (re)generation of this quarkonium state inside the QGP should also be possible [5]. It is expected that deeply bound states can start to (re)generate at high temperatures and do not have to wait until T c , when light particles hadronize.
Quarkonium suppression has been intensively investigated in experiments at both the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC). Semiclassical transport equations that account for static screening, dissociation and recombination [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], as well as statistical hadronization models [21,22], have been used to study quarkonium production in heavy ion collisions and achieved great success in phenomenology. Anisotropic effects in the QGP have also been explored [23][24][25]. Recently, we derived the semiclassical transport equations using nonrelativistic QCD in the potential regime and the open quantum system formalism. The validity of this approach relies on a hierarchy of scales: M M v M v 2 , T , where M is the heavy quark mass, v is the velocity of the quarks in the quarkonium and T is the temperature of the plasma [26]. Long after quarkonium suppression was proposed as a diagnostic for the existence of the QGP, several other observables such as jet quenching and elliptic flow have been studied which more convincingly establish the existence of the QGP (see Refs. [27][28][29][30] for recent reviews on each topic). We no longer need to use quarkonium as a probe to answer the question whether the QGP is formed in heavy ion collisions. Then the natural question to ask is what we can learn about the QGP from quarkonium measurements. Since now the hot medium effects contain static screening, dissociation and recombination, a simple answer seems implausible.
In deep inelastic scattering (DIS), electrons are shot onto protons to probe the inner structure of proton. The reason why this works is the possibility of constructing factorization theorems in certain kinematic limits (for a general discussion of factorization, see Ref. [31]). These factorization theorems express physical observables as convolutions of perturbatively calculable cross sections with parton distribution functions (PDF) that can be expressed as matrix elements of operators within the proton state. Thus measurements of these cross sections determine specific correlation functions within the proton. Furthermore, different observables can probe different kinds of parton distributions of the proton. Studying these distributions such as the transverse momentum dependent (TMD) PDF will be among the central scientific goals of the future Electron-Ion Collider.
This then leads us to ask: what correlation functions are measured when we study the production of quarkonium within the QGP? The same question could be asked of quarkonium production within, say, cold nuclear matter. We will use effective field theory (EFT) and the open quantum system framework to answer this question. For processes involving light-like partons, Soft-Collinear Effective Theory has been widely applied to study factorization in various processes [32][33][34][35][36][37][38][39][40][41][42]. In our case since the quarkonium is nonrelativistic we will use potential nonrelativistic QCD (potential NRQCD or pNRQCD) [43][44][45]. This EFT has been used to study static screening and dissociation [46][47][48]. Also, it has been used in a Lindblad equation to define a transport coefficient of quarkonium [49]. New developments of nonrelativistic EFT for quarkonium can be found in Refs. [50,51]. The open quantum system framework describes the dynamics of a subsystem, interacting with an environment. When the environment is traced over, the subsystem evolves as an open system. The open quantum system framework has been recently applied to study quarkonium transport in the QGP [26,[52][53][54][55][56][57][58][59][60][61][62][63][64][65] and provides new insight to our understanding of quarkonium transport. The open quantum system approach can be thought of as an extension of the Schrödinger equation with a complex potential [66,67] and it takes into account both correlated and uncorrelated recombination consistently. The difference between correlated and uncorrelated recombination is discussed in Ref. [68].
We will derive the semiclassical Boltzmann equation for quarkonium in the QGP by applying a Wigner transform (a Gaussian smearing is required for sustaining positivity) to the Lindblad equation, which is the evolution equation for the open system. The interaction between quarkonium (subsystem) and the thermal QGP (environment) is weak according to the power counting of pNRQCD under the assumed hierarchy. Therefore the density matrix of the total system can be assumed to factorize into the density matrix of the subsystem and that of the environment. Furthermore, in the weak coupling (between the subsystem and the environment) limit, the time evolution of the subsystem can be shown to be Markovian. We will work to leading order in the power counting parameter of pNRQCD but resum interactions to all orders in the strong coupling constant, which can be written in terms of Wilson lines. We will show that the reaction rates in the Boltzmann equation factorize into a quarkonium dipole transition function and a chromoelectric gluon distribution function of the thermal QGP. For the differential reaction rate, the chromoelectric gluon distribution function is momentum dependent and has spatially separated chromoelectric fields connected via a staple-shaped Wilson line extending to the time infinity. For the inclusive reaction rate, the function is momentum independent and the Wilson line collapses into a straight timelike line of finite length. The structures of the Wilson lines are very similar to the case of gluon TMDPDF and gluon PDF, respectively, except that here the Wilson lines lie along a timelike direction rather than the lightcone. The momentum independent chromoelectric gluon distribution function has been studied in Ref. [49] as one quarkonium transport coefficient and in Refs. [69,70] as the heavy quark diffusion coefficient. A recent lattice calculation of the heavy quark diffusion coefficient can be found in Ref. [71]. What is new in this paper is the definition and discussion of the momentum dependent chromoelectric structure function. We will discuss one application of the factorized rate to constrain the in-medium real potential between a QQ pair. A point of emphasis is that the factorization of the reaction rate crucially depends on the factorization of the density matrix into the subsystem density matrix and the environment density matrix. This is generally believed to be true in the weak coupling limit, where the weak coupling is between the subsystem and the environment. The subsystem and the environment themselves can be strongly coupled.
This paper extends our earlier work [26] which derived semiclassical transport equations from pNRQCD and the open quantum system formalism. Major improvements compared to Ref. [26] include: • We resum the A 0 interaction and the interaction mediated by Coulomb modes between the octet QQ state and the thermal QGP. This allows us to define the chromoelectric gluon distribution function nonperturbatively and construct the factorization. The derivation of Ref. [26] only works for a weakly coupled QGP while the derivation we show here is also valid for a strongly coupled QGP. A formalism compatible with a strongly coupled QGP is crucial for phenomenological studies of the future experiments at RHIC and LHC, especially the one carried out by the sPHENIX collaboration since the plasma temperature achieved at RHIC is not very high.
• We carry out a systematic semiclassical expansion for the recombination term. For dissociation, no semiclassical expansion is needed. The semiclassical expansion is a gradient expansion. We work out the leading quantum correction to the semiclassical Boltzmann equation.
The paper is organized as follows: In Section 2, we will discuss the hierarchy of scales in the problem and briefly explain pNRQCD with a focus on the power counting. Then in Section 3, a short introduction to the open quantum system framework will be given. Derivation of the semiclassical Boltzmann equation will be elucidated in detail in Section 4. Factorization of the reaction rates will also be discussed there. We will derive the leading quantum correction to the Boltzmann equation in Section 5. Finally, we will summarize and draw conclusions in Section 6.

Separation of Scales and Potential NRQCD
We consider the following hierarchy of scales: M M v M v 2 , T, Λ QCD , where M is the heavy quark mass, v is the typical relative velocity between the heavy quark pair inside quarkonium, T is the temperature of the medium and Λ QCD is the nonperturbative scale of QCD. In vacuum, M M v M v 2 is the standard separation of scales for nonrelativistic heavy quarks [72]. For both charmonium and bottomonium, M v 2 ∼ 500 MeV. In current heavy ion collision experiments, T ∼ 500 MeV is roughly the highest temperature achieved in the early stage of the medium expansion. Naively, we would expect M v 2 T to be approximately valid during the whole lifetime of QGP. But due to the static plasma screening, the binding energy decreases as the temperature increases, and M v 2 ∼ 500 MeV is probably no longer true except for the bottomonium ground state at T ∼ 500 MeV. So we do not specify the hierarchy between M v 2 and T . However, we still believe M v T is a relevant hierarchy for quarkonium transport. The reason is the following: when the temperature is high and the quarkonium size is large, this hierarchy M v T may be violated and the interaction between the heavy quark pair is significantly suppressed. Loosely bound quarkonium states at such a high temperature, even if they exist in the medium, are no longer well-defined bound states, the dissociation rate is very large and formation is ineffective and can be neglected. The state of the heavy quark and antiquark will be more like two open heavy quarks and their dynamics is governed by the transport of open heavy flavors rather than the transport of quarkonium. Only when the QGP expands and the temperature drops, does the quarkonium formation become effective and transport become applicable. At the end of the QGP expansion when hadronization starts to occur, the temperature is about T c ∼ 150 MeV and every quarkonium state should regain their vacuum properties. Therefore, the hierarchy M M v M v 2 , T, Λ QCD should be valid for ground and lower excited quarkonium states in the later stage of the QGP expansion, when the formation of these states becomes effective.
Under the separation of scales M M v M v 2 , T, Λ QCD , one can construct an effective field theory of QCD, pNRQCD, by systematic nonrelativistic expansion in v and multipole expansion in rT ∼ T M v , where r is the typical size of a quarkonium state. 2 In our assumed hierarchy of scales, the effective field theory can be constructed perturbatively at the scale M v. The Lagrangian can be written as where higher order terms in the expansion are neglected. We will work at leading power in the nonrelativistic and multipole expansions throughout the paper. Here E = E A T A (A is the adjoint color index) represents the chromoelectric field and The gluon and light quark parts are just QCD with momenta M v. The matching coefficients are V A = V B = 1 at lowest order in the coupling constant at the scale M v. The composite fields for the quarkonia are the color singlet S(R, r, t) and color octet O(R, r, t) where R denotes the center-of-mass (c.m.) position and r the relative coordinate. The composite fields here are constructed by projecting onto the proper color space, a heavy quark field and a heavy antiquark field at the same time, connected by a spatial Wilson line. We will assume the medium is invariant under translation so the existence of the medium does not break the separation into the c.m. and relative motions. H s and H o denote the color singlet and octet Hamiltonians of the relative motion. At leading order in the nonrelativistic expansion, H s,o = − ∇ 2 r M + V s,o (r), where V s,o are the singlet and octet potentials and are attractive and repulsive, respectively. So only a color singlet QQ pair can be bound. Both potentials are spin independent. Thus the pNRQCD Lagrangian is invariant under heavy quark spin symmetry and we will ignore spin quantum numbers in this paper. In the v expansion, the Fock state |QQg of quarkonium, in which the QQ pair is a color octet state, is suppressed by at least v 2 in probability with respect to the leading Fock state |QQ in which the QQ pair is a color singlet state. Therefore, at leading order in the v expansion, which is the order we are working, quarkonium can only be a color singlet QQ pair. Furthermore, the QGP is a deconfined phase of QCD, where light quarks and gluons are liberated. Thus it is reasonable to assume no gluon can be bound with an octet QQ pair.
For our construction that is at leading power in v and rT but resummed to all orders in α s (M v 2 , T ), we do the following redefinitions: where the Wilson line in the fundamental representation is defined by where W [(R,t),(R,t 0 )] denotes a Wilson line in the adjoint representation that has a straight path from (R, t 0 ) to (R, t). So far t 0 is an arbitrary constant, but later we will follow arguments given in Ref. [73] and show the results are independent of the choice of t 0 . Now we can simplify the octet kinetic term The octet kinetic term can then be rewritten as The new dipole interaction between the singlet and octet is given by If we define the dipole interaction term can be written as For later convenience, we introduce the "bra-ket" notation: for any function f of r. At the current leading power calculation, we will use 2Nc . The quantization of the fields is given by where E is the eigenenergy of a state that will be explained below. The operators a p rel (p cm ) act on the Fock space to annihilate (create) composite particles with the c.m. momentum p cm and the corresponding quantum numbers in the relative motion. These quantum numbers can be nl for bound singlet states, p rel for unbound singlet states and color A and p rel for unbound octet states. When we compute the matrix elements, we will average over the polarizations of non-S wave quarkonium states. So in our notation, we omit the quantum number m l of the bound singlet state. In the octet channel no bound state exists because of the repulsive octet potential. The corresponding wavefunctions of the relative motion are |ψ nl and |ψ p rel for color-singlet states and |Ψ p rel for color-octet states. They can be obtained by solving the equations of motion of the free composite fields, which are Schrödinger equations. The eigenenergies are E = −|E nl | and E = p 2 rel M for the bound and unbound states, respectively, with higher order terms in v neglected. Here E nl is the binding energy of the bound state |ψ nl .
The dissociation and recombination of quarkonium occur via the dipole interaction between the color singlet and octet states. As explained above, quarkonium is treated as a color singlet QQ pair in this work, consistent with the leading power (in v) calculation. The dipole vertex scales as rT where r ∼ 1 M v is the typical quarkonium size and T is the typical energy and momentum scale of excitations in the plasma (the excitation energy and momentum comes from the derivative in the chromoelectric field). In our assumed hierarchy M v T , the dipole vertex between the singlet and octet states scales as rT 1. Thus, the interaction between quarkonium and the medium is weak. This argument is about the weak coupling between quarkonium and the QGP and is valid even if quarkonium and the QGP are strongly coupled themselves.

Open Quantum Systems
In this section, we briefly introduce the framework of open quantum systems. We consider a total system consisting of a subsystem and an environment (thermal bath). The total Hamiltonian is given by where H S is the subsystem Hamiltonian, H E is the environment Hamiltonian, and H I contains the interactions between the subsystem and the environment. The interaction Hamiltonian is assumed to be factorized as follows: α where α denotes all quantum numbers. (For local quantum field theory, this is generally true and α includes the spatial coordinates.) The operators O Here ρ E is the density matrix of the environment. Each part of the Hamiltonian is assumed to be Hermitian.
The von Neumann equation for the time evolution of the density matrix in the interaction picture is given by We will omit the superscript "(int)" in the following. The symbolic solution is given by where the evolution operator is

4)
and T is the time-ordering operator. The starting time t = 0 is arbitrary in the Markovian limit. So for later convenience, we will shift the starting time to −t/2 and obtain We will assume the subsystem and the environment are weakly interacting. This is valid for quarkonium inside the QGP with the hierarchy of scales M M v M v 2 , T, Λ QCD because the heavy quark pair interacts with the medium chromoelectric field via the dipole interaction, which scales as rT ∼ T M v 1. Then we can assume the initial density matrix factorizes which is generally true for weakly coupled systems (factorization breaking terms come at higher orders in the coupling). The environment density matrix is assumed to be in thermal equilibrium thus ρ E = 1 Z e −βH E where β = 1/T , and is time-independent. We will use O (E) T to denote Tr E [O (E) ρ E ] from now on, where the subscript T indicates the environment is a thermal bath at temperature T . If we expand the interaction to second order in perturbation (which corresponds to leading power in rT ) and take the partial trace over the environment degrees of freedom, we obtain the Lindblad equation: Each term is defined as where {|a } forms a complete set of states in the Hilbert space of the subsystem. We can choose them to be the set of eigenstates of H S . To calculate these in pNRQCD at finite temperature, we need the following dictionary: where i = x, y, z is the spatial component and the superscript A denotes the adjoint color index. The complete set of states |a in pNRQCD is where the label 1 means the state is a color singlet while A means the state is in a specific color octet state. Since the redefined octet field is dressed with a Wilson line, Eq. (3.17) contains a Wilson line for the octet state created by the original octet field |p cm , p rel , A (t), where c B † p rel (p cm ) is creation operator for the redefined octet field. With these preparations, we are ready to derive the semiclassical transport equation and construct factorized reaction rates.

Transport Equation and Factorized Rates
We will first sketch the derivation of the semiclassical transport equation and then explain it in detail. The semiclassical transport equation for quarkonium can be obtained from the Lindblad equation by first making the Markovian approximation, which is valid when the environment correlation time is much smaller than the subsystem relaxation time. The Markovian approximation is generally true when the subsystem interacts weakly with the environment. The environment correlation time is given by ∼ 1/T while the subsystem relaxation time can be estimated as where we use the dipole interaction between quarkonium and the medium as the relaxation process. It can be seen that the Markovian approximation is valid in our power counting. Under the Markovian approximation, the Lindblad equation (3.7) in the Schrödinger picture when t is small turns to [26] f nl x, k, if a Wigner transform is applied to the subsystem density matrix 3 Dividing (4.1) by t and taking the limit t → 0, we obtain the Boltzmann transport equation at t = 0. Since the starting time is arbitrary, we interpret a similar equation at an arbitrary time t as The Boltzmann equation describes the phase-space evolution of the quarkonium state with the quantum number nl. To solve this equation, one needs to couple it with transport equations for unbound singlet and octet QQ's, since the recombination term depends on the distribution of unbound QQ's (see Section 4.3). In Eqs. (4.1, 4.3), the free streaming term − k 2M · ∇ x f nl (x, k, t) comes from the von-Newmann evolution of the density matrix in the Schrödinger picture, i.e., −i[H S + a,b σ ab L ab , ρ S ]. This has been explained in Ref. [26] in detail, and will not be repeated here. We will explain the dissociation C − nl and recombination C + nl collision terms in the following. In the derivation of C + nl , a semiclassical expansion will also be used.

Dissociation
We will first work out the dissociation term C − nl from a,b,c,d γ ab,cd (t)L † cd L ab ρ S (−t/2). When we sandwich it between k 1 , nl, 1| and |k 2 , nl, 1 , as in the Wigner transform (4.2), we find the state |d in L † cd must be |k 1 , nl, 1 and |c = |a . Since at the order we are working, the only vertex that couples the color singlet state and the environment is the singlet-octet dipole interaction, we must have |c = |a = |p cm , p rel , A (t/2) where the A denotes the color of the octet state. By the same argument, we find |b = |k 3 , n l , 1 . So we need to compute Under the Markovian approximation, we can set the upper limits of the two time integrals to infinity t → +∞. 4 Then the octet state |p cm , p rel , A (t/2) can be thought of as an "asymptotic" outgoing state, which is defined at t → +∞ by the original octet field creation operator O A . Using Eqs. (2.4, 2.15, 3.15, 3.17), we find Plugging everything into Eq. (4.4) gives Together with the redefined chromoelectric fields (2.8, 2.9), the term inside the partial trace over the environment degrees of freedom can be written as (in the Markovian limit where E i = E A i T A and · · · T = Tr E (· · · ρ E ). In the last line we defined the chromoelectric gluon distribution function g E++ i 1 i 2 of the thermal QGP. In the definition of g E++ i 1 i 2 , the chromoelectric fields at different spacetime points are dressed with timelike Wilson lines extending to infinity. A spatial Wilson line connecting the open ends of the timelike Wilson lines is needed to restore gauge invariance, however, it cannot be generated from the field redefinition applied above. The spatial Wilson line at the infinite time is generated from resumming offshell Coulomb modes which are exchanged between the color octet pair and the medium. Details of the calculation are shown in Appendix A. The timelike and spatial Wilson lines together form a staple shape, as shown in Fig. 1. The shape of the Wilson lines is similar to the case of the gluon TMDPDF, though the time direction here is along the real time, rather than the lightcone time. If we assume the thermal QGP is invariant under spacetime translation, we can define We can plug the Fourier transform (4.9) into (4.7) and carry out the integrals over R i and t i where i = 1, 2. The spatial integrals lead to (2π) 6 δ 3 (k 1 −p cm +q)δ 3 (k 3 −p cm +q). The two time integrals in the Markovian limit t → +∞ give ∝ δ(E nl − E p + q 0 )δ(E n l − E p + q 0 ). If we assume quarkonium states with different quantum numbers nl are non-degenerate (they have different binding energies), then the summation over n l will fix n = n and l = l. Finally using the following integral we find (4.7) can be simplified to Defining the quarkonium dipole transition function and making the Wigner transform (4.2) (by setting k 1 = k + k 2 , k 2 = k − k 2 and shifting p cm → p cm + k 2 ), we find (4.11) turns to So far, we only consider the a,b,c,d γ ab,cd (t)L † cd L ab ρ S (−t/2) term in the Lindblad equation (3.7). The other term in the anti-commutator a,b,c,d γ ab,cd (t)ρ S (−t/2)L † cd L ab can be shown to give the same result (4.13). Their sum will cancel the factor of 1 2 in Eq. (3.7). So in the Markovian limit, after the Wigner transform, the anti-commutator term in the Lindblad equation leads to tC − nl (x, p, −t/2), as previously shown in Eq. (4.1). It should be pointed out that in the derivation of the dissociation collision term, no semiclassical approximation is made.

Recombination
To derive the recombination term C + nl from the Lindblad equation (3.7), we need to sandwich a,b,c,d γ ab,cd (t)L ab ρ S (−t/2)L † cd between k 1 , nl, 1| and |k 2 , nl, 1 and apply the Wigner transform (4.2). We find the state |a in L ab is |k 1 , nl, 1 and |c in L † cd is |k 2 , nl, 1 . Since at the order we are working, the only vertex that couples the color singlet and the environment is the singlet-octet dipole interaction, we must have |b = |p 1cm , p 1rel , A 1 and |d = |p 2cm , p 2rel , A 2 where the A i denotes the color of the octet state. We need to compute (4.14) Now we note that the octet states |p i cm , p i rel , A i (i = 1, 2) are defined at −t/2 by the original octet field O A i , since they sandwich ρ S (−t/2). Similar to Eqs. (4.5, 4.6), we can show where E p i = (p i rel ) 2 /M . Plugging into (4.14) and using the Wilson lines dressed on the chromoelectric fields gives where the function The newly defined function (g E−− i 2 i 1 ) A 2 A 1 is different from the previously defined chromoelectric gluon distribution of the thermal QGP in two aspects: First, (g E−− i 2 i 1 ) A 2 A 1 has open color indexes and thus one may worry that it is gauge dependent. However, the timelike Wilson line together with the spatial Wilson which is explained in Appendix A connects the chromoelectric field to a point at R = ∞, t = −∞. So the combination W [(∞,−∞),(R,t)] E(R, t) transforms as an adjoint representation at R = ∞, t = −∞, which can be set to be a trivial transformation by a choice of the global gauge. Therefore, both WE and EW are gauge invariant objects independently in the definition of (g E−− i 2 i 1 ) A 2 A 1 . Second, the end point of the Wilson lines along the time axis is t = −∞ in (g E−− i 2 i 1 ) A 2 A 1 rather than t = +∞. The physical interpretation of the different end points of the time axis is as follows: In quarkonium dissociation, the color octet state is a final state and the Wilson line resums the A 0 interaction between the octet state and the medium, which is not suppressed by the nonrelativistic expansion. Since only final state interactions are involved, the Wilson lines go to t = +∞. In quarkonium recombination, the color octet state is an initial state and the Wilson line resums the A 0 interaction before recombination occurs, which is an initial state interaction. But one should keep in mind that the density matrix for the incoming state may be off-diagonal in color space and still contribute to recombination.
Recombination from the state with an off-diagonal color density matrix is genuinely a quantum effect with no classical analog. To further simplify the recombination term (4.17) and derive the recombination term in the Boltzmann equation, we need to make a semiclassical approximation, which will be explained in the next subsection.

Semiclassical Approximation in Recombination
As discussed above, we will make semiclassical approximation to write the recombination term (4.17) as a collision term in the Boltzmann equation (4.3). First we approximate the subsystem density matrix by its diagonal piece in the color space where the superscript (8) indicates the density matrix is a color octet state. With this approximation, we can contract the color indexes in (g E−− i 2 i 1 ) A 2 A 1 and define (4.20) The function g E−− i 2 i 1 (t 2 , t 1 , R 2 , R 1 ) is another chromoelectric gluon distribution function of the thermal QGP, similar to the previously defined g E++ i 1 i 2 (t 1 , t 2 , R 1 , R 2 ). The only difference is the orientation of the Wilson line. For g E++ i 1 i 2 (t 1 , t 2 , R 1 , R 2 ), the Wilson line goes to t → +∞ while for g E−− i 2 i 1 (t 2 , t 1 , R 2 , R 1 ), the Wilson line comes from t → −∞. A spatial Wilson line connecting the end points of the timelike Wilson lines is also needed for gauge invariance and can be generated from resumming Coulomb modes, as in the case of g E++ . Detailed calculations of the spatial Wilson line at infinite time can be found in Appendix A. The timelike and spatial Wilson lines in the definition of g E−− i 2 i 1 (t 2 , t 1 , R 2 , R 1 ) form a staple shape, as plotted in Fig. 1.
Using the assumption of translational invariance, we have Plugging everything into (4.17) and integrating over R 1 and R 2 leads to When applying the Wigner transform to (4.22), we set k 1 = k + k /2 and k 2 = k − k /2. Changing variables p 1cm → p cm + p cm /2 and p 2cm → p cm − p cm /2, we find Applying the Wigner transform, we find where f (8) QQ (x cm , p cm , x rel , p 1rel +p 2rel 2 , −t/2) is the phase space distribution function of a color octet QQ pair with center-of-mass and relative positions and momenta x cm , p cm , x rel , p 1rel +p 2rel 2 . If the color reaches thermal equilibrium, statistically we will have where f QQ is the distribution function of an unbound QQ pair that can be either a color singlet state or an octet state. Now we take a crucial step in the derivation of the semiclassical transport equation: the semiclassical expansion, also known as the gradient expansion. A general discussion of the gradient expansion in the derivation of semiclassical transport equations can be found, for example, in Ref. [74]. We will expand f (8) QQ around some x rel = x 0 and assume the distribution varies slowly as x rel changes f (8) QQ (x cm , p cm , x rel , where higher order terms in the gradient expansion are omitted. In this section, we will focus on the leading term in the gradient expansion. The next-to-leading term, which will be discussed in Section 5, corresponds to a quantum correction to the semiclassical Boltzmann equation. In practice, we want to choose x 0 such that quantum corrections are minimized. We will choose x 0 = 0 when we compute the correction in the next section. The derivation shown in Ref. [26] uses the gradient expansion implicitly, by assuming the distribution function of a QQ pair is uniform in the relative position. This assumption of uniformity is exactly the leading term in the gradient expansion (4.26). The argument given in Ref. [26] relies on a large diffusion rate for open heavy quarks and the angular dependence of the octet state wavefunctions |Ψ p rel (see Eq. (D8) of Ref. [26]). Thus it is not obvious how to generalize the derivation in Ref. [26] to incorporate quantum corrections.
Here the derivation is based on a gradient expansion and higher order corrections can be worked out systematically.
Plugging the leading term in the gradient expansion back into (4.24), we find the integral over x rel can be done trivially which gives (2π) 3 δ 3 (p 1rel − p 2rel ). Now we can carry out the time integrals in (4.23) in the Markovian limit when p 1rel = p 2rel ≡ p rel . The time integrals in the Markovian limit have been explained in Section 4.1. We can show after the Wigner transform (4.24), Eq. (4.23) turns to where E p = (p rel ) 2 /M and we defined the collision term for recombination C + nl in the Boltzmann equation (4.3). The structure of (4.27) is very similar to that of (4.13). We will discuss these two collision terms in the next subsection.

Factorization of Reaction Rates
In the previous subsections, we have derived the collision terms in the semiclassical Boltzmann equation for dissociation and recombination. From Eqs. (4.13, 4.27), we have Then the dissociation and recombination rates can be defined. General expressions for the reaction rates can be found in Refs. [75,76]. We will first study the dissociation rate of a quarkonium state with position x and momentum k, which can be written as , (4.30) where the rate depends on position and time via the dependence of the QGP temperature on position and time. 5 The QGP temperature determines the chromoelectric gluon distribution function. Using (4.28), we find The summation over i 1 , i 2 can be further simplified. So far we have not written out the dependence on the third component of the orbital angular momentum m l explicitly. In practice, we will average over m l since current heavy ion experiments do not measure m l . Temporarily restoring the m l dependence in the bound state wavefunction, we obtain (note that in the integrand, the only dependence onp rel is in the dipole transition function) where dΩ p rel = d cos θ p rel dφ p rel . Defining we can write the dissociation rate as where G E++ is the integrated chromoelectric gluon distribution function: In the integrated chromoelectric gluon distribution function G E++ , the spatial index i is summed over and W denotes a Wilson line in the adjoint representation. The Wilson lines connecting R 1 and R 2 shown in Fig. 1 overlap at the same position and only the part between t 1 and t 2 remains after cancellation. The differential rate can be written as Eqs. (4.34, 4.36) are important results of this paper. They show the inclusive and differential dissociation rates of quarkonium factorize into the quarkonium dipole transition function d nl and the chromoelectric gluon distribution function of the QGP. In the inclusive rate, the gluon distribution function is momentum independent while in the differential rate, it is momentum dependent. The connection between the Wilson line structures in the definitions of the momentum independent and momentum dependent chromoelectric gluon distribution functions is very similar to the relation between the gluon PDF and the gluon TMDPDF. Through the use of this factorization theorem, experimental measurements of quarkonium nuclear modification factors probe the chromoelectric gluon distribution function of the QGP. The centrality dependence of the quarkonium nuclear modification factor probes the momentum independent distribution while the transverse momentum dependence and measurements of the azimuthal angular anisotropy may be able to probe the momentum dependent distribution, since both the differential dissociation and recombination rates depend on the momentum dependent distribution.
One application of the factorization formula (4.34) is to combine it with lattice QCD calculations to constrain the real part of the in-medium potential of quarkonium. The thermal width extracted from lattice QCD calculations of the spectral functions contains both the dissociation rate and the diffusion rate. In the diffusion process, the singlet QQ pair exchanges some momentum with the medium but does not break up. The diffusion process is suppressed with respect to dissociation when the temperature is small, i.e., rT 1 where r ∼ 1/(M v) is the typical size of a quarkonium state. In our power counting, the dissociation amplitude scales as rT while the diffusion amplitude scales as (rT ) 2 [76]. So at leading order in the multipole expansion, the thermal width is equal to the dissociation rate R − . Then if the chromoelectric gluon distribution function G E++ can be calculated in lattice QCD, we can combine the two lattice calculations and use the factorization formula (4.34) to constrain the quarkonium dipole transition function d nl . The dipole transition is between the bound state and the unbound scattering state. Their wavefunctions can be solved by using parametrized in-medium real potentials. So we can calculate d nl with different parametrizations of the real potential to compare with the one constrained by lattice QCD calculations. This method can indirectly constrain the inmedium real potential of quarkonium. It may also be used to test the consistency between the real and imaginary parts of the potential calculated in lattice QCD [77]. In practice, one may first carry out the above analysis for Υ(1S) at low temperature, where the power counting parameter is small and the framework presented here is under good theoretical control. Recent lattice QCD calculations of the thermal width and other developments for bottomonium at finite temperature can be found in Refs. [78][79][80].
Next we will study the recombination term C + . Using (4.32) and integrating over p cm , we find (4.29) becomes QQ (x cm , k − q, x 0 , p rel , t) , (4.37) which factorizes into three pieces: dipole transition function, chromoelectric gluon distribution function and the octet QQ distribution function. Eq. (4.37) should be thought of as a differential recombination process because the final state momentum, i.e., the quarkonium momentum, k, is not integrated over. Integrating over k leads to where G E−− is the integrated chromoelectric gluon distribution function and n QQ is a density. They are given by QQ (x cm , k, x 0 , p rel , t) . (4.40) We note that the integrated chromoelectric gluon distributions G E++ and G E−− are the same. Converting Eqs. (4.37, 4.38) into differential and inclusive recombination rates of a heavy quark requires knowledge of the relation between the two-particle QQ distribution and the one-particleQ distribution. We will not pursue writing out the recombination rates explicitly here. Relevant formulas can be found in Ref. [76].
Finally, we comment on the scale dependence of each component in the factorization formula. The chromoelectric gluon distribution function has a natural scale T , the plasma temperature. So we need to compute the dipole transition function d nl at the scale T . The EFT pNRQCD is constructed by matching at the scale M v, we need to solve the renormalization group equation for d nl from M v to T . It has been shown that at one loop, no extra renormalization is needed for the dipole interaction vertex beyond the renormalization of the strong coupling constant α s [81]. We believe this is true to all orders due to the reparametrization invariance of the Lagrangianψ(iD 0 − D 2 2M )ψ for a single heavy quark field ψ, from which the leading pNRQCD Lagrangian is derived (see the derivation of pNRQCD in [45]).

Quantum Correction to Semiclassical Transport
In this section, we work out the leading quantum correction to the semiclassical Boltzmann transport equation. For dissociation, no semiclassical expansion is applied. For recombination, we make two semiclassical approximations. The first one is (4.19), where we assume the octet state density matrix is diagonal in the color space. We want to point out that this assumption is not necessary if one works at leading order in the coupling constant. The reason why we have to make this semiclassical assumption is the open color indexes in (4.18). If we only keep the leading terms (in the coupling constant) in (4.18), we can set all the Wilson lines to be unity. Then (4.18) becomes which is proportional to δ A 1 A 2 up to higher order corrections (in the coupling constant). Therefore, at leading order, only the diagonal entries of the color density matrix contribute to recombination. But in general, off-diagonal entries can contribute. To derive the recombination collision term in semiclassical Boltzmann equation, we have to approximate the octet density matrix to be diagonal in color.
The second semiclassical approximation is the gradient expansion (4.26). So far, we only take the leading term in the gradient expansion. We now work out the recombination term from the next-to-leading term: (x rel − x 0 ) · ∇ x 0 f (8) QQ (x cm , p cm , x 0 , p 1rel +p 2rel 2 , t). For simplicity, we will set x 0 = 0. In practice, one wants to choose a x 0 such that the gradient expansion converges fastest. With the next-to-leading term, (4.24) becomes (between the subsystem and the environment) limit, the total density matrix factorizes into the subsystem density matrix and the environment density matrix. We demonstrated how the Lindblad equation for quarkonium as an open system turns into a Boltzmann equation after taking the Markovian limit and applying a Wigner transform (a Gaussian smearing is required for maintaining positivity). We justified the Markovian approximation using our power counting. The derivation is valid for a strongly coupled QGP at leading power of the nonrelativistic and multipole expansions since we resummed relevant interactions to all orders in the coupling constant at the scale M v 2 and T . Reaction rates in the Boltzmann equation factorize into the quarkonium dipole transition function and the chromoelectric gluon distribution function of the QGP. The factorization originates from the factorization of the total density matrix. The structures of the Wilson line in the chromoelectric gluon distribution function are different for inclusive and differential reaction rates. The relation between the two is similar to that between the gluon PDF and the gluon TMDPDF, except that the time axis here is along the real time rather than the lightcone time. In the recombination, we also made semiclassical approximations. One semiclassical approximation assumes the initial state density matrix is diagonal in color space, while the second semiclassical approximation keeps only the lowest term in the gradient expansion. Finally we worked out the leading quantum correction to the semiclassical Boltzmann transport equation by computing the next-to-leading term in the gradient expansion. The factorization in the transport equation allows us to use experimental measurements on quarkonium suppression in heavy ion collisions to probe the chromoelectric gluon distribution functions of the QGP. The chromoelectric gluon distribution functions are defined nonperturbatively here so in principle, they can be computed by using lattice QCD or the AdS/CFT correspondence. It would be interesting to investigate how much perturbative and nonperturbative calculations differ for the distribution function. Once the distribution function is determined nonperturbatively, one can combine it with the thermal width of quarkonium extracted from lattice QCD calculations to constrain the quarkonium in-medium real potential indirectly. In practice, one may choose Υ(1S) at low temperature for a well-controlled power counting. Furthermore, the differential reaction rates depend on a new momentum dependent chromoelectric distribution function defined by two electric fields connected via staple-shaped Wilson lines. It will be interesting to explore what other physics processes are sensitive to this correlation in the QGP. Finally, one can implement the quantum correction to the semiclassical transport equation in phenomenological studies and investigate the importance of quantum corrections. The framework developed here can be easily generalized to study quarkonium transport in cold nuclear matter by replacing the thermal QGP density matrix with a density matrix describing cold nuclear matter, which is relevant for quarkonium production in eA collisions in the future Electron-Ion Collider.

A Gauge Link at Infinite Time
In the main text, the Wilson lines along the time axis arise due to field redefinitions. Here we show how to obtain the Wilson line along the spatial direction at infinite time. The derivation is similar to that of Ref. [82] where the gauge link at infinite lightcone time was derived for the TMDPDF. For our purposes here, we first discuss modes that contribute to this gauge link.
The effective theory pNRQCD is based on separation of scales: M M v M v 2 T . We will distinguish modes by their momentum scaling p µ = (p 0 , p 1 , p 2 , p 3 ). The hard, soft and ultrasoft modes scale as In the derivation of pNRQCD, the hard and soft modes are integrated out, so this theory describes the dynamics of the ultrasoft modes. The Wilson line along the time axis discussed in the main text resums the interaction between the c.m. motion of the octet and the ultrasoft A 0 field. In fact, another mode with a scaling that falls between the soft and ultrasoft modes can influence the ultrasoft dynamics. It has the momentum scaling and thus is named the Coulomb mode. The Coulomb modes that mediate the interaction between the QQ pair has been included in the pNRQCD Lagrangian (2.1) as potentials. However, the Coulomb mode that couples the c.m. motion of the octet with the medium has not been considered yet. As we will show below, it is this Coulomb mode that leads to the gauge link at infinite time. The relevant part of the Lagrangian density is where the new term compared with Eq. (2.1) is the c.m. kinetic energy D 2 R /(4M ). For ultrasoft gauge fields, the c.m. kinetic energy can be omitted since it is at higher order in the v expansion. However, for Coulomb gauge fields, the c.m. kinetic energy is at the same order as iD 0 , which is at leading order. In the adjoint representation, the c.m. kinetic energy can be written as We will first focus on the interaction given by the term linear in A and discuss the term quadratic in A later. For the interaction linear in A, we will resum a series of Feynman diagrams. The n-th diagram is depicted in Fig. 2. The outgoing octet field is on-shell and ultrasoft: p ∼ M (v 2 , v 2 , v 2 , v 2 ). All gluon fields interacting with the octet are Coulomb: q j ∼ M (v 2 , v, v, v), so q j + p = q j + O(M v 2 ) for the spatial components of the momentum. The amplitude can be written as (we want to keep the A fields explicitly for later convenience) Here R 0 is the starting position of the octet field sitting on the left end of Fig. 2. We set R 0 = 0 for simplicity for the moment and discuss the case with R 0 = 0 later. (We can also keep t 0 in the phase where t 0 is the starting time of the octet field. But t 0 will become irrelevant in the following derivation.) Shifting q n → q n + p and applying the trick used in Ref. [82], we find for the first propagator and a similar expression for a generic propagator Then the integral over q 0 m can be thought of as a Fourier transform, with the conjugated time 4M m j=1 λ j , since the gluon lines are incoming with the phase defined by e −iq·x . In the large mass limit, the conjugated time 4M m j=1 λ j → +∞, since the λ j 's are positive. The octet field is outgoing in Fig. 2 so the diagram corresponds to final-state interactions in dissociation. If the octet field is incoming, as in recombination, a different sign will be obtained in the conjugated time, which leads to −∞ in the large mass limit. After the Fourier transform in the large mass limit, the expression (A.7) becomes M n = g n n j=1 d 3 q j (2π) 3 q n · A(t = +∞, q n ) (q n ) 2 − i (2q n + q n−1 ) · A(t = +∞, q n−1 ) (q n + q n−1 ) 2 − i × · · · × (2 n j=2 q j + q 1 ) · A(t = +∞, q 1 ) ( n j=1 q j ) 2 − i d 3 q j (2π) 3 d 3 R j e −iq j ·R j q n · A(t = +∞, R n ) (q n ) 2 − i × 2(q n + q n−1 ) · A(t = +∞, R n−1 ) + i∇ n−1 · A(t = +∞, R n−1 ) (q n + q n−1 ) 2 − i × · · · × 2( n j=1 q j ) · A(t = +∞, R 1 ) + i∇ 1 · A(t = +∞, R 1 ) ( n j=1 q j ) 2 − i . (A.12) One can show n j=1 e −iq j ·R j = e −iq n ·(Rn−R n−1 ) e −i(q n +q n−1 )·(R n−1 −R n−2 ) × · · · × e −i( n j=2 q j )·(R 2 −R 1 ) e −i( n j=1 q j )·R 1 . (A.13) Then we can do a change of variables, k m = n j=m q j to simplify the momentum integrals: M n = g n n j=1 d 3 R j A(t = +∞, R n ) · i∇ n d 3 k n (2π) 3 e −ikn·(Rn−R n−1 ) (A.14) Using the standard integral in the derivation of Coulomb potential from the propagator, we find M n = (ig) n n j=1 d 3 R j A(t = +∞, R n ) · ∇ n 1 4π|R n − R n−1 | × 2A(t = +∞, R n−1 ) · ∇ n−1 1 4π|R n−1 − R n−2 | + 1 4π|R n−1 − R n−2 | ∇ n−1 · A(t = +∞, R n−1 ) × · · · × 2A(t = +∞, R 1 ) · ∇ 1 1 4π|R 1 | + 1 4π|R 1 | ∇ 1 · A(t = +∞, R 1 ) .
In this derivation we chose R 0 = 0 and thus the Wilson line starts at R 0 = 0. For a nonvanishing R 0 , we need to restore the factor exp i n j=1 q j · R 0 , (A. 26) in (A.7). Then the last delta function will become δ 3 (R 1 − R 0 ) and the Wilson line will start at R 0 . The gauge link in the amplitude starts from the spatial position of the octet field and extend to spatial infinity. In the complex conjugate of the amplitude, the gauge link comes from the spatial infinity and stops at the spatial position of the other octet field. This leads to the Wilson lines at infinite time shown in Fig. 1.