Heavy quarkonium suppression in a fireball

We perform a comprehensive study of the time evolution of heavy-quarkonium states in an expanding hot QCD medium by implementing effective field theory techniques in the framework of open quantum systems. The formalism incorporates quarkonium production and its subsequent evolution in the fireball including quarkonium dissociation and recombination. We consider a fireball with a local temperature that is much smaller than the inverse size of the quarkonium and much larger than its binding energy. The calculation is performed at an accuracy that is leading-order in the heavy-quark density expansion and next-to-leading order in the multipole expansion. Within this accuracy, for a smooth variation of the temperature and large times, the evolution equation can be written as a Lindblad equation. We solve the Lindblad equation numerically both for a weakly-coupled quark-gluon plasma and a strongly-coupled medium. As an application, we compute the nuclear modification factor for the $\Upsilon(1S)$ and $\Upsilon(2S)$ states. We also consider the case of static quarks, which can be solved analytically. Our study fulfils three essential conditions: it conserves the total number of heavy quarks, it accounts for the non-Abelian nature of QCD and it avoids classical approximations.


I. DILEPTON EMISSION FROM QUARKONIUM
The main way in which heavy quarkonium is detected in heavy-ion collisions is through its decay into a lepton pair. Electromagnetic interactions are slow compared with the strong interactions that drive the dynamics of the fireball; therefore, the bulk of these decays will happen after freezeout. An observation that supports this understanding is that the positions of the peaks in the dilepton emission spectrum are the same in pp and AA collisions [1].
The decay rate into leptons in thermal equilibrium was computed in [2]. Here we generalize that result to the case of a medium that is not in thermal equilibrium. The Hamiltonian of the system can be written as H ¼ H QCD þ H EW , where H QCD is the QCD Hamiltonian, and H EW is the part of the Hamiltonian that includes leptons and the electroweak interaction (for the following use, only the electromagnetic part of the electroweak interaction is relevant). At some early time t 0 we have a system that only contains quarks and gluons. It is described by a density matrix ρðt 0 Þ ¼ P nm ρ nm ðt 0 Þjnihmj, where jni are eigenstates of H QCD . At a much later time, the system is made by an arbitrary number of quarks and gluons and a lepton pair, l þ l − . Hence, up to higher orders in the electromagnetic coupling, α, the state describing such a system is the product of a QCD state, jni, and a lepton-pair state, jl þ l − i: jn; l þ l − i ≡ jnijl þ l − i. Naming k 1 and k 2 the momenta of the two outgoing leptons, the differential emission rate is At leading order in α, we can write hj; l þs ðk 1 Þl −s 0 ðk 2 Þjni where j0i EW is the vacuum of the electroweak theory, a s k 1 and b s 0 k 2 are the lepton and antilepton annihilation operators, l is the lepton field, A μ the photon field, J μ the quark electromagnetic current, and e Q the fraction of electron charge e carried by the heavy quark. Summing over the polarizations, we obtain dR ¼ − e 4 e 2 Q L μν ðk 1 ; k 2 Þ jk 1 where L μν ðk 1 ; For ρ ¼ e −H QCD =T , where T is a temperature, we recover the results of [2]. In general, the medium formed in a heavy-ion collision will be out of thermal equilibrium and characterized by some correlation lengths that change with time. In this work, in order to keep close contact with and take advantage of existing studies of the medium in thermal equilibrium, we will often assume that the medium evolves in time in a quasistatic fashion. This means that the medium is locally in thermal equilibrium and that at each time we can define a temperature T. This assumption allows us to take over results obtained in thermal equilibrium, with the only difference being that now the temperature depends on time. Our explicit model for the time dependence of the temperature will be discussed in Sec. II. A quantitative characterization of a temperature varying slowly (quasistatically) with time in the context of heavy quarkonium dissociation will be given at the end of Sec. III B. Nevertheless, it is worth stressing that the evolution equations that we will derive in Sec. III B do not rely on the quasistatic approximation and that in some of the following discussions, including the rest of this introduction, we may simply understand T as the inverse of a correlation length characterizing the system.
In order to write more explicitly Eq. (3), it is convenient to take advantage of the heavy-quark mass, M, being much larger than the typical momentum of the particles in the fireball, which is proportional to T 1 : One can therefore use nonrelativistic QCD (NRQCD) to describe the heavy-quark dynamics [3,4]. In NRQCD, heavy quarks are represented by a Pauli field ψ that annihilates the heavy quark and a Pauli field χ that creates the heavy antiquark. Up to corrections of order α s ðMÞ and T=M, the NRQCD electromagnetic current reads In the quarkonium rest frame, it is v ¼ ð1; 0Þ. J 0 does not contribute to the emission of the lepton pair with an invariant energy around 2M. This contribution comes only from J i . The Pauli matrix in J i projects onto the subspace of quarkonia with spin 1. In terms of the NRQCD heavyquark fields, the emission rate can be written as The right-hand side of the previous equation will depend on different energy scales, some coming from the thermal plasma, which depend on the temperature T, and some from the nonrelativistic nature of the heavyquark-antiquark bound state. These last ones are the typical heavy-quark-antiquark distance, a 0 , and the typical quark-antiquark binding energy, E. They respect the nonrelativistic hierarchy: The physics of heavy quarkonium in the fireball is going to depend on the relation between these two sets of scales. We will assume in the following that Under the condition (10), we can use potential NRQCD (pNRQCD), which provides a valid description of the quarkonium at energy scales below 1=a 0 [5][6][7][8][9]. In this effective field theory (EFT), the heavy-quark-antiquark system can be described in terms of a color-singlet field S and a color-octet field O instead of the fields ψ and χ of NRQCD. At leading order, where S l ≡ TrfSσ l g= ffiffi ffi 2 p . The result reflects the expectation that the emission of dileptons should be proportional to the number of singlet states with spin 1.
The integration of the temporal components t and t 0 may be split into where t F is the time of freeze-out. The argument that most of dilepton production from quarkonium happens in the vacuum amounts to having After freeze-out, the system evolves as it would in the vacuum. Hence for t, t 0 larger than t F , we have where ϕ n (ϕ m ) is the wave function of the nth (mth) state, E n (E m ) is the corresponding energy (including the decay width as an imaginary phase), and P (P 0 ) the center-of-mass momentum; ρ s ðt 0 ; tÞ will be defined in full generality later on [see (23)] and will be the central object of our study. Therefore under the approximation (14) and in the quarkonium rest frame, we have where we have defined q ¼ k 1 þ k 2 . The above equation tells us that the emission rate of dileptons with the energy of a given quarkonium state n [q 0 ¼ 2M þ ReðE n Þ þ q 2 =4M] will be proportional to the projection of ρ s ðt F ; t F Þ into that state at freeze-out time. The problem reduces to find ρ s ðt F ; t F Þ given some initial conditions at t ¼ 0, ρ s ð0; 0Þ. Within this framework, the yield of quarkonium nS states in AA collisions normalized with respect to the yield in pp collision is given by Note that the ratio above depends in principle on the centerof-mass momentum of the quarkonium state q. We will eventually neglect this dependence and write jnSi for jn; qi. The paper is structured in the following way. In Sec. II we briefly describe the evolution of the medium in a simple model. In Sec. III we derive the general evolution equations for the heavy-quark-antiquark densities in the framework of pNRQCD. This section contains the main theoretical results of the paper. In Sec. IV we write the Lindblad equation and solve it for a quarkonium in a weakly coupled plasma that fulfills the hierarchy 1=a 0 ≫ T ≫ E ≫ m D , where m D is the Debye mass (for details see Appendix B), while in Sec. V we do the same for a quarkonium in a strongly coupled plasma that fulfills the hierarchy 1=a 0 ≫ T ∼ m D ≫ E (the weakly coupled case 1=a 0 ≫ T ≫ m D ≫ E is also addressed as a particular case; see Appendix C). Numerical results are presented for the bottomonium states ϒð1SÞ and ϒð2SÞ. The analytical solution of the Lindblad equation in the case of static sources is derived in Appendix A. Finally, in Sec. VI we draw some conclusions and discuss possible developments. A concise version of the evolution equations and their solution in the case of a strongly coupled plasma has been presented in [10].

