Quantum and classical dynamics of heavy quarks in a quark-gluon plasma

We derive equations for the time evolution of the reduced density matrix of a collection of heavy quarks and antiquarks immersed in a quark gluon plasma. These equations, in their original form, rely on two approximations: the weak coupling between the heavy quarks and the plasma, the fast response of the plasma to the perturbation caused by the heavy quarks. An additional semi-classical approximation is performed. This allows us to recover results previously obtained for the abelian plasma using the influence functional formalism. In the case of QCD, specific features of the color dynamics make the implementation of the semi-classical approximation more involved. We explore two approximate strategies to solve numerically the resulting equations in the case of a quark-antiquark pair. One involves Langevin equations with additional random color forces, the other treats the transition between the singlet and octet color configurations as collisions in a Boltzmann equation which can be solved with Monte Carlo techniques.


JHEP06(2018)034
Heavy quarkonia, bound states of charm or bottom quarks, constitute a prominent probe of the quark-gluon plasma produced in ultra-relativistic heavy ion collisions, and are the object of many investigations, both theoretically and experimentally. Recent data from the LHC provide evidence for a sequential suppresion, with the most fragile (less bound) states being more strongly suppressed [1], while there are indications that charm quarks are sufficiently numerous to recombine, an effect that is seen to counterbalance the expected suppression [2]. These findings are in line with general expectations. The dissociation of quarkonium was suggested long ago [3] as resulting from the screening of the binding forces by the quark gluon plasma. Recombination is a natural phenomenon to expect [4,5] whenever the number of heavy quarks is sufficiently large, which seems to be the case of charm quarks at the LHC. However, in order to go beyond these qualitative remarks, and extract precise information about the dynamics, we have to address a rather complicated many-body problem. Even leaving aside the production mechanisms of heavy quarkonia in hadronic collisions, which are not fully understood yet (see e.g. [6]), the description of the interactions of heavy quarks with an expanding quark-gluon plasma is indeed complicated for a number of reasons. Many effects can contribute, among which: screening affecting the binding potential, collisions with the plasma constituents, absorption of gluons of the plasma by the bound states. It should be added to this that the bound states do not exist as objects "deposited" in the plasma: it takes time before a newly created quark-antiquark pair can be considered as a bound state, and during this time it is interacting with the plasma. This is a feature that is often forgotten in many models that attempts to describe the data (for representatives of recent phenomenological analyses, see e.g. [7][8][9] and references therein). Thus, most models emphasize static or stationary aspects (even when the expansion of the plasma is taken into account): this is the case of potential models [10], spectral function calculations [11], or kinetic approaches based on rate equations [9,12]. Clearly a fully time-dependent, out of equilibrium treatment is called for. Such a treatment should also establish contact between the dynamics of heavy quark-antiquark pairs and their bound states, and that of isolated heavy quarks in a quark-gluon plasma (see e.g. [13,14]). In short, there is a need for a general, simple, and robust formalism, where all the relevant effects can be treated within the same framework. In this respect, the observation that the collisions could be taken into account by an imaginary potential is a significant one [15][16][17][18][19].
In recent years it has been recognized that techniques from the theory of open quantum systems (see e.g. [20,21]) could offer a fruitful perspective on this problem. A system of heavy quarks in a quark gluon plasma falls indeed in the category of typical problems JHEP06(2018)034 addressed by this theory: a small system, weakly interacting with a large "reservoir", the quark-gluon plasma. This point of view has emerged explicitly or implicitly in a number of recent works: derivation of a master equation, and corresponding rate equations [22], use of the influence functional method [23,24], solution of a stochastic Schrödinger equation [25][26][27], or direct computation of the evolution of the density matrix [28,29].
The present work follows similar lines. Its initial motivation was to generalize the results of [24] to the non-Abelian case (QCD). Part of that generalization is straightforward, and relies on the same approximations as those used in the case of the abelian plasma. To some extent, this program has already been considered in the recent work by Akamatsu [30], albeit using a formalism slightly different from that used in [24] and in the present paper. However, color degrees of freedom modifies the picture in a very substantial way. The reason is that, in a collision involving one gluon exchange for instance, color changes in a discrete way, in contrast to position or momentum which vary continuously. Thus, while we can treat the motion of the heavy quarks within a semi-classical approximation, there is no such semi-classical limit for the color dynamics (except perhaps in the large N c limit). It follows that the derivation of Fokker-Planck or Langevin equations made in the abelian case needs to be reconsidered, which we do in this paper. We shall see that the complete dynamics, including the color degrees of freedom, can still be described by Fokker-Plack and Langevin equations, but only in very specific circumstances.
This paper focusses on conceptual issues. It is organized as follows. In section 2 we derive the quantum master equation for the reduced density matrix of a system of heavy quarks and antiquarks immersed in a quark-gluon plasma, in thermal equilibrium. This equation, whose structure is close to that of a Lindblad equation, is used as a starting point of all later developments. In section 3 we rederive from it the results that we had previously obtained for the abelian plasma [24] using a path integral formalism. In particular we recover, after performing a semi-classical approximation, the Fokker-Planck and Langevin equations that describe the random walks of center of mass and relative coordinates of a quark-antiquark pair. This section on the abelian plasma paves the way for the treatment of the non abelian case discussed in section 4. The equations that we present there, before we do the semi-classical approximation, are fully quantum equations. But they are difficult to solve in general. Thus, in section 5 we look for additional approximations that allow us to obtain solutions in some particular regimes, in order to start getting insight into the general solution. In particular, we explore two ways of implementing the semi-classical approximation. In the first case, we restrict the dynamics to stay close to a maximum entropy color state, where the colors of the heavy quarks are random. In this case the dynamics is described by a Langevin equation with a new random color force. The method used in this case is easily extended to the case of an arbitrary number of quark-antiquark pairs, and allows us to address the question of recombination. However, it is based on a perturbative approach that breaks down for some values of the parameters. Another strategy focuses on the case of a single quark-antiquark pair. The transition between singlets and octets are treated as "collisions" in a kinetic equation that we solve using Monte Carlo techniques. The last section summarizes our main results, and presents a brief outlook. Several appendices at the end gather various technical material.

JHEP06(2018)034 2 Equation for the density matrix of heavy quarks in a quark-gluon plasma
Our description of the heavy quark dynamics in a quark-gluon plasma is based on the assumption that the interaction between the heavy quarks and the quark-gluon plasma is weak, and can be treated in perturbation theory (with appropriate resummations). The generic hamiltonian for such a system reads where H Q describes the dynamics of the heavy quarks in the absence of the plasma, H pl is the hamiltonian of the plasma in the absence of the heavy quarks, and H 1 is the interaction between the heavy quarks and the plasma constituents. The heavy quarks are treated as non relativistic particles, and the spin degree of freedom is ignored: the state of a heavy quark is then entirely specified by its position and color. As we have mentioned already, we shall consider H 1 to be small and treat it as a perturbation. In Coulomb gauge, and neglecting magnetic interactions, this interaction takes the form A a 0 (r)n a (r), (2.2) where n a denotes the color charge density of the heavy particles. For a quark-antiquark pair, for instance, this is given by 1 n a (x) = δ(x −r) t a ⊗ I − I ⊗ δ(x −r)t a , where we use the first quantization to describe the heavy quark and antiquark, and the two components of the tensor product refer respectively to the Hilbert spaces of the heavy quark (for the first component) and the heavy antiquark (for the second component). In eq. (2.3), t a is a color matrix in the fundamental representation of SU(3) and describes the coupling between the heavy quark and the gluon. The coupling of the heavy antiquark and the gluon is described by −t a , witht a the transpose of t a . We are looking for an effective theory for the heavy quark dynamics, obtained by eliminating the plasma degrees of freedom. In previous works, this was achieved explicitly by constructing the Feynman-Vernon influence functional [31], using the path integral formalism (see e.g. [24,30]). In the present paper, we shall proceed differently, by writing directly the equations of motion for the reduced density matrix of the heavy quarks. Although the derivations presented here are self-contained, we emphasize that the main approximations that are implemented in the present section are quite common in various fields, and belong to what is commonly referred to as the theory of open quantum systems (see e.g. [20]).
We assume that the system contains a fixed number, N Q , of heavy quarks (and, in general, an equal number of antiquarks). We call D the density matrix of the full system,

