Bottomonium suppression in an open quantum system using the quantum trajectories method

We solve the Lindblad equation describing the Brownian motion of a Coulombic heavy quark-antiquark pair in a strongly coupled quark-gluon plasma using the highly efficient Monte Carlo wave-function method. The Lindblad equation has been derived in the framework of pNRQCD and fully accounts for the quantum and non-Abelian nature of the system. The hydrodynamics of the plasma is realistically implemented through a 3+1D dissipative hydrodynamics code. We compute the bottomonium nuclear modification factor and compare with the most recent LHC data. The computation does not rely on any free parameter, as it depends on two transport coefficients that have been evaluated independently in lattice QCD. Our final results, which include late-time feed down of excited states, agree well with the available data from LHC 5.02 TeV PbPb collisions.


Introduction
The aim of relativistic heavy-ion collisions is to recreate and study the quark-gluon plasma (QGP), a primordial state of matter that existed microseconds after the Big-Bang. In such a state quarks and gluons are not confined within hadrons and, in the limit of high temperatures, they behave almost as free particles. The study of the QGP is very challenging due to its short life time: we can only infer its properties from the way it affects particles that are detected at the end after freeze-out.
Among the hard probes of the QGP there is quarkonium suppression. The original idea was put forward by Matsui and Satz in [1] who assumed the quarkonium interaction to be screened in the hot QGP leading to a suppression in the number of quarkonia. Quarkonia would eventually be detected through their decays into muon-antimuon pairs. In this scenario, the observation of quarkonium suppression is a signal of QGP formation and, by measuring its strength, it is possible to learn information about the medium. Such suppression is quantified in the quarkonium nuclear modification factor R AA . It is this quantity, among others, that is measured in heavy ion collision experiments at LHC and at RHIC.
The screening mechanism underlying the idea of Matsui and Satz originates from the chromoelectric screening of the medium. It sets in at the momentum scale given by the Debye mass, m D . Hence, the typical screening distance for a heavy quark-antiquark pair (QQ) is of order r ∼ 1/m D or larger. Sequential screening as a function of the radius of the quarkonium state is a consequence of this mechanism. This paradigm was challenged when the quark-antiquark potential in the medium was first calculated in weak coupling in the screening regime r ∼ 1/m D [2]. The calculation showed that, in addition to the screening of the real part, the potential also posseses an imaginary part. The imaginary part also leads to quarkonium dissociation, which turns out to happen at a temperature lower than the screening one [3][4][5][6], i.e., at least in weak coupling, quarkonium has already dissociated when reaching the screening temperature. Since then, many phenomenological studies solving the Schrödinger equation with a complex potential have appeared [7][8][9][10][11][12][13][14][15][16].
Nonrelativistic Effective Field Theories (NREFTs) allow one to appropriately define the heavy quark-antiquark potential and supply a scheme for the systematic calculation of quarkonium properties. They exploit the separation of energy scales characteristic of nonrelativistic bound states. At zero temperature, the energy scales are the heavy-quark mass, m, the inverse of the Bohr radius of the bound state 1/a 0 ∼ mv, and the binding energy E ∼ mv 2 , where v 1 is the relative quark velocity in the bound state. The EFT that is obtained by integrating out degrees of freedom associated with the scale m is Non Relativistic QCD (NRQCD) [17,18] and the EFT obtained by integrating out gluons with momentum or energy scaling like the inverse of the Bohr radius is potential NRQCD (pN-RQCD) [19][20][21]. At leading order in v, the equation of motion of pNRQCD is the quantum mechanical Schrödinger equation for a nonrelativistic bound state. Differently from a pure quantum mechanical treatment of the bound state, however, pNRQCD provides an unambiguous field theoretical definition of the potential. The potential encodes contributions coming from modes with energy and momentum above the scale of the binding energy. Moreover, pNRQCD adds systematically to the leading order Schrödinger equation higher order corrections. The first one in a weak coupling regime is carried by chromoelectric dipole terms.
In the last decade, pNRQCD has been applied also to study quarkonium at finite temperature [3,4,22]. At finite temperature more scales are relevant, for instance the temperature T and, at weak coupling, the Debye mass m D . Nevertheless, at leading order in v the equation of motion of pNRQCD is still a Schrödinger equation, which describes the real time evolution of the QQ pair in the medium. The potential encodes now also thermal contributions if there are thermal modes associated with energy scales larger than mv 2 . The thermal part of the potential has a real part (well described in the weak coupling regime by the singlet free energy in Coulomb gauge [23]) and an imaginary part. The real part of the potential is screened only if m D is of the order of the inverse of the Bohr radius. If it is smaller, then the potential gets at most thermal corrections that are power like in T . It is in this regime, and not in the screened one, that dissociation happens at weak coupling due to the imaginary part of the potential being as large as the real one [2]. At one loop, the imaginary part is a consequence of two distinct phenomena: Landau damping [2][3][4], an effect that exists also in QED, and singlet-octet transitions, which are specific of QCD [3]. The Landau damping originates from the inelastic scattering of the heavy quark or antiquark with the partons in the medium [24], while the singlet to octet transition originates from the gluodissociation of quarkonium [25]. Which phenomenon dominates depends on the ratio between the scales m D and E. These findings in the EFT, mostly in the weak coupling regime, have inspired several subsequent nonperturbative calculations of the static potential at finite T . In particular, the Wilson loop at finite T has been computed on the lattice [26,27] finding hints of a large imaginary part. These calculations are rather challenging and refinements of the extraction methods are currently in development [28,29].
Quarkonium scattering in the medium, quarkonium dissociation into an unbound color octet QQ pair, and the inverse processes of QQ pair generation call for an appropriate framework to describe the quarkonium non-equilibrium evolution in the QGP: the open quantum system framework (OQS) (see [30] for a review and [31] for a seminal paper). The system is in non-equilibrium because through interaction with the environment color singlet and octet QQ states continuously transform into each other although the total number of heavy quarks is conserved. In [32][33][34], an OQS framework rooted in pNRQCD has been developed that is fully quantum, conserves the number of heavy quarks and takes into account both the color singlet and the color octet QQ degrees of freedom. In this framework, the QGP plays the role of the environment characterized by a scale πT and the quarkonium is the system characterized by the scale E. The inverse of the energy can be identified with the intrinsic time scale of the system, τ S ∼ 1/E, and the inverse of πT with the correlation time of the environment, τ E ∼ 1/(πT ). If the medium is in thermal equilibrium, or locally in thermal equilibrium, we may understand T as a temperature, otherwise it is just the inverse of the correlation time of the environment. The medium can be strongly coupled. The evolution of the system is characterized by a relaxation time τ R that is proportional to the inverse of the QQ self-energy in pNRQCD, i.e. τ R ∼ 1/[a 2 0 (πT ) 3 ]. Under the condition that the quarkonium has a small radius (i.e. that it is a Coulombic bound state) such that 1/a 0 πT, Λ QCD , and that πT E, a set of master equations governing the time evolution of the heavy QQ pairs in the medium has been derived in [32,33]. The equations express the time evolution of the density matrices of the color singlet and color octet QQ states. They account for the mass shift of the heavy QQ pair induced by the medium, the decay widths induced by the medium, the generation of QQ color singlet states from QQ color octet states interacting with the medium and the generation of QQ color octet states from QQ (color singlet or octet) states interacting with the medium. At leading order the interaction between a heavy QQ field and the medium is encoded in pNRQCD in a chromoelectric dipole interaction.
The master equations are, in general, non Markovian. They become Markovian if τ R τ E , while the condition τ S τ E qualifies the regime as a quantum Brownian motion. Under these conditions, the master equations assume a Lindblad form [35,36]. If we further consider the medium isotropic and the quarkonium at rest with respect to the medium, in the large time limit the Lindblad equation for a strongly coupled medium turns out to depend, remarkably, on only two transport coefficients: the heavy quark momentum diffusion coefficient, κ, and its dispersive counterpart γ. Both coefficients have a field theoretical definition: they are given by time integrals of gauge invariant correlators of chromoelectric fields [32][33][34]. These coefficients can be evaluated nonperturbatively within QCD at finite T regularized on the lattice. 1 Recently, a quenched lattice QCD evaluation of κ in an unprecedented large window of temperatures has been completed displaying a significant temperature dependence [37]. Similar studies are ongoing for γ. Unquenched values for both parameters have been obtained in [34]. Once the Lindblad equation has been solved and evolved up to freeze out time, the result can be used to compute observables like R AA and v 2 by projecting on the quarkonium states of interest. These observables can be compared with data from the LHC [32,33]. We are, therefore, in a position to determine much needed information about the properties of the QGP by solving evolution equations for the out of equilibrium dynamics of the heavy quarkonium in the QGP that have been derived in a controlled and systematic fashion from QCD. However, for this purpose we need to develop an efficient algorithm to solve the resulting dynamical equations and we need to couple them properly to the hydrodynamical evolution of the medium.
In this work, we revisit the solution of the evolution equations found in [32,33] with the aim of determining R AA and comparing with experimental data. We extend this previous work by relaxing several approximations that were implemented in the numerical solution of the Lindblad equation due to the high computational cost associated with it. The main goal of this paper is to present a new method for solving the Lindblad equation that substantially increases numerical efficiency through massive parallelization. In addition, the new framework allows for realistic hydrodynamical evolution of the medium.
Solving the Lindblad equation is challenging, the reason being that, for a Hilbert space of dimension N , one needs to compute the evolution of a density matrix with N 2 entries. In the case of quarkonium, the Hilbert space has infinite dimensions. For numerical studies, quarkonium is simulated using a finite size lattice. The main approximations that were used in [32,33] were the following. (i) A lattice of size 40 times the Bohr radius of the Υ(1S), a 0 , with a spacing of 0.1a 0 was used. Note that doubling the size of the lattice makes solving the Lindblad equation four times more expensive. (ii) An expansion in spherical harmonics was used and only S-wave and P -wave states were considered. (iii) For each centrality window only a simulation with the average temperature of the window was performed. Moreover, a temperature profile given by boost-invariant ideal Bjorken evolution was used for its simplicity and its analytic closed form.
In the present study, we relax these approximations by using the Monte Carlo wave- 1 It is remarkable that the OQS/pNRQCD framework allows to use input coming from a lattice QCD calculation of a quantity in thermal equilibrium to describe the out of equilibrium evolution of quarkonium in the medium.
function method [38]. Other methods exist and have been used to simplify the computation of the Lindblad equation in the context of quarkonium suppression [39,40]. However, for reasons that we will explain later, we believe that the Monte Carlo wave-function method is more suitable for quarkonium studies, especially when color degrees of freedom are taken into account. The method reduces to simulating the evolution of an ensemble of vector states that, on average, behave like the density matrix. The vector states evolve either according to an effective non-Hermitian Hamiltonian or by quantum jumps at random times. In our case, the Hamiltonian does not mix states with different color or orbital angular momentum. Therefore, the most numerically demanding task of the algorithm is to solve a Schrödinger equation in one dimension and this even when all possible orbital angular momenta are taken into account. In summary, the Monte Carlo wave-function method allows to decrease the lattice spacing and increase the volume without increasing the numerical cost as much as it would require directly solving the Lindblad equation and, additionally, without truncating the spherical harmonics as was the case in [32,33].
We believe that our method may be useful also for similar phenomenological applications. The semiclassical limit of the Lindblad equation has been studied in [41] and the relevance of correlated versus non correlated noise in [40]. In [42][43][44], using the same pNRQCD and OQS framework that we use here and a specific scale hierarchy, transport equations and, in particular, a semiclassical Boltzmann equation, have been obtained for the evolution of quarkonium in medium. For the differential reaction rate, the information about the QGP is contained in a chromoelectric gluon correlator that involves also stapleshaped Wilson lines. Similar correlators show up at T = 0 in the gluon parton distribution functions, the gluon transverse momentum dependent parton distribution functions and in the quarkonium production cross section expressed in pNRQCD [45]. Other applications of the OQS framework that do not use the EFT approach can be found in [41,[46][47][48][49].
The outline of the paper is as follows. In section 2, we review the Monte Carlo wavefunction method and how it is implemented in order to compute the quarkonium evolution. In section 3, we discuss the initial conditions. In section 4, we present our results and compare with what was previously obtained by directly solving the Lindblad equation. In section 5, we outline the hydrodynamical evolution used for the QGP background. In section 6, we compare with the data collected at LHC and finally, in section 7, we give our conclusions.

The Monte Carlo wave-function method
In this paper, we focus on the evolution of a heavy quark-antiquark pair that follows the GKSL or Lindblad equation [35,36] derived in [32,33] in the regime in which all the thermally induced energy scales are much smaller than the inverse Bohr radius and much larger than the binding energy E. The general form of the Lindblad equation is where ρ is the reduced density matrix, H is a Hermitian operator and the operators C n 's are called collapse operators. The specific values of H and C n in our case will be discussed later.
In [32,33], eq. (2.1) was solved numerically by expanding in spherical harmonics, keeping only S-wave and P -wave states and discretizing the radial component on a lattice. This is a very computationally demanding procedure that limits phenomenological applications. This is a common problem when solving the Lindblad equation, which is due to the fact that the size of the density matrix scales with the square of the number of degrees of freedom. In our case, for example, this means that doubling the lattice size implies multiplying by four the computational cost. In the literature, several techniques have been developed to tackle this problem, and they are named master equation unraveling (see [50] and references therein). One such unraveling, the Quantum State Diffusion method [51], has been already applied to the study of quarkonium, see for example the recent papers [39,40]. However, we believe that a more robust method is provided by the Monte Carlo wave-function (MCWF) method [38]. The reason is that it makes a more efficient use of the symmetries of the evolution of quarkonium, particularly when color degrees of freedom are taken into account. For the case of the equations derived in [32,33], it allows to take into account all possible orbital angular momenta using only a one dimensional lattice in the numerical computation.
Let us now review the basis of the MCWF method. First, we define a partial decay width operator as Then, the total decay width is simply given by Γ = n Γ n . We can also define an effective Hamiltonian (which is non Hermitian whenever Γ = 0) We can now rewrite the Lindblad equation as Note that in QCD H eff does not mix states with different color or angular momentum. Now, we consider the evolution of ρ(t) during a small time step dt the density matrix ρ is a Hermitian semi-positive definite matrix with trace equal to one. Then ρ(t) = n p n |Ψ n (t) Ψ n (t)|, where p n ≥ 0 and n p n = 1. Because the Lindblad equation is linear in the density matrix, if we are able to solve the evolution for an initial condition that is a pure state then we can solve the evolution for any initial condition.
We can understand the previous equation in the following way.
(a) With probability 1 − Ψ(t)|Γ|Ψ(t) dt the wave-function follows the evolution This evolution can be computed by performing the operation (1 − iH eff dt)|Ψ(t) and then normalizing the result.
(b) With probability Ψ(t)|Γ n |Ψ(t) dt the wave-function makes what is called a quantum jump, and transforms into C n |Ψ(t) Similarly to the case of no jump, the evolution can be computed by performing the operation C n |Ψ(t) and then normalizing the result.
We see that the evolution of a density matrix that fulfils the Lindblad equation is equivalent to making the previously discussed operations to a wave-function with the corresponding probability. This can potentially reduce the numerical cost of the simulation of the evolution because, although the algorithm involves averaging over many trajectories, the cost of computing each trajectory scales like N (the number of entries in the wave-function), while directly solving the Lindblad equation scales like N 2 . There is an additional overhead due to the requirement of sampling a large enough number of trajectories, which however does not vary strongly with N , see e.g. appendix E, where the box size N has been varied by a factor two with the same number of trajectories 2 . How many trajectories are a large enough number may depend on the state in question. Generally more trajectories are necessary for reliably sampling final states whose l value is further away from the l value of the initial state. This property is shared with the Quantum State Diffusion method, however, the advantages of the MCWF will become clearer when we review the specific form of the Lindblad equation derived in [32,33].
2.2 Quarkonium evolution in the regime 1/a 0 T, m D E Here we review the evolution equation obtained from pNRQCD in [32,33] in the regime 1/a 0 T, m D E for Coulombic quarkonia. The evolution equation has the Lindblad form of eq. (2.1). In our case, ρ is the reduced density matrix of the heavy quark-antiquark pair. We assume it to be block diagonal in the color degrees of freedom such that The symmetries of the evolution equations ensure that if we start with a density matrix with this structure at some given time t 0 it will keep the same structure during all the evolution; ρ s and ρ o are formally operators in a continuous Hilbert space labeled by the relative coordinate r. The probability for the quarkonium to be in a specific color state is p s,o = Tr(ρ s,o ). The Hamiltonian is given by , where h s,o is the in vacuum Hamiltonian of a singlet or an octet in pNRQCD [19,20]: p is the relative momentum of the QQ pair, C F = (N 2 c − 1)/(2N c ) and N c is the number of colors. We assume the quarkonium under examination to be a Coulombic bound state, hence the potential in (2.9) is the Coulomb potential in the color singlet (attractive) and color octet (repulsive) representation. The coefficient γ is obtained by matching pNRQCD with QCD; it can be expressed as the imaginary part of the integral of a chromoelectric correlator: (2.10) where · · · stands for the in medium average and T for time ordering. The field E(t, 0) should be understood as in [32,33], i.e., as Ω † (t) E(t, 0) Ω(t), where now E(t, 0) stands for the usual chromoelectric field in QCD and Ω(t) is a Wilson line going from −∞ to t: . The Wilson lines ensure that γ is gauge invariant.
Regarding the collapse operators, there are six of them: C 0 i and C 1 i (the subindex i corresponds to the spatial directions and may assume the values 1,2,3). They read The coefficient κ is the heavy quark momentum diffusion coefficient [53,54]. Like γ it is obtained by matching pNRQCD with QCD; it can be expressed as the real part of the integral of a chromoelectric correlator: (2.13) The chromoelectric fields are defined as in the case of γ, hence also the above expression of κ is gauge invariant. Both γ and κ depend on the medium and, because the medium is evolving, on time. In the case of a medium that is in local thermal equilibrium, we can define a temperature and encode in it the time dependence of γ and κ. How the temperature depends on time relies on the hydrodynamical description of the medium. Under the assumption of local thermal equilibrium, lattice determinations of γ and κ at different temperatures coupled to the hydrodynamical evolution of the medium can, therefore, describe how these coefficients evolve in time.
The collapse operators define two partial decay widths (see eq. (2.2)) and Their sum gives the total decay width, Let us now discuss the general structure of the density matrix regarding the orbital angular momentum. We can always expand the density matrix in terms of spherical harmonics: ρ lm;l m = dΩ(r) dΩ(r ) Y lm (r) ρ Y l m * (r ) . (2.17) As stated above, we assume the density matrix to be block diagonal with respect to color (see eq. (2.7)). We further assume that the same happens regarding the orbital angular momentum quantum numbers l and m. Hence, we only need to consider the case l = l and m = m. Moreover, spherical symmetry enforces that all polarizations are equally possible.
the reduced density matrix can be written as ρ = l ρ l and we can understand Tr(ρ l ) as the probability that the value of the angular momentum squared is l(l + 1). The main difference with the studies in [32,33] is that in this work, thanks to the application of the MCWF method, we do not need to truncate the sum of the quantum number l, while in previous works only S-wave and P -wave states were considered. In summary, the MCWF algorithm applied to the case of quarkonium evolution in the 1/a 0 T, m D E regime consists of the following steps: 1. Write the density matrix at the initial time as a sum of pure states The symmetries of the problem ensure that each of these pure states will have a well defined color and angular momentum.
2. With probability p n take |Ψ n (t 0 ) as the initial condition.
3. For the first time step, follow the recipe given in section 2.1 to determine whether there is no jump or, if there is a jump, which collapse operator has to be applied.
(a) If there is no jump, the evolution is given by a Schrödinger equation. Since neither |Ψ n (t 0 ) nor H eff mix different colors or angular momenta, this is a problem that can be solved numerically using a 1D lattice.
(b) If there is a jump, the collapse operators can change color or angular momenta. However, in the case of the Lindblad equation considered here this is done in such a way that the resulting state also has a well defined color and angular momentum (although different from the previous one). The probability of transition between different angular momentum states are discussed in appendix A. Regarding color, a color singlet QQ state always jumps to a color octet QQ state. In the case that the state that jumps is an octet, it has a 2/(N 2 c − 2) chance to jump to a singlet. 3 4. Repeat step 3 for each time step until the end of the evolution.

The waiting time approach
In the numerical implementation of the algorithm, we use an approach that we call the waiting time approach (see section IIID of [55] and references therein). Let us consider the case in which the initial wave-function is |Ψ(t 0 ) . The jump rate at time t is 19) which means that at the time t the norm of the state evolving according to H eff is equal to the probability of no jumps up to the time t. With this in mind we can make a more efficient algorithm to compute one trajectory.