II. TIME EVOLUTION OF THE THERMAL MEDIUM
We will consider a medium that is infinite, homogeneous, and isotropic in space but that changes with time. These approximations are appropriate for very large nuclei in central collisions. For the description of these kinds of systems we can use Bjorken's evolution [11]. These assumptions about the medium could be relaxed by describing the medium using relativistic hydrodynamics (see [12] for a review); however, this is not the largest source of uncertainty in our calculation (higher-order corrections, uncertainties related to the hierarchy of scales, and the assumed Coulombic nature of some excited quarkonium states appear to have a larger impact).
According to [11], the effective temperature of the system evolves in time following where T 0 and t 0 are, respectively, the initial temperature and initial (proper) time, and v s is the velocity of sound in the medium. In a deconfined plasma at a very high temperature, v 2 s ¼ 1=3. As values of T 0 and t 0 for central collisions at the LHC, we use T 0 ¼ 475 MeV and t 0 ¼ 0.6 fm. These values are taken from [13].
We will study collisions with different centralities. As we are assuming the plasma to be homogeneous and isotropic, the only effect that a difference in centrality will produce will be to change the initial value of the energy density and hence T 0 . For this we assume the initial energy density to be proportional to T 4 0 (this is consistent with the sound velocity that we are using). Moreover, we assume the initial energy density to depend on the number of participants. From (2.6) of [14] we obtain where b is the impact parameter; A the number of nucleons; and σ the nucleon scattering cross section, which, according to [15], is taken as σ ¼ 64 AE 5 mb ¼ 164 AE 13 GeV −2 and where ρ A is the nucleon density distribution in the nucleus. Following [15], we approximate ρ A by a Woods-Saxon profile (assuming that the nucleon density is proportional to the charge density) taken from [16]: where which, in the case of lead at the LHC, is 207.
Experimentally, the results are given in terms of centrality bins. When possible [15], we characterize each centrality bin used in CMS by its mean impact parameter. There are some centrality bins whose mean impact parameter is not given in [15]. In these cases we perform the computation ourselves, following [17]. In all calculations in the rest of the paper, we will describe the different centralities by the different T 0 given in Table I. For all the centrality bins given in this table except 20%-90%, the value of T 0 corresponds to T 0 ðhbiÞ [which is similar to hT 0 ðbÞi]. In the case of 20%-90% centrality, the value given is hT 0 ðbÞi due to the fact that the bin is so large that it is not well represented by the average impact parameter. In Fig. 1 we show the time evolution of the temperature according to (18) for the most and least central collisions of Table I. The evolution starts at t 0 ¼ 0.6 fm and ends at about 4 fm for the most central collisions and at about 1.1 fm for the most peripheral ones, when the fireball has cooled down to a temperature of about 250 MeV, which is the smallest temperature, but still larger than the crossover temperature to the quark-gluon plasma, where we expect the framework of our calculation to be safely realized. Indeed, we will need T ≫ E to be always fulfilled, and, in the weakly coupled plasma case, we will also need T large enough to trust our perturbative calculations. In the strong coupling case, the last limitation does not apply, but in practice the lack of lattice data below 1.5T c for the relevant nonperturbative quantities forces us to stick to the abovementioned temperature. We recall that the crossover temperature to the quark-gluon plasma, T c , is expected to be at about 150 MeV [18][19][20]. Henceforth, we will compute the quarkonium nuclear modification factors R AA at the earlier time when a temperature of about 250 MeV is reached. Clearly the last piece of the evolution from 250 MeV to the freeze-out temperature is missing. Nevertheless, we expect it to have a modest impact on our results.

III. TIME EVOLUTION OF THE QUARKONIUM DENSITY MATRIX FOR 1=a 0 ≫ T
The typical LHC temperature of the fireball goes from around 475 MeV to the freeze-out temperature. The energy scale induced by the temperature is πT, which is at most about 1.5 GeV. Since the fireball expands very fast at the beginning of the thermal evolution, it will reach very soon lower energies. The inverse of the Bohr radius, 1=a 0 , of the ϒð1SÞ is in the vacuum about 1.3 GeV. The conclusion is that, for most of the evolution of the fireball (except perhaps for a very short time at the beginning), the condition (10) is fulfilled at least by the bottomonium ground state. Moreover, 1=a 0 ≈ 1.3 GeV ≫ Λ QCD implies that the bound state is Coulombic, i.e., described by a Coulomb potential.

