Bottomonium production in heavy-ion collisions using quantum trajectories: Differential observables and momentum anisotropy

We report predictions for the suppression and elliptic flow of the $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$ as a function of centrality and transverse momentum in ultra-relativistic heavy-ion collisions. We obtain our predictions by numerically solving a Lindblad equation for the evolution of the heavy-quarkonium reduced density matrix derived using potential nonrelativistic QCD and the formalism of open quantum systems. To numerically solve the Lindblad equation, we make use of a stochastic unraveling called the quantum trajectories algorithm. This unraveling allows us to solve the Lindblad evolution equation efficiently on large lattices with no angular momentum cutoff. The resulting evolution describes the full 3D quantum and non-abelian evolution of the reduced density matrix for bottomonium states. We expand upon our previous work by treating differential observables and elliptic flow; this is made possible by a newly implemented Monte-Carlo sampling of physical trajectories. Our final results are compared to experimental data collected in $\sqrt{s_{NN}} = 5.02$ TeV Pb-Pb collisions by the ALICE, ATLAS, and CMS collaborations.


I. INTRODUCTION
Ultra-relativistic nucleus-nucleus (AA) collisions performed at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and the Large Hadron Collider (LHC) at the European Organization for Nuclear Research (CERN) have provided unprecedented insight into the behavior of matter at extreme energy and baryon number densities the likes of which previously only existed in the very early Universe [1,2]. The goal of these experiments is to produce and study a color-ionized, or deconfined, quark-gluon plasma (QGP), a state of matter in which the degrees of freedom are quarks and gluons rather than the hadronic degrees of freedom observed at low energies in which quarks and gluons are confined. In order to determine the properties of the QGP, experimentalists at RHIC and LHC measure a variety of observables in AA, pA, and pp collisions including the spectra of produced hadrons, their azimuthal momentum correlations, photon production, dilepton production, etc. * nora.brambilla@ph.tum.de † miguelangel.escobedo@usc.es ‡ mstrick6@kent.edu § antonio.vairo@tum.de ¶ vandergriend@tum.de * * johannes.weber@physik.hu-berlin.de An observable of particular interest is the ratio of the number of heavy quarkonia observed in an AA collision to the number observed in a pp collision (scaled by the number of binary collisions); this defines the nuclear modification factor R AA of the particular quarkonium species. It was predicted decades ago that due to Debye screening at distances larger than approximately the inverse of the Debye mass, the inter-quark potential of heavy quarkonium in a color-ionized QGP becomes more short range, and consequently, the measured rates of heavy quarkonium bound state production in AA collisions would be suppressed relative to the rates in pp collisions in which no QGP is generated [3,4]. Since these early papers, there has been considerable progress in understanding the dynamics of heavy quarkonia in the QGP. A paradigmatic shift in our theoretical understanding of heavy quarkonium suppression occurred in 2007 with the findings of thermal corrections to the real part of the in-medium potential related to screening and a nonzero imaginary part related to the in-medium dissociation rate due to Landau damping [5]. Subsequent works extended this to include the effect of non-abelian singlet-octet transitions using the effective field theory (EFT) potential non-relativistic QCD (pNRQCD) [6][7][8][9]. In the interim, the existence of a large in-medium decay width has been taken into account in phenomenological calculations of R AA which use complex potential models [10][11][12][13][14][15][16]. Nonrelativistic EFTs, arXiv:2107.06222v1 [hep-ph] 13 Jul 2021 and especially pNRQCD, allow for a systematic and nonperturbative exploitation of the separation of scales inherent in heavy quark bound states.
In order to fully understand the dynamics of inmedium heavy quarkonium, a careful consideration of in-medium scattering including both dissociation and recombination is necessary. The formalism of open quantum systems (OQS) allows for a rigorous treatment of a quantum system (here the heavy quarkonium) coupled to an external environment (here the QGP) and thus provides a useful framework for treating heavy quarkonia in medium [17][18][19][20]. In the present work, we utilize a set of evolution equations describing the in-medium evolution of heavy quarkonium realizing the hierarchy of scales 1/a 0 πT ∼ m D E where a 0 is the Bohr radius of the bound state, T is the medium temperature, m D ∼ gT is the Debye screening mass, and E is the binding energy of the bound state. In this regime, the evolution equations take the form of a Lindblad equation describing the Markovian quantum Brownian motion of a heavy quarkonium in the QGP [21][22][23].
In this work, we extend Ref. [24] wherein the Lindblad equation was solved numerically using the quantum trajectories algorithm which represents a quantum unraveling of the Lindblad equation. The numerical code, developed for and presented in Ref. [24], is called QTraj and was used to make phenomenological predictions for the nuclear suppression of Υ(1S), Υ(2S), and Υ(3S) states in 5.02 TeV Pb-Pb collisions as a function of the number of participating nucleons N part . The quantum trajectories algorithm requires averaging over a set of stochasticallygenerated quantum evolutions. Due to the associated computational costs, in Ref. [24], the temperature evolution of the plasma was simplified by using an average temperature profile per centrality class computed from the average of Monte-Carlo sampled physical trajectories in that centrality class. In this work, we compute the QGP survival probability for each physical trajectory and bin the results as is done experimentally. This has been made possible by efficiency and scalability improvements to the QTraj code [25]. As a result of these improvements, we are able to present predictions for R AA and associated double ratios as functions of both N part and p T . In addition, due to the large number of physical trajectories now considered, we are able to make statistically significant predictions for the elliptic flow v 2 of the Υ(1S), Υ(2S), and Υ(3S) states as functions of both N part and p T . We compare our results to experimental data collected by the ALICE, ATLAS, and CMS collaborations.
The structure of this work is as follows: in Sec. II, we review the derivation of the Lindblad equation describing in-medium heavy-quarkonium dynamics in a stronglycoupled QGP and the quantum trajectories algorithm as implemented in the QTraj code; in Sec. III, we present our numerical results and compare to experimental data; in Sec. IV, we present our conclusions and an outlook for the future; in App. A, we present a table of QTraj pre-dictions for the centrality-integrated R AA of the Υ(1S), Υ(2S), and Υ(3S); finally, in App. B, we investigate the role of quantum jumps in heavy-quarkonium dynamics and their effect on experimental observables.

