Quantum dissipation of a heavy quark from a nonlinear stochastic Schr\"odinger equation

We study the open system dynamics of a heavy quark in the quark-gluon plasma using a Lindblad master equation. Applying the quantum state diffusion approach by Gisin and Percival, we derive and numerically solve a nonlinear stochastic Schr\"odinger equation for wave functions, which is equivalent to the Lindblad master equation for the density matrix. From our numerical analysis in one spatial dimension, it is shown that the density matrix relaxes to the Boltzmann distribution in various setups (with and without external potentials), independently of the initial conditions. We also confirm that quantum dissipation plays an essential role not only in the long-time behavior of the heavy quark but also at early times if the heavy quark initial state is localized and quantum decoherence is ineffective.


Introduction
The study of phases of strongly interacting matter under extreme conditions [1,2] is attracting broad interest in recent years, extending far beyond nuclear matter. Fruitful interdisciplinary collaborations, among others, with the field of ultracold atomic gases (for a review see e.g. [3]) have broadened the scope of how to explore the origins of the universe.
One example of extreme conditions is very high temperatures, which are particularly interesting, since they bear direct relevance to the birth and early history of our universe. The properties of the hottest matter ever created on the Earth are being investigated at current collider facilities, such as the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) and upcoming facilities, such as NICA and FAIR. There, two heavy nuclei are collided at ultra-relativistic energies in order to create a novel deconfined state of matter known as the quark-gluon plasma (QGP) [4].
The transition from hadronic matter to the QGP is characterized by the liberation of colored degrees of freedom (quarks and gluons) otherwise confined inside hadrons. This qualitative picture is confirmed by lattice QCD calculations of, for example, the QCD entropy, which show a significant rise around the pseudocritical temperature [5][6][7][8].
In analogy with the Debye screening phenomenon in electromagnetic plasmas, the liberated colored particles may rearrange themselves around a test color charge such that the test charge is screened with a finite screening length [9]. Compared to the confining stringlike force between color sources in the vacuum, the in-medium force in a high temperature QGP is hence drastically modified and becomes short ranged [10][11][12].

JHEP07(2018)029
The in-medium modification of the QCD force has been expected to have dramatic consequences on the behavior of bound states of heavy quarks and antiquarks, so called heavy quarkonium [13,14]. In turn the dynamics of quarkonium states observed in heavyion collisions promises a direct window on the phase structure of the bulk matter, in which they are immersed.
However, the enhanced production rate of charm quark pairs at LHC complicates an interpretation of J/ψ yields at LHC [26,27]. The reason is that now another process to create J/ψ needs to be taken into account, i.e. the non-negligible probability that initially uncorrelated charm quarks end up forming bound states at the phase boundary [28], i.e. in the late stages of the collision at freezeout. The successful prediction of J/ψ yields at the LHC by means of the statistical model of hadronization [29] is seen as support for the existence of such a production mechanism.
At present it is still not clear at which collision energy the two effects, i.e. the suppressed and enhanced yields of quarkonium become comparable. In order to more clearly interpret the collected experimental results we thus need to better understand the dynamics of quarkonia in the QGP. I.e. the development of a unified quantum mechanical description for the real-time equilibration of quarkonia is called for.
In the open quantum system formulation [65], we distinguish the subsystem of interest (quarkonium) from the environment (QGP), which is made possible by a hierarchy of time scales in each sector. The dynamics of the environment is fast enough, so that we can trace it out and replace its coupling to the subsystem by the average response to a slowly changing external source.
In this way a master equation for the reduced density matrix of a heavy quark pair can be obtained, which can be expressed as arising from the contribution of three kinds of forces: the screened potential, thermal fluctuations, and dissipation [48]. The first two forces combine into a fluctuating potential force (stochastic potential [46,47,50]), while the last one is an irreversible force. So far there exist only a few numerical analyses of the quantum dissipation of quarkonium in the QGP [55][56][57]. However, quantum dissipation is an essential ingredient to understand the long-time behavior of a heavy quark pair in the QGP. It is particularly important to know how and when quantum dissipation influences the time evolution of the heavy quark pair because the lifetime of the QGP in heavy-ion collisions is not long enough for the heavy quark pair to get fully equilibrated.
In this paper we consider as a first step the physics of a single heavy quark immersed in a hot medium. We numerically solve the corresponding master equation and study its JHEP07(2018)029 equilibrium solution as well as the effects of dissipation. The central conceptual result of this study is a stochastic unravelling prescription for the master equation, in which the wave function is evolved in terms of a stochastic Schrödinger equation and the mixed states of the density matrix emerge from an ensemble average.
Our results are obtained following the "Quantum State Diffusion" approach developed by Gisin and Percival [66,67] and by subsequently solving the corresponding nonlinear stochastic Schrödinger equation. Note that by applying the Quantum State Diffusion approach to a heavy quark master equation derived in the Lindblad form [48,68] (see section 2 for the definition) we obtain for the first time in this context a non-linear Schrödinger equation with a clear connection to the underlying microscopic theory.
For simplicity, we consider a heavy quark in the QGP in one spatial dimension. It is then shown that the master equation possesses a steady state solution consistent with the Boltzmann distribution ρ eq ∝ e −βH . Furthermore, we observe that the approach to equilibrium depends on initial conditions and cannot be captured by a single decay rate. Finally, we analyze the effect of quantum dissipation by comparing with simulations in which the dissipative terms are dropped. Their effect, as expected, is essential at later times but interestingly already plays an important role at rather early times if the initial wave function is well localized and decoherence by thermal fluctuations is ineffective.
This paper is organized as follows. In section 2, we introduce the Quantum State Diffusion approach applicable to a general Lindblad master equation. We then apply this approach to the Lindblad equation for a single heavy quark in the quark-gluon plasma at high temperature and derive the corresponding nonlinear stochastic Schrödinger equation. In section 3, we solve this nonlinear stochastic Schrödinger equation numerically and analyze the relaxation process of the heavy quark. In addition we study the effect of quantum dissipation on the time evolution of a heavy quark, before we summarize our work in section 4.