JHEP06(2018)034
and D Q the reduced density matrix for the heavy quarks. The latter is defined as the partial trace of the full density matrix over the plasma degrees of freedom, that is (2.4) In order to make contact with the work of ref. [24], we recall that a typical question addressed there was the following: given a set of heavy quarks at position X i at time t i , where X denotes collectively the set of coordinates of the quarks and antiquarks (temporarily ignoring color), what is the probability P (X f , t f |X i , t i ) to find them at position X f at time t f ? This probability is given by that is, it can be obtained as a specific element of the reduced density matrix. In [24] a representation of this quantity was obtained in terms of a path integral which is difficult to evaluate in general. 2 However, in the regime where a semi-classical approximation is valid, the dynamics that it describes is equivalent to that of a Fokker-Planck equation which can be easily solved numerically, in particular by solving the associated Langevin equation. Two approximations are involved in the construction of the influence functional such as presented in [24,30]. The first one is the weak coupling approximation for the interaction of the heavy quarks with the plasma, the second assumes that the response of the plasma to the perturbation caused by the heavy quarks is fast compared to the characteristic time scales of the heavy quark motion. An additional approximation, to which we refer to as a semi-classical approximation, leads, as we have just mentioned, to Fokker Planck and Langevin equations. The last two approximations exploit the fact that the mass of the heavy quark is large, i.e., M T . Thus, when the heavy quark is not too far from thermal equilibrium, its thermal wavelength λ th ∼ 1/ √ M T is small compared to the typical microscopic length scale ∼ 1/T . Under such condition, the density matrix can be considered as nearly diagonal (in position space), motivating a semi-classical approximation: indeed the off-diagonal matrix elements X|D Q |X die off when |X − X | λ th . The typical heavy quark velocity is of the order of the thermal velocity ∼ T /M 1, and the dynamics of the heavy fermions is much slower than that of the plasma. The typical frequency for the plasma dynamics is the plasma frequency which, for ultra-relativistic plasmas, is of the order of the Debye screening mass m D . During a time t ∼ m −1 D , the heavy quark moves a distance which is small compared to the size of the screening cloud, ∼ m −1 D . Thus, over a time scale characteristic of the plasma collective dynamics, the heavy quark positions are almost frozen (they are completely frozen in the limit M → ∞). One can also recognize that the collisions of the heavy particles with the light constituents of the plasma involve the

Equation for the density matrix
The density matrix obeys the general equation of motion In order to treat the interaction between the plasma and the heavy particles using perturbation theory, we move to the interaction representation. We set H = H 0 + H 1 , with H 0 = H Q + H pl and define where D I (t), the interaction representation of the density matrix, satisfies the equation Here, H 1 (t) denotes the interaction representation of H 1 . The evolution operator in the interaction representation, where the symbol T denotes time ordering. Similarly, eq. (2.8) can be integrated formally using the Schwinger-Keldysh contour [33,34]: where the operator T C orders the operators H 1 (t C ) along the contour parameterized by the contour time t C , with the operators carrying the largest t C coming before those with smaller t C (see figure 1). The upper branch of the contour, with time running from t 0 to t, represents the evolution operator U I (t, t 0 ), the lower branch of the contour, with time running from t to t 0 , represents the operator U † I (t, t 0 ). As can be seen in eq. (2.10), in the expansion of D I (t) in powers of H 1 (t), the operators H 1 (t) that sit on the left of D(t 0 ) live on the upper branch of the contour, while those that appear on the right of D(t 0 ) live on the lower branch (they come later along the contour).
To proceed further, we assume that, at the initial time t 0 , the density matrix factorizes

JHEP06(2018)034
We also assume that at time t 0 , the plasma is in a state of thermal equilibrium, so that its density matrix D I pl (t 0 ) = D pl (t 0 ) is a canonical density matrix, where β = 1/T , with T the equilibrium temperature. This factorization of the density matrix allows for a simple calculation of the trace over the plasma degrees of freedom. Let us then examine perturbation theory at second order in H 1 , with H 1 given by eq. (2.2). Performing the trace over the plasma degrees of freedom is immediate, thanks to the specific structure of H 1 and the factorization of the density matrix at t = t 0 . One obtains where, in the last three lines, we have used the convention that t 1 , t 1 run on the upper part of the contour, while t 2 , t 2 run on the lower branch. Note that the linear term vanishes since the plasma is color neutral (so that A a 0 (x) 0 = 0). Here the notation · · · 0 stands for the average with the plasma equilibrium density matrix, that is (2.14) Similarly the correlators of the gauge fields are diagonal in color, i.e. they are proportional to δ ab . These correlators are the exact correlators in the plasma (the fields are in the interaction representation, which corresponds to the Heisenberg representation when considering the plasma alone). They are written as The apparent inversion of the order of times in the last correlator results from the relation x) 0 which follows from the cyclic invariance of the trace.
It is convenient to represent the evolution of the density matrix by a diagram such as that in figure 2, where the upper and lower parts of the diagram may be associated to JHEP06(2018)034 Graphical representation of the evolution of the density matrix from t 0 to t. The horizontal lines represent the evolution operators U I (t, t 0 ) (upper branch) or U † I (t, t 0 ) (lower branch). When D Q is the density matrix of a single heavy quark, these horizontal lines may be interpreted as the associated propagators of the heavy particle. When D Q is the density matrix of a heavy quark-antiquark pair, a single horizontal line is replaced by a pair of lines, associated with the propagator of the pair (see figure 6 below). . These diagrams are in one-to-one correspondence with the terms in the last three lines of the right-hand side of eq. (2.13) for the single particle density matrix D I Q (t).
the corresponding upper and lower parts of the Schwinger-Keldysh contour. The explicit expression that this diagram represents is where α of β represent the set of quantum numbers that are necessary to specify the state of the heavy particles (essentially the position and color). The diagrammatic interpretation of eq. (2.13) is then given in figure 3 (for the case of a single particle density matrix).
In order to implement our further approximations, it is convenient to consider the time derivative of the density matrix. This can be obtained by taking the derivative of eq. (2.13) above (see eq. (B.1) in appendix B). But it is more instructive to return to eq. (2.8), and JHEP06(2018)034 rewrite it as ]. (2.17) This exact equation is obtained by formally integrating eq. (2.8) and inserting the solution back into the equation. Perturbation theory at second order in H 1 is recovered by replacing D I (t ) by D I (t 0 ) in the double commutator in the right hand side. One may then proceed to the average over the plasma degrees of freedom, as we did before, and get the following equation for the reduced density matrix D Q We have used the fact that the linear term vanishes in a neutral plasma, and the sum over the color index a is implicit. At this point, we can improve on strict perturbation theory. To do so we notice that the integral over t in eq. (2.17) is in fact limited to a region near t t: this is because ∆(t−t ) dies out when t − t m −1 D , and we are interested in the evolution of the density matrix over time scales that are much larger than m −1 D . Thus, noticing also that the difference D I (t) − D I (t ) involves in any case an extra power of H 1 , we replace D I (t ) by D I (t) in the right hand side of eq. (2.17), 3 turning the equation into an equation for D(t) which is now local in time. We shall furthermore exploit the fact that the density matrix approximately factorizes at all times, as does the density matrix at the initial time t 0 . Again, this is consistent with the weak coupling approximation since the correction to the factorized form necessarily involves additional powers of H 1 . The latter approximation allows us to perform the trace over the plasma degree of freedom, in the same way as we did earlier for strict perturbation theory. The resulting equation is in fact identical to eq. (2.18) in which we replace in the right hand side D I Q (t 0 )by D I Q (t). It is convenient for the following to write this equation in the Schrödinger picture. A simple calculation yields where we have set t − t = τ . This equation has the same physical content as the influence functional derived in [24], and it is based on analogous approximations. It relies on a weak coupling approximation, but goes beyond strict second order perturbation theory, in particular by implementing partial resummations. This equation contains in its right hand side a time integral that we shall simplify thanks to our last approximation: in line with the fact that only small values of τ are JHEP06(2018)034 relevant, it consists in replacing e −iH Q τ 1 − iH Q τ , and keep terms up to linear order in τ in the integrals. More precisely, we write and the time-dependence of n x (t) being given by the Heisenberg representation, n a x (t) = e iH Q t n a x e −iH Q t . We get At this point we use the values of the time integrals given in appendix A. These involve the zero frequency part of the time-order propagator ∆(ω = 0) = ∆ R (ω = 0, r) + i∆ < (ω = 0, r), which we identify with the real and imaginary part of a "potential". More precisely, we set This terminology stems from the fact that V (r) + iW (r) plays the role of a complex potential in a Schrödinger equation describing the relative motion of a quark-antiquark pair: the real part represents the screening corrections, and adds to the familiar interaction arising in leading order from one-gluon exchange, the imaginary part accounts effectively for the collisions between the heavy quarks and the plasma constituents [15,16]. After a simple calculation that uses the properties The first line of the right hand side of this equation describes a hamiltonian evolution, that is, it can be written as the commutator on the left hand side, with H Q replaced by 1 2 xx V (x − x )n a x n a x . It follows that we can shift the "direct", one-gluon exchange potential initially contained in H Q into V , and keep in H Q only the kinetic energy of the JHEP06(2018)034 heavy quarks. This is the strategy that was followed in [24] and that we shall adopt in this paper. In this way the potential V (r) becomes the screened Coulomb potential, and H Q represents only the kinetic energy of the heavy particles (see also the discussion after eq. (3.29) below).
The equation (2.24) is our main equation. It is a fully quantum mechanical equation. It is a Markovian equation for the reduced density matrix D Q (t). We shall write this equation in the following way with L = L 0 + L 1 + L 2 + L 3 , and The structure of eq. (2.25) is close to that of a Lindblad equation [35], but eq. (2.25) is not quite a Lindblad equation: while the operator L 2 can be put in the Lindblad form, this is not so of the operator L 3 , unless one adds extra, subleading terms (see the discussion in appendix B). For a recent discussion of the Lindblad equation for an abelian plasma, in a formalism not too different from the present one, see [29]). The notation is, at this stage, symbolic and just expresses the fact that the right hand side of eq. (2.24) is a linear functional of the density matrix. It will acquire a more precise meaning as we proceed. We may however make the following observation. When taking matrix elements between localized states, specified by the coordinates of the heavy particles, the density operators n x play the role of projection operators, and are diagonal in the coordinate representation. Thus the same matrix elements, as far as the coordinates are concerned, will appear on the left and the right. The operator L will then appear as a differential operator acting on this matrix element (in fact L 1 and L 2 are simply multiplicative, as we shall see).
It is convenient to associate a diagrammatic representation to the various contributions that we shall calculate. The relevant diagrams will preserve the topological structures of those already introduced, but because of the various approximations that we have performed, they cannot be calculated with standard rules. As an illustration, we display in figure 4 diagrams corresponding to the time derivative of the single particle density matrix (diagrams corresponding to the two particle density matrix are displayed in figure 6 below). All interactions in eq. (2.24) have become instantaneous. For this reason, we draw these as vertical gluon lines, or as tadpole insertions, located anywhere between t − τ and t. Note that terms where the two densities sit close together in eq. (2.26), like in [n a x n a x , D Q ], correspond to diagrams where the two ends of the gluon is hooked on the upper (or lower) part of the diagram, while a term such as n a x D Q n a x corresponds to a gluon joining the JHEP06(2018)034 . These diagrams illustrate generic processes taken into account in eq. (2.26), in the case of the single particle density matrix. Note that there is another diagram with a tadpole insertion, on the lower line, not drawn. Depending on the operator considered, the propagator of the gluon corresponds to V , W or involves spatial derivatives of W . Note that since we treat the heavy quarks and antiquarks as non relativistic particles, the direction of the arrows in such a diagram does not refer to the nature (quark or antiquark) of the heavy particle: rather, it is correlated to the SK contour time, forward (to the right) in the upper branch, backward (to the left) in the lower branch.
upper and lower parts of the diagram. Since, as we shall see, in QCD these two types of terms correspond also to different color structures, we shall find it convenient to split the operators L i into two components, x D Q n a x . Note that L 1 has only contributions of type a, i.e., L 1 = L 1a . In the rest of this paper, we shall deal only with the heavy quark reduced density matrix. We shall then drop the subscript Q in order to simplify the notation and write simply D in place of D Q .

