Quantum Brownian motion of a heavy quark pair in the quark-gluon plasma

TL;DR

This paper develops a quantum master equation approach to study how dissipation influences the real-time evolution and thermalization of heavy quarkonium in the quark-gluon plasma, linking QCD to phenomenological models.

Contribution

It introduces a novel non-equilibrium master equation derived from quantum field theory for quarkonium dynamics in QGP, connecting QCD with non-linear Schrödinger models.

Findings

Dissipative effects enable quarkonium to thermalize within QGP lifetime.

The quantum evolution can be represented as a stochastic non-linear Schrödinger equation.

Dissipation impacts quarkonium dynamics at early stages in heavy ion collisions.

Abstract

In this paper we study the real-time evolution of heavy quarkonium in the quark-gluon plasma (QGP) on the basis of the open quantum systems approach. In particular, we shed light on how quantum dissipation affects the dynamics of the relative motion of the quarkonium state over time. To this end we present a novel non-equilibrium master equation for the relative motion of quarkonium in a medium, starting from Lindblad operators derived systematically from quantum field theory. In order to implement the corresponding dynamics, we deploy the well established quantum state diffusion method. In turn we reveal how the full quantum evolution can be cast in the form of a stochastic non-linear Schr\"odinger equation. This for the first time provides a direct link from quantum chromodynamics (QCD) to phenomenological models based on non-linear Schr\"odinger equations. Proof of principle simulations in one-dimension show that dissipative effects indeed allow the relative motion of the constituent quarks in a quarkonium at rest to thermalize. Dissipation turns out to be relevant already at early times well within the QGP lifetime in relativistic heavy ion collisions.

Tables

Table 1. Table 1: Numerical setup, parametrization of V ​ ( x ) 𝑉 𝑥 V(x) and D ​ ( x ) 𝐷 𝑥 D(x) , and center-of-mass momentum

$\Delta x$	$\Delta t$	$N_x$	$\gamma$	${\ell }_{\mathrm{corr}}$	$\alpha$	$m_D$	$x_c$	$P_{\mathrm{CM}}$
$1/M$	$0.1M{\left(\Delta x\right)}^2$	254	$T/\pi$	$1/T$	$0.3$	$2T$	$1/M$	0

Equations

$$\displaystyle\partial_{t}\rho_{\rm{S}}(t)=-i[H_{\rm{S}},\rho_{\rm{S}}(t)]+\sum_{i}\Bigr{(}2L_{i}\rho_{\rm{S}}(t)L_{i}^{\dagger}-L_{i}^{\dagger}L_{i}\rho_{\rm{S}}(t)-\rho_{\rm{S}}(t)L_{i}^{\dagger}L_{i}\Bigr{)},$$

$${\rho }_{\mathrm{S}}\equiv {\mathrm{Tr}}_{\mathrm{E}}{\rho }_{\text{tot}},$$

$${\partial }_t{\rho }_{\mathrm{S}}=\mathcal{L}{\rho }_{\mathrm{S}},$$

$$\displaystyle\partial_{t}\rho_{\rm{S}}(t)=-i[H_{\rm eff},\rho_{\rm{S}}(t)]+\sum_{i}\Bigr{(}2L_{i}\rho_{\rm{S}}(t)L_{i}^{\dagger}-L_{i}^{\dagger}L_{i}\rho_{\rm{S}}(t)-\rho_{\rm{S}}(t)L_{i}^{\dagger}L_{i}\Bigr{)}.$$

$$H_{\mathrm{eff}}$$

$$L_{\vec{k}a}$$

$$\frac{2g^2}{N_c^2-1}{\int }_0^{\infty }dt⟨A_0^a(t,\vec{r})A_0^a(0,\vec{0})⟩=D(\vec{r})-iV(\vec{r}),$$

$$D(\vec{r})=g^{2}T\int \frac{d^{3}k}{(2\pi {)}^3}\frac{\pi m_{\mathrm{D}}^2{e}^{i\vec{k}\cdot \vec{r}}}{k(k^2+m_{\mathrm{D}}^2{)}^2},V(\vec{r})=-\frac{g^2}{4\pi r}{e}^{-m_{\mathrm{D}}|\vec{r}|},$$

$$\mathrm{Re}[V_{Q\bar{Q}}](\vec{r})=V(\vec{r}),\mathrm{Im}[V_{Q\bar{Q}}](\vec{r})=D(\vec{r})-D(0).$$

$$\vec{R}=\frac{{\vec{x}}_Q+{\vec{x}}_{\bar{Q}}}{2},\vec{r}={\vec{x}}_Q-{\vec{x}}_{\bar{Q}},\vec{P}={\vec{p}}_Q+{\vec{p}}_{\bar{Q}},\vec{p}=\frac{{\vec{p}}_Q-{\vec{p}}_{\bar{Q}}}{2},$$

$$H_{\mathrm{eff}}$$

$$L_{\vec{k}a}$$

$$-\sqrt{\frac{\tilde{D}(\vec{k})}{2L^3}}{e}^{i\vec{k}\cdot \vec{R}}\left[1-\frac{\vec{k}}{4MT}\cdot \left(\frac{1}{2}\vec{P}-\vec{p}\right)\right]e^{-\frac{i\vec{k}\cdot \vec{r}}{2}}(1\otimes t^{a\ast }).$$

$${\mathcal{H}}_{\mathrm{Q}\bar{\mathrm{Q}}}={\mathcal{H}}_R\otimes {\mathcal{H}}_r,$$

$${\rho }_r\equiv {\mathrm{Tr}}_R({\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}).$$