A. pNRQCD for an open quantum system
Under the condition (10) and the assumption that the bound state is Coulombic, 1=a 0 ≫ Λ QCD , we can describe the quarkonium evolution in the fireball using pNRQCD. The Lagrangian of pNRQCD at next-to-leading order in the multipole expansion is [5][6][7] where r is the distance between the heavy quark and the heavy antiquark, and S ¼ S1 c = ffiffiffiffiffiffi N c p and O ¼ ffiffi ffi 2 p O a T a stand for the heavy-quark-antiquark fields in a colorsinglet and color-octet configuration, respectively. The operator h s ¼ p 2 =M þ V s is the color-singlet Hamiltonian. The potential V s is the color-singlet potential, which at leading order reads is the Casimir of the fundamental representation and N c ¼ 3 is the number of colors. The operator h o ¼ p 2 =M þ V o is the color-octet Hamiltonian. The potential V o is the color-octet potential, which at leading order reads V o ¼ α s ð1=a 0 Þ=ð2N c rÞ. We have made manifest that the strong coupling in the potentials is evaluated at a scale that is of the order of the inverse Bohr radius. We have set equal to 1 the Wilson coefficients of the dipole operators (corrections are suppressed by powers of α s and are beyond our aimed leading-order accuracy). The term L light is the QCD Lagrangian with light quarks.
In (22) there is a covariant derivative acting on the octet field O. This can be eliminated by means of suitable field redefinitions: (22) and to rename the fields O and E into O 0 and E 0 . The field Ω can be chosen to be a Wilson line going from −∞ to t: In the following, we will adopt these field redefinitions and we will understand the octet field and the chromoelectric field as the redefined ones. We drop, however, the superscript 0 to simplify the notation.
As discussed in Sec. I, we want to compute Trfρðt 0 ÞS l † ðt; r; 0ÞS l ðt; r 0 ; 0Þg. The experimental fact that the number of bottom quarks found in heavy-ion collisions is much smaller than that of lighter quarks implies that Trfρðt 0 Þg is dominated by contributions to the density matrix, ρ, coming from states made of light quarks and gluons. This is clearly the case in thermal equilibrium, where contributions coming from the heavy quarks are suppressed by a factor e −M=T . A consequence of this is that Trfρðt 0 ÞS l † ðt; r; 0ÞS l ðt 0 ; r 0 ; 0Þg ≪ Trfρðt 0 ÞS l ðt 0 ; r 0 ; 0ÞS l † ðt; r; 0Þg and similarly for the octet field.
We may look at the heavy quarks as an open quantum system that interacts with the (slowly) evolving medium of the fireball made of light quarks and gluons. The computation can be done in the close-time-path formalism (see [21] for the thermal equilibrium version and [22] for the nonequilibrium case). The formalism consists in rewriting field correlators by allowing the time to evolve from a path that goes from −∞ to ∞ and then from ∞ − iϵ to −∞ − iϵ. The fields are then ordered along the path. To make the ordering manifest, we call the fields on the upper branch type 1 and those on the lower branch type 2, and we identify them by a corresponding index. Hence, we write Trfρðt 0 ÞS † ðt; r; RÞSðt 0 ; r 0 ; R 0 Þg ¼ hS 1 ðt 0 ; r 0 ; R 0 ÞS † 2 ðt; r; RÞi ≡ hr 0 ; R 0 jρ s ðt 0 ; tÞjr; Ri; ð23Þ Trfρðt 0 ÞO a † ðt; r; RÞO b ðt 0 ; r 0 ; R 0 Þg where in each equation the last equality defines the colorsinglet and the color-octet density respectively. Since we do not have a preferred direction in color space, the octet density is taken as diagonal in color. We further assume that the heavy quarks are comoving with the medium. Under this simplifying assumption, the evolution does not depend on R and we drop it from the arguments of the densities. To simplify the notation, we have also dropped the dependence on the heavy-quarkantiquark distance; nevertheless, the densities should be understood as operators that depend on it.
We emphasize that the situation that we consider here is different from the one studied in [8,9], where the system formed by the plasma plus the heavy-quark states was assumed to be in thermal equilibrium. Under that condition the density of heavy-quark states is exponentially suppressed and, when computing heavy-quark correlators, we only need to include heavy-quark fields living on the upper branch of the closed-time path. This is not the case here, where the number of heavy-quark states is not the equilibrium distribution. As a technical remark, we further observe that the 12 propagator does not select a specific time ordering, while the 11 and 22 propagators describing heavy-quark fields living on the upper and lower branch, respectively, select instead the forward and backward propagation: they are proportional to θðt 0 − tÞ and θðt − t 0 Þ, respectively.

B. Evolution equations
We assume the following simplified model for the evolution of the quarkonium in the medium.
(a) From t ¼ 0 (which could also be taken as a time in the infinite past, t ¼ −∞) until t ¼ t 0 , the heavy quarks evolve as in the vacuum. Therefore up to corrections of relative order a 2 0 E 3 , which are negligible with respect to thermal corrections as long as T ≫ E, the color-singlet and color-octet densities evolve as The singlet and octet Hamiltonians have to be understood as operators in the relative-distance space like the densities. The initial conditions, ρ s;o ð0; 0Þ, of the densities will be discussed in suddenly the medium appears and the heavy quarks start interacting with it. We model the medium as a fireball that follows Bjorken's time evolution (see Sec. II). Since the number of heavy quarks in the medium is relatively small, we organize the computation as an expansion in the heavyquark densities, ρ s and ρ o . The 12 propagators are proportional to these densities. We compute them by keeping only Feynman diagrams with a single 12 propagator of heavy quarks, which amounts to considering only terms that are linear in the heavy-quark densities. Diagrams that contain two 12 correlators are quadratic in the densities of heavy quarks, and so on. Since we compute only diagrams that are linear in the heavy-quark densities, we consistently ignore the density dependence in the 11 and 22 propagators too. The initial conditions ρ s;o ðt 0 ; t 0 Þ are determined from the evolution starting at t ¼ 0 and ending at t ¼ t 0 computed in (a). We calculate now corrections to ρ s ðt; tÞ for t > t 0 at order r 2 . The relevant diagrams are shown in Fig. 2. From the pNRQCD Lagrangian it follows that at zeroth order in the multipole expansion, ρ s ðt; tÞ is just given by the tree-level diagram The expression hE a;i ðt; 0ÞE a;j ðt 2 ; 0Þi, like similar expressions below, stands for the in-medium light degrees of freedom average of the correlator of two chromoelectric fields. We recall that the fields in the previous expression are the redefined ones and hence gauge invariant by themselves. In this work, we assume that the medium is isotropic in space and color, and locally and instantaneously in thermal equilibrium at a temperature T. The third diagram in Fig. 2 is the complex conjugate of this one: ð30Þ Finally the last diagram in Fig. 2 reads where H.c. stands for Hermitian conjugate.
Taking into account all the diagrams in Fig. 2 and deriving both sides with respect to t, we get the evolution equation for ρ s : A similar computation leads to the evolution equation for ρ o : FIG. 2. Diagrams contributing at order r 2 to ρ s . A single line stands for a singlet propagator, a double line for an octet propagator, and a curly line for gluons. The vertices (circle with a cross) are the chromoelectric dipole vertices of the pNRQCD Lagrangian (22). The numbers 1 or 2 near the vertices mean insertions of fields from the upper or lower branches of the closed-time path, respectively. In the second diagram we also write explicitly the time variables of the propagators according to Eq. (28).
Both Eqs. (33) and (34) (33) and (34) induces new terms at order r 4 only; hence this substitution is consistent with the accuracy of those equations, which is of order r 2 . Furthermore, promoting e −ih s;o ðt−t 0 Þ ρ s;o ðt 0 ; t 0 Þe ih s;o ðt−t 0 Þ to ρ s;o ðt; tÞ in the right-hand sides of Eqs. (33) and (34) makes the equations Markovian and, in particular, leads to the standard evolution equations for the density matrix in the case of unitary evolution . Finally, it has the advantage of providing a set of equations that hold for arbitrarily large times. One should notice in fact that the validity of the expansion in Fig. 2  ðt; tÞ sensitive to two different time scales. One will come from the energy exponentials e AEih s;o ðt−t 0 Þ and will be of the order of the inverse of the binding energy, 1=E, whereas the other will come from ρ I s;o ðt; tÞ and will be of the order of 1=ðα s a 2 0 Λ 3 Þ, due to the pNRQCD power counting. Λ is generically the next relevant scale in the system (T or E or a combination of them); hence The function Ξ so accounts for the production of singlets through the decay of octets, while the functions Ξ os and Ξ oo account for the production of octets through the decays of singlets and octets, respectively. There are two octet production mechanisms that can eventually be traced back to the two octet chromoelectric dipole vertices in the pNRQCD Lagrangian (22). Finally, we note that the conservation of the trace of the sum of the densities, which amounts to the conservation of the number of heavy quarks, requires the functions Σ s , Σ o , Ξ so , Ξ os , and Ξ oo to satisfy The above equations, relating the singlet and octet decay widths to the corresponding production matrix elements, represent nothing else than the optical theorem in the problem at hand. They are fulfilled by the expressions in (29), (32), and (35)- (37), which becomes apparent if they are written like in Appendix D.
In the following we will assume t − t 0 to be larger than any other time scale that appears in the problem. This is indeed so for any time of the order of magnitude of the freeze-out time. It also amounts to assuming that the time during which the subsystem is observed is much larger than the time scale of any correlation between the subsystem and the environment. Under this assumption, we will approximate Furthermore, we will approximate This is an exact equality in vacuum and in thermal equilibrium. In the out-of-equilibrium case that we are considering here, the equality is broken by the time dependence of the temperature. Nevertheless, if the evolution of the temperature is quasistatic, which is our case at large times [cf. (18) 1=T × dT=dt ∼ 1=t ≪ E], the time dependence of the temperature may be neglected at leading order. One should point out that this approximation may be problematic at early times.

