Fate of in-medium heavy quarks via a Lindblad equation

TL;DR

This paper models the quantum dynamics of heavy quarks in a quark-gluon plasma using a Lindblad equation, exploring dissociation, recombination, and decoherence phenomena.

Contribution

It derives Lindblad equations for heavy quarks from gauge theory and performs simulations to study their in-medium behavior.

Findings

Bound-state dissociation depends on the imaginary part of the potential.

Recombination and decoherence are influenced by medium interactions.

Simulations reveal the importance of quantum effects in quark dynamics.

Abstract

What is the dynamics of heavy quarks and antiquarks in a quark gluon plasma? Can heavy-quark bound states dissociate? Can they (re)combine? These questions are addressed by investigating a Lindblad equation that describes the quantum dynamics of the heavy quarks in a medium. The Lindblad equations for a heavy quark and a heavy quark-antiquark pair are derived from the gauge theory, following a chain of well-defined approximations. In this work the case of an abelian plasma has been considered, but the extension to the non-abelian case is feasible. A one-dimensional simulation of the Lindblad equation is performed to extract information about bound-state dissociation, recombination and quantum decoherence for a heavy quark-antiquark pair. All these phenomena are found to depend strongly on the imaginary part of the inter-quark potential.

Equations

$${\hat{\rho }}_{\mathrm{tot}}(t_0)=\hat{\rho }(t_0)\otimes {\hat{\rho }}_{\mathrm{env}}(t_0),$$

$${\hat{\rho }}_{\mathrm{tot}}(t)={\mathrm{e}}^{-i{\hat{H}}_{\mathrm{tot}}(t-t_0)}{\hat{\rho }}_{\mathrm{tot}}(t_0){\mathrm{e}}^{i{\hat{H}}_{\mathrm{tot}}(t-t_0)},$$

$${\hat{H}}_{\mathrm{tot}}=\hat{H}\otimes {\mathbb{I}}_{\mathrm{env}}+\mathbb{I}\otimes {\hat{H}}_{\mathrm{env}}+{\hat{H}}_{\mathrm{int}},$$

$$P(\psi ,t|{\psi }_0,t_0)$$

$$P(\psi ,t|{\psi }_0,t_0)=\text{mbox}{Tr}_{\mathrm{env}}⟨\psi |{\mathrm{e}}^{-i{\hat{H}}_{\mathrm{tot}}(t-t_0)}{\hat{\rho }}_{\mathrm{tot}}(t_0){\mathrm{e}}^{i{\hat{H}}_{\mathrm{tot}}(t-t_0)}|\psi ⟩$$

$$P(q,q^{\prime },t|q_0,q_0^{\prime },t_0)=\text{mbox}{Tr}_{\mathrm{env}}\left\{⟨q|{\mathrm{e}}^{-i{\hat{H}}_{\mathrm{tot}}(t-t_0)}{\hat{\rho }}_{\mathrm{env}}(t_0)|q_0⟩⟨q_0^{\prime }|{\mathrm{e}}^{i{\hat{H}}_{\mathrm{tot}}(t-t_0)}|q^{\prime }⟩\right\}.$$

$$\rho (t,q,q^{\prime })=\int \mathrm{d}{q}_0\int \mathrm{d}{q}_0^{\prime }P(q,q^{\prime },t|q_0,q_0^{\prime },t_0)\rho (t_0,q_0,q_0^{\prime }).$$

$$\mathrm{i}\hslash \frac{d{\hat{\rho }}_{\mathrm{tot}}}{dt}(t)=\left[{\hat{H}}_{\mathrm{tot}},{\hat{\rho }}_{\mathrm{tot}}(t)\right],$$

$$\mathrm{i}\hslash \frac{d\hat{\rho }}{dt}(t)$$

$$\dot{\hat{\rho }}=-\frac{\mathrm{i}}{\hslash }[\hat{H},\hat{\rho }]+\frac{1}{2\hslash }\sum _{\mu }\left([{\hat{L}}_{\mu }\hat{\rho },{\hat{L}}_{\mu }^{†}]+[{\hat{L}}_{\mu },\hat{\rho }{\hat{L}}_{\mu }^{†}]\right),$$

$$\dot{\hat{\rho }}=\frac{\mathrm{i}}{\hslash }\left(-{\hat{H}}_{\mathrm{eff}}\hat{\rho }+\hat{\rho }{\hat{H}}_{\mathrm{eff}}^{†}-\mathrm{i}\sum _{\mu }{\hat{L}}_{\mu }\hat{\rho }{\hat{L}}_{\mu }^{†}\right),$$

$$⟨q|\hat{\rho }(t+\Delta t)|q^{\prime }⟩={\int }_{q_0}{\int }_{q_0^{\prime }}⟨q_0|\hat{\rho }(t)|q_0^{\prime }⟩(⟨q|q_0⟩⟨q_0^{\prime }|q^{\prime }⟩+\frac{\Delta t}{\hslash }(-\mathrm{i}⟨q|{\hat{H}}_{\mathrm{eff}}|q_0⟩⟨q_0^{\prime }|q^{\prime }⟩$$

$$+\mathrm{i}⟨q|q_0⟩⟨q_0^{\prime }|{\hat{H}}_{\mathrm{eff}}^{†}|q^{\prime }⟩+\sum _{\mu }⟨q|{\hat{L}}_{\mu }|q_0⟩⟨q_0^{\prime }|{\hat{L}}_{\mu }^{†}|q^{\prime }⟩))$$

$$⟨q|q_0⟩={\int }_p{\mathrm{e}}^{\frac{i}{\hslash }p(q-q_0)},{\int }_p\equiv \frac{1}{(2\pi \hslash {)}^d}\int \mathrm{d}p,$$

$$O(q,q^{\prime })={\int }_{p}O\left(\frac{q+q^{\prime }}{2},p\right){\mathrm{e}}^{\frac{i}{\hslash }p(q-q^{\prime })},$$

$$O(q,p)\equiv {\int }_{q^{\prime }}⟨q-\frac{q^{\prime }}{2}|\hat{O}|q+\frac{q^{\prime }}{2}⟩{\mathrm{e}}^{\frac{i}{\hslash }p{q}^{\prime }}.$$

$${\int }_{q_0}{\int }_{q_0^{\prime }}{\int }_p{\int }_{p^{\prime }}⟨q_0|\hat{\rho }(t)|q_0^{\prime }⟩{\mathrm{e}}^{\frac{\mathrm{i}}{\hslash }\left(p(q-q_0)-p^{\prime }(q^{\prime }-q_0^{\prime })\right)}[1+\frac{\mathrm{i}\Delta t}{\hslash }\cdot$$

$$\cdot \left(-H_{\mathrm{eff}}\left(\bar{q},p\right)+H_{\mathrm{eff}}^{\ast }\left({\bar{q}}^{\prime },p^{\prime }\right)-\mathrm{i}\sum _{\mu }{L}_{\mu }\left(\bar{q},p\right)L_{\mu }^{\ast }\left({\bar{q}}^{\prime },p^{\prime }\right)\right)],$$

$$⟨q|\hat{\rho }(t+\Delta t)|q^{\prime }⟩={\int }_{q_0}{\int }_{q_0^{\prime }}{\int }_p{\int }_{p^{\prime }}⟨q_0|\hat{\rho }(t)|q_0^{\prime }⟩\exp [\frac{\mathrm{i}}{\hslash }p(q-q_0)-\frac{\mathrm{i}\Delta t}{\hslash }{H}_{\mathrm{eff}}\left(\bar{q},p\right)$$

$$-\frac{\mathrm{i}}{\hslash }{p}^{\prime }(q^{\prime }-q_0^{\prime })-\frac{\mathrm{i}\Delta t}{\hslash }{H}_{\mathrm{eff}}^{\ast }\left({\bar{q}}^{\prime },p^{\prime }\right)+\frac{\Delta t}{\hslash }\sum _{\mu }{L}_{\mu }\left(\bar{q},p\right)L_{\mu }^{\ast }\left({\bar{q}}^{\prime },p^{\prime }\right)].$$

$$\rho (t,q,q^{\prime })={\int }_{q_0}{\int }_{q_0^{\prime }}P(q,q^{\prime },t|q_0,q_0^{\prime },t_0)\rho (t_0,q_0,q_0^{\prime }),$$

$$P(q,q^{\prime },t|q_0,q_0^{\prime },t_0)=\underset{N\to \infty }{\lim }{\int }_{p_1}{\int }_{q_1}\cdots {\int }_{q_{N-1}}{\int }_{p_N}{\int }_{p_1^{\prime }}{\int }_{q_1^{\prime }}\cdots {\int }_{q_{N-1}^{\prime }}{\int }_{p_N^{\prime }}$$

$$\exp \left[\frac{\mathrm{i}\Delta t}{\hslash }\sum _{i=1}^N\left(p_i{\dot{q}}_i-H_{\mathrm{eff}}({\bar{q}}_i,p_i)-p_i^{\prime }{\dot{q}}_i^{\prime }+H_{\mathrm{eff}}^{\ast }({\bar{q}}_i^{\prime },p_i^{\prime })-\mathrm{i}\sum _{\mu }{L}_{\mu }({\bar{q}}_i,p_i)L_{\mu }^{\ast }({\bar{q}}_i^{\prime },p_i^{\prime })\right)\right],$$

$$P(q_f,q_f^{\prime },t|q_0,q_0^{\prime },t_0)={\int }_{(q_0,t_0)}^{(q_f,t)}\mathcal{D}q\mathcal{D}p{\int }_{(q_0^{\prime },t_0)}^{(q_f^{\prime },t)}\mathcal{D}{q}^{\prime }\mathcal{D}{p}^{\prime }\exp \left[\frac{\mathrm{i}}{\hslash }S[q,p;q^{\prime },p^{\prime }]\right],$$

$$S[q,p;q^{\prime },p^{\prime }]=$$

$$\delta S[q,p;q^{\prime },p^{\prime }]=0$$

$$\dot{q}=\frac{\partial H_{\mathrm{eff}}(q,p)}{\partial p}+\mathrm{i}\sum _{\mu }\frac{\partial L_{\mu }(q,p)}{\partial p}{L}_{\mu }^{\ast }(q^{\prime },p^{\prime }),$$

$$\dot{p}=-\frac{\partial H_{\mathrm{eff}}(q,p)}{\partial q}-\mathrm{i}\sum _{\mu }\frac{\partial L_{\mu }(q,p)}{\partial q}{L}_{\mu }^{\ast }(q^{\prime },p^{\prime }),$$

$${\dot{q}}^{\prime }=\frac{\partial H_{\mathrm{eff}}^{\ast }(q^{\prime },p^{\prime })}{\partial p^{\prime }}-\mathrm{i}\sum _{\mu }{L}_{\mu }(q,p)\frac{\partial L_{\mu }^{\ast }(q^{\prime },p^{\prime })}{\partial p^{\prime }},$$

