Quarkonium in-Medium Transport Equation Derived from First Principles

We use the open quantum system formalism to study the dynamical in-medium evolution of quarkonium. The system of quarkonium is described by potential non-relativistic QCD while the environment is a weakly coupled quark-gluon plasma in local thermal equilibrium below the melting temperature of the quarkonium. Under the Markovian approximation, it is shown that the Lindblad equation leads to a Boltzmann transport equation if a Wigner transform is applied to the system density matrix. Our derivation illuminates how the microscopic time-reversibility of QCD is consistent with the time-irreversible in-medium evolution of quarkonium states. Static screening, dissociation and recombination of quarkonium are treated in the same theoretical framework. In addition, quarkonium annihilation is included in a similar way, although the effect is negligible for the phenomenology of the current heavy ion collision experiments. The methods used here can be extended to study quarkonium dynamical evolution inside a strongly coupled QGP, a hot medium out of equilibrium or cold nuclear matter, which is important to study quarkonium production at heavy ion, proton-ion, and electron-ion collisions.


I. INTRODUCTION
Heavy quarkonium production in hadron colliders has been studied extensively in both theory and experiment. In proton-proton collisions, the production process factorizes into a short-distance process of producing a heavy quark antiquark pair and a long-distance coalescence into a bound state [1]. In heavy ion collisions, the production process is complicated by the existence of a hot nuclear environment, the quark-gluon plasma (QGP). By comparing the quarkonium production at proton-proton and heavy ion collisions, one can study the properties of the hot medium produced during the collision (with the modification of the initial hard production due to heavy nuclei properly included). Static screening has been studied since the pioneering work of Ref. [2], which provides a partial understanding of the the suppression of quarkonia in heavy ion collisions. For a more complete understanding, a theoretical description of quarkonium dynamics that also accounts the dynamical screening and recombination inside the hot nuclear medium is needed.
There have been several approaches to address the question. First, statistical hadronization models have been used to describe charmonium production [3,4]. In these models it is assumed that the charm quark evolves unbound inside the hot medium due to the Debye screening. During the evolution, the charm quark equilibrates kinematically but not chemically, because the total number of charm quarks is fixed by the initial hard scattering.
The annihilation of charm quarks is negligible during the lifetime of the QGP. Thermal production is also negligible because of the large quark mass, compared with the medium temperature. Charmonium is produced from coalescence of charm quarks and anti-quarks with thermal momenta at the transition hyper-surface of QGP to a hadron gas. Although the model has some phenomenological success, it is limited to the study charmonium with low transverse momentum. The kinematic thermalization assumption is never justified for charmonium at large transverse momentum and bottomonium.
Debye screening of the potential is also accounted for when solving the bound state wavefunction. In many studies, the dissociation rate is calculated from perturbative QCD while the recombination is modeled from detailed balance with an extra suppression factor accounting for the incomplete thermalization of heavy quarks. The recombination process has also been analyzed in the framework of perturbative QCD with parametrized non-thermal heavy quark momentum distributions [19]. Many studies have used potential non-relativistic QCD (pNRQCD) to study quarkonium dissociation rates inside the QGP [20][21][22]. Recombination in a pNRQCD-based Boltzmann equation has been studied in Ref. [23]. New studies construct coupled Boltzmann transport equations of both heavy quarks and quarkonia, in which the heavy quark momentum distribution is not from an assumed parametrization but rather calculated from real-time dynamics, and quarkonium dissociation and recombination are calculated in the same theoretical framework [23][24][25]. By using the coupled Boltzmann transport equations, detailed balance and thermalization of heavy quark and quarkonium can be demonstrated from the real-time dynamics of heavy quark energy loss and the interplay between quarkonium dissociation and recombination.
More recently, an approach based on open quantum systems has been studied widely [26][27][28][29][30][31][32][33][34][35]. In this approach, the system of heavy quark and quarkonium and the medium evolve unitarily together. When the environment degrees of freedom are traced out, the system evolves non-unitarily and stochastic interactions can appear. This approach is a quantum description rather than a semi-classical equation. It has the advantage that non-unitarity appears automatically after tracing out the environment, while at the same time, preserving the total number of heavy quarks (by preserving the trace of the system density matrix).
Quarkonium dissociation occurs during the non-unitary evolution but the unbound heavy quark pairs from dissociation never disappear from the system and they may recombine. This feature is never easily realized in transport models based on complex potentials. Another advantage is that the recombination effect is included systematically in this procedure.
Meanwhile, the non-unitary time evolution is generally irreversible. For a general discussion of occurrence of time-irreversible processes from time-reversible underlying theory, we refer to Ref. [36]. The combination of the open quantum system and effective field theory (EFT) has also been recently used to study different physical systems: dissipative fluids [37], deep inelastic reactions [38] and bottomonium suppression in Au-Au collisions [34].
In this paper, we demonstrate a deep connection between approaches of open quantum systems and transport equations. More specifically, we use the open quantum system formalism, EFT of QCD, and the Wigner transform to derive the Boltzmann transport equation.
Our derivation clarifies the conditions for the validity of quarkonium transport (rate) equations that are based on Boltzmann transport equations. We will justify the Markovian approximation in the open quantum system approach and the molecular chaos approximation in the Boltzmann equation. The work of Ref. [34] focuses on quantum evolution of the density matrix and neglects center-of-mass (c.m.) motions of heavy quark anti-quark pairs, thus is unable to study observables as functions of transverse momentum and rapidity of the quarkonium. In this work, we explicitly keep track of the c.m. motion and focus on deriving the semi-classical Boltzmann transport equation from the quantum evolution of the system density matrix. This paper is organized as follows: first, the open quantum system and the quantum master equation, Lindblad equation, are briefly reviewed in Section II. The Boltzmann transport equation is derived in Section III. Then quarkonium annihilation is studied similarly in Section IV. Finally, conclusions are drawn in Section V.