II. METHODOLOGY
A. Heavy quarkonium dynamics in a strongly-coupled quark-gluon plasma In this paper, we solve the Lindblad equation describing the in-medium dynamics of a heavy quarkonium that was derived using the EFT pNRQCD and the OQS formalism in Refs. [21][22][23]. The nonrelativistic nature of heavy-heavy bound states, i.e., v 1 where v is the quark-antiquark relative velocity, leads to at least three hierarchically ordered scales: the hard scale M of the heavy quark mass, the soft scale M v of typical momentum transfers, and the ultrasoft scale M v 2 associated with the binding energy E. If the bound state is Coulombic then v ∼ α s . Integrating out the hard scale M from full QCD gives rise to the EFT nonrelativistic QCD (NRQCD) [26,27]; further integrating out the soft scale M v gives rise to pNRQCD [28][29][30]. In this treatment, the small radius r of the lowest lying bound states allows for a multipole expansion in r. pNRQCD implements this expansion in the bound state radius r and in the inverse of the heavy quark mass M at the Lagrangian level and is thus ideally suited for describing low lying bottomonium states of small radius. The degrees of freedom in the resulting effective Lagrangian are composite fields made of heavy quark and heavy antiquark pairs in a color singlet or color octet configuration, and light quarks and gluons at the ultrasoft scale. Transitions between the singlet and octet fields are encoded in chromoelectric-dipole interaction terms.
The OQS formalism allows for the rigorous treatment of a quantum system coupled to an external environment (see Ref. [31] for a general introduction). The relevant time scales of the full system are a time scale τ S characterizing the system, a time scale τ E characterizing the environment, and a relaxation time τ R characterizing the interaction between the system and the environment. The scale τ S is set by the characteristic time scale of internal transitions in the system and, as such, is related to the inverse of the internal level spacing of states. The scale τ E is set by the time scale of equilibration of the environment, and the scale τ R is the characteristic time scale associated with the in-medium evolution of the reduced density matrix. Hierarchical orderings of these scales allow for simplifications of calculations and the realization of different evolution paradigms. For the system treated in this work, i.e., a bottomonium in a QGP at temperatures reached in current heavy ion collision experiments, one has which allows for the Markovian approximation, i.e., the system is insensitive to its prior evolution. Furthermore, one has which qualifies the evolution as quantum Brownian motion.
We consider a strongly coupled plasma in which the heavy-quark mass M , the Bohr radius of the quarkonium a 0 , the temperature of the medium T , the Debye mass m D ∼ gT , and the binding energy of the quarkonium E fulfill the hierarchy of scales In this regime, the system, the environment, and the relaxation time scales are given by where Σ s is the thermal self-energy of the system. The hierarchy of scales in Eq. (3) ensures that the evolution of the reduced density matrix is Markovian and exhibits quantum Brownian motion. Using pNRQCD and OQS and working in the regime specified in Eq. (3), in Refs. [21,22] a set of master equations governing the in-medium evolution of a heavy quarkonium was derived. In the limit T E, an expansion in E/T may be performed; at leading order, the evolution equations take the form of a Lindblad equation [32,33] where The singlet and octet density matrices ρ s (t) and ρ o (t) describe quarkonium in the singlet and octet configurations, respectively. The operators h s,o = p 2 /M + V s,o are the singlet and octet Hamiltonians with V s = −4α s (1/a 0 )/(3r) and V o = α s (1/a 0 )/(6r); α s (1/a 0 ) is the strong coupling at the energy scale of the inverse of the Bohr radius. Interactions with the strongly-coupled medium are encoded in the non-perturbative transport coefficients κ and γ with Ω(t) being a temporal Wilson line running from time negative infinity to time t, i.e., κ is the heavy quark momentum diffusion coefficient [34,35], and γ is its dispersive counterpart. As noted in Ref. [22], κ and γ are related to the thermal width Γ and mass shift δM of the bottomonium, respectively, and can, therefore, be extracted indirectly from unquenched lattice measurements of these quantities as done in Ref. [23]. More recently, direct quenched lattice measurements of κ have been performed across an unprecedentedly large range of temperatures allowing to detect the dependence of κ on the medium temperature [36]. Direct lattice extractions of γ (as opposed to the indirect extractions via δM in Ref. [23]) are currently in progress.

