The approach to equilibrium of a quarkonium in a quark-gluon plasma

We derive equations of motion for the reduced density matrix of a heavy quarkonium in contact with a quark-gluon plasma in thermal equilibrium. These equations allow in particular a proper treatment of the regime when the temperature of the plasma is comparable to the binding energy of the quarkonium. These equations are used to study how the quarkonium approaches equilibrium with the plasma, and we discuss the corresponding entropy increase, or free energy decrease, depending on the temperature regime. The effect of collisions can be accounted for by the generalization of the imaginary potential introduced in previous studies, and from which collision rates are derived. An important outcome of the present analysis is that this imaginary potential has a sizeable dependence on the energy of the relevant transitions.


Introduction
There is an ongoing major effort to measure the production of heavy quark bound states in heavy ion experiments (for a recent review, see [1]). The goal of such measurements is to obtain information on the medium in which these heavy quark systems evolve. However, to achieve such a goal, we need to have good control of the dynamics of heavy quarks in a plasma, which is a difficult many-body problem. Different physical effects could play a role in modifying the properties of heavy quark bound states in a quark-gluon plasma, the most prominent ones being the screening of the binding forces and the collisions of the heavy quarks with the plasma constituents. The various models used in phenomenological analysis emphasize one aspect or the other, with of course many refinements in either direction. It is important however that all aspects of the dynamics be treated on the same footing, within a coherent formalism. Only then can one get confident that we understand the processes considered, and eventually extract from the data the properties of the medium in which the bound state evolves. Important progress in this direction has occurred in the last few years. A major step forward was the recognition that the effect of the collisions could be incorporated in an imaginary potential [2][3][4][5][6], somewhat analogous to the optical potential used in nuclear physics. This imaginary potential can then be calculated, albeit not yet with the same degree of accuracy as the real potential that is responsible for binding and is screened in a plasma. Attempts to access it via lattice calculations can be found for instance in [7,8]. As for the real potential, effective field theories have been used to constrain it in some particular regimes [4][5][6]. It has also been realized that techniques borrowed from the theory of open quantum systems (see e.g. [9]) could offer a new perspective on these issues. In particular, the imaginary potential appears naturally in the construction of the operators of the Lindblad equation [10]. The stochastic potential used in some approaches in connection with a Schrödinger equation (see e.g. [11,12]) is also intimately related to this imaginary potential. As we shall see in this paper, the imaginary potential also directly enters the calculations of the relevant transition rates.
The present paper complements the study presented in Ref. [13]. There, a complete derivation of the equation of motion for the reduced density matrix has been given, under the assumption that the intrinsic dynamics of the heavy quarks is slow compared to that of the plasma. This assumption allowed us to reduce the equations of motion to equations of a Langevin type. This assumption is strictly valid in the regime of high temperature, where the effect of binding forces are small and can be incorporated in the Langevin dynamics. The results that had been obtained along the same lines in the abelian case in [14] suggest that it is in this case a reasonable approximation, even when bound states can form. However, this is not so in QCD: when a quarkonium absorbs or emits a gluon its color state changes, from singlet to octet or vice versa. Since the potential between a quark and an antiquark is attractive in the singlet channel, and repulsive in the octet channel, the absorption of a gluon leads to a significant change in the effective heavy quark hamiltonian. This conflicts with some of the assumptions underlying the derivations presented in [13], which need therefore to be revisited. More broadly, we need to address more precisely the regime of moderate temperature where the binding energy is of the order of the temperature.
We consider in this paper a simplified set up, a static quark-gluon plasma in thermal equilibrium, and study the time evolution of a single heavy quarkonium in such a medium. The paper contains three main parts. In the next section, we derive equations of motion for the reduced density matrix for a quark-antiquark pair. These equations reproduce in some limit those obtained in [13], but they lend themselves to more accurate approximations in the regime where the temperature is of the order of the binding energy of the quarkonium. These equations are simplified by integrating out the center of mass coordinates, leaving us with equations for the relative motion alone. The second part of the paper, which covers Sects. 3 and 4, presents a general discussion of how the quarkonium approaches equilibrium with the quark-gluon plasma. We shall see that different treatments can be given depending on whether the temperature is large compared to the binding energy, or comparable to it. This will lead us to consider the variation with time of an (off-equilibrium) entropy and free energy. The third part of the paper, Sect. 5, presents some numerical calculations illustrating the main features of the general equations in some simplified situations. Conclusions are summarized at the end.

The evolution equation for the density matrix
We consider a single heavy quark-antiquark pair immersed in a plasma of light quarks and gluons in thermal equilibrium at a temperature T much smaller than the mass M of the heavy quark. The condition M T ensures that we can treat the heavy quark and antiquark as non-relativistic particles. Also, since the velocity of the heavy particles is small ( T /M ), we neglect their magnetic interactions (among themselves, and with the plasma constituents) 1 . We assume then that the whole system can be described by the following Hamiltonian where H pl is the QCD Hamiltonian governing the dynamics of the plasma while H Q controls the dynamics of the heavy quark-antiquark pair in the absence of the plasma. The hamiltonian H Q reads where r and R denote respectively the relative and the center of mass coordinates of the heavy particles. The interaction potential V s,o (r) is a function of the relative coordinates, and it depends also on the color configuration of the pair. Thus, as indicated in Eq. (2), we shall often write H Q as either H s or H o , depending on whether the quark-antiquark pair is in a color singlet (H s ) or a color octet (H o ) configuration. In leading order non-relativistic limit, i.e., keeping only the color Coulomb interaction into account, we have where C F = (N 2 c − 1)/(2N c ), with N c = 3 the number of colors, and α s is the strong coupling constant, α s = g 2 /(4π) with g the gauge coupling.
The last term in Eq. (1) is the interaction between the plasma and the heavy quarks. It is of the form 2 where a A 0 denotes the (color) Coulomb field created by the plasma particles, while n A denotes the color charge density of the heavy particles, with A a color index. For a quark-antiquark pair, the color charge density is given by wherer denotes the position operator 3 , and the two components of the tensor product refer respectively to the Hilbert space of the heavy quark (for the first component) and the heavy antiquark (for the second component). In Eq. (5), 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 .
2.1. The reduced density matrix and its color structure Consider now the density matrix D of the whole system. We assume that initially, at time t 0 , this density matrix factorizes where the plasma density matrix D pl (t 0 ) is an equilibrium density matrix at temperature T = 1/β: The reduced density matrix, D Q , the objet that we are mostly concerned with, is defined by taking the trace over the plasma degrees of freedom D Q (t) = Tr pl (D(t)).
The state of a heavy quark can be characterized by a position, a color, and a spin. We ignore here the spin degree of freedom. Then the reduced density matrix D Q has matrix elements of the form where a, b andā,b are color indices in the fundamental representation and its conjugate, respectively, while r i andr i (i = 1, 2) denote respectively the coordinates of the quark and the antiquark. Factorizing the color structure, one can write D Q as follows (see [13] for more details on the color structure of D Q ).
where D s and D o are matrices in the 2 particle space, with only coordinates as entries, e.g. the matrix elements of D s are r 1 ,r 1 |D s |r 2 ,r 2 . In the formula above, T F = 1/2 and the color matrices are normalized as Tr T A T B = δ AB /2. The relation between the first and second lines of Eq. (10) follows from the following formulae where |s and |o C denote respectively color singlet and octet (normalized) states, the index C in o C being a color index that distinguishes the various members of the octet. In the limit where the mass of the heavy quark is infinite, the density matrix is diagonal in coordinate space, and similarly for D o . In this limit the density matrix depends only on the relative coordinate r 1 −r 1 , which follows from the fact that the plasma in equilibrium is invariant under translations.