Lindblad equation and quantum state diffusion
There is a particularly useful class of master equations for open quantum systems, which is Markovian and fulfills basic physical requirements: the reduced density matrix ρ is hermitian (ρ = ρ † ), correctly normalized (Trρ = 1), and positive ( α|ρ|α ≥ 0 for any state |α ) during its time evolution. Such master equations may be written in general as in the so called Lindblad form [68]. Here the evolution of the density matrix operator is described in terms of the full dimension of the system Hilbert space. This makes a direct numerical simulation computationally highly demanding, especially in realistic 3+1 dimensions. There are, however, several ways to solve the Lindblad equation by what is known as stochastic unravelling. I.e. by carrying out a stochastic evolution of the wave functions of the system instead, whose ensemble average then correctly reproduces JHEP07(2018)029 the density matrix. One such stochastic unravelling corresponds to the quantum state diffusion (QSD) approach [66]. Given the Lindblad equation (2.1), the corresponding QSD equation is a stochastic nonlinear Schrödinger equation: with complex white noises dξ n whose mean and variance are given by Here O ψ ≡ ψ|O|ψ denotes the quantum expectation value of an operator with respect to a state ψ and M(O) denotes the statistical average of O. This stochastic evolution equation is solved in the Itô discretization scheme. By using ψ(t) everywhere in eq. (2.2) we obtain the wave function at the next discrete time step ψ(t + dt). In the limit dt → 0, the QSD equation (2.2) preserves the norm of ψ in each stochastic update. The initial wave function is distributed according to the initial mixed (or pure) state of the density matrix. The density matrix is subsequently constructed from an ensemble average of wave functions ψ(t), and obeys the Lindblad equation (2.1) in the dt → 0 limit. The QSD equation can also be formulated for unnormalized wave functions φ with the same complex noise. Here O φ ≡ φ|O|φ / φ|φ denotes the quantum expectation value. The density matrix constructed by is again a solution of the master equation (2.1). In our numerical simulation, we implement eq. (2.5).

Quantum state diffusion for a heavy quark in the quark-gluon plasma
Let us now consider the theory of open quantum systems for a single heavy quark in the QGP and stochastically unravel its Lindblad equation via the QSD approach. The Lindblad equation for heavy quarks has been derived in [48] by treating the scattering JHEP07(2018)029 between heavy quarks and medium particles perturbatively, i.e. by assuming that the QCD coupling constant g is small. It is further assumed that the heavy quark mass M is much larger than the temperature T /M 1 so that there exists a time scale hierarchy between heavy quarks and medium particles.
The Lindblad master equation for a single heavy quark is given by the following operators [48] with N c and N f being the numbers of colors and quark flavors, respectively. The Lindblad operator L k describes the scattering process between a heavy quark and medium particles with momentum transfer k, taking place with rateD(k). The term ∝ e ik·x in L k describes thermal fluctuations, while the term ∝ e ik·x/2 ik·∇ 4M T e ik·x/2 describes dissipation and originates in the recoil of the heavy quark during the collision. For simplicity, we ignore the effects of internal color degrees of freedom. 1 We can calculate the parts of the QSD equation as follows. The nonlinear term is given by (2.8) where n φ and j φ denote the probability density and current: and f and g i are operators defined as The function D(x) in the equations above is the inverse Fourier transform ofD(k). The linear deterministic term reads