B. Quantum trajectories algorithm
Directly solving the Lindblad equation given in Sec. II A is computationally demanding, and previous works relied on simplifying assumptions. Specifically, Ref. [22] expanded the density matrix in spherical harmonics and introduced a cutoff at = 1, thus only considering S-and P -wave states. In Ref. [24], the quantum trajectories algorithm was utilized to solve the Lindblad equation via a computationally less intensive Monte-Carlo method. This allowed for solving of the evolution equations to all orders in while also dramatically increasing the spatial extent of the lattice and decreasing the lattice spacing compared to Ref. [22].
The quantum trajectories algorithm implements a stochastic evolution of each quantum trajectory in order to solve the Lindblad equation (frequently referred to as an unraveling of the Lindblad equation). 1 The central idea of the algorithm is to split the full evolution specified by the Lindblad equation into a diagonal contribution that leaves the quantum numbers of the system unchanged and an off-diagonal contribution that changes the quantum numbers. For this purpose, we rewrite the Lindblad equation as where The non-unitary effective Hamiltonian H eff is diagonal; its action on ρ(t) leaves the color and angular momentum state of ρ(t) unchanged but decreases its trace. The jump operators C n entering into the summation in Eq. (16) are off-diagonal and their action on ρ(t) results in a change of quantum numbers. 2 The diagonal contributions include the effect of the thermal width Γ = n C † n C n in the evolution (cf. Eq. (2.2) of Ref. [24]), and the off-diagonal terms can be mapped to quantum jumps between different states. Both of these contributions can be implemented at the level of one-dimensional wave functions rather than density matrices, thereby greatly reducing both the memory needed for the simulation and the number of computational cycles required. 3 The QTraj code implements the quantum trajectories algorithm as follows: 1. Initialize a wave function |ψ(t 0 ) at initial time t 0 which corresponds to the initial quantum state of the particle given by 2. Generate a random number 0 < r 1 < 1 and evolve the wave function forward in time with Denote the first time step fulfilling the inequality of Eq. (18) as the jump time t j . If the jump time is greater than the simulation run time t f , end the simulation at time t f ; otherwise, proceed to step 3.
3. At time t j , initiate a quantum jump: (a) If the system is in a singlet configuration, jump to octet. If the system is in an octet configuration, generate a random number 0 < r 2 < 1 and jump to singlet if r 2 < 2/7; otherwise, remain in the octet configuration.
2 This is clearly the case for C 0 as it is off diagonal in color space, i.e., it induces a singlet-octet transition (and a change of ±1 in ). C 1 is diagonal in color space but off diagonal in angular momentum space, i.e., it induces an octet-octet transition between states of angular momentum and ± 1. This can be made manifest by expanding in spherical harmonics; cf. Eqs. (83) and (84) of [22]. 3 Details concerning the QTraj implementation, including scaling studies, benchmarks, and runtime comparisons to other methods can be found in Ref. [25]. This reference accompanies the opensource release of QTraj.
(c) Multiply the wavefunction by r and normalize.