$${\dot{p}}^{\prime }=-\frac{\partial H_{\mathrm{eff}}^{\ast }(q^{\prime },p^{\prime })}{\partial q^{\prime }}+\mathrm{i}\sum _{\mu }{L}_{\mu }(q,p)\frac{\partial L_{\mu }^{\ast }(q^{\prime },p^{\prime })}{\partial q^{\prime }},$$

$$q(t_0)=q_0,q^{\prime }(t_0)=q_0^{\prime },q(t)=q_f,q^{\prime }(t)=q_f^{\prime },$$

$$H_{\mathrm{eff}}(q,p)=H(q,p)-\frac{\mathrm{i}}{2}\sum _{\mu }|{\Lambda }_{\mu }(q,p){|}^2,$$

$$S_{\mathrm{cl}}[q,p]={\int }_{t_0}^t\mathrm{d}\tau \left(\dot{q}p-H(q,p)\right),$$

$$\mathcal{F}[q,p;q^{\prime },p^{\prime }]=\exp \left[\frac{\mathrm{i}}{\hslash }\phi [q,p;q^{\prime },p^{\prime }]-\frac{1}{2\hslash }{D}^2[q,p;q^{\prime },p^{\prime }]\right],$$

Full text

††arxiv: 1705.03567

1 Introduction

Experimental and theoretical physicists have long sought to understand the properties of the elementary particles (quarks and gluons) described by QCD at extreme temperatures and densities. Experimentalists are creating these extreme conditions by means of relativistic heavy-ion collisions performed at LHC and RHIC. During these collisions, a fireball of deconfined quarks and gluons, known as the quark gluon plasma (QGP), is formed for a very brief time, which makes the investigation of the fireball a challenging task. Joint experimental and theoretical efforts have been carried out for about fourty years Shuryak:1978ij to find and analyse accessible observables for diagnosing this short-lived QGP.

One of the promising candidates for diagnostics are heavy quarks produced and partecipating in the collisions. Due to the separation of scales between the heavy quarks and the main plasma constituents, the former can be used as dynamical probes of the surrounding medium. In particular, the seminal work of Matsui and Satz Matsui:1986dk suggested that the analysis of bound states of heavy quarks (quarkonia) could lead to a better understanding of the QGP, the reason being that, at sufficiently high temperatures, colour screening weakens or even prevents a $q\bar{q}$ pair from binding in the deconfined medium, resulting in a suppression of quarkonia (such as $c\bar{c}$ and $b\bar{b}$ bound states) at the hadronization of the QGP. This paved the way for many quests for understanding the fate of quarkonia.

Some of them made use of screened potential models Mocsy:2013syh , based on quantities like the free or internal energies. The screened potential, whose short- and medium-range parts can be obtained from resummed perturbative calculations, was also used to determine stationary states of a Schrödinger equation Karsch1988 . Non-perturbative studies on the lattice followed a different line aiming at reconstructing spectral functions of quarkonia around the crossover temperature Petreczky:2012rq ; Aarts:2011sm ; Aarts:2013kaa ; Aarts:2014cda , whereas other approaches exploit the AdS/CFT correspondence PhysRevD.74.085012 ; Herzog:2006gh ; Albacete:2008dz ; Hayata:2012rw . Recent progresses showed that the inter-quark potential not only gets screened at high temperatures, but develops an imaginary part Laine:2006ns ; Burnier:2007qm , which reflects the Landau damping mechanism due to the collisions between the heavy quarks and plasma constituents Beraudo:2007ky ; Brambilla:2008cx ; Brambilla:2010vq . The fact that the inter-quark potential is actually complex poses serious issues on the validity of the approaches relying on the standard Schrödinger equation, which is a good approximation only when the state lives long enough, that is when the deviation from the unitary evolution is small. This happens when the imaginary part of the potential (decay rate) is much smaller than the mass and binding energy of respectively a heavy quark and a heavy-quark bound state111Moreover the Schrödinger equation can not describe possible bound-state recombinations.. However, the imaginary part of the potential can be used to compute spectral functions of in-medium quarkonia or in the context of the stochastic Schrödinger equation (SSE) PhysRevD.13.857 ; Gisin1989 ; 0305-4470-26-9-019 . Both the real and the imaginary part of the potential have been calculated in perturbation theory, but there have also been attempts to get these quantities using lattice Rothkopf:2011db and holographic Noronha:2009da ; Giataganas:2013lga ; Fadafan:2013coa ; 0954-3899-43-9-095001 techniques. Some reviews on recent developments on the study of heavy quarks in a quark gluon plasma can be found in Brambilla:2010cs ; Rapp:2009my ; CasalderreySolana:2011us ; Aarts:2016hap .

The fact that the inter-quark potential is screened and features an imaginary part fosters the phenomenon of quarkonia suppression, which has indeed been partially detected for $J/\Psi$ (see reviews Rapp:2008tf ; Kluberg:2009wc ) and $\Upsilon$ states Khachatryan:2016xxp in heavy-ion collision experiments. On the other hand, evidence of recombination of bound states has also been observed for $J/\Psi$ Abelev:2013ila (the phenomenon is much less important for $\Upsilon$ CMS:2017ucd ), hindering the opposite mechanism of dissociation. In fact it is not surprising that some of the heavy particles liberated from the bound states can recombine to form a bound state before the freeze-out of the QGP, especially when there are many heavy quarks and antiquarks tossing about in the hot medium.

However, in order to describe the recombination process, a dynamical treatment of the heavy particles becomes necessary, and a static description of the heavy quarks, either via lattice techniques in imaginary time or potential models, is not appropriate anymore. Formulating it as a real-time problem of probes immersed in an environment, the correct language to use is then the one of open quantum systems.

The goal of this manuscript is to provide a unifying framework to study both dissociation and recombination mechanisms within the language of open quantum systems. The main idea is to exploit the weak-coupling expansion and the separation of scales between heavy probes and environment to construct an effective theory from the underlying gauge theory, after integrating out the plasma degrees of freedom. A master equation of the Lindblad type for the density matrix of the heavy particles is derived, in which all the quantities appearing are obtained from the complex potential. This equation allows one to study numerically the real-time dynamics of the heavy quarks and antiquarks.

Similar strategies based on a Lindblad equation have been used before using a different formalism with Brambilla:2016wgg and without Akamatsu:2012vt ; Akamatsu:2014qsa numerical computations. A one-dimensional numerical simulation has been performed for the master equation and the related SSE in Akamatsu:2011se , but the master equation was not directly obtained from the gauge theory. The SSE for in-medium quarkonium has been investigated in Akamatsu:2011se ; Rothkopf:2013ria ; Rothkopf:2013kya ; Kajimoto:2017rel , but there seems to be no clear way to incorporate dissipative effects (equivalent to the classical friction) within this formalism. The Schrödinger-Langevin equation (SLE) Katz:2015qja ; Katz:2015edb includes dissipation due to friction but it is not associated to a master equation (hence to a SSE). This means that it is not clear how to derive a SLE from the underlying theory. Another strategy for studying the dynamics of in-medium quarkonium exploited a classical Fokker-Planck and Langevin equation PhysRevC.79.034907 ; Young:2009tj ; Moore:2004tg ; Akamatsu:2015kaa . A generalised Langevin equation for many heavy particles, obtained from the underlying gauge theory, has been simulated in Blaizot:2015hya . The limitation of the Langevin equation is that it describes the classical dynamics of the heavy particles, but cannot address topics like quantum decoherence and quantum bound states.

This paper extends the semiclassical results of the previous work Blaizot:2015hya to a quantum level, by deriving a quantum master equation in the Lindblad form, which ensures positivity of the density operator, for the density matrix of heavy quarks and antiquarks. This equation allows one to follow the quantum dynamics of the heavy particles and compute survival and formation probabilities of heavy-quark bound states, as well as quantum decoherence.

The flow chart in Fig.1 serves as a guide to understanding how this paper is organised. Mimicking what has been done in 0305-4470-30-11-030 , section 2 shows how to build a propagator (or conditional probability), which is related to the influence functional, starting from a generic Lindblad equation. Section 3 introduces the model used here to study heavy particles in a thermal bath of light particles. The heavy particles are regarded as being nonrelativistic, whereas the light particles of the medium are treated relativistically. As long as only unbound heavy quarks are concerned, this means that their momenta have to be much smaller than their masses. Similarly, for quarkonia the nonrelativistic condition entails that their relative velocities are much smaller than the speed of light. This paper ignores the specifics of QCD interactions (like attractive/repulsive interactions according to the singlet/octet state for example) and retains only Coulomb interactions. Within this framework, the influence functional for a single heavy quark is obtained with the help of some well-defined approximations, such as the weak coupling and small frequency (Markovian limit) expansions. In section 4 the associated Lindblad equation is then read off from the influence functional. The Lindbladian terms will result in an expression written solely in terms of the imaginary part of the potential. From the structure of the Lindblad equation, it will be shown that different types of dynamics of the heavy particle are obtained according to the ratio between the correlation length of the environment and the typical size of the heavy probe. In particular, it will be proven that the classical Langevin dynamics stems from a semiclassical approximation of the Lindblad equation (see lower part of Fig.1). Section 5 contains the computation of the Lindblad equation for a heavy quark-antiquark pair, which could in principle be generalised to $N$ heavy quarks and antiquarks. Section 6 shows how quantum decoherence and the probability of having a $q\bar{q}$ bound state at a certain time can be calculated from the density matrix. This section also contains the setup for the one-dimensional numerical simulation of the Lindblad equation for the relative motion of the $q\bar{q}$ pair, after the motion of the centre of mass has been traced out. The results are displayed in section 7, with specific focus on the linear entropy (a measure of decoherence) and bound-state probabilities, which can be used to extract dissociation and recombination probabilities for different correlation lengths of the thermal bath. Conclusions and perspectives are found in section 8, followed by appendices A, B, C.

2 Density matrix and propagator of an open quantum system

The goal of this work is to study a probe propagating through a medium in thermal equilibrium. The probe is treated as an open quantum system (also called reduced system or subsystem) that can dissipate or gain energy from the plasma environment. The environment and the probe form a closed quantum system. A master equation for the density matrix describing the probe will be derived in order to study its real-time dynamics.

Of particular interest is the probability $P(\psi,t|\psi_{0},t_{0})$ to find the reduced system in a quantum state $|\psi\rangle$ at time $t$, given that it was in a state $|\psi_{0}\rangle$ at time $t_{0}$. This quantity can be written in terms of the density operator of the subsystem, $\hat{\rho}(t)\equiv\mbox{Tr}_{\rm env}\left\{\hat{\rho}_{\rm tot}(t)\right\}\,$, where $\hat{\rho}_{\rm tot}$ is the density operator of the total (closed) system and the trace is taken over the degrees of freedom of the environment. It is assumed that the subsystem has not yet interacted with the medium at time $t_{0}$, so that the initial total density operator takes the factorized form

[TABLE]

with $\hat{\rho}(t_{0})=|\psi_{0}\rangle\langle\psi_{0}|$ and $\hat{\rho}_{\rm env}(t_{0})=\frac{1}{Z_{\rm env}}{\rm e}^{-\beta\hat{H}_{\rm env}}$ is the density operator of the environment, whose Hamiltonian and partition function are respectively $\hat{H}_{\rm env}$ and $Z_{\rm env}=\mbox{Tr}_{\rm env}\exp(-\beta\hat{H}_{\rm env})$. The total density operator at time $t$ is given by

