Quantum evolution of quarkonia with correlated and uncorrelated noise

In the quark-gluon plasma (QGP), it is well known that the evolution of quarkonia is affected by the screening of the interaction between the quark and the anti-quark. In addition, exchange of energy and color with the surrounding medium can be included via the incorporation of noise terms in the evolution Hamiltonian. For noise correlated locally in time, these dynamics have been studied in a simple setting by Kajimoto et al.[PHYS. REV. D 97, 014003(2018)]. We extend this calculation by considering non-Abelian dynamics for a three-dimensional wavefunction. We also propose a modification of the noise correlation, allowing it to have a finite correlation in time with the motivation to include long-lived gluonic correlations. We find that in both cases the results differ significantly from solutions of rate equations.


I. INTRODUCTION
The propagation of quarkonia in the QGP is influenced by several processes. Screening in the thermal medium weakens the interaction between the Q and theQ [2]. Interaction with "on shell" thermal gluons can lead to dissociation (gluo-dissociation) [3,4]. In systems with high occupation numbers of heavy quarks (for example in heavy ion runs at the LHC) recombination [5,6] of c and c may also play an important role. All these effects play a role in the determination of the experimental observable, R AA , which is the normalized (per binary collision) ratio of the observed quarkonium yields in heavy ion (AA) collisions versus the yields in pp collisions.
The large mass of the heavy quark, M , provides a natural starting point for the analysis of these effects. There is a clear separation of energy scales between the mass M of the heavy quark (∼ 1.5GeV for c and 4.5GeV for b) and the scales Λ QCD and the temperature T 500MeV.
In contrast with open heavy flavors, quarkonia are nonrelativistic bound states and have additional scales: the inverse of the size, 1/r, and the binding energy E b . If the strong coupling α = g 2 /(4π) at the scale 1/r is sufficiently smaller than 1, then the bound states are Coulombic and these additional scales can be written in terms of the velocity v ∼ α: In this case the hierarchy of scales can be written as M M v M v 2 [7]. Even with optimistic estimates of α, the approximation α 1 is not expected to be quantitatively reliable for most quarkonium states except for perhaps the lowest bb bound state. (One way to see it is that matching the observed quarkonium spectra requires a long distance piece in the QQ potential in addition to the Coulombic piece [8,9]). It is assumed more generally that a nonrelativistic treatment of quarkonia is still valid with the hierarchy M 1/r E b . This hierarchy in scales allows for application of an effective field theory (EFT) treatment of the system which is valid at the lowest energy scale E b . At the lowest order in rE b , the EFT consists of non-relativistic quarks bound by a potential [10] (see Ref. [11] for a comprehensive review). At higher order the theory features interactions mediated by gluons of wavelength 1/E b . Effects of higher order terms are suppressed by positive powers of rE b , where factors of r can be seen as arising from a long wavelength expansion of the fields. This framework is called pNRQCD.
At T = 0 the potential can be calculated using nonperturbative techniques ( [12]). At finite T , the coupling between Q (andQ) and the gluons in the thermal medium at the energy scale T , and the coupling between the medium gluons at that scale, also play a role. It is typically assumed that 1/r T, Λ QCD but the relative hierarchy between E b , T , and Λ QCD is unclear. A finite temperature version of pNRQCD [13] has been developed to analyze this system It is well known that the QGP medium formed in heavy-ion collisions such as RHIC and LHC is best described as a strongly-coupled medium. Therefore, the ultimate goal should be to use EFT methods to write observables in terms of quantities which can be calculated on lattice. As a concrete example, the singlet potential has been computed on the lattice [14][15][16].
However, non-perturbative calculations of some relevant dynamical processes is still challenging and weakcoupling calculations are still useful. An important result in weak-coupling was obtained in Ref. [17] which showed that the Wilson loop of heavy quarks which is related to the potential between quark-antiquark pair, is complex at finite T . Furthermore, it was shown [13] that pNRQCD naturally incorporates the process known as gluo-dissociation ( [3,4]) as its dynamical degree of freedom include low energy gluonic degrees of freedom (and other light degrees of freedom if any) in addition to the wavefunctions of QQ pair.
Such weak-coupling calculations have given insight into the problem and results from these calculations can be used to obtain semi-quantitative estimates for experimental observables of interest: for example R AA in heavy ion collisions.
In the remaining part of this introduction we will review aspects of the theory particularly relevant for our work to set up our calculation.