Approximate evolution equation for the reduced density matrix
The time evolution of the density matrix of the full system obeys the general equation of motion In order to treat the interaction between the plasma and the heavy particles by using perturbation theory, we move to the interaction representation with respect to the unperturbed hamiltonian H 0 = H Q + H pl and define where D I (t) satisfies the equation We can then rewrite the equation of motion (15) as This exact equation is obtained by formally integrating Eq. (15) and inserting the obtained 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. We can however improve on strict perturbation theory, with the help of two approximations. These are consistent with strict second order perturbation theory but go beyond, in particular by performing partial resummations (analogous to those in Schwinger-Dyson equations). The first approximation consists in replacing in the double commutator in the right hand side of Eq. (16) the density matrix by the factorized form This is consistent with second order perturbation theory since deviation from this form necessarily involves extra powers of H 1 . The factorization allows us to perform easily the average over the plasma degrees of freedom. We then obtain the following equation for the reduced density matrix of the heavy quarks We have used the fact that the linear term vanishes in a neutral plasma, and the sum over the color index A is implicit. Finally, we have written the correlator of the a 0 fields as Figure 1: These four diagrams are in one-to-one correspondence with the four terms in Eq. (18). The time flows as in a Schwinger-Keldysh contour, forward in the upper part, and backward in the lower part. The two upper lines represent the quark and the antiquark propagating from t to t, while the lower lines represent the same particles propagating from t to t .
The equation (18) can be given a simple diagrammatic interpretation, illustrated in Fig. 1 (see [13] for more details). The diagrams involve single gluon exchanges, represented by the correlators (19), with the gluon attached at points (t , x ) and (t, x). Contributions where the two densities are on the same side of the density matrix in Eq. (18), like in n A (t, x)n A (t , x )D I Q (t ), are associated to diagrams where the gluon joins lower or upper particle lines among themselves (the first and third diagrams in Fig. 1). Contributions where a density is lying on each side of D, as in in n A (t, x)D I Q (t )n A (t , x ), are represented by diagrams where the gluon joins upper and lower lines (the second and fourth diagram in Fig. 1).
The equation (18) contains a non trivial memory integral. However a second approximation allows us to obtain a Markovian equation. Indeed, we note that the difference D I Q (t ) − D I Q (t) involves powers of the interaction. Thus, at the order at which we are working, we can neglect this difference in the right hand side of Eq. (18), and simply substitute there D I Q (t ) with D I Q (t). At this point, the equation still contains a non trivial time integral, but it is Markovian. Moving back to the Schrödinger picture, one can write this equation as where we have made a change of variable in the time integration, and set t−t = τ . This Markovian equation is the equation that was studied in Ref. [13], and it will prove useful later on. In the rest of this section though, we shall use Eq. (18), which has a simpler diagrammatic interpretation. In fact, most of the derivations in this section are blind to this modification of the equation, as it will be clear. In the Schrödinger picture, Eq. (18) reads

Equation of motion for D s and D o
We shall write Eq. (21) as follows where L is to be understood as a linear operator acting on the density matrix. It corresponds typically to one gluon exchange processes (see Fig. 1), with the gluon being emitted at time t − τ and absorbed at time t.
Given the color structure of the density matrix (see Eq. (10)), it is convenient to view L as a matrix in the 2-dimensional space spanned by the two components D s and D o of the density matrix. Thus, we write In order to perform the color algebra needed to obtain the explicit expressions of the operators L ij , we note that both the density matrix and the heavy quark hamiltonian are diagonal in the singlet-octet basis (a property that we have already used in writing Eq. (23)). Furthermore, we note that the density operator n A (x) can connect singlet to octet states, and also various octet states among themselves. Its matrix elements are given by (see e.g. [13]) where The calculation then proceeds easily by inserting closure relations in the singletoctet basis at appropriate places in Eq. (21), for instance and using the following formulas to complete the color algebra It is then straightforward, by taking matrix elements of Eq. (21) in singlet or octet states, to obtain the expressions for the operators L ij in Eqs. (23). Thus, by taking matrix elements in singlet states, we get (with t ≡ t − τ ) In these equations, and similarly for ∆ < − (X, X ). This expression represents the combination of propagators that naturally emerges when one adds the four possible ways to hook the gluon in diagrams with a given topology (i.e. in one of the diagrams of Fig. 2). The minus sign in the last two terms finds its origin in the minus sign present in n(x) in Eq. (25) and affects the contributions where the gluon couples a quark and an antiquark. Furthermore, we have introduced the notation X = (t, X), where X represents the set of coordinates of the quark-antiquark pair, i.e., X = {x,x}, with coordinates without (with) bar giving the position of the quark (antiquark). The integral X in Eq. (29) runs over these coordinates, i.e., Finally, P X = |X X| is a projector, whose matrix elements between two localized states read  Figure 2: These four diagrams are in one-to-one correspondence with the four terms in Eqs. (28,29). In these diagrams the evolution operators concern the heavy quark pair (not the gluon). Each of these diagrams represents four similar diagrams where the gluon is hooked in four possible ways without changing the topological structure (these four contributions are summarized by the propagators ∆ < − and ∆ > − , cf. Eq. (30)). The first and third diagrams represent L ss , the second and fourth L so .
The evolution operators now depend on the color state of the propagating quark-antiquark pair. They are U o (τ ) = e −iHoτ for an octet state and U s (τ ) = e −iHsτ for a singlet state. The operator U (τ ) propagates the quark-antiquark pair forward in time, i.e., from t = t−τ to t, while its hermitian conjugate, U † (τ ) propagates the pair backward in time, from t to t − τ . The two operators are therefore attached respectively to the upper and lower pairs of lines in diagrams such as those introduced in Fig. 1. The structure of Eqs. (28, 29) may then be understood with the help of the diagrams displayed in Fig. 2.
For the operators L oj , we get and, writing L oo = L oo 1 + L oo 2 + L oo 3 , where and similarly for ∆ < + (X, X ). Note that there is no minus sign in ∆ > + (X, X ). This is because this contribution arises from products of factors m(x) in Eq. (25). As for the color factors, their origins can be easily traced back to Eqs. (24 -27).
The equations of motion that we have obtained for D s and D o are similar to those derived in Ref. [13], to within the small change discussed above (see after Eq. (20)), and the fact that in Ref. [13] the hamiltonian H 0 used in the interaction representation contains only the heavy quark and antiquark kinetic energy. At that point, in Ref. [13] a further approximation was performed, that consists in expanding the evolution operators at short time, i.e., writing U (τ ) 1 − iH 0 τ . Here we shall proceed differently in our treatment of the time integrals. The need to go beyond the approximation used in Ref. [13] is motivated in particular by the color changing transitions that take place in QCD: when a quark-antiquark pair in a singlet state absorbs a gluon, it turns into an octet state. This produces an immediate change in the effective hamiltonian of the pair, the force between the quark and the antiquark turning from attractive in the singlet state to repulsive in the octet state. It is then important to keep track of this change of the pair hamiltonian as the pair propagates during the lifetime of the exchange gluon. This is precisely what the various evolution operators do in the equations written above (see also Figs. 2 and 3), and why we have included the leading order interaction between the quark and the antiquark in the hamiltonian H Q .