Semi-classical approximation for abelian plasmas
The equation (2.24) is quite general. It holds for any system of heavy quarks and antiquarks. Depending on the system considered, the color density n a (x) and the density matrix D take different forms. In this section, we study the specific form of eq. (2.24), and the associated operators L i in eq. (2.26), for the single particle and the two particle density matrices, in the case of an electromagnetic (abelian) plasma. For simplicity, we shall continue to refer to the charged particles as quarks (positive charge) and antiquarks (negative charge). The interaction hamiltonian reads as in eq. (2.2), with n a (x) replaced by the density of charged particles.
Our goal here is twofold: i) This section is a preparation for the more complicated case of non abelian plasmas presented in the next section. Some of the results obtained here will indeed extend trivially to QCD, to within multiplicative color factors. ii) We wish to establish the relation with results obtained previously for the influence functional obtained in the path integral formalism. In particular, we shall show that we obtain, after a semi-JHEP06(2018)034 classical approximation, the same Fokker-Planck equations, and the associated Langevin equations, as derived earlier in the path integral formalism in ref. [24].

Single particle density matrix
In first quantization, the charge density n(x) of a single heavy quark is the operator We also need the matrix elements of the time derivative of the density. These can be easily obtained from the continuity equation, where the matrix elements of the current j(x) are given by One then gets We can then proceed to the evaluation of the various contributions L i in eq. (2.26), first in the case of the single particle density matrix. It is easy to show that the first line of eq. (2.26) yields Thus, the real part of the potential does not contribute to the evolution of the single particle density matrix. In terms of diagrams, this results from the cancellation of the tadpole insertions in the upper and lower branches (see the second diagram in figure 4), which represent here (unphysical) self-interactions.
Taking the matrix element of the second line of eq. (2.26), one obtains where we have set This equation illustrates the role of the collisions, captured here by the imaginary part of the potential, in the phenomenon of decoherence (the damping of the off-diagonal matrix elements of the density matrix). In contrast to what happens with the real part of the potential that we have just discussed, in the present case the two tadpole contributions add up, instead of cancelling. They are in fact needed to properly define the damping rate Γ, and insure in particular that it cancels when r → 0, so that the density (the diagonal part of the density matrix) is not affected by the collisions.

JHEP06(2018)034
A useful estimate of Γ(r) is obtained in the Hard Thermal Loop approximation [36][37][38] which gives where T is the temperature, and φ(x) a monotonously increasing function such that φ(x = 0) = 0 and φ(x → ∞) = 1 [15]. The same formula holds in the case of QCD, with α replaced by α s , the strong coupling constant, and the multiplication by appropriate color factors (see section 5). In the limit of a large separation, Γ(r) 2γ Q , where γ Q = α s T /2 is the so-called damping factor of a heavy quark (or antiquark) [39]. At small separation, interferences cancel the effect of collisions: the heavy quark pair is seen then as a small electric dipole, i.e., an electrically neutral object on the scale of the wavelengths of the typical modes of the plasma. Note that at large separation, the imaginary part of the potential itself vanishes, so that W (0) = −2γ Q . Considering finally the third line of eq. (2.26) one gets 4 The spatial derivatives originate from the time derivatives of the density (see eq. (3.3)), which involve the velocity of the heavy quark (hence the factor 1/M ). In fact, there is a close correspondence between L 3 and L 2 . Observe indeed that L 3 can be obtained from L 2 by multiplying the latter by the overall factor 1/(4M T ), and performing the following substitutions: We shall see that analogous correspondences also exist in the more complicated case of the 2 particle density matrix. At this point, we make the following change of variables and set r|D(t)|r = D(R, y, t). (3.10) The equation (2.24) becomes then d dt D(R, y, t) = LD(R, y, t), with L appearing now explicitly as a differential operator acting on the function D(R, y, t): The first term arises from the kinetic energy, i.e., it represents L 0 . Note that the other terms, which represent the effect of the collisions, vanish for y = 0, in particular thanks to the property ∇W (0) = 0. As already mentioned, this reflects the fact that the collisions do not change the local density of heavy quarks.

JHEP06(2018)034
Equation (3.11) above represents the explicit form of the operators L i in eq. (2.26) for the density matrix of a single heavy quark (in the abelian case). It has been obtained without any additional approximation beyond those leading to eq. (2.26). We may proceed further and simplify eq. (3.11) by performing a small y expansion. The variable y measures by how much the density matrix deviates from a diagonal matrix, a situation which is reached in the classical limit. Thus, the small y expansion may be viewed as a semiclassical expansion. We have where H(0) is the (positive definite) Hessian matrix of W , evaluated at y = 0, and we have used ∂ y W (y)| y=0 = 0. Note that if we stop the expansion of W (y) at quadratic order, ∇ 2 W (0) − ∇ 2 y W (y) = 0. The differential operator (3.11) reads then At this point some comments on the order of magnitude of the various terms are appropriate. It is convenient to measure the time in terms of the inverse temperature, setting for instance τ = T t. Dividing both sides of the equation by T , on gets on the left hand side ∂ τ , and the operator L/T on the right hand side is dimensionless. We shall assume in this paper that the heavy particles are initially close to rest. In interacting with the medium they ultimately thermalize, their velocity reaching values of order T /M , so that ∇ R √ M T . The variable y measures the non locality of the density matrix. When the heavy quark is not too far from equilibrium, this non locality is of the order of the thermal wavelength, that is D(R, y, t) dies out when y λ th ∼ 1/ √ M T . Thus in the first term, typically ∇ y ∼ √ M T , so that ∇ R · ∇ y ∼ M T . It follows that the term L 0 /T , where L 0 represents the kinetic energy of the heavy quark, is of order unity, while the other two terms are both of order Γ(y). The range of variation of Γ(y) is controlled by the Debye mass, i.e., it varies little on the scale of the thermal wavelength of the heavy particles. More precisely, using the HTL estimate Γ(r) ≈ αT (m D r) 2 , we get Γ(y)/T ≈ αm 2 D /(M T ) 1, the inequality following from our assumption M T , and the fact that m D T (in strict weak coupling m 2 D ≈ αT 2 ). In summary, the ratio of the last two terms in eq. (3.14) to the kinetic term is of order αm 2 D /(M T ) 1, which justifies the semi-classical expansion. To see better the physical content of eq. (3.14), we take its Wigner transform with respect to y. We define, with a slight abuse of notation, and obtain

JHEP06(2018)034
The corresponding equation for D(R, p, t) may be interpreted as a Fokker-Planck equation. The second term in eq. (3.16), proportional to ∇ 2 p can be viewed as a diffusion term (in momentum space), and is associated with a noise term in the corresponding Langevin equation (see below). It originates from the contribution L 2 . The last term, steming from the opeartor L 3 , can be associated with friction. This can be made more transparent by introducing the following definitions Then we operator above yields the following Fokker-Planck equation where v ≡ p/M is the velocity of the particle. It is easily shown that this equation can be obtained from the following Langevin equation The relation η = 2γT between the diffusion constant η and the friction coefficient γ can be viewed as an Einstein equation and expresses the fact that both noise and friction have the same origin, as can be made obvious by rewriting eq. (3.16) as follows (3.20)

The two particle density matrix
We consider now a heavy quark-antiquark pair. The charge density operator is written as where the first term refers to the quark and the second to the antiquark, the minus sign reflecting the fact that the antiquark has a charge opposite to that of the quark. The matrix elements of n(x) are given by Similarly, the matrix elements of the time derivative of the density are given by which is easily obtained from eq. (3.3). Note that we have introduced here a short notation, r ij ≡ r i −r j , that will be used often in the following. We shall also occasionally write ∇ 1 for ∇ r 1 , and introduce similar other shorthands, in order to reduce the size of some formulae.