II. LINDBLAD EQUATION IN WEAKLY-COUPLED SYSTEM
In this section we briefly review standard results in open quantum systems, which are covered in many textbooks, see, for example, Ref. [39]. Assume the Hamiltonian of the system and environment (thermal bath) is given by where H S is the system Hamiltonian, H B is the environment Hamiltonian, and H I contains the interactions between system and environment. The interaction Hamiltonian is assumed to be factorized as follows: where α denotes all quantum numbers.
separately. Here ρ B 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 and T is the time-ordering operator. We assume the interaction is a weak perturbation and expand the evolution operator to the second order in H I .
We shall assume the initial condition is given by where the environment density matrix is assumed to be time-independent. We define Then by taking the partial trace over the environment we can obtain the evolution equation Using O (B) α = 0 and inserting complete sets of the system we obtain Finally defining the Lindblad operator L ab ≡ |a b| and we obtain the Lindblad equation up to second order in perturbation theory The relation θ(t) = (1 + sign(t))/2 has been used in the derivation. It will be shown in the next section that for quarkonium, the commutator term is a loop correction of the real part of the Hamiltonian. The anticommutator term describes the dissociation of quarkonium, which can also be thought of as an imaginary part of the potential. The second term on the right hand side of Eq. (12) represents the recombination contribution. A direct conclusion from Eq. (12) is the conservation of probability: Trρ S (t) = Trρ S (0). This implies the unbound heavy quark antiquark pair from quarkonium dissociation stays as active degrees of freedom of the system and may recombine later in the evolution.
The form of the Lindblad equation is valid up to all orders in the perturbative expansion [40]. So the higher-order terms neglected here can also be written in the form of the Lindblad equation. The Lindblad equation cannot be written in the form of a von Neumann equation because the evolution is non-unitary. The time-irreversibility can be seen by noting that the relative entropy of the system with respect to a steady state under the partial trace is monotonically decreasing [40]. The partial trace over the environment can be thought of as an average over different environment configurations. Though the dynamics involving each configuration is governed by a time-reversible theory with a unitary evolution, after average, the dynamics becomes time-irreversible and non-unitary.