2.
We evolve the wave-function with the effective Hamiltonian H eff until the norm is equal to the random number or smaller.
3. We compute the effect of the jump. 3 If |Ψo is the wave-function of an octet state, then Ψo|Γ 0 |Ψo 4. We repeat the process with the wave-function resulting from the jump.
In this way, we only need to generate a few random numbers per trajectory.

Algorithm for updating the wave-function between quantum jumps
For the evolution of the wave-function between jumps we use a split-step pseudospectral method [13,56,57]. Due to the spherical symmetry of the underlying potentials we can compute the evolution of a state with a given angular momentum quantum number l in a color representation c using where u(r, t) = rR(r, t) with R(r, t) being the radial part of the wave-function. The Hamiltonian operator H l,c contains the angular momentum term and a potential that depends on whether the pair is in a color-singlet or a color-octet configuration. For the normalization of the various states, we take with L being the upper bound on the radius in the simulation.
To enforce the boundary condition on u(r, t) at the origin, we use real-valued Fourier sine series in a domain r ∈ (0, L] to describe both the real and imaginary parts of the wave-function. To perform the update specified in eq. (2.20) we split the Hamiltonian into kinetic and potential contributions H l,c = T + V l,c and use the Baker-Campbell-Hausdorff theorem to approximate, using a split-step decomposition, the time evolution operator as The resulting time evolution steps are 1. Update in configuration space using a half-step: ψ 1 = exp(−iV ∆t/2)ψ 0 .
6. Repeat until the next jump is triggered.
The discrete sine transforms (DST) above can be implemented using standard routines for Fast Fourier Transforms (FFTs). By repeating this procedure, we can evolve the wavefunction forward in time in a manner that is manifestly unitary for real-valued potentials. For the final results reported herein, we take the momentum-space operator to be T = p 2 /m where m is the mass of the heavy quark; however, for testing we also consider a kinetic energy operator that uses centered discrete differences to compute the second derivative.
In this second case, the momentum-space representation of the kinetic energy operator is of the form T discrete = 2[1 − cos(p∆r)]/(m∆r 2 ). We note that besides the improved performance observed, one major benefit of the DST algorithm is that, when dealing with the derivative operator, T , the derivatives are effectively computed using all points in the lattice, not just a fixed number of them. As a result, the evolution obtained using the DST algorithm is more accurate than the one one would obtain using, for example, a Crank-Nicolson (CN) scheme with a three-point secondderivative [13]. Using the same sized derivative stencil as the default DST algorithm, the CN scheme scales as O(N 2 ) where N is the number of lattice points in the wave-function, whereas the split-step DST evolution scales as O(N log N ) [13]. In order to maximize the calculation speed, we use the CUDA Fast Fourier Transform (CUFFT) library which leverages massively parallel graphics processing units (GPU) to compute the necessary FFTs more efficiently [58]. As an additional optimization, we have implemented batched FFTs in order to simulate hundreds of quantum trajectories simultaneously, which maximizes GPU utilization. When run in computing environments containing multiple high-end GPUs, e.g. NVIDIA Tesla P100, the resulting quantum trajectory code allows us to efficiently simulate millions of quantum trajectories in parallel. The resulting code is called QTraj.