Continue from step 2.
The procedure for the calculation of the jump time t j in step 2 is known as the waiting time approach and reduces the number of random numbers to be generated compared to the standard quantum trajectories approach (see Sec. III.D of Ref. [37] and references therein). The probabilities in step 3 correspond to the branching fractions into a state of different angular momentum and/or color and are calculated via the relation Each evolution of the wave function from time t 0 to t f is called a quantum trajectory. In practice, a large number of quantum trajectories must be generated and averaged over, and, as the number of trajectories considered increases, the average converges to the solution of the Lindblad equation. This equivalence can be explicitly proven by writing |ψ(t + δt) as a superposition of a jumped state and a state evolved with H eff . For details of this proof, see Sec. III.A of Ref. [37].

C. Simulation details
In order to solve Eq. (7), we must specify the values of the transport coefficients κ and γ. For the former, we make use of recent quenched lattice measurements of κ carried out in Ref. [36] which provide κ(T ) over a large range of temperatures. All results reported in this work are carried out using three temperature-dependent parameterizations ofκ(T ) = κ(T )/T 3 which are given by the lower, central, and upper bounds of the "fit" curve of Fig. 13 of [36]. We denote these three parameterizationsκ L (T ),κ C (T ), andκ U (T ), respectively. For γ, we perform simulations with three temperature-independent values ofγ = γ/T 3 = {−3.5, −1.75, 0}. These values are taken from the relation δM (1S) = (3/2)a 2 0 γ where δM (1S) is the in-medium mass shift of the Υ(1S) state as detailed in Ref. [23]. We note that the lattice studies of Refs. [38,39] used in Ref. [23] favor larger absolute values of (the negative parameter)γ, while more recent lattice studies [40,41] favor δM (Υ(1S)) 0 and thuŝ γ 0.
For the mass, we take M = m b = m Υ(1S) /2 = 4.73 GeV with m Υ(1S) from [42]. 4 The strong coupling α s is calculated by solving where α s is evaluated at the inverse of the Bohr radius using the 1-loop running with N f = 3 flavors, and Λ [43]. The resulting value of the strong coupling constant is α s = 0.468.
For the initial state radial wave-function we use a Gaussian-smeared delta function multiplied by a power of r appropriate for the initial angular momentum state , i.e., with r ψ (t 0 ) normalized to one when summed over the entire (one-dimensional) lattice volume. Narrower initial states (smaller c) require a significantly larger number of trajectories to obtain similar statistical errors. We take the width of the Gaussian to be c = 0.2 to balance accuracy and computational effort; while this choice may cause relative systematic uncertainties of about 10% or 15% for the excited S-wave states, the S-wave ground state is unaffected (below 5% level) by changes of c within a factor of two [25]. We employ a radial lattice of NUM = 4096 lattice sites and a radial volume of L = 80 GeV −1 , corresponding to a radial lattice spacing of a ≈ 0.0195 GeV −1 . Systematic errors due to the finite lattice spacing or volume are of the same order as those due to the smeared initial state; the former is more significant for the ground state, the latter for the excited states. The real time integration employed for deterministic evolution between jumps is discretized with a time step of dt = 0.001 GeV −1 ; this time discretization leads to a quantitatively similar level of systematic errors as the other sources [25].
We expand upon our work reported in Ref. [24] by Monte-Carlo generating independent physical trajectories through the quark-gluon plasma rather than using a single path-averaged temperature evolution in each centrality bin. In Ref. [24], in each centrality bin, a pathaveraged temperature evolution was computed from the average of approximately 132000 Monte-Carlo generated physical trajectories and used to compute the survival probability. In the present work, due to increased code efficiency/scalability and access to large-scale computational resources, we sample approximately 7 -9 ×10 5 independent physical trajectories for each choice ofκ(T ) andγ, with approximately 50-100 quantum trajectories per physical trajectory. To generate each physical trajectory, we sample the bottomonium production point in the transverse plane using the nuclear binary collision overlap profile N bin AA (x, y, b), the initial transverse momentum of the state p T from an E −4 T spectrum, and the initial azimuthal angle φ of the state's momentum uniformly in [0, 2π). We bin the results for the survival probability as a function of centrality, p T , and φ. This allows us to make predictions for differential observables such as R AA as a function of p T and elliptic flow.
We use the same medium evolution as Ref. [24] that is modeled using a 3+1D dissipative relativistic hydrodynamics code, which makes use of the quasiparticle anisotropic hydrodynamics (aHydroQP) framework [44][45][46]. The code uses a realistic equation of state fit to lattice QCD measurements [47] and was tuned to soft hadronic data collected in 5.02 TeV collisions using smooth optical Glauber initial conditions in Ref. [48]. The resulting hydrodynamic parameters provide an excellent description of the experimentally observed hadronic spectra/multiplicities, extracted femtoscopic radii, and identified hadron elliptic flow with an initial central temperature of T 0 = 630 MeV at τ 0 = 0.25 fm/c and a constant specific shear viscosity of η/s = 0.159. The anisotropic hydrodynamics framework allows for an accurate description of both the early-time evolution of the quark-gluon plasma and the evolution near the transverse edges of the plasma where deviations from equilibrium are large. This is due to an all orders resummation in the inverse Reynolds number [49]. As a result, aHydroQP reliably describes even the very early stages of the collision, when non-equilibrium corrections are large, in addition to extreme cases of the flow profile, such as Gubser flow where non-equilibrium corrections are large both at early and late times [50][51][52][53][54][55][56][57][58][59][60][61][62].
In our simulations, the wave-function is initialized at time τ = 0 fm/c and evolved in the vacuum until the interaction with the medium is initialized at τ = 0.6 fm/c. To ensure that the hierarchy of scales of Eq. (3) is fulfilled and our evolution equations are valid, we evolve the state in the vacuum when the temperature falls below T f = 250 MeV. In this temperature region, the hierarchy of scales given in Eq. (3) is no longer fulfilled as πT is no longer significantly greater than the binding energy E. Hence, in this temperature region, the medium effects are ignored, and the quantum state is evolved using the vacuum potential. As this particular value of T f is somewhat arbitrary, in Ref. [24], a set of simulations were performed varying T f by ±25 MeV; the uncertainty from this variation was found to be similar in magnitude to that obtained from variation ofκ(T ) andγ. We note that the most recent lattice quantum chromodynamics (LQCD) calculations find that the pseudocritical temperature for the QGP phase transition is approximately T pc 158 MeV [63,64]. 5 All results reported in this work are obtained using T f = 250 MeV.  QTraj -  [66], ATLAS [67], and CMS [68] collaborations. The bands represent theoretical uncertainties as in Fig. 1.