C. Lindblad equation
Equations (38) and (39) form a system of equations that can be solved if the properties of the medium and the initial conditions for ρ s and ρ o are known. This set of equations is the main result of this work. Although it is possible to solve these equations numerically for given initial conditions, it is indeed very challenging and computationally expensive. For this reason we have chosen here to focus on cases where these equations can be simplified to a Lindblad equation [23,24]. The Lindblad equation is well known in the fields of quantum optics and quantum information. It was studied in relation with quarkonium in [25]. From a mathematical point of view, the Lindblad equation follows from requiring the time evolution of the density matrix of the open quantum system to be Markovian, to preserve the trace, the equation to be linear in the density, and the corresponding linear operator to be a completely positive map. It has the following form: where H is a Hermitian operator and the operators C n are called collapse operators. In our case, ρ is the matrix

D. Expansion in spherical harmonics
We solve the Lindblad equation using numerical libraries available in the literature [26,27] and putting the system on a lattice. However, this lattice is three dimensional, making the number of entries for operators in the Lindblad equation still prohibitively large for our present means. As a way to deal with this practical difficulty, we do an expansion in terms of spherical harmonics. We define Since there is no preferred direction in space, during the entire evolution only the components with l ¼ l 0 and m ¼ m 0 are nonzero. Moreover, by the same argument, all polarizations are equally possible; therefore, all the information can be encoded in The normalization of ρ l has been chosen to ensure that P l ρ l ¼ Trfρg. The crucial approximation that simplifies the numerics of the calculation is to consider only l ¼ 0 and l ¼ 1. We have checked that the results change very little if we also include l ¼ 2. The reason is that we are interested in the suppression factor of S-wave states. Via a chromoelectric dipole transition, these states can decay to or be generated by P-wave states only. Hence, they are affected by states with l larger than 1 only indirectly. Under this approximation the density matrix can be written as The corresponding Hamiltonian is E. Initial conditions for ρ s and ρ o The computation of the production cross section of quarkonium in pp and pA and its extrapolation to AA collisions is a nontrivial problem and an active topic of research (see [28] and references therein). An additional difficulty is that results for production are usually written as expectation values for the different quarkonium states, while what we need is the reduced density matrix, which contains more information (for example, the relative phase between the different states). Moreover, the thermalization procedure from the collision time until the hydrodynamic regime and the way in which quarkonium is affected by the medium between these two times are still largely unknown.
Facing these difficulties, we choose to make a naive assumption about the form of the initial conditions. We impose that the relation between ρ s and ρ o at the initial time be controlled by just one parameter δ. The fact that the creation of heavy quarks requires high energies tells us that singlets and octets will be formed in a configuration similar to a Dirac delta, which implies an S-wave state; this assumption can also be found in [29]. Because our evolution equations are linear in the densities and we are interested in R AA , which is a ratio, we do not need to care about the absolute size of ρ s and ρ o but only about their relative size. The production cross section of singlets in S-wave states was computed at leading order in [30] and that of octets in [31]. Since the production of singlets is α s suppressed compared to that of octets, we use as an initial condition that at collision time (t ¼ 0) where j0i is an eigenstate of r with eigenvalue r ¼ 0. The normalization N is fixed by Trfρ s g þ Trfρ o g ¼ 1.
In the following we will perform the computations using different values of δ (1, 0.1, and 10). Between the collision of the two ions at t ¼ 0 and the beginning of the hydrodynamic evolution at t ¼ 0.6 fm, we assume that quarkonium evolves as if it were in vacuum. We use the bottom mass M ¼ 4.8 GeV. We obtain the Bohr radius for the 1S bottomonium state from the condition 1=a 0 ¼ MC F α s ð1=a 0 Þ=2. Using the MS scheme and one-loop running with Λ MS ¼ 250 MeV, we get 1=a 0 ¼ 1.334 GeV. The pNRQCD singlet and octet Hamiltonians are given after (22) and for each angular momentum component in (51). For the singlet and octet potentials V s and V o , we use the expressions given below (22) in accordance with the tree-level matching of pNRQCD with NRQCD. To compute the evolution of ρ s and ρ o (in vacuum and in the medium) we use a lattice of size 40a 0 and spacing a 0 =10. For the numerical computation we make use of the library Qutip [26,27].