A. Theory Overview
In Ref. [17] the QQ system was analyzed in weakcoupling in the regime where the relevant energy scales satisfy the hierarchy E b 1/r T . With these assumptions it was proved that at late times, the time evolution equation for a thermal averaged correlator for a static QQ pair, Ψ QQ ( r, t)Ψ QQ ( 0, 0) satisfies a Schrödinger like equation. The evolution kernel has an imaginary piece with the formal structure of an imaginary potential which arises due to the Landau damping of the gluons exchanged between the Q and theQ due to thermal gluons.
From the complex potential one can calculate the thermal width of quarkonia in the medium. Interpreting the inverse width for a quarkonium state as its decay rate one can solve the rate equation to find the fraction of quarkonia that survive in the medium during its evolution. Thus one has a theoretical calculation for R AA for various quarkonium states [17,[28][29][30][31][32][33].
In Refs. [13,[55][56][57] the calculation was extended by considering different hierarchies of the energy scales (between 1/r, E b , T , Λ QCD ), and additional processes like gluo-dissociation, within the weak coupling approximation using pNRQCD. Boltzmann equations in weak coupling have been written down and solved in Refs. [43,58,59]. In Refs. [60][61][62], a Lindblad equation was derived and used to obtain a Boltzmann transport equation and compute R AA .
However most calculations of QQ described above ignore the coherence of the quarkonium wavefunction on the time scale of the medium evolution. Therefore one requires a formalism which tracks the full quantum evolution of the QQ state.
The correct way to dynamically interpret the results obtained in [17] is to look at the evolution of the QQ density matrix by treating the QQ pair as an open quantum system [63,64]. The complex potential corresponds to the decoherence of a QQ state. In addition, another process -dissipation (which is required for heavy-quarks equilibration but is expected to be small for tightly bound quarkonia [63,64], however see Ref. [65]) can also be naturally derived in this formalism [64]. This approach to quarkonium dynamics was introduced in various physical regimes in Refs. [66][67][68]. It was developed in the weak coupling regime in Refs. [1,63,64,66], in the pNRQCD framework in Refs. [69,70], and more recently in Ref. [71].

B. Summary
In this paper we follow the formalism developed in Refs. [63,64]. In the weak coupling regime the authors derived equations for the evolution of the density matrix for the QQ system in contact with a thermal medium. It undergoes decoherence, which refers to processes where interactions with the environment convert a pure quantum state of the system to a mixed state. In this context it refers to scatterings with the medium gluons. If the typical energy scale of the system (here E b , which is inverse of the system time scale) is much smaller than the environment relaxation rate, then the system evolution during a typical interaction can be taken to be slow. Formally taking the system frequency to be much smaller than gT , in Refs. [17,64] a Markovian master equation in Lindblad form was derived. Then, the evolution is only controlled by two parameters -the temperature T and value of the strong coupling g. These evolution equations can be naturally solved by introducing appropriate noise fields, solving the resulting stochastic Schrödinger equations, and taking the ensemble average [72]. Refs. [17,64] derived the corresponding stochastic Schrödinger equation with a noise term which is correlated locally in time.
In Ref. [1], the authors solved a simplified version of these equations for one dimensional wavefunctions and ignoring the color structure. We expand their implementation into a more general setup with a simplification which we argue from the viewpoint of pNRQCD. We implement stochastic Schrödinger equations which keep track of the color, angular-momentum and radial wavefunction in position space for the quarkonia pair. This is the main technical advance presented in our paper.
In Sec. III, we propose a modification to the stochastic Schrödinger equation which can incorporate finitefrequency processes. We argue that the process of absorption or emission can be described if the noise field is allowed to be correlated in time with a finite time scale. This makes the system evolution non-Markovian due to memory-effects in the bath degrees of freedom. The modification can be checked by comparing the results at early time with classical decay approach.
The brief outline of the paper is as follows. In Sec. II, we extend the calculation of [1] to a realistic three dimensional case while keeping the complete color structure of the QQ pair. We also make an expansion in small r for the noise fields. A small r expansion is justified as long as r m D 1. In Sec. III we extend the stochastic equation used in Sec. II, by allowing the noise fields to be correlated in time. This allows us to perform a quantum calculation of gluo-dissociation. Our main results for the above two cases are presented in Sec. IV. We also make a comparison with simple rate-equation like approaches which has been traditionally used in phenomenological approaches.
Finally, in the appendix we provide a comparison between a r expanded and a calculation without making the expansion ("un-expanded") for a simple one-dimensional colorless system, for which results were available from [1].

II. DECOHERENCE IN SMALL r LIMIT
In this section, we briefly review the evolution equations for quarkonia in the QGP [63,64] and simplify them using the approximation r 1/T .
A. Master equation for the quarkonium density matrix ρ QQ (t) The QQ "system" continuously exchanges energy with the thermal "environment". The density matrix (ρ) of the QQ is obtained by tracing out the environmental degrees of freedom. In general the process of tracing out the environmental degrees is complicated. However tractable evolution equations for the system can be obtained under some simplifying equations. The starting point of our calculation is the evolution equation for the QQ density matrix (Eq. 1) derived in Refs. [63,64] using the following approximations.
1. All interactions are governed by a single coupling constant g and it was assumed that g 1 The evolution equation was derived keeping terms up to O(g 2 ) 2. It was assumed that E b gT . Physically this corresponds to assuming that the thermal gluons relax [on a rough time scale ∼ 1/(gT )] on a shorter time scale than the natural time scale for the system oscillations [∼ 1/(E b )]. Then each exchange with the environment can be treated as independent and hence the density matrix evolution is Markovian: the operator governing the evolution of ρ does not depend on the history and is local in time. Given that the scales T and E b are not well separated, it is worth scrutinizing this assumption further, and we will do this in Sec. III.
3. It was assumed that M is much greater than any other scale in the system. Then the Hamiltonian for the fermionic part can be expanded in powers of 1/M [73]. Only the leading order terms in 1/M were retained.
4. Under the further assumption that M T , dissipation terms are smaller than terms leading to the decoherence of the wavefunction which is the regime that we will focus on here. Using these approximations, a master-equation for the QQ pair was derived in Lindblad form [64,74] ∂ ∂t Here r corresponds to the relative separation between the QQ in the "ket" space and s is the separation in the "bra" space. M/2 is the reduced mass of the QQ system. ρ 1, 8 = ρ 1, 8 (t, r, s) are the singlet and octet components of the QQ density matrix in position space. V ( r), V ( s) correspond to the potential between Q andQ. N c is the number of color degree of freedom and C F = (N 2 c − 1)/2N c . We consider the QQ pair at rest in the medium and hence the center-of-mass coordinates R, S do not play a role and we have suppressed the dependence on them.
D( r, s) are terms related to decoherence of the QQ state [64], The function D( r) is related to the imaginary part of gluonic self-energy. It reflects the scattering rate of offshell (ω | k|) longitudinal gluons. For r ∼ 1/m D , the most important contributions are captured by the Hardthermal-loop(HTL) approximations [75]. In this approx-imation, Here m D is the Debye mass for which we use the one-loop result m D = N c /3 + N f /6 gT . N f here is the number of light flavors. It is easy to see that D( r) approaches 0 as r increases beyond 1/m D . This master equation satisfies the necessary physical constraints of linearity, positivity and trace-preservation. Techniques of quantum-state diffusion methods [72] can then be applied to numerically simulate the evolution of such a master equation. For example, Eq. 1 can be simulated using the stochastic evolution in the following manner [64].
One starts from the pure state at the initial time t 0 (although mixed states can easily be used [70]). One introduces noise fields θ a (t, r) which are picked from an ensemble which is specified by the expectation values, where .. means taking the stochastic average over the noise fields.
For each member of the ensemble θ a (t, r), ψ is evolved using the Schrödinger equation, and the density matrix can be obtained by taking a stochastic average of the outer product ρ(t, r, s) = |ψ(t, r) ψ(t, s)| .
Master equations of similar form have also been solved in [62,69,70,76] with different implementations. A simplified version of Eq. 1 was simulated in [1], where the system was assumed to be one-dimensional and the colorstructure of QQ pair was neglected (Abelian dynamics).
To incorporate these effects, we first simplify the stochastic evolution equation (Eq. 1) (and therefore its corresponding master equation) by expanding the decoherence terms in small r, s. This approximation is motivated by the hierarchy between the inverse size of the states and the temperature 1/r T . This allows us to extend the calculation to a three dimensional system while keeping all the color structure of QQ pair intact without a high computational cost. The calculation is three dimensional in the sense that we allow for transitions between different angularmomentum states (l = 0, 1). Transitions which change the angular-momentum by two units or more are suppressed by O(r 2 T 2 ). (See Ref. [69] for a similar analysis.) To check the accuracy of r approximations, we performed a similar expansion for Abelian dynamics in one dimension, for which results are known from Ref. [1]. The comparison is presented in Appendix A. Without expanding in r we were able to match the results of Ref. [1], thereby testing our implementation. Then we analyze conditions on the wavefunctions for which the r expansion is accurate. The main conclusion from the analysis is that this is a good approximation for the lowest two bound states of Bottomonia, and we focus on these states in the three dimensional calculation.
B. Small r expansion and the momentum diffusion coefficient κ We start with the stochastic evolution equation for a quark-antiquark QQ pair in its rest frame (5). The decoherence terms, in the density matrix for QQ pair, are all expanded in small r expansion.
This simplifies the calculation in two ways. First, the noise field θ a ( r, t) correlated in space, is replaced by just two different noises which are r independent and only depend on time. We only need the noise field at the center-of-mass coordinate R, and its first derivative at R. This makes generation of the stochastic noise much cheaper computationally. Second, the r expansion allows one to compute transitions between different angular momentum states, thus facilitating a three-dimensional calculation.
The r expanded stochastic evolution operator up to O( r 2 ) for a l = 0 initial state is, where the noise field θ( r, t) was defined in Eq. 4.
are operators in the color-space of QQ pair. The subscript i refers to the spatial index, and we refer ∇θ at r = 0 as θ for notational convenience.) The noises appearing in Eq. 7 can be generated as random-fluctuations correlated locally in time as, Note that the factors of g have been absorbed in the definition of the correlation function (Eq. 4). The above Hamiltonian evolution is written for a three dimensional system. Since, V ( r) is rotationally invariant, we can separate the radial part of the three dimensional wavefunction from its angular part. The wavefunction in position space can be written as where ψ(r) is the radial wavefunction and ϑ is the wavefunction in angular momentum space, with β being the polar angle and φ azimuthal angle. We also define the normalized color states for QQ octet and singlet wavefunction as, The indices l, k denotes the color states of a single quark or antiquark. T F = 1/2 is the index for SU (3) Lie group, for the fundamental representation.
Finally, we project the evolution operator in the Eq. 7 into the color and angular momentum space of QQ pair, This Hamiltonian acts on the wavefunction given in the form Here, ψ S (r, t) and ψ O a (r, t) denote radial wavefunctions for QQ pair in singlet and octet states respectively and the index a runs from 1 to (N 2 c −1) for different coloroctet states. l denotes the angular momentum states, which take the values l = 0, 1. The Hamiltonians for the singlet and octet states are One can check that the color factors between singlet and octet states are same as those obtained in pNRQCD [10,11]. Under the approximations considered, the correlation functions (Eqs. 8) are the most important quantities which control the suppression pattern. In Ref. [1], the correlation function (Eq. 4) was approximated by a gaussian function with a width l corr ∼ gT . Here we simply use the HTL form (Eq. 3). (In Appendix A we use the gaussian form since we wanted to compare with Ref. [1].) − ∇ 2 D(0) at one loop HTL is divergent. This problem is well known in the perturbative calculations of momentum-diffusion coefficients for a heavy-quark [78]. This is not physical as the problem arises from the use of the HTL form for D( r) for very short distances where it is not valid. For the momentum diffusion coefficient, the contribution from the scales above gT is important. The resulting ultraviolet divergence in the softmomentum region which is regulated by a cutoff of the order of gT is cancelled by the infrared divergence coming from the upper momentum sector k ∼ gT − T (see [13,78,79]). Since the constant κ is closely related to the physical observable such as the flow patterns of heavy quarks inside the QGP medium, it has been investigated extensively. In our calculation, we use the weak coupling result for κ at the leading-order (LO) including the UV contributions [77,80], given in Eq. 15, The value of κ has also been calculated on lattice [81,82] for a pure SU (3) theory and it was seen to be larger compared to its LO estimates from perturbative calculations,1.8 κ Lattice /T 3 3.4. Including light quarks in the calculation might modify this value further.
Intuitive understanding of the noise field can be gleaned by looking at the non-perturbative expression for the diffusion constant in terms of the correlation function of color-electric fields at different times [80], where E a i (t) is the color electric field and W (t; 0) is the gauge link in the fundamental representation.
The small r expansion of D( r) gives us exactly the same quantity which one uses in these sorts of perturbative calculation (see [64,77,80,81]).
Comparing Eqs. 16, 14 with with Eq. 8 we see that ∇θ a can simply be interpreted as E a in the temporal gauge except for the factors of g which has been absorbed in the definition of noise correlations.
It was shown in Ref. [13] that the same structure of the electric-field correlator appears in the calculation of the gluo-dissociation rate, where the gauge link connecting the two fields is in the adjoint space. In temporal gauge the two different correlator defined perturbatively becomes same. Therefore we have also plotted the relevant correlator (see Eq. 20 and denoted as κ GD on the plot) on the same plot (Fig. 1). Therefore one can argue that noise-fields here can be thought of as the electric field present in the pNRQCD lagrangian [10,11]. We can extend the definition of the θ a (t) correlator from being uncorrelated in time to have a finite correlation in time to include on-shell processes. This modification ensures the gluonic emission and absorption processes are include in our calculation. This we do next.

III. GLUO-DISSOCIATION
In this section, we describe our implementation of the quantum calculation of the process called gluodissociation [3,4,13,56,83] in literature. At finite temperature, a singlet bound state can absorb a gluon from the medium and jump to one of the excited state. This process changes the color state of the quarkonia to a color-octet state. In perturbation theory, the short distance potential for an octet state is repulsive and thus it is typically assumed that this transition destroys the bound state.
The decay rate from this process assuming T 1/r was first calculated in Refs. [3,4]. More recently, the decay rate was computed for T E b in Ref. [13] and 1/N c corrections for T 1/r were computed in Ref. [56].
The process of gluo-dissociation is naturally described in pNRQCD [10,11]. This EFT is valid at energies (this could be an energy scale like T , m D or E b depending on the hierarchies between these scales) 1/r. The degrees of freedom in the theory are light degrees of freedoms like gluons and light quarks, and the singlet and octet wavefunctions of the QQ. (For a detailed study of pNRQCD at finite T see [13].) Starting with the lagrangian, where Ψ S ( r, t) = Ψ( r, t) ⊗ |S and Ψ O ( r, t) = Ψ a ( r, t) ⊗ |O a are the three dimensional wavefunctions of QQ pair in color-singlet and octet states, where color states |S and |O a has been defined in Eq. 10. ∂ 0 denotes the time derivative and D 0 = ∂ 0 − igA a 0 t a is the covariant derivative acting on the octet states. The Hamiltonian for singlet and octet states are The interaction vertices are same as in Eq. 11. F a µν ( R, t) is the field strength tensor for the long wavelength gluons and q( R, t) represents the light quarks. Because of the multipole expansion of pNRQCD, all the light degrees of freedom are function of center-of-mass coordinate ( R) only. V A and V B are the coefficients of dipole-interactions. At leading order in α they are 1.
The wavefunctions Ψ S , Ψ O in Eq. 17 are different from the wavefunctions ψ S , ψ O in Eq. 12 as they are for a three dimensional system right now. One can project out the above lagrangian in the color and angular momentum space of QQ to get back to an equation similar to Eq. 11. The singlet to octet transition rate in a thermal medium at a uniform, time independent, temperature T to first order in perturbation theory is The integration is over the set of continuum of octet states. We will focus on l = 0 singlet initial states and therefore the final octet states have l = 1. ∆E = E(ψ O ) − E(ψ S ).G(ω, T ) is given by the thermal expec- Where φ ab (t) is the gauge link connecting the two electric-fields in adjoint representation [13]. Looking at Eq. 20, we see that it has the same structure as that of the correlator in Eq. 16 with two important differences. First, the gauge link is adjoint in Eq. 20 and fundamental in Eq. 16. Second, Eq. 20 has an additional factor of e −iωt corresponding to the fact that during the gluo-dissociation process, the QQ state absorbs energy ω from the gluon.
The absence of e −iωt in Eq. 16 can be traced to the hierarchy between the energy scales assumed in the derivation of Eq. 1. E b gT implies that the relaxation time scale for the thermal gluons is much shorter than the system time scales. Therefore, the electric field correlator can be taken to be local in time on long time scales. Physically it corresponds to the assumption that there are no long-lived (compared to 1/E b ) gluonic degrees of freedom in the medium.
On relaxing this assumption, ω can no longer be taken to be zero in Eq. 20, and the electric field correlator has a finite correlation in time. In the calculation of the decay rate this does not cause any technical complication (Eq. 19). However, in a quantum calculation which follows the density matrix evolution of the QQ, the steps  show that the evolution of survival probability is non trivial. Comparing this with Fig. 2, we see that their behavior is different for two different medium effects.
involved in deriving a Markovian evolution in the form Eq. 1 can no longer be followed.
As an illustrative example, consider a regime when the relaxation rate of the thermal gluons is small compared to T , and let T be comparable to E b . Concretely, one scenario where this can be realized in the weak coupling regime when the gluonic screening mass, ∼ gT , and the relaxation rate, ∼ g 2 T , [75] are both much smaller than T . Then, at leading order in g, the electric field correlator can be written as, where a, b are color indices and i, j are spatial indices. In this case the thermal gluons can not be integrated out from the influence functional to obtain an interaction term which is local in time. To make progress on the quantum implementation in presence of a correlated electric field, we start from the stochastic Schrödinger Eq. 11. Following the interpretation in Sec. II of the noise field θ a (t) as g E a (t), the correlation function of stochastic noise is given as θ a (t) = 0 .
The correlator θ i a (t) θ j a (t) is given by Eq. 21. The density matrix at any given time can be obtained by taking the noise average Eq. 6. The evolution equation for the density matrix thus obtained can not be written in a Markovian form as the correlations between the noise terms are not local in time.
A more rigorous approach to obtaining a time evolution equation for the density matrix would involve deriving the influence functional without making an expansion in ω, and using the full gluonic propagator. Here we have used the lowest order form for the electric-field correlator (Eq. 21). At one-loop the spectral function of gluons changes drastically (see fig. 1) for the particles with momenta less than T . The finite thermal mass and decay width is important and can not be ignored. These corrections will change the spectral density and also the analysis of non-Markovian regime. Finally spontaneous emission processes need to be included. We leave these considerations for future work. In spirit, our calculation is similar to what was done in Ref. [66], before it was made theoretically concrete in subsequent works Refs. [63,64].
At this point we would like to make a comment about an alternative approach to deriving the quantum evolution equations for quarkonia in the QGP. Open quantum treatment of quarkonia starting from the pNRQCD lagrangian has been performed in [69,70]. The authors derived a general evolution equation for the QQ density matrix including gluo-dissociation processes. Furthermore, in two different physical regimes, they were able to write the density matrix equations in Lindblad from. The first case was the strong coupling regime, in which case the static limit of the electric field correlator was considered. The second case was the weak coupling limit, g 1, where the leading order form for the electric-field correlator was used just as in Eq. 21. The evolution equations were simplified using the hierarchy E b T . In this case an expansion in (V o − V s )/T is possible, and the selfenergy correction and the gluo-dissociation rate can be simplified. We do not make this assumption. As a result we can not write a simple equation for the density matrix evolution and prove its validity in the weak coupling regime.
However, we believe that this is a good first step towards incorporating non-Markovian effects in the QQ evolution equations in the presence of long-lived gluonic degrees of freedom. This can also be confirmed by looking at the classical decay picture with the quantum one at early time. We confirm in Sec. IV that for small time, both results follow each other and the two approach start diverging at late times (the details depend on the initial states chosen, see IV).
In a medium evolving with time we can modify the generation of noise to incorporate the dependence of T on time. Assuming that temperature change is slow enough, we can approximate the physical picture as follows.
1. The entire evolution of the QQ pair is divided into time blocks. During each block we take the T to be constant and equal to the mean temperature in the block. The division has to be done while keeping in mind that the time blocks we choose are large enough to include the finite correlation for the dominant gluons (which here are of the order ∼ T ).
3. Then for each block, the noise is generated using the equilibrium correlation function. They are stitched together using a linear interpolation functions, λ β (t) in Eq. 22, which are normalized to one.
We have used 3 time-blocks for the results presented in Sec. IV. We have also checked that using 5 blocks gives the same results. The stitching was done as given below.
Here, θ a i,β (t) denotes the i'th component of noise θ a (t) generated in β'th block and G(t, T ) is given by Eq. 21.
Having described these two different decay mechanisms and our implementation, we proceed to the next section where we present our main results.

IV. RESULTS
The main results of the paper are presented here. We calculate the survival probability of the vacuum states by doing a three dimensional quantum evolution of the density matrix using the stochastic Schrödinger equation defined in the Eq. 11, for two different physical cases (decoherence and gluo-dissociation). The survival probability is defined as where |ψ 0 is the vacuum wavefunction for 1S or 2S states. The evolved wavefunction is The quantity P (t) is related to the observed suppression number R AA of quarkonium states at RHIC and LHC. Typically, phenomenological calculations of R AA in the literature (see [18,[31][32][33] use a classical rateequation approach. We call these approaches "classical", since a full quantum evolution of the density matrix is not done. Quarkonia has a finite decay width inside a thermal medium due to different physical processes, such as inelastic scattering with thermal particles, gluodissociation etc. The width can be calculated in perturbation theory at desired order. Suppose there were initially N ψ numbers of quarkonia in some state, labelled as ψ here. The number of surviving quarkonia after a finite time t in the state ψ, in the classical approach is given by The value of Γ ψ (t) also depends on the choice of wavefunction one uses to calculate the width. For example in [31][32][33], instantaneous value of the width was used to calculate N ψ (t) by solving the three-dimensional Schrödinger equation at each time step using a complex potential. Its time dependence comes from the fact that dissociation rate depends on T . The quantity P (t) defined in Eq. 23 is equivalent to N ψ (t) defined above in the sense that starting from N QQ ψ pair in state ψ at time t 0 , the number of surviving pair after time t is P (t) × N . Therefore, from here on we use P (t) to denote both the quantum and classical survival probability.
In our calculation, we use the vacuum wavefunction at each time step to calculate the value of the width. The expression for decay width for two different cases (subscript 'dc' for decoherence and 'gd' for gluo-dissociation) is (26) where the quantities κ andG(ω, t) have been defined The production cross section of quarkonium states in heavy nuclei relative to proton-proton collisions is still an active area of research (see [84] and references therein). Different initial states have been used to calculate the survival probability of quarkonia in medium. For example in [70] initial states were chosen to be a delta function in position space in l = 0 state. In [1,65] initial states were chosen to be eigenstates of the vacuum Cornell potential. To investigate the effects of size and shape of initial wavefunction, we choose first two lowest lying eigenstates of the Coulomb and Cornell potential for our calculation.
For the eigenstates of the Coulomb potential, we used the following parameters These value of the α is determined by the self consistency equation where a 0 is the radius of the ground state of bottomonium. For the initial states of the Cornell potential we used the following parameters where r 0 = 1.2fm is the string breaking parameter (threshold for heavy-light meson production) as determined in [85]. These parameters were taken from [86] (although we use M = 4.8GeV whereas the mass used in [87] was M 4.6GeV).
To implement the evolution of the temperature as the QGP medium cools down, we use the expression for a Bjorken expanding medium. These parameters were taken from [88].
which were also used in [70].
The survival probability as a function of time t for 1S and 2S states for the decoherence case has been presented in Fig. 2. The same results for the case of gluodissociation has been presented in Fig. 3. 1. The average radii r avg = r for 1S states of Coulomb and Cornell potentials are r avg (Coulomb, 1S) = 0.26fm, r avg (Cornell, 1S) = 0.16fm. We find that for both the casesdecoherence and gluo-dissociation ( Fig. 2  2. For 2S states the Cornell initial wavefunctions is much narrower than Coulomb one: r avg (Coulomb, 2S) = 0.95fm, r avg (Cornell, 2S) = 0.29fm. The difference in P (t) between Coulomb 2S and Cornell 2S originates from the huge difference in their average radii. However, the evolution pattern is not very intuitive. Figs. 2 and 3, we see that despite being a much narrower state, P (t) for Cornell 2S is different but of the same order of magnitude as Coulomb 2S. In the evolution, we have taken the potential to be screened Coulomb, which is closer in form to the Coulomb potential. Just the difference in the evolution potential to the potential used to calculate the eigenstate, leads to a rapid change in the wavefunction for the Cornell 2S state. (This effect is very prominent in particular for decoherence.) On the other hand, we expect decoherence and gluo-dissociation to be more important for the initially wider Coulomb state. The competition between these is subtle. Such large effects for 2S arise from the QQ wavefunctions becoming broad with time very quickly, and suggest that other effects that we have ignored here (in particular dissipation) could play an important role and need to be studied further. For the eignestates of the Cornell potential it would also be natural to evolve using a non-perturbative potential obtained from the lattice. The calculation of non-perturbative forms for both the real and imaginary parts for the potential at finite temperature is an active area of research [14][15][16], and we leave this exercise for future.

From the
4. Comparing the P (t) for 1S states, for two different cases, we find that gluo-dissociation has a much stronger effect on quarkonium decay than decoherence.

5.
A direct comparison of our results with the "strong coupling" results of Refs. [69,70] is not possible as we do not work in that regime, but the closest comparison that we can consider is between our gluo-dissociation results and the "weak coupling" results of Refs. [69,70]. The main difference is that 2S shows substantially larger suppression than 1S in our calculation which is not seen in Ref. [69,70]. This might be because of the following reasons (g and parameters in the Bjorken expansion are taken to be the same) (c) It could also be due to a difference in the choice of the initial state.
We have presented our comparison of survival probability P (t) between the classical and quantum approach in Figs. 4 and 5 for decoherence and gluo-dissociation respectively. We only present our results for 1S states as for 2S states, as discussed above, additional effects might play an important role.
We note from the Figs.4 and 5 that the two approaches give different results for both decoherence and gluodissociation. For both Cornell and Coulomb 1S, the quantum decay probability is substantially larger than its classical counterpart. One can understand this as follows. The wavefunction gets wider as it evolves in time, therefore at a later time the decay rate will be much higher than what it was at early times in the quantum scenario. If the medium evolution is quasi-static (the time scale over which width becomes constant is very small compared to other time scales) and higher order contributions are not important, one would expect that both quantum and classical approaches should give similar results.
The numerical details of our calculations are as follows. We took a lattice of spatial extent x ∈ [−2.56, 2.56]fm with 200 lattice points. The states were evolved under the stochastic Hamiltonian given in Eq. 11 from t 0 to t max = 5fm with time-steps of size dt = t max /400. We used 500 instances to perform the stochastic averaging.

V. CONCLUSIONS AND OUTLOOK
In this paper we studied the quantum evolution of the density matrix of QQ inside a hot quark-gluon plasma. We have used the techniques of stochastic Hamiltonian evolution to simulate the density matrix evolution [72]. The evolution is unitary and therefore the number of QQ pairs is conserved. The main technical advancement in this paper over the work done in [1,65] is that the wavefunction is three dimensional and the evolution of the complete color structure of the QQ pair was done.
Our starting point was the recoilless master equation derived for a QQ pair inside a weakly coupled medium in [64]. The master equation describes the process of decoherence in a regime where effects of dissipation can be ignored. (See Ref. [65] for a more quantitative estimate of dissipation effects.) We argued that for a medium at temperature T which satisfies the relation rm D 1, for a quarkonium state of size r, an expansion of the stochastic Hamiltonian in r is justified. In Sec. II B, we derived a small r expanded version of the stochastic evolution operator derived in [64] in recoilless limit. We tested this expansion in Appendix A for a one dimensional colorless system for which results are available from [1]. For the un-expanded case our results matched the results of [1]. We checked that the expansion in r gives accurate results for the lowest 3 states of the Cornell potential.
This expansion allows one to solve the equations for three-dimensional wavefunctions by including transitions between different angular momentum states. r expansion also makes the generation of noise much cheaper computationally. Finally, in the r expansion, we can relate the correlator of stochastic noise to the momentum-diffusion coefficient κ which can be expressed as correlator of color electric field [63,80]. A similar correlator was derived for the process gluo-dissociation in [13].
Since the hierarchy between E b , T and m D is not very clear for the realizable temperatures at RHIC and LHC, a Markovian evolution of the density matrix may not be well justified. Therefore in Sec. IV we proposed a modification of the stochastic noise correlator from Eqs. 8 to Eq. 21 to include on-shell gluons in our calculation. The main idea was to implement a stochastic evolution equation which gives us the same decay rate as when calculated in pNRQCD for gluo-dissociation at leading order in g. A quasi-static medium evolution was assumed to perform the calculation for the Bjorken expanding medium.
Finally, in the Sec. IV we made a comparison of the survival probability P (t) when calculated in a classical rate-equation approach versus in a quantum approach. Typically, most phenomenological calculations of R AA for quarkonium states in literature have implemented a rate-equation approach. We call it "classical", since a quantum evolution of the density matrix was not done. We found that the two approaches do not always produce the same results for the survival probability P (t). The difference depends on the initial states chosen (we considered eigenstates of the Coulomb potential and the Cornell potential as examples) and also the form of potential one uses to evolve. Our main results were presented in the Sec. IV where we compared the survival probability between classical and quantum approaches separately for decoherence and gluo-dissociation.
For the 1S state we found that for the choice of parameters given in Sec. IV, gluo-dissociation gives a substantially larger suppression compared to decoherence. For both the cases -decoherence and gluo-dissociation -we found that for 1S states, the survival probability is very similar for the two choices for the initial wavefunction. However, it should be noted that the effects of potential change on Coulomb and Cornell states are quite different. For both cases, we found that the 2S states are highly suppressed relative to 1S states.
For the 1S state we also found that the quantum calculation shows larger suppression than the classical calculation. Finally, the dependence of classical survival probability on the initial wavefunction is much stronger compared to its quantum counterpart.
Our work can be extended in several directions. Within our framework, dissipative effects [64] can be included. Their effects have been in studied in a recent work [65] for one dimensional Abelian dynamics, and our implementation can extend it to three dimensional wavefunctions with color dynamics.
We would also like to put our formalism on a stronger theoretical footing. A simple conceptual advance would involve going to higher order in in g in Eq. 21. This introduces screening and a width for the gluons, thereby relaxing a severe approximation in our calculation of gluodissociation (Sec. III). Eventually it would be very useful to derive quantum evolution equations for QQ from first principles, only assuming M 1/r E b . Refs. [69,70] have derived Lindblad equations by making choices about the hierarchy between the scales E b , T and gT . But since these scales are not well separated, and it would be interesting if evolution equations can be derived without making these approximations.
Eventual connection with phenomenology would require effort in other directions. The initial production of quarkonia from QQ and the dynamics of quarkonia in the pre-thermalized medium, t 1fm have to be investigated to get a better understanding of the initial wavefunctions and remove an important systematic uncertainty. To make a connection to the observed experimental R AA for quarkonium states, in addition, one needs to take as a background thermal system, a realistic three dimensional hydrodynamic simulation.
We hope to make progress on these directions in future.

VI. ACKNOWLEDGEMENTS
We would like to acknowledge discussions with Dibyendu Bala and Sourendu Gupta. A.T. would like to thank Shiori Kajimoto for valuable discussions. We would, in particular, like to acknowledge many conversations with Saumen Datta. RS would like to thank the hospitality of INT during Program INT-19-1a and IIT Guwahati during WHEPP 2019, where part of the work was done. We acknowledge support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0200. In this section we numerically investigate a few examples to test under what conditions the expansion of the stochastic Hamiltonian in r is a good approximation.
We reproduce results Ref. [1] obtained without making an expansion in r, checking our implementation. We then verify that for wavefunctions smaller than the noise correlation length, the expansion works quite well.
Ref. [1] investigated two kinds of systems. First they considered states propagating in a time independent thermal system at temperature T . The initial states for this system were taken to be eigenstates of the screened Coulomb interaction (the real part of the singlet potential). Second, they considered states propagating in a Bjorken expanding medium. The initial state for this system was considered to be the eigenstate of a vacuum potential of a Cornell form.
For this section we consider parameters (coupling constants, temperatures, and parameters in the Bjorken expansion) which are the same as those considered in Ref. [1] for easy comparison. These are different from those used to obtain our final results in Sec. IV.
(A1) and D( r) was modelled as a gaussian,

Time independent background
The strength γ was taken to be which includes the factor C F which carries the imprint of color in the Abelian dynamics. In weak coupling l corr ∼ gT (Eq. 3). However, since the hierarchy between T and gT (g ≈ 1.7 for g 2 C F 4π = 0.3) is unclear, Ref.
[1] considered a range of values of l corr varying from 0.04fm to 0.96fm. Here we compare our results for two values within this range, The stochastic fields are correlated over length l corr . Therefore, if the hierarchy l corr r is satisfied, one can expand the stochastic fields present in the above equations around r = 0. Then up to leading order in r/l corr , we expect the system to interact with the environment as a dipole. The stochastic Schrödinger equation in this approximation can be written as: In this section we would like to examine how reliable such an expansion is in practice. EFTs like pNRQCD and its various versions are based on these assumptions and this exercise provides with a simple check.
Since we are in a one-dimensional system, let us make our naming scheme clear. We label the lowest lying ground state for the given Hamiltonian by n 1 = 0. The first and second excited states are respectively labelled as n 1 = 1, 2. The averaged squared radius is defined by r rms = r 2 . For the Debye screened potential used in [1], with parameters given in Table I the value of r rms for the first two bound states is r rms (n 1 = 0) = 0.11fm (A7) r rms (n 1 = 1) = 0.66fm.
(A8) A convenient dimensionless quantity to characterize the separation of scales is ξ = r rms /l corr which we want to be much smaller than 1. In this section we take T = 0.4GeV to be constant in time.
Taking l corr = 1/T gives the value l corr = 0.5fm. For g corresponding to α eff in Table I, the second value of l corr is 1/gT = 0.3fm.
Taking l corr = 1/T , we find that ξ is indeed smaller than 1 for the ground state. For the n 1 = 1 state the ratio is not small and one expect that r expansion will FIG. 6: (color online) Survival probability for expanded vs un-expanded case for time-independent medium (l corr ∼ 0.3fm). For l corr = 1/gT , both n 1 = 0 (red solid and blue dash dotted lines) and n 1 = 1 (black solid and pink dash dotted lines) start differing as soon as t = 1.5fm by more than 5%.
We see that for l corr = 1/T , the expanded vs unexpanded results are very close for n 1 = 0 state. For l corr = 1/T at final time the expanded result is 3.8% smaller than the un-expanded case. For l corr = 1/gT at final time the expanded result is 16% smaller than the un-expanded case.
For n 1 = 1 state we see that r expansion breakdowns and therefore is not reliable. For l corr = 1/T at final time the expanded result is 40% smaller than the un-expanded case. For l corr = 1/gT at final time the expanded result is 56% smaller than the un-expanded case. We conclude that for the realistic system with initial states of similar size, r expansion is not a reliable tool.
To understand the implications for the three dimensional calculation, we note that eigenstates in the one dimensional problem and the three dimensional problem are related. For a rotationally invariant potential, the radial part of a three dimensional Schrödinger equation is equivalent to a one dimensional Schrödinger equation with the additional constraint that the wavefunction is 0 at the origin for l = 0 states, l being the quantum number for orbital angular momentum. The n 1 = 0 state is finite at the origin and does not correspond to any three dimensional state. The n 1 = 1 state corresponds to the n = 0 three dimensional, l = 0 state. Therefore our results show that for an eigenstate of the screened Coulomb potential in three dimensions with value of α similar to table I , the r expansion is invalid and a full three dimensional simulation is necessary. However, it is reasonable to argue that in a rapidly evolving plasma, it is not well motivated to use the eigenstate of the screened Coulombic state as the initial state of evolution anyway. It is more appropriate to start the evolution with a narrow state, which is often taken in the literature to be the state in the vacuum. This is the system we analyze in the next section.
Finally, we comment on some technical aspects of our implementation. We took a lattice of spatial extent x ∈ [−2.56, 2.56]fm with 512 lattice points. The states were evolved under the stochastic Hamiltonian given in Eq. (A1 and A5) for t max = 5fm with time-steps of size dt = t max /5000. We used 1000 ensembles to perform the stochastic averaging. FIG. 8: (color online) Survival probability for a Bjorken expanding medium. For Cornell states, n 1 = 0 (red solid and blue dash dotted lines) and n 1 = 1 (pink solid and cyan dash dotted lines) expanded vs un-expanded results are very close to each other. Only for n 1 = 2 state (black solid and sky-blue dash dotted lines), we see a big difference around t 2fm between the two cases, which narrows down as the wavefunctions evolves further in time.

Non-equilibrium QGP
The Bjorken expansion is a well studied model for QGP dynamics and in this section we present the same comparison in the case of a Bjorken expanding medium. The temperature changes with time according to the relation, (A11) Although the full dynamics of the evolution of medium can be quite complicated, for a simple calculation the picture of the Bjorken expansion is a good check. The initial temperature was chosen to be T 0 = 0.4 and t 0 = 1fm to match Ref. [1]. The initial states were chosen as the first three bound states of the one dimensional Cornell potential, V Cornell (r) = σr − α eff r , σ = 0.2GeV, α eff = 0.3 (A12) and the parameters of the model are same as in Table I.
After solving for the eigenstates the Cornell potential in one dimension, we obtain the following values for r rms r rms (n 1 = 0) = 0.098fm r rms (n 1 = 1) = 0.27fm r rms (n 1 = 2) = 0.53fm (A13) In this section we only show the results for l corr = 1/T . As both the wavefunction width and the temperature change with time in this system, it is not convenient to quote a single number ξ to check whether the expansion in r will be accurate. We simply compare the results with and without the approximation to see how well it works. The result is presented in the Fig. 8. At t = 5fm the difference between the expanded and un-expanded cases is given below, P (t = 5fm)(n 1 = 0) = 0.89,P (t = 5fm)(n 1 = 0) = 0.89, P (t = 5fm)(n 1 = 1) = 0.16,P (t = 5fm)(n 1 = 1) = 0.16, P (t = 5fm)(n 1 = 1) = 0.047,P (t = 5fm)(n 1 = 1) = 0.044 (A14) For the Bjorken case we find that the r expansion is a much better approximations. The largest different between the expanded vs un-expanded case is for n 1 = 2 case, which at t = 5fm is 6%. This originates from two facts 1. The Cornell wavefunctions are much narrower compared to the Debye screened potential ones. Even for the n = 2 wavefunction, the value of ξ ∼ 1 at the earliest time.
2. As the medium cools down, the l corr grows larger. At late times, it will be harder for the medium to resolve the details of QQ. For e.g. at t = 0, the typical correlation lengths are of the order l corr = 0.5fm which grows up to l corr = 0.9fm by t = 5fm.
Translating this to three dimensions, the n 1 = 1 and the n 1 = 2 states correspond to 1S and 2S states of a realistic three dimensional system. We conclude that at least for the initial states chosen from the Cornell like potential evolving in a Bjorken expanding medium, r expansion can be used reliably.
Similarly, for the lowest three eigenstates of the three dimensional Coulomb potential with larger value of α eff , such as those used in [70], one can use the r expansion to simplify the calculation.