Initial production
In this section, we discuss the initial conditions that we use in order to solve the Lindblad equation. The authors of [32,33] assumed that the initial density matrix is diagonal in color and consistent of only S-wave states. We are going to assume the same here, except when studying the evolution with initial P -wave states. In [32,33], it was additionally assumed that the initial wave-function was a delta function in coordinate space. The physical reason is the following. The production of a heavy quark-antiquark pair is a process that involves energies of the order of m or higher. In coordinate space this means that the process is localized in a region with a radius of the order of 1/m or smaller. This is tiny compared to the size of the typical wave-function of a bound state, which is related with the Bohr radius a 0 ∼ 1/(mv). This is the same reason why in NRQCD studies of production and decay the leading order result for S-waves are only sensitive to the wave-function of the bound state at the origin.
However, using a delta function for the initial state in the MCWF method is problematic. The reason is that a wave-function proportional to a delta function cannot be properly normalized. A solution consists in regularizing the delta function on a finite lattice [59]. However, this solution is still problematic when using the MCWM. The reason is that the delta function regularized on a finite lattice contains momenta much higher that the ones that are relevant for the bound state physics. The wave-function of high energy modes grows faster than that of low energy ones and, since the decay width increases with distance, higher energy modes undergo more quantum jumps. Starting with a delta function regularized on a finite lattice generates many trajectories of high momentum that jump several times and have no impact on quarkonium physics. This leads to the fact that the number of trajectories that need to be simulated to obtain precise results becomes extremely large. The situation gets worse as the lattice spacing is reduced since even higher momentum modes enter in the delta function.
To overcome the above issues, we use an initial wave-function in which very high momentum modes are cut off, but whose size is still much smaller than the Bohr radius a 0 . This can be achieved using a Gaussian with a small width instead of a delta function. To fix the notation, we consider the Gaussian wave-function where N can be computed by demanding that the state is normalized. In figure 1 we compare the results of [32,33] for R AA with the ones that are obtained using the same code and parameters but with a Gaussian instead of a delta function regularized on a finite lattice as initial condition. 4 Note that the agreement increases as c is reduced, but it also makes it harder to implement in the MCWF method. We find that c = 0.2 is a good compromise, and this is the value that we will use in the following. c=0.5 c=0.2 c=0.1 c=0.05 Figure 1: Relative difference in R AA , computed as in [32,33], when using a Gaussian regulated delta function and a finite lattice regulated delta function defined according to [59].
In this work, we also study the survival probability of P -wave states that are produced at collision time. The same arguments that we used in the S-wave case apply also here. Production takes place in a very small region compared to the size of the bound state, but this time the angular distribution is that of a P -wave state. Then, at a resolution comparable with the size of the quarkonium, the wave-function after production behaves as a derivative of a delta function. This has the same problems as the delta function in the S-wave case and the solution is analogous: use instead of the derivative of a delta function the derivative of a Gaussian, which happens to be proportional to the Gaussian multiplied by r.