IV. QUARKONIUM IN A WEAKLY COUPLED
PLASMA: 1=a 0 ≫ T ≫ E ≫ m D In this section, we derive the Lindblad equation and solve it for a weakly coupled quark-gluon plasma in a particular thermodynamical regime. A plasma is weakly coupled if T ≫ m D , where m D is the Debye mass. Because in perturbation theory m D ∼ gðTÞT, the condition is indeed fulfilled if the coupling gðTÞ is small. If there is no dynamical scale in between T and m D the evolution equations are a particular case of the strongly coupled case, T ∼ m D , that we will address in more generality in the next section. Explicit expressions can be found in Appendix C.
Here we consider the particular thermodynamical regime where the binding energy is in between T and m D . This means that we consider a system that fulfills the hierarchy of scales: 1=a 0 ≫ T ≫ E ≫ m D ; Λ QCD . In this situation and for thermal equilibrium, the modifications to the heavy quarkonium dynamics have been studied in detail [32]. Both the energy levels and the decay widths get thermal corrections. The thermal corrections to the width are mostly due to gluodissociation [33] (the dissociation of a quarkonium through the scattering with a real-time-like gluon from the medium), which at temperatures such that E ≫ m D is a dissociation mechanism more important than the dissociation by inelastic parton scattering [34] (the dissociation of the quarkonium through scattering with partons in the medium). In [35] it has been argued that this hierarchy of scales might be realized by the ϒð1SÞ state produced in heavy-ion collisions at the LHC.
Whereas it is clear that the real and imaginary parts of the singlet potential computed in [32] must be related to Σ s , it may not be obvious that they can be taken to be the same. In Appendix B we show that up to a redefinition of the density matrices, this is indeed the case. An additional motivation for this redefinition is to obtain a system that can be solved using a Lindblad equation. We show in Appendix D that, in general, (38) and (39) are not positive, and hence cannot be written in the Lindblad form. However, if we consider a medium that changes from thermal equilibrium in a quasistatic way, and an evolution time large enough so that Σ s only depends on the medium through the temperature (see the discussion at the end of Sec. III B), then (38) and (39) simplify and can be brought to the Lindblad form by the same redefinitions that led to the thermal equilibrium results of [32]. Taking over those results by just allowing a slow change of temperature with time, we see that they read at leading order where we have distinguished between the strong coupling running with the temperature scale, μ T ¼ πT; the one running with the scale of the energy, μ E ≳ −E 0 ¼ 1=ðMa 2 0 Þ; and the one coming from the Coulomb potential running with the inverse of the Bohr radius.
The thermal modifications to the octet potential were not calculated in [32]. The calculation is very similar to the one for the singlet. We present it in Appendix B 3 for completeness. The result reads at leading order The different Ξ's must be computed directly from (32), (36), and (37). They must also fulfill (42) and (43), which follow from the conservation of the number of heavy quarks. Those equations, for instance, implement the fact that the decay width of the singlet is related with Ξ os because any singlet that decays contributes to the creation of an octet. As discussed in Appendix B, the final expressions for the Ξ's can be cast, after a suitable redefinition of the density matrices, in the following form: which are suitable to be written in a Lindblad form. The Hermitian operator H entering the Lindblad equation is Furthermore, the Lindblad equation is made in this case by a set of nine collapse operators C 0 i , C 1 i , and C 2 i , which are To solve the Lindblad equation numerically, we proceed along the lines discussed in Sec. III D and perform an expansion in spherical harmonics, keeping only the S-and P-wave terms. The Hamiltonian reduces to The nine collapse operators above combine into three upon projection, The operators ∂=∂r and ∂=∂r þ 2=r come from p i acting on S-wave and P-wave states, respectively.

A. Results
In Fig. 4 (left plot), we show the results that we obtain using μ E ¼ 900 MeV (chosen to be the closest scale to −E 0 where perturbation theory may still apply) and the initial conditions as defined in Sec. III E with δ ¼ 1. We see that the suppression is slightly larger for ϒð1SÞ than for ϒð2SÞ and strong for both. This is in fact consistent with the leading order (LO) behavior of the decay widths calculated in [32], but clashes with experimental observations. This can be due to the fact that one (or several) of the assumed hierarchies of energy scales may not be fulfilled (e.g., if T ∼ m D , the plasma is not weakly coupled; see the next section) or that α s is not small enough. In particular, for our choice of μ E , α s ðμ E Þ is rather large. Hence, it could happen that perturbation theory is still reliable at the scale T, but not at the scale E. If we repeat the computation using μ E ¼ πT, we obtain the result shown in Fig. 4 (right plot). We see now that the suppressions of ϒð1SÞ and ϒð2SÞ are similar, but still very strong, both features are difficult to reconcile with experimental observations. Finally, we investigate the sensitivity to the initial conditions. We show in the left and right plots of Fig. 5 the results for δ ¼ 0.1 (more singlets than in the LO NRQCD production) and δ ¼ 10 (more octets than in the LO NRQCD production), respectively. While the case δ ¼ 0.1 is similar to the case δ ¼ 1, for δ ¼ 10 we observe slightly less quarkonium suppression. This indicates that above some threshold the suppression pattern will be quite sensitive to the ratio between the singlets and the octets initially produced. In particular, the larger the initial fraction of the quarkantiquark color octets is, the more marked is the feedback of singlets coming from the octet decays. This is at the origin of the small kink at about 2 fm in the right plot of Fig. 5.
In this section we apply the general equations derived in Sec. III to the case in which the thermodynamical scales are smaller than 1=a 0 but larger than the binding energy E. Thermodynamical scales are T and the Debye mass m D . In a weakly coupled plasma, one assumes that T ≫ m D ∼ gT. In this section, however, we assume more generally that the plasma is strongly coupled. This means that we take T ∼ m D . As we will see, the evolution equations can be written in terms of just two real constants [10]. Because the entire information about the medium is contained in these two constants that we are not going to evaluate, we do not need to make here any special assumption about the properties and degrees of freedom of the medium (which in the previous section was explicitly taken to be a quarkgluon plasma). Hence, everything that we write in this section applies to a generic strongly coupled hot medium.
If 1=a 0 ≫ T, Λ QCD , we can describe the evolution of the heavy-quark densities in the fireball by means of the equations found in Sec. III. If T is larger than E, we may neglect the energy-dependent exponentials e AEih s;o ðt−t 2 Þ that appear in the definitions of the Σ's and Ξ's: for t − t 2 ∼ 1=T, e AEih s;o ðt−t 2 Þ ∼ 1 holds. Using also the approximations (44) and (45) (this is where the quasistatic approximation enters), we can write Ξ os ðρ s ; tÞ ¼ r i ρ s r i κðtÞ; ð72Þ The real quantity κ is the heavy-quark momentum diffusion coefficient [36,37]: where T stands for time ordering. The quantity κ is related to the thermal decay width of the heavy quarkonium. In particular for 1S states, we have [see (41)] The heavy-quark momentum diffusion coefficient has been recently computed on the lattice [38]: The above estimate has been obtained from pure SU (3) calculations at temperatures of about 1.5 T c . Higher temperatures seem to suggest a smaller κ and/or a temperature dependence of κ=T 3 [39]. Also, light quarks may modify the above values. A perturbative determination of κ at next-to-leading order can be found in [40]. At next-to-leading order, κ=T 3 turns out to be in the range 9 ≲ κ=T 3 ≲ 18, but, as pointed out in [40], the perturbative series does not converge for any realistic value of the strong coupling. In fact, at leading order [see (C4)] it turns out to deliver a negative κ for realistic LHC temperatures. At the moment no definite statement can be drawn from looking at κ in perturbation theory. The real quantity γ is defined as The quantity γ is related to the thermal mass shift of the heavy quarkonium. In particular for 1S states, we have [see (40)] So far γ has not been computed on the lattice. The only estimate we have for γ is the perturbative calculation done at leading order in [8]; see (C2). Note that in our setup both κ and γ depend on time through the evolving temperature of the plasma.
With the above functions, we can write the evolution equations (38) and (39) in the Lindblad form (46). In addition to the Hermitian Lindblad operator H, we need six collapse operators C 0 i and C 1 i , which are Following Sec. III D, the numerical computation is done by expanding the density matrix and the collapse operators in spherical harmonics, and keeping only S-and P-waves. The projected Hamiltonian reads and the six collapse operators above are combined into two,