D. Feed down
The QTraj code allows for a computationally efficient solution of the Lindblad equation describing the in-medium evolution of bottomonium states in the QGP. From this evolution, one can extract the survival probability of a state that has traversed the QGP. However, in order to compare to experimental measurements of the nuclear modification factor R AA , one must take into account the probability that an excited bottomonium state emerging from the plasma decays to a lower-lying bottomonium state in the vacuum before being experimentally detected. At the level of the cross section, the experimentally observed and direct production cross sections are related by σ exp = F σ direct where each entry of the σ vectors corresponds to a particular bottomonium state, and F is a matrix related to the branching ratios of the excited states. We consider the states {Υ(1S), Υ(2S), χ b0 (1P ), χ b1 (1P ), χ b2 (1P ), Υ(3S), χ b0 (2P ), χ b1 (2P ), χ b2 (2P )}. The entry F ij i < j is the branching ratio of state j to state i, F ii = 1, and F ij = 0 for i > j. The explicit values of F ij are taken from the Particle Data Group [69] and presented in Eq. (6.4) of Ref. [24].
The resulting nuclear suppression R AA of each state is computed using where S(c, p T , φ) is a diagonal matrix which collects the survival probabilities extracted from the QTraj evolution; c labels the centrality class, p T the transverse momentum, and φ the azimuthal angle.  [68,70] as explained in Sec. 6.4 of Ref. [24].