Equations for the relative motion
At this stage, before doing any further approximation, we shall first simplify the equations that we have obtained by eliminating the center of mass coordinates. To that aim, we define a further reduced density matrix by taking a trace over the center of mass coordinates. With r andr denoting respectively the coordinates of the quark and the antiquark, we call R = (r +r)/2 the center of mass coordinate and s = r −r the relative coordinate. The reduced density matrix is defined from the matrix elements r 1r1 |D s,o |r 2r2 which, with a slight abuse of notation, we write also as R 1 , s 1 |D s,o |R 2 , s 2 . We callD s,o the reduced density matrix obtained after taking the trace over the center of mass coordinates. That is, The derivation of the equations of motion for the reduced density matrices D s,o is presented in Appendix Appendix A. It is assumed there that the center of mass velocity is small, typically T /M . We obtain then the following equations. For the singlet, we get and In these equations (see Appendix Appendix A) whereŝ is the operator measuring the relative coordinate.
For the octet, we get where we have set C q·ŝ ≡ 2 cos (q ·ŝ/2).
The equations above give a fairly complete account of the relative motion of a heavy quark-antiquark pair in a static quark-gluon plasma in thermal equilibrium. They are however difficult to solve in full generality. In order to get some familiarity with their physical content, we consider, in the next two sections the general question of how they describe the approach to thermal equilibrium, in two distinct regimes: The first regime is that of high temperature, controlled by the increase of the entropy. In the second regime, where the temperature is of the order of the binding energy, entropy effects compete with binding forces; there, the relevant quantity to look at is a non equilibrium free energy which will be seen to decrease monotonously as the equilibrium is approached. We analyze first the Abelian case, and in the following section we consider QCD.
From now on, in order to alleviate the notation, we omit the tilde inD s,o since we shall be dealing only with the reduced density matrix of the relative motion.

Entropy and free energy in the abelian case
A derivation completely analoguous to the one done in the previous section gives, for the case of an abelian plasma, with This equation is identical to the equation for D s (Eq. (23)) in which we replace C F by unity, set D o = D s = D, and use the expressions (38) and (39). Recall that D is the reduced density matrix of the quark-antiquark pair, after taking the trace over the center of mass coordinates. The operator S q·r is an operator in the space of the relative coordinates (cf. Eq. (40)), withr denoting the relative coordinate operator. Finally, in H Q , the potential is the ordinary Coulomb potential, V (r) = −α/r (analogous to V s in Eq. (3) from which it is deduced via the substitutions C F → 1 and α s → α, with α the fine-structure constant).
Our goal now is to perform, approximately, the integration over the time τ , in order to simplify the right hand side of the equation of motion. The strategies to do so differ, depending on whether we are in the high temperature limit of not, that is whether binding energy effects play an important role or not.