JHEP06(2018)034
The various coordinates that are used in the evaluation of the two particle density matrix elements r 1 r 2 |D Q |r 1 r 2 .
It will be also convenient at a later stage to change variables. Thus, we define the center of mass and relative coordinates, as well as the set of coordinates that generalize those introduced in eq. (3.9) for the single particle case, The latter are useful to derive the semi-classical approximation. In this approximation, Y → 0, y → 0, and R and r become respectively the center of mass and the relative coordinates. These various coordinates are illustrated in figure 5.
We now turn to the specific evaluation of the matrix elements of eq. (2.24) in the case of a quark-antiquark pair. Consider first the matrix element of the free evolution, governed by the hamiltonian We have that is Turning now to the operator L 1 , a simple calculation yields JHEP06(2018)034 Figure 6. Graphical illustrations for typical contributions to the operators L i for the two-particle density matrix. In the first two diagrams, the gluon line represents either V or W , while in the last two, only W and its spatial derivatives are involved (the hamiltonian evolution, involving the real part of the potential, does not connect the upper and lower parts of the diagrams). In the last diagram, the gluon line connects two particles with the same charge, and contribute to the quantity called W a . In the third diagram, the gluon line connects two particles with opposite charges, and contributes to the quantity called W a . When W is involved in the first diagram, it represents a contribution to W c , and finally the tadpole insertion in the second diagram is associated to V (0), or to W (0) and its spatial drivatives. The defiitions of W a,b,c are given in eq. (3.32).
Note the cancellation of the self-interaction terms, as was the case for the single particle density matrix. The real part of the potential produces just a phase in the evolution of the density matrix. This can be understood as a hamiltonian evolution, with here H → −V , the minus sign resulting from the fact that the two interacting heavy particles have opposite charges. As we have mentioned earlier (see the discussion after eq. (2.24)), the structure of the equation makes it possible to include in the potential V both the "direct" interaction between the heavy quarks, by which we mean the interaction that exists in the absence of the plasma, as well as the "induced" interaction that results from the interaction of the heavy quarks with the plasma constituents. The latter is responsible for the screening phenomenon. In the HTL approximation, we have where the first term cancels the constant contribution hidden in the screened Coulomb potential (the second term) at short distance. Thus as r → 0, V (r) reduces to the Coulomb potential, V (r) ∼ α/r. Note that the function V (r) thus defined corresponds to the interaction potential of two particles with identical charges. By taking the matrix element of the second line of eq. (2.26), we obtain The various terms in this expression correspond to the various ways the exchanged gluon can be hooked on the upper and lower lines. To simplify the bookkeeping of the various contributions, and the writing of the equations, we define the following quantities, which JHEP06(2018)034 will often appear in forthcoming formulae: These quantities correspond to the diagrams in figure 6, where W plays the role of the propagator: W a connects particles with the same charge in the bra and the ket in eq. (3.31), while W b connects particles with opposite charges; W c connects particles within the bra, or within the ket. In the infinite mass limit, is defined in eq. (3.6). As was the case for the single particle density matrix, the collisions tend to equalize the coordinates (here the relative coordinates) in the ket and in the bra, bringing the density matrix to the diagonal form. The structure of the entire damping factor takes actually a form more complicated than in the case of the single particle density matrix. The combination of terms in the right hand side of eq. (3.31) can indeed be written Note that the entire contribution vanishes when r 1 = r 1 and r 2 = r 2 : Γ 11 → Γ(0) = 0, and similarly for Γ 22 while the other terms mutually cancel. This is of course related to the fact that the collisions do not change the local density of heavy particles, as we have already discussed. For future reference, we write L 2 as a sum of two contributions (as explained at the end of section 2), and write The diagonal elements (r 1 = r 1 , r 2 = r 2 ) of L 2a and L 2b mutually cancel, as we have seen. Finally, we turn to the 1/M corrections, which are more involved. The calculations are straightforward, but lengthy. However, as we shall see, the results are simply related to those obtained for L 2 . Again, we split L 3 into two contributions, L 3 = L 3a + L 3b . We obtain where we have used ∇W (0) = 0, and introduced the following shorthand notation with analogous definitions for W a , W b , W ± (to be used later). In this formula, and recall that ∇ 1 stands for ∇ r 1 and W 11 for W (r 1 − r 1 ).

JHEP06(2018)034
The second contribution to L 3 reads Note the analogy between eqs. (3.34) and eqs. (3.35) and (3.38). The latter follow from the former after the replacements and the multiplication by the overall factor 1/(4M T ). This correspondence is analogous to that already observed after eq. (3.8), and its origin is the same. Until now, the expressions that we have obtained are an exact transcription of the operators L i in eq. (2.26) to the case of a (abelian) quark-antiquark pair. At this point it is useful to go to the coordinates (3.25), i.e., r 1 r 2 |D|r 1 r 2 → D(R, Y , r, y), and perform a semi-classical expansion similar to that which leads to eq. (3.14) for the single particle density matrix. We obtain then After performing the Wigner transform with respect to the variables Y and y, we obtain We note that the operators for the relative coordinates are independent of the center of mass coordinates. It is then easy to identify the operators for the relative coordinates,

JHEP06(2018)034
and determine the elements of the corresponding Langevin equation. The relative velocity is given by 2p/M = p/(M/2). Thus Similarly, The term L 2 corresponds to the noise term. We write it as Finally L 3 corresponds to the friction, and we write it as Friction and noise are related by an Einstein relation The Langevin equation associated with the relative motion is then of the form Note that for an isotropic plasma, we have (cf. eq. (3.17)) One can repeat the same for the center of mass coordinates. Here we set v = P /(2M ). We get The Langevin equation associated with the center of mass motion is then These two equations (3.49) and (3.54) are identical to those obtained in [24] (see eqs. (4.69) there). Note that while the Langevin equation for the relative motion does not depend on the center of mass, this is not so for the Langevin equation describing the center of mass motion, which depends on the relative coordinate: this is because, as we have already emphasized, the effect of the collisions on the system depends on its size, with small size dipoles being little affected by the typical plasma field fluctuations.

QCD
We turn now to QCD. Much of the calculations are similar to those of the QED case, with however the obvious additional complications related to the color algebra. As we did in QED, we shall consider successively the case of the single particle density matrix, and that of a quark-antiquark pair.

Single quark density matrix
For a single quark, the color charge density can be written as (see eq. (2.3)) with matrix elements r, α|n a (x)|r , α = δ(r − r )n(r) α|t a |α , where n(r) is the density of heavy quarks, that is, the number of heavy quarks located at point r irrespective of their color state. The reduced density matrix of a single quark can be written as follows (see appendix D) It depends on 9 real parameters, and contains a scalar as well as a vector (octet) contributions. In fact, since we assume the plasma to be color neutral, we need consider only the scalar part of the density matrix (see however appendix G), that is the quantity r|D 0 |r . With D having only a scalar component, i.e., D = D 0 I, one can perform immediately the sum over the color indices in eq. (2.24), using The result is then identical to that obtained in QED, to within the multiplicative factor C F : there is no specific effect of the color degree of freedom on the color singlet component of the density matrix, aside from this multiplicative color factor. The resulting Fokker-Planck equation is then essentially identical to that first derived by Svetitsky long ago [40], which has been used in numerous phenomenological applications.

Density matrix of a quark-antiquark pair
The color density of a quark-antiquark pair is given by eq. (2.3). The color structure of the reduced density matrix D is discussed in appendix D. We shall use two convenient representations. In the first one, to which we refer as the (D 0 , D 8 ) basis, the density matrix takes the form where D 0 and D 8 are matrices in coordinate space (product of coordinate spaces of the quark and the antiquark), e.g. r 1 , r 2 |D 0 |r 1 , r 2 , and similarly for D 8 . The second representation is in terms of components D s and D o defined by (see appendix D) where |s s| denotes a projector on a color singlet state, and |o c o c | a projector on a color octet state with given projection c. We shall refer to this basis as the (D s , D o ) basis, or as the singlet-octet basis. The relation between the two basis is given by the following equations The advantage of the singlet-octet basis is that it involves states of the quark-antiquark pair with well defined color (singlet or octet), which is not the case in the (D 0 , D 8 ) basis. The latter will play a role when we address the issue of equilibration of color. Then the matrix D 0 , which represents a completely unpolarized color state, or a maximum (color) entropy state, plays an essential role. Because of the existence of two independent functions of the coordinates to describe the density matrix of the quark-antiquark pair, the equation of motion for D takes the form of a matrix equation where D can be viewed as a two dimensional vector, e.g.
and L is a 2 × 2 matrix in the corresponding space. The derivation of the equations for the various components of the reduced density matrix in the two basis is straightforward, but lengthy. The results are listed in appendix F. In this section we shall give brief indications on how to obtain the equations in the singletoctet basis: the color algebra is then transparent, and most of the equations can be related to the corresponding ones of the abelian case.