[TABLE]

with the Hamiltonian

[TABLE]

where $\hat{H}$ is the Hamiltonian of the isolated subsystem and $\hat{H}_{\rm int}$ accounts for the interaction between the subsystem and the medium. The desired probability can now be written as

[TABLE]

where $\left\{\phi_{\rm env}\right\}$ forms a basis of the Hilbert space of the medium and $\rho(t,q,q^{\prime})=\langle q|\hat{\rho}(t)|q^{\prime}\rangle$ is the density matrix of the subsystem. Using (2.2), this probability can equivalently be expressed as

[TABLE]

where the propagator (or conditional probability) has been defined as

[TABLE]

Equations (2.4) and (2.5) imply that the density matrix is related to the propagator through the relation

[TABLE]

In the next paragraph the path integral expression for the propagator, hence for the density matrix, will be derived in the Markovian limit through the well-known Trotter decomposition. In order to obtain the path integral, one needs to know the equations of motion satisfied by the reduced density operator. The density operator of the closed system obeys the Liouville-von Neumann equation

[TABLE]

which is equivalent to a Schrödinger equation, describing reversible quantum dynamics. Once the trace over the environmental degrees of freedom is taken, the reduced density operator does not obey the previous equation anymore, being the energy of the open quantum system not conserved.

Instead of eq.(2.8), one gets a more complicated equation for $\hat{\rho}(t)\,$, known as master equation:

[TABLE]

where $\hat{\mathcal{D}}$ is an operator (acting on the Hilbert space of the subsystem) that describes the irreversible part of the dynamics of the subsystem propagating through the medium.

2.1 Path integral solution of the Lindblad equation

Here a path integral expression is derived for the reduced density matrix in the Markovian limit, that is when memory effects of the medium can be largely neglected due to the fact that the intrinsic time scale of the subsystem is much larger than the typical bath correlation time. The path-integral construction mimics the one found in 0305-4470-30-11-030 . The same derivation is presented here for self-containment and to set up the language used throughout this work.

It has been shown CIS-338355 ; 8ab9645a29cc4e23a4837264fa594328 that the most general Markovian master equation can be written in the Lindblad form (also known as Bloch equation)

[TABLE]

or, equivalently,

[TABLE]

with $\hat{H}_{\rm eff}\equiv\hat{H}-\frac{\mathrm{i}}{2}\sum_{\mu}\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}$ and $\left\{\hat{L}_{\mu}\right\}_{\mu}$ forms an arbitrary linear basis of operators acting on the Hilbert space of the subsystem. Starting from the discretised version of the last expression, one finds a relation between $\hat{\rho}(t+\Delta t)$ and $\hat{\rho}(t)$ in the position basis, that is

[TABLE]

where the completeness relationship $\int_{q}|q\rangle\langle q|={\mathbbm{1}}$ has been inserted, together with the notation $\int_{q}\equiv\int\mathop{}\!\mathrm{d}{q}$. Using that

[TABLE]

and the identity

[TABLE]

where the Wigner transform is defined as

[TABLE]

The right hand side of eq.(2.1) then reads

[TABLE]

with the midpoints $\bar{q}\equiv\frac{q+q_{0}}{2}$ and $\bar{q}^{\prime}\equiv\frac{q^{\prime}+q_{0}^{\prime}}{2}$. Therefore, up to first order in $\Delta t\,$, eq.(2.1) re-exponentiated becomes

[TABLE]

The next step consists in deriving the final density matrix $\rho(t,q,q^{\prime})$ from the initial density matrix $\rho(t_{0},q,q^{\prime})\,$. To accomplish this, the time interval $[t_{0},t]$ is partitioned in $N$ small steps of size $\Delta t\,$. Then, after using eq.(2.1) iteratively, one is left with

[TABLE]

with the propagator

[TABLE]

and $(q_{N},p_{N})=(q,p)\,$, $(q_{N}^{\prime},p_{N}^{\prime})=(q^{\prime},p^{\prime})\,$. One should notice that the appearance of the midpoints $\bar{q}_{i}=\frac{q_{i}+q_{i-1}}{2}$ is a result of the use of the Wigner representation (2.15). On the other hand the velocities are defined through the retarded rule $\dot{q}_{i}=(q_{i}-q_{i-1})/\Delta t\,$.

In the continuous limit the propagator (2.1) can be formally written as

[TABLE]

with the phase space action functional

[TABLE]

The appearance of the primed and unprimed coordinates can be traced back to the way the density operator evolves with time (see eq.(2.2)), that is with a time evolution from $t_{0}$ to $t$ and one in the opposite direction. The unprimed coordinates live in the former interval, whereas the primed ones live in the latter branch. These two oriented time branches form the well-known Keldysh contour. Observe that the primed and unprimed coordinates decouple in the absence of Lindblad operators.

2.2 Equations of motion and quantum decoherence

When the semiclassical limit $\hbar\to 0$ of expression (2.19) is taken, the variational principle

[TABLE]

gives the complex equation of motions corresponding to the stationary paths of the subsystem, that is,

[TABLE]

with the boundary conditions

[TABLE]

The equations of motion for one heavy quark propagating through a medium will be derived in section 3.2. Notice the remarkable fact that dissipative dynamics can be derived from a variational principle. Following 0305-4470-30-11-030 , the propagator (2.18) can be rewritten in a more transparent way, which separates the reversible dynamics generated by $H(q,p)$ from the irreversible contribution coming from the Lindbladian matrix elements $L_{\mu}(q,p)\,$. Denoting the Wigner transform of $\hat{L}_{\mu}^{\dagger}\hat{L}_{\mu}$ as $|\Lambda_{\mu}(q,p)|^{2}\,$, which implies that

[TABLE]

the propagator can be expressed in terms of the classical action of the isolated subsystem,

[TABLE]

and the in-medium piece

[TABLE]

where

[TABLE]

is a phase functional modifying the classical action of the isolated system and

[TABLE]

is called decoherence distance functional, since it measures the decoherence between the phase space histories 0305-4470-30-11-030 . Notice that the phase functional $\varphi$ has to vanish when one takes $(q^{\prime},p^{\prime})=(q,p)\,$, and $D$ is symmetric under $(q^{\prime},p^{\prime})\leftrightarrow(q,p)$ exchange.

The propagator (2.18) can be finally written in terms of the functionals $\varphi$ and $D$ as

[TABLE]

In the following sections the phase functional $\varphi$ and the decoherence functional $D$ will be calculated first for one heavy quark and then for a heavy $q\bar{q}$ pair propagating through a quark gluon plasma.

3 Heavy quarks in a thermal medium: the model

The general setup used in this paper is the same discussed in the previous work Blaizot:2015hya . Only the key features of the model (such as physical scales and approximations) are presented here, whereas all the other details and calculations are found in Blaizot:2015hya .

In order to describe the real-time dynamics of heavy quarks and antiquarks propagating through a plasma of thermalised light particles, an effective theory for the heavy probes is obtained by tracing out the plasma degrees of freedom.

The model considered here consists of an abelian plasma of massless charged particles interacting with positive- and negative-charged heavy particles, which are called heavy quarks and antiquarks respectively. Strictly speaking, this picture applies only to electromagnetic interactions, being non abelian features of QCD not taken into account. However, the language of QCD has been used throughout this manuscript due to the similarities between the QED and QCD plasmas at very high temperatures. In this model $N$-nonrelativistic222Consequently their number is fixed. heavy quarks and antiquarks have been considered, whose mass $m$ is much larger than the temperature $T=1/\beta$ of the plasma ($k_{B}$ and $c$ have been set to $1$). The heavy particles are described using first quantisation, whereas the light particles of the thermal bath are considered within a quantum field theory framework. The heavy probes interact amongst themselves and with the plasma particles via Coulomb interactions. Magnetic interactions among the heavy probes are negligible in the non relativistic limit, being suppressed by powers of the velocity, or $p/m\,$. Magnetic interactions are also not taken into account for the light particles of the medium since the heavy quarks do not induce magnetic excitations. Therefore the whole system undergoes only Coulomb interactions. In the Coulomb gauge, the Hamiltonian of the total system reads

[TABLE]

where $\alpha^{i}=\gamma_{0}\gamma^{i}$ is a Dirac matrix, and $j_{\rm tot}^{0}=j^{0}+j_{\psi}^{0}$ is the total charge density, with

[TABLE]

the charge density of the heavy quarks and antiquarks, and

[TABLE]

the density of the charged light quarks and antiquarks of the plasma, where $g$ is the gauge coupling. The spatial coordinates ${\bm{q}}_{j}$ and $\bar{\bm{q}}_{j}$, with $j=1,\cdots,N$, refer to the heavy quarks and antiquarks respectively, and ${\bm{p}}_{j}$, $\bar{\bm{p}}_{j}$ are the corresponding momenta. The plasma is supposed to be electrically neutral, that is, it contains the same number of light quarks and antiquarks.

3.1 Influence functional for $N$ heavy quarks and antiquarks

In Blaizot:2015hya it has been shown in details how to get the path integral expression of the propagator (2.6) for the heavy particles, after integrating out the plasma degrees of freedom. The conditional probability can be written as

[TABLE]

where the $2N$-dimensional vector ${\bm{Q}}=\left({\bm{q}}_{1},\cdots,{\bm{q}}_{N},\bar{\bm{q}}_{1},\cdots,\bar{\bm{q}}_{N}\right)$ (and similarly for ${\bm{Q}}^{\prime}$) collectively denotes the spatial coordinates of all the heavy particles. The coordinate ${\bm{Q}}$ lives in the first (upper) branch of the Keldysh contour in the complex time plane. This branch runs from $\tau=0$ to $\tau=t\,$. The coordinate ${\bm{Q}}^{\prime}$ lives in the second (lower) branch, which runs backwards in time from $\tau=t$ to $\tau=0\,$. The free action of the heavy particles is

[TABLE]

The influence functional $\Phi$ is obtained after performing a perturbative expansion up to order $\alpha_{s}=g^{2}/(4\pi)$ (non-linear interactions are ignored) and considering almost instantaneous interactions between heavy and light particles, exploiting the fact that one is interested in studying the dynamics of the heavy probes over a time scale much larger than the typical one of the particles of the medium333Refer to the small frequency expansion of the correlation functions in Blaizot:2015hya for further details and the derivation of the influence functional.. The influence functional can be split into three terms as $\Phi=\Phi_{{}_{QQ}}+\Phi_{{}_{\bar{Q}\bar{Q}}}+\Phi_{{}_{Q\bar{Q}}}\,$, where

[TABLE]

and similarly for $\Phi_{{}_{\bar{Q}\bar{Q}}}$ with the substitution $\{{\bm{q}}_{i}\}\rightarrow\{\bar{\bm{q}}_{i}\}$. The mixed quark-antiquark contribution reads

[TABLE]