The entropy increase at high temperature
In this subsection we focus on the high temperature regime, i.e., the regime where the temperature is much higher than the binding energy. In this regime, one can ignore the binding energy, and the approach to equilibrium is dominated by entropy effects. We shall indeed show that the equations of motion predict a monotonous increase of the entropy. These equations of motion lead naturally to quantum Brownian motion, which was studied more extensively in [13]. In fact the equation for the reduced density matrix takes the form of a Lindblad equation, where the effect of the collisions are accounted for by an imaginary potential.
The temperature enters the equations through the gluon propagator, and limits the range of the τ -integration to τ 1/m D , with m D T . 4 The evolution operators U (τ ) contribute phase factors of the form e −iτ H , where H is the sum of a kinetic and a potential energy V . Collisions change the kinetic energy by an amount q 2 /M + P q/M , where q ∼ m D is the typical momentum of a gluon exchanged in a soft collision, and P ∼ αM the typical momentum of the heavy quark or antiquark in its relative motion in a Coulomb bound state. The contribution to the phase factor coming from the change of the kinetic energy is then of order m D /M 1 for the first term, and of order α for the second. In either case, these are small contributions that can be safely neglected. As for the binding force acting on the heavy quark, this can be estimated as follows: for a Coulomb bound state, we have parametrically, α/r ∼ αp ∼ M α 2 . Thus, when M α 2 m D , which can be satisfied at sufficiently high temperature, the potential energy V , and the binding energy, are small compared to the Debye mass, and one can safely set U (τ ) ≈ 1 in Eq. (46). 5 Along the same line, we assume that during the time interval between t − τ and t, the density matrix does not vary significantly so that we can replace D(t − τ ) by D(t) in the right hand side of Eq. (46) (alternatively, we could use the form (20) of the evolution equation for D). With these approximations the evolution equation greatly simplifies, and reads The integrand in the right hand side can be written as follows The contributions of the first and the second lines are qualitatively different: as we shall see shortly, the first line yields a correction to the real part of the potential, and corresponds to hamiltonian evolution, while the second line involves the imaginary part of the potential, and accounts in particular for dissipation. To see that, we perform the integration over q 0 , and then the integration over τ in Eq. (47). In the long time limit, t − t 0 1/m D , we can let the upper limit of the τ -integration go to infinity. The τ -integrations yield then (see [13]) The correction δV to the real part of the potential provides a contribution that adds up to the hamiltonian in the left hand side ( where we have used δV (−q) = δV (q). This is the screening correction to the real part of the potential. In the Hard Thermal Loop (HTL) approximation, this is given by (see e.g. [3]) so that the total potential in H Q is the screened potential V = (α s /r) exp(−m D r) (to within an irrelevant constant term). The remaining terms can be rearranged as follows This equation has the structure of a Linblad equation [15]. It can be written as where the Lindblad operators L q take the form ( At this point, we could rely on the theorem derived in [16] for the Lindbald equation, in order to show that the entropy is a monotonically increasing function of time. This theorem requires that , which clearly holds in our case since L † q = L q . However in order to highlight the difference between the present high temperature regime, and the low temperature regime to be discussed in the next subsection, we shall proceed with an explicit and elementary derivation.
From the definition of the von Neumann entropy, and using the fact that TrD(t)) is independent of time, a property that can be verified explicitly on Eq. (54), one easily obtains At this point, it is convenient to use a representation in which D(t) is diagonal, viz.
with p n (t) ≥ 0, and the states n are function of time. We then obtain This expression is manifestly positive, which implies that the entropy (56) indeed increases with time. As already emphasized, this proof is less general than the use of the theorem in [16] (it relies in particular on the property L † q = L q ).
It is interesting to relate the rate of entropy increase to our function W (r), or equivalently Γ(r). A crude estimate can be obtained as follows. First we note that therefore the combination (p n (t)−p m (t)) log pn(t) pm(t) that appears in Eq. (59) is a priori of the same size as S. Also, having in mind the approach to equilibrium of a system initially in a bound state, we may write (cf. Eq. (55)) where a 0 is a vector whose modulus coincides with the Bohr radius of the bound state, and A nm are constants of order unity. We can then estimate the order of magnitude of the change of entropy as follows where we have used W (q) = −g 2 ∆ > (0, q), and Γ(r) = W (r)−W (0), with W (r) the imaginary part of the potential. This estimate relates the rate of entropy increase to the imaginary part of the potential, that is to a typical collision rate, at a scale determined by the size of the bound state.

Free energy minimization in the Abelian limit
Now we look at the regime in which the temperature and the binding energies are of the same order of magnitude. It is no longer legitimate to approximate the evolution operators by unit operators, as we did in the previous subsection. In this case, the effects of entropy and binding compete. Then, a relevant object to look at is the non equilibrium generalization of the free energy. To analyze the time dependence of this quantity one is naturally led to expand on a complete set of eigenstates of the hamiltonian. The equations of motion lead in this case to rate equations, with rates that identify with those obtained from Fermi golden's rule. One can still express the effect of the collisions, at least partially, through the imaginary part of a potential, provided one takes into account the energy dependence of this potential.
We start again from Eq. (45), with the right hand side written as in Eq. (20), that is, In order to handle more easily the evolution operators U (τ ) = e −iH Q τ and U † (τ ), we introduce at appropriate places projectors on eigenstates of H Q , P n = |n n|, and assume for simplicity absence of degeneracy. Note that, in contrast to the previous subsection, the states n are now independent of time. We get The integration over τ can then be performed, using with → 0 + . We get In order to perform the integration over q 0 we note that where we have set E kn ≡ E k − E n , and the symbol P in front of the integral denotes the principal value. Similarly, We where the last equality follows from the rotational invariance of the plasma. We can then rewrite Eq. (66) as In the case where the typical transitions involve energy differences that are small compared to the Debye mass, which controls the decay of ∆(q 0 ) with q 0 , we can ignore the energies E kn , and perform freely the sums over n and k, which eliminates the projectors. Using the identities one then easily recovers the result of the previous section, i.e., Eq. (53).
We return now to Eq. (69). In order to minimize the effects of the principal parts and focus on the dissipative part of the equation, we use the eigenstates of H Q instead of H Q , and accordingly subtract the corresponding contribution of δV in the right hand side of the equation. We get where the energies E n are the eigenvalues of H Q . A this point, it is convenient to consider the explicit matrix elements of D and write the equation in a Liouvillian form In order to simplify the writing, we set i|S q·r |j → S ij in the following. We then obtain In this equation, P denotes the principal part of the integral from which the contribution 1/q 0 is subtracted (cf. Eq. (71)).
Assuming that the Liouvillian in the right hand side of Eq. (72) can be treated as a perturbation, we expect the effect of this perturbation to be dominant when it connects pairs of states with comparable energy differences, that is, when |E ij | |E kl |. In particular, one expects the diagonal elements of the density matrix, i.e. the occupation probabilities of the various levels, for which E ij = 0 to decouple from the non diagonal ones. 6 We now restrict ourselves to these diagonal matrix elements, and ignore possible degeneracies. It is easy to verify that the principal values then cancel. We get, for the case i = k, This is nothing but the decay rate Γ k→i calculated according to Fermi's golden rule. That is, where ∆ > (E ki , q) plays the role of the density of available states for the transition k → i. Similarly, for the case i = k we get 6 A more formal discussion of this issue is presented in Appendix Appendix B.
At this point we denote by p n (t) = n|D(t)|n the probability to find the system in the eigenstate n of H Q . The equation (72) yields then dp n dt = From the property and Eq. (75), it follows that Thus, in equilibrium where dp n /dt = 0, the detailed balance relation p k Γ k→n = p n Γ n→k implies that is p n ∝ e −En/T . In other words, the system reaches thermal equilibrium at the temperature of the plasma.
In order to see globally how the equilibrium is achieved, we look at the free energy F , defined in terms of the density matrix as in equilibrium, viz.
where in the last equality, the states n are the eigenstates of H Q . Taking the time derivative of this equation, and using Eqs. (77) and (78), we get Since all the terms in the sum are positive, this equation shows that the free energy is a monotonously decreasing function of time (at least in the large time limit). An alternative and more formal proof can be obtained by using Lemma 1 of Ref. [17]. Note that if we use this evolution equation to compute the derivatives of S and E separately, one obtains expressions that do not have necessarily a well defined sign. We may estimate the rate of change in the free energy, using a similar argument as that used for the entropy. We get where ∆E is a quantity that represents an average value for the binding energy differences. In the last line Γ(∆E, a 0 ) is a damping rate that summarizes the effect of the collisions. This can be viewed as a definition that generalizes that of Γ(a 0 ) to include an energy dependence (i.e. Γ(a 0 ) = Γ(∆E = 0, a 0 ), see last section). Eq. (83) is essentially the same as that giving the entropy increase, Eq. (62), with ∆ > (0, q) replaced by ∆ > (∆E, q). We shall discuss in the last section how large is the correction due to this energy dependence.

Entropy and free energy in a non-Abelian theory
In this section we repeat the analysis of the previous section in the case of QCD. The generalization is straightforward except for obvious complications related to the color algebra, and the existence of several components of the density matrix.

Entropy increase
As in the abelian case, in the high temperature limit, the evolution equations for the reduced density matrix are a set of Lindblad equations, in which the effect of collisions is taken into account via an imaginary potential. We shall verify that the entropy increases monotonously as the equilibrium is approached.
Following the same reasoning as in Sect. 3.1 we get the following simplified expressions for the various operators L: and In the previous equations we have used L q as defined in Eq. (55) whilē In the QCD case the entropy can be written where the factor N 2 c − 1 is due to the normalization chosen in Eq. (10). By using the explicit expression of the operators L ij just given above, we can write the derivative of the entropy as Using exactly the same reasoning as in Sect where, in the last expression, |o, m o, m| actually stands for Then, a simple calculation allows us to write the first line of the right hand side of Eq. (92) as In conclusion, all the terms contributing to the derivative of the entropy are positive. This implies that in the regime where the temperature of the quarkgluon plasma is large in comparison to the typical binding energies, the equations for the reduced density matrix yield a monotonous increases of the entropy as the quarkonium approaches thermal equlibrium. The rate of entropy change can be estimated in the same way as we did for the abelian case in Sect. 3.1.

Free energy minimization
With consider now the regime of moderate temperatures, and will proceed to the calculation of the free energy. We shall first write the necessary rate equations describing the evolution of the populations of the various states. Although many of these states belong to a continuum, we write the summations over states as discrete sums, as in Eq. (94), instead of integrations, since we focus here on the general structure of the equations. The contimuum states will be explicitly dealt with in the examples treated in the next section.
The probabilities p s n and p o m fulfil the following evolution equations dp s and dp o Note that in order to obtain Eq. (97), we had to combine Eq. (84) and (85) The first equation gives the transition rate Γ o→s from one particular member of an octet state to a singlet state (the factor 1/(2N c ) follows from Eq. (24)).
In the second equation, giving Γ s→o , all members of the considered octet are summed over (producing a factor N 2 c − 1)). Similarly, we have, for the octet to octet transitions Γ (2) o,m→o,k = o,m→o,k = All these transition rates are those which control the evolution of the populations, according to the equations written above. With these ingredients we can compute the evolution of the free energy. Expanding on the basis of the eigenstates of H Q , and dropping the on the energies in order to simplify the notation, one gets Using analogous manipulations as the ones we used in section 3.2 we get Again the physical interpretation is straightforward: each of the elementary transitions makes the free energy decrease separately, in a way that is very similar to what happens in the Abelian limit.