Comparisons between QuTiP and QTraj
Previous studies [32,33] made use of the open-source QuTiP 2 Python package [60,61] to solve the Lindblad equation with the Hamiltonian H of eq. (2.8) and the collapse operators C n i of eqs. (2.11) and (2.12) truncated at l = 1 in the spherical harmonic expansion discussed in section 2.2. The newly developed QTraj code that we present in this work possesses several distinct advantages over the previously utilized QuTiP code. 5 Foremost among these is, for a system simulated with N discrete points, the reduced size in memory from O(N 2 ) to O(N ), due to calculating with the wave-function rather than with the density matrix. It should be also added that in QTraj many single calculations of O(N ) must be performed and averaged to arrive at a final result with sufficiently small statistical error; the overall computational complexity 6 may, therefore, exceed O(N ). However, as each trajectory is independent, this can be counterbalanced by running trajectories in an embarrassingly parallel setup. Other advantages of the QTraj code include working to infinite order in the orbital quantum number l and the use of an all-points derivative rather than forward-backward finite differences. Compared to the QuTiP based code developed for refs. [32,33] QTraj also includes the capability to couple to realistic 3+1D hydrodynamical backgrounds.
We can use the QuTiP code as a benchmark for the QTraj code. In order to do so, first we set up the QTraj code as the QuTiP code by implementing a cutoff in l, a forwardbackward finite difference derivative, T discrete , and Bjorken evolution. Finally, we compare results obtained in this way with the QTraj code with the results obtained with the QuTiP code.
We perform two series of tests to establish the agreement of the results obtained from the new QTraj code with the results obtained from QuTiP. Our procedure is to simulate the evolution of the state from t = 0 fm to t = 0.6 fm in the vacuum and from t = 0.6 fm to t = 2.95 fm in a medium of initial temperature T 0 = 425 MeV undergoing Bjorken evolution. On the QuTiP side, we solve the Lindblad equation and use the time evolved density matrix to calculate the expectation values of the 1S, 2S, 3S, 1P , and 2P Coulombic states and of the radius r. On the QTraj side, we run 196608 trajectories and calculate the corresponding expectation values; finally we compare with the QuTiP results. We perform these tests on a lattice of spatial extent L = 51.2 a 0 and lattice spacing a s = 0.1 a 0 . We take κ/T 3 = 2.6 and γ = 0. In QuTiP, we work with angular momentum cutoffs l max = 1 and l max = 2; in QTraj, we implement an angular momentum cutoff l max and run simulations 5 We note that a Monte Carlo solver is also implemented in the qutip.mcsolve function of QuTiP; it still requires, however, the implementation of a cutoff in l in contrast to the QTraj code. 6 A systematic scaling study of computational complexity and number of trajectories needed will be included in a forthcoming publication by some of the authors [62]. with l max = 1, 2, and ∞. There is a subtlety to consider when working with a cutoff in angular momentum in the quantum trajectories algorithm. The implementation of a cutoff in l in QTraj leaves the evolution of states of orbital angular momentum l < l max unaffected, i.e., they proceed according to the prescriptions of section 2.2. However, a state of angular momentum l = l max must be evolved with a reduced width cf. eq. (2.2). The reason for this is that a state of angular momentum l max can only jump down to a state of angular momentum l max − 1. Hence, the total width of the state is reduced by an amount equal to the probability of jumping down by one unit in angular momentum with respect to the width of the state without cut off on the maximal orbital angular momentum. This probability is given by P lmax d = l max /(2l max + 1), see appendix A. It is therefore the reduced width, Γ = Γ n with Γ n given by eq. (4.1), that enters the effective Hamiltonian describing the evolution of the state with l = l max . Once a jump is triggered, this is deterministically a jump down. The color evolution remains unchanged. Under the above conditions, we test the two programs using identical initial conditions.
In the first test, we run both programs with a Coulombic 1S wave-function as initial condition and plot the results in figure 2. We observe excellent agreement for all measured quantities at identical angular momentum cutoffs l max .
In the second test, we run both programs with a Gaussian of width c = 0.2 as initial condition, see section 3. We plot the results in figure 3. In this case, we do observe some disagreements at late times that cannot be fully accounted for by the statistical errors of QTraj. The position space wave-function of a narrowly peaked Gaussian quickly expands, and for our simulation parameters reaches the outer edge of the box before the end of the simulation. The reflection behavior exhibited at the spatial boundary of the lattice, which may be deduced also from the shrinking of the average radius after an initial expansion in the last plot of figure 3, differs between QuTiP and QTraj. We attribute the late time (small) discrepancies in the expectation values between the two programs to this phenomenon. Since we find that a Gaussian wave-function evolved in QuTiP hits the edge of the box at approximately t = 1.3 fm, we take t > 1.3 fm as the region where finite size effects become significant. When running the QTraj code to simulate our final results for the quarkonium nuclear modification factors in section 6, we will avoid finite size effects by increasing the spatial extent of the lattice. We note that while increasing the extension of the lattice is rather cheap for QTraj, it is computationally unfeasible for QuTiP.
In summary, we find that the results of the QuTiP and QTraj programs are in excellent agreement when run in their regions of validity with identical parameters. They agree on both the extracted overlaps and the expectation value of r including their dependence on l max . This gives us confidence to proceed with the use of QTraj for our calculations.