The equations in the singlet-octet basis
The contribution to L 1 . When written in terms of D s and D o the two equations decouple. This is because the product n a x n a x is a scalar in color space, and hence does not change the color state of the pair. In other terms, the one-gluon exchange does not change the color state (singlet or octet) of the heavy quark-antiquark pair. The color algebra is then trivial, and yields The diagrams contributing here are the first two diagrams in figure 6, and the equations above are analog to that obtained in QED. In fact we have where L QED 1 is given by eq. (3.29). Note the well known property that the interaction is attractive (and proportional to C F ) in the singlet channel, and repulsive (and proportional to 1/(2N c )) in the octet channel.
The contribution to L 2 . In this case we have two types of contributions. The first ones involve products of the color charges, making up a color scalar. These contribute to L 2a , which is diagonal in color. The second type of contribution involves transitions from singlet to octet or from octet to octet. We denote these contributions by L 2b .
Consider first L 2a . We have, for the singlet with L QED 2a given in eqs. (3.34). For the octet Consider next L 2b . The corresponding contributions involve transitions from singlet to octet intermediate states, or, when considering the octet density matrix, transitions from octet to singlet and also octet to octet transitions that generate an additional diagonal contribution. We have where L QED 2b is given in eqs. (3.34). Similarly, for the octet to octet transition, we have which has no counterpart in QED.

JHEP06(2018)034
The contribution to L 3 . The calculation of these contributions is more involved, but we can relate simply the results to those obtained earlier for the operators L 2 . Let us first list the results. We have and It is easy to verify that the equations giving L 2 and the corresponding equations for L 3 are related via the same substitutions that are discussed after eq. (3.38).

Semiclassical approximation
The formulae listed in the previous subsection are an exact transcription of the main equation, eq. (2.24) for a quark-antiquark pair. Analogous formulae can be written for a pair of quarks. They are given in appendix H. We shall be mostly concerned in this paper, in particular for the numerical studies presented in the next section, by the semi-classical limits of these equations. These can be obtained easily by using the formulae given in appendix F.3. In this subsection, we just reproduce the few equations that will be used in the next section. We consider first the equation for the component D s , which we write as dD s dt = (D s |L|D). (4.21) We obtain

JHEP06(2018)034
One reason why we display this equation is that it reduces to the QED equation when D s = D o . Thus, if color quickly equilibrates, an assumption that we shall exploit in the next section, the dynamics becomes analogous to that of the abelian case. In this case, color degrees of freedom play a minor role, and the motion of the heavy particles can be described by a Fokker-Planck equation or the associated Langevin equation.
As we have already emphasized, the component D 0 corresponds to the maximum (color) entropy state, where all colors are equally probable. This state plays an important role in the calculations to be presented in the next section. Thus, it is useful to write the corresponding equations of motion, or equivalently the operators L of the (D 0 , D 8 ) basis, in the semi-classical limit. We have: and We shall also need the operators L 08 and L 88 at leading order in y. These are given by

Numerical studies
The equations for the time evolution of the reduced density matrix that we have obtained in the previous sections are difficult to solve in their original form, that is, for the operator L given in section 4.3, or appendix F, for a quark-antiquark pair. We shall not attempt to solve them directly in the present paper. In the case of QED, we have seen that an additional approximation, the semi-classical approximation, allows us to transform these equations into Fokker-Planck, or equivalently, Langevin equations, describing the relative and center of mass motions of the heavy particles as simple random walks. In QCD, the presence of transitions between singlet and octet color states complicates the situation, since such transitions are a priori not amenable to a semi-classical description. The purpose of this section is to present numerical studies that illustrate two possible strategies to cope with this problem, namely preserve as much as possible the simplicity of the semi-classical description of the heavy particle motions, while taking into account the effects of color transitions. To simplify the discussion we shall ignore the center of mass motion in most of this section. The new feature in QCD, as compared to QED, namely the transitions between distinct color states, is best seen in the infinite mass limit, where the relative motion is entirely JHEP06(2018)034 frozen. Then the only dynamics is that of color: as a result of the collisions with the plasma constituents the colors of the heavy quarks and antiquarks can change in time. The corresponding equations of motion for the density matrix are easily obtained from the formulae listed in appendix F. They read, for a quark-antiquark pair, where r is the (fixed) relative coordinate. These equations exhibit the decay widths in the singlet (2C F Γ(r)) and the octet ((1/N c )Γ(r)) channels, respectively. These two coupled equations acquire a more transparent physical interpretation in the (D 0 , D 8 ) basis, where they take a diagonal form The first equation merely reflects the conservation of the trace of the density matrix.
Recall also that D 0 is associated with the maximum color entropy state, where all colors are equally probable (see eq. (D.13)): this component of the density matrix represents an equilibrium state that remains unaffected by the collisions. The second equation indicates that D 8 ∝ D s − D o decays on a time scale (N c Γ(r)) −1 . Thus, at large times only D 0 survives, that is, the collisions drive the system to the maximum entropy state. Note that the distance between the quark and the antiquark enters the rate ∝ Γ(r) at which equilibrium is reached. When |r| m D , Γ(r) ≈ 2γ Q , where γ Q is the damping factor of one heavy quark (or antiquark): at large separation, the quark and the antiquark equilibrate their color independently (with a rate N c γ Q -see appendix G). On the other hand, when |r| m D the two quarks screen their respective colors, hindering the effect of collisions in equilibrating color. The first strategy that we shall explore in order to treat the relative motion of the heavy quarks semi-classically and at the same time take into account the color transitions just discussed, rests on the following observation. There is one instance where one can recover a situation analogous to that of QED: this is the regime where the color exchanges are fast enough to equilibrate the color on a short time scale (short compared to the typical time scale of the relative motion). We have for instance observed in the previous section that the component D s of the density matrix, eq. (4.22), obeys the same equation as the QED density matrix when D o = D s (to within the multiplicative color factor C F ), which corresponds indeed to the maximum entropy state. In this case, one can explore the dynamics in the vicinity of this particular color state, treating the color transitions in perturbation theory. One can then derive Langevin equations which contain an additional random force arising from the fluctuations of the color force between the heavy particles. This perturbative approach is easily generalized to the case of a large number of quarkantiquark pairs.

JHEP06(2018)034
The second strategy is based on an analogy between the equations (5.1), and their generalizations to include the semi-classical corrections, and a Boltzmann equation, with the right hand side being viewed as a collision term. That is, the changes of color that accompany the singlet-octet transitions are then treated as collisions rather than an additional random force in a Langevin equation. This strategy allows us to overcome some of the limitations of the perturbative approach.

Langevin equation with a random color force: single quark-antiquark pair
We shall now examine the corrections to eq. (5.2) that arise in the semi-classical approximation, i.e. taking into account corrections to the infinite mass limit. Note first that the kinetic energy of the heavy quarks leaves L as a diagonal operator. In the (D 0 , D 8 ) basis, this operator reads The semi-classical corrections brings L to the form   One can diagonalize this new operator L, and in particular find the eigenvalue that corresponds to the maximum entropy state in the limit where y → 0. This is given by usual perturbation theory By using the explicit expressions given above for the various coefficients, one deduces that the corresponding eigenvector, D 0 , fulfills the equation

JHEP06(2018)034
Performing the Wigner transform we obtain The comparison with the Fokker-Planck operators given in eqs. (3.43) to (3.47), allows us to write the corresponding Langevin equation: where and .

(5.12)
As compared to the QED case, eq. (3.49), we note a number of new features. The most noteworthy is the presence of two random forces. The force ξ is the familiar stochastic force and is related to the drag force as indicated in eq. (5.11). The second random force, Θ, has a different nature: it originates from the fact that the force between a quark and an antiquark is a function of their color state. Now, D 0 represents a state close to the maximum entropy state, that is, a state in which the probability to find the pair in an octet is approximately N 2 c − 1 times bigger than the probability to find it in a singlet. At the same time, the force between heavy quarks in a singlet state is N 2 c − 1 times bigger than that in an octet state, and it has opposite sign. The net result is that, on average, the force between the heavy quark and the heavy antiquark is zero. But this is true only in average. There are fluctuations, and these give rise to Θ. The vanishing of the average force between the quark and the antiquark explains the absence of the force term in the Langevin equation, as compared to the QED case. Note that this picture is valid as long as transitions between singlet and octet states are fast compared with the rest of the dynamics. This is no longer the case when the size of the pair becomes too small: then, Γ(r) becomes small, reducing the energy denominator in eq. (5.12), i.e., amplifying the effect of the random color force. This, as we shall see, can lead to unphysical behavior.

Many heavy quarks and antiquarks
The discussion of the previous subsection can be generalized to a system containing many heavy quarks and antiquarks. We call N Q the number of heavy quarks, and for simplicity we assume that it is equal to the number of heavy antiquarks. The density matrix can be written as a product of density matrices of the individual quarks and antiquarks, generalizing the construction of eq. (D.4) for the quark-antiquark density matrix. One can then write where the dots represent all the scalar combinations that can be formed with products of n ≤ 2N Q color matrices t a , with coefficients corresponding to the components of D. We