A. Results
In the bottomonium case, the time evolutions of R AA for 30%-50% centrality and 50%-100% centrality are shown in the left and right plots of Fig. 6, respectively. Note that, in the left plot, the R AA for the 2S state becomes insensitive to κ at large times, an indication that it reaches a steady state before the quark-gluon plasma vanishes.
We have taken κ=T 3 in the range (76), while we have set γ ¼ 0 and δ ¼ 1. We have no a priori knowledge for either γ, since a nonperturbative determination of this parameter is missing, or δ, since we ignore the precise initial conditions. Hence, we have scanned for several values of γ and δ. Here, we have restricted ourselves to γ < 0, the sign of the perturbative result (C2). We find that the CMS results of [1] prefer small values of γ, which is the rationale for choosing γ ¼ 0. We further find that the results are rather insensitive to δ, in contrast to the weakly coupled case discussed in the previous section. The choice δ ¼ 1 assumes the initial ratio of octets over singlets to be just 1=α s ðMÞ.
In Table II we show our predictions for the centrality bins studied by CMS at 2.76 TeV in [41]. Results in this table are corrected for feed-down effects using the method of [42] with the updated feed-down fractions from [43]. The reason why feed-down is taken into account in the table and not in the time evolution plots is that it takes place after freeze-out. In order to be on the safe side regarding the condition 1=a 0 ≫ T ∼ m D ≫ E, we focus only on centralities between 30% and 70%. All our determinations are summarized and compared with the CMS data in Fig. 7. The theoretical error band accounts only for the lattice uncertainty in κ.
ϒ suppression in heavy-ion collisions has also been studied by the Alice Collaboration [44]. They have only considered the centrality bins 0%-20% and 20%-90%. The initial temperature for the 0%-20% centrality bin is too high for our present study. Regarding the centrality bin 20%-90%, the average initial temperature happens to be very similar to the one in the centrality bin 50%-70% and, therefore, our prediction is approximately the same. Analyses for LHC data at 5.02 TeV are under way [45].
In order to check that the hierarchy of scales assumed at the beginning of Sec. V is maintained during the evolution of the ϒ states in the fireball, we have computed the time evolution of the matrix element hnSjf1=r; ρ s g=2jnSi and of the binding energy E nS for bottomonium n ¼ 1 and n ¼ 2 states. This is shown in Fig. 8. We see that both h1Sjf1=r; ρ s g=2j1Si=h1Sjρ s j1Si and h2Sjf1=r; ρ s g= 2j2Si=h2Sjρ s j2Si remain close to their initial, in-vacuum, values: 1.334 and 0.884 GeV, respectively. For our choice of parameter γ ¼ 0, the Hamiltonian H coincides with the in-vacuum one [see (79)], and so do the binding energies: E 1S ≈ 0.37 GeV and E 2S ≈ 0.04 GeV. This implies that the in-vacuum hierarchy between the inverse sizes of the bound states and their binding energies is preserved during the evolution in the plasma. Furthermore, we also plot πT as a function of time for the 30%-50% centrality class. We observe that the hierarchy between the bound-state scales and the thermal scale is preserved for late time evolution.
Many effects that have not been considered in the present analysis or considered in a simplified form FIG. 6. Time evolution of R AA for bottomonium with κ=T 3 in the range (76), γ ¼ 0 and δ ¼ 1, for 30%-50% centrality (left plot) and for 50%-100% centrality (right plot). (e.g., the hydrodynamical evolution) have the potential to quantitatively impact the nuclear modification factor calculated here. Inside the framework presented here, the results depend on the initial conditions and on just the parameters κ and γ. The impact of a value of κ outside the range in (76) or of a positive value of γ is shown for two illustrative examples in Fig. 9. The result suggests that there may exist values of κ, γ, and δ that reproduce the present data. The fact that a value of κ lower than the lattice prediction is needed may be explained if most of the quarkantiquark pairs are moving with respect to the plasma [46,47], at least in the weak-coupling case. As mentioned before, the lattice results of [39] also seem to point to a lower value of kappa at higher temperatures. Note finally that a positive value of γ means that the medium is very different from a weakly coupled quark-gluon plasma, since the latter has a negative gamma. This reinforces the need for a lattice evaluation of γ.
From the comparison between strongly and weakly coupled results, it can be seen that there is more suppression in the weakly coupled case than in the strongly coupled one, which might seem surprising. In order to understand this, one has to take into account the several differences between the two cases. In the weakly coupled case, we restrict ourselves to very central collisions (in order to guarantee high temperatures and make the weakly coupled case plausible), while in the strongly coupled scenario, which is more phenomenologically oriented, we focus on less central collisions since there our approximations are more likely to be fulfilled. Another important difference is that the leading-order thermal corrections in the weakly coupled case are approximately linear with the temperature, while in the strongly coupled one they are cubic. In a system that expands following Bjorken evolution and with a sound velocity close to the ideal gas case, linear corrections tend to have a larger impact than cubic FIG. 8. Time evolution of hnSjf1=r; ρ s g=2jnSi (continuous line), the binding energy E nS (dotted line), and the temperature (dashed line) for bottomonium n ¼ 1 (left plot) and n ¼ 2 states (right plot), with κ=T 3 ¼ 2.6, γ ¼ 0, δ ¼ 1, and for 30%-50% centrality collisions.
FIG. 9. R AA as obtained using κ=T 3 ¼ 0.25 and γ ¼ 0 in the left plot and κ=T 3 ¼ 2.6 and γ=T 3 ¼ 6 in the right plot (dots) compared with the CMS data of [41] (triangles). Upper (red) entries refer to the ϒð1SÞ, and lower (green) entries to the ϒð2SÞ. The vertical dashed lines highlight the window in which we expect the approximation 1=a 0 ≫ T ∼ m D ≫ E to be valid.
ones. This can be seen by looking at our analysis of the static limit in Appendix A.