III. DERIVATION OF BOLTZMANN EQUATION
In this section we will derive the Boltzmann transport equation by applying the Lindblad equation, Eq. (12), to the Wigner transform of the density matrix describing heavy quark antiquark pairs that can be bound or unbound. The system in vacuum can be described by pNRQCD [42,45]. The effective theory can be constructed from QCD by a non-relativistic expansion assuming the separation of scales: where M is the heavy quark mass and v is the velocity of the heavy quark antiquark inside a quarkonium. The quarkonium size is roughly given by r ∼ 1/(M v). The environment is a weakly-coupled QGP in local thermal equilibrium, ρ B = 1 Z e −βH B , where Z = Tr B e −βH B . Thus the correlations in Eq. (7) can be calculated in real-time thermal field theory. We review different definitions of thermal correlations (Green's functions) in Appendix A. We will use free thermal Green's functions of gauge fields. Our derivation can be extended by using resummed thermal propagators. Resummed thermal propagators and pNRQCD have been used to investigate static heavy quark antiquark pairs at finite temperature [46]. The plasma provides two extra scales: the temperature T and the Debye mass m D (we use units in which k B = 1). Here we will focus on the case where quarkonium exists as a well-defined bound state in a QGP that is below the melting temperature of the quarkonium, so We do not consider cases with MeV for charmonium and bottomonium, and the temperatures realized in current heavy ion experiments are smaller than this. For our choice of scaling both dissociation and recombination are possible.
PNRQCD can be constructed by matching with NRQCD at the scale M v. The matching can be done perturbatively if M v Λ QCD or non-perturbatively. In either case, the quarkonium interacts with gluons from the QGP via a dipole interaction at lowest order.
As will be seen below, the dipole interaction scales as rT ∼ T M v v, which is small in the assumed separation of scales. We assume a perturbative matching throughout the paper.
The dipole interaction is not running at one loop level [43,44], which means the coupling constant in the dipole term is set at the scale of M v, no matter of the scale of the scattering.
To make calculations easier here, we follow Ref. [45] and use a slightly different notation for the pNRQCD Lagrangian density.
The degrees of freedom, in the standard pNRQCD Lagrangian, are the color singlet S(R, r, t) and octet O a (R, r, t) where R and r are the center-of-mass (c.m.) and relative positions of the heavy quark antiquark pair. Here we define the "bra-ket" notation via for any function f of r. We use the "bra-ket" notation so that we no longer need to write the integral over r explicitly, which simplifies notations in the derivation. Summations over color indexes are assumed and higher order terms in the velocity expansion are neglected.
for later use. The covariant derivative on the octet field has been written out explicitly We will work up to the leading order (LO) in v 2 and s,o ∼ M v 2 by the virial theorem. When the medium is static in the rest frame of quarkonium, the quarkonium exchanges gluons with the medium with energy T . In our power counting, T ∼ M v 2 and hence the c.m. kinetic energy, p 2 cm 4M , is O(M v 4 ) and is therefore neglected. If the medium moves with respect to the quarkonium at a velocity v med , the gluon energy is boosted by We assume v med = 0 in the following but generalization to v med = 0 can be easily done by boosting the gluon distribution function.
We do keep track of the c.m. momentum so that momentum is conserved. The singlet and octet composite fields are given by where E is the eigenenergy of a state in the whole Hilbert space. The whole Hilbert space factorizes into two part: one part for the c.m. motion and the other for the relative motion. The operators a p rel (p cm ) and c 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 state, p rel for unbound singlet state and color a and p rel for unbound octet state. When we compute the square of matrix elements, we will average over the polarizations of non-S wave quarkonium states. In our notation, we omit the quantum number m 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 , |ψ p rel and |Ψ p rel . 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 separately with higher order terms in v neglected. Here E nl is the binding energy of the bound state |ψ nl . The annihilation and creation operators in the Fock space satisfy the following commutation relations: and all other commutators vanish.
The interaction part of the Hamiltonian of the theory is given in Eq. (13)  As discussed above, this is true in both perturbative and non-perturbative constructions of the pNRQCD. When rT ∼ 1, static screening effect of the potential is too strong to support the quarkonium bound state.
To use the Lindblad equation derived in Section II, we write The sum over α means The complete set used to construct the Lindblad operators are where 1 denotes the singlet while a is the color index of an octet. The unbound singlet state will not be used in our current calculation because at the order we are working, an unbound singlet cannot form a bound singlet by radiating out one gluon, only unbound octet can do so.
We are interested in the bound state evolution. Therefore our basic strategy is to study the time evolution of k 1 , n 1 l 1 , 1|ρ S (t)|k 2 , n 2 l 2 , 1 by sandwiching Eq. (12) between k 1 , n 1 l 1 , 1| and |k 2 , n 2 l 2 , 1 . To obtain the evolution equation of the semi-classical phase space distribution function, we will take the Wigner transform of the density matrix We will extract the linear dependence on t of γ ab,cd and σ ab terms in Eq. (12) The part inside the curly bracket gives the loop correction of the potential, which can be calculated as usual by the standard quantum field theory perturbative technique: computing the loop shown in Fig. 1 where other terms in the Lindblad Eq. (12) we obtain Here if we restore the c.m. kinetic energy, we can write where the c.m. velocity of the quarkonium is defined as v = k 2M . Now we proceed to compute the contributions from the other two terms in the Lindblad Eq. (12) omitted in Eq. (27). The − a,b,c,d The dissociation rate Γ disso of a quarkonium with momentum k and position x can be defined k,t) . The dissociation rate derived here is the same as calculated in Ref. [21] by taking the imaginary part of Fig. 1. The same dissociation term has been used in the Boltzmann transport equation in Ref. [23].
The a,b,c,d γ ab,cd L ab ρ S (0)L † cd term gives t a,i where f QQ (x, p cm , r, p rel , a, t = 0) is the two-particle distribution function of a heavy quark antiquark pair in color octet a with the c.m. position x and momentum p cm and relative position r and momentum p rel . Unlike in the dissociation term, one of the integrals over the wavefunctions of the relative motion involves the two-particle distribution function of QQ.
Now putting Eqs. (27), (28) and (29) together we finally infer the Boltzmann transport where the dissociation C nl (x, k, t) terms are defined in Eqs. (28) and (29). Both terms C Taking T ∼ M v 2 as previously assumed we find we must have p rel M v/(v 2 α 2 s ) which is clearly satisfied for p rel ∼ M v. For p rel large enough that this condition is not satisfied, the contribution to the integral involving ψ nl |r|Ψ p rel is negligible for such large p rel . This agrees with the intuition: a heavy quark antiquark pair with large relative momentum cannot form a bound state. So we can take f QQ (x, p cm , r, p rel , a, t) out of the wavefunction integral with the awareness that the contribution from r a B should vanish.
Furthermore, we make the molecular chaos assumption and write where x, r, p cm , p rel are the c.m. and relative positions and momenta of the heavy quark antiquark pair with positions x 1 , x 2 and momenta p 1 , p 2 . The factor 1 9 accounts the probability of the color state of QQ being in a specific octet state a. The molecular chaos assumption is valid when the rate of decorrelation between the heavy quark and antiquark is much larger than the relaxation rate of the system. The former is given by D −1 ∼ α 2 s T with α s at the scale T or m D while the later has been estimated above and is ∼ α s v 2 T with α s at the scale M v. In NRQCD, v ∼ α s (M v), so the molecular chaos assumption is valid.
Combining these two assumptions gives where the sum over color index a has been carried out. This is the recombination term used in the Boltzmann equation in Ref. [23]. In order to take the spin multiplicity into account, one must further insert a factor g s into C open effective field theory in order to conserve probability can be found in Ref. [38]. The annihilation is too slow to be of much interest for phenomenology but we study it as an interesting example of how Lindblad-type operators enter the time evolution equation for the density matrix. In our case, we first restore the standard pNRQCD notation of singlet field, S(R, r, t) ≡ r|S(R, t) , i.e., we project the wavefunction of the relative motion onto the relative position space. Then we can add two new terms in the density matrix evolution equation where the evolution term is explicitly trace-preserving.
As above, we are interested in the bound state and will sandwich the density matrix between two bound quarkonium states and then do a Wigner transform.
where we have used the Markovian approximation and written the delta function in energy as t. In the summation over n 3 and l 3 , only n 3 = n, l 3 = l contributes due to the delta function in energy (we assume no degeneracy in the bound state eigenenergy beyond that implied by rotational invariance).
We can define the annihilation rate of a quarkonium state nl, which is non-zero even in vaccum and should be distinguished from the dissociation rate inside QGP. For S-wave, one may set Γ(r) = Γδ 3 (r) and then Γ S = Γ|ψ S (0)| 2 , i.e., the annihilation rate depends on the wavefunction of the relative motion at the origin.
The other term in the anticommutator will give the same result. Under a Wigner transform, these two terms lead to in the Boltzmann equations Typically Γ nl ∼ 10 keV, so for a QGP with a lifetime ∼ 10 fm ∼ 0.05 MeV −1 , the effect from quarkonium annihilations is negligible on the in-medium evolution. It is justified to assume that the total number of heavy quarks is conserved during the in-medium evolution.
where we have inserted an identity d 3 k 3 δ 3 (k 3 ) = 1 and written δ 3 (k 3 ) as a spatial integral over x . It should be noted that the integral over k 3 is already a Wigner transform on the density matrix of the second particle with momentum k 3 and position x . If we further apply a Wigner transform on the density matrix of the first particle and properly reshuffle labels, we obtain the contribution of this term in the Boltzmann equation It involves the two-particle distribution function f (x, k, nl; x , k , n l ; t) of two quarkonium states nl and n l with positions x, x and momenta k, k respectively. When the second quarkonium with the quantum number n l annihilates, it leads to an increase in the oneparticle distribution function of quarkonium with quantum number nl. Therefore, this term together with the term in Eq. (37) guarantees the conservation of probability in the oneparticle distribution of quarkonium. However, as mentioned earlier, the annihilation effect is negligible in current heavy ion collision experiments.