JHEP06(2018)034
shall not need the explicit form of these extra components. As for D 0 , this is clearly the maximum entropy state, where all colors of individual quarks and antiquarks are equally probable and uncorrelated. It corresponds to the following density matrix (cp. eq. (D.13)): where the sum runs over all the colors of the quarks (α 1 , · · · α N Q ) and the antiquarks (ᾱ 1 , · · · ,ᾱ N Q ). We want to construct for this system the analog of eq. (5.8) which describes how the maximum entropy state is modified by the semiclassical corrections. Our starting point is the main equation, eq. (2.24). To proceed, it is useful to have in mind the diagrammatic representation of the density matrix that we have introduced earlier. As compared to the diagrams displayed in figure 6, in the present case, the diagrams will contain 2N Q lines in the upper part, and 2N Q lines in the lower part. All the interactions that we are dealing with involve a single gluon exchange, represented by one gluon line joining quark or antiquark lines in various ways. The evolution equation for D is still described by an operator L, which is a matrix in the space of all the independent components. For our perturbative calculation, we need only to consider diagonal (to order y 2 ) and non diagonal (to order y) elements of this matrix, the non diagonal elements involving the maximum entropy state as one of their entries.
Consider first the diagonal elements. We have first the kinetic energy, trivially given by The leading order diagonal element that involves the interaction can be obtained in the infinite mass limit. It represents the decay of the components of D that are connected to D 0 by one gluon exchange. It is given by where r kl = r k − r l , r k and r l denoting the coordinates of the quark or the antiquarks to which the gluon is attached. The factor N c can be understood from the same argument as that given after eq. (5.2): when the separation r kl is large, the color of the quark (or antiquark) at r k and r l relax independently at a rate N c γ Q . At the order of interest, we need also the diagonal element for the maximum entropy state, including the semi-classical corrections up to order y 2 and y M . This is given by We turn now to the non-diagonal elements. To leading order, these involve solely the real part of the potential. Diagrammatically, the corresponding exchanged gluon connects only either the upper or the lower lines among themselves, but does not connect upper with lower lines.

JHEP06(2018)034
Consider first the non-diagonal elements that bring the maximum entropy state to another state. If the pair connected by the exchanged gluon is formed by quark k and antiquark l then the element is (cp. eq. (5.5)) while if it is formed by quark (antiquark) l and quark (antiquark) k then it is (cp. eq. (H.9)) We also need the non-diagonal elements that bring the system back to the maximum entropy state. If this pair is formed by quark k and antiquark l then the element is (cp. eq. (5.5)) − i y kl · F (r kl ) , (5.20) while if it is formed by quark (antiquark) l and quark (antiquark) k then it is (cp. eq. (H.10)) iy kl · F (r kl ) .

(5.21)
With this information, it is straightforward to construct the generalization of eq. (5.8) for an arbitrary number of quark-antiquark pairs. The corresponding equations will be presented later.

Simulation of a Langevin equation with a random color
In this subsection we present numerical results obtained by simulating the equations of the previous subsections. We examine successively the evolution of the relative coordinate of a single heavy quark-antiquark pair, and then that of fifty pairs initially produced in a thin layer. The first case will help us to understand the range of applicability of the perturbative method, while the second will illustrate how the Langevin equation may account for recombination.
In these calculations, we use the standard QCD running coupling constant α s determined at one loop order for N f = 3 massless flavors and Λ QCD = 250 MeV. The screening Debye mass is given by its HTL approximation, m 2 D = (2π/3)(6 + N f )α s T 2 , with α s evaluated at the scale 2πT , with T the temperature. Further details on the parameters of the calculation will be given as we proceed.
We should emphasize here that the numerical simulations to be presented in this section are meant to illustrate the main physical content, as well as the limitations, of the equations obtained in this paper. Although the numbers that go into the calculations are appropriate for the physics of quarkonia in a quark-gluon plasma, we make no attempt to develop a phenomenological discussion.

A single heavy quark-antiquark pair
The Langevin equation for the relative motion is given in eq. (5.10). The information about the medium is encoded in the functions V (r) and W (r) which we estimate using the HTL approximation. Note that the resulting value of ∆W (r) and V (r) diverge as r → 0 (for different reasons though, see e.g. [24]). In [24] the divergence of ∆W (r) was regulated with the help of an ultraviolet cut-off. Here, we shall use a simpler procedure and choose as a free parameter of our simulation, for which we choose the value γ = 0.19T 2 . For the real part of the potential, we proceed as in [24], that is we define it as the Fourier transform of a Yukawa potential integrated in momentum space up to a cutoff Λ = 4 GeV. The coupling constant that appears in V (r) is evaluated at a scale corresponding to the inverse Bohr radius of the bottomonium (1348 MeV). The spatial dependence of W (r) is obtained, as already mentioned, from the HTL approximation, and is of the form W (r) = W (0) + α s T φ(m D r) [15], with α s evaluated at the scale 2πT , and φ(x) a monotonously increasing function such that φ(x = 0) = 0 and φ(x → ∞) = 1. At small separation, i.e., for m D r 1, φ(x) can be approximated by φ(x) = x 2 3 (− log x + 4 3 − γ E ). In figure 7 we show a set of ten random trajectories produced for bottomonium (with mass M b = 4881 MeV) at a temperature of T = 350 MeV. We take as initial distribution of relative distances that obtained from the wave function of the 1S state (which we plot as an histogram with 1000 events in figure 8 for further comparison). We see that many of the trajectories are clearly unphysical since they involve supraluminal velocities: this is because some of the random kicks can occasionally be very strong due to the amplification produced by the small denominator in eq. (5.12) at small r. A more systematic comparison can be done by looking at the distribution of relative distances after some time t. This is shown in figure 9 for t = 5 fm. The histogram reveals that there remains at this time only  a tiny probability to find the pair within a relative distance corresponding to the size of the bound state: the random color force is clearly too efficient in suppressing the bound state!

Many heavy quark-antiquark pairs
In spite of its shortcomings, the perturbative method remains interesting as it allows us to treat an assembly with an arbitrary number of quark-antiquark pairs, and address in particular the issue of recombination. The relevant equations can be constructed along the lines developed in the previous subsection. We need to take into account not only the relative coordinates but also the center of mass motions. Moreover the random force does not only act between heavy quarks and antiquarks but also between two heavy quarks and two heavy antiquarks (the equations for the density matrix of a pair of two heavy quark are given in appendix H). The resulting equations are given by where the noises have vanishing means and correlators given by In these equations N Q is the number of heavy quarks, equal to the number of heavy antiquarks. The indices a or b are color indices in the 3 representation while the same with hat are a color indices in the3 representation. The nature of the color index specifies whether a given quantity refers to a quark or an antiquark. Thus, r a represents the coordinate of a quark, while râ represents the coordinate of an antiquark. We use also compact notation, such as r ab = r a − r b to denote the relative distance between a quark of color a and a quark of color b, or r ab = r a − rb to denote the relative distance between a quark of color a and an antiquark of colorb. Finally, for functions of coordinates, we set Fâ b = F (r a − r b ), Γ ab = Γ(r a − r b ), and so on. In figure 10 we plot some random trajectories of fifty pairs of quarks and antiquarks. The parameters are different from the ones used in the previous section, now the temperature is T = 250 MeV and the cut-off for V (r) is Λ = 1500 MeV. We keep the same value of JHEP06(2018)034  Figure 11. Histogram of the number of bound states formed in 300 simulations with the same initial conditions as in figure 10 γ. This new choice of parameters makes the problem of the violent hard kicks less severe, at the cost of having a cutoff Λ unrealistically small (it is of the order of the inverse of the Bohr radius of the ground state). The system is prepared in the following way: in a square of size 2.5 fm we chose fifty random points; in each point we put a quark-antiquark pair following a probability for the relative coordinate given by 2 πr sin(Λr) , (5.29) where Λ is the same cut-off as for the real part of the potential (this distribution becomes a Dirac delta as Λ → ∞). The fifty quark-antiquark pairs then evolve for a time t = 5 fm/c, according to the stochastic equations displayed above. As can be observed by looking at figure 10 some of the trajectories remain close enough to allow for "recombinations" into bound state.
To quantify the phenomenon, we perform a statistical analysis of how many bound states are observed at the end of the evolution, starting from the previous initial condition. To define a bound state, we follow the procedure of ref. [24], but with slightly different parameters. We declare a heavy quark-antiquark pair to be bound if the quark and the antiquark remain at a distance smaller than 0.5 fm during a time bigger than 0.1 fm/c. This procedure can become ambiguous when many quarks and antiquarks are found in a small region, for example in the case in which two quarks and two antiquarks are contained in a sphere of radius smaller than 0.25 fm. In this situation we count the number of bound states by choosing the combination that yields the maximum number. The results obtained after 300 simulations are shown in figure 11. We can see that from the 50 initial bound states, on average 17 remain after a time of 5 fm/c spent inside the medium. This is to be contrasted for figure 9 which would suggest that all pairs should become unbound if they were evolving independently of each other. Of course, we should recall that different parameters have been used in the two calculations. However, repeating the simulation of JHEP06(2018)034 the previous subsection with the present parameters, one finds that about 10 to 15% of the bound states would survive after a time t = 5 fm/c. This is about half of what the present calculation suggests. We may therefore take this result as evidence for recombination, in line with what was found in ref. [24].