$${\mathrm{Tr}}_R\left[H_{\mathrm{eff}},{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}\right],{\mathrm{Tr}}_R\left(L_{\vec{k}a}{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}{L}_{\vec{k}a}^{†}\right),{\mathrm{Tr}}_R\left(L_{\vec{k}a}^{†}{L}_{\vec{k}a}{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}\right),{\mathrm{Tr}}_R\left({\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}{L}_{\vec{k}a}^{†}{L}_{\vec{k}a}\right).$$

$${\mathrm{Tr}}_R[H_{\mathrm{eff}},{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}]=\left[H_{\mathrm{eff}}^{(r)},{\rho }_r\right].$$

$$L_{\vec{k}a}=e^{i\vec{k}\cdot \vec{R}}{C}_{\vec{k}a}(O_r,\vec{P}),L_{\vec{k}a}^{†}=C_{\vec{k}a}^{†}(O_r,\vec{P})e^{-i\vec{k}\cdot \vec{R}},O_r=\left\{\vec{r},\vec{p},t^a\otimes 1,1\otimes t^{a\ast }\right\},$$

$${\mathrm{Tr}}_R\left(L_{\vec{k}a}{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}{L}_{\vec{k}a}^{†}\right)$$

$${\mathrm{Tr}}_R\left(L_{\vec{k}a}^{†}{L}_{\vec{k}a}{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}\right)$$

$${\mathrm{Tr}}_R\left({\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}{L}_{\vec{k}a}^{†}{L}_{\vec{k}a}\right)$$

$${\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}=|{\vec{P}}_{\mathrm{CM}}⟩⟨{\vec{P}}_{\mathrm{CM}}|\otimes {\rho }_r.$$

$${\mathrm{Tr}}_R\left(L_{\vec{k}a}{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}{L}_{\vec{k}a}^{†}\right)$$

$${\mathrm{Tr}}_R\left(L_{\vec{k}a}^{†}{L}_{\vec{k}a}{\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}\right)$$

$${\mathrm{Tr}}_R\left({\rho }_{\mathrm{Q}\bar{\mathrm{Q}}}{L}_{\vec{k}a}^{†}{L}_{\vec{k}a}\right)$$

$$\frac{d}{dt}{\rho }_r(t)$$

$$H_{\mathrm{eff}}^{(r)}$$

$$L_{\vec{k}a}^{(r)}$$

$$-\sqrt{\frac{\tilde{D}(\vec{k})}{2L^3}}\left[1-\frac{\vec{k}}{4MT}\cdot \left(\frac{1}{2}{\vec{P}}_{\mathrm{CM}}-\vec{p}\right)\right]e^{-\frac{i\vec{k}\cdot \vec{r}}{2}}(1\otimes t^{a\ast }).$$

$${\rho }_{\mathrm{S}}=\mathrm{M}[|\psi (t)⟩⟨\psi (t)|],$$

$$|d\psi ⟩$$

$$\displaystyle=-iH_{\rm{eff}}\ket{\psi(t)}dt+\sum_{i}\Bigr{(}2\langle L^{\dagger}_{i}\rangle_{\psi}L_{i}-L^{\dagger}_{i}L_{i}-\langle L^{\dagger}_{i}\rangle_{\psi}\langle L_{i}\rangle_{\psi}\Bigr{)}\ket{\psi(t)}dt$$

$$\displaystyle\hskip 109.54327pt+\sum_{i}\Bigr{(}L_{i}-\langle L_{i}\rangle_{\psi}\Bigr{)}\ket{\psi(t)}d\xi_{i}\,,$$

$$\mathrm{M}[d{\xi }_i]=\mathrm{M}[\Re (d{\xi }_i)\Im (d{\xi }_j)]=0,$$

$$\mathrm{M}[\Re (d{\xi }_i)\Re (d{\xi }_j)]=\mathrm{M}[\Im (d{\xi }_i)\Im (d{\xi }_j)]={\delta }_{ij}dt.$$

$$|d\psi (t)⟩$$

$${\rho }_{\mathrm{S}}(t)$$

Full text

††preprint: APS/123-QED

I Introduction

Over the past decades, properties of nuclear matter in extreme conditions have been vigorously studied both experimentally and theoretically.

At modern collider facilities, such as the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC), heavy nuclei are collided at ultrarelativistic energy to create a multiparticle system in the collision center, endowed with an extremely high energy density. A multitude of measurements suggest that the energy densities reached are high enough that a new phase of nuclear matter, the “quark-gluon plasma (QGP)” is realized [1, 2, 3].

On the theory side significant progress has been made in understanding the thermodynamic, i.e. static properties of the QGP at the temperatures reached in heavy ion collisions at which QCD dynamics is still nonperturbative. Lattice QCD simulations have played a central role in this regard shedding light on e.g. the crossover transition temperature [4, 5], the physics of static screening in QCD [6, 7] and the equilibrium spectral properties of heavy quarkonium [8, 9, 10, 11, 12, 13, 14]. On the other hand our understanding of the dynamical properties of QCD is much less developed due to the notorious sign problem preventing lattice QCD simulations from being performed directly in Minkowski time. Perturbation theory may provide some input on dynamical properties but is applicable only at much higher temperatures than achievable in experiments.

One way to handle dynamical evolution is to turn to effective field theories (EFTs), such as relativistic hydrodynamics for bulk matter or potential non-relativistic QCD for the bound states of heavy quarks, so-called quarkonium [15]. An EFT allows us to systematically simplify the theory description by focusing on the relevant degrees of freedom at a certain energy scale only. This opens up the possibility to use quantities, computable on the lattice (e.g. the equation of state), to implement a dynamical evolution of strongly interacting matter in a heavy ion collision. For the case of heavy-quark pairs it has become possible to define and derive the concept of an in-medium potential [16, 17, 18, 19, 20, 21], which summarizes how the quarkonium interacts with the surrounding QGP. One central open question is how to use this in general complex valued potential to implement the microscopic evolution of the quarkonium system. In this paper we will provide a concrete example of how this may be achieved.