Hydrodynamic background evolution
In order to compute the in-medium survival probability of a given quantum trajectory one must specify the temperature evolution of the QGP. For this study we make use of a 3+1D dissipative hydrodynamics code that is based on the quasiparticle anisotropic hydrodynamics (aHydroQP) framework for dissipative relativistic hydrodynamics [63][64][65]. The aHydroQP framework has been shown to well-reproduce a variety of experimental soft-hadronic observables such as the total charged hadron multiplicity, identified hadron spectra, integrated and identified hadron elliptic flow, and HBT radii at both RHIC and LHC nucleus-nucleus collision energies [66][67][68][69][70][71]. The aHydroQP framework is an extension of the originally formulated conformal anisotropic hydrodynamics [72][73][74]  The gray band denotes the area where the wave-function hits the edge of the lattice, i.e., finite size effects become relevant, which occurs at t = 1.3 fm (see text for discussion). We note that for the S-wave states we report their overlaps normalized to their initial overlaps and denote this quantity R AA .
effects of non-conformal transport coefficients such as the bulk viscosity. The resulting code uses a realistic equation of state determined from lattice QCD measurements [75] and self-consistently computed second-and higher-order transport coefficients. Due to the resummation to all orders in the inverse Reynolds number, the anisotropic hydrodynamics framework can be applied at early times after an AA collision, which is a period of time when strong non-equilibrium effects are present [76][77][78][79][80][81][82][83][84][85][86][87]. This allows us to more accurately describe the dynamics of the QGP during the entire evolution of the produced bottomonium states since they are created at early times in the QGP's lifetime.
For this paper we use the aHydroQP tuning recently reported in ref. [71]. We assume a smooth optical Glauber initial condition that provides the initial spatial energy density profile of the QGP as a function of the impact parameter. The study of ref. [71] concluded that the best fit to soft hadron observables is achieved using an initial central temperature of T 0 = 630 MeV and a constant specific shear viscosity of 4πη/s = 2. 7 To determine the temperature experienced by bottomonium states produced in the QGP, we assume that the initial transverse spatial distribution for bottomonium production is proportional to the binary overlap profile of the two colliding nuclei, N bin AA (x, y), and then use Monte-Carlo sampling to generate the initial production points. We take the initial transverse momentum (p T )-distribution to be proportional to p T /(p 2 T + M 2 ) 2 for all states, where M is the average mass of all states being considered. We then Monte-Carlo sample the p T for each particle generated. We take into account the approximate boost-invariance of the QGP and assume all bottomonia to have zero momentum rapidity, y = 0. Finally, we sample the initial azimuthal angle φ from a uniform distribution between 0 and 2π. We then averaged the temperature obtained along each Monte-Carlo generated path using approximately 132000 samples per centrality bin. The resulting path-averaged temperature evolution is plotted in figure 4 for each centrality bin used in this work. Note that at late times (t 3 fm) one can see the onset of 3D expansion for centralities 50%, 7 The initial longitudinal proper time corresponding to this initial central temperature is 0.25 fm [71]. which results in more rapid cooling of the QGP. 8 In table 1, for each centrality class considered, we list the average impact parameter, the average number of participating nucleons, the initial central temperature, and the path-averaged initial temperature. We note that since many bottomonium states are created away from the center of the collision, the path-averaged temperature is always lower than the central temperature. We evolve the quantum wave-packets using the in vacuum potential starting at t = 0 fm and turn on the in medium complex potential at t = 0.6 fm. Finally, when the averaged temperature drops below T f = 250 MeV, we again use the in vacuum potential for the wave-function evolution. Hence the hydrodynamical evolution of the medium does not play any role between t = 0 fm and t = 0.6 fm and after bottomonium freeze out. By neglecting medium effects below T f we are ignoring physical effects that may be relevant specially for excited states. The inclusion of physical effects in the regime T ∼ E is beyond the scope of this paper. It would involve solving a more complicated master equation [33] in which the information of the medium cannot be encoded in two transport parameters, as it is the case in the regime that we study in this work. In Appendix D, we try to assess the size of the physical effects we are leaving out by varying the value of T f .

Results and discussion
In this section, we present the QTraj results for the survival probability of various bottomonium states and quantify the effect of quantum jumps on this observable. Using the resulting survival probabilities, we then include the effect of late-time feed down of excited states and compare the QTraj results with available LHC data for R AA and double ratios of various states.

Parameters
We look at the quarkonium evolution in the QGP in the regime 1/a 0 T, m D , Λ QCD E; we further assume that the quarkonium is mostly Coulombic. 9 This has been discussed in section 2.2. From that discussion it follows that the evolution depends on four parameters. One parameter is the bottom quark mass m b , which has to be understood as the pole mass. We take the bottom quark mass to be m b = 4.881 GeV. Another parameter is the strength of the Coulomb potential. We may trade this parameter for the Bohr radius a 0 , which we take to be a 0 = 0.742 GeV −1 = 0.146 fm to reproduce the value used in [32,33]. This value is the solution of the self consistency equation a 0 = 2/(C F α s (1/a 0 )m b ) when the strong coupling α s in the MS scheme is taken at one loop accuracy. The strong coupling α s at the scale of the inverse of the Bohr radius is then α s (1/a 0 ) = 0.414. Finally, the evolution equations also depend on the two coefficients γ (see eq. (2.10)) and κ (see eq. (2.13)). The coefficients γ and κ have mass dimension three. Hence it is convenient to definê The coefficient γ has been taken equal to zero in [32,33]. More recently, in [34] γ has been extracted using 2+1 flavor lattice QCD data for the quarkonium mass shift in a thermal bath. Note that γ and κ are flavor independent. It was found that the lattice data for the J/ψ and Υ(1S) thermal mass shifts at the temperatures 251 MeV and 407 MeV (Υ(1S) case only) fall inside the interval −3.8 <γ < −0.7. In this work, we will takê γ = −1.75 ± 1.75 (6.2) to encompass also the value used in [34]. On general grounds we expect thatγ depends on the temperature. Indeed, the data in [34] seem to suggest thatγ is closer to zero at higher temperatures. Nevertheless, our present knowledge ofγ is clearly insufficient to parameterizeγ in terms of the temperature. Hence, we will assumeγ constant (in the temperature, and, therefore, also in time) in the present analysis. The coefficient κ has been recently computed at temperatures in the range 1.1 T /T c 10 4 , where T c = 155 MeV, from pure SU(3) gauge lattice data in [37]. Because of the wide range of temperatures, it has been possible to parameterize the temperature dependence ofκ. The change ofκ with the temperature is well parameterized by its nextto-leading order expression when the coefficient of the term proportional to the ratio of the Debye mass over T is fitted. The parameterization ofκ as a function of the temperature with the corresponding error band is shown in figure 5. In this work, we will takeκ and its uncertainties according to figure 5. We will callκ C (T ) the central line,κ U (T ) the upper 9 In this work, we investigate the Υ(1S), Υ(2S), Υ(3S), χ b (1P ) and χ b (2P ) under these assumptions. In particular, we remark that the Υ(1S) is commonly treated as a Coulombic bound state, whereas the Υ(2S), χ b (1P ) and more critically the Υ(3S) and χ b (2P ) have been investigated as Coulombic bound states, for instance, in [88][89][90][91][92][93]. boundary andκ L (T ) the lower one. Recall that, once coupled with the hydrodynamical evolution of the QGP, the temperature dependence ofκ translates into a time dependence.
We emphasize that in the regime 1/a 0 T, m D , Λ QCD E all parameters describing the quarkonium evolution in the QGP are determined independently of the quarkonium nuclear modification factors that we eventually compute. Indeed, m b and α s are parameters of the QCD Lagrangian, and the coefficients γ and κ are determined from first principle (lattice) computations in QCD.

Numerical setup
For all results reported in this section we use a lattice size of N = 4096 points with L = 80 GeV −1 ≈ 108 a 0 . The temporal step size for evolving the wave-function between jumps is ∆t = 0.001 GeV −1 ≈ 2 × 10 −4 fm. We initialize the wave-function at t = 0 and evolve it with the vacuum potential until t = 0.6 fm. From this point forward in time, in between quantum jumps, we evolve the wave-function with the potential appropriate for the considered QQ state labeled by its integer orbital angular momentum l ≥ 0 and color state (singlet or octet). We terminate the evolution when the path-averaged temperature, T average , in a given centrality bin drops below T = T f = 250 MeV. We then compute the final quantum mechanical overlaps of the vacuum eigenstates with the QTraj evolved wave-function to obtain the survival probability of the state.
For all results presented in this section we initialize the wave-function using a smeared Gaussian delta function (3.1) of the form u l (r, t = 0) = N r l+1 exp(−r 2 /(ca 0 ) 2 ) , with c = 0.2 and N being fixed by the normalization condition (2.21).