III. RESULTS
In this section, we present our final results for the nuclear modification factor R AA and the elliptic flow v 2 of the Υ(1S), Υ(2S), and Υ(3S). The theoretical uncertainties, which are indicated as shaded bands, come from varying the values of the parametersκ(T ) andγ as detailed in Sec. II C, while statistical errors are indicated by narrow bands around the individual lines which are, in many cases, smaller than the respective line widths. In App. B, for a subset of observables, we present comparisons between QTraj simulations run with the full evolution including jumps as detailed in Sec. II B and results obtained by evolving the wave function using only the effective Hamiltonian H eff without applying the jump operators. The full QTraj results presented in this section are obtained from approximately 50-100 quantum trajectories per 7 -9 ×10 5 physical trajectories for each combination ofκ(T ) andγ. The H eff results were obtained by sampling approximately 10 6 physical trajectories for each combination ofκ(T ) andγ. We compare our results with experimental data collected by the ALICE [66,72], ATLAS [67], and CMS [68,71,73]

collaborations.
A. Nuclear modification factor RAA In Fig. 1, we plot the results of our QTraj simulations for the nuclear modification factor R AA of the Υ(1S), Υ(2S), and Υ(3S) as a function of the number of participating nucleons N part . In the left panel of Fig. 1, the shaded bands indicate the variation in our QTraj results for R AA when varyingκ while holdingγ fixed at its central value; the dashed lines correspond to the lower boundκ(T ) =κ L (T ); and the dot-dashed lines correspond to the upper boundκ(T ) =κ U (T ). In the right panel of Fig. 1, the shaded bands indicate the variation in our QTraj results for R AA when varyingγ while holdinĝ κ fixed at its central value; the dashed lines correspond to the lower boundγ = −3.5; and the dot-dashed lines correspond to the upper boundγ = 0. As can be seen from this figure, the central values of these two parameters provide a good description of the N part dependence of R AA for all three states considered. Comparing the left and right panels of Fig. 1, one sees that the uncertainty associated with the variation ofγ (right panel) is larger than the one associated with the variation ofκ (left panel).
In Fig. 2, we present our results for R AA [1S], R AA [2S], and R AA [3S] as a function of transverse momentum p T . The bands, line styles, and panels represent the same variation as in Fig. 1. We observe that, within uncertainties, our results are in agreement with the experimental data. In fact, the dependence of R AA on p T is very mild. This behavior is seen both in our results and in experimental measurements. Our results for Υ(1S) show a greater sensitivity to variation ofγ than toκ; however, the opposite is true for the excited states.
In Figs. 3 and 4, we present our results for the double ratio of R AA [2S] and R AA [3S], respectively, to R AA [1S] as a function of N part . As in Fig. 1, the left and right panels correspond to the variation overκ(T ) andγ, and the line styles for the bounds are the same. We note that the data from the CMS collaboration in Fig. 4 give only an upper bound on R AA [3S]. We observe good agreement between our QTraj results and the experimentally measured values of the double ratios across the entire range of N part . Our results show a much larger dependency on γ than onκ. This suggests that this measurement can potentially constrain the value ofγ, for which there are much less lattice QCD data than forκ. Unfortunately, at the moment, the experimental uncertainties are of the order of the effect of theγ variation.
In Fig. 5, we plot the double ratio of R AA [2S] to R AA [1S] as a function of p T . The notation and parameter variation are the same as in the previous plots. What we observe in this figure confirms what we saw in previous plots. The dependence of this double ratio with p T is very mild. And similarly to what we observed in the double ratio versus the number of participants, varyinĝ κ has almost no influence while varyingγ is significant. Regarding the comparison with experimental data, we see a reasonable agreement within reported uncertainties with some tension with the data seen at large p T .