The functions $V$ and $W$ come from the correlation functions describing the interactions between the heavy probes amongst themselves and with the bath. They are respectively the real and the imaginary part of the complex potential experienced by the heavy particles, derived as

[TABLE]

where $x=|{\bm{x}}|\,$, $m_{{}_{\scriptsize D}}\propto gT$ is the Debye mass and $\Lambda$ is a momentum cutoff coming from the Hard Thermal Loop approximation (see discussion in Blaizot:2015hya ). In the next paragraph, all the terms appearing in the influence functional will be related to the Lindblad matrix elements introduced in section 2.2. This is the point where this project takes another path compared to the one followed in the previous work Blaizot:2015hya . In the latter, a semiclassical approximation of $\Phi$ was performed in order to get a generalised Langevin equation for $N$ heavy quarks and antiquarks. Instead, here a quantum master equation will be derived from the influence functional (see Fig.1). The master equation allows one to deal with quantum objects like bound states and study related quantum phenomena. These aspects could not be studied with a classical Langevin equation.

3.2 Conditional probability for a single heavy quark

For a single heavy quark in the plasma the influence functional simplifies to

[TABLE]

In order to calculate the phase functional $\varphi$ and the decoherence distance functional $D$ for a single heavy quark, the conditional probability (3.32) has to be rewritten in phase space as in (2.2). To achieve this, the following identity is useful:

[TABLE]

Now one defines the functions

[TABLE]

and neglect terms of order $g^{4}$ when computing $f^{2}$ and $h^{2}\,$. It is then easy to prove that the propagator can be expressed in phase space as in (2.2), where the classical action is

[TABLE]

and the in-medium contribution

[TABLE]

One can show that the Lindblad equation that generates the propagator (2.2) with (3.41) violates the conservation of probability ($\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}t}\mbox{Tr}\,\hat{\rho}(t)=0\,$) to order $g^{4}$. Hence, with hindsight, the following term, which is negligible both in perturbation theory and in the non relativistic limit ($T/m\ll 1\,$), is added ad hoc to the in-medium part $\mathcal{F}\,$, in order to exactly preserve probability:

[TABLE]

where ${\bm{z}}\equiv\frac{m_{{}_{\scriptsize D}}}{\hbar}({\bm{q}}-{\bm{q}}^{\prime})$ is a dimensionless parameter. The corrected functional then reads

[TABLE]

Finally eq.(3.43), (2.25), (2.26) and (2.27) give the phase and distance decoherence functionals for the case of a single heavy quark in the medium:

[TABLE]

and

[TABLE]

Notice that $\left.\frac{\partial W}{\partial{\bm{x}}}\right|_{{\bm{x}}=0}=0$ gives the required property $\varphi[{\bm{q}},{\bm{p}};{\bm{q}},{\bm{p}}]=0\,$. Moreover, the fact that $W$ is an even function makes the decoherence distance satisfy the symmetric property $D^{2}[{\bm{q}},{\bm{p}};{\bm{q}}^{\prime},{\bm{p}}^{\prime}]=D^{2}[{\bm{q}}^{\prime},{\bm{p}}^{\prime};{\bm{q}},{\bm{p}}]\,$. Equations (3.44) and (3.45) show explicitly how the Lindbladian terms are related to the imaginary part of the potential. Combining the first lines of (3.44) and (3.45), one also obtains the general relation444The fact that here $D$ does not depend on the momenta can be understood, for example, in the case of a Lindblad operator of the form $\hat{L}=\sigma\hat{f}(\hat{\bm{q}})+\mathrm{i}g^{2}\hat{\bm{p}}\>,$ with $\sigma,g^{2}\in\mathbbm{R}\,$, $\hat{f}^{\dagger}=\hat{f}$ and negligible o$(g^{4})$ terms. This Lindblad operator gives $|\Lambda({\bm{q}},{\bm{p}})|^{2}=|\Lambda({\bm{q}})|^{2}+\,\mbox{o}(g^{4})$ and $\mbox{Re}\left[L({\bm{q}},{\bm{p}})L^{*}({\bm{q}}^{\prime},{\bm{p}}^{\prime})\right]=\mbox{Re}\left[L({\bm{q}})L^{*}({\bm{q}}^{\prime})\right]+\,\mbox{o}(g^{4})\,$. This is why eq.(3.2) assumes $\Lambda_{\mu}$ to be momentum independent.

[TABLE]

where the definitions $\varphi\equiv\int_{t_{0}}^{t}\!\!\mathop{}\!\mathrm{d}{\tau}\varphi_{\tau}\,$, $D^{2}\equiv\int_{t_{0}}^{t}\!\!\mathop{}\!\mathrm{d}{\tau}D^{2}_{\tau}$ have been used. Relation (3.2) will be useful to derive the Lindblad equation for the density matrix.

Moreover, from (3.2) and minimising the action, one obtains the classical equations of motion (2.2) for the single quark case. It is interesting to observe that the Lindbladian terms modify the equations of motion from the standard $m\dot{\bm{q}}={\bm{p}}$ and $m\dot{\bm{q}}^{\prime}={\bm{p}}^{\prime}$ to

[TABLE]

This will turn out to be a crucial point when a Fokker-Planck-like equation for the Wigner function of a $q\bar{q}$ pair is derived from a Langevin equation (see Appendix C).

4 Lindblad equation for one heavy quark

This section shows how to get the Lindblad equation for the density matrix $\rho(t,{\bm{q}},{\bm{q}}^{\prime})$ starting from the Lindblad equation for the density operator $\hat{\rho}\,$. The first step is to project eq. (2.11) into position basis:

[TABLE]

The intermediate steps of the calculation are shown in Appendix A. Only the main results are reported here. For a heavy quark in the plasma (without external potential) the first two terms on the right-hand side read

[TABLE]

and

[TABLE]

The last term in (4.49) gives

[TABLE]

Collecting all the results in (4.49), (4.50), (4.51) and (4), one gets

[TABLE]

Using relations (3.44) and (3.45) for $\varphi_{\tau}$ and $D_{\tau}$ respectively, together with integration by parts, the Lindblad equation for the density matrix finally reads

[TABLE]

It is immediate to check that probability is preserved by the Lindblad equation:

[TABLE]

being $\left.\partial_{\bm{q}}W({\bm{q}})\right|_{{\bm{q}}=\mathbb{0}}=0\,$. Notice that (4.54) is the same Lindblad equation derived in Akamatsu:2012vt (see their eq.$(66)\,$) using a different (variational) approach. It is important to observe that the Lindbladian terms in (4.54) have little effect on the diagonal elements of the density matrix. On the other hand, it will be shown that the off-diagonal terms become suppressed at long times. The disappearence of the off-diagonal peaks of the density matrix, which are due to quantum interference, is called quantum decoherence. This phenomenon makes the system undergo a transition from quantum to classical dynamics Zurek2007 .

4.1 Resolution of the quantum probe by the environment

Different types of dynamics of the heavy quark in the bath arise according to the ability of the medium to resolve the subsystem. This ability depends on the correlation length of the environment ($l_{env}$) relevant to the heavy quark-environment system, which corresponds to the typical spatial region in which $W({\bm{x}})$ is significantly different from zero. The environment is expected to be able to resolve, hence affect, the system if and only if its correlation length is of the order of the typical size of the quantum system ($l_{sys}$). This can also be understood from the fact that the typical momentum kicks that the medium gives to the probe are $p_{env}\sim h/l_{env}\,$. This means that a very large correlation length allows the environment to distribute only very small momentum kicks555A possible exception to this argument has to be considered if there are correlated fields in the background instead of random fluctuations., therefore the subsystem is barely perturbed. In the limit of $p_{env}=0\,$, the system continues following its reversible dynamics dictated by the Schrödinger equation.

When $l_{\rm env}\gtrsim l_{\rm sys}$, it will be shown that the system undergoes a Brownian motion, with the bath providing a white noise disturbing the probe. In the $l_{env}\ll l_{sys}$ regime the magnitude of the momentum kicks that the medium gives to the system is much bigger than the typical energies of the latter, and the dynamics of the system is then largely affected by the presence of the environment. In this scenario all bound states are expected to dissociate quickly. Of particular interest is the case with $l_{env}/l_{sys}\lesssim 1\,$, for which the dynamics of the distinct bound states is affected differently by the medium. In this case the fate of the bound states can be regarded as a hadronic calorimeter. Below, this qualitative physical description will be made more quantitative.

4.2 Emergence of different types of dynamics

An inspection of the Lindblad equation will bring to the conclusion that there are four distinct types of dynamics that the system can exhibit:

• •

$l_{env}/l_{sys}\gg 1$

Reversible dynamics;

• •

if $l_{env}/l_{sys}\gtrsim 1$

Open-system dynamics.

The dynamics becomes classical (Langevin-like) after a decoherence time;

• •

$l_{env}/l_{sys}\lesssim 1$

The bath acts like a bound-state sieve;

• •

$l_{env}/l_{sys}\ll 1$

The bath does not discriminate amongst bound states; they all decay at the same rate.

In order to prove these statements let us first perform the following change of variables:

[TABLE]

where the ${\bm{y}}$ coordinate describes the off-diagonal elements of the density matrix. The Lindblad equation (4.54), with an additional external (real) potential $V_{\rm ext}\,$, then reads

[TABLE]

The external potential has been introduced in order to study the dynamics of bound states. Notice again that the diagonal elements $\rho(t,{\bm{r}},\bm{0})$ are almost unaffected by the Lindbladian terms.

It is immediate to see that, if $l_{env}/l_{sys}\gg 1\,$, $W({\bm{x}})$ is constant in the region where the density matrix is different from zero. This means that all Lindbladian terms go to zero in the region where the subsystem is localised. Therefore the subsystem behaves like a closed quantum system obeying a standard Schrödinger equation.

In the scenario where $W$ slightly varies over the region where $\rho$ has support, the former can be Taylor expanded for small $|{\bm{x}}|/l_{env}$. The expansion of $W$ up to second order in its argument is called semiclassical approximation and reads

[TABLE]

where $\mathcal{H}(0)$ is the positive-definite Hessian matrix of $W({\bm{y}})$ evaluated at ${\bm{y}}=0$ and $\left.\partial_{\bm{y}}W({\bm{y}})\right|_{{\bm{y}}=0}=0$ has been used. Assuming that the semiclassical approximation is good enough also for $V_{\rm ext}\,$, the Lindblad equation reads

[TABLE]

Using the fact that $\mathcal{H}(0)$ is a diagonal matrix, and introducing the friction coefficient $\gamma$ and the momentum diffusion coefficient $\kappa$ respectively as (see Blaizot:2015hya ) 666It will be shown below that the momentum diffusion coefficient $\kappa$ gives the strength of the momentum kicks that the system receives from the medium. The associated (position) diffusion coefficient is $D=\frac{2T^{2}}{\kappa}=\frac{T}{m\gamma}$. Observe that $D\sim 1/\left(g^{4}T\ln\left(1/g^{2}\right)\right)$, like all transport scales in perturbation theory Arnold:2000dr . This result comes from the fact that $\gamma\sim\frac{g^{4}T^{2}}{m}\ln(1/g^{2})$, where the expression (3.1) for $W$ has been used together with $\Lambda^{2}/m_{{}_{D}}^{2}\sim 1/g^{2}\gg 1\,$.