VI. CONCLUSIONS AND OUTLOOK
In the paper we present a systematic description in an effective field theory framework (pNRQCD) of heavyquark-antiquark systems as open quantum systems interacting with an environment made of light quarks and gluons, the fireball formed in heavy-ion collisions. While this work derives its original inspiration from [25], it is quite different from all previous studies on the subject [48][49][50][51][52][53][54][55][56], as the derived evolution equations fulfill three essential conditions: they conserve the total number of heavy quarks (i.e., Trfρ s g þ Trfρ o g is preserved by the evolution equations); they account for the non-Abelian nature of QCD (through gluon exchanges' color-singlet quarkonia may dissociate into quark-antiquark color-octet states, and vice versa quark-antiquark color-octet states may generate quarkonia); and, finally, they do not rely on classical approximations but rather follow from the closedtime-path formalism applied to quantum field theory. The work substantially extends, updates, and completes a previous strongly coupled analysis done in [10]. A recent study also based on pNRQCD can be found in [57].
The evolution equations in terms of quark-antiquark color-singlet and color-octet density matrices have been written in (38) and (39). The equations rely on the assumption that the typical inverse size of the quarkantiquark system, 1=a 0 , is larger than any other scale of the medium and larger than Λ QCD . This implies that the quark-antiquark interaction is mainly Coulombic and that the interaction with the medium may be multipole expanded. The evolution equations follow from the calculation of the singlet and octet density matrices at the leading nontrivial order in the multipole expansion. Since the heavy quark density is expected to be small, we only keep linear terms in them. Then we show that, at the order in the multipole expansion we are working, the time derivative of the density matrices at a given time can be written as a linear function of the density matrices at the same time. This produces an evolution equation that is reliable at arbitrarily large times and that turns out to be equivalent to the Schwinger-Dyson equations depicted in Fig. 3. The evolution equations (38) and (39) do not make any special assumption on the medium and may be valid either for a quark-gluon plasma or a different medium formed in the heavy-ion collisions. They are also valid either if the medium is in thermal equilibrium or if it is far from it (provided that no dynamical scale is larger than 1=a 0 , as mentioned above). The evolution equations preserve the number of heavy quarks, but, in general, they are not of the Lindblad form.
In order to get a Lindblad equation, we consider some specific cases where we assume that the time scales we are interested in are larger than any other time scale in the problem and that the evolution is quasistatic. For a quasistatic evolution, the environment may be locally in thermal equilibrium. Thermal equilibrium allows us to define at each time a temperature. In Sec. IV we consider the explicit case of a weakly coupled quark-gluon plasma. The temperature, T, and the Debye mass, m D , are assumed to be such that T ≫ E ≫ m D , where E is the typical binding energy of the heavy-quark-antiquark system. In Sec. V we consider a strongly coupled plasma at a temperature T ∼ m D ≫ E. In Appendix C we also discuss the case T ≫ m D ≫ E. In order to numerically solve the Lindblad equation, we further make a truncation on the number of partial waves taken into account. This approximation is of a technical nature, and it is only useful to simplify the solution in the face of limited computing time.
In the future, one may relax some of these approximations, or even solve directly the evolution equations (38) and (39).
The numerical solutions of the Lindblad equations are presented and discussed in Sec. IVA for the weakly coupled quark-gluon plasma, and in Sec. VA for the strongly coupled medium. We analyze some specific set of initial conditions, parameter, and time evolution of the thermal medium [we work with the Bjorken evolution (18)]. A more extensive study is surely due in the future. For instance, the Bjorken evolution could be substituted by more refined hydrodynamical models (see [58] for a review of the state of the art). There is also plenty of room for improvement at early times, where color glass condensate physics [59,60] may be included or higher-order NRQCD production results (see [61] and references therein) for the initial conditions of the density matrices may be incorporated. Eventually, the momentum dependence of the quarkantiquark pairs should also be addressed. This requires enlarging the Hilbert space of the density matrices and incorporating the effects of the relative motion of the pair with respect to the medium; see, for instance, [46,47].
The strongly coupled plasma case may be of particular interest at the LHC, since the temperature T or even πT may not be much larger than m D . The Lindblad equation for the strongly coupled case depends on only two parameters: the heavy-quark momentum diffusion coefficient κ, defined in (74), and γ, defined in (77). While κ has been computed on lattice QCD (although in pure gluodynamics and for a limited range of temperatures), γ has not. A first determination of γ remains therefore the most urgent missing ingredient for the computation of the suppression factor R AA . ACKNOWLEDGMENTS J. S. thanks Anton Andronic for providing invaluable information on the current codes used by the experimental collaborations to map the number of participants to average impact parameters. The work of N. B. and A. V. was supported by the Deutsche Forschungsgemeinschaft (DFG) Grant No. BR 4058/1-2 "Effective field theories for heavy probes of hot plasma" and by the DFG cluster of excellence "Origin and structure of the universe" (www .universe-cluster.de). The work of M. A. E. was supported by the Academy of Finland, Project No. 303756. J. S. was supported by the CPAN CSD2007-00042 Consolider-Ingenio 2010 program; the FPA2013-43425-P, FPA2013-4657, FPA2016-81114-P and FPA2016-76005-C2-1-P projects (Spain); and the 2014-SGR-104 grant (Catalonia).

APPENDIX A: STATIC LIMIT
In this appendix, we consider the evolution equations for static quarks under the condition Interestingly, this case can be solved analytically. Because T is much larger than the typical energy scale, the energy-dependent exponentials in (29), (32), (35), (36), and (37) can be set equal to 1. We obtain The evolution equations depend, therefore, on only one time-dependent parameter Γ s ðtÞ ¼ Σ s ðtÞ þ Σ † s ðtÞ, the color-singlet width, and they read The initial condition describes two heavy quarks at a given distance r in an arbitrary color state: The problem consists in solving the evolution equations (A4) and (A5) for the two functions ρ sr and ρ or . The solution is with uðtÞ ¼ N 2 Γ s ðtÞ . The static limit does not give us information on how the singlets are distributed in the different possible states [for example, ϒð1SÞ, ϒð2SÞ, and so on], but it does give qualitative information on how the population of singlets compares to that of octets. The crucial parameter is uðtÞ: for t − t 0 ≪ uðtÞ, the thermal medium has a small impact on the distribution of quarkonia, while for t − t 0 ≫ uðtÞ, the density approaches the large-time asymptotic value.
(a) As a first special situation, we consider a strongly coupled plasma: From (69) it follows that in this case, Γ s ðtÞ ¼ Σ s ðtÞ þ Σ † s ðtÞ ¼ κðtÞr 2 . The heavy-quark momentum diffusion coefficient κ has been defined in (74). The equations do not depend on the coefficient γ defined in (77).
We estimate the order of magnitude of uðtÞ by taking r ¼ a 0 [the Bohr radius of the ϒð1SÞ], the average temperature T ¼ 317 MeV (that we define as the average between the temperature at t ¼ 0.6 fm and T ¼ 250 MeV), and κ ¼ 2.5T 3 . We obtain u ∼ 4 fm, which is about the time the fireball temperature is above T c for central collisions (see Fig. 1). In Fig. 10 we plot ρ s as a function of t for the two extreme initial conditions.
(b) The second situation that we consider is the one of a weakly coupled quark-gluon plasma: FIG. 10. Evolution in time of the color-singlet density for static quarks in a strongly coupled plasma. The blue continuous line shows the evolution when the initial state is made only of singlets, whereas the green dashed line shows the evolution when the initial state is made only of octets. Asymptotically both curves approach 1=N 2 From (55) [32] were used in [62] to fit lattice results using α s ≈ 0.4 independently of the renormalization scale. If we proceed like that, we get u ≈ 2.16 fm. The plot of ρ s as a function of t for the two extreme initial conditions corresponding to this last choice of the coupling is shown in Fig. 11.
(c) Finally, let us consider the more general case where Γ s has an unspecified power-law dependence on the temperature: This case allows to understand, at least qualitatively, the interplay between the time evolution of the decay width and Bjorken's expansion in the time evolution of the colorsinglet and color-octet densities. Note that, if we neglect the running of the strong coupling, both previously considered situations, (a) and (b), are of this type.
Using the time dependence of the temperature in Bjorken's expansion (18), one gets which implies Inserting (A14) into (A8) and (A9) provides the solution of the evolution equations: ρ sr ðt; tÞ ¼ 1 − ρ or ðt; tÞ The solution shows three very different behaviors depending on the value of nv 2 s . (1) If nv 2 s > 1, the color-singlet density never reaches the value 1=N 2 c , which is its thermal equilibrium value. Instead it approaches the value The physical interpretation is that for nv 2 s > 1, the decrease with time of the decay width in the fireball, described by Bjorken's expansion, is so fast that the static quark-antiquark densities do not have time to equilibrate.
(2) In the case nv 2 s < 1, we are exactly in the opposite situation. In the absence of freeze-out effects that could modify the evolution, the color-singlet and color-octet densities reach their thermal equilibrium values exponentially fast. This is the situation realized by static quarks and antiquarks in the weakly coupled plasma of case (b), for which n ¼ 1.
which implies that ρ sr ðt; tÞ ¼ 1 − ρ or ðt; tÞ Like before, the color-singlet and color-octet densities reach after some time their thermal equilibrium values. Differently from the previous case, however, the falloff is powerlike instead of exponential, which means that the asymptotic values are reached slower [in fact, much slower if Γ s ðT 0 Þt 0 ≪ 1]. The physical interpretation is that we are describing a situation in which the speed of the expansion of the fireball is competing with the decrease of the decay width with the temperature. The actual time required to reach the thermal equilibrium values depends crucially on FIG. 11. Like in Fig. 10, but in a weakly coupled quark-gluon plasma with the choice α s ≈ 0.4 for all strong couplings.
the initial condition. This situation is realized by static quarks and antiquarks in the strongly coupled plasma of case (a), for which n ¼ 3.
In this appendix we consider quarkonium in a plasma that realizes the regime 1=a 0 ≫ T ≫ E ≫ m D . We aim at justifying Eqs. (54)-(60). We will start with some general remarks and then apply them to the case at hand.

Density matrices' redefinitions
The evolution equations (38) and (39) allow for the following redefinitions of the density matrices ρ s and ρ o and of the functions Σ s , Σ o , Ξ so , Ξ os , and Ξ oo that preserve the evolution equations at the given accuracy of order r 2 in the multipole expansion: This implies that at order r 2 we may neglect redefinitions of the density matrices inside the functions Ξ so , Ξ os , and Ξ oo . Moreover, whenever dρ s =dt or dρ o =dt multiplies functions of order r 2 , we can replace them, through the leading-order evolution equations, with −i½h s ; ρ s and −i½h o ; ρ o , respectively. In order for these transformations to preserve the trace of the total heavy-quark density, i.e., Trfρ s g þ Trfρ o g, the operators O s , O o , O so , O os , and O oo must be related in such a way that the transformed color-singlet and color-octet density matrices, and the transformed functions Σ s , Σ o , Ξ so , Ξ os , and Ξ oo fulfill the conditions (42) and (43). This is guaranteed if the final evolution equation is of the Lindblad form (see Sec. III C).

Computation of Σ s and Ξ so
Equations (54) and (55) were originally derived in [32] by computing the singlet self-energy in momentum space for a generic incoming energy and by dropping terms that would vanish on physical states. Such terms contribute to the wave function normalization only. We may distinguish two momentum regions in the chromoelectric correlator appearing in the self-energy. The momentum region scaling like the temperature T contributes with where ΔV ¼ h o − h s ¼ N c α s =ð2rÞ. The momentum region scaling like the energy scale E contributes with The factor of T appearing in (B9) is a consequence of the Bose enhancement. The real part of −i times (B8) gives (54), and the imaginary part of −i times (B9) gives (55). The function Σ s that appears in the evolution equation (38) has been defined in (29). It also gets contributions from the scale T and the scale E. The contribution from the scale T, which is computed by expanding the exponents and keeping terms linear in h o and h s , reads Σ s j T ¼ Eq:ðB8Þ þ i π 9 C F α s T 2 ½h s ; r 2 : The contribution from the scale E, which is computed by expanding the Bose distribution function for large T, reads Σ s j E ¼ Eq:ðB9Þ þ 2 3 C F α s T h s ; 1 2 ½h s ;r 2 þr 2 ΔV

: ðB11Þ
We see that neither (B10) nor (B11) coincides with (B8) and (B9). The difference is significant: for instance, in (B10) it gives rise to an imaginary part of −iΣ s that is larger than the one listed in (55). However, assuming that the time dependence of the temperature is a subleading effect (see the comment at the end of Sec. III B), we may reabsorb the difference in a redefinition of the color-singlet density, ρ s , and Σ s according to (B3) with This shift in the color-singlet matrix element has been taken into account when computing the results shown in Sec. IV A. Its numerical impact turns out to be small.

APPENDIX C: QUARKONIUM IN THE REGIME 1=a 0 ≫ T ≫ m D ≫ E
We consider here quarkonium in a plasma that realizes the regime 1=a 0 ≫ T ≫ m D ≫ E. This regime was studied in [8] for the static case and in [63] for muonic hydrogen. It is straightforward to combine those results to obtain the relevant expansions that enter the Lindblad equation. Under the condition 1=a 0 ≫ T ≫ m D ≫ E, all thermal contributions can be encoded in modifications of the potential. In addition, most of the contributions can be regarded as a particular case of Sec. V, where γ and κ can be computed in perturbation theory. The expected leading contributions to γ and κ, which would be of Oðα s ðμ T ÞT 3 Þ, turn out to vanish. As a consequence, the leading nonvanishing contributions to these quantities are of Oðα 2 s ðμ T ÞT 3 Þ and Oðα s ðm D Þm 2 D TÞ. They may compete in size with terms of relative size of Oðα s ðμ T ÞET 2 Þ calculated in (54) that come from higher-order terms in the E=T expansion and are not included in γ and κ. Putting all the contributions together, we obtain Imð−iΣ s Þ ¼ − r 2 2 κ; ðC3Þ where μ T ¼ πT, m 2 D ¼ 4πα s ðμ T ÞT 2 ðN c þ n f =2Þ=3, and ΔV ¼ N c α s ð1=a 0 Þ=ð2rÞ; n f is the number of active massless flavors (n f ¼ 3 in the bottomonium case); and ζ is the Riemann zeta function. Note that the term r 2 γ=2 in (C1) is suppressed by a factor a 0 T with respect to the other terms; we keep it nevertheless to maintain a closer analogy to the discussion of Sec. V.
The thermal contributions to the octet potential have not been calculated before. However, they can be easily obtained from the expressions above by following the same steps as in Appendix B 3. We obtain The quantities Ξ so , Ξ os , and Ξ oo are obtained as particular cases of (71)-(73), in which γ and κ take the values of (C2) and (C4), respectively. In addition, the Lindblad operator H is obtained from (79) by making the following replacements, It is not obvious that this regime applies to LHC temperatures. In particular, the separation of scales between T and m D does not seem to be large enough to guarantee a positive κ and hence a positive decay width. We can nevertheless compute the quarkonium suppression just due to screening, which amounts to computing R AA while ignoring the imaginary part of the potential and all the Ξ's. In Fig. 12, this suppression is computed from (17) in the case of the 1S and 2S bottomonium states. We can see that, while screening alone is able to make ϒð2SÞ disappear almost completely, ϒð1SÞ suppression is at most 10%.
If h were a positive definite matrix, then it would always be possible to redefine the operators L n i in such a way that the evolution equation would be of the Lindblad form (46). Since, however, h is not a positive definite matrix, the Lindblad theorem [23] does not guarantee that Eqs. (38) and (39) may be brought into a Lindblad form. A special case is the strongly coupled case studied in Sec. V. There L 1 i ∝ L 0 i and L 3 i ∝ L 2 i , which allows us to set to zero, after a redefinition of the operators L n i , the eigenvectors of h associated to negative eigenvalues, eventually leading to an evolution equation of the Lindblad form.