B. Elliptic flow v2
In Fig. 6, we plot our results for the elliptic flow v 2 of the Υ(1S) as a function of centrality. Again, the notation is as in previous plots. Our results agree to within uncertainties with the experimental results of the CMS collaboration, though we note the large uncertainities of the experimental results. In this case, we see that the influence ofκ andγ is similar. It is noteworthy that the more inclusive prediction (in the 10 to 90% percent centrality window) is very precise and close to the central value of the experimental results (more details below).
In Fig. 7, we plot v 2 [Υ(1S)] as a function of p T . We observe agreement to within uncertainties with the experimental results of the ALICE and CMS collaborations. In this case, the sensitivity of our results toκ andγ is similar, except for the lower momentum region, in which the sensitivity toγ is larger.
In Fig. 8, we plot our results for the elliptic flow of the Υ(2S) and Υ(3S) as a function of centrality. Our results  [67] and CMS [71] collaborations. The bands and bars represent uncertainties as in Fig. 3; we note that the CMS measurements give only an upper bound at 95% confidence level.
agree to within uncertainties with the experimental data point from the CMS collaboration, although the experimental uncertainties are at least an order of magnitude larger than our theoretical uncertainty. It is interesting to see that our model predicts very similar v 2 for both Υ(2S) and Υ(3S); v 2 for the excited states appears to be somewhat larger than for the ground state.
In the case of v 2 of the Υ(1S), we predict that it has a maximum on the order of 1.5% as a function of both centrality and transverse momentum. Our prediction for the 10-90% centrality-and p T -integrated Υ elliptic flow is v 2 [Υ(1S)] = 0.008 ± 0.003 ± 0.002, v 2 [Υ(2S)] = 0.016 ± 0.003 ± 0.002, and v 2 [Υ(3S)] = 0.015 ± 0.002 ± 0.001, where the first uncertainty corresponds to bothκ and γ variation and the second uncertainty corresponds to the statistical uncertainty due to the average over physical and quantum trajectories. We find that the 2S and 3S states have similar integrated elliptic flow, which is roughly a factor of two larger than the 1S state, 2 for 10-90%. When considering the 2S to 1S v 2 -ratio in different centrality bins, we find that, taking into account the variation over botĥ