[TABLE]

one gets an equation of the Caldeira-Leggett type for the density matrix:

[TABLE]

Notice that the semiclassical approximation preserves the trace of the density matrix. The penultimate term in eq.(4.61) has very little effect on the diagonal peaks of the density matrix, but it causes the off-diagonal ones ($\rho_{\rm off}$) to decay like Zurek2007

[TABLE]

where

[TABLE]

is the decoherence time, $\tau_{{}_{R}}=\gamma^{-1}$ the relaxation time and $\lambda_{\rm th}$ the thermal de Broglie wavelength. For macroscopic objects, the decoherence time is typically much shorter than the relaxation time. Taking the Wigner transform (see 2008AmJPh..76..937C for a good review on Wigner functions and their properties)

[TABLE]

of eq.(4.61), one gets the following Fokker-Planck-like equation for the Wigner function:

[TABLE]

Recall that the Wigner function can be negative and therefore is not a phase space probability. However, after a short time of the order of the decoherence time, the Wigner function becomes positive and can be regarded as a classical probability distribution. In this scenario it is easy to show (see schwabl2006statistical ) that the above Fokker-Planck equation can be generated by a Langevin equation

[TABLE]

where $\bm{\eta}$ is a noise vector satisfying

[TABLE]

The stationary solution of the Fokker-Planck equation is nothing but the Maxwell-Boltzmann equilibrium distribution 777Eq.(4.65) admits a normalised stationary solution if and only if the integral $\int\!\mathop{}\!\mathrm{d}{{\bm{r}}}\,\mathrm{e}^{-\frac{V_{\rm ext}({\bm{r}})}{T}}$ is finite. Notice that this is not the case for an attractive potential, that is $V_{\rm ext}({\bm{r}})<0$ for all ${\bm{r}}\in\mathbb{R}^{3}\,$, since it causes an instability.

[TABLE]

with the normalisation constant

[TABLE]

such that

[TABLE]

which is equivalent to the requirement Tr$\hat{\rho}=1\,$.

The interesting case $l_{env}/l_{sys}\lesssim 1$ will be numerically analysed later for the density matrix of a $q\bar{q}$ pair. The last case corresponds to $l_{env}/l_{sys}\ll 1\,$. In this scenario the imaginary part of the potential has support only in a very small region888Notice that $W(|{\bm{x}}|\to\infty)=0$ because of the definition $W({\bm{x}})=-\int\mathop{}\!\mathrm{d}{t}\,\Delta^{<}(t,{\bm{x}})$ and the fact that the Wightman function $\Delta^{<}$ dies off at large space separation. Moreover $W(0)<0$ since $\Delta^{<}(0,0)=\langle A_{0}(0,0)A_{0}(0,0)\rangle>0$ and $\Delta^{<}(t\to\infty,0)=0\leq\Delta^{<}(t,0)\leq\Delta^{<}(0,0)$ for each time $t\,$.. This means that, excluding the region where ${\bm{y}}={\bm{q}}-{\bm{q}}^{\prime}=0\,$, the Lindblad equation (4.54) with external potential becomes

[TABLE]

The solution of this equation for the off-diagonal (${\bm{q}}\neq{\bm{q}}^{\prime}$) elements of the density matrix is easily found to be

[TABLE]

where $W(0)<0$ and $\rho_{\rm rev}(t,{\bm{q}},{\bm{q}}^{\prime})$ is the solution of the reversible von Neumann equation, i.e. eq.(4.71) without the dissipative term. One immediately notices that the off-diagonal elements of the density matrix become suppressed as time increases, signaling a strong quantum decoherence. In sec.7 it will be shown numerically that the same phenomenon occurs also when $l_{\rm env}\sim l_{\rm sys}\,$, but decoherence manifests itself more slowly. This can be understood in terms of the environment performing frequent measurements on the subsystem, every time the former kicks the latter. When $\rho_{\rm rev}(t,{\bm{q}},{\bm{q}}^{\prime})=\langle{\bm{q}}|\psi(t)\rangle\langle\psi(t)|{\bm{q}}^{\prime}\rangle$ is expanded in the basis of the Hilbert space defined by the Hamiltonian of the closed subsystem ($H|\psi_{n}\rangle=E_{n}|\psi_{n}\rangle$), the solution of eq.(4.71) reads

[TABLE]

where the quantity $\mathrm{e}^{\frac{W(0)}{\hbar}t}|\lambda_{n}|^{2}$ is the probability associated to the $n^{\rm th}$ eigenstate of $H\,$. This means that all the bound states decay with the same rate $-g^{2}W(0)/\hbar=g^{2}T/(4\pi\hbar)\,$, which is also the rate of collisions between one heavy quark (with negligible kinetic energy) and the light quarks of the medium (see Appendix C of Beraudo:2007ky ). Therefore the environment does not discriminate between different bound states (see also Akamatsu:2011se for a similar result). This is not surprising since, in the limit $l_{env}/l_{sys}\ll 1\,$, the large momentum kicks distributed by the medium are much larger than the typical binding energies of the bound states.

Following this argument, one could think that the trace of the density operator is not preserved since probabilities of the eigenstates tend to zero at large time. This reasoning does not hold since the trace is computed on the diagonal elements of the density matrix, which are not described by the solution (4.73).

It is important to bear in mind that $l_{env}$ cannot actually be too small compared to $l_{sys}\,$. In fact, if $l_{env}$ were too small, the momentum kicks of the bath would excite modes of system that are beyond the description of this non relativistic model for the heavy probes.

4.3 Ehrenfest relations and Langevin dynamics

Once the semiclassical approximation has been performed, the stochastic (Langevin) dynamics can be also obtained via the Ehrenfest relations. It is straightforward to show that eq.(4.59) is nothing but the well-known Caldeira-Leggett master equation in the position representation 1983AnPhy.149..374C ; breuer2007theory . The master equation for the density operator reads

[TABLE]

with the friction ($\gamma$) and momentum diffusion ($\kappa$) coefficients defined in (4.60). From this equation one obtains the Ehrenfest equations of motion for the first and second moments of the position and the momentum operators (see Appendix B for the derivation of the first two lines), where the hat on the operators has been omitted to simplify the notation:

[TABLE]

Observe that, due to the non-commuting nature of the operators, this set of equations is not closed for a generic potential, hence the system can not be solved in general. The Ehrenfest relations are closely related to the stochastic differential equations

[TABLE]

with a white noise $\eta(t)\,$, formally given by the time derivative of the Weiner process $W(t)\,$, i.e. $\!\frac{\mathop{}\!\mathrm{d}W(t)}{\mathop{}\!\mathrm{d}t}=\eta(t)\>.$ The white noise satisfies $\langle\eta(t)\rangle_{\eta}=0$ and $\langle\eta(t)\eta(t^{\prime})\rangle_{\eta}=\delta(t-t^{\prime})\,$, implying that the expressions in (4.3) are equivalent to the Langevin form (4.66). It is straightforward to show that the Ehrenfest relations (4.3) have the same form of the noise average of the equations of motion for the first and second moments of the position and the momentum. In particular, the equations for the first moments,

[TABLE]

are nothing but the Langevin equation for the Brownian motion of a particle immersed in a thermal bath.

5 Lindblad equation for a heavy $q\bar{q}$ pair

The same procedure used for the single particle case is applied here to the case of a heavy $q\bar{q}$ pair propagating through the plasma. The interaction part of the influence functional appearing in the path integral (3.32) is derived from eq.(3.1) for $N=2$ and reads

[TABLE]

where the coordinates $({\bm{q}},{\bm{q}}^{\prime})$ and $(\bar{\bm{q}},\bar{\bm{q}}^{\prime})$ refer to the heavy quark and antiquark respectively. The non-interacting term for the single quark is given again by eq.(3.37) and analogously for the antiquark with the replacements ${\bm{q}}\to\bar{\bm{q}}$ and ${\bm{q}}^{\prime}\to\bar{\bm{q}}^{\prime}$. As before, the phase-space path integral for the conditional probability will be derived. For this purpose, one again makes use of the identity (3.38), this time written in the form

[TABLE]

with

[TABLE]

and similarly for the heavy antiquark with the obvious replacements ${\bm{q}}\to\bar{\bm{q}}$ and ${\bm{q}}^{\prime}\to\bar{\bm{q}}^{\prime}$. Neglecting terms of order $g^{4}$ or higher, the propagator in phase space becomes

[TABLE]

with the classical action

[TABLE]

and the in-medium contribution999Notice that the line containing the second derivatives of $W$ is a negligible o$(g^{4})$ term that has been added by hand in order to exactly preserve the trace of the density operator. The same procedure was carried out for a single quark (see eq.(3.42)).

[TABLE]

From the last expression and definition (2.25), one reads off the phase functional and distance decoherence functional for the $q\bar{q}$ pair:

[TABLE]

and

[TABLE]

Observe that $\varphi$ vanishes in the coincidence limit $({\bm{q}}^{\prime},{\bm{p}}^{\prime},\bar{\bm{q}}^{\prime},\bar{\bm{p}}^{\prime})=({\bm{q}},{\bm{p}},\bar{\bm{q}},\bar{\bm{p}})$, and $D^{2}$ is symmetric when $({\bm{q}}^{\prime},{\bm{p}}^{\prime},\bar{\bm{q}}^{\prime},\bar{\bm{p}}^{\prime})$ and $({\bm{q}},{\bm{p}},\bar{\bm{q}},\bar{\bm{p}})$ are swapped, as required by the definitions (2.26) and (2.27). Following a very similar calculation to the one performed in sec.4 for one heavy quark, one finds the Lindblad equation for a $q\bar{q}$ pair in the plasma101010Notice that the Lindblad equation is symmetric under quark-antiquark exchange, as physically expected:

[TABLE]

An interesting feature can be inferred from this result. If the quark and antiquark are taken far apart, which implies that $|{\bm{q}}-\bar{\bm{q}}|,|{\bm{q}}^{\prime}-\bar{\bm{q}}^{\prime}|\gg r_{{}_{D}}\,$, where $r_{{}_{D}}$ is the Debye radius, and consider that111111The distances $|\bar{\bm{q}}-\bar{\bm{q}}^{\prime}|$ and $|\bar{\bm{q}}^{\prime}-\bar{\bm{q}}|$ are of the order of the de Broglie wavelength of the heavy particles.

[TABLE]

one finds that eq.(5) admits a factorised solution

[TABLE]

since both the real and imaginary part of the potential vanish when their arguments are large compared to the Debye radius. The density matrix for the single heavy particle satisfies the Lindblad equation (4.54). Eventually, after an integration by parts, it is immediate to check that the Lindbald equation (5) preserves the trace, that is

[TABLE]

5.1 Tracing out the motion of the centre of mass