Bottomonium survival probabilities
We present, first, results for the survival probability of various bottomonium states. The survival probability tells us the probability to find a given quantum state after the wavefunction is evolved in the QGP. This measure does not yet take into account late-time feed down of excited states, but does allow us to assess the impact ofκ andγ. In addition, we can use the pre feed down results to more cleanly quantify the effect of quantum jumps on the resulting survival probabilities.   demonstrates, our resulting statistical errors are smaller than the systematic uncertainties coming from the choice ofκ andγ. We also find that the QTraj predictions for the survival probability are sensitive to the choice ofκ andγ, with the variation ofγ resulting in the larger variation of the survival probability for the Υ(1S). The eventual comparison of the quarkonium nuclear modification factors with the experimental data will, therefore, be in the position to validate or invalidate our independent choice of values for the coefficientsκ andγ given respectively in figure 5 and eq. (6.2).
In order to quantify the effect of quantum jumps on the S-wave survival probabilities, in figure 7 we present five panels where we compare the result of evolving the system with no quantum jumps to the full QTraj result. In the case that no jumps are allowed, this reduces to evolving the wave-function solely with the complex Hamiltonian H eff . For the full QTraj results including quantum jumps, we once again use 98304 quantum trajectories for each point in N part . In the top row of figure 7, we present the results obtained when varyingκ in the same range used in the left panel of figure 6. In the bottom row of figure 7, we present the results obtained when varyingγ in the same range used in the right panel of figure 6. From these figures we firstly note that the survival probability of the Υ(1S) is well reproduced by the H eff (no jump) evolution. For the excited states, we see a larger effect from the correct implementation of the quantum jumps. Focusing on the top row of figure 7, we see that the effect of jumps on the excited states increases with increasingκ. We also notice that increasingκ generally decreases the survival probability of the states. The importance of quantum jumps is largest in the case ofγ = 0, which is shown in the bottom left panel of figure 7. In this case, the H eff evolution (dashed lines) predicts strong suppression of the 2S and 3S states. When this occurs, the corrections from the jumps become comparatively more important and result in a flattening of the survival probability as a function of N part at a magnitude that is similar to the survival probabilities of the 2S and 3S states seen in other panels of figure 7.
One of the reasons that we see small effects of quantum jumps on the singlet S-wave survival probabilities when usingκ andγ at their central values is that, once such an initial state makes one quantum jump, it necessarily changes its angular momentum quantum number l and its color state to octet. Later quantum jumps are more likely to cause the state to jump to color octet states with higher angular momentum. For example, starting from a singlet 1S state one can only jump to an octet 1P state. Once the state is in a color octet configuration, the potential becomes repulsive and the wave-function will spread to larger radii. 10 After evolving the wave-function, a second jump can occur that causes the state to jump back down to the singlet 1S state, however, this is a disfavored transition since the probability of an octet to singlet transition is 2/(N 2 c − 2) = 2/7 (see footnote 3) and the probability of an l = 1 to l = 0 transition is P 1 d = 1/3 (see eq. (A.8)), resulting in a combined probability of 2/21 ≈ 10% for transitioning from a P -wave color octet to an S-wave color singlet. Moreover, the compact 1S color singlet state overlaps only marginally with the wide color octet state and the impact of the jump on the evolution equation is small. The other 90% of the time the state will transition to an l = 2 singlet or octet state or to an l = 0 octet state. Since for any finite l transitioning to higher l is favored (see eq. (A.8)), the result is, on average, a directed random walk towards higher l. As the state transitions to higher and higher angular momentum states, the probability to jump up or down in l approaches P ∞ d = P ∞ u = 0.5 and the state eventually does a balanced random walk in l. The net result of all this is that the survival probability for the singlet 1S state is strongly dominated by the case of no quantum jumps and that the corrections, while important to quantify, are small. For the 2S and 3S states, one sees larger effects from the quantum jumps because they have larger average radii than the 1S state and hence will have a larger overlap with color octet states after a series of jumps. Nevertheless, also in this case transitioning to higher l is favored for the same argument exposed above.
In order to assess what role the repulsive nature of the octet potential plays in the importance of quantum jumps, in figure 8 we present QTraj results obtained when promoting in the code all color octet states to color singlet ones. For ease of comparison, we have made the vertical scale in figure 8 the same logarithmic scale that has been used in figure 7. Comparing the central panel in the top row of figure 7 with figure 8, we see that the effect of quantum jumps is enhanced in particular on the excited states, when using a framework that includes only the attractive color singlet potential. This enhancement is due to the fact that the singlet-potential is attractive, which causes states to have smaller average radii than when they interact through a repulsive color octet potential. A smaller radius leads to a larger overlap with S-wave states and eventually to a larger probability to jump back to a state with lower angular momentum.  Singlet delta P -wave initial conditions In figure 9, we show the QTraj results obtained with color singlet P -wave (l = 1) initial conditions. As in figure 6, the left and right panels correspond to varyingκ andγ while holding the other coefficient fixed. We have averaged over approximately 393216 quantum trajectories at each point in N part , and the statistical errors are indicated by error bars on the central line for each state. As can be seen from these figures, there is a stronger variation withγ than withκ. In the right panel, the upper bound for the P -wave survival probabilities has been obtained when usingγ = 0, with the other two values considered resulting in much stronger P -wave suppression.

Off-diagonal overlaps and octet initial conditions
In the previous two subsections, we presented the singlet S-wave overlaps resulting from S-wave delta initial conditions and the singlet P -wave overlaps resulting from P -wave delta initial conditions. Due to the fact that the full QTraj evolution can mix different angular momentum and color states, one can also consider, for example, (i) the singlet S-wave overlaps resulting from P -wave delta initial conditions, (ii) the singlet P -wave overlaps resulting from S-wave delta initial conditions, and (iii) the singlet S-wave overlaps resulting from octet P -wave initial conditions. In all three cases, we find that these contributions to the final overlaps are small and can be ignored in phenomenological applications. We provide details concerning our findings in appendix B.

Final results including late-time feed down of excited states
Once each quantum state is evolved using QTraj, the survival probabilities are converted into particle numbers by multiplying by (a) the expected number of binary collisions in the centrality bin sampled and (b) the direct production cross section for each bottomonium state. After the number of states that survived transversal of the QGP is computed, one needs to take into account the late time feed down of excited bottomonium states. Here, we follow refs. [15,16] and introduce a feed down matrix F which collects the known information about excited bottomonium state decays available from the Particle Data Group [94].
In the case of pp collisions, one can convert the direct production cross sections into the post feed down cross sections by multiplying a vector containing them by the feed down matrix σ exp = F σ direct with where the vectors σ collect the experimentally-observed and direct cross sections for the  [95]. For the Υ(1S), Υ(2S), and Υ(3S) cross sections, we use the 5.02 TeV data obtained by the CMS collaboration in the rapidity interval |y| ≤ 2.4 [96]. In ref. [96] the left panel of figure 3 presents B × dσ/dy, where B is the dimuon branching fraction. Averaging over rapidity in the interval presented in the CMS figure, we obtain B × dσ/dy ≈ 1.44 nb, 0.37 nb, and 0.15 nb, respectively. Dividing by the branching fractions for Υ(1S), Υ(2S), and Υ(3S) → µ + µ − , which are ≈ 2.5%, 1.9%, and 2.2%, respectively [94], one obtains To compute the effect of final-state feed down in AA collisions, we first construct a vector N QGP containing the numbers of each state produced at the end of each simulation (survival probability × N bin (b) × σ direct ). We then multiply the result by the same feed down matrix used for pp feed down, i.e. N final = F N QGP . The use of the same feed down matrix for both pp and AA collisions is related to the fact that feed down occurs on a time scale that is much longer than the QGP lifetime. After the feed down is complete, we compute the post feed down R AA for each state by dividing the final number of each state produced by the average number of binary collisions in the sampled centrality class times the post feed down pp production cross-section for that state (σ i exp ), giving where i labels the state of interest, S is a diagonal matrix that contains the quantumtrajectory averaged survival probabilities for each state along the diagonal, and c indicates the centrality class considered. The resulting R i AA (c) contain the suppression factors for all states included in our analysis. Since, with respect to the Hamiltonian in our QTraj simulations, states belonging to the same spin multiplet are degenerate, we take the survival probabilities of these states to be the same, i.e. S[χ bj (nP )] = S[χ b (nP )] when constructing the survival probability matrix S.
In figure 10, we present our final results for the suppression of the Υ(1S), Υ(2S), and Υ(3S) as a function of N part compared to experimental data available from the ALICE [100], ATLAS [101], and CMS [96] collaborations. The bands shown on the QTraj results represent the variations with respect toκ (left) andγ (right) with the systematic uncertainties propagated through feed down. As with the pre feed down survival probability we see a large sensitivity to the choice ofγ andκ. Overall, comparing the central lines with the available data we find quite reasonable agreement between the QTraj predictions and the experimental data, however, for the most central collisions, QTraj seems to show a somewhat stronger suppression than is seen in the data.
It is important to note that the minimum temperature at which we consider in-medium damping of states is T f = 250 MeV. The reason for this choice is that at temperatures lower than 250 MeV the hierarchy T, m D E may not hold, and one needs to solve a different set of evolution equations [33]. As a result, QTraj predicts that R AA = 1 for N part 9.5 (see table 1), when using the path-averaged temperature. Nevertheless, damping processes do occur at low temperatures. This calls for a future extension of the in medium evolution equations and QTraj to the low temperature region, T c T T f , in order to describe more accurately data at very small N part .
Next, we turn to the double ratios constructed from Υ(2S) and Υ(3S) to Υ(1S). These are obtained by computing the ratio of Υ(nS) and Υ(1S) yields in PbPb collisions, divided by the same ratio in pp collisions; we will indicate this double ratio with [Υ(nS)/Υ(1S)] PbPb / [Υ(nS)/Υ(1S)] pp . On the experimental side, these quantities have typically smaller systematic uncertainties, which allow for tighter constraints on theoretical We compare the QTraj results with data available from the ALICE [100], ATLAS [101], and CMS [102] experiments. For QTraj, we show the theoretical uncertainty as a shaded blue band; for the experimental systematic and statistical uncertainties we use black and red error bars, respectively. As can be seen from the figure, QTraj is within 1σ of the statistical error bars. We see visible deviations at small N part , with these again being due to the fact that we do not allow in-medium breakup at low temperatures as discussed above. By comparing the left and right panels of figure 11, we see that the QTraj results depend more on the assumed value ofγ than on the value ofκ. In all cases considered, QTraj predicts that the 2S to 1S double ratio depends only mildly on N part for N part 150.
Finally, in figure 12 we present the QTraj results for the 3S to 1S double ratio and compare with experimental data from the ATLAS [101] and CMS [102] experiments. In the case of the CMS data, the results were reported as upper bounds on the 3S to 1S double ratio. In the case of ATLAS, the collaboration presented results only for an integrated 2S+3S double ratio. As can be seen from this figure, the QTraj results for the 3S to 1S ratio depend less strongly on the assumed value ofγ than in the double ratio [Υ(2S)/Υ(1S)] PbPb / [Υ(2S)/Υ(1S)] pp . Similar to the 2S to 1S double ratio, we see that QTraj predicts that the 3S to 1S ratio does not depend strongly on N part for N part 150. In the future, increased statistics from the ATLAS and CMS collaborations may allow for more constraints on the 3S to 1S double ratio.