V. CONCLUSION
In this paper, we used the open quantum system formalism where the system of heavy quarks and quarkonium is described by pNRQCD at LO in the non-relativistic expansion while the environment is a weakly-coupled thermal QGP. We derived the Boltzmann trans- The derivation here provides a theoretical justification of quarkonium transport equations inside a weakly-coupled QGP below the melting temperature of the quarkonium. It connects two main approaches of the phenomenology of quarkonium production in heavy ion collisions. One can improve the derivation by working to the next-leading-order in the coupling constant and expansion parameters in both the system sector (pNRQCD) and the environment sector (thermal QCD). In the case of a non-perturbative construction of pNRQCD, a similar derivation is possible. The connection between the complex potential calculated on lattice [52] and the transport equation is worth exploring in our framework. The derivation can be extended to the case of quarkonium evolution inside a strongly-coupled QGP, a hot medium out of equilibrium or cold nuclear matter by replacing the Green's functions of thermal QCD with those in the corresponding media. It would also be interesting to study the viscous and anisotropic corrections to the Debye screening, the dissociation rate of quarkonium [53][54][55] and the recombination in a non-thermal QGP. The effect of a turbulent plasma on the heavy quark antiquark pair or quarkonium in the early stage of heavy ion collisions can also be explored [56]. A description of the quarkonium evolution through cold nuclear matter will be useful to studies of quarkonium production at both proton-ion and electron-ion colliders. In real-time thermal field theory, the commonly used Green's functions are the ">", "<", retarded, advanced and time-ordered Green's functions. They are defined as follows: The Fourier transform is defined as where X could be >, <, R, A or T . In momentum space for a free theory To compute a,b σ ab L ab , we first note that where in the last line we flipped α ↔ β, t 1 ↔ t 2 and |a ↔ |b . We can split sign(t 1 − t 2 ) into θ(t 1 − t 2 ) and −θ(t 2 − t 1 ) in a,b σ ab L ab and just need to compute the θ(t 1 − t 2 ) term.
The −θ(t 2 − t 1 ) term is given by the Hermitian conjugate of the θ(t 1 − t 2 ) term. Therefore a,b σ ab L ab is Hermitian and this term can be thought of as a correction to the system Hamiltonian. To carry out the calculation explicitly, we first write Then we can replace C αβ (t 1 , t 2 ) in Eq. (B1), due to the θ(t 1 − t 2 ), with The term inside the square brackets is the time-ordered thermal propagator in momentum space.
We are interested in the bound state part of the density matrix, so we set |a = |k 1 , n 1 l 1 , 1 and |b = |k 2 , n 2 l 2 , 1 . Then we can compute where E p = p 2 rel M . This can be written as It should be noted that p 0 cm here does not represent the c.m. energy of the octet. In fact, it is the total energy of the composite octet particle, p 0 cm = p 2 cm 4M + . To simplify the expression, we make the Markovian approximation t → ∞. Then integrating over t 1 and t 2 will give two δ-functions in energy. Plugging Eqs. (B3) and (B5) into a,b σ ab L ab and integrating over t 1 , t 2 , R 1 and R 2 we find The two time integrals give a product of two delta functions in energy where ω i = E k 1 − p 0 cm − q 0 for i = 1, 2 and ω 1 = ω 2 = ω. We interpret one factor of 2πδ(ω) to be the time interval, so the double time integral is interpreted as follows See Ref. [57] for details. This argument also applies in the next two sections C and D. The δ-functions in energy and momentum give k 1 = k 2 = k, n 1 = n 2 = n and l 1 = l 2 = l (we assume no degeneracy in the bound state eigenenergy beyond that implied by rotational invariance). So we have We will show this term gives the dissociation term in the Boltzmann equation. We first compute the term − 1 2 γ ab,cd L † cd L ab ρ S (0) . As explained previously, we are interested in the bound state part of the density matrix and take k 1 , n 1 l 1 , 1|ρ S (t)|k 2 , n 2 l 2 , 1 , so we set |d = |k 1 , n 1 l 1 , 1 . Since at lowest order of the expansion, the transition between bound state and unbound pair only occurs via singlet-octet transition, we have |a = |c = |p cm , p rel , a 1 and double summations over |a and |c become just one summation. Similarly we have |b = |k 3 , n 3 l 3 , 1 . We need to compute We can start with The correlation needed is where we have used the expression of D < ab µν (q) in Appendix A. It should be pointed out that one can also use the expression of D > ab µν (q) and will obtain the same result due to the relation 1+n B (q 0 )+n B (−q 0 ) = 0. Now we can write the term k 1 , n 1 l 1 , 1|γ ab,cd L † cd L ab ρ S (0)|k 2 , n 2 l 2 , 1 out explicitly as Integrating over R 1 and R 2 gives two delta functions in momenta, δ 3 (k 1 − p cm + q)δ 3 (k 3 − p cm + q). Under the Markovian approximation, t → ∞, integrating over t 1 and t 2 will give another two delta functions δ(E k 1 −E p +q 0 )δ(E k 3 −E p +q 0 ). Since E k i < 0 and E p > 0, some energy has to be transferred to the bound state to break it up to an unbound state, and thus q 0 has to be positive. (If we use D > ab µν (q), here we would have δ(E k 1 −E p −q 0 )δ(E k 3 −E p −q 0 ) and q 0 is negative.) Integrating over R 1 , R 2 , t 1 , t 2 and k 3 gives n 3 ,l 3 ,a 1 ,i 1 n B (q)(2π) 5 δ(E k 1 − E p + q)δ(E k 3 − E p + q)δ 3 (k 1 − p cm + q) 2T F 3N C q 2 g 2 ψ n 1 l 1 |r i 1 |Ψ p rel Ψ p rel |r i 1 |ψ n 3 l 3 k 1 , n 1 l 1 , 1|ρ S (0)|k 2 , n 2 l 2 , 1 , where we have used for any smooth function f (q) Due to the energy δ-functions, the sum over n 3 and l 3 gives n 3 = n 1 , l 3 = l 1 (we assume no degeneracy in the bound state eigenenergy beyond that implied by rotational invariance).
Under a Wigner transform of the form Eq. (21) (where we set k 1 = k + k 2 , k 2 = k − k 2 , n 1 = n 2 = n and l 1 = l 2 = l and then do a shift in c.m. momentum p cm → p cm + k 2 ), Eq. (C7) finally leads to The other term in the anti-commutator gives the same result. So applying the Wigner transform to the − 1 2 γ ab,cd k 1 , n 1 l 1 , 1|{L † cd L ab , ρ S (0)}|k 2 , n 2 l 2 , 1 term in the Lindblad equation, yields the negative of Eq. (C8).