JHEP07(2018)029
and the linear stochastic term amounts to where the definition and the correlation of the complex noise field dζ(x) is given by

Results of numerical simulation
In this section we explicitly check the application of the QSD approach and investigate the properties of the Lindblad equation in three simple settings. We consider a single heavy quark in the QGP in one spatial dimension either with or without external potentials.
In particular, we study the equilibration of the heavy quark in each setting and discuss the importance of quantum dissipation. As external potentials we deploy either the harmonic potential or the regularized Coulomb potential: where r c = 1/M . The noise correlation function D(x) is set to have correlation length ∼ 1/m D and is approximated with a Gaussian dependence on distance The parameters of our numerical setup are summarized in table 1. The Hamiltonian time evolution is solved by the fourth-order Runge-Kutta method and that for the other parts is implemented via an explicit forward step according to the QSD equation (2.5) with dt = ∆t. We check that the discretization effects for both ∆t and ∆x are negligible by comparing with results with smaller ∆t or ∆x. In the simulation periodic boundary conditions are employed so that the noise correlation is replaced with M (dζ(x)dζ * (y)) = D(r xy )∆t, (3.3a)

Equilibration of a heavy quark
The Lindblad equation itself does not guarantee the Boltzmann distribution ρ ∝ exp(−H/T ) to be the static solution (see appendix A). In the derivation of the Lindblad equation, the fluctuation-dissipation theorem for the environment, i.e. the QGP sector, constrains the terms implementing the fluctuation and dissipation of the heavy quarks.
To be specific, the coefficient i/4M T in L k is determined from the fluctuation-dissipation theorem for a thermal QGP medium. Therefore, we expect that the equilibrium density matrix is close to the Boltzmann distribution ρ ∝ exp(−H/T ). Here we analyze how well the equilibration of the density matrix is achieved and study its equilibrium properties.

In the absence of an external potential
The initial heavy quark wave function is taken to be uniform (plane wave with zero momentum). In figure 1, we show the profiles of a normalized wavefunction at different times in one sample event. A wavefunction type typically encountered here is a localized solitonic state, arising from the nonlinearity of the evolution equation. Figure 2 on the other hand contains the time evolution of the momentum distribution of the heavy quark:  The corresponding classical dynamics of the heavy quark is a Brownian motion with a drag force dp dt = − γ M T l 2 corr p. (3.5) Its typical relaxation time is τ relax = M T l 2 corr γ = 100π/M . We can see that the momentum distribution approaches the Boltzmann distribution with temperature T = 0.1M over a time scale ∼ τ relax . Note also that at late times (t = 620/M, 1240/M ) slight deviations from the Boltzmann distribution are observed above p 1.5M . The reason lies in the poor convergence of the gradient expansion, applied in evaluating the Lindblad operators, at high momenta (p M ), when one takes into account the effects of dissipation. On the other hand, we should not rely on the nonrelativistic description for a heavy quark with such a high momentum, so this limitation is not so restrictive in practice.

In the presence of external potentials
Let us turn to the case with external potentials present next. The potential here is added by hand out of pure theoretical interest and should not be confused with the potential for a quarkonium state. One should rather think of a single heavy quark in a fictitious trap.
The initial heavy quark wave function is taken to be the ground state, the first, or the second excited states of the corresponding Hamiltonian. The time evolution of the occupation number of these levels, is shown in figure 3 for the harmonic potential and in figure 4 for the regularized Coulomb potential. Independent of the initial conditions, the occupation numbers converge to their equilibrium values. In contrast, the relaxation time depends on the initial condition. One might expect a naive relaxation process  to describe the dynamics, as motivated and applied in the rate equation approach to heavy quarks. However, from figure 3 and figure 4 it is obvious that a single decay rate cannot capture the relaxation of the eigenstate occupation. Note that actually there are cases where the occupation number shows a non-monotonic approach to equilibrium. In order to investigate the properties of the equilibrium density matrix, we show in figure 5 the equilibrium distribution of the lowest ten levels as a function of the eigenenergy for the harmonic and the regularized Coulomb potentials. In the figures, we plot the results for several different potential parameters: ω/M = 0.01, 0.04, 0.09 for the harmonic potential and α = 0.2, 0.3, 0.4 and r c = 1/M for the regularized Coulomb potential. The initial condition is chosen to be the ground state of each Hamiltonian. The equilibrium distribution is calculated at late enough time for each setup: M t = 1550, 3100, 4650 for ω/M = 0.01, 0.04, 0.09 and M t = 4650, 7750, 9300 for α = 0.2, 0.3, 0.4, respectively. In all cases, the distribution is close to the Boltzmann distribution with T = 0.1M . We also calculate the real and imaginary parts of the off-diagonal elements of the density matrix in equilibrium and check that they are consistent with zero within the statistical uncertainty.