Some illustrative calculations
In this last section, we present results of some numerical solutions of the equations for the density matrix in simplified situations. We emphasize that our main goal here is to illustrate some of the concepts that we have introduced. Thus, although the numbers are adjusted to bottomonium physics, we make no attempt to a complete phenemonological description. The first example that we treat is that of an infinitely massive quark-antiquark pair. This provides a simple illustration of the role of the energy dependence of the imaginary potential, as well as a quantitative indication of the magnitude of the effect. In principle, such a setting corresponds to that used in lattice QCD calculations, and we briefly compare with relevant lattice results. The next example involves a simplified picture of a bottomonium in a plasma, with a single bound state in the singlet channel, and octets states involving the free quark and antiquark. In this case, rate equations are complemented by a Langevin equation describing the Brownian motion of the heavy quark and antiquark in the plasma.

Infinitely massive quark-antiquark pair
The physics of singlet to octet transitions is best analyzed by ignoring completely the motion of the heavy particles and focussing on their color degrees of freedom alone. This is what we do in this subsection.
We consider an infinitely massive quark-antiquark pair, and assume that it can exist in two color states, a singlet (s) and N 2 c −1 degenerate octet (o) states. There are no continuum states, so that the problem reduces to that a two level system. The partition function reads where V s (r) and V o (r) denote the energies of the pair in a singlet and an octet state, respectively. Since the particles do not move, they have no kinetic energy, and the energy of the pair is just the potential energy, which depends on the distance r between the quark and the antiquark. We assume that V s < V o , and set ∆V = V o − V s . The free energy is given by F = E − T S = −T ln Z, with E the internal energy and S the entropy. The latter can be deduced from F by using the thermodynamic relation S = −∂F/∂T . In the low temperature limit, In this regime, the dominant contribution to the free energy is the energy of the ground state V s , the correction from the octet excited states being exponentially small, ∝ T e − ∆V T . The internal energy and entropy are given by In the opposite limit of high temperature, the free energy is entirely dominated by the entropy. A simple calculation yields indeed The factor N 2 c is just the total number of available states, one singlet and (N 2 c − 1) octet. All are present in equilibrium with the same probability. The Boltzmann factors in the partition function can all be approximated by unity, and the internal energy is simply given by It is independent of the temperature. Let p s and p o be the probabilities to find the system respectively in the singlet or a given octet state. In the infinite mass limit, these are simply the diagonal elements of the density matrix (cf. Eq. (12)), i.e., p s = D s (r) and p o = D o (r). In equilibrium, we have These probabilities depend only on the ratio ∆V /T , which controls the transition between the low and the high temperature regimes: At low temperature, T ∆V , and p s 1. At high temperature, T ∆V , In the numerical calculations, we use where m D is the HTL Debye mass calculated with a running coupling at the scale 2πT .
We are interested in the dynamics of the approach to the equilibrium state. As we have seen in the previous sections, this is controlled by a rate equation, generically of the form where Γ o→s denotes the transition rate from any one of N 2 c − 1 degenerate octet states, and similarly for Γ s→o . In a stationary state, the rate equation dp s /dt = 0 yields the detailed balance condition This is to be compared to the result that we expect when the stationary state is the state of thermal equilibrium (cf. Eq. (109)) By comparing the two equations (113) and (112) This relation is satisfied by the rates that we have obtained in the previous section. It implies in particular that their energy dependence needs to be taken into account when the temperature is of the order of magnitude of ∆V : in that case the static imaginary potential is not sufficient to fully account for the effects of collisions. We return to this issue shortly. The evolution equation (111) has the following solution, for arbitrary initial conditions, and where p eq s and p eq o are the equilibrium values, given in Eq. (108), andΓ is an effective rate defined asΓ The solution p s (t) and p o (t) obtained for p s (0) = 1 and r = 0.15 fm are plotted in Fig. 4, together with the free energy calculated from Eq. (101). For both the survival probability p s (t) and the free energy, the effective rateΓ(r) determines the time scale that controls the approach to equilibrium. In order to quantify the importance of the energy dependence of the rates, we use the explicit results that we have obtained in the previous section. From Eq. (97), we get This equation is identical to Eq. (111), with now the following explicit expressions for the rates It is easily verified that these rates satisfy Eq. (114), as stated above.
In the static limit, and at high temperature, one can express the survival probability of the singlet state in terms of an imaginary potential. We have where we have used W (q) = −g 2 ∆ < (0, q), and in the last line Γ s→o is given by Eq. (119) in which one sets ∆V = 0. This relation suggests the following definition of a generalized, energy dependent, "potential" W (ω, r), viz.
where σ(ω, q) = ∆ > (ω, q) − ∆ < (ω, q) is the (longitudinal) gluon spectral function, and we have used the relation Note that as ω → 0, W (ω, r) reduces to W (r) since, in this limit, N (ω) ∼ T /ω, and T σ(ω, q)/ω → ∆ < (0, q). It is perhaps useful to recall here a few basic properties of the gluon spectral function σ(ω, q). To be specific, we shall rely on the HTL approximation, for which an explicit expression is known (see e.g. [3]) 7 . At fixed momentum, σ(ω, q) is an increasing function of the energy (linear at small energy), in the space-like domain |ω| < |q|. For |ω| > |q| it vanishes, except for an isolated delta-function contribution corresponding to the plasmon excitation at ω q (ω 2 q ω 2 pl + 6q 2 /5, with ω pl = m D / √ 3), which exists only for |q| m D . The specific contribution of the plasmon to Eq. (121) will be ignored in this paper. 8 For ω = 0, we know σ(ω, q), and hence ∆ < (0, q), analytically. This is so that [2] When the energy is non vanishing, the expression of σ(ω, q) is more complicated. It can be obtained from the analytic propagator (see e.g. [21]) The imaginary part of Π L (obtained by setting ω → ω + iη, with ω real) determines the continuum part of the spectral function at small energy. It is given by More generally, we have (128) 7 Note that all the numerical calculations presented in this paper use this approximation. Note also that we shall be using this approximation beyond its strict regime of validity, which requires ω, q T . 8 It should of course be included in a more quantitative study. It represents a process analogous to that of gluon dissociation involving the transverse modes of the gluon [18][19][20]. The collective plasmon exists only at small momentum q m D , and its contribution to W (ω, r) − W (ω, 0) is expected to be modest in the region of interest, and taking it into account would not alter the main conclusions of this section.
Note that the temperature enters the spectral function only through the Debye mass m D , and we can set σ(ω, q) = m −2 Dσ (ω/m D , q/m D ), whereσ(ω/m D , q/m D ) is a dimensionless function. On the other hand, the statistical factor that multiplies σ(ω, q) in Eq. (122) depends only on T . Knowing the spectral function, one can then determine the energy dependent potential (121). The results of this calculation are displayed in Fig. 5 for two values of the temperature, T = m D and T = 5m D . One sees there that the dominant effect of the energy dependence is a sizeable reduction of the imaginary potential, a reduction which gets amplified as the temperature, when it is of the order of the energy, decreases. This reduction arises from the fact that as the energy ω of the transitions increases the phase space of the space-like gluons that induce such transitions decreases. The density of such gluons in momentum space is essentially the quantity ∆ < (ω, q) and its decrease with increasing ω, for a given q, results from the combination of two effects: the statistical factor suppresses transitions with ω > T , and the spectral density vanishes when ω > q.
To proceed further, it is convenient to rewrite Eq. (121) in terms of dimensionless variables, as follows The curves in Fig. 5 are obtained after integration over q, which affects the dependence on r of W (ω, r) at fixed ω. In particular, at large values of rm D , the last term in the integral (129) yields a vanishing contribution, which is the origin of the flat behavior observed in Fig. 5. Another factor determines the r dependence of the rates: the energy ω is to be set equal to ∆V (r) (see Eqs. (119)). It turns out that, after this substitution, the dominant r-dependence, in the relevant range, is captured by the function h in Eq. (129), that is where, after reinstating the appropriate color factor C F , we have set and Γ s (r) is given explicitly in Eq. (120).