Conclusions
In this paper, we have presented a new program called QTraj to solve the Lindblad equation describing the nonequilibrium evolution of a color singlet and octet heavy quark-antiquark pair of a small radius in a strongly coupled QGP characterized by an inverse correlation length larger than the typical quark-antiquark energy. The evolution equations were originally derived in [32,33] and solved there with the publicly available code QuTiP. The new code QTraj presents several distinct advantages with respect to the one previously used. It shows a sizeable reduction of the memory requirement from O(N 2 ) to O(N ) for a system simulated with N discrete points, due to the fact that it is computing the wave-function instead of directly the density matrix. The need to average over a large number of trajectories is efficiently counterbalanced by running many trajectories in an embarrassingly parallel setup. The need to average over a large number of trajectories is efficiently counterbalanced by running many trajectories in parallel. In addition, QTraj allows to account for wavefunctions of any angular momentum quantum number l and uses all-points derivatives allowing for more accurate solutions of the evolution equations. In particular, with the new code we could eliminate several of the restrictions necessary in [32,33], such as only including l ≤ 1 states and using a Bjorken time evolution for the temperature profile of the plasma.
Many new investigations have become possible that have been performed in this paper for the first time. (i) We included in the solution of the evolution equations the contributions coming from states of all angular momentum quantum numbers l, without the need of introducing a cutoff in this parameter. (ii) We studied also the evolution the in-medium of quarkonium P states. (iii) We coupled the Lindblad equation to a realistic medium evolution based on the quasiparticle anisotropic hydrodynamics (aHydroQP) framework for dissipative relativistic hydrodynamics. (iv) We studied the dependence on the initial condi-tions. (v) We thoroughly investigated the impact of color octet heavy quark-antiquark pairs on the evolution of the quarkonium in the medium. (vi) We quantified the recombination contribution with respect to the effective Hamiltonian evolution, in this way assessing the importance of quantum jumps for heavy quarkonium evolution in the QGP. (vii) We investigated extensively the dependence of the quarkonium nuclear modification factor R AA on the two parameters that characterize the quarkonium evolution in the QGP, i.e. the heavy quark momentum coefficient κ and its dispersive counterpart γ. In particular, we explored for the first time the impact of the temperature dependence of κ/T 3 , as extracted from [37], on the quarkonium evolution and on R AA .
An interesting feature of our approach is the characterization of the quarkonium evolution in the QGP in terms of only two transport coefficients κ and γ, which are, in general, temperature dependent. We have seen that the survival probabilities depend on the values of these parameters. In particular, increasing or decreasing κ results in a decreased or enhanced survival probability for all three S-wave bottomonium states below threshold, while the impact of γ is different from one state to the other. For the Υ(1S), a smaller value of γ results in a reduced survival probability, while for the Υ(2S) the reverse is true. Notable is, in particular, the γ dependence of the P -state survival probabilities. Therefore, the quarkonium survival probability can act directly as a diagnostic of the QGP through these transport coefficients.
We explored and quantified the impact of the quantum jumps, i.e. the recombination effects. We found that quantum jumps seem to only marginally affect the Υ(1S), whose survival probability can be well described using just the effective Hamiltonian evolution. For excited states, however, quantum jumps are found to give a sizeable contribution. The effect increases by increasing the value of κ. We already commented that increasing κ decreases the survival probabilities of the states: quantum jumps correct and mitigate this effect with respect to the pure effective Hamiltonian evolution.
Finally, we included the late time feed down from excited states and we compared our results to the bottomonium nuclear suppression factor measured by the ALICE, ATLAS, and CMS collaborations. It is important to emphasize that the computed bottomonium nuclear suppression factor does not depend on any free parameter as both κ and γ have been determined independently by lattice QCD in ref. [37] and [34], respectively. Hence the comparison is really between a QCD prediction and data. We obtain a good agreement with data, which is better than the one reported in [32,33], confirming the importance of developing the new program in order to lift the limitations of the previous code. The predictions of QTraj have small uncertainties that depend on the allowed values for the κ and γ transport coefficients. In particular, there is a larger dependence on γ that calls for a dedicated lattice study of this coefficient in a wide temperature range. Our formalism, as currently implemented, does not include medium effects at temperatures below 250 MeV. This is an approximation that is accurate up to corrections of relative order (a 0 T ) 2 , which are not included. As already mentioned at the end of section 5, this approximation might neglect physical effects, like thermal gluodissociation, possibly more relevant for excited states, Υ(2S) and Υ(3S), than for the Υ(1S). The reason is that our assumed hierarchy of scales is more marginally realized for the former than for the latter bottomonia. This is visible in our results where the relative error in the survival probability coming from the variation of κ and γ is indeed larger for the excited states. Although the absolute size of the medium effects at low temperatures is presumably small (specially for the case of Υ(1S)), it would be desirable to include these effects in a future work to better address collisions with a very small number of participants and to improve our description of excited states. In a forthcoming work, we will include corrections of order E/T to the currently considered evolution equations, thus extending the range of validity of our method to lower temperatures [103].
We considered also double ratios of S-wave production cross sections in pp and PbPb, which eliminated several of the systematics, both theoretically and experimentally. We predicted that the 2S to 1S and the 3S to 1S double ratios do not depend strongly on N part for N part 150 and we made quantitative predictions for these ratios.
Using QTraj we plan to explore other interesting observables in the near future such as the p T dependence of R AA and v 2 . Moreover, we plan to relax the assumption of isotropy and solve the Lindblad equation in such cases where the κ and γ coefficients become tensors instead of numbers. Longer term goals include to consider quarkonium not at rest with respect to the QGP [104], to extend these studies to quarkonia with larger radius and eventually to solve the full evolution master equations [32,33] far from equilibrium.

A Change of orbital momentum during a quantum jump
In order to determine the probability to jump to any given orbital angular momentum state, we need to study the term of the Lindblad equation that goes like n C n ρC † n . If we consider a density matrix ρ block diagonal in the quantum numbers l and m, then the result of computing n C n ρC † n is also block diagonal. For simplicity, let us consider a density matrix that only contains states with a given orbital momentum. Since the Lindblad equation is linear, it will be trivial to generalize to any other density matrix with spherical symmetry. Then we can split the original term i C x i ρ l C x i † (where x can be either 0 or 1, see notation in section 2.2) in a sum l C x,l→l ρ l C x,l→l † , where each term generates a density matrix with an orbital momentum l . The procedure that we follow in practice is the following. First, in the waiting time approach, we determine if a jump takes place by computing Ψ l |Ψ l , where we explicitly write the orbital angular momentum of the wave-function to clarify the notation. After determining whether the wave-function is affected by the collapse operator that induces singlet-octet transitions, C 0 i , or octet-octet transitions, C 1 i , the system will jump to a state with angular momentum l with probability where Y lm are spherical harmonics and we have taken into account that the system does not have any preferred polarization. We note that the previous quantity is only non-zero if l = l + 1 or l = l − 1. Now we can perform the computation using the following argument. Starting from the identity we can use the fact that the spherical harmonics can be used to expand any function on the unit sphere to deduce that and We can identify P l d as the probability to jump to a state with angular momentum l − 1 and P l u as the probability to jump to a state with angular momentum l + 1. It then holds that From this recurrence relation we determine P l d and P l u knowing that P 0 d = 0:

B Off-diagonal overlaps
In the main body of the paper, we focused on the singlet overlaps resulting from singlet initial conditions with a fixed angular momentum l. For singlet S-wave initial conditions we presented the resulting S-wave overlaps and, likewise, for singlet P -wave initial conditions we presented the resulting P -wave overlaps. In this appendix, we demonstrate why it is consistent to ignore the off-diagonal contributions corresponding to, e.g., singlet Pwave overlaps resulting from singlet S-wave initial conditions and singlet S-wave overlaps from octet P -wave initial conditions. Such off-diagonal contributions are generated during the QTraj evolution due to the quantum jumps, while without quantum jumps (H eff evolution only) such overlaps are identically zero. We will present evidence that, when including quantum jumps, the off-diagonal contributions are small enough as to be ignored for phenomenological applications.  Figure 13: Quantum-mechanical overlaps with singlet S-wave states obtained using singlet P -wave initial conditions. The bands shown in the left and right panels correspond to the variations detailed in the caption of figure 6. In both panels, the central line shows the result averaged over the corresponding variation.
In figure 13, we present the singlet S-wave overlaps resulting from singlet P -wave initial conditions as a function of N part . As in the main body of the text, the left and right panels show the variation over the assumed values ofκ andγ, respectively. In order to gauge the magnitude of these numbers, we note that for a central collision the singlet S-wave overlaps resulting from singlet S-wave initial conditions (corresponding to figure 6) are approximately 6 × 10 −3 , 1 × 10 −4 , and 2 × 10 −5 for the Υ(1S), Υ(2S), and Υ(3S), respectively. Comparing to the magnitude of the overlaps shown in figure 13, we see that the off-diagonal contribution in angular momentum is small in this case. This contribution is additionally suppressed by the fact that the P -wave initial production cross sections are down by a factor of approximately four with respect to the S-wave production cross-section. For this reason we ignore this contribution in our final phenomenological predictions.  Figure 14: Quantum-mechanical overlaps with singlet P -wave states obtained using singlet S-wave initial conditions. The bands shown in the left and right panels correspond to the variations detailed in the caption of figure 6. In both panels, the central line shows the result averaged over the corresponding variation.
In figure 14, we present the singlet P -wave overlaps resulting from singlet S-wave initial conditions as a function of N part . In order to gauge the magnitude of these numbers, we note that for a central collision the singlet P -wave overlaps resulting from singlet P -wave initial conditions (corresponding to figure 9) are approximately in between 10 −4 and 10 −7 , and in between 10 −5 and 10 −8 for the χ b (1P ) and χ b (2P ), respectively, depending on the assumed values ofκ andγ. These are very small numbers, whose effect falls well inside the range of variations considered when varyingκ andγ (see figure 9). We show that the smallness of the off-diagonal contributions is generic in appendix C, where we present the diagonal and off-diagonal overlaps obtained when the initial wave-function is taken to be a pure Υ(1S) or Υ(2S).
In figure 15, we present the singlet S-wave overlaps resulting from octet P -wave initial conditions as a function of N part . In order to gauge the magnitude of these numbers, we note that for a central collision the singlet S-wave overlaps resulting from singlet S-wave initial conditions (corresponding to figure 6) are approximately 6 × 10 −3 , 1 × 10 −4 , and 2 × 10 −5 for the Υ(1S), Υ(2S), and Υ(3S), respectively. As can be seen from figure 15, the S-wave overlaps resulting from octet P -wave initial conditions are orders of magnitude smaller than the overlaps resulting from singlet S-wave initial conditions. For this reason, we can safely ignore the off-diagonal octet-singlet contributions when considering phenomenology.

C S-wave initialization studies
In this appendix, we present results for the overlap factors versus N part obtained with pure S-wave initial conditions. Figure 16 shows results for an Υ(1S) initial condition and  Finally, in figure 18 we present a comparison of the Υ(1S) and Υ(2S) survival probabilities obtained from (i) S-wave eigenstate initial conditions and (ii) Gaussian delta initial conditions as used in the main body of the paper. We see less suppression of the Υ(1S) when using the pure Υ(1S) eigenstate initial condition than when using the Gaussian delta initial condition. For the Υ(2S), however, this pattern is reversed, and one sees a much larger difference between the results obtained using pure eigenstate and Gaussian delta initial conditions. The same pattern is observed when ignoring jumps, i.e. evolving only with the temperature-dependent H eff . The key difference between the two types of initial conditions is that the Gaussian delta is a linear superposition of many vacuum eigenstates. The results indicate, therefore, that quantum state mixing due to the temperature-dependent Hamiltonian is important and leads to a substantial reduction in the suppression of the Υ excited states.

D Effect of varying T f
As discussed at the end of section 5, we evolve the quarkonium in the medium down to a temperature T f , which we take to be 250 MeV, and in the vacumm from T f to the crossover temperature T c . The reason for the first choice is that the evolution equation set in section 2.2 holds under the condition T , m D E. Only under this condition, effects coming from the energy region E may be neglected with respect to medium effects, and medium effects can be cast in the two transport coefficients κ and γ [32][33][34]. The largest Coulombic binding energy is the one of the Υ(1S) state; its value is −1/(m b a 2 0 ), which is about 400 MeV for our choice of parameters. The condition T ≥ T f = 250 MeV guarantees reasonably well that the temperature scale, which is more properly πT , and m D ≈ 2T are larger than 400 MeV during the quarkonium evolution in the medium. The reason for the second choice is that if the temperature becomes as low as the binding energy, in medium effects do not affect the potential (real and imaginary). The potential is the vacuum one. This is an exact statement that follows from the effective field theory description of the system [3]. In this situation, in medium effects, which are of relative order (a 0 T ) 2 or smaller according to the power counting of the effective field theory, do not enter the Hamiltonian but affect the evolution equation through state dependent functions rather than transport coefficients. They may be included systematically in the evolution equation as done in [32,33]. By neglecting order (a 0 T ) 2 corrections, which are small at low temperatures, the quarkonium evolution for T ≤ T f is therefore that one of a Coulombic bound state in the vacuum.
The value that we have chosen for T f is somewhat arbitrary and a possible source of uncertainty in the determination of R AA . In figure 19, we show for T f = (250 ± 25) MeV the survival probabilities of the 1S, 2S, and 3S bottomonium states computed evolving from an initial Gaussian state with H eff ; H eff provides the bulk of the evolution as we have seen in the main body of the paper (cf. figure 7). The variations of the survival probabilities due to a 10% uncertainty in T f are of the same size as the ones due to the uncertainties in κ and γ.

E Finite size effects
As mentioned in section 4, the Gaussian initial condition used in the main body results in a rapidly spreading wave function when jumps are turned on. As a result, one should make sure that there are not significant effects due to the finite box size (L) used in the simulations. The evolution of the expectation value of r varies for each quantum trajectory and in runs in which there are many jumps, one may start to hit the edge of the box and the wave function may even reflect back from the edge of the box. However, such situations are dominated by octet configurations with rather large l and usually do not have time to reflect back to the left side of the box by the end of the run. As a result, such cases have essentially zero overlap with the low l compact bound states of interest, due to having very large l and being in the octet configuration most of the time.
In order to demonstrate that the results presented in the main body of the manuscript are not significantly affected by finite size effects, in figure 20 we present a comparison of survival probabilities for the 1S, 2S, and 3S states and quantum mechanical overlaps for the 1P and 2P states obtained using two different box sizes corresponding to L = 108 a 0 (open red circles) and L = 216 a 0 (open black squares) with fixed lattice spacing. In both cases, we used S-wave Gaussian initial conditions with zero impact parameter b = 0 fm, T f = 250 MeV,κ(T ) =κ C (T ), andγ = −1.75. We choose b = 0 fm for this test because this case has the longest evolution time and hence is the most susceptible to any finite size effects during the evolution. As shown in figure 20, the S-wave results obtained are independent of L within the statistical errors reported. Quantitatively, the ratio of the L = 206 a 0 results to the L = 108 a 0 , shown in figure 20 are {0.98 ± 0.04, 0.95 ± 0.06, 0.97 ± 0.06, 0.95 ± 0.07, 0.96 ± 0.06, 1.10 ± 0.10} computed from left to right, which demonstrates that our S-wave results for the two box sizes are consistent with unity within statistical errors stemming from the average over quantum trajectories. We have performed similar tests for P -wave initial conditions and found that, for both box sizes, the resulting b = 0 survival probabilities/overlaps are always extremely small.