The effect of dissipation
As a last consideration let us now turn off the quantum dissipation. Since quantum dissipation is described by the term ∝ e ik·x/2 ik·∇ 4M T e ik·x/2 , we can switch it off by taking the  In figure 7, we show the time evolution of the eigenstate occupation as given by the QSD equation without dissipation. The lowest three levels become equally occupied regardless of the energy gaps. 2 It is expected that not only these levels but all the levels are eventually occupied equally if the quantum dissipation is neglected. We also observe that the effect of quantum dissipation at early time strongly depends on the external potential. This dependence can be understood by analyzing the Lindblad equation as follows. The initial decay rate of an eigenstate ψ i of the Hamiltonian is given by  Using the Lindblad operators in eq. (2.7), we obtain 9) in which the first (second) line arises from thermal fluctuations (dissipation). With the Gaussian approximation (3.2) for D(x), we get for our one-dimensional model (3.10) With the values from table. 1, the effect of dissipation (the second line in eq. (3.10)) is −0.048 − 0.125 ∇ 2 ψ i /M 2 . The ground state wave function yields ∇ 2 ψ 0 /M 2 0, −0.02, −0.08 for the free case, the harmonic potential, and the regularized Coulomb potential, respectively. Therefore, the effect of dissipation in this case ranges from −0.048 to −0.038 and only slightly depends on the potential.
On the other hand, the decay rate due to thermal fluctuation is quite sensitive to the size of the wave function as can be seen from the first line in eq. (3.10). In fact, the wave function size is M ∆x 3.5, 1.9 < 10 = M l corr for the ground states of the harmonic and the regularized Coulomb potentials respectively. If the size is smaller than the correlation length l corr , the decay rate in the absence of dissipation is already comparatively small so that the relative importance of dissipation increases.

JHEP07(2018)029 4 Conclusion
In this paper, we investigate how quantum dissipation influences the time evolution of the density matrix of a heavy quark in the quark-gluon plasma (QGP). The master equation for heavy quark systems in the QGP has been obtained in the Lindblad form [48,68] and thus possesses particularly useful properties: the density matrix ρ stays hermitian, remains correctly normalized and positive during the time evolution. We solve this Lindblad equation for a single heavy quark by stochastic unraveling. Applying the approach of quantum state diffusion (QSD) [66] to our Lindblad equation, we derive a nonlinear stochastic Schrödinger equation for the heavy quark wave function.
Subsequently we solve the QSD equation for a heavy quark in simple settings, i.e. in one spatial dimension and with or without external potentials (harmonic and regularized Coulomb potentials). We found that in both cases the density matrix relaxes to ρ eq ∝ e −H/T within statistical errors. This property is expected from a constraint in the Lindblad equation introduced by the fluctuation-dissipation theorem for the QGP sector but was not explicitly guaranteed by the Lindblad equation itself (see appendix A). We also found that the relaxation process strongly depends on the initial condition so that it is not captured by a simple rate equation.
As a further topic we study the effect of quantum dissipation by switching off the dissipative terms in the QSD equation. Without the dissipative terms, the heavy quark is overheated because it only receives energy from the thermal medium which is not dissipated back. It is shown that the importance of dissipation, as compared to the thermal fluctuations, strongly depends on the wave function size. The relative importance of dissipation increases when the wave function is small.
In the future, we plan to extend our analysis to the description of heavy quarkonium in the QGP. In that case, not only thermal fluctuations but also dissipation takes place nontrivially because the collisions of a heavy quark and those of a heavy antiquark interfere with each other. It would be interesting and phenomenologically relevant to study the effects of quantum dissipation in that case. The computations for quarkonium in three spatial dimensions and in evolving fluid background for heavy-ion collisions will be one of the ultimate goals of our project.

JHEP07(2018)029
A Approximate steady state solution of the Lindblad equation First we give a few algebraic relations of the Lindblad operators in eq. (2.7): (A.1c) The steady state solution ρ eq (p) of the Lindblad equation for a heavy quark without external potential V ext = 0 satisfies In a collision with momentum transfer k, the on-shell energy of the heavy quark changes by an amount, In the momentum representation, L k L † k and L † −k L −k are given by (A.4b) Then detailed balance in reactions between p ↔ p − k dictates ρ eq (p) ρ eq (p − k) The above approximation is reliable for ∆E k /T 1. For p ∼ M v, ∆E k /T ∼ v/T l corr + 1/M T l 2 corr ∼ v +0.1 so that our approximation is valid for v 1. Therefore, it is surprising that the numerical result in figure 2 is well approximated by the Boltzmann distribution for as large as p ∼ 1.5M .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.