The suppression factor h ∆V (r)
T is plotted as a function of r for different temperatures in Fig. 6. Note that h(x) → 1 as x → 0, while, h(x) ∼ e −x as x → ∞. Thus, the strong suppression at small r originates from the fact that ∆V (r) → ∞ as r → 0, that is, h(r) ∼ exp{−∆V (r)/T }. This overwhelms the suppression already present in Γ s (r) ∼ r 2 ln(1/r), reflecting color transparency, i.e. the suppression of interactions when the size of the color dipole made by the quark and the antiquark in a color singlet vanishes. At large r, ∆V (r) → 0, and h(r) ∼ 1 − ∆V (r)/(2T ).
The setup discussed in this subsection is very close to that used in lattice QCD calculations of the potential, or free energy, of a heavy quark-antiquark pair. In particular, we may attempt a comparison with the recent results of Ref. [8]. Since the potential calculated there is reconstructed from the spectral function, it should contain, in principle, the energy dependence that we have been discussing. A comparison with the lattice results show that, as is the case with the imaginary potential that we calculate, the small r dependence is clearly different from the behavior (∼ −r 2 log(r)) expected in the absence of energy dependence: there is a strong suppression at small distances which carries on up to larger radius as the temperature decreases. Unfortunately, a more quantitative comparison between our computation and the lattice simulations is difficult, since the Coulomb approximation that we use is not accurate at large distance, and an additional effect due to the string tension cannot be excluded, as discussed in [8].
Finally, we return to the estimates of the rate of entropy or free energy variations for which expressions were derived in Sect.3. The explicit expression of Γ(r) in the HTL approximation is given in Eq. (124). Using this approximation we estimate that at T = 250 MeV, and using as Bohr radius a 0 = 1 1200 MeV −1 the time scale that characterizes the changes in the entropy is around 12 fm/c. This is of the order of the typical total lifetime of the fireball produced in heavy ion collisions. This suggests that the state of quarkonium will be substantially modified but perhaps not fully thermalized. Note also that, according to Fig.7, the energy dependence reduces the rate by about a factor 3.

