<figure><img src="image_1.png" /></figure>

HAL
open science

# Quarkonium dynamics in the quantum Brownian regime with non-abelian quantum master equations

Stéphane Delorme, Roland Katz, Thierry Gousset, Pol-Bernard Gossiaux,
Jean-Paul Blaizot

► To cite this version:

Stéphane Delorme, Roland Katz, Thierry Gousset, Pol-Bernard Gossiaux, Jean-Paul Blaizot. Quarkonium dynamics in the quantum Brownian regime with non-abelian quantum master equations. Journal of High Energy Physics, 2024, 06, pp.060. 10.1007/JHEP06(2024)060. hal-04470951

HAL Id: hal-04470951

https://hal.science/hal-04470951v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

RECEIVED: February 28, 2024

REVISED: May 9, 2024

ACCEPTED: May 10, 2024

PUBLISHED: June 12, 2024

# Quarkonium dynamics in the quantum Brownian regime with non-abelian quantum master equations

Stéphane Delorme<sup>id,a,b,*</sup> Roland Katz<sup>id,b</sup> Thierry Gousset<sup>id,b</sup> Pol Bernard Gossiaux<sup>id,b,*</sup>
and Jean-Paul Blaizot<sup>id,c</sup>

E-mail: stephane.delorme@ifj.edu.pl, rolandkatz@outlook.com,
gousset@subatech.in2p3.fr, gossiaux@subatech.in2p3.fr,
jean-paul.blaizot@ipht.fr

**ABSTRACT:** We present numerical solutions in a one-dimensional setting of quantum master equations that have been recently derived. We focus on the dynamics of a single heavy quark-antiquark pair in a Quark-Gluon Plasma in thermal equilibrium, in the so-called quantum Brownian regime where the temperature of the plasma is large in comparison with the spacing between the energy levels of the $Q\bar{Q}$ system. The one-dimensional potential used in the calculations has been adjusted so as to produce numbers that are relevant for the phenomenology of the charmonium. The equations are solved using different initial states and medium configurations. Various temperature regimes are studied and the effects of screening and collisions thoroughly analyzed. Technical features of the equations are analyzed. The contributions of the different operators that control the evolution are discussed as a function of the temperature. Some phenomenological consequences are addressed.

**KEYWORDS:** Finite Temperature or Finite Density, Quarkonium, Effective Field Theories of QCD

ARXIV EPRINT: 2402.04488

<footnote><sup>a</sup>Institute of Nuclear Physics, Polish Academy of Sciences,
ul. Radzikowskiego 152, 31-342, Kraków, Poland</footnote>

<footnote><sup>b</sup>IMT Atlantique, Nantes Université, CNRS/IN2P3, SUBATECH,
44307, Nantes, France</footnote>

<footnote><sup>c</sup>Université Paris-Saclay, CNRS, CEA, Institut de physique théorique,
91191, Gif-sur-Yvette, France</footnote>

<footnote>*Corresponding author.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Equations for the $Q\bar{Q}$ density matrix</td><td>4</td></tr><tr><td>2.1</td><td>Master equations</td><td>4</td></tr><tr><td>2.2</td><td>Various regimes and characteristic temperatures</td><td>7</td></tr><tr><td>3</td><td>Main physics results</td><td>11</td></tr><tr><td>3.1</td><td>In-medium 1S singlet initial state</td><td>11</td></tr><tr><td>3.2</td><td>1S vacuum singlet initial state</td><td>16</td></tr><tr><td>3.3</td><td>Initial compact wave packet</td><td>18</td></tr><tr><td>3.4</td><td>Charmonium in an expanding plasma</td><td>20</td></tr><tr><td>4</td><td>Discussion of more technical issues</td><td>21</td></tr><tr><td>4.1</td><td>The role of various superoperators in the evolution</td><td>22</td></tr><tr><td>4.2</td><td>Asymptotic distributions and steady states</td><td>30</td></tr><tr><td>5</td><td>Conclusion</td><td>38</td></tr><tr><td>A</td><td>Positivity and structure of the QME</td><td>40</td></tr><tr><td>B</td><td>Implementation of the higher-order operator and UV regularization</td><td>42</td></tr><tr><td>C</td><td>Color transition operators</td><td>45</td></tr><tr><td>D</td><td>Details on the numerical implementation</td><td>46</td></tr><tr><td>E</td><td>Charmonium dissociation rates</td><td>49</td></tr><tr><td>F</td><td>Evolution under $\mathcal{L}_2$ alone</td><td>52</td></tr><tr><td>G</td><td>Equilibrium limit of the momentum distribution</td><td>61</td></tr></table>

## 1 Introduction

The quarkonia, bound states of heavy quark-antiquark pairs, have been used for a long time to diagnose the properties of matter produced in ultrarelativistic heavy ions collisions [1]. The dynamics of bound states are indeed very much affected by the presence of a medium in which they propagate: they can in particular disappear due to screening or collisional effects if the temperature of the medium is high enough; if the number of $Q\bar{Q}$ pairs produced in a given collision is large, one may anticipate that these pairs recombine into quarkonia during the evolution of the system. Thus, the measurements of the properties of quarkonia produced in ultrarelativistic heavy ion collisions can provide information on the state of the medium. And indeed many experimental data provide clear evidence that quarkonia are strongly affected by the matter produced in high energy nuclear collisions [2].
Extracting significant information from these data is however not a simple task. We are facing a rather complicated many-body problem, which requires going beyond the elementary approaches which have been used in the past. In recent years, numerous efforts have been made to adapt to this complicated problem tools that have been developed in the study of open quantum systems [3–16]. Such approaches generically lead to Quantum Master Equations (QME) that capture the most important effects of the medium on the dynamics of the $Q\bar{Q}$ heavy quark pairs (for recent reviews see [17–19]). Generically two regimes can be identified, depending on whether the bound states of the $Q\bar{Q}$ pairs are an important part of the dynamics, or not. Which regime occurs depends on the temperature $T$. The so-called “Quantum Optical” regime corresponds typically to “low” temperatures, $T \ll \Delta E$ (and $T \lesssim \Delta E$ in practice), where $\Delta E$ denotes the typical energy difference between two eigenstates of the (in-medium) Hamiltonian. In this regime, the effect of the medium can be seen as a perturbation applied to the bound states of the system, modifying their binding energies and contributing (small) widths to the states. This is a regime where rate equations or Boltzmann equations can be derived from the QME [20, 21] and where other approaches, such as transport approaches based on these rate equations can also prove useful [16, 22–27]. In the other regime, which takes place at “high” temperatures, $T \gg \Delta E$ (and $T \gtrsim \Delta E$ in practice), the binding forces represent a small perturbation of the dynamics of nearly independent heavy quarks subjected to the influence of stochastic forces generated by the collisions with the plasma constituents. This regime is commonly referred to as the Quantum Brownian Motion (QBM), and this is the regime we shall be concerned with in this paper.

These two regimes are usually described by two different types of master equations, and the way to analyze the corresponding physics refers usually to two different basis. In the optical regime, the basis of localized bound states is the most appropriate one. In the QBM regime however, the analysis is better carried out in coordinate space, which plays the role of a “preferred basis”. This is because the operators that represent the effects of the collisions of the heavy quarks with the plasma constituents are local operators, and are (nearly) diagonal in the coordinate space basis. In this basis, the collisions bring the density matrix of the $Q\bar{Q}$ system to a diagonal form in a relatively short time scale, which is not necessarily the case in other basis. This entails simplifications in the description of the dynamics, suggesting in particular the usefulness of semi-classical approximations [14, 28]. The dynamics of the color degrees of freedom remain fully quantum however, and they are treated as equations that couple the different color components (singlet and octet) of the $Q\bar{Q}$ density matrix.

Common to all the derivations of the QME in the present context is the fact that the mass of the heavy quarks is large, which implies that their motion is slow on the time scale of collective plasma dynamics. This is used to obtain Markovian equations, which in most approaches take the form of Lindblad equations. In the case of quantum chromodynamics (QCD), further simplifications have been implemented, relying on the use of effective field theories. In addition to the large quark mass, such approaches focus on the dipolar nature of the microscopic interactions, and exploit various scale hierarchies that allow for systematic and controlled approximations. This is the case of the non-relativistic QCD (NRQCD) approach [6, 8, 10, 14], and the so-called potential non-relativistic QCD (pNRQCD) approach [11, 12, 16, 29–31]. Various strategies have been followed to solve the resulting equations for one single $Q\bar{Q}$ pair: by transforming the quantum master equations into Schrödinger equations with a complex potential [32–34] or into stochastic Schrödinger
equations [35, 36], by expanding them to the first spherical harmonics [11, 12, 29–31], by performing a direct resolution in one dimension and limited to the quantum electrodynamics (QED) case [13, 37] or by reducing them to semi-classical Langevin equations [10, 14, 15]. The latter is of particular interest, as it allows to treat multiple quark-antiquark pairs simultaneously, which has not been done yet with a full quantum equation, except in the abelian case [10].

In this paper, we present the full numerical solution of the QME derived in [14], in the QBM regime and for the particular case of a single $Q\bar{Q}$ pair, using a potential tuned specifically for such one-dimensional studies. It is the first time that the complete solution of these particular equations is obtained. We note however that these equations have essentially the same physics content as those studied by different means in ref. [36]. A one-dimensional setting limits some aspects of the dynamics, but as we shall see, it provides interesting information on the evolution of a $Q\bar{Q}$ in a medium in thermal equilibrium. It allows us to identify what we consider as robust features that are likely to survive in a more realistic three-dimensional setting, while being much less computationally demanding.

The paper is organized as follows. Section 2 is an introductory section where we present a brief summary of the QME that have been derived in [14]. These equations describe the evolution of the reduced density matrix of heavy quark-antiquark pairs coupled to a quark gluon plasma in thermal equilibrium. These equations are written in the coordinate space representation, and take the form of a Lindblad equation, involving a set of four Liouvillian superoperators $L_i$. We discuss the various regimes that these equations can describe as a function of the temperature of the plasma, and we provide estimates of the characteristic temperatures that delineate these various regimes, based on the parameters of the specific one-dimensional potential that we use. The main body of the paper is divided into two sections. In the first one, section 3, we present the results obtained by solving the QME for a single charm-anticharm pair with different initial conditions for the pair, and various temperatures of the plasma. Assuming a medium with fixed temperature, we first study the evolution and late stage steady state of a pair initially in a color singlet eigenstate of the screened Hamiltonian. The evolution of the density matrices, the marginal distributions and the instantaneous projections on the vacuum eigenstates are analyzed in detail. We then investigate more systematically the effect of the medium temperature by either starting from a 1S-like vacuum eigenstate or from a more compact Gaussian singlet state. The impact of the initial state on the evolution toward equilibrium is further analyzed by testing an initial P-like color octet state. Most of the analysis of this section involves a medium of fixed temperature. At the end of the section we consider the modification of the dynamics that take place in an expanding plasma, assuming a simple one-dimensional Bjorken expansion [38]. The second major section of the paper, section 4, addresses several more “technical” issues related to generic properties of the QME. In section 4.1, we inspect the role of the various superoperators of the Lindblad equation at various stages of the evolution and for the various temperature regimes introduced in section 2. A particular discussion concerns the higher order superoperator $L_4$ in ensuring the positivity of the density matrix evolution. Section 4.2 is dedicated to the study of the asymptotic distributions toward which the $c\bar{c}$ pair relaxes at late time. We discuss the similarities and differences between the QCD case and the corresponding abelian case, and compare the steady states of the QME to thermal equilibrium distributions. The last section of the paper summarizes our conclusions. Seven appendices
contain technical material and additional developments. The explicit expressions of the Lindblad superoperators for the case of a single $Q\bar{Q}$ pair are given in appendices A, B, and C, where an explicit check of the positivity in the evolution is carried out, as well as a discussion on the regularization adopted for the imaginary part of the potential. The main features of the complex one-dimensional potential from [39] are recalled in appendix D. The corresponding dissociation rates are calculated in appendix E. The following appendix F is a standalone study of the evolution under the single superoperator $\mathcal{L}_2$. The last appendix G presents the details of the deviation of the asymptotic momentum distribution obtained solving the QME from a thermal one.

## 2 Equations for the $Q\bar{Q}$ density matrix

In this section, we start by reviewing briefly the quantum master equations (QME) that were established in [14], and that we shall be solving in this paper in a simplified one-dimensional setting. We recall the assumptions under which these equations were derived and comment on their range of validity. In particular we emphasize that they are appropriate in the so-called Quantum Brownian Motion (QBM) regime.<fnref n="1" /> Although this will not be pursued here we also note that they constitute a natural step toward semi-classical approximations,<fnref n="2" /> which have been explored in details in the abelian case in [10]. In fact, the effect of the collisions, which are part of the QME through the $\mathcal{L}_2$, $\mathcal{L}_3$ and $\mathcal{L}_4$ operators, is essentially classical, as already pointed out in [42], and more thoroughly in [14].

The calculations to be presented in this paper concern the time evolution of a single $Q\bar{Q}$ pair immersed in a quark-gluon plasma which is in thermal equilibrium at a given temperature. In the second part of this section, we identify the characteristic temperatures that delineate the various regimes that can be identified in this particular context. The calculations are based on a simplified complex one-dimensional potential that has been adjusted so that numbers that come out of the calculations are in a physically relevant range for the corresponding three-dimensional systems, in particular for the physics of the charmonium [39].

### 2.1 Master equations

The equations that we want to solve describe the time evolution of the density matrix of a $Q\bar{Q}$ pair in a quark-gluon plasma in equilibrium at temperature $T$. In the position representation and in one spatial dimension, the reduced density matrix, obtained by integrating out the plasma degrees of freedom, takes the form

$$ \langle r_1, r_2 | \mathcal{D} | r'_1, r'_2 \rangle = \langle R + s/2, R - s/2 | \mathcal{D} | R' + s'/2, R' - s'/2 \rangle \quad (2.1) $$

where we have introduced center of mass ($R, R'$) and relative ($s, s'$) coordinates of the two members of the $Q\bar{Q}$ pair. In fact we shall only study here the relative motion, and reduce further the density matrix by integrating out the center of mass. Thus, in the space of relative coordinates, and with a slight abuse of notation, we set

$$ \langle s | \mathcal{D} | s' \rangle = \int_R \langle R + s/2, R - s/2 | \mathcal{D} | R + s'/2, R - s'/2 \rangle. \quad (2.2) $$

<footnote n="1">¹In particular we ignore the effect of the energy gap in the various transitions between states [40].</footnote>

<footnote n="2">²Such approximations will be discussed more thoroughly in a forthcoming paper [41].</footnote>
Thus defined, the reduced density matrix $\mathcal{D}$ remains a matrix in color space which we may conveniently decompose into singlet and octet components [14]

$$ \mathcal{D} = \mathcal{D}_s |s\rangle \langle s| + \mathcal{D}_o \sum_{i=1}^{N_c^2-1} |o^i\rangle \langle o^i|, \qquad (2.3) $$

where $|s\rangle$ denotes a color singlet state and $|o^i\rangle$ a projection of a color octet state in the Hilbert space of the quark-antiquark pair.<fnref n="3" />

The equations satisfied by $\mathcal{D}_s$ and $\mathcal{D}_o$ can be written in the form [14]

$$ \frac{d}{dt} \begin{pmatrix} \mathcal{D}_s \\ \mathcal{D}_o \end{pmatrix} = \mathcal{L} \begin{pmatrix} \mathcal{D}_s \\ \mathcal{D}_o \end{pmatrix}. \qquad (2.4) $$

The Liouville operator $\mathcal{L} = \sum_{i=0}^{4} \mathcal{L}_i$ (also referred to as a superoperator) is a 2 × 2 matrix in the space spanned by $\mathcal{D}_s$ and $\mathcal{D}_o$.<fnref n="4" /> The expressions of the superoperators $\mathcal{L}_i$, before color projection and center of mass elimination, are given by