Langevin equations coupled to random collisions
We now turn to the second strategy presented in the introduction of this section. This will allow us, in particular, to bypass the limitations of the perturbative diagonalization of the matrix in eq. (5.4), caused by the vanishing of Γ(r) at small r. Let us then go back to the small y expansion of eq. (4.22), and let us temporarily neglect the terms that go like y M or y 2 , that is we keep only the kinetic term and the force term. We get The equation for D o is not given explicitly in eq. (4.22), but it is easily derived from the material presented in appendix F. We note that only the terms proportional to Γ(r) mix singlets and octets, i.e. the terms involving the force preserve the color state of the pair. We now perform a Wigner transform with respect to the variable y, and define P s = D s and P o = (N 2 c −1)D o , the probabilities for the pair to be in a singlet or octet state, respectively. We get The right hand sides of these two equations can be interpreted as a "collision term" in a Boltzmann equation, with gain and loss terms. Note that these collision terms are opposite in the singlet and octet channels, as expected: The left hand sides of the equations (5.32) and (5.33) describe the relative motion of the pair under the influence of the color force F (r). The corresponding classical equations of motion read where the last two equations refer to the singlet and octet channels, respectively. Thus, instead of treating the singlet-octet transitions as an additional color force in a Langevin equation, we can treat these transitions as "collisions". In practice we can solve the set of equations (5.32) and (5.33) using a Monte Carlo method, deciding at each time step, according to a probability proportional to the respective decay widths, whether a JHEP06(2018)034 transition takes place or not, and then evolve the system through the time step according to the classical equations of motion (5.35). This is somewhat analogous to the Monte Carlo Wave Function method applied to a 2-level problem in ref. [41].
The equations (5.32) and (5.33) capture some of the important physics but they miss the drag forces and the stochastic forces that have to go with them in order to fulfill the fluctuation-dissipation theorem. These come from the semi-classical corrections that we have left out in writing eqs. (5.30). However, if we were to include these corrections as they appear for instance in eq. (4.22), we would introduce extra couplings between D s and D o that would lead in particular to a collision term involving derivatives of P s,o and we do not know of any efficient numerical tools to solve the resulting equation. However, if the system is not too far from the maximum entropy state, terms that go like y 2 P s − Po c −1 are small and can be safely neglected. We again rely on the assumption that color relaxes faster than the relative motion. Under these conditions, and after performing the Wigner transform, we obtain for the singlet Comparing the operator in the first line of this equation with that given in eqs. (3.43) to (3.47), we easily derive the following stochastic equations: These equations are the analogs of eq. (3.49) for the QED case. Performing completely analogous manipulations, we obtain a similar result for the octet. The Fokker-Planck equation reads and the corresponding Langevin equation is with now We shall now present the results of the simulation of these equations, for the case of a single quark-antiquark pair. We consider first the static limit, and then turn to the full equations including the semi-classical corrections.

The static limit
The study of the static limit (or infinite mass limit) offers us the possibility to test the numerical method, since the exact solution can be obtained analytically in this case. In particular, this will give us an idea of the number of iterations that are needed in order to get a good estimate. We consider a heavy quark-antiquark pair, whose relative distance is r = 0.1 fm, and in a well-defined color state, in a quark-gluon plasma at temperature T = 250 MeV.
The equations to be solved are eqs. (5.1). If the initial conditions are such that singlet or and octet states are equally probable, i.e., P s (0) = P o (0), the probability to be in an octet state at time t is and that to be in a singlet state is P s (t) = 1 − P o (t).
We can compare this result to that of a simulation using the Monte Carlo method described above. The results are plotted in figure 12, as a function of τ ≡ (2C F α s T t)/10, with a time step ∆τ = 0.02. We see that for 100 events the results of the simulation match relatively well the analytic result, although sizeable fluctuations remain. The simulations to be presented next involve 1000 events.  Figure 14. Probability to find a heavy quark-antiquark pair in an octet state at time t, after it has been prepared at time t = 0 as a J/Ψ state.

Simulation with dynamical quarks
We consider now the full eqs. (5.36) and (5.40). In figure 13 we plot the average mean distance r cc of a pair of charm quarks prepared in a J/Ψ, according to the prescriptions used in ref. [24]. That is, the radius of the pair is chosen randomly between 0 and 1/m D , and the relative momentum is chosen according to a Maxwell distribution with most probable velocity given by v 2 = 0.3. Finally one retains only pairs with binding energy bigger than 500 MeV and radius bigger than 0.1 fm. The temperature is taken to be T = 160 MeV, the charm quark mass M c = 1460 MeV, and γ/M c = 0.2 fm. These conditions differ slightly from those used earlier: the reason is that we want to check our results against those of [24] JHEP06(2018)034 in a domain where they could be compared (see figure 13). Thus, we also use a different running coupling than above, α s = 0.5/(1 + 0.76 ln(T /160)), and m D = (16πα s /3)T 2 (for two massless flavors). The cutoff on V (r) is 4GeV, while that on W (r) is 4.58 m D . The unit of time in figure 13 is the physical unit fm/c.
We compare the results of the simulation of eqs. (5.36) and (5.40) with those obtained by neglecting color rotation, i.e., the singlet-octet transitions. The latter case is equivalent to a QED simulation, and indeed our result in that case reproduce those obtained in [24] (cp. the corresponding result in figure 13 with figure 5 in [24]). As expected, we see that the bound state tends to remain bound longer if the transition to octet is not taken into account. The effect of color rotation is clearly to accelerate the melting of the bound state, although, according to the criterion used in [24], r cc r D = m −1 D , we may consider the system to remain bound at time t = 4 fm/c. This is to be contrasted with the result obtained with the Langevin equation with a color random force: in the present case, the disappearance of the bound state is a more gradual phenomenon, not amplified by unphysical violent kicks of a random color force. This gradual transition can be visualized by looking at the evolution of the probability to find the pair in an octet state, which is plotted in figure 14.
We can see that it takes a non negligible time to lose the information that the system was initially in a singlet state.

Summary and outlook
In this paper we have obtained a set of equations for the time evolution of the reduced density matrix of a collection of quark-antiquark pairs immersed in a quark-gluon plasma in thermal equilibrium. These equations are fairly general (they are valid for an arbitrary number of heavy particles), and rely on two major approximations: weak coupling between the heavy quarks and the quark-gluon plasma, small frequency approximation for the plasma response. In the weak coupling approximation, the plasma sees the heavy quarks as a perturbation, and responds linearly to it. This response is characterized by a set of correlators, expectations values of gauge fields in the equilibrium state of the plasma, and because the heavy quark motion is slow on the typical scale of the plasma dynamics, only static, or nearly static response functions are needed. These functions account for some of the dominant effects of the plasma on the dynamics of the heavy quarks: the screening of the instantaneous Coulomb interaction between the heavy quarks, and the effect of soft, low momentum transfer, collisions of the heavy quarks with the plasma constituents taken into account by an imaginary potential. The main equations that result from these two approximations alone generalize the equations that were obtained for an abelian plasma in the path integral formalism, using the Feynman-Vernon influence functional method [24]. Their structure is close to that of a Lindblad equation, and they are essentially equivalent to the equations obtained for QCD by Akamatsu [23], although the present formulation differs from his in several aspects. Recently a Lindblad equation was obtained for the evolution of the density matrix of a quark-antiquark pair, using similar approximations, but formulated in the context of a non relativistic effective theory (pNRQCD, [28,42]). This formalism puts the emphasis on the singlet-octet transitions, and the validity of the JHEP06(2018)034 employed effective theory requires specific conditions, viz. 1/r T ∼ E, with E the typical binding energy. The corresponding Lindblad equation keeps the quantum features of the problem, however at the price of a high computational cost.
In the case of abelian plasmas, a further approximation, the semi-classical approximation, leads to a Fokker-Planck equation, and a corresponding Langevin equation, which are relatively easy to solve numerically. When trying to extend this semi-classical approximation to QCD, we have to face new features related to color dynamics. In the particular case of a quark-antiquark pair, this involves the transitions between the singlet and the octet color configurations of the pair. Taking these transitions into account yields coupled equations for the two independent components of the density matrix, that are not easily solved, even when the motion of the heavy particles is treated semi-classically.
We have then explored numerically two strategies to solve approximately these coupled equations. In the first one, we assume that the color dynamics is fast compared to the motion of the heavy quarks. In this case, the collisions drive the systen quickly to a maximum entropy state where all colors are equally probable and uncorrelated. One can then use the Langevin equations to describe the dynamics in the vicinity of this maximum color entropy state, using a perturbative approach. This is sufficiently simple that it can be generalized to a system of an arbitrary number of quarks and antiquarks. However, the perturbative approach is limited by the fact that the color relaxation is slow when the size of the quark-antiquark pairs is small, which may lead to unphysical behavior for a physically relevant choice of parameters. To overcome this limitation, we have explored another strategy, which appear more promising. It consists in treating the singlet-octet transitions as collisions, viewing the corresponding equations as Boltzmann equations that we solved using Monte Carlo techniques.
Although they are fairly general, the equations that we have obtained so far do not yet capture all the relevant physics. For instance, in the particular case of a single quarkantiquark pair, the transitions between singlet and octet color states cause rapid changes in the heavy quark hamiltonian. These are not properly handled, and may be in conflict with the assumption that the dynamics of the heavy particles is slow compared to that of the plasma. In addition we have left aside the possibility of absorption or emission of real gluons from the plasma, that are responsible in particular for gluo-dissociation, known to be an important mechanism in some temperature range. These shortcomings, and further aspects of the problem, will be addressed in a forthcoming publication.

JHEP06(2018)034
and on the difference of coordinates, which we shall denote respectively by τ and x in this appendix. They are invariant under the change x → −x.

B Alternative time discretization
Our starting point is eq. (2.13), with symmetrical time integrations, of which we take the time derivative. We get