A simplified model of quarkonium evolution
We now move away from the static limit and consider a more "realistic" model for the quarkonium. This is based on the following assumptions: • We neglect the quark-antiquark interaction in the octet channel, i.e., we set V o = 0. This approximation would be justified in a large N c limit, and it was used in the original derivation of the gluon-dissociation crosssection [18]. Comparison with later derivations shows that it remains a reasonable approximation for N c = 3 [19,20].
• This implies in particular that the heavy particles behave as free particles once they are in the octet channel, and octet to octet transitions can be treated in the high temperature limit. In this case, the corresponding equation of motion for the heavy quarks will reduce to a Langevin type equation.
• We assume that a single bound state exists in the singlet channel. The survival probablity of this singlet bound state is entirely controlled by its interaction with octet (continuum) states. In fact, we also ignore continuum singlet states (these represent only about 10% of the available continuum states).
In summary the model that we consider consists in one bound state, the singlet state, and a continuum of free, octet, states. We want to see how our equations describe the approach to equilibrium of this particular system.
We start with elementary remarks concerning the system when it is in thermal equilibrium with the plasma. Recall that we ignore the center of mass motion. We write the partition function of the relative motion as follows where |E s | is the binding energy (E s < 0), Ω is the volume of the plasma, and λ T = 4π/M T is the thermal wavelength of the relative motion. At low temperature, the bound state dominates, and Z ≈ e |Es|/T , while at high temperature Z ≈ Ω(N 2 c − 1)/λ 3 T . Clearly, there is a transition temperature T c , of the order of |E s |, corresponding to the situation where these two contributions are of the same order of magnitude, Note that T c has a (weak, logarithmic) dependence on the volume, and would vanish in an infinite volume. Let p s and p o be the probabilities for the system to be in the ground state or in a continuum state, respectively. We have Clearly, when T T c , p s ≈ 1, while p o ≈ 1 when T T c . The time evolution of the probabilities are given by the simplified system of equations and where E o p = p 2 /M is the kinetic energy of the relative motion. We assume that the medium is contained in a cubic box of volume Ω. Computations are made for two different volumes, Ω = 1 fm 3 and Ω = 100 fm 3 (these values cover the orders of magnitude of the typical volumes of the fireballs produced in a heavy-ion collision). This volume factor affects the numerical results, and it has been made explicit in Eq. (137). Thus, the plane wave in the equation above is normalized so that r|o, p = e ir·p .
The first equation expresses the change in the bound state population, with a loss term caused by the singlet to octet transitions, while the gain term represents the possible transitions of any of the continuum octet states to the bound singlet. The second equation accounts in addition for the Brownian motion of the particles in the continuum. The specific form of the Langevin terms in the left hand side is taken from Ref. [13]. As a simple consistency check of these equations, one may verify that and that the steady state solution is given by Eqs. (135).
We have solved Eqs. (136, 137), taking for γ the value used in Ref. [14], but adapted to the bottomonium mass, assuming that γ goes as the inverse of the mass, that is, γ = 0.060 fm −1 . Other needed inputs are the binding energy and the wave function of the singlet ground state (that enters the computation of the matrix element s|S q·r |o, p ). We obtain these by solving the Schrödinger equation with a screened potential 9 The results are shown in table 1. As can be seen in this table, screening sub-   stantially reduces the binding energy. Note that this reduction of the binding energy, together with the corresponding modification of the singlet wave function, entail a substantial increase of the decay width at a given temperature. This is of course in line with the energy dependence of the rates that we analyzed in the previous section. At T = 400 MeV, the decay width is of the same magnitude as the binding energy, suggesting that at this temperature the singlet can hardly be considered as a bound state anymore.
To solve these equations we use the same numerical methods as in Sect. 5.4 of Ref. [13]. The most relevant difference as compared to the case treated in [13] is that, in the present case, the singlet bound state can decay into octets with different momenta. To include this feature in our simulation, we use a rejection sampling based on the differential decay width to select the momentum.
The value of the survival probability of the singlet state, p s , obtained by solving these equations is given in the following table, for two different temperatures, two different interaction times, and two different volumes. For comparison the value of the equilibrium probability is indicated in the last column. We see that at late times p s becomes very small: the system is then completely dominated by the octets, the more so the larger the volume. However, on time scales that are typical of the lifetime of the plasma in a heavy-ion collision, ∼ 5 fm/c a significant amount of singlets survive, the survival probability being essentially independent of the volume.
In order to get a feeling for the role of the Brownian motion of the free quarks, we have repeated the calculations dropping the Langevin terms in the left hand side of Eq. (142), i.e., keeping only the time derivative. We obtain then the results listed in the table below:  As can be seen from this table, the survival probability is much reduced, and it eventually vanishes at large time: the absence of an energetic penalty for the transition to an octet state allows for a rapid occupation of the large continuum phase space. This provides another indication of the importance of the energy dependence of the rates. As a final remark, we have evaluated the free energy, internal energy and the entropy in equilibrium. These are shown in the following table In all cases, the free energy is dominated by the entropy, the more so the larger T and/or Ω. Note that in the large volume Ω limit < E > goes to a constant, while F ∼ −T S goes to infinity but just as -T log Ω.

Conclusions
The equations for the reduced density matrix that we have derived in this paper describe the evolution of a quarkonium towards thermal equilibrium in both regimes of high and moderate temperatures. The high temperature regime is that where the binding energies can be ignored. Then the dynamics is well described by a Lindblad equation and the approach to equilibrium is controlled by the increase of the von Neumann entropy. In this regime, binding forces can be treated perturbatively, and the effects of the collisions accounted for by a static imaginary potential. In the regime of moderate temperature, binding energies cannot be ignored, and it is convenient to use as a reference basis, that of the eigenstates of the effective heavy quark hamiltonian. The equation for the density matrix then leads to rate equations describing transitions between these eigenstates, and the approach to equilibrium is accompanied by the decrease of a free energy. The dynamics of continuum states remains dominated by Brownian motion and is described by a Langevin equation. In this regime of moderate temperature, the effect of collisions is still captured by an imaginary potential, which enters the determination of the collision rates. An important feature of the imaginary potential is that it depends on the energy: this is because as the energy of a transition increases, the phase space of the space-like gluons that cause this transition decreases. This energy dependence is found to be numerically important and should be taken into account in phenomenological studies. As we have emphasized, this effect is expected to be much more important for QCD than it would be for QED. This is because the absorption of a gluon by a quark-antiquark pair changes the color state of the pair, and turns an attractive force into a repulsive one, or vice versa.
The last section of the paper presented numerical studies illustrating the main concepts discussed in the earlier sections. The first example is that of a pair of infinitely massive quark and antiquark. This is close to the typical set up used in lattice gauge calculations, and some comparison with recent lattice results has been attempted. It would certainly be worthwhile to extend such comparison and see whether the strong suppression arising from the energy dependence of the imaginary part of the potential at small separation can be reproduced by lattice calculations. The second example treats a simplified model of a quarkonium with a single bound state in a singlet state, and continuum octet states. Although this is an oversimplification of the realistic situation, many interesting features emerge from this study, that could be relevant in phenomenological studies. This example illustrates in particular the interplay between screening and collisions, and the importance of treating both on the same footing, as we do in this paper. 10 In this paper we have focussed on a simple question, how a quarkonium approaches equilibrium when it is in contact with a static quark-gluon plasma in thermal equilibrium at temperature T . Although we have examined only simplified models, the equations that we have derived allow in principle for a quantitative answer to this question. They should provide a consistent starting point for more elaborate phenomenological work. The formalism developed in this paper should be well suited to the study of the bottomonium, presumably to be found in the moderate temperature regime in heavy ion collisions. The case of charmonium is more intricate, and presumably calls for a mix of techniques, in particular if one wishes to address the issue of recombination. Then the approximations developed in [13] may be useful.
was supported by the Academy of Finland, project 303756.