$$ \begin{aligned} \mathcal{L}_0 \mathcal{D} &= -i [H_Q, \mathcal{D}] \\ \mathcal{L}_1 \mathcal{D} &= -\frac{i}{2} \int_{xx'} V(x-x') [n_x^a n_{x'}^a, \mathcal{D}] \\ \mathcal{L}_2 \mathcal{D} &= \frac{1}{2} \int_{xx'} W(x-x') (\{n_x^a n_{x'}^a, \mathcal{D}\} - 2 n_x^a \mathcal{D} n_{x'}^a) \\ \mathcal{L}_3 \mathcal{D} &= -\frac{i}{4T} \int_{xx'} W(x-x') \left( \dot{n}_x^a \mathcal{D} n_{x'}^a - n_x^a \mathcal{D} \dot{n}_{x'}^a + \frac{1}{2} \{\mathcal{D}, [\dot{n}_x^a, n_{x'}^a]\} \right) \\ \mathcal{L}_4 \mathcal{D} &= \frac{1}{32T^2} \int_{xx'} W(x-x') (\{\dot{n}_x^a \dot{n}_{x'}^a, \mathcal{D}\} - 2 \dot{n}_x^a \mathcal{D} \dot{n}_{x'}^a). \end{aligned} \qquad (2.5) $$

In these equations, $H_Q$ is the free quark Hamiltonian (i.e., the kinetic energy), $V$ and $W$ are respectively the real and imaginary parts of a complex potential, and $n_x^a$ is the color charge density operator at location $x$:

$$ n_x^a = \delta(x - \hat{r}) t^a \otimes \mathbb{I} - \mathbb{I} \otimes \delta(x - \hat{r}) \tilde{t}^a, \qquad (2.6) $$

where $\hat{r}$ is the position operator, and the two components of the tensor product refer respectively to the Hilbert spaces of the heavy quark for the first component and the heavy antiquark for the second component. In eq. (2.6), $t^a$ is a color matrix in the fundamental representation of SU(3) and $\tilde{t}^a$ denotes its transpose. The matrix elements of the time derivative of $n_x = \delta(x - \hat{r})$ are given by [14]

$$ \langle r | \dot{n}_x | r' \rangle = -\frac{1}{2iM} \{ [\partial_r \delta(r-r')] \cdot \partial_x[\delta(x-r') + \delta(x-r)] \}, \qquad (2.7) $$

where $M$ is the mass of the heavy quark. Thus the presence of a time derivative of the density in eqs. (2.5) brings in a factor $1/M$ and two spatial derivatives in the final expressions of

<footnote n="3">³Contributions $\propto |o^i\rangle\langle s|$ have been discarded as we consider only states that are invariant under color rotations (see discussion in [14], appendix D).</footnote>

<footnote n="4">⁴The matrix elements of $\mathcal{L}$ in this two-dimensional space will be denoted $\mathcal{L}^{ss}$, $\mathcal{L}^{so}$, $\mathcal{L}^{os}$ and $\mathcal{L}^{oo}$.</footnote>
the corresponding coordinate space representation of the superoperators (see appendix C for details).

In writing eqs. (2.5) we have used a specific time discretization discussed in appendices A and B of [14]. This discretization allows us to deduce the sum of operators $\mathcal{L}_2 + \mathcal{L}_3 + \mathcal{L}_4$ by making, in the expression of $\mathcal{L}_2$ given in eq. (2.5), the following substitutions

$$n_x^a \rightarrow \left(n_x^a - \frac{i}{4T} \dot{n}_x^a\right), \quad n_{x'}^a \rightarrow \left(n_{x'}^a - \frac{i}{4T} \dot{n}_{x'}^a\right). \qquad (2.8)$$

The operators $\mathcal{L}_3$ and $\mathcal{L}_4$ that result from this substitution collect the terms that are respectively linear and quadratic in the time derivative of the density. This procedure has the advantage to put the evolution equation for the density matrix into a Lindblad form. This form is explicit in the expression of $\mathcal{L}_2$ (and that of $\mathcal{L}_4$), but not in the expression of $\mathcal{L}_3$. In [14] it was argued that the operator $\mathcal{L}_4$, which contains two derivatives of the density, is subleading (by a factor $1/MT$) and it was not considered in the various applications discussed there. This operator is however necessary to guarantee the positivity of the evolution (i.e. making sure that the density matrix remains a positive definite matrix at all times or, equivalently, that the probability to find the system in any arbitrary state is non-negative) and it will be included in the present work. The question of positivity and the role of this operator is discussed in section 4.1 and in appendix A. The explicit expressions of the operators $\mathcal{L}_i$ used in this paper are given in appendix C in terms of derivatives of the complex potential.

In eqs. (2.5) the operator $\mathcal{L}_0$ describes the free motion of the quarks (i.e. their kinetic energy), while the operator $\mathcal{L}_1$ describes the attractive or repulsive effects of the real potential $V$. Together, the operators $\mathcal{L}_0$ and $\mathcal{L}_1$ describe a unitary evolution of the $Q\bar{Q}$ system. In contrast, the other operators, which involve the imaginary potential $W$, imply a non-unitary evolution. The operator $\mathcal{L}_2$ captures the effect of the random kicks due to collisions of the heavy quarks with the plasma constituents which lead to diffusion in momentum space and collisional decoherence. The operator $\mathcal{L}_3$, linear in the velocity, accounts for friction and dissipation. The operator $\mathcal{L}_4$, brings in corrections which are quadratic in the velocity and which guarantee the positivity of the evolution.<fnref n="5" />

In deriving these equations the fact that the mass of the heavy quarks is large, i.e. $M \gg T$, was exploited. This allowed to treat the heavy quarks non relativistically, and to simplify their mutual interaction, as well as their interactions with the medium constituents (essentially keeping the Coulomb interaction and ignoring the magnetic interactions). Aside from screening effects, which are included in the real potential $V$, the other interactions of the heavy quarks with the medium constituents can be viewed as collisions. These involve essentially one-gluon exchanges that take place over a typical time $\tau_D \sim 1/m_D$, where the Debye screening mass $m_D$ characterizes the plasma collective dynamics, such as the screening effects included in $V$. An important step in the derivation rests on the fact that we are interested on the evolution of the $Q\bar{Q}$ pair on times scales that are large compared to $\tau_D$. Under such circumstances, the low-frequency part of the plasma response is filtered out, and one can treat the collisions as if they were essentially instantaneous (hence their representation by an imaginary “potential”).

<footnote n="5">⁵See remark after table 2 in appendix B.</footnote>
## 2.2 Various regimes and characteristic temperatures

We shall now elaborate more on the specific features of the complex potential that we use in this paper [39]. As already mentioned, its form and parameters are adjusted so as to produce numbers (for binding energies, damping rates, etc.) that are of comparable magnitudes with those one would obtain in a more realistic situation based on a three dimensional potential. Note that the constraints of the fitting procedure make the potential well defined only above some non-vanishing temperature ($T \gtrsim 120$ MeV), and still appears to present some non-physical artifact below $T = 175$ MeV. As we shall restrict our studies to $T \gtrsim 200$ MeV, which corresponds to the range of validity of the master equations recalled in the previous subsection — as we shall verify shortly —, these non-physical features will be fairly innocuous.

The potential to be used departs from that obtained in weak coupling, including HTL resummations, by taking into account in particular string breaking at large distance or other non perturbative information obtained from lattice calculations. However, its main component remains inspired by the form of the weak coupling potential in three dimensions, whose general expression, for a quark-antiquark pair, reads<fnref n="6" /> [42] (see also [43])

$$V_{\text{HTL}}(r, T) + iW_{\text{HTL}}(r, T) = -\alpha \left[ m_D + \frac{e^{-m_D r}}{r} \right] + i\alpha T \phi(m_D r), \quad (2.9)$$

where

$$\phi(x) = 2 \int_{0}^{\infty} dz \frac{z}{(z^2 + 1)^2} \left( 1 - \frac{\sin(zx)}{zx} \right). \quad (2.10)$$

The function $\phi(x)$ increases monotonously from 0 to 1 as $x$ runs from 0 to $\infty$, and it behaves as $\phi(x) \sim x^2 \ln 1/x$ at small $x$. This particular small distance behavior is characteristic of the underlying dipolar transitions associated with gluon emission or absorption.

The complex potential (2.9) depends on the plasma properties only through the Debye screening mass $m_D$. As mentioned in the previous subsection, in the regime where the master equations (2.5) are valid, only the static response of the plasma is needed, and this is indeed entirely controlled by the Debye mass. At weak coupling $m_D \sim gT$ where $g$ is the gauge coupling, supposed to be small, so that $m_D \ll T$. However, since at the temperatures of interest $g$ can be of order unity, the magnitude of the Debye mass can be of the same order as the temperature. We shall therefore set $m_D = CT$, with $C$ a coefficient of order unity (in practice, $C \simeq 2$).

The real part of the potential $V_{\text{HTL}}(r, T)$ in eq. (2.9) exhibits screening at a scale characterized by the Debye mass $m_D$. At very short distance, i.e. for $r$ such that $m_D r \ll 1$, screening plays no role and the potential indeed identifies with a Coulomb potential, $V_{\text{HTL}}(r) \simeq -\alpha/r$. At large distance, the potential goes to a constant, $V_{\text{HTL}}(r) \to -\alpha m_D$. This constant can be associated to a screening correction to the heavy quark mass. It controls the edge of the continuum where a bound state becomes degenerate with the energy of two independent quarks with mass $M - \alpha m_D/2 \equiv M(T)$.

<footnote n="6">⁶As written the complex potential of eq. (2.9) is that of an abelian theory. In QCD, the actual potential depends on the color state of the $Q\bar{Q}$ pair on which it is acting. To get, for instance, the potential acting on a color singlet $Q\bar{Q}$ pair, one needs to multiply $V + iW$ in eq. (2.9) by $C_F = (N_c^2 - 1)/2N_c$. The precise color factors entering the superoperators of eq. (2.5) are given in appendix C. See also appendix D.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. Charmonium mass spectrum of the S states for the one-dimensional potential (solid lines) obtained from the expectation values of the Hamiltonian compared to the three-dimensional case (dashed lines).</figcaption></figure>

These generic features of the weak HTL potential are also found for the real potential of ref. [39] that will be used for the numerical simulations. They are well illustrated in figure 1 for this specific 1D-potential. The Debye mass parameterization is taken directly from [44] and satisfies $m_D \approx CT$ with $C \approx 2$ for $T > 200$ MeV. Screening manifests itself by the decrease of the binding energy of the various bound states with increasing temperature.<fnref n="7" /> The binding energy $E_n(T)$ is defined as the eigenvalue of the non-relativistic Hamiltonian

$$H = \frac{P^2}{M} + V(T, r) - V(T, r \to \infty), \quad H|\psi_n\rangle = E_n|\psi_n\rangle \quad (2.11)$$

where the subtraction of $V(T, r \to \infty)$ removes the screening correction to the heavy quark mass. The total mass, $M_n$ of the $Q\bar{Q}$ system in an eigenstate $|\psi_n\rangle$ of the non-relativistic Hamiltonian $H$, is given by

$$ \begin{aligned} M_n &= 2M + \langle \psi_n | \frac{P^2}{M} + V(T, r) | \psi_n \rangle \\ &= E_n + 2M + V(T, r \to \infty) = E_n + 2M(T) \end{aligned} \quad (2.12) $$

When the binding energy vanishes, $M_n = 2M + V(r \to \infty) = 2M(T)$, generalizing the expression given above for the effective mass $M(T)$. At the temperature for which the binding energy vanishes, the mass of the $Q\bar{Q}$ system is equal to twice the effective mass of two heavy quarks at the same temperature (and in the same overall color state). This temperature marks the disappearance of the bound state. It will be referred to as the “melting temperature”,<fnref n="8" /> and denoted $T_{melt}$. Note that $T_{melt}$ depends on the state considered: it decreases as the vacuum binding energy of the state decreases. Thus, in figure 1, the 2S melts at a lower

<footnote n="7">⁷In fact, for the particular potential that we use, the parameter $V(T, r \to \infty)$ in eq. (2.12) is kept constant up to $T = 175$ MeV. This results in a small artifact: the binding energies slowly increases with the temperature up to $T = 175$ MeV (see figure 2), beyond which the binding energies decrease as expected.</footnote>

<footnote n="8">⁸In order to avoid dealing with too loosely bound states, $T_{melt}$ will be operationally defined as the temperature at which $E_{binding} \lesssim 1$ MeV.</footnote>
temperature than the 1S ($T_{\text{melt}}$ is approximately 440 MeV for the 1S and 250 MeV for the 2S). As explained in appendix D we indeed refer to the bound states of the 1D potential by using the three-dimensional conventions, S and P waves corresponding to 1D-even and 1D-odd states, although there is of course no angular momentum in one dimension. Note that the screening correction to the heavy quark mass does not enter the calculation of the binding energy, and plays therefore no role in the determination of $T_{\text{melt}}$.

The real part of the potential is involved in the unitary part of the evolution. In terms of the Liouville operators acting on the density matrix, this unitary evolution is driven by $\mathcal{L}_0+\mathcal{L}_1$. In the low temperature regime, where the bound states can be clearly identified, these bound states form a convenient basis to analyze the density matrix.<fnref n="9" /> In this low temperature regime, the motion of the heavy quarks in the bound state is fast compared to the plasma collective dynamics. Let us characterize the time scale of the relative motion by $\tau_{\text{BS}} \sim 1/E_{\text{BS}}$, where $E_{\text{BS}}$ is the binding energy. As recalled earlier, the time scale of plasma dynamics is $\tau_D \sim 1/m_D$, with $m_D$ the Debye mass. The condition $\tau_{\text{BS}} \sim \tau_D$, or equivalently $E_{\text{BS}} \sim m_D$, characterizes the transition region between the “low temperature” and the “high temperature” regimes. We shall call $T_{\text{QBM}}$ the corresponding temperature (see figure 2). Below $T_{\text{QBM}}$, the bound states remain a dominant feature of the dynamics, which can be understood in terms of small modification of the bound state properties and the various transitions between these bound states. In this paper we shall be interested in the high temperature regime, i.e. in the regime where $T > T_{\text{QBM}}$. Then, $\tau_{\text{BS}} > \tau_D$, and the motion of the heavy quarks is slow on the scale of the plasma dynamics. The perturbations caused in the plasma by the motion of the heavy quarks are nearly instantaneously damped. This high temperature regime is often referred to in the literature as the QBM regime. This is the regime for which the equations (2.5) are supposed to be valid. In the present case, with the one-dimensional potential that we use, $T_{\text{QBM}} \simeq 210$ MeV for a $J/\Psi$ state, as can be read from the left panel of figure 2. The transition to the QBM regime is of course not a sharp one, and we could as well define $T_{\text{QBM}}$ as the temperature at which the binding energy and the temperature are approximately equal. From figure 2, one sees that the corresponding temperature is slightly higher, about 240 MeV. In this paper we shall stick to the value defined as above, i.e. 210 MeV, keeping in mind that, at this temperature, effects of bound states can still be visible. We shall find that it is indeed the case with the simulations performed at $T= 200$ MeV.

We turn now to the effects of collisions which are responsible for the non-unitary part of the evolution. In weak coupling regime, the rate of such collisions is of the order of $\alpha T$, the order of magnitude of the imaginary potential of eq. (2.9). Defining the collision time $\tau_R \sim 1/\alpha T$ and the plasma time $\tau_D \sim 1/m_D$, we get for the ratio

$$ \frac{\tau_R}{\tau_D} = \frac{C}{\alpha} \gg 1, \quad (m_D = CT). \qquad (2.13) $$

Thus, independently of the value of the temperature, the collision time is large compared to $\tau_D$ (even in strict weak coupling): relaxation processes governed by collisions are slow on the time scale of the plasma collective dynamics. However, the damping rate of the bound state depends strongly on the size of the $Q\bar{Q}$ pair, in particular it vanishes for small size,

<footnote n="9">⁹Note however that the continuum states play an important role. Depending upon the initial conditions these continuum modes may be populated and be part of the unitary evolution described by $\mathcal{L}_0 + \mathcal{L}_1$.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 2. Left: various energy scales (in GeV) obtained from our 1D-potential construction [39], in the case of a c$\bar{c}$ pair in a 1S state. An estimation of the temperatures $T_{\text{melt}}$, $T_{\text{QBM}}$ and $T_{\text{damp}}$ is also provided (see text for details). Right: same for a b$\bar{b}$ pair.</figcaption></figure>

due to the dipolar nature of the interaction. To take this into account, we calculate the decay width $\Gamma$ of the 1S eigenstate of the temperature dependent potential (see appendix E) and plot it in figure 2. One sees that the decay width increases with temperature. This increase is the combination of several effects: one is the increase of the overall magnitude of the imaginary potential, proportional to $\alpha T$, another is the temperature dependence of the function $\phi(m_D r) \simeq m_D^2 r^2$ of the imaginary potential (see eq. (2.10)). For a parametric estimate which takes approximately these two effects into account, we may use

$$ \Gamma \approx \alpha T (a_0 m_D)^2 \approx 4C^2 T^3 / (\alpha M^2), \qquad (2.14) $$

where $a_0$ denotes the typical size of the bound state. Still another effect is the increase of this bound state size that accompanies the decrease of the binding energy.<fnref n="10" /> Since the size of the bound state increases with the temperature (see e.g. figure 27), the rate $\Gamma$ varies more rapidly with temperature than predicted by eq. (2.14), which appears to be the case as illustrated in figure 33.

The point where the decay width becomes equal to the binding energy defines another characteristic temperature $T = T_{\text{damp}}$, with $T_{\text{QBM}} < T_{\text{damp}} < T_{\text{melt}}$. This temperature marks the beginning of an overdamped regime where one expects the $Q\bar{Q}$ spectral function to become structureless.

In summary, we have identified three important temperatures. The first one is $T_{\text{QBM}}$, defined as the temperature at which the binding energy equals the Debye mass. This temperature marks the onset of the Quantum Brownian Motion regime. The next temperature is $T_{\text{damp}}(> T_{\text{QBM}})$ defined as the temperature at which the binding energy equals the decay width. Above that temperature, it is expected that the $Q\bar{Q}$ spectral function becomes structureless. Finally, above the temperature $T_{\text{melt}}$ the binding energy vanishes. Note that the two temperatures $T_{\text{QBM}}$ and $T_{\text{melt}}$ depend only on the Debye mass and the real part of the potential. Only $T_{\text{damp}}$ depends on collisions, i.e. on the imaginary part of the potential. Finally, although the results of the next section will concern solely the physics of charmonium, we present in the right panel of figure 2 results corresponding to the bottomonium. It can

<footnote n="10">¹⁰Note that at low temperature ($T \lesssim T_{\text{QBM}}$) the calculation overestimates the decay width since it ignores the effect of the energy gap between levels. See [40] for a recent discussion of this effect.</footnote>
be seen that the same pattern is observed and the same characteristic scales emerge, to
within an overall rescaling of the temperature.

3 Main physics results

In this section, we study the time evolution of the density matrix $\langle s | \mathcal{D} | s' \rangle$ of a single c$\bar{c}$ pair, in the form given in eq. (2.2), where the center of mass degrees of freedom are traced out and $s$ and $s'$ refer to relative coordinates.<fnref n="11" /> We examine various contexts. In section 3.1, the initial state is chosen to be the 1S, color singlet, eigenstate of the in-medium Hamiltonian (2.11), and we provide an overview of various aspects of its time evolution in a medium at temperature $T = 300$ MeV. At this temperature, the 1S state is the only bound state. We also investigate the changes in the evolution that occur at the smaller temperature $T = 200$ MeV at which three bound states exist, like in the vacuum. In order to study the overall evolution of the system at the higher temperature $T = 600$ MeV, at which the 1D potential does not sustain any bound state, we choose, in section 3.2, a vacuum 1S state as initial state. We then explore more completely the dependence of the medium temperature $T$, by varying $T$ from 200 to 600 MeV. In section 3.3, we consider the evolution of a compact wave packet which mimics the wave packet created in the gluon fusion picture of c$\bar{c}$ production [45]. Finally, in section 3.4, in order to get closer to the environment of heavy ion collisions, we relax the fixed medium temperature assumption and consider the evolution of a c$\bar{c}$ pair in a time-dependent environment with a temperature following (adiabatically) a Bjorken profile [38].

This section aims at providing an overall physical picture of the fate of a single c$\bar{c}$ pair
in a hot QGP medium whose temperatures are chosen so as to illustrate the various regimes
introduced in the previous section. We recall that the one-dimensional potential used in this
work was constructed with the aim to yield numerical estimates that are phenomenologically
relevant. Hence, the various time and length scales that appear in this section are meant to
be as realistic as can possibly be in a one-dimensional setting. The details of the potential
and of the numerical implementation are presented in appendix D. Let us recall here that
the one-dimensional grid used in the calculations extends from -10 fm to +10 fm with respect
to the origin of the relative distances.

3.1 In-medium 1S singlet initial state

As a first step in our analysis, we consider as initial density matrix the projector on the
1S ground state of the in-medium real potential in the color singlet channel, at the plasma
temperature $T = 300$ MeV. Starting from an eigenstate of the in-medium potential at the
plasma temperature leads to the cancellation of the contribution of $\mathcal{L}_0 + \mathcal{L}_1$ at the initial time
and allows us to focus on the role of the imaginary potential. Besides, only one singlet bound
state exists for $T = 300$ MeV, which simplifies somewhat the analysis. Figure 3 shows the time
evolution of the $s, s'$ matrix elements of the components $\mathcal{D}_s$ and $\mathcal{D}_o$ of the density operator (see
eq. (2.3)). We see that the octet channel, initially absent, is quickly populated (within 1 fm/c)
as a P state (see bottom left panel) via dipolar transitions from the 1S state induced by the
operator $\mathcal{L}^{\text{os}}$. Since the potential is (weakly) repulsive in the octet configuration, the relative

<footnote n="11">¹¹Note that the center of mass momentum does not enter explicitly the expressions of the operators $\mathcal{L}_i$, except $\mathcal{L}_4$ (see appendix C).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time evolution (from left to right: $t = 1, 5$ and $10$ fm/c) of the singlet density matrix $\langle s | \mathcal{D}_s | s' \rangle$ (top panels) and octet density matrix $\langle s | \mathcal{D}_o | s' \rangle$ (bottom panels) in a medium at temperature $T = 300$ MeV. The initial state is an “in-medium” 1S color singlet. For better readability, the color scale is different for each time and the octet densities have been rescaled by a numerical factor, 10 for $t = 1$ fm/c and 5 for $t = 5$ fm/c.</figcaption></figure>

distance between c and $\bar{c}$ increases with time, leading to a spreading of the density along the $s = s'$ axis. Since an octet state can be locally converted back into a singlet state via the $\mathcal{L}^{so}$ operator, the same spreading occurs for $\mathcal{D}_s(s, s)$ as well, in spite of the attractive nature of the potential in the singlet channel. In fact, in the singlet channel, a peak remains visible in the central region of the corresponding plots, i.e. for $s$ and $s'$ within the range of the singlet potential. This peak reflects the survival of bound states until the late stage of the evolution.

One also observes in figure 3 that the reduced density matrix $\langle s | \mathcal{D} | s' \rangle$ tends to become diagonal at late time. This is the result of collisional decoherence induced by the operators $\mathcal{L}_2$, which brings the system closer and closer to a classical system with a diagonal density matrix $\langle s | \mathcal{D} | s' \rangle \propto \delta(s - s')$ (see appendix F). The operator $\mathcal{L}_3$ introduces friction that eventually balances the effect of $\mathcal{L}_2$ so that the off-diagonal elements of $\langle s | \mathcal{D} | s' \rangle$ remain significant within a range of the order of the thermal wavelength, $|s - s'| \lesssim \lambda_{\text{th}} = \frac{1}{\sqrt{MT}} \approx 0.3$ fm. The two features of the reduced density matrix $\langle s | \mathcal{D} | s' \rangle$ that we have just discussed, namely the spreading along the $s + s'$ axis, and the shrinking along the $s - s'$ axis (decoherence) are found generically in all our simulations (in agreement with other results in the literature, see e.g. [34, 36]).

Color equilibration is another feature of the present QME. As a global measure, we consider the singlet and octet weights, defined respectively as tr $\mathcal{D}_s$ and tr $\mathcal{D}_o$. These weights are plotted on the left panel of figure 4, which shows that they become equal after a time $t \approx 10$ fm/c. Note that since $1 = \text{tr } \mathcal{D} = \text{tr } \mathcal{D}_s + 8\text{tr } \mathcal{D}_o$, the equality tr $\mathcal{D}_s = \text{tr } \mathcal{D}_o$ implies that there are, at that time, eight times more octet pairs than singlet pairs in the system,<fnref n="12" />

<footnote n="12">¹²A word of caution about the wording is perhaps needed here in order to avoid confusion. There is only one pair in the system, but the results would be identical for an arbitrary number of independent pairs. It is often convenient to think about the system as a collection of pairs, as done implicitly here and repeatedly throughout, but of course what is referred to is a probability for one pair.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 4. Left: time evolution of the singlet weight (tr $D_s$, solid) and the octet one (tr $D_o$, long-dashed). The dashed line represents $\exp(-\Gamma_{1S}t)$. Same initial conditions as in figure 3. Right: spatial densities $\rho_s(r,t)$ (solid) and $\rho_o(r,t)$ (dashed) for various times: 0 fm/c (black), 2 fm/c (blue), 5 fm/c (orange), 10 fm/c (green) and 20 fm/c (red). The medium temperature is $T = 300$ MeV.</figcaption></figure>

and according to figure 3, these pairs are to be found in continuum states. At very early time, the evolution of the density matrix is dominated by that of its major component, the in-medium 1S state. Indeed an exponential fit at very early time ($t \lesssim 0.5$ fm/c) yields a decay rate $\Gamma_{1S} \approx 0.1$ GeV, which is in good agreement with the value estimated in appendix E for the in-medium 1S state (see figure 33). One notes however that the actual behavior of the singlet weight rapidly deviates from this initial exponential law, in agreement with expectations based on the analysis presented in appendix F (see eq. (F.7)) of the evolution induced by the operator $\mathcal{L}_2$ alone.

A more differential information can be obtained from the diagonal elements of the density matrix, i.e., the local densities

$$ \rho_{s/o}(r,t) = \mathcal{D}_{s/o}(r,t). \qquad (3.1) $$

These are plotted in the right panel of figure 4. Color equilibration manifests itself as $\rho_s(r,t)$ becoming equal to $\rho_o(r,t)$, which occurs rapidly at large relative distance: $\rho_s \approx \rho_o$ for $|r| \gtrsim 2$ fm already at early time. Between 5 and 10 fm/c, the densities reach the box boundaries ($\pm 10$ fm/c) and level off toward their asymptotic values ((1/9)×(1/20) fm$^{-1} \approx 5.5 \times 10^{-3}$ fm$^{-1}$). At small distance, the dipolar nature of the imaginary potential leads to a reduced equilibration rate, which delays color equilibration as illustrated in the figure. We then arrive at the physical picture of a decaying singlet bound state in a growing and equilibrating background of unbound $c\bar{c}$ pairs: at $t = 20$ fm/c the distribution at large relative distance reaches a stationary pedestal, while a trace of a bound state stills remains at small relative distance and is also seen in the asymptotic distribution, as discussed in subsection 4.2.

Further insight into the equilibration process can be gained from the momentum distributions $\rho_{s,o}(p,t)$, defined as

$$ \rho_{s/o}(p,t) = \frac{1}{2\pi} \int dy e^{-ipy} \int dr \mathcal{D}_{s/o} \left( r + \frac{y}{2}, r - \frac{y}{2}, t \right), \qquad (3.2) $$

and displayed in figure 5. One can see there that the momentum densities evolve gradually toward the thermal distribution $\exp(-p^2/2M_{\text{red}}T) = \exp(-p^2/MT)$, which reflects the
<figure><img src="image_8.png" /><figcaption>Figure 5. Time evolution of the singlet (left) and octet (right) momentum densities for $t = 0$ fm/c (black), 1 fm/c (blue), 2 fm/c (orange), 5 fm/c (green) and 20 fm/c (red). The initial condition (1S in-medium state) and the temperature ($T = 300$ MeV) are the same as in figure 3. Each curve has been rescaled to unity at $p = 0$ for better readability. The classical Maxwell distribution $\exp(-p^2/MT)$ — with $M = m_c = 1.469$ GeV — is also shown for comparison (dashed). A fair agreement with the actual distribution at $t = 20$ fm/c holds for $p \lesssim \lambda_{th}^{-1} \approx 0.7$ GeV/c.</figcaption></figure>

statistical population of the continuum momentum states. Note that for $t = 20$ fm/c, the thermal (Maxwell-Boltzmann) distribution is well reproduced up to momenta of the order of $\sqrt{MT} \simeq 0.7$ GeV, corresponding to a thermal wavelength $\sim 0.3$ fm. However the tail at large $p$ deviates systematically from the Maxwell-Boltzmann distribution. This behavior is further discussed in section 4.2.2 and appendix G where a more detailed analysis of the late-time momentum distribution is given.

A different perspective on the time evolution of the density matrix can be obtained by projecting the singlet density operator $D_s$ on various charmonium bound states, that we denote by $\Phi_n$. Recall that we have two types of bound states: in-medium bound states, i.e., eigenstates of the in-medium screened potential, and vacuum bound states. The probability to find the system at time $t$ in one of these bound states is given by

$$p_n(t) = \langle \Phi_n | D(t) | \Phi_n \rangle = \iint ds \, ds' \Phi_n^*(s) D_s(s, s', t) \Phi_n(s'). \quad (3.3)$$

This is of course a conditional probability that depends on the choice of the initial density matrix, the survival probability corresponding to the choice $D(0) = |\Phi_n\rangle\langle\Phi_n|$. The evolution of the survival probability of the in-medium 1S state (the only bound state for $T = 300$ MeV) is shown in figure 6 (left panel).

At early time, the evolution of this in-medium eigenstate is dominated by the imaginary potential, that is, by collisions. These naturally decrease the survival probability as continuum states are being populated. As already mentioned, at very early time, $t \lesssim 0.5$ fm/c, this initial depletion of the 1S state is well captured by an exponential decay law, with the decay rate $\Gamma_{1S} \approx 100$ MeV. At later times, say up to $t \lesssim 6$ fm/c, the survival probability of the initial state still follows an approximate exponential decay law, but with a smaller effective rate of $\approx 77$ MeV, indicated by the dotted curve. We do not have a simple explanation for the emergence of this approximate exponential law, besides the fact that it results from the combined contributions of several operators $\mathcal{L}_i$. We show in appendix F that the action of $\mathcal{L}_2$ alone predicts a definite deviation from the exponential $\exp(-\Gamma_{1S}t)$ toward an algebraic
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Left: time evolution of the probability $p_{1S}$ of the in-medium 1S state (solid line). The initial condition and temperature ($T = 300$ MeV) are the same as in figure 3. The dashed curve represents the singlet weight $\text{tr} \mathcal{D}_s$ (from figure 4). Right: the probabilities of vacuum eigenstates $p_{1S}$ (orange), $p_{1P}$ (green) and $p_{2S}$ (red). On both panels, the horizontal long-dashed lines indicate the respective asymptotic values (see section 4.2.3).</figcaption></figure>

law (see e.g. eq. (F.26) and figure 37). Beyond 6 fm/c, the trend progressively changes and one starts observing the beginning of an asymptotic regime where the $c\bar{c}$ relative motion reaches some form of equilibrium. After 20 fm/c, about 10% of the singlet weight can be attributed to the surviving in-medium 1S state, a value which has been checked to scale like the inverse of the box size $L$, provided that $L$ is significantly larger than the potential range (see the discussion in section 4.2.1, and in particular figure 24).

On the right panel of figure 6, we display the projection of the density matrix on the three vacuum eigenstates, i.e., the probabilities $p_{1S}$, $p_{1P}$ and $p_{2S}$, the initial state being as before the in-medium 1S state. The motivation for this analysis is that if the medium were to freeze-out at say, 10 fm/c, these projections would yield the actual populations, at that particular time, of the various (vacuum) states. The in-medium 1S state has initially a projection to the vacuum 1S state nearly equal to unity and a small component (of the order of 5%) on the vacuum 2S. Because of the large overlap of the in-medium and vacuum 1S states, the evolution of the vacuum $p_{1S}$ does not differ significantly from the in-medium $p_{1S}$ displayed on the left panel. The probability $p_{2S}$ has a more complex evolution than the $p_{1S}$, reflecting the effective couplings of the three vacuum eigenstates, in addition to their decay to continuum states. Such a coupling is obvious for the 1P state, which is not present in the initial state. It is only populated through the chain of transitions (1S singlet $\rightarrow$ odd octet $\rightarrow$ even octet $\rightarrow$ 1P singlet) found at the 3rd order in time-dependent perturbation theory (this explains in particular the slow growth (as $t^3$) of the probability $p_{1P}$). We note that $p_{2S}$ and $p_{1P}$ reach a common value when $t \gtrsim 10$ fm/c, and then slowly converge, as well as $p_{1S}$, to their respective asymptotic values. That $p_{2S}$ and $p_{1P}$ reach a common value before reaching their asymptotic limit can be understood from the fact that the density matrix becomes rapidly diagonal, as shown in figure 3. If it were completely diagonal, $\mathcal{D}(s, s') = \rho(s)\delta(s-s')$, then the projection on any state would be almost the same, since $\langle\Psi|\mathcal{D}|\Psi\rangle = \int dr\rho(r)|\Psi(r)|^2$ and $\rho(r)$ is nearly constant, except in a small region around the origin. Another important factor that contributes is the fact that the 1P and 2S states become rapidly delocalized, and are therefore subject to frequent transitions toward octet states. A more detailed discussion of
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left: same as figure 4 (right panel) for $T = 200$ MeV, with one extra time ($t = 60$ fm/c, purple). Right: the probabilities $p_{1S}$ (orange), $p_{1P}$ (green) and $p_{2S}$ (red) corresponding to the in-medium bound states. The horizontal dashed line indicates the asymptotic value of $p_{1S}$. The asymptotic values of $p_{1P}$ and $p_{2S}$ are already reached before $60$ fm/c. The black dashed line illustrates the exponential decay law with $\Gamma_{1S}$ evaluated from eq. (E.2).</figcaption></figure>

the late time behavior and the corresponding asymptotic state is presented in section 4.2.3.

It is finally instructive to contrast the evolution of the in-medium 1S state at temperature $T = 300$ MeV with the corresponding evolution taking place at $T = 200$ MeV. According to figure 2, this temperature lies on the borderline of the QBM regime, but as can be seen in figure 7, the results are qualitatively the same as the ones displayed in figure 4 for $T = 300$ MeV, with a lowering of the transition probabilities (and hence a dilatation of the time scale), and a stronger peak at small $r$ remaining at late time. For $T = 200$ MeV, the 1D potential still supports three bound states, with the following binding energies: 0.474 GeV (1S), 0.221 GeV (1P), and 0.103 GeV (2S) (see appendix D). The probabilities $p_{1S}$, $p_{1P}$ and $p_{2S}$, corresponding to these in-medium bound states, are displayed on the right panel. One notes that the decay of $p_{1S}$ follows more accurately the exponential law than at higher temperature, with the value $\Gamma_{1S} = 20$ MeV corresponding to $T = 200$ MeV. One notes also that $p_{1P}$ and $p_{2S}$ reach their common asymptotic value before the end of the evolution, as already observed for $T = 300$ MeV.

## 3.2 1S vacuum singlet initial state

In this subsection, we consider as initial state a 1S vacuum state, and analyze its evolution in terms of the projections of the density matrix on the various vacuum bound states, i.e., in terms of the vacuum probabilities $p_{1S}$, $p_{1P}$ and $p_{2S}$, for temperatures in the range $T = 200-600$ MeV.

We start by considering the evolution of the singlet density at $T = 600$ MeV, a temperature at which the in-medium potential does not sustain any bound state. As shown in figure 8 the qualitative features of the evolution are quite similar to those observed for $T = 200$ MeV in figure 7 (left panel) and for $T = 300$ MeV in figure 4 (right panel). The variations with the temperature can be understood simply in terms of the collisions: at $T = 600$ MeV, the system suffers harder and more frequent collisions, the $\bar{c}\bar{c}$ relative distance increases faster and the central peak disappears earlier, as compared to the $T = 300$ and $200$ MeV cases. In fact beyond $t \simeq 5$ fm/c, no peak survives at small $r$ that would signal the presence of a bound state, in line with the fact that the in-medium potential is not binding at $T = 600$ MeV.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Time evolution of the singlet density $\rho_s(r) = \langle r | D_s | r \rangle$ for an initial 1S vacuum state and two different medium temperatures: $T = 300$ MeV (solid curves) and $T = 600$ MeV (dashed curves).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Up-left: time evolution of the probabilities $p_{1S}$ (orange), $p_{1P}$ (green) and $p_{2S}$ (red) of the vacuum singlet states in a medium at temperature $T = 200$ MeV (solid lines). The initial state is the 1S vacuum state. The horizontal long-dashed lines represent the asymptotic values. Up-right, low-left, and low-right: same for $T = 300, 400$ and $600$ MeV.</figcaption></figure>

The time evolution of the vacuum probabilities $p_n$ are shown in figure 9 for four different values of the temperature. At early times, $p_{1S}$ exhibits an exponential decay, while the growth of $p_{1P}$ and $p_{2S}$ reflects the generation of the respective populations via various transitions, as already discussed. After this initial transient stage, $p_{1S}$ and $p_{2S}$ appear to decrease with similar
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Half-life of an initial in-medium (resp. vacuum) 1S-state estimated from figures 6 and 7 as well as from other evolutions not shown (resp. figure 9), in solid red line (resp. long-dashed blue line), as a function of the temperature and compared to the estimation made using the decay rate of either the in-medium 1S state (red dashed) or the vacuum 1S state (blue dot-dashed) evaluated in appendix E (see figure 33).</figcaption></figure>

rates. As the temperature increases, the same pattern develops, while the three probabilities become rapidly of comparable magnitudes, eventually reaching a nearly common value at late time, for the reasons invoked earlier about the diagonal structure of the density matrix and the delocalized nature of the various states. Finally, note that the oscillatory behavior of the 2S weight found at early time for the lowest temperature ($T = 200$ MeV) disappears for temperature $\approx T_{\text{damp}} \approx 300$ MeV and beyond, in agreement with the interpretation of $T_{\text{damp}}$ given in section 2 as marking the onset of an overdamped regime.

The time dependence of the various probabilities in figure 9 appears to follow, for $T \gtrsim 300$ MeV, a common pattern, with the plots at different temperatures being related by a simple scaling of the time scale with the decay rate. In fact, because the exponential decay with the appropriate in-medium or vacuum decay width as calculated in appendix E, holds only at very early time, and can be affected with transient effects, a more accurate characterization of the decay of the (dominant) 1S state is provided by its half-life $\tau_{1/2}$, that is the time it takes for $p_{1S}$ to drop by a factor 2. The half-life derived from the QME is plotted in figure 10. We found that the plots of figure 9 that correspond to temperatures $T \gtrsim 300$ MeV approximately overlap once plotted as a function of the reduced time $t/\tau_{1/2}$.

### 3.3 Initial compact wave packet

As a somewhat more realistic initial condition, we now investigate the case where the initial system is described by a compact wave packet. If one imagines the $c\bar{c}$ pairs as being produced in a gluon fusion process one expects an initial mixture of singlet and (dominantly) octet pairs. We shall examine separately these two color components of the initial density matrix,
<figure><img src="image_14.png" /><figcaption>Figure 11. Left: rms radius of $D_s(r, r, t)$ for the 1S initial state (dashed curves) and the compact initial wave packet (solid curves), for three values of the medium temperature, $T = 200, 300$ and $600$ MeV (bottom to top, i.e. blue, orange, green). Right: same, focusing on the first $3$ fm/c of the evolution.</figcaption></figure>

respectively $D_s$ and $D_o$, for which we choose the following simple Gaussian forms<fnref n="13" />

$$D_s(s, s') \propto \exp\left(-\frac{s^2 + s'^2}{2\sigma^2}\right) \qquad D_o(s, s') \propto s s' \exp\left(-\frac{s^2 + s'^2}{2\sigma_o^2}\right). \qquad (3.4)$$

The parameter $\sigma$ of $D_s$ is adjusted so that the rms radius of the density $D_s(r, r)$ is $\frac{\sigma}{\sqrt{2}} \approx 0.12$ fm (i.e. less than one half of the 1S rms radius (see table 3)), leading to $p_{1S} \approx 0.74$ and $p_{2S} \approx 0.12$. For $D_o$, $\sigma_o$ is chosen such that the rms of the momentum distribution, has the same value as for the initial singlet Gaussian state, rms($p$) = 0.83 GeV/c, corresponding to $\sigma_o \approx 0.29$ fm.<fnref n="14" />

The time evolution of the rms radius of $D_s(r, r, t)$ is displayed in figure 11 and compared with that of the rms of the vacuum 1S initial state. One can see, in particular in the right panel, that the size of the distribution corresponding to the 1S state does not change much for the first half fm/c. This is because the 1S state is, as we have already indicated, nearly an eigenstate of the in-medium Hamiltonian. When collisions set in the system starts to expand, the faster the higher the temperature. When the initial condition corresponds to a compact wave packet one observes an initial rapid growth due its natural spreading induced by the operator $L_0$. This initial growth, of purely quantum mechanical origin, is essentially independent of the temperature, and lasts for less than about 1 fm/c, at which point the collisions take over and lead to smooth evolution.

The role of $L_0$ in the early stage of the evolution can also be seen in the probabilities $p_n(t)$, as illustrated in the right panel of figure 12. The early transient evolution of the compact initial singlet state is dominated by quantum mechanics, and transitions between the three states, as revealed for instance by the rapid population of the 1P state originally absent. A distinctive feature of the compact-state evolution seen on the right panel of figure 12, and absent from figure 9, is the early time evolution of $p_{1S}$ and $p_{2S}$ which appears to be mostly quadratic in time during the first 0.5 fm/c. Such a behavior is known to originate from the

<footnote n="13">¹³The octet state is taken here in a “P wave”, one of the possible quantum numbers for the $c\bar{c}$ final state resulting from the gluon fusion process at LO.</footnote>

<footnote n="14">¹⁴The specific choice of these values of $\sigma$ and $\sigma_o$ is made so that the two wave packets have initially the same kinetic energy content.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 12. Left: spatial densities $\rho_s(r)$ (solid) and $\rho_o(r)$ (dashed) for various times: 0 fm/c (black), 2 fm/c (blue), 5 fm/c (orange), 10 fm/c (green) and 20 fm/c (red), corresponding to an initial compact octet state. Right: time evolution of the populations of the first three vacuum eigenstates starting either from a compact P-octet wave packet (dashed curves) or a compact S-singlet wave packet (solid curves). The medium temperature is $T = 300$ MeV.</figcaption></figure>

unitary operators $\mathcal{L}_0$ and $\mathcal{L}_1$ (see e.g. [46]). After this transient regime, the probabilities follow a regular pattern analogous to the ones illustrated on figure 9 and eventually reach their asymptotic values at late time.

We turn now to the evolution corresponding to the octet wave packet initial condition. Except for the very early stage, the overall evolution of the densities is similar to that of a singlet initial state (see figure 3). In figure 12 (left panel), the singlet and octet densities are shown for $t = 0, 2, 5, 10$ and $20$ fm/c. One can see that, thanks to the dipolar transitions, within less than two fm/c the singlet density develops a peak at small relative distance, while the peak corresponding to the initial octet density disappears. After that, the evolution of the singlet and octet density is similar to that observed in the right panel of figure 4. The early time behavior of the densities is also reflected in a fast increase of the probabilities $p_{1S}$, $p_{1P}$ and $p_{2S}$, as can be observed on the right panel of figure 12. In particular, one observes a faster growth of P singlet states in the case of an initial P octet state (behaviors in $t^2$ in time-dependent perturbation theory vs $t^3$ for initial singlet state). Generically, S (resp. P) states appear to be less (resp. more) populated when starting from a P octet initial state. Note finally that at late time the various probabilities reach their asymptotic values, independently of the initial conditions.

## 3.4 Charmonium in an expanding plasma

In the previous subsections, we have assumed a medium with a fixed temperature. In order to get closer to potential applications in heavy ion collisions, we now consider an expanding medium with a time-dependent temperature given by the commonly used (see e.g. [32, 34]) Bjorken expansion

$$T(t) = T_0 \left( \frac{1}{\tau_0 + t} \right)^{1/3}, \qquad (3.5)$$

with $T_0 = 600$ MeV the initial medium temperature, and $\tau_0 = 1$ fm/c. We limit the evolution to time $t \le 20$ fm/c in order to fulfill the condition $T(t) > T_{QBM}$ for all times ($T(20$ fm/c) $\approx$
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Left: spatial densities for several times — 0 fm/c (black), 2 fm/c (blue), 5 fm/c (orange), 10 fm/c (green) and 20 fm/c (red) — as a function of the relative distance for a singlet compact initial state; full (resp. dashed) lines are for the singlet (resp. octet) density. Right: same for the P octet compact initial state.</figcaption></figure>

217 MeV). As for the initial state, we consider the two compact states discussed in the previous subsection, namely the singlet and the octet wave packets in eq. (3.4).

The time evolution of the densities $\rho_s(r)$ and $\rho_o(r)$ is shown on figure 13. A question of phenomenological interest is whether the peak near the origin in the singlet density $\rho_s$ survives at the end of the evolution, as it is correlated with significant weights of the vacuum eigenstates. Figure 13 shows that this is clearly the case: the magnitude of the peak at the end of the expansion, is about five times larger than the “equilibrium” value reached in a static medium of constant temperature $T = 200$ MeV, as observed on figure 7. A fully quantitative analysis is delicate because of effects coming from the finite size of the box used in the calculation, but the main point is that a significant peak survives at the end of the expansion. This suggests that, in the present setting, the early hot stage of the QGP may not last long enough to eradicate the peak that arises at small $r$ from the presence of a bound state, as would be the case if the medium, instead of expanding, would remain at its initial temperature $T = 600$ MeV (see figure 8). This conclusion is reinforced by the fact that color has not completely relaxed at the end of the evolution either, as can be inferred from the densities shown in the left panel of figure 13 for $\lesssim 1$ fm. In fact, after $20$ fm/c, still $17\%$ of the $c\bar{c}$ pairs are found in singlet states whereas the equilibrium value would be $1/9 \simeq 11\%$.

The evolution of the densities obtained when one starts with the P octet density matrix, resembles that displayed in the left panel of figure 12. One observes fast transitions toward the singlet channel, building a peak near the origin, as can be seen in the right panel of figure 13. Color relaxation is reached within $\approx 2$ fm/c for $r \gtrsim 2$ fm. The later cooling preserves the peak which appears nearly frozen to a constant magnitude. The final value of $\rho_s(0) \approx 0.03$ fm$^{-1}$ is sizable (compare with figure 12) and of the same order of magnitude as that of the peak obtained at the end of the evolution of the singlet initial state (left panel). This argues in favor of considering both octet and singlet initial states in phenomenological applications.

## 4 Discussion of more technical issues

In this section, we return to some of the main features observed in the analysis of the previous section. In particular, we clarify certain technical aspects of the QME that we use, the physical roles of the various superoperators in the evolution of the density matrix, the
interplay of collisional decoherence and dissipation in driving the system to a steady state at late time, and the specific role of the superoperator $\mathcal{L}_4$ in ensuring the positivity of the evolution. We also discuss characteristics of the asymptotic steady state that is reached at the end of the evolution in a static medium with a constant temperature. We analyze the origin of the small observed deviations of the corresponding momentum distribution from a thermal distribution. We compare more broadly the late-time steady state with thermal equilibrium, and point out interesting differences between abelian and QCD plasma.

# 4.1 The role of various superoperators in the evolution

## 4.1.1 Global overview

As can be seen on the formulae given in appendix C, the non-unitary superoperators $\mathcal{L}_2$, $\mathcal{L}_3$ and $\mathcal{L}_4$ are multiplied by powers of the thermal wavelength, $\lambda_{th} = 1/\sqrt{MT}$, with $\mathcal{L}_2$ independent of $\lambda_{th}$, $\mathcal{L}_3 \propto \lambda_{th}^2$ and $\mathcal{L}_4 \propto \lambda_{th}^4$. This suggests a natural hierarchy for the magnitudes of these superoperators (see also eq. (E.3)), which has been exploited to some extent in the derivation of eqs. (2.5) in the QBM regime. However, the superoperators $\mathcal{L}_3$ and $\mathcal{L}_4$ contain also gradients, of the imaginary potential, and of the density matrix. These have to be taken into account in order to quantify the action of each superoperator on a given density matrix. Thus, in order to gauge (numerically) the relative importance of the action of the various superoperators $\mathcal{L}_k$ on the density matrix, we define the following averages (for $k = 0, \dots, 4$)

$$ \langle \mathcal{L}_k \rangle = \frac{\sum_{i,i'} |\mathcal{L}_k^{ss}[\mathcal{D}_s] + \mathcal{L}_k^{so}[\mathcal{D}_o]|_{i,i'} + (N_c^2 - 1) |\mathcal{L}_k^{os}[\mathcal{D}_s] + \mathcal{L}_k^{oo}[\mathcal{D}_o]|_{i,i'}}{\sum_{i,i'} |\mathcal{D}_s|_{i,i'} + (N_c^2 - 1) |\mathcal{D}_o|_{i,i'}} \quad (4.1) $$

where $i$ (resp $i'$) is a discrete index running over all the values of $s$ (resp. $s'$) coordinates on the grid of the numerical simulations. The linear combinations $\mathcal{L}_k^{ss}[\mathcal{D}_s] + \mathcal{L}_k^{so}[\mathcal{D}_o]$ and $\mathcal{L}_k^{os}[\mathcal{D}_s] + \mathcal{L}_k^{oo}[\mathcal{D}_o]$ inside the $|\cdots|$ have been chosen as they respectively represent the contribution from the superoperator $\mathcal{L}_k$ to the time derivatives of $\mathcal{D}_s$ and $\mathcal{D}_o$ according to eq. (A.1). The particular weightings $|\dot{\mathcal{D}}_s| + (N_c^2 - 1) |\dot{\mathcal{D}}_o|$ and $|\mathcal{D}_s| + (N_c^2 - 1) |\mathcal{D}_o|$ are privileged in the numerator and in the denominator as they correspond to those appearing in the conserved trace. The time evolution of the various $\langle \mathcal{L}_k \rangle$ is illustrated in figure 14 for medium temperatures $T = 200, 300, 400$, and $600$ MeV. This allows us to follow the gradual change of the dynamics from $T \approx T_{QBM}$ to beyond $T_{melt}$. At moderate temperature ($T = 200$ MeV, upper left panel of figure 14), the hierarchy between the various superoperators is well recognized. The decay width $\Gamma$ is small compared to the binding energy $E_{BS}$ of the initial state, or equivalently $\tau_R > \tau_{BS}$ (from figure 2 one gets approximately $\Gamma \approx 50$ MeV, $E_{BS} \approx 500$ MeV so that $\tau_R \approx 4$ fm/c and $\tau_{BS} \approx 0.4$ fm/c). In these conditions, the collisions represent a perturbation of the dynamics which is dominated at early time by $\mathcal{L}_0$ and $\mathcal{L}_1$. Such dynamics lead to a slow evolution of all the $c\bar{c}$ properties (size, internal energy, ...) and of the probabilities to find the system in a given eigenstate. The $c\bar{c}$ system evolves then quantum mechanically for several “periods” $\tau_{BS}$ before the effects of the collisions encoded in $\mathcal{L}_2$ and $\mathcal{L}_3$ take over after $t \approx 10$ fm/c.

With increasing temperature ($T = 300$ MeV, upper right panel of figure 14), $\mathcal{L}_2$ plays an increasing role already at early time, and $\langle \mathcal{L}_2 \rangle$ becomes comparable in magnitude with $\langle \mathcal{L}_0 \rangle$, while the potential becomes more and more screened and plays a negligible role after
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Mean values of the $\mathcal{L}_k$ superoperators (defined according to eq. (4.1)) as a function of time for a medium temperature of 200 MeV (top left panel), 300 MeV (top right panel), 400 MeV (bottom left panel) and 600 MeV (bottom right panel). The initial state is the singlet 1S vacuum state.</figcaption></figure>

a few fm/c. As a consequence, the size of the $c\bar{c}$ pair grows faster and the time over which $\mathcal{L}_0 + \mathcal{L}_1$ governs the evolution shrinks to $t \lesssim 2$ fm/c. Note that $T = 300$ MeV corresponds approximately to the temperature $T_{\text{damp}}$ at which the magnitude of the decay width is comparable to the binding energy (or $\tau_R \approx \tau_{\text{BS}} \approx 1.3$ fm/c).

For even larger temperature (lower panels), collisional effects dominate right from the start ($\tau_R \ll 1$ fm/c) and naturally contribute to enhance the contribution of $\langle \mathcal{L}_0 \rangle$ which remains of comparable magnitude as $\langle \mathcal{L}_2 \rangle$ during the first few fm/c. In contrast, binding effects induced by $\mathcal{L}_1$ are completely negligible. The physical picture is clearly that of a pair with negligible internal interaction, whose constituents separate rapidly and are driven toward thermal equilibrium by stochastic forces.

Note that in all the plots, but this is more obvious at the highest temperature, the magnitudes of $\langle \mathcal{L}_2 \rangle$ and $\langle \mathcal{L}_3 \rangle$ reach a common value. This reflects the effect of friction (encoded in $\mathcal{L}_3$) in taming the effect of $\mathcal{L}_2$, eventually leading to equilibrium with the heat bath where $\langle \mathcal{L}_2 \rangle \approx \langle \mathcal{L}_3 \rangle$. As can be inferred from the discussion in appendix F, $\mathcal{L}_2[\mathcal{D}(s, s', t)]$ generically decreases with time, as $\exp(-\kappa(s-s')^2t)$, with $\kappa$ a constant. This particular
dependence on $s - s'$ is responsible for collisional decoherence, whose effect on the density matrix $\langle s|D|s' \rangle$ is to decrease the range of the off-diagonal elements from a typical bound state value, $|s - s'| \approx \frac{1}{M\alpha_S}$, to $|s - s'| \approx 0$ where the density matrix reaches a diagonal form, characteristic of a classical behavior. The superoperator $\mathcal{L}_3$ prevents this limit to be reached, restricting instead the off-diagonal elements of $\langle s|D|s' \rangle$ to remain non negligible over a distance of the order of the thermal wavelength when equilibrium is reached and $\langle \mathcal{L}_2 \rangle \approx \langle \mathcal{L}_3 \rangle$ (see figure 15, and also figure 22).

As for $\langle \mathcal{L}_1 \rangle$, its general behavior is easy to understand. Under the effect of collisions, the quark and the antiquark gain kinetic energy ($\mathcal{L}_2$ term). Then, due to the action of $\mathcal{L}_0$, they move away from the screened range of the potential and are therefore less affected by the attractive forces present in the singlet channel; besides, the color equilibration resulting from these collisions populate the octet channel and convert the attractive singlet potential into a slightly repulsive potential. Finally, one may also observe that the shrinking of the $|s - s'|$ range of the off-diagonal elements of the density matrix, caused by $\mathcal{L}_2$, contributes to an additional weakening of $\langle \mathcal{L}_1 \rangle \propto \langle V(s) - V(s') \rangle$.

One sees also in figure 14 that, after a short transient time where kinetic energy increases because of collisions with the heat bath particles, $\langle \mathcal{L}_0 \rangle$ decreases regularly with time. Since $\mathcal{L}_0 \propto y\partial_r$, with $y = s - s'$ and $r = (s + s')/2$, this regular decrease can be associated to that of the average value of $y$, while the gradient $\partial_r$ remains small and also decreases as the distance between the c and the $\bar{c}$ increases. Finally, figure 14 suggests that the effect of $\langle \mathcal{L}_4 \rangle$ is marginal. A specific analysis of the role of this operator is presented later in this subsection.

### 4.1.2 More on $\mathcal{L}_0 + \mathcal{L}_1$, $\mathcal{L}_2$ and $\mathcal{L}_3$

Figure 15 illustrates the impact of the superoperators $\mathcal{L}_3$ and $\mathcal{L}_4$ on the evolution of the singlet density operator $D_s$. With $\mathcal{L}_2$ alone (in addition to $\mathcal{L}_0 + \mathcal{L}_1$), the density matrix becomes quickly diagonal in coordinate-space, and by the end of the evolution, corresponding to $t = 10$ fm/c, most of the central peak reflecting the initial presence of a bound state is lost. The dynamics is dominated by collisions without dissipation, leading to collisional decoherence which makes the density matrix diagonal in coordinate space. As already discussed, the effect of $\mathcal{L}_3$ is to tame the effect of $\mathcal{L}_2$, and one can see that indeed in figure 15: at 10 fm/c the “diagonal” of the density matrix keeps a width of the order of the thermal wavelength ($\lambda_{th} \approx 0.3$ fm), and the correlation at small $r$ remains clearly visible.

The role of $\mathcal{L}_2$ is further illustrated in figure 16, which displays the singlet momentum space density obtained at $t = 20$ fm/c. It can be seen that the density obtained with $\mathcal{L}_2$ alone and with the combination $\mathcal{L}_0 + \mathcal{L}_1 + \mathcal{L}_2$ are pretty similar and both deviate significantly from the equilibrium Maxwell-Boltzmann distribution, already at small momenta: overall the distribution is much harder than the equilibrium one, a situation sometimes qualified of “overheating”. The full evolution yields a distribution which agrees with the thermal equilibrium distribution up to momenta of order $\sqrt{MT}$, in agreement with the results displayed in figure 5 (for the initial in-medium 1S state). This agreement with the equilibrium distribution is to a large extent due to the superoperator $\mathcal{L}_3$, which brings in friction, and slows down the diffusive evolution caused by $\mathcal{L}_0 + \mathcal{L}_1 + \mathcal{L}_2$ (see also [35, 36]). Finally, adding $\mathcal{L}_4$ does not change the picture in any visible way, except at large momenta where the presence of $\mathcal{L}_4$ is needed to ensure that the distribution remains positive.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Time evolution of the singlet density matrix $D_s(s', s)$, for three different cases: including superoperators up to $\mathcal{L}_2$ (top panels), up to $\mathcal{L}_3$ (middle panels) and up to $\mathcal{L}_4$ (bottom panels). From left panel to right panel: $t = 0.1$, 5 and 10 fm/c. For better readability, the color scale is different for each time. The initial state is a $c\bar{c}$ bound vacuum 1S singlet state, and the medium temperature is $T = 300$ MeV.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Momentum space singlet density obtained after evolving the vacuum 1S state up to 20 fm/c, in a medium at temperature $T = 300$ MeV. The impact of the superoperators $\mathcal{L}_3$ and $\mathcal{L}_4$ is illustrated by implementing various truncations, as indicated on the figure. The red line represents the evolution obtained with the full superoperator $\mathcal{L} = \mathcal{L}_0 + \mathcal{L}_1 + \mathcal{L}_2 + \mathcal{L}_3 + \mathcal{L}_4$. The dashed green curve shows the result of the evolution with $\mathcal{L}_2$ alone ($\mathcal{L} \mapsto \mathcal{L}_2$). The black dotted curve represents the Maxwell-Boltzmann distribution $\exp(-\frac{p^2}{MT})$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Left: time-evolution of rms(p) of the singlet component for various truncated QME as well as for the full QME (same conventions as for figure 16). Right: time evolution of the rms radius of the singlet spatial density for the same situation and the same selections of superoperators as in figure 16. The saturation for $t \gtrsim 10$ fm/c is due to the finite size of the grid.</figcaption></figure>

The impact of the various superoperators is further analyzed on figure 17, illustrating the evolution of $\sqrt{\langle p^2 \rangle}$ and $\sqrt{\langle r^2 \rangle}$. Starting from the 1S vacuum state, the unitary part $\mathcal{L}_0 + \mathcal{L}_1$ does not modify significantly the momentum distribution and leads to a moderate spatial expansion. The superoperator $\mathcal{L}_2$ steadily converts heat bath energy into $c\bar{c}$ kinetic energy but locally conserves the probability. In fact the evolution with $\mathcal{L}_2$ alone leads to an early-time reduction of $\langle r^2 \rangle$ (see eq. (F.11) and the corresponding discussion in appendix F). When considering $\mathcal{L}_0$ on top of $\mathcal{L}_2$, this effect is counterbalanced by the enhanced spatial diffusion triggered by collisions. The combined action of $\mathcal{L}_0$, $\mathcal{L}_1$ and $\mathcal{L}_2$ hence eventually yields a fast enhancement of both $\sqrt{\langle p^2 \rangle}$ and $\sqrt{\langle r^2 \rangle}$, which is tamed by the friction brought in by $\mathcal{L}_3$.

### 4.1.3 $\mathcal{L}_4$ and the positivity issue

The role of the superoperator $\mathcal{L}_4$ in insuring the positivity of the evolution equation for the density matrix is discussed in appendix A (see also appendix E and in particular the left panel of figure 32). In this section, we present a more direct analysis, by considering the spectrum of eigenvalues of the density matrix $\langle s|D|s' \rangle$, after an evolution of 5 fm/c starting from the 1S vacuum state, in a medium of temperature $T = 200$ MeV. The low value of the temperature is chosen on purpose, in order to maximize the possible effect of the $\mathcal{L}_4$ superoperator. At $t = 0$ fm/c, all eigenvalues vanish, except the one corresponding to the initial pure 1S state. As can be seen in figure 18, the evolution without the $\mathcal{L}_4$ term generates negative eigenvalues, a clear signal that the evolution does not preserve the positivity of the density matrix. Note however, that the negative eigenvalues are small, of the order of $10^{-4}$ at most in absolute magnitude, and their accumulated effect on the bulk of the dynamics is not significant, as we have seen before in several occasions. After adding the $\mathcal{L}_4$ term, one finds that the distribution of eigenvalues is shifted so that most eigenvalues end up positive, except for a very limited set of negative ones, whose magnitude $\sim 10^{-12}$ make them indistinguishable from the numerical noise affecting the evolution. Focusing on the distribution of the largest positive eigenvalues one finds that the impact of $\mathcal{L}_4$ remains quite moderate. Note finally that the impact of the negative eigenvalues is expected to be maximum at the low temperature chosen for the present analysis, which is at the boarder
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Distribution of eigenvalues of the singlet density matrix $D_s$ after 5 fm/c starting from the 1S vacuum state. The orange (resp. blue) — light gray (resp. mid-gray) in B&W printing — histogram represents the distribution obtained without (resp. with) the $\mathcal{L}_4$ superoperator. The temperature of the medium is $T = 200$ MeV.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Time evolution of the probabilities of the first three vacuum eigenstates (from top to bottom: $p_{1S}$, $p_{2S}$ and $p_{1P}$). These are obtained by evolving the density matrix first with $\mathcal{L}_0 + \mathcal{L}_1$ only (dot-dashed lines), then adding $\mathcal{L}_2$ (solid lines), $\mathcal{L}_3$ (long-dashed lines) and finally $\mathcal{L}_4$ (thick lines); the 1S evolution under $\mathcal{L}_2$ alone is also represented, by a dashed orange line. The medium temperature is $T = 300$ MeV, and the initial state is the vacuum 1S state. The horizontal dashed lines represent the asymptotic values of the survival probabilities.</figcaption></figure>

line of the validity of the master equations. At higher temperatures the impact of $\mathcal{L}_4$ on the evolution is found to be quite small, and in most cases negligible.

### 4.1.4 Probabilities and preferred basis

We turn now to a different way to analyze the results of the calculations, namely an analysis based on the probabilities $p_n = \langle \Phi_n | D_s | \Phi_n \rangle$ where $|\Phi_n\rangle$ is one of the vacuum eigenstates (see eq. (3.3)). In figure 19, we display the time evolution of these probabilities and, as we did before for other quantities, we consider the impact of the various superoperators on
<figure><img src="image_23.png" /><figcaption>Figure 20. The matrix element $\langle 2S | D_s | 1S \rangle$ as a function of time for various initial conditions.</figcaption></figure>

the evolution. When only the unitary superoperators $\mathcal{L}_0$ and $\mathcal{L}_1$ are involved, the system exhibits some kind of Rabi oscillations, mostly visible in the plot corresponding to $p_{2S}$. These oscillations find their origin in the fact that the initial state, the vacuum 1S state, is not an eigenstate of the in-medium Hamiltonian. This is also, partially, the cause of the rapid rise of $p_{2S}$ at early time. When $\mathcal{L}_2$ is added, collisions cause a rapid decrease of the various probabilities before these equilibrate to a common value for $t \gtrsim 20$ fm/c. If one tries to interpret the corresponding steady state as a thermal equilibrium state, and extract an effective temperature from the momentum density $\rho_s(p)$ at $t = 20$ fm/c, one finds a huge temperature, $T \approx 1.9$ GeV (hence the terminology of “overheating” alluded to earlier). When the $\mathcal{L}_3$ friction term is included, the decrease saturates earlier to larger values, as the c$\bar{c}$ system equilibrates at the temperature of the system. Again, one can note that the addition of the $\mathcal{L}_4$ terms leads only to tiny differences when $T > 200$ MeV, which confirms that those terms are, for practical purpose, unimportant when $T \gtrsim T_{QBM}$. As a final remark, consider the time evolution of the probability $p_{1S}$ under $\mathcal{L}_2$ alone. This is represented in figure 19 by the orange-dashed curve, where the evaluation is made with the same imaginary potential $W$ as for other curves. While it coincides with the full evolution for $t \lesssim 1$ fm/c, it deviates at later time. We note in particular that combining $\mathcal{L}_2$ with $\mathcal{L}_0 + \mathcal{L}_1$ yields a much more rapid suppression of $p_{1S}$. This can be attributed to the fact that, thanks to the operator $\mathcal{L}_0$, the c$\bar{c}$ pair, while gaining kinetic energy through the collisions with the plasma constituents can also expand spatially, as already observed on figure 17, resulting in a substantial suppression of the 1S probability at intermediate time.<fnref n="15" />

The analysis of the results in terms of eigenstates of the Hamiltonian complements that carried out earlier in terms of the densities, or more generally in terms of the coordinate space representation of the density matrix. On may argue that the coordinate space representation plays here the role of a “preferred basis” in the sense where this notion is used in the field of open quantum systems (see e.g. [47]). It is in this basis that we see most clearly the effect of decoherence induced by the collisions and leading the density matrix to a diagonal form. The

<footnote n="15">¹⁵In appendix F, the evolution under $\mathcal{L}_2$ alone is studied. In the absence of $\mathcal{L}_0$, diffusion occurs only in momentum space, and leads to broadening of the momentum distribution.</footnote>
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Orange (resp. red) solid line: initial time derivative of the probability to observe a vacuum-1S (resp. 2S) from a initial Gaussian state (3.4) under the evolution of $\mathcal{L}_2$ alone, according to eq. (F.33); the $\kappa$ value is tuned to reproduce $\Gamma_{1S}^{\text{vac.}} = 45$ MeV at $T = 300$ MeV. The 3 vertical lines correspond to the 3 initial states studied in section 3: vacuum 1S ($\sigma = \sigma_{1S} \equiv 0.387$ fm), in-medium 1S state ($\sigma \equiv 0.55$ fm), and compact Gaussian state ($\sigma = 0.165$ fm), while the square symbols represent the results from the full QME at $T = 300$ MeV for these IS (plus one extra simulation performed for $\sigma/\sigma_{1S} = 0.75$). The dashed lines represent the usual decay law $\dot{p}_{nS}(\sigma, 0) = -\Gamma_{nS} p_{nS}(\sigma; 0)$, thus simply reflecting the variation of $p_{nS}(\sigma; 0)$ with $\sigma$.</figcaption></figure>

representation in terms of “states”, meant to be localized states representing bound states, hides the fact that a dominant aspect of the physics is precisely the dissolution of bound states into continuum states: as can be read on figure 19, at $t = 10$ fm/c, $p_{1S} \lesssim 5\%$ and $p_{2S} \simeq 1\%$, meaning that most of the “states” at that time are continuum states. These represent weakly interacting c and $\bar{c}$, and are clearly better described in the coordinate representation. There is another issue. In the coordinate representation, decoherence occurs quickly. This is not the case, for instance, in the basis of vacuum eigenstates. To illustrate the point, consider figure 20, where the matrix element $\langle 2S|\mathcal{D}_s|1S\rangle$ between the first two vacuum-eigenstates is displayed as a function of time for $T = 300$ MeV for the three initial conditions considered in section 3. By looking at figure 20, one sees that this non-diagonal matrix element evolves on appreciable time scales (and even develops when not present in the initial condition). At late time, this non-diagonal element however vanishes. This reflects the fact that the density matrix approaches a unit matrix at late time, as we have seen already.

A further illustration of the interplay between quantum mechanical effects and the collisions is provided by the evolution of the 1S and 2S vacuum states, starting from an initial Gaussian wave packet with parameter $\sigma$ (see eq. (3.4) and appendix F), and a medium temperature $T = 300$ MeV. In figure 21, we compare the time derivatives of $p_{1S}$ and $p_{2S}$ at $t = 0$ obtained by solving the full QME with the values obtained from a calculation with $\mathcal{L}_2$ alone. The values $\sigma = 0.55$ fm, 0.387 fm and 0.165 fm are representative of the situations described respectively in figure 6 (right panel) for an initial in-medium 1S, figure 9 (upper right panel) for the vacuum 1S and figure 12 (right panel) for a compact singlet state. Also shown in figure 21 are the curves corresponding to simple decay laws of the individual states, $\dot{p}_{nS}(\sigma, 0) = -\Gamma_{nS} p_{nS}(\sigma; 0)$ (dashed lines). If the density matrix was of the form
$D^\sigma(t) = p_{1S}(t)|1S\rangle\langle1S| + p_{2S}(t)|2S\rangle\langle2S|$, the probabilities $p_{1S}$ and $p_{2S}$ would decay according to these laws. Figure 21 shows that this is clearly not the case. Even in cases where the initial density matrix is nearly in this diagonal form, non-diagonal elements develop in time and alter the evolution of the probabilities in a non trivial way. This is discussed in details in appendix F in the case where the evolution is driven by $L_2$ alone. The comparison displayed in figure 21 indicates that this approximate evolution with $L_2$ alone provides an accurate representation when $\sigma$ is not too different from $\sigma_{1S}$. The effect of coherence of the initial Gaussian wave packet is particularly visible for small values of $\sigma$. For $\sigma \lesssim 0.5$, one sees very little evolution, in particular for the 2S. This is because the small size of the (coherent) wave packet hinders the effect of the dipolar transitions. Also, the 2S exhibits a non-uniform pattern, increasing for $\sigma \gtrsim 0.5$ and decreasing for $\sigma > \sigma_{1S}$. Such a non-uniform behavior can be attributed to the destruction of the coherence in the initial overlap by collisions with the medium constituents. More details are given in appendix F.

## 4.2 Asymptotic distributions and steady states

We have seen in section 3 that, at late time, the $c\bar{c}$ density matrix appears to reach a steady state, independent of the initial conditions. This steady state presents features of a thermal equilibrium state,<fnref n="16" /> with however some (small) deviations (see figure 5). The goal of this section is to provide some understanding of these deviations and, more broadly, to better characterize the steady state. To do so, we shall examine the fixed point solutions of the QME, that is the solutions $D_{\text{asymp}}$ of the equation

$$ \mathcal{L} D_{\text{asymp}} = 0. \qquad (4.2) $$

It will be interesting to contrast the solution of the full QCD equations with that of their abelian version, which we shall refer to (somewhat abusively) to QED, and which is obtained by setting $D_s = D_o$ in the QCD equation for $D_s$, and take $D_s$ to be proportional to the abelian density matrix, which we shall denote $D_{\text{QED}}$. In order to ease the comparison of the numerical results, we keep the color factor $C_F$ in the expression of the singlet potential when solving the equation for $D_{\text{QED}}$. Note finally that the normalization $\text{tr}D_s + (N_c^2 - 1)\text{tr}D_o = 1$ entails, when $D_s = D_o$, that $\text{tr}D_s = 1/N_c^2$. As $D_{\text{QED}}$ is normalized so that $\text{tr}D_{\text{QED}} = 1$, one should therefore confront $D_{\text{QED}}/N_c^2$ and $D_s$ when comparing QCD and QED.

### 4.2.1 Generic structure of $D_{\text{asymp}}$

A generic feature of the asymptotic solution, common to QED and QCD, is illustrated in figure 22. The nearly diagonal structure of the density matrix is the result of the collisions of the $c\bar{c}$ pair with the surrounding medium constituents. We have associated this phenomenon to collisional decoherence and it is further discussed in appendix F where the role of the operator $L_2$ is thoroughly discussed, as well as the important role of the coordinate space basis as the preferred basis in which decoherence is the fastest. The final width is eventually controlled by the thermal wavelength (in both cases of QCD and QED), as indicated in the figure.

<footnote n="16">¹⁶To our knowledge, there is no guarantee that the late time behavior of a system described by a Lindblad equation in the QBM regime always corresponds to a thermal equilibrium distribution expressed in terms of the energy levels like the canonical Gibbs-Boltzmann distribution. In fact, this subject is still debated nowadays (see e.g. [48] and refs. therein). See also [49, 50].</footnote>
<figure><img src="image_25.png" /><figcaption>Figure 22. Left: contour plots of $D_{QED}(s, s')/N_c^2$ (levels: 0.001, 0.002, ..., 0.005, 0.01, 0.02, ..., 0.05 fm$^{-1}$); the skewed blue line across the main diagonal represents $2\lambda_{th}$ for $T = 200$ MeV. Right: same for $D_s(s, s')$ (levels: 0.001, 0.002, ..., 0.005, 0.01, 0.015 fm$^{-1}$).</figcaption></figure>

The peak seen near the origin is stronger in the case of QED than in the case of QCD. It reflects the correlations induced by the singlet attractive potential, which leads to an increase of the density near the origin. To quantify this effect, one can consider the Wigner transform of the density matrix

$$W_{s/QED}(r, p) = \int dy e^{-ipy} D_{s/QED}(r + \frac{y}{2}, r - \frac{y}{2}). \qquad (4.3)$$

A semi-classical approximation (see e.g. [14]) can provide useful guidance at sufficiently high temperature. This approximation, in its leading order, yields the following simple form for the Wigner transform for the asymptotic distribution

$$W_{sc}(r, p) \propto \exp\left(-\frac{p^2/M + V(T, r)}{T}\right). \qquad (4.4)$$

The corresponding density, $\rho_{sc}(r) \propto \exp(-V(T, r)/T)$, is shown in figure 23 for the two temperatures, $T = 200$ MeV and $T = 300$ MeV. In the abelian case, the semi-classical expression (4.4) indeed compares well with the exact asymptotic density, except at very small $r$.<fnref n="17" /> At $T = 200$ MeV, there is a significant difference between the abelian and the non-abelian cases. In the QCD case, the singlet to octet transitions, that are particularly efficient at large $r$, bring the pedestal of continuum states near its equilibrium value ($\simeq (1/9)(1/20)$ fm$^{-1}$). In the QED case, the pedestal is smaller. This can be attributed to the attraction of the potential which maintains a significant fraction of pairs concentrated at small distance. The picture changes somewhat at $T = 300$ MeV. At this temperature, collisions are sufficiently strong in QED to balance the effect of the attractive potential, and the pedestal is equilibrated in both cases. Correlatively the magnitude of the peak is significantly reduced.

<footnote n="17">¹⁷The agreement with the semi-classical expressions holds provided the potential $V$ is bounded from below, which is the case for the one-dimensional potential used in our calculations (see left panel of figure 30).</footnote>
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Left: asymptotic density associated to the steady state of the QM equations for $T = 200$ MeV: $D_s$ (solid curve), $D_{QED}$ (dashed blue), and semi-classical density $\propto \exp(-V(T, r)/T)$ (dot-dashed black), properly normalized. Right: the same for $T = 300$ MeV.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 24.** The density $\rho(r)$ for various sizes of the box ($L = 5$ fm, 10 fm and 20 fm); QED case. For positive values of $r$, a multiplicative factor $L/20$ has been applied to illustrate the scaling of the central peak with the pedestal.</figcaption></figure>

The density in the steady state depends on the size $L$ of the box where the system is confined. One expects in particular that in a box of infinite length, the weight of the bound state vanishes. The dependence on $L$ of the density $\rho(r)$ is illustrated in figure 24 in the abelian case. As $L$ increases, the shape of the peak remains essentially unaltered. This is because this peak is dominantly composed of localized bound states that do not depend too much on the box size. We may write $\rho(r) = c_1 e^{-V(r)/T}$, with $c_1$ a normalization constant. That is, the density is essentially constant outside the region of the peak (assuming that there $V(r) = 0$). The integral of the density is composed of two contributions: the peak contribution $\sim c_1 a_0 e^{V_0/T}$ — where $a_0$ is the width of the peak and $-V_0$ is some average value of the potential over the peak — and the continuum contribution $\sim c_1(L - a_0)$. Now, the normalization is kept fixed as the size of the box increases, $\int dr \rho(r) = 1 \simeq c_1 [L - a_0 + a_0 e^{V_0/T}]$. It follows that when $L$ is large enough ($L \gg a_0$), $c_1 \sim 1/L$ and the weight of the peak becomes indeed a negligible contribution to the total weight. Note however that the relative magnitude of the peak with respect to the pedestal is mostly determined by the potential $V(r)$ and is essentially independent of $L$. This is indeed illustrated in figure 24.
<figure><img src="image_28.png" /><figcaption>Figure 25. The Wigner distribution $W_s(r,p)$ for $r = 1$ fm and $T = 300$ MeV, showing the emergence of the Maxwell distribution, becoming gradually peaked at $p = 0$ as time increases: $t = 0$ fm/c (black), $1$ fm/c (purple), $2$ fm/c (blue), $5$ fm/c (orange), $10$ fm/c (green), and $20$ fm/c (red). For early times, the Wigner distribution presents oscillations reaching negative values that cannot be displayed in the semi-log plot.</figcaption></figure>

### 4.2.2 Asymptotic momentum distribution

At large relative distance, where the potential becomes essentially constant, the Wigner transform predicts, in the semi-classical approximation (4.4), that the momentum distribution is a Maxwell distribution. This is indeed the qualitative behavior that one observes for the full quantum calculation in figure 25, even for distances as small as $r = 1$ fm from the origin. The standard deviation $\sqrt{\langle p^2 \rangle_r - \langle p \rangle_r^2}$ is shown as a function of $r$ in the left panel of figure 26. One sees that the classical expectation $p_{\text{th}} \approx \sqrt{MT/2} \approx 0.5$ GeV/c is reached at late times, for $r \gtrsim 2$ fm. The relaxation is somewhat slower for the region $r \lesssim 2$ fm/c which remains under the influence of the potential for a longer time. The influence of the temperature on the asymptotic state is illustrated on the right panel of figure 26. Focusing on the large values of $r$ ($\gtrsim 2$ fm/c), one observes small but definite deviations with respect to the Maxwell-Boltzmann distribution, which do not disappear at the highest temperature considered. Their origin is clarified in appendix G where it is shown that the variance of $p$ at large relative distance is given by

$$ \text{var}(p) = \frac{MT}{2 + \gamma} \quad \text{with} \quad \gamma = \frac{\tilde{W}^{(4)}(0)}{16MT \tilde{W}''(0)}. \qquad (4.5) $$

The deviations from the expected value $MT/2$ are encoded in the coefficient $\gamma$ which involves 2nd and 4th derivatives of the imaginary potential at the origin.<fnref n="18" /> These derivatives depend on the prescription adopted to control the singular behavior of $W(r)$ near $r = 0$.<fnref n="19" /> For the choice made in appendix G, we have $\frac{\tilde{W}^{(4)}(0)}{\tilde{W}''(0)} \sim -\frac{3q_{\text{max}}^2}{20\ln(q_{\text{max}}/m_D)}$, with $q_{\text{max}} \gtrsim m_D$ (see

<footnote n="18">¹⁸This coefficient $\gamma$ should not be confused with the mass-shift introduced in some QME approaches.</footnote>

<footnote n="19">¹⁹For well-behaved $W$, one would expect $\frac{\tilde{W}^{(4)}(0)}{\tilde{W}''(0)} \sim m_D^2$, resulting in a subleading value of $\gamma$ in the QBM regime.</footnote>
<figure><img src="image_29.png" /><figcaption>**Figure 26.** Left: time evolution of the standard deviation of the momentum distribution as a function of $r = (s + s')/2$ for a fixed medium temperature $T = 300$ MeV. Right: same in the large time limit at for three different medium temperature $T = 200$ MeV ($t = 60$ fm/c), $300$ MeV ($t = 20$ fm/c) and $600$ MeV ($t = 20$ fm/c). The black dotted lines correspond to the classical result $\sqrt{MT}/2$ while the black dashed lines correspond to the corrected value $\sqrt{MT}/(2 + \gamma)$ (see text for details).</figcaption></figure>

eq. (G.5)).<fnref n="20" /> The dashed lines in figure 26 show that the resulting estimate of $\gamma$ is in good agreement with the numerical simulations.

Besides the deviations just noted for var($p$), one should also mention that the equilibrium distribution exhibits large deviations from the Maxwell-Boltzmann distribution when $p \gg \sqrt{MT}$, as seen e.g. on figure 16, and figure 28 below (see also the discussion after eq. (G.3) in appendix G). The asymptotic solution in the semi-classical approximation is, to within a normalizing factor, $D_s(y) = \tilde{y}^\eta K_\eta(\sqrt{\eta}\tilde{y})$, where $\eta = \frac{1}{2} + \frac{1}{|\eta|}$ and $\tilde{y} = \sqrt{MT}y$. The Wigner transform is $\propto (p^2 + \eta)^{-(\frac{1}{2}+\eta)}$. When $\gamma \to 0$, one recovers the correct MB distribution, but for finite $\gamma$, the decay is algebraic, in line with that of the Wigner transform of the QME, as illustrated for instance in figure 25.

We thus conclude that the QME used in the present paper lead to steady state momentum distributions that deviate slightly from thermal distributions. Analogous deviations have been pointed out recently by several authors [51, 52], who have also suggested methods to cure them. We shall not explore these here, since these deviations remain quite small in the present setting, and they have no impact on the main physics results presented in this paper.

### 4.2.3 Analysis in terms of in-medium eigenstates

We now turn to the analysis of $D_{asymp}$ in terms of its projections $p_n = \langle n | D_{asymp} | n \rangle$ on the eigenstates $|n\rangle$ of the in-medium Hamiltonian. The spectrum of the first three eigenstates of the in-medium Hamiltonian are displayed in the left panel of figure 27 for three values of

<footnote n="20">²⁰The correction $\gamma$ depends on a choice of regularization and is cut-off dependent. In [10], a fixed cutoff with a large value, $q_{max} = 4$ GeV, is adopted. As a consequence, $\gamma$ decreases $\propto \frac{q_{max}^2}{MT \ln(q_{max}/CT)}$ and the Maxwell-Boltzmann distribution is restored at “large” temperatures $\sim \frac{q_{max}^2}{M}$. The drawback of such a prescription is that it leads to large deviations for $T \simeq T_{QBM}$, up to 50% for c-quarks (depending on the precise value of $m_D$). In the present work, the definition adopted for the regularized $\tilde{W}$ is equivalent to take $q_{max} = \sqrt{MT}$, which yields a milder temperature dependence.</footnote>
<figure><img src="image_30.png" /><figcaption>**Figure 27.** Left panels: eigenstates (the first 3 as well as the 10th one) of the unitary $c\bar{c}$ Hamiltonian (with the T-dependent potential); central and right panels: eigenstates of the asymptotic $D_{QED}$ density matrix in the abelian case (central) and of the asymptotic $D_s$ for QCD (right), for $T = 200$ Mev (top), $T = 300$ Mev (middle) and $T = 400$ Mev (bottom). The eigenstates of $D_{QED}$ and of $D_s$ are generically complex-valued. For the selected temperatures, eigenstates of $D_s$ have a small imaginary part which has not been represented. See text for comments.</figcaption></figure>

At $T = 200$ MeV, one can identify three localized states that correspond to the three bound states. At the higher temperatures, $T = 300, 400$ MeV, only the 1S state remains localized, and its size is seen to increase with temperature. The 2S and 1P states are completely delocalized already at $T = 300$ MeV. And so are the higher eigenstates, which are continuum states spread over all the available volume of the simulation box. For comparison, we have also plotted the 10th eigenstate which displays a regular plane-wave behavior, with a slight distortion near the origin which gradually disappears as the temperature increases.

The Hamiltonian eigenstates are compared to the eigenstates of the abelian (middle column) and non-abelian (right column) density matrices. A striking similarity with the Hamiltonian eigenstates can be observed in the abelian case, suggesting that, in this case, the asymptotic density matrix is diagonal in the Hamiltonian eigenbasis. This is not so in the QCD case, as can already be observed at $T = 200$ MeV: while the 1S looks quite similar to the 1S eigenstate of the Hamiltonian, the 2S and 1P states are not localized states, but are completely delocalized; they are best viewed as continuum states. We note also that the spreading of the 1S state is much more rapid with increasing temperature than in the abelian case, while at $T = 400$ MeV, an apparent shift of the 1S and the 2S weights to higher values of $s$ is observed.

These features have a strong impact on the interpretation, in terms of Boltzmann weights, of the probabilities $p_n$ to find a given in-medium eigenstate in the asymptotic density
<figure><img src="image_31.png" /><figcaption>Figure 28. Probability to find in-medium eigenstates of the screened Hamiltonian in the asymptotic density operator $\mathcal{D}$ for either the abelian case (left) or QCD (right), for both $T = 200$ MeV (top) and $T = 300$ MeV (bottom).</figcaption></figure>

matrix. If $\mathcal{D}_{\text{asymp}}$ was close to a thermal equilibrium distribution, one would expect $p_n$ to be proportional to the Boltzmann factor, $p_n \propto \exp(-E_n/T)/Z$. This is indeed what one observes in the abelian case, for all the states with $E \lesssim 1$ GeV including the bound states, as illustrated in the left panel of figure 28. The QCD case presents a similar pattern, except for the lowest eigenstates at $T = 200$ MeV which significantly depart from the general trend. A similar behavior can be observed in the results of [35, 36], although a detailed comparison is difficult, given the differences in the potentials and in the numerical techniques employed to solve the QME. We note that deviations are observed only at the lowest temperature $T = 200$ MeV, and the 1S bound state exhibits thermal behavior at $T = 300$ MeV. We do not have a clear interpretation of this observation. Clearly $T = 200$ MeV is at the borderline of the QBM regime, and this could be the signal that some physical effects are missed by the present QME [19, 40]. However, one may also note that at $T = 200$ MeV, the eigenstates of $H_{Q\bar{Q}}$ differ significantly for those of the density matrix $\mathcal{D}_s$, in contrast to the QED case (see figure 28). Clearly this issue deserves further investigation.

The variation with temperature of the in-medium probabilities $p_n$ is illustrated in figure 29, for both the QCD and the abelian cases. Also shown in this figure is the analytic continuation of the probabilities (dashed lines) when these states cease to be bound. One expects that when the temperature is high enough all states are fully delocalized in the box. Then a simple calculation allows us to get the value of $p_n$. We have indeed

$$p_n = \int ds \, s' \psi_n^*(s) \mathcal{D}_s(s, s') \psi_n(s'), \quad (4.6)$$
<figure><img src="image_32.png" /><figcaption>**Figure 29.** Left: the probabilities $p_{1S}$ (orange), $p_{1P}$ (green) and $p_{2S}$ (red) to observe asymptotically the $c\bar{c}$ pair in the corresponding in-medium singlet state, as a function of the medium temperature $T$ (solid lines); for each state, the transition from a solid to a dashed curve corresponds the melting temperature. The (black) dashed-dotted line is the approximation (4.9). Right: the same in the abelian approximation.</figcaption></figure>

where $\psi_n(s)$ is the wave function of one of the delocalized state. When the temperature is large enough, we expect the density matrix to become diagonal quickly, to within the thermal wavelength $\lambda_{th}$ which is much smaller than the typical range of variation of $\psi_n(s)$, which, for the states considered, is of the order of a fraction of the size of the box, i.e. $\lambda_{th} \ll L$. For instance at $T = 300$ MeV, $\lambda_{th} \simeq 0.3$ fm, while $L = 20$ fm. Equivalently the typical momenta that control the dependence of $D(s, s')$ on $|s - s'|$ are large compared to those in $\psi_n$ (the typical momentum components of $\psi_n$ are of the order of a few times $\pi/L$). Under these conditions, we can approximate

$$ p_n = \int dr dy \psi_n^*(r + y/2) \psi_n(r - y/2) D_s(r, y) \\ \approx \int dr dy \psi_n^*(r) \psi_n(r) D_s(r, y). \qquad (4.7) $$

For $D_s(r, y)$ we can use the following factorized expression, which reproduces the Maxwell-Boltzmann momentum distribution with $\langle p^2 \rangle = MT/2$,

$$ D_s(r, y) = \rho_s(r) e^{-\frac{y^2}{4\lambda_{\text{th}}^2}} \qquad (4.8) $$

Noticing that in the considered regime the density is approximately constant $\rho_s(r) \simeq (1/N_c^2)(1/L)$, and assuming $\phi_n(r)$ to be normalized to unity, one gets

$$ p_n \approx \frac{2\sqrt{\pi} \lambda_{\text{th}}}{N_c^2 L}. \qquad (4.9) $$

As can be seen in figure 29, this simple formula accounts quite well for the values of the probabilities at temperatures $T \gtrsim 300$ MeV $\approx T_{\text{damp}}$. In fact the same agreement is found for the abelian case, but at higher temperature. This may be attributed to the fact that, in QCD, color singlet-octet transitions accelerate the delocalization of states, as is visible for instance in figure 27.
## 5 Conclusion

We have solved numerically the non-abelian quantum master equations derived in ref. [14] in the quantum Brownian regime. These can be put in the form of a Lindblad equation. This was used to study the dynamics of a heavy charm-anticharm pair in a Quark-Gluon Plasma in thermal equilibrium. The operators entering these equations were simplified in order to cope with various artificial divergences induced by successive derivatives of the complex potential. The resulting “minimal” set of operators conserves the basic property of a Lindblad equation, in particular the positivity of the evolution. The specific role of the operator $L_4$ in guaranteeing this property was emphasized. The equations were solved in a one-dimensional setting, exploiting the complex potential specifically developed for one-dimensional studies in [39]. This potential is adjusted so as to yield numbers for physical quantities that are comparable to those of a three dimensional potential. This allowed us, along the paper, to make remarks of phenomenological interest, although a detailed phenomenological study is beyond the scope of this paper.

The evolution of the reduced density matrix of the $c\bar{c}$ pair was studied for different initial conditions of the pair (size and color), and various temperatures of the quark-gluon plasma. The case of an expanding plasma was also briefly analyzed. A generic feature of the evolution is the fact that collisions lead to a rapid squeezing of the coordinate space density matrix along its diagonal, a characteristic effect of decoherence making part of the semi-classical dynamics. It also underlines the role of coordinate space as a “preferred basis” for the analysis of the dynamics. The color dynamics, captured by coupled equations that connect the singlet and octet sectors of the reduced density matrix, involve a subtle coupling between coordinate space and color: this is because of the coordinate space dependence of the color-dipolar interactions that control the overall dynamics. Also, the continuum states play an important role. They are quickly populated, and it is through them that color equilibrates the fastest. In contrast, since dipolar interactions vanish at vanishing relative distances, color equilibrates on a much longer time scale in the region occupied by bound states, which enhances the survival probability of the bound states. Two subregimes were identified in the quantum Brownian regime. For $T \gtrsim T_{\text{damp}}$ (overdamped subregime), one observes the dominance of the collisional decoherence over the binding effects, while for $T \lesssim T_{\text{damp}}$ the vacuum bound states gain in significance as the temperature decreases. Much of the qualitative features of the overdamped regime can be understood from the dominant role in this regime of a specific operator of the QME, the operator $L_2$.

We have shown that at late time the density matrix evolves toward a steady state that is close to a thermal equilibrium state if the temperature is high enough. Because of dissipation the typical range of the non-diagonal matrix elements of the reduced density matrix ends up being of the order of the thermal wavelength. The equilibrium state is that of a bound pair in equilibrium with continuum states that occupy the box where the simulation is done. The momentum distribution for the continuum states is close to a Maxwell-Boltzmann distribution with a small deviation that depends on the details of the complex potential. For the lowest temperature available in our calculations, one finds that the populations of the bound states deviate significantly from the Boltzmann factors. This is to be contrasted with the abelian situation where the populations of the various states are nearly perfectly
given by the Boltzmann weights. This feature, and more broadly the nature of the steady state, deserve further study.

The analysis in terms of projections on charmonium eigenstates as a function of time provided further insight. By projecting on vacuum eigenstates one gets access to the probability to find a particular charmonium state in the plasma as a function of time, a quantity which is often calculated in phenomenological studies. The projections onto the in-medium screened Hamiltonian provides a complementary view of the system properties. This is complicated by the fact that such bound states quickly disappear in favor of continuum states as the temperature increases. The 1S state, however, survives for a longer time, and its evolution dominates that of the density matrix. Thus the evolution of the latter can, to a large extent, be understood in terms of the evolution of the 1S state. In fact, the very first instants of the evolution depend much on the initial state, and can exhibit quantum mixing before such effects start to be hindered by collisions. The later stages of the evolution, on the other hand, have a universal character, driven essentially by the generic evolution of the coordinate space density matrix which, when analyzed in terms of projections on vacuum eigenstates, yields indeed a similar evolution of both the 1S and the 2S.

There are several aspects of the present work that are worth investigating further. One of them concerns the nature of the steady state, as mentioned above. Another one is related to the semi-classical dynamics and its interplay with color dynamics. As pointed out in [14] this interplay hinders the straightforward application of semi-classical approximations that have proven to be successful in the abelian case (see e.g. [10]). Having at our disposal complete solutions of the QME constitutes a benchmark that could be very helpful for deepening our understanding of this issue. At the lowest temperature that we have explored, we have seen that the eigenstates of the asymptotic QCD density matrix differ from those of the in-medium Hamiltonian, in contrast to what happens in the analogous abelian situation. This points possibly to missing ingredients in the QME, which would be needed for instance to cope with the energy gaps between various bound states when these bound states are playing a prominent role in the dynamics [15]. Ways to improve the present QME in order to address this issue are known and implementing them in the present context would represent an interesting extension of the present study. Finally the tools that have been developed for the present study may constitute a useful playground for a number of more phenomenologically oriented studies.

## Acknowledgments

We wish to thank Aoumeur Daddi Hammou and Joerg Aichelin for valuable discussions. The authors thank the Région Pays de la Loire and Subatech for support. S.D. is supported by the Centre national de la recherche scientifique (CNRS) and Région Pays de la Loire and acknowledges the support of Narodowe Centrum Nauki under grant no. 2019/34/E/ST2/00186. R.K. was under contract No. 2015-08473. PBG is supported by the European Union's Horizon 2020 research and innovation program under grant agreement No 824093 (STRONG-2020).
<table><caption>**Table 1.** Structure of the various superoperators $\mathcal{L}_{2b}$, $\mathcal{L}_{3b}$ and $\mathcal{L}_{4b}$ in terms of their spatial derivative content, generically written as $D_o$, $\partial D_o$ and $\partial^2 D_o$. $W^{(i)}$ denotes the $i$th derivative of $W$. The subscripts “cst”, “lin” and “quad” refer to the number of derivatives acting on $D_o$, respectively 0, 1, and 2.</caption><thead><tr><th></th><th>D<sub>o</sub></th><th>∂D<sub>o</sub></th><th>∂<sup>2</sup>D<sub>o</sub></th></tr></thead><tbody><tr><td>L<sub>2b</sub></td><td>L<sub>2b,cst</sub>: ×W</td><td>0</td><td>0</td></tr><tr><td>L<sub>3b</sub></td><td>L<sub>3b,cst</sub>: ×W<sup>(2)</sup></td><td>L<sub>3b,lin</sub>: ×W<sup>(1)</sup></td><td>0</td></tr><tr><td>L<sub>4b</sub></td><td>L<sub>4b,cst</sub>: ×W<sup>(4)</sup></td><td>L<sub>4b,lin</sub>: ×W<sup>(3)</sup></td><td>L<sub>4b,quad</sub>: ×W<sup>(2)</sup></td></tr></tbody></table>

## A Positivity and structure of the QME

In this appendix, we sketch the proof that, once the term $\mathcal{L}_4$ is included, the equations used in this paper preserve the positivity of the density matrix.<fnref n="21" /> To that aim, we first note that the Liouville superoperator $\mathcal{L}$ of eq. (2.4) can be written explicitly in terms of sub-operators that represent the various color transitions:

$$ \begin{aligned} \frac{d}{dt} D_s &= \mathcal{L}^{ss} D_s + \mathcal{L}^{so} D_o \\ \frac{d}{dt} D_o &= \mathcal{L}^{os} D_s + \mathcal{L}^{oo} D_o. \end{aligned} \qquad (A.1) $$

In turn, each of these sub-operators can be decomposed in terms of the $\mathcal{L}_i$'s defined in eq. (2.5). We shall be concerned in this appendix only with the $\mathcal{L}_i$'s which are proportional to the imaginary part $W$ of the complex potential and its spatial derivatives, namely $\mathcal{L}_2$, $\mathcal{L}_3$ and $\mathcal{L}_4$. These are the operators that are involved in the positivity constraint that we want to study ($\mathcal{L}_0$ and $\mathcal{L}_1$ are unitary operators which preserve positivity). Recall that the operators $\mathcal{L}_2$, $\mathcal{L}_3$ and $\mathcal{L}_4$ have ordered dependence on the mass of the heavy quark, $\mathcal{L}_2 \propto \mathcal{O}(M^0)$, $\mathcal{L}_3 \propto \mathcal{O}(M^{-1})$ and $\mathcal{L}_4 \propto \mathcal{O}(M^{-2})$. The same ordering holds for the color transition operators introduced above. As an illustration consider the operator $\mathcal{L}^{so}$ which realizes transitions from octet to singlet states. Using the convention of ref. [14], we write this operator as follows:

$$ \mathcal{L}^{\text{so}} = C_F (\mathcal{L}_{2b} + \mathcal{L}_{3b} + \mathcal{L}_{4b}), \qquad (A.2) $$

where the superoperators of type $\mathcal{L}_b$ correspond to one-gluon exchange diagrams in which the gluon propagator connects quark lines that are located on both sides of the Schwinger-Keldysch contour. Referring to eq. (2.5), and to the explicit expressions listed in appendix C, one sees that $\mathcal{L}_{2b}$ just involves $W$ while $\mathcal{L}_{3b}$ and $\mathcal{L}_{4b}$ contain both derivatives of $W$ and derivatives acting on $D_o$. The derivative content of the various operators is summarized in table 1.

For the physical interpretation of the various superoperators,<fnref n="22" /> we note that, in the large mass limit, the sole effect of the collisions is to shift the momentum of the heavy quarks and to rotate their colors. The dominant superoperator in this limit is $\mathcal{L}_2$. Finite-mass corrections, which involve in particular recoil corrections that account for the energy conservation in

<footnote n="21">²¹Although the proof holds for any dimension, we focus on the 1D case when it comes to specific expressions.</footnote>

<footnote n="22">²²See also [19], section 3.2.1.</footnote>
individual collisions, yield the superoperators $\mathcal{L}_3$ and $\mathcal{L}_4$. The superoperator $\mathcal{L}_3$ contributes a term $\mathcal{L}_{3b,cst}$ which has a structure similar to $\mathcal{L}_{2b,cst}$, namely it is independent of the heavy quark momenta, $\mathcal{L}_{3b,cst} \sim \frac{q^2}{8MT}$; there is also in $\mathcal{L}_3$ a contribution that depends explicitly on the heavy quark momentum $\vec{p}_Q$, $\mathcal{L}_{3b,lin} \sim -\frac{\vec{q}\vec{p}_Q}{4MT}$, which leads to dissipation. Similar considerations apply to the superoperator $\mathcal{L}_4$.

To proceed, it is convenient to express $W(r)$ in terms of the following spectral representation

$$W(r) = \int_{-\infty}^{+\infty} dq \rho(q)(1 - \cos qr) = \int_{-\infty}^{+\infty} dq \rho(q)(1 - e^{-iqr}), \quad (A.3)$$

where the spectral density $\rho(q)$ is an even function of $q$. Its explicit form is actually not needed for the argument of this appendix, aside from the fact that it is a positive definite function of the (one-dimensional) momentum $q$.<fnref n="23" />

The positivity property is in principle guaranteed by the Lindblad structure of the master equations (2.5). However, it is instructive to examine concretely how this is realized in the present case. To that aim, we shall look at the transition rates and verify that they are all positive definite. We shall do that by focusing on the set of superoperators (A.2) which describe the octet to singlet transitions. In the basis that, at time $t$, diagonalizes the density matrix, the octet and singlet density operators at time $t$ take the form $\mathcal{D}_o = \sum_l p_{o,l}|l\rangle\langle l|$, and $\mathcal{D}_s = \sum_m p_{s,m}|m\rangle\langle m|$. One can then show that the time variation of the singlet population $p_{s,n}(t) = \langle n|\mathcal{D}_s(t)|n\rangle$ in this specific basis is given by the following rate equation

$$\frac{dp_{s,n}}{dt} = \sum_{l} t_{nl} p_{o,l} - \Gamma_{n} p_{s,n} \quad \text{with} \quad \Gamma_{n} = \sum_{l} t_{ln}, \quad (A.4)$$

where the gain (resp. loss) term stems from the $\mathcal{L}^{\text{so}}$ (resp. $\mathcal{L}^{\text{ss}}$) superoperator.<fnref n="24" /> For the sake of illustration, we show how the transition coefficients $t_{nl}$ relate to $\mathcal{L}^{\text{so}}$ in the case of the superoperator $\mathcal{L}_2 = \mathcal{L}_{2b} = 2(W(\frac{s+s'}{2}) - W(\frac{s-s'}{2}))$ (see appendix C). Using the spectral decomposition (A.3) of $W$, one obtains

$$\mathcal{L}_{2b} = 2 \int dq \rho(q) \left( e^{-iq \frac{s-s'}{2}} - e^{-iq \frac{s+s'}{2}} \right) = 4 \int dq \rho(q) \sin \left( \frac{qs}{2} \right) \sin \left( \frac{qs'}{2} \right). \quad (A.5)$$

The contribution from $C_F \mathcal{L}_{2b} \mathcal{D}_o$ to $\dot{p}_{s,n}$ thus reads

$$\frac{dp_{s,n}}{dt} = C_F \langle n | \mathcal{L}_{2b} \mathcal{D}_o | n \rangle = C_F \sum_l \int dq \rho(q) \left| 2 \int ds \psi_{s,n}^*(s) \sin\left(\frac{qs}{2}\right) \psi_{o,l}(s) \right|^2 p_{o,l}. \quad (A.6)$$

On recognizes in eq. (A.6) the form of the gain term in eq. (A.4), provided the transition coefficient $t_{nl}$ is defined as $t_{nl} = C_F \int dq \rho(q)|T_{n,l}(q)|^2$, where $T_{n,l}(q)$ is the following transition amplitude

$$T_{n,l}(q) = 2 \int ds \psi_{s,n}^*(s) \sin\left(\frac{qs}{2}\right) \psi_{o,l}(s).$$

<footnote n="23">²³In the case of the HTL imaginary potential of eq. (B.1), $\rho(q) = \frac{\alpha_s T}{2m_D} \left[ \frac{\pi}{2} - \tan^{-1} \left( \frac{|q|}{m_D} \right) - \frac{m_D |q|}{m_D^2 + q^2} \right]$, which is indeed > 0.</footnote>

<footnote n="24">²⁴The detailed calculations leading to these results can be found in [53]. Let us mention that off-diagonal elements are generated from the dynamics, that is $\langle n|\frac{dD}{dt}|n' \neq n \rangle \neq 0$, but they only correspond to the “rotation” of the local basis $\{|n(t)\rangle\}$. In this local basis $\mathcal{D}$ stays diagonal by construction.</footnote>
For the full $\mathcal{L}_{\text{so}}$ superoperator, the expression of $t_{nl}$ reads

$$t_{nl} = C_F \int dq \rho(q) |(1 - (q\lambda_{\text{th}})^2/8)T_{n,l}(q) + T_{n,l}^{(3)}(q)|^2 \quad (A.7)$$

where $\lambda_{\text{th}} = 1/\sqrt{MT}$, while a new transition amplitude $T_{n,l}^{(3)}(q)$ needs to be introduced at the $\mathcal{L}_{3b}$ level:

$$T_{n,l}^{(3)}(q) = \frac{q\lambda_{\text{th}}^2}{2} \int ds \psi_{s,n}^{\star}(s) \cos\left(\frac{q}{2}s\right) \psi_{o,l}^{\prime}(s).$$

Clearly, the expression (A.7) is always positive iff $\rho(q) > 0$ for all $q$.

In eq. (A.7), six terms contribute to the integrand, which can be mapped to those enumerated in table 1. The leading fluctuations term for an infinite mass heavy quark (recoilless limit) corresponds to $|T_{n,l}(q)|^2$ originating from $\mathcal{L}_{2b}$, as shown above. The cross terms, $(T_{n,l}^{(3)}(q)T_{n,l}^{\star}(q) + \text{cc})$ and $-(q\lambda_{\text{th}})^2/4|T_{n,l}(q)|^2$, correspond respectively to $\mathcal{L}_{3b,\text{lin}}$ and $\mathcal{L}_{3b,\text{cst}}$. The last 3 terms $|T_{n,l}^{(3)}(q)|^2$, $-(q\lambda_{\text{th}})^2/8(T_{n,l}^{(3)}(q)T_{n,l}^{\star}(q) + \text{cc})$ and $(q\lambda_{\text{th}})^4/64|T_{n,l}(q)|^2$, correspond respectively to the higher order superoperators $\mathcal{L}_{4b,\text{quad}}$, $\mathcal{L}_{4b,\text{lin}}$, and $\mathcal{L}_{4b,\text{cst}}$.

A very similar decomposition applies for the $\mathcal{L}^{\text{os}}$, $\mathcal{L}^{\text{ss}}$, and $\mathcal{L}^{\text{oo}}$ blocks, including superoperators of type $\mathcal{L}_a$ (in particular $\mathcal{L}^{\text{ss}}$ will generate the loss terms for the singlet state probabilities in eq. (A.4)). This demonstrates the equivalence between the action of the superoperator $\mathcal{L}_2 + \mathcal{L}_3 + \mathcal{L}_4$ and the gain-loss kinetic rate equation (A.4) which preserves positivity.

## B Implementation of the higher-order operator and UV regularization

While the superoperator $\mathcal{L}_3$ involves at most second derivatives of $W(r)$ with respect to $r$, the superoperator $\mathcal{L}_4$ involves third and fourth order derivatives at finite values of $r$, as well as a fourth derivative at the origin for $\mathcal{L}_{4a}$. Because of the specific analytic form of $W(r)$, such derivatives entail possible UV divergences if the integrals entering the definition of $W$ are not properly regulated.

To appreciate the difficulty, let us recall that the one-dimensional imaginary potential that we use derives from the three-dimensional HTL expression given in eqs. (2.9) and (2.10), which we recall here for convenience:

$$W(r, T) = \alpha_S T \phi(m_D r), \quad (B.1)$$

where $m_D$ is the Debye mass, $r$ is the distance between the quark and the antiquark, and

$$\phi(\tilde{r}) = 2 \int_0^\infty dz \frac{z}{(z^2+1)^2} \left(1 - \frac{\sin(z\tilde{r})}{z\tilde{r}}\right), \quad (B.2)$$

with $\tilde{r} = m_D r$ and $z = q/m_D$. The derivatives of $\phi(\tilde{r})$ at the origin are divergent, starting from the second one which suffers from a logarithmic UV divergence.<fnref n="25" /> Since the range of $q$ values that lead to such divergences is in fact outside the domain of validity of the HTL approximation, a standard procedure consists in regulating the integrals by imposing some

<footnote n="25">²⁵At finite $r$ all derivatives are well defined, although starting from the fourth one, they cannot be evaluated by derivation under the integral sign.</footnote>
<table><caption>**Table 2.** Same as table 1 for the “minimal set” of superoperators associated to the regularization of the imaginary potential defined by eq. (B.3).</caption><thead><tr><th></th><th>D<sub>o</sub></th><th>&part;D<sub>o</sub></th><th>&part;<sup>2</sup>D<sub>o</sub></th></tr></thead><tbody><tr><td>L<sub>2b</sub></td><td>L<sub>2b,cst</sub>: &times; W&#771;</td><td>0</td><td>0</td></tr><tr><td>L<sub>3b</sub></td><td>0</td><td>L<sub>3b,lin</sub>: &times; W&#771;<sup>(1)</sup></td><td>0</td></tr><tr><td>L<sub>4b</sub></td><td>0</td><td>0</td><td>L<sub>4b,quad</sub>: &times; W&#771;<sup>(2)</sup></td></tr></tbody></table>

cutoff [13]. In our work, we have followed an alternate method that has the advantage of simplifying the structure of the equations. This method is based on the observation that the transition probability (A.7) suffers from the same UV divergences. To cure the divergences in eq. (A.7), one proceeds as follows: one assumes that higher contributions in the gradient expansion would generate a prefactor $(1 - \frac{(q\lambda_{\text{th}})^2}{8})$ in front of $T_{n,l}^{(3)}(q)$ and one then substitutes this common prefactor by $(1 + \frac{(q\lambda_{\text{th}})^2}{8})^{-1}$, so that

$$ \rho(q) \left| \left( 1 - \frac{(q\lambda_{\text{th}})^2}{8} \right) T_{n,l}(q) + T_{n,l}^{(3)}(q) \right|^2 \rightarrow \frac{\rho(q)}{\left( 1 + \frac{(q\lambda_{\text{th}})^2}{8} \right)^2} \left| T_{n,l}(q) + T_{n,l}^{(3)}(q) \right|^2, $$

which is UV-convergent and still positive. Up to $O((q\lambda_{\text{th}})^4)$, the new expression is similar to the former one, apart from different numerical coefficients for the terms $L_{4b,cst}$ and $L_{4b,lin}$ which are anyhow expected to be subdominant at high temperature. This procedure leads to a regularized potential $\tilde{W}$:

$$ \tilde{W}(r) = \int dq \frac{\rho(q)}{\left(1 + \frac{(q\lambda_{\text{th}})^2}{8}\right)^2} (1 - e^{-iqr}) = \int dq \tilde{\rho}(q) (1 - e^{-iqr}), \quad (B.3) $$

where the modified spectral density $\tilde{\rho}(q) = \rho(q)\left(1 + \frac{(q\lambda_{\text{th}})^2}{8}\right)^{-2}$ contains all the corrections discussed above. The corresponding regularized transition rates then read

$$ \tilde{t}_{nl} = C_F \int dq \tilde{\rho}(q) |T_{n,l}(q) + T_{n,l}^{(3)}(q)|^2. \quad (B.4) $$

Note that the substitution of the original $t_{nl}$ by the regularized $\tilde{t}_{nl}$ is equivalent to selecting a reduced number of terms in the original $L_b$ superoperator ($L_{2b}$, $L_{3b,lin}$ and $L_{3b,quad}$, together referred to the “minimal set” of superoperators warranting positivity), and expressing them in terms of $\tilde{W}$ (see table 2). Let us emphasize that the present construction of the minimal set allows to recover the original $L_{3b}$ term while displaying a lighter structure.<fnref n="26" /> We also note that the definition (B.3) amounts essentially to take as UV regulator $\Lambda = \lambda_{\text{th}}^{-1}$ in the original expression of the potential $W$ (see the discussion in appendix G).

<footnote n="26">²⁶The simplification brought by the use of the minimal set allows a physical illustration of the role of $L^{(4)}$. In a semi-classical expansion, it yields in particular a contribution $\propto \tilde{W}^4(0)y^2\partial^2 D$, or in terms of the Wigner transform $W(r = \frac{s+s'}{2}, p)$, a contribution $\propto \partial_p^2(p^2W)$. This represents a correction $\propto O(p^2)$ to the diffusion coefficient appearing in the Fokker-Planck equation obtained from $L_2$ and $L_3$ [14].</footnote>
Up to now, we have discussed the $\mathcal{L}_b$ type of superoperators, present in the $\mathcal{L}^{\text{so}}$ and $\mathcal{L}^{\text{os}}$ blocks, as well as in $\mathcal{L}^{\text{oo}}$. We now discuss the case of the $\mathcal{L}_a$ type of superoperators, present in the $\mathcal{L}^{\text{ss}}$ and $\mathcal{L}^{\text{oo}}$ blocks and obtained by considering exclusively diagrams for which the gluon propagator connects quark lines located on a single side of the Schwinger-Keldysch contour. The regularization discussed above, $W \to \tilde{W}$, also leads to a simplification of the $\mathcal{L}_a$ superoperators. The $\mathcal{L}_{2a}$ contribution to the ss block is readily defined as (see appendix C)

$$ \mathcal{L}_{2a}^{\text{ss}} = C_F(2\tilde{W}(0) - \tilde{W}(s) - \tilde{W}(s')). \qquad (B.5) $$

This leads to the following contribution to $\Gamma_n$, the dissociation of a singlet state $|n\rangle$:

$$ \Gamma_{n,2} = 2C_F \langle n | \tilde{W}(s) - \tilde{W}(0) | n \rangle, \quad \text{with} \quad \tilde{\Gamma}(s) = \tilde{W}(s) - \tilde{W}(0). \qquad (B.6) $$

In order to construct $\mathcal{L}_{3a}$ consistently within the minimal set prescription we start from the transition probabilities $t_{nl}$ defined in appendix A and rely on detailed balance ($t_{nl} = t_{nl}$)$^{27}$ to derive the decay rates $\Gamma_n$ governing the loss term in the rate equation (A.4). This allows us to calculate the various contributions of $\Gamma_n$ by selecting the corresponding terms in the expression of $t_{nl}$ and ultimately express these contributions as “$\mathcal{L}_a$-type” superoperators acting on the density operator. For instance, one explicitly recovers the expression (B.6) by just retaining the $|T_{nl}|^2$ contributions in $t_{nl}$. Following this method, one defines

$$ \Gamma_{n,3} = C_F \sum_l \int dq \, \tilde{\rho}(q) (T_{ln}^*(q) T_{ln}^{(3)}(q) + cc) \qquad (B.7) $$

where

$$ T_{ln}(q) = 2 \int ds' \psi_l^*(s') \sin\left(\frac{qs'}{2}\right) \psi_n(s') \quad \text{and} \quad T_{ln}^{(3)}(q) = \frac{\lambda_{\text{th}}^2}{2} \int ds \psi_l^*(s) \cos\left(\frac{qs}{2}\right) q\psi_n'(s). $$

Injecting these transition elements into the expression of $\Gamma_{n,3}$, one obtains a threefold integral. However using the summation on $l$ and the completeness relation leads to a factor $\delta(s - s')$ which is used to simplify the final result:

$$ \begin{aligned} \Gamma_{n,3} &= C_F \lambda_{\text{th}}^2 \int dq \, \tilde{\rho}(q) \int ds \, ds' \, \delta(s-s') \left[ \sin\left(\frac{qs'}{2}\right) \cos\left(\frac{qs}{2}\right) q\partial_s + \right. \\ & \qquad \left. \sin\left(\frac{qs}{2}\right) \cos\left(\frac{qs'}{2}\right) q\partial_{s'} \right] \psi_n^*(s') \psi_n(s) \\ &= \frac{C_F \lambda_{\text{th}}^2}{2} \int dq \int ds \, \tilde{\rho}(q) \sin(qs) q\partial_s |\psi_n(s)|^2. \end{aligned} \qquad (B.8) $$

Expressing $\sin(qs)q$ as $-\partial_s \cos(qs)$, one finds

$$ \begin{aligned} \Gamma_{n,3} &= -\frac{C_F \lambda_{\text{th}}^2}{2} \int ds \, \partial_s \left( \int dq \, \tilde{\rho}(q) \cos(qs) \right) \nabla |\psi_n(s)|^2 = \frac{C_F \lambda_{\text{th}}}{2} \int ds \, \tilde{W}'(s) \partial_s |\psi_n(s)|^2 \\ &= \frac{C_F \lambda_{\text{th}}}{2} \left( \int ds \, \psi_n^*(s) \tilde{W}'(s) \partial_s \psi_n(s) + cc \right), \end{aligned} \qquad (B.9) $$

<footnote>$^{27}$Because of the specific structure of the $\mathcal{L}_{3a}$ superoperator, the procedure put forward for the $\mathcal{L}_b$ superoperators can not be blindly applied. This is the reason for the strategy used here.</footnote>
corresponding to the following (1D$^{28}$) contribution to $\mathcal{L}_{3a}^{ss}$,$^{29}$

$$ \mathcal{L}_{3a}^{ss} = -\frac{C_F \lambda_{\text{th}}^2}{2} (\tilde{W}'(s) \partial_s + \tilde{W}'(s') \partial_{s'}) . \qquad (B.10) $$

Following the same method, one obtains

$$ \begin{aligned} \Gamma_{n,4} &= \frac{C_F \lambda_{\text{th}}^4}{4} \int dq \, \tilde{\rho}(q) \int ds \, ds' \, \delta(s-s') \cos\left(\frac{qs}{2}\right) \cos\left(\frac{qs'}{2}\right) q^2 \partial_s \psi_n(s) \partial_{s'} \psi_n^*(s') \\ &= \frac{C_F \lambda_{\text{th}}^4}{4} \int ds \int dq \, \tilde{\rho}(q) q^2 \cos^2\left(\frac{qs}{2}\right) \psi_n'(s) \psi_n'^*(s) \\ &= \frac{C_F \lambda_{\text{th}}^4}{4} \int ds \, (\tilde{W}''(s) + \tilde{W}''(0)) \psi_n'^*(s) \psi_n'(s), \end{aligned} \qquad (B.11) $$

while the associated $\mathcal{L}_{4a}^{ss}$ superoperator is$^{30}$

$$ \mathcal{L}_{\text{ss},4a} = \frac{C_F \lambda_{\text{th}}^4}{16} \left[ \tilde{W}''(0) (\partial_s \partial_s + s \leftrightarrow s') + (\partial_s (\tilde{W}''(s) \partial_s) + s \leftrightarrow s') \right]. \qquad (B.12) $$

### C Color transition operators

This appendix provides the complete list of the color transition superoperators $\mathcal{L}^{ss}$, $\mathcal{L}^{so}$, $\mathcal{L}^{os}$ and $\mathcal{L}^{oo}$ corresponding to the “minimal set” discussed in appendix B. Details on their derivation can be found in [53]. The superoperators below are written in the coordinate representation and are supposed to act on the color components $\langle s | \mathcal{D}_s | s' \rangle$ and $\langle s | \mathcal{D}_o | s' \rangle$ of the density matrix. The color factors appropriate to the various channels are properly included. Note that the non-diagonal superoperators $\mathcal{L}^{so}$ and $\mathcal{L}^{os}$ (superoperators of type $\mathcal{L}_b$ in the terminology of [14]) are proportional to the imaginary potential $W$ and its derivatives. The diagonal superoperators $\mathcal{L}^{ss}$ and $\mathcal{L}^{oo}$ (where $\mathcal{L}^{ss}$ is of type $\mathcal{L}_a$, while both $\mathcal{L}_a$ and $\mathcal{L}_b$ types are present in $\mathcal{L}^{oo}$) contain in addition kinetic energy contributions as well as the real part $V$ of the potential. In each superoperator, the terms proportional to $W$ contribute to $\mathcal{L}_2$, the terms proportional to $\hbar^2 c^2 / MT$ to $\mathcal{L}_3$ and the terms in $\hbar^4 c^4 / M^2 T^2$ to $\mathcal{L}_4$. Note the presence in $\mathcal{L}_4$ of terms proportional to $p_{\text{tot}}^2$, where $p_{\text{tot}}$ is the center of mass momentum. These terms survive the “elimination” of the center of mass, and they limit the application of the formalism to small $p_{\text{tot}}$ [15].

$$ \begin{aligned} \mathcal{L}^{ss} &= \frac{i\hbar^2 c^2}{M} [\partial_s^2 - \partial_{s'}^2] - iC_F [V(s) - V(s')] + C_F [2\tilde{W}(0) - \tilde{W}(s) - \tilde{W}(s')] \\ &\quad - \frac{\hbar^2 c^2 C_F}{2MT} [\tilde{W}'(s)\partial_s + \tilde{W}'(s')\partial_{s'}] + \frac{\hbar^4 c^4 C_F}{64M^2T^2} [4\tilde{W}''(0)(\partial_s^2 + \partial_{s'}^2) - 2p_{\text{tot}}^2\tilde{W}''(0)] \\ &\quad + \frac{\hbar^4 c^4 C_F}{64M^2T^2} [4\tilde{W}''(s)\partial_s^2 + 4\tilde{W}''(s')\partial_{s'}^2 + 4\tilde{W}'''(s)\partial_s] \\ &\quad + \frac{\hbar^4 c^4 C_F}{64M^2T^2} [4\tilde{W}'''(s')\partial_{s'} + p_{\text{tot}}^2\tilde{W}''(s) + p_{\text{tot}}^2\tilde{W}''(s')], \end{aligned} $$

<footnote>$^{28}$For higher dimensions, the expression reads $\mathcal{L}_{3a}^{ss} = -\frac{C_F \lambda_{\text{th}}^2}{2} (\nabla \tilde{W}(\vec{s}) \cdot \nabla_s + \nabla \tilde{W}(\vec{s}') \cdot \nabla_{s'})$; details of the derivation can be found in [53], section 4.4.2.</footnote>

<footnote>$^{29}$Note that this expression is neither the $\mathcal{L}_{3a}$ nor the $\mathcal{L}'_{3a}$ obtained in [14].</footnote>

<footnote>$^{30}$In higher dimensions, this is $\mathcal{L}_{4a}^{ss} = \frac{C_F \lambda_{\text{th}}^4}{16} [\partial_i \partial_j \tilde{W}(0) (\partial_{s,i} \partial_{s,j} + s \leftrightarrow s') + (\partial_{s,i} (\partial_i \partial_j \tilde{W}(s) \partial_{s,j}) + s \leftrightarrow s')]$.</footnote>
$$
\begin{aligned}
\mathcal{L}^{\text{so}} = C_F & \left[ -2\tilde{W}\left(\frac{s-s'}{2}\right) + 2\tilde{W}\left(\frac{s+s'}{2}\right) \right] \\
& + \frac{\hbar^2 c^2 C_F}{2MT} \left[ -\tilde{W}'\left(\frac{s-s'}{2}\right) (\partial_s - \partial_{s'}) + \tilde{W}'\left(\frac{s+s'}{2}\right) (\partial_s + \partial_{s'}) \right] \\
& + \frac{\hbar^4 c^4 C_F}{8M^2 T^2} \left[ \tilde{W}''\left(\frac{s-s'}{2}\right) \partial_s \partial_{s'} + \tilde{W}''\left(\frac{s+s'}{2}\right) \partial_s \partial_{s'} \right] \\
& + \frac{\hbar^4 c^4 C_F}{32M^2 T^2} \left[ p_{\text{tot}}^2 \tilde{W}''\left(\frac{s-s'}{2}\right) - p_{\text{tot}}^2 \tilde{W}''\left(\frac{s+s'}{2}\right) \right],
\end{aligned}
$$

$$
\begin{aligned}
\mathcal{L}^{\text{os}} ={}& \frac{1}{2N_c} \left[ -2\tilde{W}\left(\frac{s-s'}{2}\right) + 2\tilde{W}\left(\frac{s+s'}{2}\right) \right] \\
& + \frac{1}{2N_c} \frac{\hbar^2 c^2}{2MT} \left[ -\tilde{W}'\left(\frac{s-s'}{2}\right) (\partial_s - \partial_{s'}) + \tilde{W}'\left(\frac{s+s'}{2}\right) (\partial_s + \partial_{s'}) \right] \\
& + \frac{1}{2N_c} \frac{\hbar^4 c^4}{8M^2 T^2} \left[ \tilde{W}''\left(\frac{s-s'}{2}\right) \partial_s \partial_{s'} + \tilde{W}''\left(\frac{s+s'}{2}\right) \partial_s \partial_{s'} \right] \\
& + \frac{1}{2N_c} \frac{\hbar^4 c^4}{32M^2 T^2} \left[ p_{\text{tot}}^2 \tilde{W}''\left(\frac{s-s'}{2}\right) - p_{\text{tot}}^2 \tilde{W}''\left(\frac{s+s'}{2}\right) \right],
\end{aligned}
$$

$$
\begin{aligned}
\mathcal{L}^{\text{oo}} ={}& \frac{i\hbar^2 c^2}{M} [\partial_s^2 - \partial_{s'}^2] + i\frac{1}{2N_c} [V(s) - V(s')] + 2C_F \tilde{W}(0) + \frac{1}{2N_c} [\tilde{W}(s) + \tilde{W}(s)] \\
& - \frac{N_c^2 - 2}{N_c} \tilde{W}\left(\frac{s-s'}{2}\right) - \frac{2}{N_c} \tilde{W}\left(\frac{s+s'}{2}\right) + \frac{1}{2N_c} \frac{\hbar^2 c^2}{2MT} [\tilde{W}'(s)\partial_s + \tilde{W}'(s')\partial_{s'}] \\
& - \frac{N_c^2 - 2}{2N_c} \frac{\hbar^2 c^2}{2MT} \tilde{W}'\left(\frac{s-s'}{2}\right) (\partial_s - \partial_{s'}) - \frac{1}{N_c} \frac{\hbar^2 c^2}{2MT} \tilde{W}'\left(\frac{s+s'}{2}\right) (\partial_s + \partial_{s'}) \\
& + \frac{\hbar^4 c^4 C_F}{64M^2 T^2} [4\tilde{W}''(0) (\partial_s^2 - \partial_{s'}^2) - 2p_{\text{tot}}^2 \tilde{W}''(0)] \\
& - \frac{1}{2N_c} \frac{\hbar^4 c^4}{64M^2 T^2} [4\tilde{W}''(s)\partial_s^2 + 4\tilde{W}''(s')\partial_{s'}^2 + 4\tilde{W}'''(s)\partial_s] \\
& - \frac{1}{2N_c} \frac{\hbar^4 c^4}{64M^2 T^2} [4\tilde{W}'''(s')\partial_{s'} + p_{\text{tot}}^2 \tilde{W}''(s) + p_{\text{tot}}^2 \tilde{W}''(s')] \\
& + \frac{N_c^2 - 2}{2N_c} \frac{\hbar^4 c^4 C_F}{64M^2 T^2} [8\tilde{W}''\left(\frac{s-s'}{2}\right) \partial_s \partial_{s'} + 2p_{\text{tot}}^2 \tilde{W}''\left(\frac{s-s'}{2}\right)] \\
& - \frac{1}{N_c} \frac{\hbar^4 c^4 C_F}{64M^2 T^2} [8\tilde{W}''\left(\frac{s+s'}{2}\right) \partial_s \partial_{s'} - 2p_{\text{tot}}^2 \tilde{W}''\left(\frac{s+s'}{2}\right)].
\end{aligned}
\qquad (C.1)
$$

D Details on the numerical implementation

The equations are solved numerically for the charmonium system (with the charm mass $m_c = 1.469$ GeV), using the Crank-Nicolson scheme. Both the s and s' variables (which correspond to the relative distance between the quark and antiquark) are discretized along $N = 501$ points, from -10 fm to +10 fm (box size $L = 20$ fm), with a spatial step $\Delta s = 0.04$ fm. We choose a time step $\Delta t = 0.01$ fm/c (lower values of $\Delta t$ were found to have no effect on
<figure><img src="image_33.png" /><figcaption>**Figure 30.** Real (left panel) and regularized imaginary (right panel) parts of the one-dimensional potential as a function of the medium temperature *T* (in GeV) and distance *x* between the quark and the antiquark of the charmonium [39]. The color factor *C*<sub>*F*</sub> is that appropriate for the singlet channel. The dashed line in the right panel represents the probability density corresponding to the 1S charmonium state.</figcaption></figure>

the results). The boundary conditions are such that $\mathcal{D}(s, s') = 0$ when either $s$ or $s'$ is on the boundary of the box, that is when $s = \pm 10$ fm or $s' = \pm 10$ fm.

We use a complex potential (shown in figure 30) that was specifically developed for one-dimensional studies [39]. The real and imaginary parts of the complex potential were tuned to reproduce at best the temperature-dependent mass spectra and decay widths of quarkonium states obtained with a three-dimensional potential model inspired by lattice QCD results [44].<fnref n="31" />

The real potential admits three eigenstates (referred to as vacuum eigenstates) which are obtained by solving the Schrödinger equation for the relative c̄c motion. The ground state is referred to as 1S, while the first (odd) and second (even) excited states are referred to respectively as 1P and 2S.

The flattening of the 1D potential at large distance is adjusted to match the corresponding behavior of the 3D potential and to mimic the temperature dependence of the charm pole-mass, as well as the effect of the “string breaking” mechanism.

Table 3 shows that the properties of the S states of the charmonium as calculated from a 3D-potential are well accounted for by the 1D potential. Note that the 3D and the 1D calculations both assume that temperature effects are irrelevant below a temperature $T_0$. Thus, the vacuum properties referred to in this paper correspond in fact as well to those calculated at this temperature $T_0$. For the P-state, the quantitative agreement with the 3D model is not as good, which is probably due to the absence of centrifugal energy in the 1D setup (see [39]).

Besides, one has $2m_c + V_{SB} \approx 2m(D^+) -$ where $V_{SB} \approx 0.8352$ GeV corresponds to the saturation for the string breaking 1D-potential — so that the liberation threshold for a c̄c pair can roughly be seen as the creation of a $D^+ + D^-$ pair for $T = T_0$. The binding energies at larger temperatures can be read directly from figure 31. With increasing temperature,

<footnote n="31">³¹Note that in [39], as well as in [44], the $C_F$ color factor associated to the singlet state was absorbed in the definition of the complex potential. In the present work, we make the color factors explicit, as they differ between the singlet and octet parts of the $\mathcal{L}$ superoperator, as recalled in appendix C; this holds both for the real part $V$ and the imaginary part $W$ of the potential.</footnote>
<table><caption>**Table 3.** Masses (in GeV) and root mean square radii (in fm) for the “vacuum” S states (at $T_0 \approx 0.126025$ GeV). $m^{\text{PDG}}$ is the experimental mass given by the particle data group (for the 1P, it corresponds to a weighted average between $\chi_{c0}$, $\chi_{c1}$, and $\chi_{c2}$), $m_{3D}^{\text{EV}}$ (resp. $\sqrt{\langle r^2 \rangle}$) is the mass (resp. root mean square radius) calculated from the expectation values of the Hamiltonian (the square radius operator) with the three-dimensional potential suggested in [44] (no spin-orbit term). $m_{1D}^{\text{EV}}$ and $\sqrt{\langle x^2 \rangle}$ are the corresponding values obtained with the one-dimensional potential.</caption><thead><tr><th>At T = T<sub>0</sub></th><th>J/ψ(1S)</th><th>χ<sub>c</sub>(1P)</th><th>ψ(2S)</th></tr></thead><tbody><tr><td>m<sup>PDG</sup> [54]</td><td>3.0969</td><td>3.525</td><td>3.6861</td></tr><tr><td>m<sup>EV</sup><sub>3D</sub></td><td>3.0964</td><td>3.509</td><td>3.6642</td></tr><tr><td>m<sup>EV</sup><sub>1D</sub></td><td>3.0981</td><td>3.4532</td><td>3.6858</td></tr><tr><td>|m<sup>PDG</sup> - m<sup>EV</sup><sub>1D</sub>|</td><td>0.0012</td><td>0.072</td><td>0.0003</td></tr><tr><td>E<sub>binding</sub></td><td>0.675</td><td>0.320</td><td>0.087</td></tr><tr><td>3D √⟨r<sup>2</sup>⟩</td><td>0.431</td><td>0.689</td><td>0.943</td></tr><tr><td>1D √⟨x<sup>2</sup>⟩</td><td>0.271</td><td>0.543</td><td>0.856</td></tr></tbody></table>

<figure><img src="image_34.png" /><figcaption>**Figure 31.** Binding energy as a function of the temperature for the 3 vacuum bound states 1S, 1P and 2S, considered to be constant up to $T_0 \approx 126$ MeV. Vertical long-dashed lines indicate the melting temperatures, corresponding, for practical purpose, to $E_{\text{binding}} = 1$ MeV. Note that the unnatural behavior of the binding energies for $T < 200$ MeV (darkened area) is an artifact of the specific potential used in this paper, and should not be regarded as physical. All our numerical simulations are done for $T \gtrsim 200$ MeV.</figcaption></figure>

the potential well gets shallower (see figure 30): the melting temperatures of the 1S, 1P and 2S charmonium state are defined when $E_{\text{binding}} < 1$ MeV, respectively $T_{\text{melt}} \approx 440$ MeV, $\approx 280$ MeV and $\approx 250$ MeV.

The imaginary part $W$ of the potential is taken to be:

$$W(x, T) = \alpha W_C(|x|, T) + \beta W_S(|x|, T), \qquad (\text{D.1})$$

with $\alpha = 1.7$, $\beta = 0.8$ for the charmonium system, and where $W_C$ and $W_S$ are the Coulombic
and string-like parts, defined as:

$$ 
\begin{aligned}
W_C(x, T) &= \frac{\tilde{\alpha}_s T}{C_F} \phi(m_D x), \quad \phi(r) = 2 \int_0^\infty dz \frac{z}{(z^2+1)^2} \left(1 - \frac{\sin(zr)}{zr}\right) \\
W_S(x, T) &= \frac{\sigma T}{C_F m_D^2} \chi(m_D x), \quad \chi(r) = 2 \int_0^\infty dp \frac{2 - 2 \cos(pr) - pr \sin(pr)}{\sqrt{p^2 + \Delta_D^2 (p^2 + 1)^2}},
\end{aligned}
\quad (D.2) 
$$

where $\tilde{\alpha}_S = 0.513 \pm 0.0024$, $\sigma = (0.412 \text{ GeV})^2$, and $\Delta_D \approx 3.0369$. The $\alpha$ and $\beta$ coefficients were chosen to reproduce at best the decay widths calculated with the three-dimensional potential. At small distance, the imaginary potential exhibits a harmonic behavior$^{32}$ which can be linked to the heavy quark diffusion coefficient. At large distance, the potential saturates toward twice the reaction rate of individual c quarks [43].

The imaginary potential given at eq. (D.1) suffers from the same UV divergences as the HTL imaginary potential recalled in appendix B and is regularized as discussed in that appendix. The right panel of figure 30 illustrates the ensuing result, as well as the probability $\rho_{1S}$ for the vacuum 1S state, which allows to appreciate the action of W on the ground state.

## E Charmonia dissociation rates

In this appendix, we evaluate the various transition rates that are referred to in the main text. Let $|n\rangle$ be a normalized singlet state, and let us assume that the initial density matrix $D_s(0)$ is the projector on this particular state, i.e. $D_s(0) = |n\rangle\langle n|$. The variation with time of the probability $p_n = \langle n|D_s(t)|n\rangle$ is given in eq. (A.4), from which we get, taking into account that $p_o(t=0) = 0$,

$$ \left. \frac{dp_n}{dt} \right|_{t=0} = \langle n | \mathcal{L}_a^{ss} \mathcal{D}(0) | n \rangle = -\Gamma_n, \quad (E.1) $$

which defines the dissociation rate $\Gamma_n$. In this equation, $\mathcal{L}_a^{ss}$ is the $\mathcal{L}_a$ part of the Lindblad superoperator $\mathcal{L}^{ss}$ that involves the imaginary potential W, viz.

$$ \mathcal{L}_a^{ss} = \mathcal{L}_{2a} + \mathcal{L}_{3a} + \mathcal{L}_{4a}. $$

We call $\Gamma_{n,2}$, $\Gamma_{n,3}$ and $\Gamma_{n,4}$ the corresponding contributions to $\Gamma_n$. These are estimated in appendix B, see eqs. (B.6), (B.9), and (B.11). Their sum can be written as

$$ \Gamma_n = C_F \int dq \, \tilde{\rho}(q) \int ds \left| \left( \sin \left( \frac{qs}{2} \right) + \frac{\lambda_{\text{th}}^2}{4} \cos \left( \frac{qs}{2} \right) q \partial_s \right) \psi_n(s) \right|^2. \quad (E.2) $$

It is clearly positive. Incidentally, one can take advantage of this explicit expression to evaluate the hierarchy between the 3 contributions, starting by the $\Gamma_{n,3}/\Gamma_{n,2}$ ratio. As compared to the evaluation of $\Gamma_{n,2}$ — the contribution involving the $\sin(\cdots)^2$ term in the integrand —, the calculation of $\Gamma_{n,3}$ (the cross term) generates: a) one extra derivative of the wave function, b) one extra q factor which is equivalent — integrating by parts — to another derivative of the wave function, and c) an extra $\lambda_{\text{th}}^2$ factor. Parametrically, each

<footnote>$^{32}$Up to some logarithmic prefactor.</footnote>
<figure><img src="image_35.png" /><figcaption>Figure 32. Left: decay rates of the vacuum 1S state (aka $J/\psi$) as a function of the temperature. The respective contributions of $\mathcal{L}_{2a}$ (blue, dashed), $\mathcal{L}_{2a} + \mathcal{L}_{3a}$ (orange, long-dashed) as well as the full result using $\mathcal{L}_{2a} + \mathcal{L}_{3a} + \mathcal{L}_{4a}$ (red, full) are shown. The calculations are done using the regularized imaginary potential of eq. (D.1), following the regularization procedure outlined in appendix B. Right: same calculation for the lowest in-medium 1S state of the T-dependent potential.</figcaption></figure>

extra derivative of the wave function $\psi'_n$ generates a factor inversely proportional to the rms radius $\sqrt{\langle r^2 \rangle_n}$ of the state. Hence, the ratio $\Gamma_{n,3}/\Gamma_{n,2}$ scales like

$$ \frac{\Gamma_{n,3}}{\Gamma_{n,2}} \sim \frac{\lambda_{\text{th}}^2}{\langle r^2 \rangle_n}. \qquad (\text{E.3}) $$

Since $\lambda_{\text{th}}$ decreases with increasing temperature, while $\langle r^2 \rangle_n$ increases with $T$ and $n$, these ratios are decreasing functions of $T$ and $n$. They are moreover expected to be less than unity in the QBM regime.<fnref n="33" /> A similar analysis results in $\Gamma_{n,4}/\Gamma_{n,2} \propto (\Gamma_{n,3}/\Gamma_{n,2})^2$. For all temperatures lying in the QBM regime, one thus expects $\Gamma_{n,3}$ to be small and $\Gamma_{n,4}$ subleading.

The calculations of the rates of the various singlet states are displayed as a function of the medium temperature in figures 32, 33 and 34. Figure 32 shows the contributions of the individual superoperators to both the vacuum (left) and the in-medium (right) 1S ($J/\psi$) state. In the case of the vacuum 1S state, one sees that the $\mathcal{L}_{3,a}$ term generates a significant reduction as compared to the dominant $\mathcal{L}_{2,a}$ term. The $\mathcal{L}_{4,a}$ term is completely subdominant for $T \gtrsim 0.25$ GeV. For low temperature ($T \lesssim 0.2$ GeV, not displayed), it is needed for preserving the positivity of the decay rate, although one may note that the positivity violation occurs only at the lower edge of the QBM regime. The in-medium eigenstate has a larger decay rate. This is related to the fact that it is less bound, and correlatively it has a larger spatial spread than the vacuum states.

On the left panel of figure 33, we compare the decay rates of the three vacuum-eigenstates of the 1D potential and find as anticipated that the deepest bound state has the smallest decay rate and is thus is the most robust against dissociation. With increasing temperature, $W$ reaches its asymptotic value for smaller and smaller distance (see figure 30, right), and all decay rates converge toward $2C_F W(\infty, T) \propto T$ which corresponds to twice the damping rate of an individual c quark.

<footnote n="33">³³For the 1S state, one has indeed $\frac{\lambda_{\text{th}}^2}{\langle r^2 \rangle_{1S}} \sim \lambda_{\text{th}}^2 \langle p^2 \rangle_{1S} \sim \frac{E_{\text{kin}}^{1S}}{T} \lesssim 1$ when $T$ becomes larger than the typical energies associated to the bound state.</footnote>
<figure><img src="image_36.png" /><figcaption>**Figure 33.** Left: full decay rate ($\mathcal{L}_{2a} + \mathcal{L}_{3a} + \mathcal{L}_{4a}$) for 1S (plain, orange), 1P (long-dashed, green) and 2S (dashed, red) states of the vacuum potential; same method as figure 32. On the top, the limiting value $2C_F W(\infty, T)$ is represented by a black dotted line. Right: same calculation with in-medium states.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Figure 34.** Fraction of the decay rate due to the Coulomb-like part of the imaginary potential (same conventions as for figure 33).</figcaption></figure>

On the right panel of figure 33, the decay rates of the in-medium states are illustrated. As these states are broader than the corresponding vacuum states, their decay rate are significantly larger — as already noticed on figure 32 — and even saturate at finite $T$ to $2C_F W(\infty, T)$ when the bound state ceases to exist. For $T \gtrsim 200$ MeV, one observes an approximate scaling $\Gamma \propto T^4$. This increase with $T$ slightly exceeds the law (2.14) due to the aforementioned growth of the states with $T$.

Finally, figure 34 illustrates the fraction of the decay width that stems from the Coulomb-like part of the imaginary potential (D.1). As the string-like part behaves as $x^4$ at small distance $x$, it is less important for 1S state than for the 2S or the 1P. But even at high temperature this strong contribution remains significant for all the states. Note that while this contribution is quantitatively significant, as shown in figure 34, we have not identified any qualitative impact on the physical quantities that we have evaluated.
## F Evolution under $\mathcal{L}_2$ alone

In this appendix, we analyze the evolution of the system under the sole action of the superoperator $\mathcal{L}_2$, which dominates the dynamics at large $T$ ($T \gtrsim T_{\text{damp}}$). This will give us the opportunity to comment on the role of coordinate space as a preferred basis (in the sense where this terminology is used in the theory of open quantum systems [47]): the superoperator $\mathcal{L}_2$ is a local operator in the coordinate space representation, which results in a rapid “diagonalization” of the density matrix as time increases; this illustrates the effect of collisional decoherence in inducing a partial “classicalization” of the system. The color degree of freedom remains quantum though, and this is manifest in the fact that the singlet and octet components of the density matrix remain coupled in their equations of motion.

The explicit equations for $\mathcal{D}_s$ and $\mathcal{D}_o$ read (with $\dot{\mathcal{D}} = d\mathcal{D}/\partial t$)

$$ \begin{pmatrix} \dot{\mathcal{D}}_s \\ \dot{\mathcal{D}}_o \end{pmatrix} = \begin{pmatrix} C_F (2W(0) - W_c) & -C_F W^- \\ -\frac{1}{2N_c} W^- & 2C_F W(0) + \frac{1}{2N_c} W_c - \left(\frac{N_c^2-2}{2N_c} W_a + \frac{W_b}{N_c}\right) \end{pmatrix} \begin{pmatrix} \mathcal{D}_s \\ \mathcal{D}_o \end{pmatrix}, \quad (\text{F.1}) $$

where$^{34}$ $W_a = 2W(\frac{s-s'}{2})$, $W_b = 2W(\frac{s+s'}{2})$, $W^- = W_a - W_b$, and $W_c = W(s) + W(s')$. Introducing $\Gamma(r) = W(r) - W(0)$, as well as the analogously defined quantities $\Gamma_a, \Gamma_b$ and $\Gamma^- = \Gamma_a - \Gamma_b$, we can rewrite eqs. (F.1) as

$$ \begin{pmatrix} \dot{\mathcal{D}}_s \\ \dot{\mathcal{D}}_o \end{pmatrix} = \begin{pmatrix} -C_F (\Gamma(s) + \Gamma(s')) & 2C_F \left(\Gamma\left(\frac{s+s'}{2}\right) - \Gamma\left(\frac{s-s'}{2}\right)\right) \\ \frac{1}{N_c} \left(\Gamma\left(\frac{s+s'}{2}\right) - \Gamma\left(\frac{s-s'}{2}\right)\right) & -\frac{1}{2N_c} (\Gamma(s) + \Gamma(s')) + \mathcal{L}_{o \leftrightarrow o} \end{pmatrix} \begin{pmatrix} \mathcal{D}_s \\ \mathcal{D}_o \end{pmatrix}, \quad (\text{F.2}) $$

where we have introduced the $\mathcal{L}_{o \leftrightarrow o}$ superoperator:

$$ \mathcal{L}_{o \leftrightarrow o} = \frac{1}{N_c} (\Gamma(s) + \Gamma(s')) - \left( \frac{N_c^2 - 2}{2N_c} \Gamma_a + \frac{1}{N_c} \Gamma_b \right), $$

which leads to transitions exclusively happening within the octet component. It can be verified that $\mathcal{L}_{o \leftrightarrow o}(s' = s) = 0$, so that $\mathcal{L}_{o \leftrightarrow o}$ does not modify the octet weight (tr$\mathcal{D}_o$) nor the singlet one (tr$\mathcal{D}_s$).

Eqs. (F.2) do not only preserve the trace of the density matrix, but also the total density of heavy quark pairs, $\rho_s(r, t) + (N_c - 1)^2\rho_o(r, t)$, where $\rho_s(r, t) = \mathcal{D}_s(r, r, t)$ and $\rho_o(r, t) = \mathcal{D}_o(r, r, t)$. It can indeed be easily verified that this total density remains constant in time, for any $r$. To see that, let us set $r = s = s'$ in eqs. (F.2), and note that $\Gamma(0) = 0$. Then these equations simplify into

$$ \begin{pmatrix} \dot{\rho}_s \\ \dot{\rho}_o \end{pmatrix} = \Gamma(r) \begin{pmatrix} -2C_F & 2C_F \\ \frac{1}{N_c} & -\frac{1}{N_c} \end{pmatrix} \begin{pmatrix} \rho_s \\ \rho_o \end{pmatrix}. \quad (\text{F.3}) $$

The property mentioned above then follows easily. We derive also the following equations

$$ \frac{d\rho_8}{dt} = -N_c\Gamma(r)\rho_8, \quad \rho_8(r,t) = \mathcal{D}_8(r,r,t), \quad \mathcal{D}_8 \equiv \frac{1}{N_c^2}(\mathcal{D}_s - \mathcal{D}_o), \quad (\text{F.4}) $$

which indicate that the system is driven to local color equilibrium ($\rho_s = \rho_o$) at a rate $N_c\Gamma(r)$ which depends on the (fixed) size $r$ of the $Q\bar{Q}$ pair.

<footnote>$^{34}$We use here the notation of ref. [14].</footnote>
To proceed further, we approximate the imaginary part of the potential W by a harmonic expansion, $W(x) \approx W(0) + \frac{\kappa}{2}x^2$.<fnref n="35" /> The densities are then given by

$$ \rho_s(r,t) - \rho_o(r,t) = e^{-N_c r^2 \tau} \rho_0(r), \quad \rho_s(r,t) = \left[ \frac{N_c^2 - 1}{N_c^2} e^{-N_c r^2 \tau} + \frac{1}{N_c^2} \right] \rho_0(r), \qquad (F.5) $$

where $\rho_0(r) = \rho_s(r, t = 0)$ is the initial singlet density ($\rho_o(r, t = 0) = 0$), and where we have introduced the reduced time $\tau = \kappa t/2$. At late time, the densities tend to a common value that is independent of time, $\rho_s^{\text{asymp}}(r) = \rho_o^{\text{asymp}}(r) = \frac{1}{N_c^2}\rho_0(r)$. This corresponds to a maximum entropy state for the color degrees of freedom: all the color states are equally populated, with the densities varying with $r$ according to the initial density profile $\rho_0(r)$. A similar result holds if the initial density is an octet density, i.e., $\rho_o(r, t = 0) = \rho_0(r)$, and $\rho_s(r, t = 0) = 0$. One finds $\rho_o(r, t) = \frac{1}{N_c^2}[N_c^2 - 1 + e^{-N_c r^2 \tau}]\rho_0(r)$.

From the densities $\rho_s(r, t)$ and $\rho_o(r, t)$, one can calculate the evolution in time of the singlet and octet weight, respectively $\text{tr}\mathcal{D}_s(t)$ and $\text{tr}\mathcal{D}_o(t)$. In order to get explicit expressions, we choose as initial (singlet) density the normalized Gaussian

$$ \rho_0(r) = \frac{1}{\sqrt{\pi}\sigma^2} e^{-r^2/\sigma^2}. \qquad (F.6) $$

We then obtain

$$ \text{tr}\mathcal{D}_s(t) = \frac{N_c^2 - 1}{N_c^2} \frac{1}{\sqrt{1 + \sigma^2 N_c \tau}} + \frac{1}{N_c^2}. \qquad (F.7) $$

This is a much slower decay than the exponential relaxation of the density toward its equilibrium value, as indicated in eq. (F.5). To understand in more details how this behavior emerges, let us expand $\text{tr}\mathcal{D}_s(\tau)$ in leading order in $\tau$. The leading order contributions to $\mathcal{D}$ read

$$ \mathcal{D}_s^{(1)}(\tau) = [1 - C_F(s^2 + s'^2)\tau] \mathcal{D}_s(0) \quad \text{and} \quad \mathcal{D}_o^{(1)}(\tau) = \frac{1}{N_c} ss'\tau \mathcal{D}_s(0). \qquad (F.8) $$

From these equations, or from the direct expansion of (F.7), it follows that

$$ \text{tr}\mathcal{D}_s^{(1)} = 1 - \frac{1}{2}\sigma^2\frac{N_c^2 - 1}{N_c}\tau = 1 - 2C_F\langle r^2 \rangle_0\tau = 1 - \Gamma_0 t, \qquad (F.9) $$

where $\Gamma_0 = C_F\kappa\langle r^2 \rangle_0$, with $\langle r^2 \rangle_0 = \sigma^2/2$, plays the role of a damping rate.<fnref n="36" /> This leading order result, however, does not exponentiate into $e^{-\Gamma_0 t}$. At second order we have indeed

$$ \text{tr}\mathcal{D}_s^{(2)}(t) = 1 - 2C_F\langle r^2 \rangle_0\tau + \left(2C_F^2 + \frac{C_F}{N_c}\right)\langle r^4 \rangle_0\tau^2. \qquad (F.10) $$

The term in $\tau^2$ obviously differs from $(1/2)\Gamma_0^2 t^2$ coming from the expansion of $e^{-\Gamma_0 t}$. There are two effects responsible for this. One is a “correction” coming from the coupling between the singlet and octet sectors, represented by the term $C_F/N_c$ which adds to the term $2C_F^2$ in the coefficient of $\tau^2$. Such a coupling term is absent in the term linear in $\tau$ if $\rho_o = 0$.

<footnote n="35">³⁵This approximation is a priori valid only when $x \lesssim m_D^{-1}$. We use it here for all values of $x$ in order to get simple analytical expressions.</footnote>

<footnote n="36">³⁶This is in agreement with the definition (E.1), with here $\mathcal{L}_{2a} = -2C_F\Gamma(r)$ and $\Gamma(r) = \kappa r^2/2$.</footnote>
<figure><img src="image_38.png" /><figcaption>**Figure 35.** Same as figure 4 for the evolution performed solely with $\mathcal{L}_2$. $\Gamma_0$ in the harmonic approximation of W has been set to 100 MeV, corresponding to the decay-width of a 1S in-medium state — $\langle r^2 \rangle_0 \approx 0.125 \text{ fm}^2$ — at $T = 300$ MeV. The slope of the dashed black line is minus this decay rate.</figcaption></figure>

at the start of the evolution, which is the considered case. The second effect comes from the fact that $\langle r^4 \rangle_0 \neq \langle r^2 \rangle_0^2$, where $\langle \cdot \rangle_0$ denotes an average over the initial density $\rho_0(r)$. For the Gaussian density (F.6), $\langle r^4 \rangle_0 = 3\langle r^2 \rangle_0^2$. Thus, even ignoring the back-contribution of the octet sector, there is a deviation from the exponential behavior that stems simply from the fact that the average of the exponential is not the exponential of the average. The comparison between the leading order terms and the exact behavior (F.7) reveals the intricate coupling between the color dynamics and the evolution of the singlet and octet densities. Note in particular how the color prefactor in eq. (F.9) contributes to change the time scale of the evolution, from $(N_c\kappa\langle r^2 \rangle_0)^{-1}$ to $(C_F\kappa\langle r^2 \rangle_0)^{-1}$. Note also that, leaving aside the color dynamics that couple strongly the singlet and octet sectors, the calculation above requires only “classical” information, meaning that it involves only the diagonal components of the coordinate space density matrices $D_s$ and $D_o$, i.e. the densities of singlet and octet $Q\bar{Q}$ pairs. We shall consider later in this appendix quantities that are sensitive to the non-diagonal matrix elements $D(s, s')$.

A plot of the singlet and octet weights as a function of time is given in figure 35. The octet weight is obtained from the trace of the full density matrix: $\text{tr}D = 1 = \text{tr}D_s + (N_c^2 - 1)\text{tr}D_o$. The qualitative behavior is quite similar to that observed in the left panel of figure 4. However it appears that color equilibration is taking here a much longer time. This can, to a large extent, be attributed to the fact that, in contrast to the situation in figure 4, the role of the continuum states is here hindered by the fact that the entire density stays localized in the region defined by the initial density $\rho_0(r)$. In comparison with the system considered in section 3.1, a $Q\bar{Q}$ pair in a box of size 20 fm, the typical values of $r$ involved in the present calculation are “small”, meaning that the damping rate is small, hence the longer time scale observed for color equilibration.

Another characterization of the variation of the densities with time is provided by the variation of the expectation value of $r^2$ (see e.g. figure 17). For the singlet density, we obtain

$$ \frac{\langle r^2 \rangle_s}{\langle r^2 \rangle_0} = \frac{1 + \frac{N_c^2 - 1}{(1 + N_c \sigma^2 \tau)^{3/2}}}{1 + \frac{N_c^2 - 1}{(1 + N_c \sigma^2 \tau)^{1/2}}}. \qquad (F.11) $$
The variation of $\langle r^2 \rangle_s$ with time is a genuine non-abelian effect, i.e. an effect which results from the coupling between the singlet and octet sectors. Indeed, in an abelian theory the density would remain constant under the action of $\mathcal{L}_2$ alone. To see that, consider the abelian version of eqs. (F.2)

$$ \dot{D}_{\text{QED}} = \left[ -C_F (\Gamma(s) + \Gamma(s')) + 2C_F \left( \Gamma\left(\frac{s+s'}{2}\right) - \Gamma\left(\frac{s-s'}{2}\right) \right) \right] D_{\text{QED}}(s, s', t). \quad (\text{F.12}) $$

The solution, in the harmonic approximation, reads simply

$$ \mathcal{D}_{\text{QED}}(s, s', t) = e^{-\frac{C_F \kappa t}{2} (s-s')^2} \mathcal{D}_{\text{QED}}(s, s', 0), \quad (\text{F.13}) $$

which shows that the main effect of collisions in the abelian case is momentum broadening (see below). It follows from this equation that $\rho_{\text{QED}}(r, t) = \mathcal{D}_{\text{QED}}(r, r, t) = \rho_0(r)$ is a constant, and so is therefore $\langle r^2 \rangle$. The equation $\rho_{\text{QED}}(r, t) = \rho_0(r)$ reflects the local conservation of the total number of $Q\bar{Q}$ pairs. The same conservation law in the case of QCD applies to the combination $\rho_s + (N_c^2 - 1)\rho_0$, as we have already emphasized. What happens in the QCD case is that the singlet to octet transitions start to populate the octet states at large $r$ (where $\Gamma(r)$ is the largest). This naturally entails a variation of the singlet density, reflected in the decrease of $\langle r^2 \rangle_t$ according to the formula above. At small time, we have approximately $\langle r^2 \rangle_s/\langle r^2 \rangle_0 \simeq 1 - 2\Gamma_0 t$. At late time however, the populations equilibrate, as we have seen earlier, and in particular $\rho_s(r, t) \simeq \rho_0(r)/N_c^2$. When this situation is reached, $\langle r^2 \rangle_s \approx \langle r^2 \rangle_0$ again. Note that the variation of $\langle r^2 \rangle_s$ concerns only the singlet density. Because the total density of $Q\bar{Q}$ pairs is locally conserved, it remains independent of time, and so does the corresponding rms radius. As for the octet, its rms radius starts to grow at early time, reaches some maximum and then decreases toward $\langle r^2 \rangle_0$ for $t \to \infty$.

We turn now to the solution for the general (non-diagonal) case. This will allows us to exhibit the general structure of $\mathcal{D}(s, s')$ and calculate some quantities which require non-diagonal information, such as the average momentum square, or survival probabilities. In the harmonic approximation, eqs. (F.2) read

$$ \left( \frac{d\mathcal{D}_s}{d\tau} \middle/ \frac{d\mathcal{D}_o}{d\tau} \right) \approx \left( \begin{matrix} -C_F (s^2 + s'^2) & 2C_F s s' \\ \frac{1}{N_c} s s' & -\frac{1}{2N_c} (s^2 + s'^2) - \frac{(N_c^2-4)(s-s')^2}{4N_c} \end{matrix} \right) \left( \begin{matrix} \mathcal{D}_s \\ \mathcal{D}_o \end{matrix} \right). \quad (\text{F.14}) $$

These equations admit explicit solutions. We focus on solutions corresponding to an initial condition such that $\mathcal{D}_o(s, s', t = 0) = 0$. One then gets:

$$ \mathcal{D}_s(s, s', \tau) = U^{ss}(s, s', \tau) \mathcal{D}_s(s, s', 0) \quad (\text{F.15}) $$

with

$$ U^{ss}(s, s', \tau) = e^{-\left(\frac{N_c}{8}(s+s')^2 + \frac{N_c^2-4}{4N_c}(s-s')^2\right)\tau} \left( \cosh(\alpha\tau) - \left(\frac{N_c^2-2}{2}(s+s')^2 + (s-s')^2\right) \frac{\sinh(\alpha\tau)}{4N_c\alpha} \right), \quad (\text{F.16}) $$

where

$$ \alpha(s, s') = \frac{1}{8} \sqrt{N_c^2 (s + s')^4 - 16 s s' (s - s')^2}. \quad (\text{F.17}) $$
<figure><img src="image_39.png" /><figcaption>**Figure 36.** Illustration of $U^{ss}(s, s', \tau)$ for $\tau = 5 \text{ fm}^{-2}$. For $\tau = 0$, $U^{ss}(s, s', \tau) = 1$, while at late time $U^{ss}(s, s', \tau)$ is given by eq. (F.18). The damping of the non-diagonal elements is obvious. Near the origin, i.e. for small values of $s$ and $s'$, the damping is slower, as anticipated from the color transparency phenomenon encoded in the harmonic approximation of $W$. This manifests itself explicitly here by the fact that $U^{ss}(s, s', \tau)$ depends on $s, s'$ and $\tau$ only through the combinations $s\sqrt{\tau}$ and $s'\sqrt{\tau}$. At late time the peak conserves its height, but it extends on a small region of size $1/\sqrt{\tau}$ and its weight becomes negligible. Right: same for the $U^{00}$ evolution superoperator in the $D_0 - D_8$ basis.</figcaption></figure>

In the long time limit, $U^{ss}(s, s', \tau)$ shrinks along the $s-s'$ direction and becomes

$$ U^{ss}(s, s', \tau) \approx \frac{1}{N_c^2} \exp \left( - \left( \frac{N_c}{2} r^2 + \frac{N_c^2 - 2}{4N_c} (s-s')^2 - \alpha \right) \times \tau \right) \approx_{\text{finite } r} \frac{1}{N_c^2} e^{-\frac{N_c^2-1}{4N_c} (s-s')^2 \tau} \quad (\text{F.18}) $$

The factor $(s-s')^2\tau$ in the exponent of this expression suggests that at late time, the density matrix $D_s(s, s', t)$ becomes diagonal, $\propto \rho_0(s)\delta(s-s')$, an expected feature of classicalization.<fnref n="37" /> This can be indeed verified in figure 36, which reveals, in addition to the diagonal structure the presence of a peak in the vicinity of the origin. To better encompass this overall structure, it is convenient to use the following approximation of $U^{ss}$:

$$ U^{ss}(s, s', \tau) \approx \frac{N_c^2 - 1}{N_c^2} e^{-N_c \left(r^2 + \frac{y^2}{4}\right) \tau} + \frac{1}{N_c^2} e^{-\frac{N_c^2-1}{4N_c} y^2 \tau}, \quad (\text{F.19}) $$

where we have set $y = s - s'$ and $r = (s + s')/2$. When $s' = s$, this expression coincides with that of the singlet density given in eq. (F.5), and it is equal to unity for $\tau = 0$ as it should. We have checked numerically that this approximate expression provides a reasonably accurate representation of $U^{ss}$.<fnref n="38" /> The formula (F.19) clearly exhibits the peak near $s = s' = 0$. The existence of this peak, which is also present in the density $\rho_s(r, t) = D_s(r, r, t)$, reflects the fact that the damping rate $\Gamma(r)$ vanishes for $r = 0$, so that there are not singlet to octet transitions at this point. The height of the peak stays

<footnote n="37">³⁷For $r = 0$, one has $\alpha = y^2/4$ and $U^{ss'} \approx \frac{1}{N_c^2} \exp(-\frac{N_c^2+N_c-2}{4N_c} (s-s')^2)$, confirming that the shrinking along the $s-s'$ direction acts $\forall r$.</footnote>

<footnote n="38">³⁸Although the $N_c$ dependence of the leading order expansion in $y$ of the expression (F.19) is not correctly reproduced, the ensuing error is quantitatively negligible (see eq. (F.27) and figure 37).</footnote>
therefore constant as time increases ($U^{ss}(0, 0, \tau) = 1$), but it extends over a smaller and smaller spatial region so that its contribution to physical quantities becomes eventually negligible at late time. Note finally that a peak exists only in the singlet density, or in $D_s$. As already emphasized, the total density of $Q\bar{Q}$ pairs is independent of $r$ and presents no peak at small $r$. To illustrate this point, it is convenient to consider an alternative representation of the density matrix [14], namely $D_0 = \frac{1}{N_c^2}(D_s + (N_c^2 - 1)D_o)$, for which $D_0(s, s', t) = U^{00}(s, s', t)D_0(s, s', 0) + U^{08}(s, s', t)D_8(s, s', 0)$. The evolution superoperator $U^{00}(s, s', t)$ is plotted in the right panel of figure 36. As expected, it reveals a uniform ridge in the $r$ direction, with a non trivial structure only in the $y$ direction.

With the full non-diagonal density matrix, there are several interesting quantities that can be calculated. One is the average momentum squared of the full distribution of heavy quarks. This can be obtained from the relation

$$ \langle p^2 \rangle = - \partial_y^2 \int dr \left[ D_s \left( r + \frac{y}{2}, r - \frac{y}{2}, t \right) + (N_c^2 - 1) D_o \left( r + \frac{y}{2}, r - \frac{y}{2}, t \right) \right] \Big|_{y=0} \quad (F.20) $$

where both singlet and octet components are taken into account. It turns out that an accurate estimate of $\langle p^2 \rangle$ can be obtained from the leading order expansion of $D$ at small time, which are given in eqs. (F.8). The singlet contribution writes

$$ \frac{d\langle p^2 \rangle_s}{dt} = -\frac{\kappa}{2} \int dr \partial_\tau \partial_y^2 \left[ 1 - C_F(s^2 + s'^2)\tau \right] D_s(0) \Big|_{y=0} = \frac{C_F \kappa}{2}, \quad (F.21) $$

where we have used $s^2 + s'^2 = 2r^2 + y^2/2$ and $\int dr D_s(0)|_{y=0} = \text{tr} D_s(0) = 1$. As for the octet contribution, one has (with $ss' = r^2 - \frac{y^2}{2}$)

$$ \frac{d\langle p^2 \rangle_o}{dt} = -\frac{N_c^2 - 1}{N_c} \frac{\kappa}{2} \int dr \partial_\tau \partial_y^2 [ss'\tau] D_s(0) \Big|_{y=0} = \frac{C_F \kappa}{2}. \quad (F.22) $$

Adding both contributions leads to $\frac{d\langle p^2 \rangle}{dt} = C_F \kappa$, which is also the result found in the abelian case. The linear relation $\langle p^2 \rangle = C_F \kappa t = \hat{q}t$ remains approximately valid till late time. It represents one of the important effects of collisions, namely the diffusion in momentum space, also referred to as momentum broadening. Since the momentum $p$ is conjugate to the variable $y = s - s'$ under Fourier transform, this growth of $\langle p^2 \rangle$ is correlated with the squeezing of the non-diagonal matrix elements of $D$ in coordinate space, a phenomenon that we have associated with collisional decoherence.

A further manifestation of collisional decoherence is seen in the time dependence of the survival probability of initial states, as can be exemplified by taking, as initial state, the 1S state

$$ \psi_{1S}(s) = \frac{1}{\sqrt{\pi \sigma_{1S}}} \exp \left( -\frac{s^2}{2\sigma_{1S}^2} \right), \quad (F.23) $$

and let $D_{1S}(s, s') = \psi_{1S}(s)\psi_{1S}^*(s')$ be the density matrix associated to this initial state. We have, for QED,

$$ p_{1S}^{\text{QED}}(t) = \int ds \, s' \psi_{1S}^*(s) \langle s | D_{1S}^{\text{QED}}(t) | s' \rangle \psi_{1S}(s') = \int ds \, s' \rho_{1S}(s) e^{-\tau(s-s')^2} \rho_{1S}(s'), \quad (F.24) $$
where we have used eq. (F.13) for $\langle s | \mathcal{D}_{1S}^{\text{QED}}(t) | s' \rangle$. Note that the squeezing along the $y = s - s'$ direction caused by the exponential factor entails a power law decrease of $p_{1S}$ at late time, $p_0 \propto \int dr \rho_{1S}^2(r) \times \int dy e^{-\tau y^2} \propto \tau^{-1/2}$. The explicit evaluation yields

$$p_{1S}^{\text{QED}}(t) = \frac{1}{\sqrt{1 + C_F \kappa \sigma_{1S}^2 t}}. \qquad (\text{F.25})$$

Note that we can write $C_F \kappa \sigma_{1S}^2 t = \hat{q}t/(2\langle p^2 \rangle_{1S})$, with $\hat{q} = C_F \kappa$ and $\langle p^2 \rangle_{1S} = 1/(2\sigma_{1S}^2)$. Thus, in the abelian case, the decay of the survival probability $p_{1S}^{\text{QED}}(t)$ is tightly connected with momentum broadening. The characteristic time scale for the decay is $2\langle p^2 \rangle_0/\hat{q}$, and the “suppression” of the initial state becomes significant when the average momentum squared gained by collisions, $\hat{q}t$, becomes comparable to the intrinsic momentum squared in the initial wave function. The survival probability $p_{1S}^{\text{QED}}(t)$ is displayed as a function of time in the left panel of figure 37. For the parameters used, we have $\hat{q} \approx 0.13 \text{ GeV}^2/\text{fm}$ and $\langle p^2 \rangle_{1S} \approx 0.07 \text{ GeV}^2$. It follows that the linear behavior predicted by eq. (F.25) holds only for a time of the order of $t = \langle p^2 \rangle_{1S}/\hat{q} \approx 0.5 \text{ fm/c}$.

While in the abelian case, the decay of the survival probability arises solely from momentum space diffusion, or momentum broadening, in QCD, an additional factor comes from the color dynamics. In order to get a simple estimate of $p_{1S}(t)$ in the QCD case, we use the approximation (F.19) and get:

$$p_{1S}(t) \approx \frac{N_c^2 - 1}{N_c^2} \times \frac{1}{(1 + \bar{\tau})} + \frac{1}{N_c^2} \times \frac{1}{\sqrt{1 + \hat{\tau}}}, \qquad (\text{F.26})$$

where we have set $\hat{\tau} = \Gamma_{1S}t$, with $\Gamma_{1S} = C_F \kappa \langle r^2 \rangle_{1S}$, and $\bar{\tau} = \frac{1}{2}N_c \kappa \langle r^2 \rangle_{1S}t$. This estimate is compared to the exact result obtained with the full expression (F.16) in figure 37. As one can see, the agreement is quite good.<fnref n="39" /> figure 37 also shows that the time dependence of $p_{1S}$ follows that of the abelian case, with however the additional damping that comes from the singlet-octet transitions.

The early time behavior $p_{1S}(t) \approx 1 - \Gamma_{1S}t$ coincides with the leading expression of the singlet weight, eq. (F.9). This is because, at early time, the density matrix is dominated by the initial singlet state. Similar results hold for other choices of initial states. Let us consider the states 1P and 2S, with the following harmonic wave functions

$$\psi_{1P}(s) = \frac{\sqrt{2}}{\pi^{1/4}} \sigma_{1S}^{3/2} s e^{-\frac{s^2}{2\sigma_{1S}^2}}, \qquad \psi_{2S}(s) = \frac{1}{\sqrt{2\sqrt{\pi}\sigma_{1S}}} \left(1 - \frac{2s^2}{\sigma_{1S}^2}\right) e^{-\frac{s^2}{2\sigma_{1S}^2}}. \qquad (\text{F.28})$$

$$p_{1S}(t) \approx 1 - \Gamma_{1S}t \left(1 + \frac{1}{2N_c^2}\right), \qquad (\text{F.27})$$

whereas the correct behavior is $p_{1S}(t) \approx 1 - \Gamma_{1S}t$, as can be easily deduced from the exact expression at linear order in $\tau$. The additional term in eq. (F.27) is a spurious artifact of the approximation (F.19). The comparison with the exact calculation shown in figure 37 indicates that this spurious contribution is of little quantitative consequence.

<footnote n="39">³⁹As mentioned earlier, a drawback of the approximation (F.19) is that it does not treat accurately the color dynamics at early time and in particular it does not yield the correct $N_c$ dependence of the small time expansion of the survival probability $p_{1S}$. We have indeed from eq. (F.26)</footnote>
<figure><img src="image_40.png" /><figcaption>**Figure 37.** Left: same as figure 6 (left panel) for the evolution performed solely with $\mathcal{L}_2$. $\sigma_{1S}$ has been taken as 0.55 fm while $\kappa \approx 500$ MeV fm$^{-2}$ has been adjusted to get $\Gamma_{1S} = 100$ MeV $= \Gamma_{1S}^{\text{in-medium}}$ for the “exact” decay width at $T = 300$ MeV. The dashed black line indicates the initial linear decay. The thick-solid orange line represents the QCD survival probability obtained numerically, and the dashed orange line is the approximation eq. (F.26), while the blue dot-dashed line illustrates the equivalent abelian law (F.25). Right: survival probabilities of the 1S, 1P and 2S states, calculated from eqs. (F.26) and corresponding laws for excited states; same conditions. The dashed lines represent the linear approximation $1 - \Gamma_{nS/P} t$.</figcaption></figure>

These wave functions are orthogonal and also orthogonal to the 1S wave function (F.23). It is easy to verify, using eq. (F.9) that the early time linear behavior of $p_{2S}$ and $p_{1P}$ is the same as that of $p_{1S}$, once the proper sizes ($\langle r^2 \rangle_{2S} = 5\langle r^2 \rangle_{1S}$ and $\langle r^2 \rangle_{1P} = 3\langle r^2 \rangle_{1S}$) are used in the respective expressions of the decay constant, i.e. $\Gamma_{1P} = 3\Gamma_{1S}$ and $\Gamma_{2S} = 5\Gamma_{1S}$. As can be seen in figure 37, this expected linear behavior is well reproduced, while at late time $p_{2S}$ and $p_{1P}$ behave similarly as $p_{1S} \propto t^{-1/2}$.

We have seen that at large time, the density matrix becomes rapidly diagonal in coordinate space (see eq. (F.13)), a property that we have associated to collisional decoherence. To emphasize the dependence on the choice of basis, and the fact that quantum coherence may survive in other basis, we return to the case of an initial wave packet considered in section 3.3 and write the corresponding initial singlet density matrix $D_s^\sigma(0) = |G_\sigma\rangle\langle G_\sigma|$ where $\sigma$ characterizes the width of the Gaussian. Depending on the value of $\sigma$, the probabilities $p_{1S}$ and $p_{2S}$ have different values at $t = 0$. These are given explicitly by

$$p_{1S}(\sigma; t=0) = \frac{2\sigma\sigma_{1S}}{\sigma^2 + \sigma_{1S}^2}, \qquad p_{2S}(\sigma; t=0) = \frac{\sigma\sigma_{1S}(\sigma_{1S}^2 - \sigma^2)^2}{(\sigma_{1S}^2 + \sigma^2)^3}, \qquad (F.29)$$

and their dependence on $\sigma$ is illustrated in figure 38. Note that $p_{2S}$ vanishes for $\sigma = \sigma_{1S}$. In this case, $|G_\sigma\rangle = |1S\rangle$ is orthogonal to $|2S\rangle$. We shall be interested in the range of $\sigma$ values delineated in figure 38 by $\sigma_{\text{compact}}$ and $\sigma_{\text{in-medium}}$. In this range, the initial density matrix is well accounted for by its matrix elements in the subspace spanned by the $|1S\rangle$ and $|2S\rangle$ states. For instance, in the range considered, $p_{1S} + p_{2S} \gtrsim 0.9$ (and equals 1 for $\sigma = \sigma_{1S}$), while the non diagonal element $\langle 2S|D_s^\sigma(t=0)|1S\rangle$ is given by

$$\langle 2S | D_s^\sigma(0) | 1S \rangle = \frac{\sqrt{2\sigma\sigma_{1S}(\sigma_{1S}^2 - \sigma^2)}}{(\sigma_{1S}^2 + \sigma^2)^2}. \qquad (F.30)$$
<figure><img src="image_41.png" /><figcaption>**Figure 38.** Initial 1S and 2S weights in a Gaussian wave packet — eq. (F.6) — as a function of the rescaled width $\sigma/\sigma_{1S}$. The 3 vertical lines illustrate the 3 initial states studied in section 3.3: vacuum 1S ($\sigma = \sigma_{1S} \equiv 0.387$ fm), in-medium 1S state ($\sigma \equiv 0.55$ fm), and compact Gaussian state ($\sigma = 0.165$ fm).</figcaption></figure>

It is then interesting to ask to which extent the evolution of the density matrix can be captured by that of these two states.

The evolution of the initial density matrix is given by eq. (F.15), from which we can extract the probability $p_{2S}(t) = \int_{ss'} \psi_{2S}^*(s) D_s^\sigma(s, s', t) \psi_{2S}(s')$. To within the normalization constants we get

$$p_{2S}(t) \propto \iint ds ds' e^{-\left(\frac{1}{2\sigma^2} + \frac{1}{2\sigma_{1S}^2}\right)(s^2 + s'^2)} \left[ \underbrace{1 - \frac{2(s^2 + s'^2)}{\sigma^2}}_{B} + \underbrace{\frac{4s^2 s'^2}{\sigma^4}}_{C} \right] U^{ss}(s, s', \tau). \quad (F.31)$$

As a result of the evolution induced by $U^{ss}(s, s', \tau)$, both $\langle (s - s')^2 \rangle$ and $\langle (s + s')^2 \rangle$ shrink with time, producing a differential evolution of the three contributions labeled A, B and C in eq. (F.31). Depending on the value of $\sigma$ this can lead to a reduction or to an increase of $p_{2S}$. Such an effect can be anticipated from figure 38, by noticing that the shrinking of $\langle (s - s')^2 \rangle$ and $\langle (s + s')^2 \rangle$ with time can be seen as an effective decrease of $\sigma$, producing a decrease of $p_{2S}$ if initially $\sigma > \sigma_{1S}$ and an increase in the opposite case.

Further insight can be gained by focusing on the short time evolution. At order $O(t)$, using eq. (F.8), one obtains

$$p_{1S}(\sigma; t) = \frac{2\sigma\sigma_{1S}}{\sigma^2 + \sigma_{1S}^2} \left(1 - C_F \frac{\sigma^2\sigma_{1S}^2}{\sigma^2 + \sigma_{1S}^2} \kappa t\right). \quad (F.32)$$

Thus, whatever the value of $\sigma$, $p_{1S}$ decreases linearly at small time. Note however that the decay constant is the harmonic average of the decay constants $\Gamma_{1S}$ and $\Gamma_\sigma$. It reduces to $\Gamma_{1S}$ only when $\sigma \to \sigma_{1S}$. For the 2S, we get

$$p_{2S}(\sigma; t) = \frac{\sigma\sigma_{1S}}{(\sigma_{1S}^2 + \sigma^2)^3} \left[ (\sigma_{1S}^2 - \sigma^2)^2 - C_F \frac{\sigma^2\sigma_{1S}^2(\sigma_{1S}^2 - \sigma^2)(\sigma_{1S}^2 - 5\sigma^2)}{\sigma_{1S}^2 + \sigma^2} \kappa t + \dots \right]. \quad (F.33)$$

In the relevant range of $\sigma$-values, the time derivative of $p_{2S}$ is negative for $\sigma \ge \sigma_{1S}$ and positive for $\sigma \le \sigma_{1S}$, confirming the previous heuristic argument. This result can also be
understood from a different perspective, relying on the expression of the density matrix on the {|1S⟩, |2S⟩} subspace. By projecting the density matrix at any time t on the basis {|1S⟩⟨1S|, |2S⟩⟨2S|, |1S⟩⟨2S|, |2S⟩⟨1S|}, one can then follow the time evolution of corresponding matrix elements. At leading order in t the expression of $p_{1S}$ is reproduced, while for $p_{2S}$ one finds

$$ 
\begin{aligned}
p_{2S}(\sigma; t) &= p_{2S}(\sigma; 0) - C_F \frac{\sigma \sigma_{1S}^3 (\sigma_{1S}^2 - \sigma^2) (3\sigma_{1S}^2 - 7\sigma^2)}{2 (\sigma^2 + \sigma_{1S}^2)^3} \kappa t \\
&\approx p_{2S}(\sigma; 0) + C_F \frac{\sigma_{1S}^2 - \sigma^2}{4} \kappa t
\end{aligned} 
\qquad (F.34) 
$$

This relation is not the same as (F.33). However, both relations yield the same time dependence when $\sigma \simeq \sigma_{1S}$, as indicated in the second line of eq. (F.34) above. The difference can be attributed to the fact that the initial wave packet $|G_\sigma\rangle$ contains an admixture of a $|3S\rangle$ state, and this can feed the 2S (but not the 1S). This admixture does not affect $p_{1S}$, and is negligible when $\sigma \simeq \sigma_{1S}$. Once this is included, one indeed reproduces (F.33).

Consider now the non diagonal matrix element $\langle 2S|D_s(t)|1S\rangle$, and choose as initial state $D_s(0) = |1S\rangle\langle 1S|$. The density matrix at any time t is then obtained from the approximate expression (F.19). The resulting expression is given by

$$ \langle 2S | D_s^{(1S)}(t) | 1S \rangle = \frac{N_c^2 - 1}{N_c^2} \frac{\bar{\tau}/\sqrt{2}}{(1+\bar{\tau})^2} + \frac{1}{N_c^2} \frac{\hat{\tau}/\sqrt{2}}{2(1+\hat{\tau})^{3/2}}, \qquad (F.35) $$

where the “(1S)” superscript recalls the initial state. The expression (F.35) shows that the matrix element $\langle 2S|D_s^{(1S)}(t)|1S\rangle$ evolves in time with the same time scales as the 1S-survival probability. In particular, it decreases at late time as $\tau^{-1/2}$, what is in fact true irrespective of the initial condition. This behavior is illustrated on figure 39. The comparison with the QED calculation shown in the right panel indicates that the same basic mechanisms are at work in the non-abelian case. The slow power law decay shows that in the basis of the vacuum eigenstates, effects due to quantum coherence can survive on a long time scale. This is to be contrasted with the exponential squeezing of the non-diagonal elements of the density matrix in coordinate space.

## G Equilibrium limit of the momentum distribution

In this appendix, we study the equilibrium limit of the momentum distribution. To keep the presentation simple, we consider the abelian case. We emphasize however that identical results are obtained in QCD (as long as the energy gap between bound states can be ignored). We assume that once equilibrium is reached (and the temperature is large enough), the bound states play a minor role, the quark and the antiquark occupying continuum states. In this regime, one expects a semi-classical approximation to be valid, and accordingly we treat the quantity $y = s' - s \sim \lambda_{th} = \frac{1}{\sqrt{MT}}$ as being small. Going back to the various Lindblad superoperators in appendix C, one then expands those of the form $\tilde{W}^{(n)}(\frac{y}{2})$ in powers of y, while others such as $\tilde{W}^{(n)}(s)$, $\tilde{W}^{(n)}(s')$ or $\tilde{W}^{(n)}(\frac{s+s'}{2})$, where $s \sim s' \sim r \gg \lambda_{th}$ are typically suppressed by powers of $\lambda_{th}/r$ and can be neglected. Besides, we only retain the derivatives of $D$ with respect to $y$ since $D$ is quasi uniform along $r \gg \lambda_{th}$ once equilibrium
<figure><img src="image_42.png" /><figcaption>Figure 39. Left: 1S survival eq. (F.26) as a function of the dimensionless time $\hat{\tau} = \Gamma_{1S}t$ (solid red), and the matrix element $\langle 2S|\mathcal{D}_s(t)|1S\rangle$, for the initial conditions $\mathcal{D}_s(t=0) = |1S\rangle\langle 1S|$ (long-dashed blue). Right: same for the QED case (which can be easily obtained by setting $N_c = 1$ in eq. (F.35)).</figcaption></figure>

is reached. The following expressions for the minimal set of Lindblad superoperators are then obtained, up to quadratic terms in $y$:

$$ \mathcal{L}_2 \approx -\frac{\tilde{W}''(0)}{4} y^2 \mathcal{D}, \qquad (G.1) $$

$$ \mathcal{L}_3 \approx -\frac{\tilde{W}''(0)}{2MT} y \partial_y \mathcal{D}, \qquad (G.2) $$

$$ \mathcal{L}_4 \approx -\frac{\tilde{W}^{(4)}(0)}{64M^2T^2} y^2 \partial_y^2 \mathcal{D}. \qquad (G.3) $$

If one solely considers the effect of $\mathcal{L}_2$ and $\mathcal{L}_3$, the unique non trivial solution of $(\mathcal{L}_2 + \mathcal{L}_3)\mathcal{D}_{asymp} = 0$ is $\mathcal{D}_{asymp} \propto \exp(-\frac{MTy^2}{4})$, from which on deduces that the Wigner transform of $\mathcal{D}_{asymp}$ is proportional to $\exp(-\frac{p^2}{MT})$, with the expected average $\langle p^2 \rangle = \frac{MT}{2}$. By also including the contribution of $\mathcal{L}_4$, one finds that the equation $(\mathcal{L}_2 + \mathcal{L}_3 + \mathcal{L}_4)\mathcal{D}_{asymp} = 0$ is equivalent to a modified Bessel equation. This leads to a more complicated solution for the Wigner transform, implying deviations from the Maxwell-Boltzmann distribution at large momenta. Note however that for the mere evaluation of $\langle p^2 \rangle$, knowing the small-$y$ behavior is enough since $\langle p^2 \rangle = -\frac{\mathcal{D}''(0)}{\mathcal{D}(0)}$. Writing $\mathcal{D}_{asymp}(0) = 1 - \frac{a_2}{2}y^2 + o(y^2)$ and injecting this expression in the $(\mathcal{L}_2 + \mathcal{L}_3 + \mathcal{L}_4)\mathcal{D}_{asymp} = 0$ equation, one obtains

$$ a_2 = \frac{MT}{2+\gamma} \quad \text{with} \quad \gamma = \frac{\tilde{W}^{(4)}(0)}{16MT\tilde{W}''(0)}. \qquad (G.4) $$

We have checked that these deviations from the Maxwell-Boltzmann distribution, represented by the factor $\gamma$ in the above equation,<fnref n="40" /> are not due to the restriction to the minimal set of superoperators that we have used, but they hold more broadly for the QME used in this paper and in [14].<fnref n="41" />

<footnote n="40">⁴⁰We may also interpret this correction as an effective modification of the temperature, $T \mapsto T_{\text{eff}} = 2T/(2+\gamma)$.</footnote>

<footnote n="41">⁴¹Such deviations are indeed also present in the equilibrium limit of the Lindblad equation with the $\mathcal{L}_2$ and $\mathcal{L}_3$ superoperators defined in appendix B of [14], where the action of $\mathcal{L}_3$ on $\mathcal{D}$ is evaluated as $\approx -\frac{1}{2MT} (W''(0)y\partial_y\mathcal{D} + \frac{W^{(4)}(0)}{8}y^2\mathcal{D})$, with an extra term $\propto W^{(4)}(0)$ as compared to eq. (G.2).</footnote>