JHEP06(2018)034
At this point we use the values of the time integrals of the correlators that are given in appendix A, and obtain As mentioned in the main text (see the discussion after eq. (2.25)) the structure of this equation is close to that of a Lindblad equation. 6 To make this more obvious, let us Fourier transform the variables x and x . One obtains easily with the shorthand notation q = d 3 q/(2π) 3 , and q the variable conjugate to x in the Fourier transform. The second line of this equation has the structure of a Lindblad operator, but this is not so for the last two lines corresponding to the operator L 3 . However, it is easy to see that the substitution n q → n q +(i/4T )ṅ q in the second line, which obviously preserves its Lindblad structure, generates all the terms in the last two lines with, in addition, terms that are quadratic in the time derivative, and are therefore suppressed with respect to the other terms by a power 1/M T . Thus, to within these terms, one may consider eq. (B.9) as a Lindblad equation. 7 Note that the substitution mentioned above works only with the present time discretization, and is not immediately applicable to that used in the main text.

C Comparison between the two discretizations
The two time discretizations differ solely in their respective contributions to L 3 , and only in that part of it that we called L 3a . In this appendix, we examine this difference in the 6 Recall that one of the virtues of the Lindblad equation is to maintain the positivity of the density matrix [35]. 7 The necessity of additonnal terms quadratic in velocities in order to obtain the Lindblad equation is also discussed in [23].

JHEP06(2018)034
case of QED. The generalization to QCD is straightforward. In the main text, we obtained (eq. (3.35)) while the discretization presented in appendix B yields (eq. (B.8)) By taking the difference one obtains It is easy to verify that this difference does not contribute to the matrix elements of the single particle density matrix. Let us then consider the two particle density matrix. A straightforward calculation yields Changing to the variables of eq. (3.25), and performing the small y expansion, one gets Using the expression of L 3 from eq. (3.40), and L 3 − L 3 = L 3a − L 3a , we obtain After taking a Wigner transform this yields This yields a modified Fokker-Planck equation for the relative coordinate, as compared to that used in the main text, eq. (3.47). However, a simple change of variables allows us to recover the Langevin equation (3.49). To this end let us set Then, to within terms that are suppressed in the semi-classical approximation, we have and a simple calculation shows that, in terms of these variables, the Langevin equation corresponding to the operator L , with L 3 given in eq. (C.7) and L 1,2 = L 1,2 , is identical to eq. (3.49).

JHEP06(2018)034
D Color structure of the density matrix The density matrix of a color quark is a 3 × 3 matrix in color space, which can be written as follows where t i = λ i /2, with λ i the Gell-Mann matrices. We use the standard normalization Tr The density matrix (D.1) depends on 9 real parameters, and contains a scalar as well as a vector (octet) contributions. The density matrix associated to an antiquark may be written as A representation of the density matrix of a quark-antiquark pair may be obtained as the tensor product Using the identity one can write D as where |s denotes a color singlet and |o c a color octet, with projection c. The states |s and |o are normalized to unity s|s = 1, o c |o d = δ cd . We have The relations between the coefficients in the two bases are easily obtained. They are given by Note that the component D 0 corresponds to a completely unpolarized system, and can be written as The same density matrix in the singlet-octet basis corresponds to D s = D o .

E Some useful formulae and matrix elements
In this appendix, we list a number of useful formulae, as well as some matrix elements that facilitate the derivation of the equations presented in the main text. We start with relations involving color matrices in the fundamental representation. Using the relations and it is easy to establish the following formulae

JHEP06(2018)034
We also need We consider now matrix elements in the singlet-octet basis. We have The following matrix elements of the color charge density, or its time derivative, are also useful. We have singlet-octet matrix elements, where n(x) is the QED charge density We can write the formula above more simply as We have also octet-octet matrix elements, 10) or, more simply,

JHEP06(2018)034
F The equations of motion for the density matrix of a heavy quarkantiquark pair In this appendix we present the equations of motion for the matrix elements of the reduced density matrix of a heavy quark-antiquark pair. We consider the two representations of the density matrix that correspond to the (D 0 , D 8 ) and the (D s , D o ) basis. In the main text we have indicated how to obtain the equations for D s and D o from the corresponding equations of the abelian case. Depending on the choice of basis, the color algebra proceeds differently, and the relevant formulae are listed in appendix E. Independently of the way we proceed, the results should eventually lead to formulae for the various components of D Q which satisfy the relations (D.12). This constitutes a useful check of the results. The equations to be presented below correspond to the coordinate space matrix element r 1 r 2 |D|r 1 r 2 . Thus, for instance, D 0 stands for D 0 (r 1 , r 2 ; r 1 , r 2 ), and similarly for D 8 , or D s and D o . We use the notation (D s |L|D) for the contribution of the operator L i to the time derivative of D s , that is dD s dt = (D s |L i |D).
(F. 1) and similarly for the other components of the density matrix. Also we use the compact notation introduced in the main text, e.g., W 12 = W (r 1 − r 2 ), ∇ 1 = ∇ r 1 , ∇ 12 = ∇ r 1 − ∇ r 2 , and so on, as well as the quantities W a,b,c,± defined in the main text (see eqs. (3.32)). The various calculations are straightforward, but lengthy and somewhat tedious. They will not be presented here, we just list the final results. Those results, listed in the next two sections for the (D 0 , D 8 ) basis and the (D s , D o ) basis, respectively, are an exact transcription of eq. (2.26) to the density matrix of a quark-antiquark pair, without additional approximation. In the last subsection of this appendix, we list a number of formulae that are useful to implement the semi-classical approximation.
F.1 (D 0 , D 8 ) basis In the (D 0 , D 8 ) basis the contributions of the operators L i to time derivative of the density matrix are given by: Note that in the infinite mass limit, the contribution (D 0 |L 2 |D) vanishes, while the second contribution reduces to (D 8 |L 2 |D) = −N c Γ(r)D 8 . These results have a simple interpretation discussed in the main text.

F.2 (D s , D o ) basis
In the (D s , D o ) basis the contributions of the operators L i to time derivative of the density matrix are given by: Note the analogy with the equation (F.7): replace W (0) → ∇ 2 W (0), W c → ∇ 2 W c + ∇W c · W c and W − → ∇ 2 W − + ∇W − · ∇ − (the color algebra is the same). Also this equation is identical to that in QED when D o = D s .

F.3 Equations in the semi-classical approximation
To derive the equations in the semi-classical approximation, the following formulae are useful.  ∇W a · ∇ a ≡ ∇W 11 ∇ 11 + ∇W 22 ∇ 22 = 2Y · H(0) · ∇ Y + 2y · H(0) · ∇ y ∇W b · ∇ b ≡ ∇W 12 ∇ 12 + ∇W 21 ∇ 21 = 2∇W (r) · ∇ r + 2Y · H(r) · ∇ Y , ∇W c · ∇ c ≡ ∇W 12 · ∇ 12 + ∇W 1 2 · ∇ 1 2 = 2∇W (r) · ∇ r + 2y · H(r) · ∇ y , (F. 13) G A single heavy quark in the static limit As an illustration of the color dynamics we study the case of a single heavy quark in the infinite mass limit. The density matrix can be written as and we leave open the possibility of a vector component (D 8 ). In the infinite mass limit, D 0 and D 8 obey the equations The first equation reflects the conservation of the trace of the density matrix. The second equation indicates that the states that have a preferred direction in color space decay exponentially in time. To illustrate this behaviour imagine that at t = 0 we have a density matrix corresponding to a heavy quark with a specific color (such that D 11 = 1 and the rest of the components are 0). The evolution of this matrix can be written as

H The case with two heavy quarks
We consider here a system formed by two heavy quarks (rather than a heavy quarkantiquark pair). Eq. (2.24) is also fulfilled in this case, but the color structure of the where P3 and P 6 denote respectively the projectors on the representation3 and 6 of SU(3). In writing eq. (H.1) we stick to the notation of eq. (D.5), although in the present case there is no octet state involved in D 8 . As for D 0 it conserves the interpretation of the maximal entropy state. Calculating, as we did for the quark-antiquark pair, the matrix element r 1 r 2 |D|r 1 r 2 of the two-quark density matrix, we obtain for the first line of eq. (2.24) This set of equations has two eigenvalues with opposite signs (i(V 12 − V 1 2 )(1 ± N c )/(2N c )). The positive sign corresponds to the color configuration3 and represents an attractive interaction, while the minus sign corresponds to the configuration 6 for which the interaction is repulsive. The equations for D3 and D 6 read Note that the attraction between two quarks in the3 channel is N c − 1 times weaker than the attraction of a quark-antiquark pair in the singlet channel. For N c = 3 this is only a factor 2, which suggests that the probability to find heavy di-quarks in a plasma may not be too different from that of finding bound heavy quark-antiquark pairs. We turn now to the second line of eq. (2.24), which yields Finally, for the third line of eq. (2.24) we obtain JHEP06(2018)034 (H.8) At this point we may use the formulae listed in appendix F in order to perform the small y expansion. We obtain and [Y · H(0) · ∇ Y +y · H(0) · ∇ y +Y · H(r) · ∇ Y −2y · H(r) · ∇ y ] D 8 2N c . (H.10) This pair of equations forms a system that we can diagonalize perturbatively, following the procedure of section 5.1. The relevant coefficients that enter eq. (5.7) are easily identified on the equations above. We then find that the evolution of the component of the density matrix that is close to the maximum entropy configuration is given by Apart from keeping the dynamics of the center of mass explicit, this equation is very similar to eq. (5.8).