Since the motion of the centre of mass of the $q\bar{q}$ pair does not tell anything about the fate of bound states, one would like to trace out this collective motion and find a Lindblad equation for the density matrix describing only the dynamics of the relative $q\bar{q}$ coordinate. To accomplish this, one starts off by defining the centre of mass and relative coordinates respectively as

[TABLE]

and similarly for the primed coordinates. The density matrix describing the dynamics of the relative $q\bar{q}$ coordinate is defined as the trace of the density matrix with respect to the centre of mass coordinate, that is

[TABLE]

In the coincident limit ${\bm{q}}_{\rm cm}^{\prime}\to{\bm{q}}_{\rm cm}\,$, the Lindblad operator does not depend on ${\bm{q}}_{\rm cm}$ anymore, therefore one can carry out the trace over the centre of mass coordinate and get the Lindblad equation for the relative motion:

[TABLE]

Performing a change of variables analogous to (4.56) but with ${\bm{y}}$ halved, that is

[TABLE]

one finally obtains

[TABLE]

Observe that, as in the one particle case, the diffusive term (second line of (5.94)), responsible for quantum decoherence, affects only the off-diagonal elements of the density matrix. Moreover, using integration by parts, it is immediate to check that the trace of the density matrix is preserved.

5.2 Fokker-Planck and Langevin dynamics in the classical limit

In section 4.2 it has been shown that the semiclassical limit of the master equation for a single heavy particle produces a classical stochastic dynamics. The same happens when one takes the semiclassical limit of eq.(5.94) by expanding up to ${\bm{y}}^{2}$ (the expansion is actually in the dimensionless quantity ${\bm{y}}/l_{\rm env}$). The result is

[TABLE]

with the Hessian matrix $\mathcal{H}({\bm{r}})$ defined through $\mathcal{H}({\bm{r}})_{ij}=\partial^{2}W({\bm{r}})/\partial r_{i}\partial r_{j}\,$. Taking the Wigner transform121212Notice the extra factors of $2$ in the definition of the Wigner function. They stem from the fact that a halved ${\bm{y}}$ coordinate has been used in the change of variables (5.93).

[TABLE]

of eq.(5.95), one finds the following Fokker-Planck-like equation for the Wigner function:

[TABLE]

As observed in section 4.2, this becomes a proper Fokker-Planck equation when the Wigner function turns into a classical probability distribution in phase space. This transition happens on time scale of the order of the decoherence time, which is much smaller than the typical relaxation time of the subsystem. Appendix C shows how equation (5.97) can be obtained from the following Langevin equation for the relative motion of the $q\bar{q}$ pair, with position-dependent friction and noise131313Also known as multiplicative noise. (see also Sec. 4.3 of Blaizot:2015hya ):

[TABLE]

where $\mu=m/2$ is the reduced mass of the subsystem and the friction matrix is

[TABLE]

The multiplicative noise $\bm{\eta}({\bm{r}},t)$ is white and autocorrelated, being

[TABLE]

where the momentum diffusion matrix

[TABLE]

satisfies the Einstein relation $\kappa({\bm{r}})=2\mu\gamma({\bm{r}})T$ required by the fluctuation-dissipation theorem applied to the Brownian motion of a particle in a thermal bath schwabl2006statistical .

Notice that a Langevin equation is much easier to be evaluated numerically than a Fokker-Planck equation. However, the information about the quantum $q\bar{q}$ bound states at $t=0$ can not be encoded in the initial conditions of a Langevin equation, which is classical. On the other hand, the quantum initial conditions are easily implemented in the Fokker-Planck-like equation for the Wigner function. One could think of the Langevin equation as the correct formalism to describe the dynamics of the $q\bar{q}$ bound states after a typical time scale of the order of the decoherence time, which is much shorter than the typical relaxation time of the subsystem when $l_{\rm env}\lesssim l_{\rm sys}$, as it will be displayed in sec.7.

6 1D simulation for a heavy $q\bar{q}$ pair

The target of this section is to analyse dissociation, (re)combination and quantum decoherence of bound states of a heavy $q\bar{q}$ pair in a medium. To achieve this, the Lindblad equation (5.94) will be solved numerically using Mathematica.

This part presents a simple model to investigate numerically the effect of the Lindbladian terms on the dynamics and fate of quarkonia. These terms depend on the imaginary part of the potential ($W$), which originates from the Landau mechanism due to the collisions between the subsystem and the medium. The real part of the potential ($V$) is found instead in the reversible part of the Lindblad equation. Surely the screening of the real part of the potential reduces the number of available bound states in the system, however, this is not the effect that this work aims to focus on. The screening mechanism is important when one wants to consider phenomenological implications, but this is not the case of this exploratory analysis. Here the focus is on understanding how the imaginary part of the potential influences the dynamics of the quantum system.

6.1 Bound-state probabilities and linear entropy

The solution of the Lindblad equation is used to compute the following two quantities for studying bound-state dissociation/formation and decoherence respectively:

• •

$P(\psi,t|\psi_{0},t_{0})=\int_{q}\!\int_{q^{\prime}}\psi(q^{\prime})\psi^{*}(q)\rho(t,q,q^{\prime})\>$;

• •

$S_{{}_{L}}=\mbox{Tr}\hat{\rho}-\mbox{Tr}\hat{\rho}^{2}=1-\mbox{Tr}\hat{\rho}^{2}\>$.

The first quantity was introduced previously in eq.(2.4) and represents the probability of finding the subsystem in the state $|\psi\rangle$ at time $t\,$, given that it was origianally in the state $|\psi_{0}\rangle$ at $t=0\,$141414The initial condition has to be specified when solving for the differential equation of the density matrix.. The second object is the so-called linear entropy, which is the leading term of the Mercator series of the von Neumann entropy $S=-\mbox{Tr}[\hat{\rho}\ln\hat{\rho}]\approx\mbox{Tr}[\hat{\rho}({\mathbbm{1}}-\hat{\rho})]=S_{{}_{L}}\,$. The computation of the standard entropy requires the diagonalisation of the density matrix, which can be done only when all the probabilities $\{P_{n}(t)\}_{n}$ of the eigenstates $\{\psi_{n}\}_{n}$ of the Hilbert space of the subsystem are known. This is usually an impossible task, this is why the linear entropy is computed as a proxy of the von Neumann one.

Another reason to consider the linear entropy is that it is a handy quantity for diagnosing quantum decoherence. Indeed a pure state has $S_{{}_{L}}=0\,$, being $|\psi(t)\rangle\langle\psi(t)|=\hat{\rho}(t)=\hat{\rho}^{2}(t)\,$, whereas a mixed state has $0<S_{{}_{L}}\leq 1\,$, since

[TABLE]

If a pure state of the subsystem is plunged into the bath, the linear entropy will increase with a rapidity dependent on the ability of the medium to resolve the quantum system, that is according to the ratio $l_{\rm env}/l_{\rm sys}$ that has been discussed in section 4.2. If this ratio is very small, the environment performs “invasive” measures on the subsystem (the momentum kicks of the bath are large), inducing a sudden raise of the linear entropy. In the opposite regime, the environment can not resolve the system and consequently the linear entropy remains close to zero.

It is interesting to see how the off-diagonal elements of the density matrix affect the time evolution of the linear entropy in the semiclassical limit ($l_{\rm env}\gtrsim l_{\rm sys}$). Using that $\dot{S}_{{}_{L}}=-2\,\mbox{Tr}\left[\dot{\hat{\rho}}\hat{\rho}\right]$ and inserting eq.(5.95), after some algebra one gets

[TABLE]

where it has been used that $\left.\partial_{\bm{r}}^{2}W({\bm{r}})\right|_{{\bm{r}}=\mathbb{0}}-\partial_{\bm{r}}^{2}W({\bm{r}})\approx 0$ when $l_{\rm env}\gtrsim l_{\rm sys}\,$. It is clear that in this regime the linear entropy rises as long as the density matrix has off-diagonal elements. This makes the role of off-diagonal terms in the process of quantum decoherence more transparent.

6.2 Model for simulations

Instead of using the expressions in (3.1) for $V$ and $W\,$, two similar expressions are used here for computational convenience. The real part is taken to be the Pöschl-Teller potential (see Fig.2)

[TABLE]

where $\omega$ is a unit of energy, $\mu=m/2$ is the reduced mass of the $q\bar{q}$ pair and $j\in\mathbb{N}\,$. This potential has been chosen since it admits $j$ bound states,

[TABLE]

with $n=0,1,\dots,j-1,j\,$, $C_{n}^{j}$ a dimensionless normalisation factor and $\mbox{P}_{j}^{n}(x)$ the associated Legendre polynomial. Therefore the maximum number of vacuum ($T=0$) bound states can be simply tuned by changing the natural number $j\,$, which is set to $2$ here. The imaginary part is taken to be (see Fig.2)

[TABLE]

where the correlation length of the environment is $l_{\rm env}=\hbar/(\sigma\,T)$ and $\sigma$ is a dimensionless parameter that will be tuned to probe the different physical regimes discussed in section 4.2. This form of the correlation length resembles the one of the Debye radius $r_{{}_{D}}=m_{{}_{D}}^{-1}=\hbar/(\alpha g\,T)$ appearing in the original expression (3.1) of $W\,$, which is plotted in Fig.2 as $\tilde{W}(x)=-\frac{Tl_{\rm env}}{|x|}\int_{0}^{5}\mathop{}\!\mathrm{d}{z}\frac{\sin\left[\frac{2z|x|}{l_{\rm env}}\right]}{(z^{2}+1)^{2}}\,$. Comparing with the expression derived from the gauge theory, a factor $\frac{1}{4\pi}$ has been omitted in front of $\tilde{W}(x)\,$. This has been done since the model is one dimensional anyway, and the form (3.1) is strictly valid only in three dimensions.

In the numerical simulations, the bound-state probabilities ($P_{0}(t)$ and $P_{1}(t)$ for the ground and excited state respectively) and the linear entropy will be extracted from the solution of the following one-dimensional Lindblad equation for the relative motion of a $q\bar{q}$ pair:

[TABLE]

This equation is the one-dimensional analog of eq.(5.94) with $g^{2}=1$ 151515Notice that the perturbative approximation is still valid since the expansion is actually in powers of $\alpha_{s}=\frac{g^{2}}{4\pi}\,$. Moreover $T=0.8$ GeV and $m=1.2$ GeV have been used here, so that the nonrelativistic limit is valid to a good degree of accuracy. and $V(x),\,W(x)$ defined as in (6.104),(6.106) respectively. The first part of the analysis assumes that initially the $q\bar{q}$ pair is in one of the two bound states of the Pöschl-Teller potential, that is

[TABLE]

where $n=0,1$ for the ground state and excited state respectively. If the bound states melt, the density matrix describes scattering states, whose density matrix does not go to zero at infinity. To avoid this numerical issue, the system is inserted in a harmonic box via the following modification of the real part of the potential161616The bound-state probabilities and linear entropy computed in the following section do not depend on the chosen value of $\epsilon$, as long as it is small enough to allow the ground and excited states of $V(x)$ to be eigenstates also of the modified potential $V(x)+\epsilon x^{2}\,$.:

[TABLE]

with $\epsilon$ very small, so that the Pöschl-Teller potential gets significantly modified only for $|x|\gtrsim 2$ fm, where the wavefunctions of the bound states are basically zero. Therefore, for numerical purposes, the eigenstates of the Pöschl-Teller potential are also eigenstates of the Pöschl-Teller potential plus the harmonic contribution. In addition, the new potential admits an infinite set of more energetic bound states, which can be interpreted as localized “scattering” states of the unmodified Pöschl-Teller potential. This trick allows one to use vanishing boundary conditions for the density matrix when solving (6.107) in a finite box.

In the cases analysed below, it has been checked numerically that the density operator is positive-semidefinite, that is $\rho(t,q,q)\geq 0$ for all times $t$ and positions $q$ (equivalently $\rho(t,r,y=0)\geq 0$ $\forall t,r\,\in\mathbb{R}\,$), as required by a master equation in the Lindblad form. Moreover, eq.(6.107) preserves the trace of the density operator to the value of one with an error of the order of one part in a thousand.

7 Numerical results

7.1 Dissociation of bound states

The probabilities of having one of the two bound states at time $t\,$, given an initial density matrix corresponding to the ground (excited) state, are shown in the top left (right) panel of Fig.3. The time evolution of the linear entropy for the two different initial density matrices is shown just in the corresponding panels below. The analysis of the outcomes is divided in three cases, according to the correlation length of the environment.

• •

$l_{\rm env}>l_{\psi_{1}}$ regime

In this scenario the correlation length of the bath is bigger than the size of the excited state, which in turns is larger than the size of the ground state. This means that the environment can “measure” the excited state but can very hardly resolve the ground state. This has indeed been found on the left panels of Fig.3, where the initial ground state (see circles) remains basically unaffected as time passes ($P_{0}(t)\gtrsim 90\%$ $\forall t\,$). On the other hand, there is a small probability of having a transition from the ground state to the excited state($P_{1}(t)\lesssim 10\%$ $\forall t\,$), but there is basically no change of dissociation since $P_{0}(t)+P_{1}(t)\approx 1=P_{0}(t)+P_{1}(t)+P_{\rm scatt}(t)$ at all times, where $P_{\rm scatt}$ is the probability of having some “scattering” states. Not only the medium cannot give large enough momentum kicks to disturb the ground state, but also cannot resolve it, hence quantum decoherence should manifest itself as a very mild effect. This is indeed seen in the small values taken by the linear entropy, indicating that the initial state remains approximately pure.

When starting off with the excited state (circles on the right panels of Fig.3), the linear entropy increases quite rapidly up to a stationary value of $S_{{}_{L}}\approx 0.8\,$, meaning that the measurements that medium makes on the subsystem bring relatively quickly the excited state from a pure to a mixed state. This is the signal of the transition from a quantum to a classical system. The kicks distributed by the bath can (with roughly the same probability of $40\%$ in equilibrium, i.e. for $t\gtrsim 5$ fm$/c$) either bring about a feed-down mechanism of the excited state to the ground state, or melt the bound state. The final probability of remaing in the excited state is only about $20\%\,$.

• •

$l_{\psi_{0}}<l_{\rm env}<l_{\psi_{1}}$ regime

In this scenario the environment can resolve very well the excited state and also start affecting the ground state. From the diamond-shaped symbols on the left panels of Fig.3 one can see that an initial ground state has survival probability $P_{0}(t)\geq 80\%\,$, which is high but less than in the previous case. The chances of promoting the ground state to the excited state are slightly higher than before but still fewer than $10\%\,$. Similarly, the linear entropy grows more than before but stays below the value of $0.4\,$. On the other hand, if the system is initially in the excited state (diamonds on the right panels of Fig.3), the environment strongly perturbs the bound state, causing the linear entropy to increase abruptly and the excited state to disappear very rapidly, reaching soon the equilibrium value of about $10\%$ for the survival probability $P_{1}$. One might think that the disappearance of the excited state is simply a signal of its melting due to the kicks received from the bath. However, the feed-down mechanism is more important here than in the previous case, in which the excited state receives lighter kicks. In fact the probability of ending up with the ground state grows rapidly from zero to the equilibrium value of $60\%\,$, which leaves a melting probability of $30\%\,$.

• •

$l_{\rm env}<l_{\psi_{0}}$ regime

Here the environment can fully resolve both bound states, causing a sudden increase of the linear entropy for any initial density matrix (see triangles in lower panels of Fig.3). If the system is in the ground state at $t=0\,$, its survival probability will reduce quite quickly to the stationary value of about $55\%\,$. The ground state may then form an excited state with $P_{1}\approx 15\%$, or melt with a probability of approximately $30\%\,$. The case corresponding to the excited state at $t=0$ is basically identical to the previous scenario. The only difference is that the survival probability of the excited state is slightly higher than before, and the linear entropy reaches a lower equilibrium value of $0.5\,$. Observe that, in this case, the equilibrium values for the probabilities and linear entropy for the two different initial density matrices are very similar. This is expected since in the limit $l_{\rm env}\to 0$ the momentum kicks coming from the bath have such a much larger value than the typical binding energies of the subsystem that, at long time scales, it does not matter which particular initial state the subsystem was prepared in.

In this paragraph only the initial conditions corrisponding to the ground and excited states have been considered. However, one can start off with an initial density matrix corresponding to any linear combination of ground and excited state. This is an important feature if one wants to consider a realistic experimental initial condition in which both the ground and excited states of a system are occupied.

7.2 Role of scattering states

It is clear that the presence of the scattering states in the Hilbert space of the subsystem is crucial for preserving probability, that is $\mbox{Tr}\hat{\rho}(t)=1\,$. However, one might think that the presence of the scattering states is not that important in certain cases. After all, the previous paragraph has shown that, for certain correlation lengths of the medium, the probability of bound states dissociation or, equivalently, the formation of scattering states is small. For example, the case $l_{\psi_{0}}<l_{\rm env}<l_{\psi_{1}}$ shows that $P_{\rm scatt}(t)=1-P_{0}(t)-P_{1}(t)\lesssim 0.1\,(0.3)\>\forall t$ when starting off with the ground (excited) state. Consequently one might hope that for this scenario the bound states probabilities $P_{0}(t)$ and $P_{1}(t)$ could be approximately computed without solving the full Lindblad equation for the density matrix, but expanding the density matrix in terms of the bound states only

[TABLE]

and solving for the coefficient $\lambda_{nm}(t)$ by taking the time derivative

[TABLE]

and using the Lindblad equation (6.107) with the replacement $q=r+y$ and $q^{\prime}=r-y\,$. One then solves the linear first-order differential equations for the probabilities $\lambda_{00}(t)=P_{0}(t)$ and $\lambda_{11}(t)=P_{1}(t)\,$, which are decoupled from the equations for the mixed coefficients $\lambda_{01}(t)$ and $\lambda_{10}(t)\,$. Observe that the ansatz (7.2) does not take scattering states into account, hence does not preserve the trace of the density matrix. This technique has been used in the literature (e.g. Akamatsu:2011se ) for computing approximate bound states survival probabilities.

Fig.4 shows that, if ansatz (7.2) is assumed, the time evolution of the bound states probabilities are quite different from the correct ones computed from the full Lindblad equation, which implicitly knows about the whole Hilbert space of the subsystem, which comprises the scattering states. Notice that the evolution of the probability of the excited state $P_{1}(t)$ with the ansatz (7.2) is similar to the correct one computed from the Lindblad equation. The situation is drastically different for the ground state, whose probability takes values significantly lower than the exact ones. This happens because (7.2) implies that all the states that “leak out” from the potential well have no chance of going back inside it, since the ansatz does not take into account scattering states.

The next subsection is going to show more explicitly that scattering states indeed roll back into the potential well to form the ground state with a non negligible probability. The feed-down mechanism from a scattering state to an excited state is instead a marginal effect, shedding light on why the excited-state probabilities calculated either via solving the full Lindblad equation or using the anstaz (7.2) are quite similar.

7.3 Recombination

This paragraph aims to show how (re)combination171717The term recombination is used when the scattering state originally came from a melted $q\bar{q}$ bound state. of a scattering state into a bound state takes place. In order to implement vanishing spatial boundary conditions for the density matrix when solving the Lindblad equation (6.107), one has to consider normalised scattering states (with positive energies) defined as

[TABLE]

where $\delta=\sqrt{2}\sqrt{\langle\hat{x}^{2}\rangle}$ is the localisation parameter and $p=\langle\hat{p}\rangle$ is the expectation value of the momentum. In principle, the localised scattering state can be written as a superposition of elements of the Hilbert state of the $T=0$ Hamiltonian with the unbounded potential (6.109), but this can not be done in practice since the full Hilbert space of the subsystem is unknown. An interesting case to analyse is the evolution of a thermal scattering state, that is a scattering state with size equal to the thermal de Broglie wavelength, i.e. $\sqrt{\langle\hat{x}^{2}\rangle}=\lambda_{\rm th}=h/p_{\rm th}\,$, with thermal momentum $p_{\rm th}=\sqrt{\mu\,T}\,$, $\mu=m/2\,$. The temperature $T=0.8$ GeV has been used again in the numerical analysis.

Fig.5 shows the probability of the thermal scattering state to become a bound state (either the ground or excited state) for different $l_{\rm env}\,$. When $l_{\rm env}>\sqrt{\langle x^{2}\rangle}_{\psi_{{}_{\rm scatt}}}\,$, both $P_{0}$ and $P_{1}$ do not vary with time (see square-shaped symbols). This entails that the recombination probability, defined as $P_{\rm rec}(t)=\sum_{i=0,1}\left(P_{i}(t)-P_{i}(0)\right)$, remains practically zero. A very similar behaviour is found for $l_{\rm env}\sim\sqrt{\langle x^{2}\rangle}_{\psi_{{}_{\rm scatt}}}$ (see circle-shaped symbols), with the only difference that here the probability of getting the ground state slightly increases with time, whereas the probability of getting the excited state slightly diminishes.

For these two correlation lengths of the environment the linear entropy grows relatively rapidly (see Fig.6). The scattering state is indeed a superposition of eigenstates of the Hamiltonian of the subsystem with the modified potential (6.109), and the medium disturbs all the eigenstates with wavelengths either bigger or of the same order of $l_{\rm env}$. This effect becomes obviously stronger when $l_{\rm env}$ becomes smaller than the size of the scattering state, causing the linear entropy to raise abruptly up to values very close to the maximum value of one. This means that the subsystem has been strongly perturbed by the medium and the quantum information is quickly lost. The fact that the dynamics becomes basically classical does not prevent the subsystem to form bound states, as can be seen from the finite recombination probability inferred from Fig.5 (see diamonds and triangles). In fact the probability of the scattering state to form the ground state, in the case with $l_{\rm env}=0.12$ fm, grows from the $t=0$ overlap value of approximately $9\%$ to $42\%$ at $t=5$ fm/c. On the other hand, the overlap with the excited states decreases of about $6\%\,$. Observe that the long-time value of the linear entropy for the two smallest values of $l_{\rm env}$ is determined very accurately by $P_{0}(t)$ only. For example, for $l_{\rm env}=0.12$ fm one has $0.8=S_{{}_{L}}\approx 1-P_{0}^{2}=0.82$ at $t=5$ fm/c.