Appendix A. Elimination of the center of mass coordinates
In this appendix, we perform the Fourier transform of our main equation, and eliminate the center of mass coordinate.
In order to proceed with the Fourier transform, we note that the correlator ∆ < − (X, X ) defined in Eq. (30) depends on the difference of times, τ = t − t , and a priori on 4 coordinates. Because of translation invariance, it is in fact function of only three differences of coordinates. To make this more explicit, we consider ∆ < − (X, X ) as an operator in the two particle space, with matrix elements The coordinates Y and y play an important role in the semi-classical approximation (see [13]). Thus, in the large mass limit, r 1 ≈ r 2 andr 1 ≈r 2 , Y → 0, y → 0 and r becomes equal to the relative coordinate. It is convenient to express the Fourier transforms in terms of these variables. We have, for instance (with the shorthand notation q = and for ∆ < − (τ ; Y , y, r), The variable q has the meaning of the momentum of the exchange gluon. We can also take a Fourier transform with respect to the time variable (q = (q 0 , q)) and write (with now q = We shall sometimes write, with a slight abuse of notation, ∆ − (q, s 1 , s 2 ) in place of ∆ − (q, y, r). To summarize, we can write r 1 ,r 1 |∆ < ± (τ, X, X )|r 2 ,r 2 as where q = (q 0 , q) has the interpretation of the exchanged gluon four-momentum, and These relations allow us to perform the partial trace over center of mass coordinate, or equivalently over the center of mass momentum. We illustrate the procedure with the first term of Eq. (23), and more precisely the first contribution to Eq. (28). Consider then the matrix element between localized states of the heavy quark antiquark pair (in a color singlet state). Taking advantage of the fact that the projectors are diagonal in coordinate space, we can rewrite this as (omitting to indicate time variables to simplify the writing) where X i = (R i , s i ). We note then that the evolution operators factorize into a center of mass contribution which depends only on the kinetic energy of the center of mass, and a part related to the relative motion that involves the potentials V s or V o . We set U s = U cmŨs and similarly for U o . We have andŝ is the relative coordinate operator. At this point it is (almost) trivial to trace out the center of mass degrees of freedom. This amounts to set P 1 = P 2 and to integrate over P 1 . Note that the commutator in Eq. (23) yields The first term will not contribute when taking the trace (with P 1 = P 2 ). As for the second term, it yields s 1 |[H s ,D s ]|s 2 .
To proceed further we need to analyze the phase factor in Eq. (A. 16). This is the product of τ by the recoil energy where we have set P 1 = P 2 = P . The quantity ∆E recoil is the recoil energy of the heavy quark system, after absorption or emission of a gluon with momentum q. The range of the τ -integration in Eq. (A.16) is limited by the propagator ∆ > (τ, q) to be of the order of the inverse Debye mass m D T . On the other hand, the collisions of the heavy particles with the light constituents of the plasma involve the exchange of soft gluons, with |q| m D . It follows that typically, qτ ∼ 1, and the recoil energy is a small fraction of the Debye mass, τ q 2 /M (m D /M ). A similar estimate holds for the term τ P · q/M (T /M ), where we have assumed that P T (we consider pairs that are initially at rest. If the center of mass momentum is high then we need to consider "hot wind" effects, which is beyond the scope of this paper [24,25]). Since we assume that both m D M and T M , one can safely ignore the phase factor. A similar reasoning can be made for all the contributions to the main equations. We then obtain the equations that are listed in the main text.

Appendix B. Multiple-scale analysis
In this appendix, we discuss the solution of Eq. (72) within perturbation theory, paying particular attention to the secular terms. We first rewrite Eq. The difficulty with this naive expansion is that the condition D 1 D 0 is not always satisfied. In particular, this condition is violated at late times if the following quantity is not equal to zero. Let us then set This expression makes explicit the secular term, growing linearly with time, at the origin of the breakdown of naive perturbation theory. The problem can be handled by multiple-scale analysis (see e.g. chapter 11 of [26]). One introduces a "slow" time τ = t, and consider D as a function of t and τ , treated (artificially) as independent variables. We have, as before, D(t, τ ) = D 0 (t, τ ) + D 1 (t, τ ) + · · · (B.8) The leading order equation reads The next to leading order equation involves the derivative of D 0 (t, τ ) with respect to τ , viz.
We can now use ∂D0 ∂τ in order to cancel the secular contribution of F[D 0 (t, τ )], that is, we demand that the following equation By construction, D 1 (t, τ ) non longer contains a secular term, and can be considered a genuine perturbative quantity. A similar result could have been obtained by applying renormalization group techniques, as discussed recently in Ref. [27]. The separation of the secular term requires the solution of Eq. (B.12). By projecting this equation on the eigenvectors of the operator H, assuming that all energy levels are discrete, we obtain This is non-zero only if E n − E m = E n − E m . Thus the evolution described by Eq. (B.12) only connects pairs of states whose energy differences E mn = E m −E n are identical. It follows in particular that the evolution of the populations of the various quarkonium states, i.e., of the diagonal elements, for which E mn = 0, decouples from that of the non-diagonal ones (at leading order and assuming absence of degenerate states). One can also evaluate similarly the matrix elements of D 1 (t, τ ). One gets where, by construction, E n m = E nm . Assuming that the particles are confined in a volume Ω ∼ L 3 , the lowest values of the energy denominators are of order L −1 . Thus, if Γ denotes the leading order decay rate, D 1 will remain a small perturbative correction as long as ΓL 1. In the example treated in Sect. 5.2, this condition is well satisfied for Ω = 1 fm 3 , but only marginally for Ω = 100 fm 3 .