IV. CONCLUSIONS AND OUTLOOK
In this paper, we presented a comprehensive set of predictions for the suppression and elliptic flow of Υ(1S), Υ(2S), and Υ(3S) in 5 TeV Pb-Pb collisions and compared our predictions to experimental data from the AL-ICE, ATLAS, and CMS experiments. To make our predictions, we numerically solved the 3D non-abelian Lindblad equation for the quarkonium reduced density matrix that emerges when OQS methods are applied within the pNRQCD effective field theory for a strongly-coupled QGP. The numerical solution was realized by mapping the solution of the Lindblad equation to a 1D Schrödinger equation with a non-Hermitian Hamiltonian that is subject to stochastic quantum jumps. Using the resulting quantum trajectories algorithm, we were able to simu-late the full 3D evolution of the wave-function, including the possibility of internal transitions between different color and angular momentum states.
To describe the interaction with the hot and threedimensionally expanding QGP, we made use of a realistic dissipative hydrodynamics simulation called anisotropic hydrodynamics. The initial conditions and transport coefficients used in the 3+1D aHydro code were tuned to reproduce soft observables such as identified pion, proton, and kaon p T -spectra, multiplicities, and elliptic flow. To compute R AA , we produced a large ensemble of physical quarkonium trajectories by Monte-Carlo sampling both the initial production points and transverse momentum vectors. We then computed the survival probability along each of these physical trajectories by averaging over ensembles of stochastically generated quantum trajectories. Based on the Monte-Carlo sampling of physical   trajectories, we could compute both the N part -and p Tdependence of R AA and the elliptic flow of the states. This extends our prior work where, due to the high computational demand of solving the Lindblad equation, we used a trajectory-averaged temperature evolution in each centrality bin [24]. Our final predictions also include the effect of late-time feed down of bottomonium states, the calculation of which is based on known experimental measurements of bottomonium production cross-sections and branching ratios in pp collisions. We find that the primary effect of computing the survival probability on a trajectory-by-trajectory basis is to increase both Υ(2S) and Υ(3S) R AA , which helps to bring our predictions for both R AA of these states, and the corresponding double ratios, into better agreement with available experimental data than the trajectory-averaged results presented in [24]. Associated with this paper, the QTraj code used to generate the results will be released under a public GPL license. We present the details of the code, along with examples, and benchmarks in a separate work with a more computational focus [25]. Due to the stochastic quantum trajectories algorithm and Monte-Carlo sampling of the physical trajectories, the results of our simulation had an associated statistical uncertainty. For each parameter set considered, the statistical uncertainty computed was reported in each figure based on an ensemble size of approximately 10 5 −10 6 physical trajectories. With these large ensemble sizes, the statistical uncertainty in the determination of R AA was on the order of the line width in the plots, while there remained somewhat larger statistical uncertainties in our predictions for v 2 . We estimated our theoretical uncertainties by varying the relevant transport coefficientsκ andγ in the range indicated by lattice measurements of these quantities. We found that R AA [Υ(1S)], the 2S to 1S double ratio, and 3S to 1S double ratio had a larger variation withγ than withκ, with the doubleratios rather strongly depending onγ but notκ. This observation offers some hope that, with increased statistics for both 1S and 2S R AA , one can constrainκ andγ based on experimental data.
In the case of the elliptic flow, we found similar variation in our predictions under variation ofκ and γ. We found reasonable agreement between our predictions for v 2 [Υ(1S)] and available experimental data and made predictions for the elliptic flow of the Υ(2S) and Υ(3S), finding that the differential suppression of these states results in a larger elliptic flow, as can be expected from the fact that their survival probabilities are smaller (stronger medium interactions). When considering the centrality dependence of the ratio of the elliptic flow of the 2S and 1S states, our approach predicts 2 v 2 [Υ(2S)]/v 2 [Υ(1S)] 4, with the maximum occurring in the 30-50% centrality bin. This prediction can hopefully soon be tested by experimentalists.
Turning to the future, one limitation of the framework used herein is that it relied on an assumed strict ordering of the binding energy and temperature, namely T E. As a result, at low-temperatures, the framework used herein becomes potentially unreliable. For this reason, we used a lower temperature of T f = 250 MeV for bottomonium interactions with the medium. In our previous work, it was shown that the variation of R AA when varying T f by 10% was on the same order as the theoretical uncertainty associated with the variation of the fundamental transport coefficientsκ andγ. That said, it seems necessary to include sub-leading corrections in E/T in order to gauge their impact on in-medium bottomonium dynamics [19]. Another interesting prospect is that, at low temperatures, one could interface QTraj output to codes based on a semi-classical approach in which one instead solves in-medium Boltzmann equations, see e.g. [74][75][76].   In this appendix, we present comparisons between the full Lindblad evolution including the effects of quantum jumps and evolution in which we only evolve the system with the complex Hamiltonian H eff . This will help us to assess the role played by quantum jumps and their final effect on experimental observables.
In Fig. 9, we plot a comparison of the QTraj results for R AA [2S] as a function of p T implementing the full evolution with jumps to those obtained using only H eff . Each panel presents results obtained using different values of κ(T ) andγ. We observe agreement to within uncertainties with the experimental data for all values ofκ(T ) and γ and between the full and H eff evolution for all values exceptγ = 0 (lower left panel). We note that the difference between the full Lindblad evolution and the H eff evolution is much smaller than the uncertainty obtained by varyingκ andγ. Therefore, until more precise determinations of κ and γ are available, the error made by ignoring jumps when computing R AA is negligible.  [67], and CMS [71] collaborations. We compare the results obtained using the full QTraj algorithm (blue) with results obtained using evolution with H eff with no jumps (orange). The bands represent uncertainties as in Fig. 1.
In Fig. 10, we present a comparison of results obtained using full evolution with jumps against results obtained using only H eff evolution for the double ratio R AA [Υ(2S)]/R AA [Υ(1S)] as a function of p T . As in the case of R AA [Υ(2S)], we observe the largest effect of the jumps in the caseγ = 0 andκ =κ C (lower left panel in Fig. 9) and in the caseγ = −1.75 andκ =κ U (upper right panel in Fig. 9) with agreement to within the re-    [73]. We compare the results obtained using the full QTraj algorithm (blue) with results obtained using evolution with H eff with no jumps (orange). The parameter variation is as in Fig. 9.
ported statistical uncertainties for the other values. We note that the error induced by ignoring the jumps is of the order of the uncertainty obtained by varying κ but much smaller than uncertainty obtained by varying γ. In summary, the uncertainty on the prediction of the double ratio R AA [Υ(2S)]/R AA [Υ(1S)] as a function of p T is driven by γ and a precise value of this quantity can potentially constrain the transport coefficient.
In Fig. 11, we plot a comparison of full and H eff evolution results for v 2 [Υ(1S)] as a function of p T ; the panels correspond to separate variation ofκ(T ) (left panel) orγ (right panel), while the other parameter is kept fixed. We observe again agreement to within uncertainties with the available experimental data. It is interesting to note that v 2 seems to be the only observable, within our obtained accuracy, in which the effect of the jumps competes with the uncertainties associated with the variation of κ and γ. Therefore, v 2 appears to be the observable most sensitive to quantum jumps and might provide, in the future, an observable that cannot be explained with purely H eff evolution.