Eventually Fig.7 showcases some direct results obtained by solving the Lindblad equation (6.107) with $l_{\rm env}=0.25$ fm, starting off with the thermal scattering state defined above. It is clear that quantum decoherence manifests itself at the very early stages via the suppression of the off-diagonal elements ($y\neq 0$) of the real part of the density matrix. After barely $\Delta t=0.2$ fm/c the density matrix considerably changes (explaining the sudden raise of the linear entropy in Fig.6), and at $t=1$ fm/c approaches a diagonal form. The imaginary part of $\rho$ has not been displayed since it goes close to zero very quickly, leaving the density matrix basically real.

By considering just the time evolution of the diagonal elements of the density matrix, one can check that positivity is indeed preserved, that is $\rho(t,x,0)\geq 0$ for all positions $x$ and times $t$ ($\,\rho(t,x,0)\in\mathbb{R}$ $\forall x,t\,$), as can be seen in Fig.8181818Note that positivity of the density operator comes automatically from the Lindblad form (2.11). However, equation (6.107) is an almost exact Lindblad equation, due to the very small terms neglected in perturbation theory. Therefore it can happen that in certain regimes, in which the approximations used throughout this work are not valid anymore, the density operator sligthly violates positivity..

In order to see how the dissociation and recombination mechanisms depend on the mass of the heavy quarks, one can find the numerical results for the bound-state probabilities for a heavier quark with mass ($m=4.7$ GeV) in Appendix D.

8 Summary and Outlook

It has been shown how to derive a Lindblad equation for non relativistic heavy quarks and antiquarks propagating out of equilibrium in a thermalised quark gluon plasma, within the framework of open quantum systems and starting from the underlying gauge theory. To achieve this, an abelian model for the plasma has been used, together with some well-defined approximations, including the perturbative expansion and the Markovian limit. All the Lindbladian terms in the master equation depend on the imaginary part of the inter-quark potential, whose shape gives the correlation length of the environment. Different types of dynamics of the heavy probes emerge according to the value of the ratio between the correlation length of the environment and the typical size of the probes. The smaller the correlation length of the medium, the higher the resolution of the quantum probes, giving rise to quantum decoherence and to classical dynamics. Numerical simulations of the one-dimensional Lindblad equation for a heavy quark-antiquark pair allow one to study quantitatively the mechanisms of melting and formation of bound states within a unique framework. It has been found that the interaction of bound states with the medium can make the bound state dissociate or undergo a feed-down mechanism to a bound state with lower energy. Similarly, a scattering state (that can be interpreted as a melted bound state) has a non-zero probability of (re)combining into a bound state.

In principle, dissociation and formation of bound states can also be analysed starting from the same initial density matrix by solving the Lindblad equation for more than one heavy $q\bar{q}$ pair, which can be straightforwardly derived with the same technique exploited here for the two-particle case. The issue with this procedure is that a numerical simulation of such equation becomes very challenging for more than two particles, already in one dimension. A solution to this problem would be to solve the stochastic Schrödinger equation PhysRevD.13.857 ; 0305-4470-25-21-023 ; Gisin1989 associated to the Lindblad equation. The advantage of this approach is that the dimension of a wave vector is the square root of that of a density matrix. This strategy has been exploited by Akamatsu:2011se in one dimension in the framework of quarkonium in a quark gluon plasma, but only the recoilless limit (no friction) has been taken into account.

Another interesting and numerically feasible approach to study the real-time dynamics of many particles in the plasma is the Langevin equation, which has been proven to be the correct description when the system has undergone a quantum to classical transition. However, in order to be able to use the classical stochastic equation, one needs to understand how to implement the initial quantum state correctly, since the Langevin formalism is strictly valid only after a time scale of the order of the decoherence time. It has been found that the Langevin equation derived from the gauge theory for a $q\bar{q}$ pair (see also Blaizot:2015hya for more particles) features a space-dependent friction coefficient, which encodes physical information about polarisation effects of the medium induced by the relative position of the heavy probes. Apart from Blaizot:2015hya , this feature has not yet been considered in the literature when solving the Langevin equation for quarkonia. It would be interesting to see whether this space-dependent friction for a heavy $q\bar{q}$ pair could be reproduced for strongly coupled plasmas in the context of the AdS/CFT correspondence, as it has already been done for a single heavy quark PhysRevD.74.085012 ; Son:2009vu 191919This idea has been suggested by Prem Kumar. Observe that single-quark diffusion and heavy $q\bar{q}$ pairs are treated differently in holography..

Eventually, in order to envisage some phenomenological applications, a natural development of this project would be to obtain a Lindblad equation in QCD using the same formalism presented here (see Akamatsu:2012vt ; Akamatsu:2014qsa ; Brambilla:2016wgg for other approaches). One would expect to be able to implement non abelian features within the present scheme, this time considering the weak coupling approximation more carefully, being the non-linear QCD couplings less strongly suppressed than the QED ones. This route has been followed in the preliminary work prep , in which the QED case can be obtained as a restriction of the QCD case.

Acknowledgments

I am very grateful to Jean-Paul Blaizot and Pietro Faccioli for fruitful discussions on theoretical aspects of this paper and for having invited me at IPhT (Saclay) and Trento University respectively. A special acknowledgment goes to Giovanni Garberoglio, who carefully checked all the numerical results, and to Gert Aarts, who provided interesting comments and did a meticulous job in proofreading the manuscript. I am also thankful to Prem Kumar, Jamie Mc Donald, Guy Jehu and Markus Muller for valuable remarks and suggestions, and to my supervisor Chris Allton, for giving me the freedom to pursue this project.

Appendix A Derivation of the Lindblad equation for one heavy quark

In order to obtain the Lindblad equation (4.54) for one heavy quark one has to compute the matrix elements in (4.49). The first one on the right hand side reads ($\int_{\bm{p}}\equiv\int\frac{\mathop{}\!\mathrm{d}{{\bm{p}}}}{(2\pi\hbar)^{3}}$)

[TABLE]

where it has been used that $|\Lambda_{\mu}|^{2}$ is momentum independent up to order $g^{2}$. The same calculation applies to the second term on the right hand side of (4.49) and gives expression (4.51) as a result. The computation of the last term in (4.49) reads

[TABLE]

Collecting all the results as shown in section 4, one finally gets to the Lindblad equation (4.54).

Appendix B Derivation of the Ehrenfest equations of motion

Here the proof of the first two equations of motion in (4.3) is shown. The derivation of the other relations is very similar and is therefore omitted here. The Ehrenfest equation of motion for the first momentum of the position operator is obtained through the following steps (hats on operators have been removed for ease of notation):

[TABLE]

where in the last line the cyclicity of the trace has been exploited to cancel the friction and diffusion terms. The fact that the position operator commutes with the potential operator has been also used. The commutator $[q,p]=\mathrm{i}\hbar$ then gives

[TABLE]

The second relation of (4.3) follows from a similar calculation. One has

[TABLE]

since

[TABLE]

The first contribution in (B) becomes

[TABLE]

where the relation $\left[p,V_{\rm ext}(q)\right]=-\mathrm{i}\hbar V^{\prime}_{\rm ext}(q)$ has been used. The second term of (B) reads

[TABLE]

Expressions (B), (B.119) and (B) together give the second Ehrenfest relation (4.3).

Appendix C From the Langevin equation with multiplicative noise to the Fokker-Planck equation

A stochastic differential (Langevin) equation for a multi-component process can be written in the general form

[TABLE]

with the white noises satisfying

[TABLE]

In 2007cond.mat..1242G (see Sec. 5.3.1) it is shown that the Fokker-Planck equation associated to the Markov process described by (C.121) is

[TABLE]

From now on $\lambda$ is set to $1/2\,$. The goal of this section is to prove that the Langevin equation (5.98) with multiplicative noise generates the Fokker-Planck-like equation (5.97). The Langevin equation (5.98) can be cast (up to o($g^{2}$) terms) in the multivariate form (C.121) as

[TABLE]

with the usual white noise $\eta(t)$ and the function $G({\bm{r}})$ defined via

[TABLE]

It is interesting to see that the first line of eq.(C) does not correspond to the classical $\dot{\bm{r}}={\bm{v}}\,$, but it contains an extra perturbative quantum correction that can be traced back to the modified classical equations of motion (2.2) derived from minimisation of the action containing the Lindbladian terms. In order to avoid a cumbersome notation, only the one-dimensional case is considered here. The proof in $d>1$ dimensions works in the same way. Contrasting eq.(C) with eq.(C.121) and using ${\mathbf{S}}=(S_{1},S_{2})\equiv(r,v)\,$, one obtains the following coefficients for the multivariate Langevin equation:

[TABLE]

Substituting these terms in (C), one gets

[TABLE]

whence, using (C.125) and $v=p/\mu\,$,

[TABLE]

which has exactly the same form of the Fokker-Planck-like equation (5.97) for the Wigner function with $m=2\mu\,$.

Appendix D Mass dependence of dissociation and recombination

In this paragraph the dissociation and recombination mechanisms for two different values of heavy-quark mass are compared. This analysis is inspired by the experimental evidence that the two phenomenons are different for charmonium and bottomonium. In particular recombination effects are very small for bottomonium. Here the values $m_{c}=1.2$ GeV and $m_{b}=4.7$ GeV have been used for the mass of the charm and bottom quark respectively.

Contrasting Fig.9 with Fig.3 of section 7, it is clear that the ground state of the heavier quarkonium is less affected by the medium than the lighter quarkonium (see bound-state probabilities on the left panel of Fig.9, where the initial density matrix corresponds to the ground state).

This can be understood from the fact that for the heavier quarkonium all the correlation lengths of the environment ($l_{\rm env}$) are greater than or equal to the size of the ground state ($l_{\psi_{0}}$), hence the system is mildly perturbed by the medium. The same behaviour is found when one starts off with the excited state and $l_{\rm env}>l_{\psi_{1}}$ (see circles and diamonds on the right panel of Fig.9). However, when $l_{\rm env}<l_{\psi_{1}}$ (see triangles on the right panel of Fig.9), the bound-state probabilities evolve in time very similarly to the ones in Fig.3 for the lighter quarkonium. In fact, in this case the dissociation process seems slightly more significant for the heavier quarkonium than for the lighter one.

Fig.10 shows the time evolution of the bound-state probabilities when the initial state corresponds to a thermal scattering state ($T=0.8$ GeV), as defined in eq.(7.112), for a fixed $l_{\rm env}$. If on one hand the probability of getting the excited state at time $t$ (right panel) does not depend much on the mass of the heavy quark, on the other hand the probability of the scattering state of getting trapped into the ground state is largely dependent on the mass of quarkonium. In fact the left panel of Fig.10 shows that the scattering state with the heavier mass does not tend to form the ground state, whereas the one with the lighter mass clearly does.