An understanding of the dynamical evolution of heavy quarkonium is a key factor in achieving insights into the properties of the hot matter created in heavy-ion collisions. The original proposal by Matsui and Satz, stating that the survival probability of quarkonia ($J/\psi$) may be used as a probe of QGP formation, points out one essential feature of deconfinement: Debye screening [22]. In the deconfined QGP phase, light quarks and gluons that carry color are liberated and move about, screening the color charges of the heavy quarks. However, in addition to the screened potential force, the QGP constituents exert two other kinds of forces on quarkonia, namely the drag force and random force (kicks). The latter two are both qualitatively and quantitatively different from the force acting in the vacuum.

Many different phenomenological approaches are currently used to describe quarkonium data. Some are based on Schrödinger equation with the in-medium complex potential [23, 24, 25, 26] and others based on kinetic Boltzmann equations with chemical reactions between a quarkonium and an unbound heavy quark pair in the QGP [27, 28, 29, 30, 31, 32, 33]. It remains, however, an open question how to derive the dynamical evolution of heavy quarkonium systematically from QCD.

As heavy quarkonium constitutes a quantum mechanical bound state immersed into a strongly interacting medium, we are forced to deploy a genuinely quantum mechanical description of its dynamics. This question of the evolution of a small quantum system “S” coupled to an environment “E” has been studied in detail in condensed matter physics. In that context the open-quantum-system approach has been developed and applied to the description of quarkonium in the QGP as well [34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. The ultimate goal then is to eventually simulate a quarkonium as an open quantum system in the QGP and extract information on the in-medium forces from experimental data of $\Upsilon$ and $J/\psi$ measured at the LHC and RHIC [2, 49, 50, 51, 3, 52, 53, 54, 55, 56, 57].

In the open quantum system approach [58], dynamical information of the small system “S” is encoded via the reduced density matrix $\rho_{\rm{S}}(t)$. The master equation for the evolution of the reduced density matrix can be derived by integrating out the environment degrees of freedom from the density matrix of the overall system. This process requires a set of approximations, whose validity rests on energy and time scale hierarchies between the system and the environment.

In this work, we introduce a novel master equation specifically for the relative motion of quarkonium, which turns out to be of “quantum Brownian motion” type. I.e., it exploits the fact that the system evolution is much slower than the environment. A general master equation for the quantum Brownian motion of quarkonia in the Lindblad form [59] has been first derived in [38]:

[TABLE]

where $H_{\rm{S}}$ denotes the Hamiltonian of the small system and $L_{i}$ the so called Lindblad operators. This master equation describes a Markovian time evolution and one can prove that it preserves the positivity of the reduced density matrix.

The initial derivation of the Lindblad equation for quantum Brownian motion of quarkonia has led to heightened interest and activity in finding Lindblad type master equations for quarkonium in various different regimes of energy and time scale separation. One pertinent example is the quantum optical master equation for gluo-dissociation dynamics. It describes the process in which a singlet bound state (quarkonium) absorbs a real gluon and turns into an octet unbound state. It is an interesting application and at the same time requires careful examination of the validity of open system approach111 In particular, the standard approximation schemes in the literature, such as the gradient expansion (quantum Brownian motion) and the rotating wave approximation (quantum optics) [58], are not applicable to the transition from a deeply bound singlet state to continuum of octet spectra by gluon absorption. It is expected that classical description by Boltzmann equation turns out to be applicable at time scale longer than the decoherence time scale as also discussed in the single heavy quark case [76]. .

In this paper, we derive the explicit Lindblad equation for the relative motion of heavy quarkonium and present numerical simulations of its quantum Brownian motion using the quantum state diffusion (QSD) method. With the explicit inclusion of dissipative effects we overcome a central limitation of the previously deployed stochastic potential approach [37, 39], which corresponds to the lowest order gradient expansion of the full Lindblad equation. The present study follows in the footsteps of our recent paper [40], where we analyzed the quantum Brownian motion of a single heavy quark. The main outcome was that quantum dissipation is important not only on long time scales such as the heavy quark equilibration but also in the early stages, if the initial heavy quark wave function is localized. The present paper extends this analysis to quarkonia. With a heavy quark pair, a scattering between a heavy quark and a gluon interferes with that between a heavy antiquark and the gluon. Therefore quantum Brownian motion for a heavy quark pair acquires nontrivial correlation between the pair, in addition to the potential force. Using simulations in one spatial dimension via the QSD method we find that

• •

The effective coupling of quarkonia to the QGP depends on the dipole size.

• •

Study of the late time steady state of the reduced density matrix reveals that it is indeed consistent with a thermal Boltzmann distribution.

• •

Inspection of the early stages reveals that the effects of quantum dissipation are also important there.

This paper is organized as follows. In Sec. II we briefly review the foundations of the open quantum system approach and introduce a novel Lindblad master equation for the relative motion of quarkonium. The explicit form of the Lindblad operators for the relative motion is provided and its physical meanings are discussed. We also discuss the quantum state diffusion approach, a method of implementing the Lindblad master equation in terms of an ensemble of stochastically evolving wave functions. In Sec. III we present the numerical setups and then in Sec. IV we present our numerical results from a one-dimensional simulation for different medium conditions and study the effect of dissipation. Section V provides a summary and outlook.

II Lindblad equation for a quarkonium in the QGP

II.1 Basics of open quantum system

In the open quantum system approach, we consider a composite system which consists of a small system “S” and an environment “E”. The unitary evolution of the degrees freedom of “S” and “E” together can be described by a hermitian Hamiltonian $H_{\text{tot}}=H_{\rm{S}}\otimes I_{\rm{E}}+I_{\rm{S}}\otimes H_{\rm{E}}+H_{\rm int}$, which induces time translations in the total density matrix of states $\rho_{\text{tot}}$.

When we are interested in the dynamics of the small system alone, it is desirable to work with an effective description written solely in terms of the system degrees of freedom. Therefore we turn to the reduced density matrix, defined by

[TABLE]

and its time evolution, usually referred to as the master equation,

[TABLE]

which is the central object of the open quantum system analysis. Here $\mathcal{L}$ refers to a superoperator that maps one reduced density matrix into another at different times and effectively describes how the environment couples to the system.

Positivity ($\forall{\alpha},\bra{\alpha}\rho_{\rm{S}}\ket{\alpha}\geq 0$), hermiticity ($\rho_{\rm S}=\rho_{\rm S}^{\dagger}$) as well as unitarity (${\rm Tr}_{\rm S}\rho_{\rm S}=1$) are the basic properties that allow a physical interpretation of $\rho_{\rm S}$ as a density matrix. It can be proven that if a Markovian master equation respects these conditions, it can be expressed in a particular form, the Lindblad form [59] 222 The Caldeira-Leggett master equation [74] is not in the Lindblad form; nevertheless it describes quantum Brownian motion well. The Lindblad master equation is, however, conceptually more useful. :

[TABLE]

Note that the Hamiltonian $H_{\rm eff}$ is not necessarily the same as the Hamiltonian $H_{\rm{S}}$ when the system is isolated.

In our study, the system “S” and the environment “E” correspond to a quarkonium and the QGP (light quarks and gluons), respectively, and the Lindblad operators $L_{i}$, as explained below, are labeled by a continuous momentum variable $\vec{k}$ and discrete color $a$, i.e. $i=(\vec{k},a)$.

II.2 Lindblad equation for quarkonium relative motion

We consider a system, where the energy density of the medium is high enough for the QCD coupling to be small and at the same time the temperature is still low enough, so that the heavy quark mass represents the largest energy scale $m_{Q}\gg T\gg\Lambda_{\rm QCD}$. In such a scenario, which may be realized in heavy-ion collisions at LHC, using well controlled expansions in $1/m_{Q}$ and the strong coupling $g$, the effective Hamiltonian and the Lindblad operators in a finite volume $L^{3}$ have been derived in the influence functional formalism in Ref. [38]:

[TABLE]

where $(\vec{x}_{Q},\vec{p}_{Q})$ and $(\vec{x}_{\bar{Q}},\vec{p}_{\bar{Q}})$ denote the position and momentum operators of the heavy quark and antiquark, respectively, and $t^{a}$ represent color matrices in the fundamental representation. The in-medium real-time dynamics is hence governed by two quantities, the real functions $V(\vec{r})$ and $D(\vec{r})$ (and its Fourier transform $\tilde{D}(\vec{k})$), which are defined by gluon two-point functions 333 In comparison to [38], we adopt $8T^{2}A=D$ for simplicity and opposite sign convention for $V(\vec{r})$ and $D(\vec{r})$ for later purposes.

[TABLE]

and given explicitly using the gluon self-energy in the hard-thermal loop (HTL) approximation,

[TABLE]

with the Debye screening mass $m_{\rm D}^{2}=(g^{2}T^{2}/3)(N_{\rm c}+N_{\rm f}/2)$ for QCD with $N_{\rm f}$ light flavors. The Lindblad operator describes collisions between heavy quarks and plasma particles. The operator $\exp[i\vec{k}\cdot\vec{x}_{Q}](t^{a}\otimes 1)$ shifts the heavy quark momentum by $\vec{k}$ and rotates its color. The operator with $\vec{k}\cdot\vec{p}_{Q}/4MT$ describes the recoil of the heavy quark in each collision and accounts for its dissipation. Note that the Lindblad operator without this dissipative term corresponds to the master equation of the stochastic potential model [37, 39].

The functions $V$ and $D$ in Eqs. (6) and (7) are closely related to the real and imaginary part of the complex in-medium potential $V_{Q\bar{Q}}$ derived from the Wilson loop at high temperatures [16, 17]. Within the HTL approximation we have

[TABLE]

Equation (5) thus provides an explicit microscopic answer to the question of how the real and imaginary part of the quarkonium in-medium potential govern the quarkonium real-time dynamics. While the real part drives the von-Neumann like part of the dynamics, it is the imaginary part that encodes how fluctuations induce the decorrelation of the quarkonium state from its initial state [37, 39], as it evolves over time. It is this decay in correlations and not the annihilation of the static quark pair what is encoded in the damping of the thermally averaged Wilson loop, described by its complex potential $V_{Q\bar{Q}}$.

Let us concentrate on the relative motion of the quarkonium and attempt to integrate out the center-of-mass motion. As we will see below, a complete decoupling of the center-of-mass momentum in general is not possible and it will affect the relative motion. Thus the procedure needs to be performed carefully. Introducing the center-of-mass and the relative coordinates,

[TABLE]

we can rewrite the effective Hamiltonian and the Lindblad operators as

[TABLE]

The Hilbert space for a quarkonium consists of a direct product of two Hilbert spaces:

[TABLE]

where $\mathcal{H}_{R}$ is the Hilbert space for the center-of-mass coordinate and $\mathcal{H}_{r}$ is that for the relative coordinate and for the color space. The density matrix for the latter is obtained by tracing out the center-of-mass coordinate

[TABLE]

To obtain the master equation for $\rho_{r}$, we need to calculate

[TABLE]

By expressing the Hamiltonian as $H_{\rm eff}=H_{\rm eff}^{(R)}\otimes 1+1\otimes H_{\rm eff}^{(r)}$, the first term is

[TABLE]

Next, by writing the Lindblad operators as

[TABLE]

we obtain

[TABLE]

At this point it is not yet possible to express the above equations as explicit expressions in terms of $\rho_{r}={\rm Tr}_{R}(\rho_{\rm Q\bar{Q}})$. In order to proceed, let us make the assumption that the center-of-mass momentum is fixed $\vec{P}=\vec{P}_{\rm CM}$. To be specific, we assume the reduced density matrix is

[TABLE]

We then obtain

[TABLE]

which tell us that for constant $\vec{P}_{\rm CM}$ the quantities $C_{\vec{k}a}(O_{r},\vec{P}_{\rm CM})$ take the role of Lindblad operators for $\rho_{r}$.

Let us summarize the Lindblad master equation for the relative coordinates and color space of quarkonium

[TABLE]

In this derivation, the center-of-mass momentum $\vec{P}_{\rm CM}$ is given as an external parameter and can depend on time. This explicit dependence provides a way to incorporate the relative motion of quarkonium traversing the QGP in the real-time description when the quarkonium velocity $\vec{P}_{\rm CM}/2M$ is small. We find that this effect on the quarkonium relative motion is mild (see Appendix A for details) and we fix $\vec{P}_{\rm CM}$ at $\vec{0}$ in all of the simulations in Sec. IV .

Let us mention how the above master equation differs from previous proposals in the literature. In a recent study [42], the effects of dissipation have been considered in a quarkonium master equation; however, the author did not derive a self-consistent Lindblad equation, and as a result additional terms needed to be added by hand. A Lindblad-like master equation in a weakly coupled setting has been derived in [43] with a focus of further simplifying the dynamics using semi-classical approximations. While our master equation relies on the weak-coupling approximation, we will avoid any further semi-classical approximations and instead implement the full quantum time evolution.

II.3 Quantum state diffusion

In quantum mechanics, we can either consider the time evolution of a system based on its density matrix or go over to a mixed state description in terms of wave functions. The change from one to the other is called stochastic unraveling of the master equation. Since evolving the density matrix in the position basis incurs high numerical cost, it is often advantageous to carry out simulations in the language of wave functions.

While the lowest order gradient expansion of the quarkonium Lindblad equation can be stochastically unravelled into unitary time evolution of wave functions based on a linear Schrödinger equation with a stochastic potential [37, 39], the full Lindblad equation requires a more sophisticated treatment. As has been worked out in detail in the quantum physics community, any Lindblad master equation may be unravelled stochastically via the “quantum state diffusion (QSD)” method [63].

In the QSD method, the density matrix $\rho_{\rm{S}}$ is obtained from the ensemble average of wave functions,

[TABLE]

and the wave functions evolve according to the following nonlinear stochastic Schrödinger equation,

[TABLE]

where $\rm M$ represents taking the ensemble average and $d\xi_{i}$ is complex white noise satisfying

[TABLE]

The nonlinearity arises from the terms containing the expectation value of the Lindblad operator with respect to the wave function $\langle L_{i}\rangle_{\psi}$. The above QSD equation can be shown to be equivalent to the following nonlinear stochastic Schrödinger equation for unnormalized wave functions

[TABLE]

which we employ in our numerical simulation. When we calculate the occupation number of some specific state $\phi_{i}(\vec{x})$ of a quarkonium, we thus calculate

[TABLE]

We would like to emphasize that with Eq. (23b) the full quantum time evolution of the in-medium quarkonium system has been cast into the form of a stochastic nonlinear Schrödinger equation. As our derivation originates in the influence functional treatment of the QCD path integral and deploys a well defined set of approximations, it is hence possible for the first time to provide a direct link between QCD and previous phenomenological models that introduce a non-linear Schrödinger equation in an ad-hoc manner.

III Numerical setup

Let us introduce our simulation prescription for the relative motion of a quarkonium in the QGP. The evolution equation for the wavefunction from the QSD method consists of two parts: The effective Hamiltonian term, as well as additional terms related to the Lindblad operators (including the stochastic ones). They represent a stochastic integro-differential equation, whose explicit form for three dimensions is shown in Appendix B.

For simplicity, in this study, we simulate the full dissipative dynamics in one spatial dimension and ignore the heavy quark colors. We solve the Hamiltonian term via the 4th order Runge-Kutta method and the Lindblad terms via a simple forward Euler time step.

Our goal is to showcase how dissipative effects influence the evolution of quarkonium states. Thus we wish to compare to simulations in which these effects are absent. Within the QSD framework this can be achieved by discarding all terms (except for the kinetic energy) that are not finite in the $T/M\to 0$ limit. This leads to the following evolution equation:

[TABLE]

which is equivalent to the master equation of the stochastic potential model [39] in spite of rather different appearances. The former is nonlinear while the latter is linear in the wave function.

The parametrization for the functions $V(x)$ and $D(x)$ are given as follows:

[TABLE]

From perturbative calculations (7), we have

[TABLE]

so that we choose $\alpha=0.3$, $m_{\rm D}=2\ell_{\rm corr}^{-1}=2T$ 444 $m_{\rm D}\simeq{2}/{\ell_{\rm corr}}$ in (27) is obtained by equating the full width half maximum of $D(r)$ in (7) and that in (26). , and $\gamma=T/\pi$. Since the one-dimensional Coulomb force is singular at the origin, we regularize the potential by a cutoff $x_{c}=1/M$ corresponding to the validity of the nonrelativistic approximation $|p|<M$. In the simulations in the next section, the temperature of the QGP is chosen to be $T/M=0.1$-$0.3$ or $T(t)=T_{0}\cdot[t_{0}/(t_{0}+t)]^{1/3}$ with $t_{0}=0.84\ {\rm fm}$ and $T_{0}=0.47\ {\rm GeV}$ for the case of Bjorken expansion. The center-of-mass momentum is set to $P_{\rm CM}=0$ as the dependence of the quarkonium dynamics on its values is found to be small (see Appendix A for the details).

To simulate quarkonium physics with this setup, we discretize space and time by $\Delta x=1/M$ and $\Delta t=0.1M(\Delta x)^{2}$. The spatial discretization $\Delta x$ is chosen much finer than the typical medium length scales $m_{\rm D}^{-1}\sim\ell_{\rm corr}/2\sim 1/2T$. We also take the system volume $L=N_{x}\Delta x=254\Delta x$ much larger than the medium length scales $m_{\rm D}^{-1}\sim\ell_{\rm corr}/2\sim 1/2T$. These parameters and the lattice setup are summarized in Table 1.

Finally, let us remark on the boundary conditions used for the noise field $d\xi(x)$. To simulate on a finite size lattice $[-L/2,L/2]$, we impose periodic boundary conditions on the wave function $\psi(x)$. Then the boundary condition for the noise field is given by

[TABLE]

requiring that the noise field $d\xi(x)$ obeys a periodicity of $L/2$. Correspondingly, the noise correlation function $D(x-y)$ for $-L/4\leq x,y\leq L/4$ should be interpreted as $D(r_{xy})$:

[TABLE]

In this way, we redefine the function $D(x)$ for a finite system with an interval $-L/2\leq x\leq L/2$, which is all that is needed for solving the QSD equation.

IV Numerical results

We study and simulate how quantum dissipation and center-of-mass motion of a quarkonium affect the dynamical evolution of its relative motion (see Appendix A for the latter). We first show in Sec. IV.1 how quantum dissipation influences the time evolution. We numerically confirm that in the presence of quantum dissipation the quarkonium system approaches a steady state at late times that is independent of the initial conditions from which the evolution commenced. We furthermore find that the late time distribution is very close to the Boltzmann distribution, whose slope is well within $10\%$ of the input temperature of the surrounding medium. The quarkonium system thus appears to approach genuine thermal equilibrium at late times, assisted by the interplay of quantum fluctuations and dissipation.

We then show in Sec. IV.2 how the evolution depends on temperature and heavy quark mass by comparing with simulations with $T/M=$ 0.1 and 0.3. By rescaling the evolution time $t$ by the heavy quark relaxation time, we find that a similar relaxation behavior shows up in both of the cases. Finally in Sec. IV.3, we present the quarkonium evolution in a time-dependent background with temperature $T(t)$, which decreases according to Bjorken expansion. To make closer contact with experimental data, we project the quarkonium wave functions onto the eigenstates in Cornell potential.

IV.1 Equilibration with quantum dissipation

We first simulate at $T/M=0.1$ and compare the two cases where the initial condition is either given by the ground state or the 1st excited state of the Hamiltonian

[TABLE]

Note that the Hamiltonian $H_{\rm Debye}$ here is different from the effective Hamiltonian $H_{\rm eff}$ in (19b).

In Fig. 1, we plot the time evolution of the occupation number $N_{i}(t)$ of the $i$-th eigenstate for the ground state ($i=0$) and the 1st ($i=1$) excited state of $H_{\rm Debye}$ for two different initial conditions $N_{0}(0)=1$ and $N_{1}(0)=1$. We find that each of the occupation numbers relaxes to a constant value irrespective of the initial condition and hence confirm that the dissipative effects lead to a steady state. This behavior is indicative that the relative motion of a quarkonium in the QGP becomes equilibrated.

Next we plot in Fig. 2 the occupation numbers of the lower levels at late times $Mt=4650$, well within the steady state regime, as a function of the eigenenergy $E_{i}$ of the Hamiltonian $H_{\rm Debye}$. We also show the results for the same analysis with $T/M=0.3$ at time $Mt=900$. The distribution can be fitted by the Boltzmann distribution $\propto\exp(-E_{i}/T_{\rm{fit}})$ with $T_{\rm fit}/M=0.099\pm 0.004$ for $T/M=0.1$ and $T_{\rm{fit}}/M=0.288\pm 0.013$ for $T/M=0.3$. The fitting range is limited to the eigenstates $\phi_{i}$ with velocity $\frac{\sqrt{\langle p^{2}\rangle_{\phi_{i}}}}{M/2}<0.5$. 555 We expect that the steady state of the Lindblad evolution is the Boltzmann distribution only when the velocity is small. See Appendix C for more details.

This corresponds to fitting the lowest 21 levels for both $T/M=0.1$ and $0.3$.

As we have seen above, the relative motion is equilibrated with quantum dissipation. To see the importance of the quantum dissipation, we switch off the $\mathcal{O}(T/M)$ terms, i.e. keeping only those terms explicit in (25), and compare with the full simulation. The comparison is made at $T/M=0.1$ and is shown in Fig. 3. We can see the clear differences not only in the asymptotic behavior on long time scales but also in the initial behavior. When the dissipation is switched off, the initial decay of the ground state occupation is faster and approaches a much smaller value than in the case with dissipation. Physically, the drag force prevents the heavy quark pair from dissociating in the QGP and balances with the thermal fluctuations to maintain the system in equilibrium. Since we observe clear dissipative effects in the initial decay, we conclude that it may not be ignored even within the finite QGP lifetime $\sim 10$ fm. We will come back to this issue in a slightly more realistic setup in Sec. IV.3.

IV.2 Dependence of the temperature and the heavy quark mass

Here we study how the time evolution of the occupation numbers depends on the temperature. We compare the results with $T/M=0.1$ and 0.3 starting from the ground state of $H_{\rm Debye}$ in Fig. 4 (a). For bottomonium, this corresponds to comparing $T\simeq 0.47$ GeV and 1.41 GeV, respectively. As we can see in the figure, the relaxation takes place much faster at higher temperature $T/M=0.3$ because of two physical effects: (i) The heavy quark damping rate is larger, and (ii) the ground state wave function is more extended and receives decoherence more easily. To cancel out the first effect, we plot the time scale in the unit of heavy quark damping rate $\tau_{\rm{eq}}\equiv MT\ell_{\rm{corr}}^{2}/\gamma=\pi(M/T)^{2}/M$ in Fig. 4 (c). There is still some difference between $T/M=0.1$ and 0.3, which we ascribe to the other effect: The decoherence of the ground state wave function. The decoherence rate for a wave function of size $\ell_{\psi}\equiv\sqrt{\langle x^{2}\rangle_{\psi}}$ is estimated as $\tau_{\rm dec}^{-1}=D(0)-D(\ell_{\psi})$, which amounts to $\tau_{\rm dec}\simeq 456/M$ and $19/M$ for the ground state at $T/M=0.1$ and 0.3, respectively. By rescaling to $\tau_{\rm eq}$, we get $\tau_{\rm dec}/\tau_{\rm eq}=1.45$ and 0.55 for $T/M=0.1$ and 0.3, respectively, which qualitatively explains the reason why the initial decay for $T/M=0.3$ is faster even after rescaling to $\tau_{\rm eq}$666 To be strict, the inclusion of the dissipation changes the initial decay rate from the estimate by the decoherence rate as we saw in the previous section. .

We can also interpret $T/M=0.1$ and 0.3 as a bottomonium and a charmonium at $T=0.47$ GeV, respectively. As shown in Fig. 4 (b), we find the relaxation of a bottomonium proceeds more slowly than that of a charmonium again with the same two physical effects as above.

IV.3 Phenomenological implication to heavy ion collision experiments

In order to relate these simulations to quarkonia in heavy ion collisions, we take account of the expansion of the QGP by solving the QSD equation in a time-dependent environment undergoing Bjorken expansion:

[TABLE]

We define occupation numbers by projecting the wave functions onto the eigenstates of the Hamiltonian with the vacuum Cornell potential,

[TABLE]

so that the quarkonium yield is related to these occupation numbers at kinetic freezeout.

Note that we use the Debye screened potential for $V(x)$ in the QSD evolution equation. Since the initial quarkonium wave function is not understood well and still an open issue, we adopt the eigenstates of $H_{\rm Cornell}$ in this study. We take the ground state and the 1st excited state as the initial conditions and show in Fig. 5 how the occupation numbers of these states evolve in time. From Fig. 5, we find that for both bottomonium and charmonium the occupation number of the 1st excited state decays faster than that of the ground state, which supports the phenomenological idea of sequential melting. Also the charmonium ground state/1st excited state decays faster than the bottomonium ground state/1st excited state. It is expected because bottomonium is more localized so that the decoherence is more ineffective and because bottomonium is heavier and thus the relaxation time is longer.

To see the effect of quantum dissipation, we also simulate without the dissipative terms as we did in Sec. IV.1 and plot in Fig. 5. Since the temperature at around 18 fm is about the transition temperature 170 MeV, 777 The typical QGP lifetime is $\sim 10~{}{\rm fm}$ in full three dimensional hydrodynamic simulations. the dissipative effects on the relative motion of a quarkonium in the QGP are up to 20% effect for the ground states and marginally effective for the 1st excited state. The reason why dissipation affects the ground state more is that the decoherence is ineffective for localized states and that the relative importance of the dissipation enhances as we found in [40]. By comparing with Fig. 3, one may notice that the dissipative effects appear to be much less important here in Fig. 5. However, when making a comparison, one needs to consider the temperature decrease for the latter because the decoherence and damping proceed much slower at lower temperatures.

V Summary

In this paper, we study the quantum Brownian motion of a heavy quark pair in the quark-gluon plasma (QGP). Treating the heavy quark pair as an open quantum system in the QGP, we derive and solve the master equation for the relative motion among the pair. In this master equation, three different forces govern the evolution of the heavy quarks: Debye screened potential force, thermal noise, and drag force. Thus the classical counterpart of our quantum mechanical description is given by the Langevin equation of two Brownian particles interacting with each other through the Debye screened potential. It is for the first time that a full quantum mechanical simulation including dissipation is performed for heavy quarks in the QGP.

The master equation of the heavy quark pair takes the Lindblad form (19), which ensures the basic physical properties of the reduced density matrix: positivity, hermiticity, and unitarity. Among the three different forces stated above, the Lindblad operator is responsible for the thermal noise and the drag force. The thermal noise is represented by a momentum shift operator while the drag force by the recoil of heavy quarks during each microscopic collision. The Lindblad operator for the relative motion of the heavy quark pair depends on the center-of-mass momentum, but its dependence is found to be small so that we only considered the static case ($P_{\rm CM}=0$) in our simulations.

We solved the Lindblad master equation by a stochastic unravelling method called quantum state diffusion (QSD). Using this method, any Lindblad master equation is shown to be equivalent to a nonlinear stochastic Schrödinger equation by which we can correctly produce a mixed state ensemble for the density matrix. In our numerical simulation in one dimension, we first checked basic properties of the master equation (Sec. IV.1). The occupation number of eigenstates relaxes toward a value independent of the initial condition, and the steady state distribution is consistent with the Boltzmann distribution. We also studied the effect of quantum dissipation by comparing a simulation without dissipation. The quantum dissipation delays the relaxation toward equilibrium, which is consistent with our intuitive classical picture that the drag force prevents a heavy quark pair from dissociating. We next simulated the temperature and heavy quark mass dependences on the time evolution (Sec. IV.2). The time evolution strongly depends on these quantities but is found to scale to a considerable extent with the heavy quark damping time. Finally, we take into account the expansion of the QGP in relativistic heavy ion collisions and solved the master equation with a time-dependent temperature (Sec. IV.3). By simulating the bottomonium and charmonium yields in the QGP undergoing the Bjorken expansion, we found that charmonium dissociates faster than bottomonium and that the excited states decay faster than the ground states. We also found that the quantum dissipation delays the ground state dissociation compared to the case without dissipation while the excited state dissociation is totally insensitive to the quantum dissipation. This difference comes from the fact that for an extended state the decoherence due to thermal fluctuation is the driving force of the dissociation while for a localized state (such as the ground state) the decoherence is ineffective and as a result the dissipation becomes as important as the decoherence.

Through this analysis, it becomes clear that quarkonia probe fundamental length scales of the QGP: the screening length $1/m_{\rm D}$ and the correlation length $\ell_{\rm corr}$ as shown in (26). The screening length has been evaluated via the real-part of the heavy quark potential (see e.g. [68, 69, 70]) and the heavy quark free energies (see e.g. [71, 72]) in various lattice QCD simulations. On the other hand, the correlation length is related to the imaginary part, whose determination from the lattice is much less robust. Using heavy quark observables in heavy-ion collisions such as single leptons and open heavy flavor mesons, the single heavy quark damping rate has been determined phenomenologically with some accuracy. The corresponding quantity in the quarkonium Lindblad equation is given by $\tau_{\rm eq}=MT\ell_{\rm corr}^{2}/\gamma$. Therefore, the suppression pattern of quarkonium yields provides a way to determine the correlation length of colored excitations $\ell_{\rm corr}$, a fundamental dynamical quantity of the QGP which has not yet been calculated precisely on the lattice. To perform such a comprehensive study in the future, we need to implement a more realistic computation, such as a three-dimensional simulation on a three-dimensional hydrodynamic background and include the color of the heavy quarks, to name a few.

Acknowledgements.

T.M. is grateful to Shiori Kajimoto for valuable discussions. The work of Y.A. is supported by JSPS KAKENHI Grant Number JP18K13538. M.A. is supported in part by JSPS KAKENHI Grant Number JP18K03646. This work has utilized computing resources provided by UNINETT Sigma2 - the National Infrastructure for High Performance Computing and Data Storage in Norway under project NN9578K-QCDrtX ”Real-time dynamics of nuclear matter under extreme conditions”

Appendix A Dependence of the center-of-mass motion

Here we show the effects of the center-of-mass motion on the relative motion of a quarkonium in the QGP at $T/M=0.1$. We set the center-of-mass momentum at different values $P_{\rm{CM}}/M=0,0.4,0.6,0.8$, and $1.0$, which are equivalent to the center-of-mass velocity $v_{\rm{CM}}=0,0.2,0.3$, $0.4$, and $0.5$, respectively. The time evolution for the ground state occupation numbers is shown in Fig. 6. The result shows that the effect of $P_{\rm CM}$ is rather small. In Fig. 7, we also show the occupation number for each eigenstate at late enough time ($Mt=4650$) when the system has reached a (non-equilibrium) steady state. Again, the result is quite insensitive to the value of $P_{\rm CM}$.

Appendix B Explicit form of QSD equation

We present here the explicit form of the QSD equation (23a) with heavy quark colors ignored:

[TABLE]

with

[TABLE]

Appendix C Equilibrium distribution in the classical limit

The explicit form of the quantum master equation for a heavy quark and antiquark pair has been given in [38]. By considering up to the 2nd order in the derivative expansion for the time coarse graining, we obtain the influence functional $S_{\rm IF}=S_{\rm pot}+S_{\rm fluct}+S_{\rm diss}+S_{\rm L}$, where each term corresponds to a specific power counting in the heavy quark velocity $v$ as $S_{\rm pot}\sim S_{\rm fluct}\sim v^{0}$, $S_{\rm diss}\sim v$, and $S_{\rm L}\sim v^{2}$. Let us ignore $S_{\rm L}$, which is justified when the velocity is small $v\ll 1$, and take the classical limit in the master equation. For simplicity, we do not consider color degrees of freedom here. Then, in contrast to the fully quantum system considered in the main text, the master equation for the center-of-mass motion can be integrated. The result is

[TABLE]

with

[TABLE]

Using the coordinates

[TABLE]

it is written as

[TABLE]

In the classical limit, we can prove that the equilibrium distribution is the Boltzmann distribution. By performing the Wigner transformation

[TABLE]

and taking the classical limit $\hbar\to 0$ 888 When the dimension of $\hbar$ is recovered, the kinetic term becomes $\nabla_{R}\nabla_{r}\hbar$, potential $V/\hbar$, fluctuation $F_{1}/\hbar^{2}$, and dissipation $F_{2}$ is unchanged. , we obtain the classical kinetic equation:

[TABLE]

Using the approximation $\vec{F}_{2}(x)=\vec{\nabla}F_{1}(x)/4MT$ employed in (48), the equilibrium distribution is

[TABLE]

If we do not use this approximation, the correction is estimated as

[TABLE]

and the effective temperature slightly deviates from the environment temperature,

[TABLE]

which gives a rough explanation of the difference between $T_{\rm fit}$ and $T$ in Fig. 2.