Bottomonium suppression and elliptic flow using Heavy Quarkonium Quantum Dynamics

We introduce a framework called Heavy Quarkonium Quantum Dynamics (HQQD) which can be used to compute the dynamical suppression of heavy quarkonia propagating in the quark-gluon plasma using real-time in-medium quantum evolution. Using HQQD we compute large sets of real-time solutions to the Schrödinger equation using a realistic in-medium complex-valued potential. We sample 2 million quarkonia wave packet trajectories and evolve them through the QGP using HQQD to obtain their survival probabilities. The computation is performed using three different HQQD model parameter sets in order to estimate our systematic uncertainty. After taking into account final state feed down we compare our results to existing experimental data for the suppression and elliptic flow of bottomonium states and find that HQQD predictions are good agreement with available data for RAA as a function of Npart and pT collected at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV. In the case of v2 for the various states, we find that the path-length dependence of ϒ(1s) suppression results in quite small v2 for ϒ(1s). Our prediction for the integrated elliptic flow for ϒ(1s) in the 10−90% centrality class, which now includes an estimate of the systematic error, is v2[ϒ(1s)] = 0.003 ± 0.0007 ±0.00130.0006\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {}_{0.0013}^{0.0006} $$\end{document}. We also find that, due to their increased suppression, excited bottomonium states have a larger elliptic flow. Based on this observation we make predictions for v2[ϒ(2s)] and v2[ϒ(3s)] as a function of centrality and transverse momentum.


Introduction
It is expected that nuclear matter undergoes a phase transition to a primordial state called the quark-gluon plasma (QGP) at temperatures exceeding T c ∼ 155 MeV [1,2]. During the phase transition hadronic states breakup and the appropriate degrees of freedom for describing the system become (quasi-) quarks and gluons. The temperature at which a given hadronic bound state is expected to breakup scales with the binding energy of the state being considered. Hadronic states with low binding energies such as the π are expected to disassociate at temperature close to T c , however, states with higher binding energy such as the J/ψ and Υ have disassociation temperatures well above T c . Since they are not fully disassociated, one can use the production of heavy quarkonium bound states as a probe of the conditions generated in the QGP, such as its initial central temperature, flow profile, etc.
Experimentally, the QGP is produced using relativistic heavy-ion collision experiments at Brook-haven National Laboratory's (BNL) Relativistic Heavy Ion Collider (RHIC) and the European Organization for Nuclear Research's (CERN) Large Hadron Collider (LHC).
One key observables for LHC heavy-ion experiments is the production of heavy quarkonium states, with early results measured at the Super Proton Synchrotron (SPS) in the mid-90s already hinting at the QGP-induced suppression of the J/ψ. The expectation that one would observe suppression of heavy-quarkonium states in a thermal medium was based on the early work of Karsch, Mehr, and Satz (KMS) [3][4][5][6] who pointed out using a non-relativistic potential model that in-medium screening of the gluonic interactions would cause heavy-quark bound states to disassociate. Such a non-relativistic treatment is justified by the fact that, as the mass of the heavy-quark increases, its velocity inside the bound state decreases [7][8][9][10]. This understanding can be made more formal using effective field JHEP03(2021)235 theory methods to integrate out different energy/momentum scales, resulting in potential non-relativistic QCD (pNRQCD) [11][12][13].
One key difference from the early days of KMS is that it is now understood that the in-medium heavy-quark potential is complex-valued, with the imaginary part being related to the in-medium breakup rate of heavy-quark bound states [14][15][16][17][18][19][20][21][22][23]. The imaginary part of the potential was shown to be related to gluon dissociation or parton free dissociation of the states in refs. [24,25]. As a consequence of the imaginary part of the potential one has a non-Hermitian Hamiltonian and the time evolution of the system becomes non-unitary. This can be understood in the context of open quantum systems in which the heavy-quark bound state is coupled to a thermal heat bath [26][27][28][29][30][31][32][33][34][35]. Such imaginary contributions also appear in the context of kinetic transport models since in these calculations one also includes the possibility of in-medium breakup [36][37][38][39][40][41][42][43].
In this paper, we follow up our previous work from ref. [44] by including more details of the calculation, in particular with respect to the numerical method for solving the timedependent Schrödinger equation and our method for including the effect of late-time feed down and we update our treatment of feed down to include a 9 × 9 feed down matrix. In addition, in this paper we provide estimates of uncertainties associated with the choice of model parameters used in ref. [44]. As in ref. [44] we will focus on bottomonium states and present a model called Heavy Quarkonium Quantum Dynamics (HQQD) in which we solve the time-dependent Schrödinger equation with a complex in-medium potential for a large set of Monte-Carlo-sampled bottomonium wave-packet trajectories. For each trajectory, the quantum states are in a quantum linear superposition and one can extract the survival probability of a given state by computing the quantum-mechanical overlap of the state's vacuum eigenstate with the in-medium evolved quantum wave-function. Since we do not include an explicit time-dependent noise contribution in the potential, the resulting wavefunctions correspond to the averaged wave-function (thermal expectation value) [31]. 1 This approximation allows us to straightforwardly go beyond previous phenomenological works which made use of the adiabatic approximation instead of real-time solutions [20,[45][46][47][48][49][50][51][52][53]. In a previous paper [54], we made a preliminary investigation of the effects of relaxing the adiabatic approximation, finding that there were potentially important effects on the survival probability of the states.
Herein, we turn this approach into a more complete phenomenological framework, which can be used for comparisons with experimental data. To do this, we make use of the output of a 3+1D anisotropic hydrodynamics code which has been tuned to reproduce a large set of soft hadronic observables in √ s N N = 5.02 TeV collisions [55][56][57][58][59][60][61][62]. After computing each state's survival probability, we then take into account late-time feed down of excited states using vacuum branching ratios available from the Particle Data Group [63]. We find that our HQQD results are in quite reasonable agreement with available data given current uncertainties, however, some quantitative differences remain which motivate going beyond the methods used herein to more fully include the effects of in-medium thermal noise, initial production in octet states, and singlet-octet transitions.

JHEP03(2021)235
The structure of this paper is as follows. In section 2 we introduce the theoretical tools necessary to derive the in-medium time-dependent Schrödinger equation obeyed by heavyquark bound states. In section 3 we introduce the phenomenological potential model used in HQQD. In section 4 we provide details of the split-step pseudo-spectral method used herein. In section 5 we provide details concerning the implementation of late-time feeddown and the computation of R AA [Υ] and v 2 [Υ]. In section 6 we present our final results, providing estimates of both the statistical errors coming from sampling and systematic errors associated with varying the model parameters. In section 7 we give our conclusions and an outlook for the future.

Relation of noise and the imaginary part of the potential
To begin we review the stochastic potential model with two noise sources for a heavy quark and anti-quark described in refs. [28,31]. Let us consider a heavy quark located at x = R + r 2 and a heavy anti-quark at x = R − r 2 , and the strength of the potential between them is not only weakened by screening but also fluctuates. Therefore, the Hamiltonian consists of a screened potential V (r) and stochastic noise term Θ(r, t) Here, Θ(r, t) is the sum of noise terms for the heavy quark, θ(x, t), and anti-quark, θ(x , t), which have the following correlations From eq. (2.2), we write Using eq. (2.3) in eq. (2.4) and assuming that the noise correlation D-function is an even function i.e., D(−r) = D(r), one obtains the correlation of Θ(r, t)

JHEP03(2021)235
Next we expand the unitary time evolution operator e −i∆t H(r,t) in terms of the infinitesimal time ∆t: where in going from the second to third lines we have dropped the terms of O (∆t) 3/2 and higher. Here in eq. (2.7) we use Θ(r, t) = θ (x, t) − θ (x , t). From eq. (2.7) we extract the effective Hamiltonian as where in eq. (2.8), we use θ( ∆t , in eq. (2.9) we use the fact that x = R + r 2 , x = R − r 2 ⇒ x − x = r, and, in eq. (2.10), we have inserted eq. (2.1). Since Θ(r, t) = 0, the averaged effective Hamiltonian can be written as Therefore, the stochastic Schrödinger equation for the noise averaged quarkonium wave function is,

Phenomenological potential
In vacuum we take the heavy-quark potential to be given by a Cornell potential with a finite string breaking distance For the in-medium potential we use an internal energy based potential. The free energy of a static heavy quark-antiquark pair in an isotropic plasma [3][4][5][6]46] where a = g 2 C F /4π and σ is the string tension. We take the in-medium gluonic Debye mass to be m D = λ 4πN c (1 + N f /6)α s T 2 /3. Here we will consider the cases λ ∈ {0.8, 1, 1.2}, with the central value being the known high-temperature limit of the gluonic Debye mass. The internal energy, U , of the states which has an entropy contribution is given by where we have used S = − ∂F ∂T . Eq. (3.3) gives us the internal energy based real part of the Karsch-Mehr-Satz (KMS) potential We note that, in equilibrium, there are arguments to support the use of the internal energy [5,6,64], however, since bottomonium states are not expected to be in thermal equilibrium with their surroundings, it is unclear what the correct choice should be. In order to make sure that we have the correct equilibrium limit we choose to use the internal energy prescription. To match smoothly onto the zero temperature limit we use . (3.5) In the limit T → 0, eq. (3.5) reduces to eq. (3.1).

JHEP03(2021)235
For the imaginary part of the potential we take the result of the leading-order resummed perturbative QCD calculation of Laine et al. (3.6) with φ(r) ≡ 1 − 2 ∞ 0 sin(z)/(z 2 +r 2 ) 2 [14]. We evaluate the strong coupling α s at the scale µ = 2πT and use three-loop running [63] with Λ M S = 344 MeV. This value of Λ M S is chosen in order to reproduce the lattice result for the running coupling α s (5 GeV) = 0.2034 [65]. The resulting final complex-valued potential is of the form (3.7)

Numerical method for solving the Schrödinger equation
The complex potential is a purely radial potential hence the general solution in spherical coordinates can be written as where Y m are spherical harmonics. Making a change of variables to u m (r, t) ≡ rR m (r, t), we obtain where u(r, θ, φ, t) = rψ(r, θ, φ, t). This change of variables allows us to write the Hamiltonian in one-dimensional formĤ where V eff, (r, t) = V (r, t) + ( +1) 2mr 2 andp = −i d dr . By applying the time evolution operator to u, one obtains [54] u(r, θ, φ, t + ∆t) = exp(−iĤ∆t)u(r, θ, φ, t) where N is a normalization constant, and the summation limits are implicit. Using which tells us that each of the different states can be updated independently using the corresponding Hamiltonian. Finally, for the normalization of the various states, we take ∞ 0 dr u * (r, t)u (r, t) = 1 .

JHEP03(2021)235
We will use the time evolution operator (4.5) to evolve the wave function given a particular initial condition. To enforce the fact that the function u (r, t) must vanish at the origin, we use real-valued Fourier sine series to describe both the real and imaginary parts of the wave function. Hence, the correct boundary conditions at r = 0 are satisfied automatically. The resulting time evolution steps are [54] • Step-1: Update in configuration space using a half-step: ψ 1 = exp(−iV ∆t/2)ψ 0 .
• Step-2: Perform Fourier sine transformations (F s ) on real and imaginary parts separately: Step-3: Update in momentum space using:ψ 2 = exp −i p 2 2m ∆t ψ 1 . • Step-4: Perform inverse Fourier sine transformations (F −1 s ) on real and imaginary parts separately: • Step-5: Update in configuration space using a half-step: The discrete sine transforms (DST) above can be implemented using standard routines for Fast Fourier Transforms. By repeating this procedure, we can evolve the wave function forward in time in a manifestly unitary manner which is generally faster and more accurate than traditional finite-difference based evolution schemes. We note that besides the improved performance observed, one major benefit of the DST algorithm is that the derivatives are computed using all points in the lattice, not just at just fixed number of points. As a result, the evolution obtained using the DST algorithm is more accurate than with, for example, a Crank-Nicolson (CN) scheme using a three-point derivative [54]. In practice, we use the CUDA CUFFT library for the implementation [66]. The resulting code is massively parallelized.

Computation of R AA including feed-down
Using the complex potential eq. (3.7) along with eq. (3.5) and eq. (3.6), we then numerically solve the time-dependent Schrödinger equation on a discrete lattice. The DST method used is manifestly unitary for real-valued potentials and based on a split-step pseudospectral method described above [54,67,68]. For the discretization we use N = 4096 points in r with L = r max = 19.7 fm with a corresponding lattice spacing of a 0.0048 fm. We compute the in-medium suppression for = 0 and = 1 states, separately, following eq. (4.5).
We choose a Gaussian initial wave-function [44] u (r, with ∆ = 0.04 fm. This initial conditions mimics the local (delta function) production of bottomonium states in space resulting from their initial hard production. This type of initial state is a quantum superposition of many eigenstates of the Schrödinger equation for a fixed . Once the wave-function evolves forward in time, we obtain the survival probability of a given vacuum state by computing the overlap of the in-medium quantum

JHEP03(2021)235
wave-function with the vacuum basis states. These overlaps decay in time because of the non-Hermitian nature of the system's Hamiltonian which is fundamentally related to the in-medium breakup of bottomonium states.
Here we solve the 3+1D Schrödinger equation for a large set of trajectories using a realistic 3+1D hydrodynamics background tuned to data [60,62]. We use smooth optical Glauber initial conditions which provide the initial energy density profile as a function of the impact parameter. For the background evolution we used an initial central temperature of T 0 = 630 MeV at τ 0 = 0.25 fm/c and a constant specific shear viscosity of 4πη/s = 2, which correspond to the results obtained in a recent aHydroQP fit to soft observables at 5 TeV [62]. We evolve the quantum wave-packets using the vacuum potential starting at τ = 0 fm/c and turn on the in-medium complex potential at τ = τ med . Finally, when the temperature on a given trajectory drops below the QGP transition temperature, T QGP = 155 MeV, we again use the vacuum potential for the wave-function evolution. In the results section we will present our findings when varying τ med ∈ {0.25, 0.4, 0.6} fm/c Due to the fact that each sampled wave-packet experiences a different temperature along its trajectory while propagating through the QGP, we numerically solve the timedependent Schrödinger equation for a large set of bottomonium trajectories (2 million). 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. Here, we consider all bottomonia to have zero rapidity, y = 0, in order to simplify the calculation in a manner that takes advantage of the approximate boost-invariance of the QGP. We assume that the transverse momentum (p T )-distribution is proportional to p T /(p 2 T + M 2 ) 2 for all the states, where M is the average mass of all states being considered and, once again, we Monte-Carlo sample the p T for each particle generated. Finally we sample the initial azimuthal angle φ from a uniform distribution between 0 and 2π. Once we have the initial position, transverse momentum, and azimuthal angle, we then record the aHydroQP result for the QGP temperature along the (assumed) straight line trajectory followed by the quantum wave-packets.
Once each quantum state propagates along its trajectory, the survival probabilities are converted into particle number by multiplying by 1. The expected number of binary collisions in the centrality bin sampled and 2. The primordial (direct) production cross section for each bottomonium state.
At this point one needs to take into account the late time feed down of excited bottomonium states. This can be done by introducing a feed down matrix F which collects the known information about excited bottomonium state decays available from the Particle Data Group [63]. In the case of proton-proton collisions one can take the primordial cross sections for production and convert this into the post feed down cross sections (abundances) JHEP03(2021)235 Table 1. The experimentally observed (σ exp ) and primordial (σ primordial ) pp → bottomonium production cross sections of the various bottomonium states. Data were collected from refs. [63,[69][70][71][72][73][74][75]. The primordial (direct) cross-sections were obtained from the experimentally measured post feed down cross sections using, σ primordial = F −1 σ exp . by applying a feed down matrix to the primordial cross sections σ exp = F σ primordial with the vector σ collecting the observed and primordial (direct) cross sections for the Knowing the experimental values for the production cross-sections σ exp (listed in table 1), one can compute the primordial cross sections via σ primordial = F −1 σ exp . The resulting values for σ primordial are listed in table 1. We extract the cross sections of Υ(1s), Υ(2s), and Υ(3s) from the 5.02 TeV data obtained by the CMS collaboration in the rapidity interval |y| ≤ 2.4 [71]. The left panel of figure 3 of ref. [71] presents B × dσ/dy, where B is the dimuon branching fraction. Once we take the average over rapidity in the interval presented in that 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% [63], respectively [63], and obtains The χ b -cross sections can be obtained by through knowledge of σ[Υ(1s)], the ratios , and the ratios of the 1p and 2p cross-sections to 1s cross sections provided in ref. [72]. We assume σ[χ b2 (np)]/σ[χ b1 (np)] = 1.176 [73] for both the 1p and 2p states. This value is consistent with available experimental data [74]. To obtain the χ bj (np) cross sections, the values of R Υ(1s) are taken from the lowest p T bins of tables 5 and 6 of ref. [72] (measured at √ s = 7 TeV and √ s = 8 TeV, respectively) and extrapolated to √ s = 5 TeV. Using this extrapolation, σ[χ bj (np)] for j, n = 1, 2 can be extracted by using eq. (1) of ref. [72]. The χ b0 (np) cross sections are taken to be one quarter of the

JHEP03(2021)235
average of the σ[χ b1 (np)] and σ[χ b2 (np)] cross-sections [73]. This gives Note that herein we ignore possible p T -dependence in the branching ratios which could potentially be important for quantitative understanding of the p T -dependence of post feed down R AA [75].
To compute the effect of final-state feed-down in heavy-ion 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) × σ primordial ). We then multiply the result by the same where we have discarded terms for which the coefficients were zero in (5.2). Finally, we compute R AA for each state by dividing the final number produced after feed down by the average number of binary collisions in the sampled centrality class times the post feed-down pp production cross-section for each state (σ exp ). The momentum anisotropic flow of bottomonium states can be obtained in a similar manner. In this case, we can exploit the fact that the aHydroQP background always has the reaction plane aligned with the x and y axes, i.e. Ψ RP = 0, to simply average cos(nφ) over all particles JHEP03(2021)235 Table 2. Decay modes of Υ(2s), χ b (1p), and Υ(3s) and the calculation of different elements of the feed-down matrix (5.2) using the experimentally measured branching ratios of the various bottomonium states [63]. Table 3. Decay modes of χ b (2p) and the calculation of different elements of the feed-down matrix (5.2) using the experimentally measured branching ratios of the various bottomonium states [63]. produced after feed down is taken into account. In practice we compute the following average over all sampled trajectories in a given centrality and p T bin, v n ≡ cos(nφ) . For both R AA and v 2 , we can directly compute the statistical uncertainty associated with the sum over the sampled quantum wave-packet trajectories. To obtain an estimate of our systematic uncertainty we will vary the HQQD model parameters used in ref. [44].

Results
In order to estimate the systematic uncertainty coming from the choice of model parameters we will vary the medium initialization time (τ med ) and the pre-factor of the peturbative Debye mass (λ) in the set where the factor λ accounts for higher-order corrections to the perturbative Debye mass with λ = 1 giving the known leading-order result from hard-thermal-loop resummation [76,77]. We present the HQQD predictions for the suppression of Υ(1s), Υ(2s), and Υ(3s) states as a function of N part in figure 1 (left). For the left panel of figure 1 we applied a transverse momentum cut of p T < 30 GeV in HQQD. The HQQD results are compared with experimental data obtained by the ALICE [78], ATLAS [79], and CMS [71] collaborations, shown as circles, squares, and triangles, respectively. From figure 1 (left), we see that, for the central choice of model parameters corresponding to τ 0 = 0.4 fm/c and λ = 1, HQQD does a quite reasonable job in describing the N part dependence of R AA [Υ(1s)]. The shaded bands in figure 1 (left) show the variation of the results obtained when taking the limits indicated in eq. (6.1), with the first and third cases mapping to stronger and weaker suppression compared to the central line, respectively. In the case of the Υ(2s), one sees that using the central parameters, HQQD seems to under-predict the observed R AA [Υ(2s)].

JHEP03(2021)235
Turning next to figure 1 (right), we present the HQQD predictions for the transversemomentum dependence of R AA [Υ]. To obtain the HQQD predictions for this panel, we averaged over centrality with a weight function w(c) = exp(−c/20), with c ∈ [0, 100]. This weight function allows us to mimic the experimentally observed distribution of the number of Υ states versus centrality [81]. As can be seen from figure 1 (right), HQQD predicts a very weak dependence of R AA [Υ] on p T , with only a small decrease in R AA [Υ] at momentum less than the mass scale of the bottomonium states. One can understand the increased suppression at low-p T in terms of the average effective lifetime of the quantum wave-packets in the QGP. For low-momentum wavepackets, their effective lifetime in the QGP is longer due to their lower velocities and, hence, they experience stronger suppression than high-momentum wave-packets which can escape from the QGP more quickly.
In figure 2 we plot the HQQD prediction for the Υ(2s) to Υ(1s) double ratio as a function of N part (left panel) and p T (right panel). We compare these predictions with data reported by the ATLAS and CMS collaboration in refs. [79,80]. The band shown in the figure provides an estimate of our systematic uncertainty. As can be seen from both of the panels of this figure, we see further evidence that HQQD is predicting too much Υ(2s) suppression compared to the experimental data. It is unclear at this moment time what is the source of this discrepancy. We will return to this question in the conclusions where we discuss the possible impact of including stochastic in-medium color and angular momentum transitions in a more complete manner.
In order to gain further insight into the spatial-distribution of bottomonium suppression in the QGP one can construct "tomographic" plots showing the starting positions of the sampled bottomonium states, color-coded by the final, time-evolved R AA for that state. From such plots one can learn where in the QGP bottomonium states are most strongly suppressed and where those with only weak-suppression emerge from. A priori, one would expect trajectories that experience the strongest suppression to be those initially produced in the central of the fireball and trajectories with the weakest suppression to be those initially produced on the periphery of the collision. This expectation is indeed observed in our simulations. In figure 3 we show such a tomographic plot for the Υ(1s) (left) and the Υ(2s) (right). For both panels in this figure, we considered a Pb-Pb collision with impact parameter of b 9.2 fm. Comparing the left and right panels we firstly see the increased overall suppression of the Υ(2s) relative to the Υ(1s). In both panels, we see that states with R AA close to one come from initial production near the surface, while those experiencing the strongest suppression come from the center of the fireball as expected.
In order to better visualize the probability density for particles in each R AA -class, in figure 4 we have plotted each R AA -class in a separate panel with the panels corresponding to 0.1 intervals in R AA in the range 0.4 ≤ R AA ≤ 1. In each panel we show the normalized probability distribution function (PDF) for each R AA -class. As in figure 3 we see that the suppression of the states is correlated with the initial geometry and that states with JHEP03(2021)235  weak suppression typically come from the surface of the fireball. Statistically, one can see, for example, that there is little overlap between the PDFs for particles with R AA ≥ 0.7 and R AA ≤ 0.5. By comparing the production of particles in these two classes, one could investigate potential differences in bottomonium production in the central and peripheral regions of the QGP. For now, however, such plots are mostly useful in making sure that our understanding of suppression dynamics is correct.
In figure 5, we present the HQQD predictions for centrality dependence of the elliptic flow of Υ(1s), Υ(2s), and Υ(3s) states. In this figure, we take p T < 50 GeV and compute v 2 in 10 equally spaced centrality bins from 0-100%. The bands show the statistical uncertainty associated with the mean values extracted in each bin. The upper (squares), central (circles), and lower (triangles) lines correspond to the three different cases of varying the model parameters specified in eq. (6.1), respectively. For each of the three choices of model parameters shown in figure 5, one sees an ordering of the elliptic flow of the states with the Υ(3s) state having the largest flow and the Υ(1s) the smallest. This result is expected since excited states are more suppressed and hence experience a larger inmedium path-length dependence. For each of the states shown in figure 5 one notices that the oscillations in v 2 (to be discussed later) depend on the choice of the model parameters. This can be attributed to the change in the effective path-length of the quantum wavepackets travel through the medium when changing τ med . Finally, we note that v 2 is more sensitive to these oscillations because it explicitly involves, for example, differences between the survival probability along the short and long sides of the QGP. In R AA one sums over all angles and averages and, as a result, these oscillations are smoothed out in the sum over quantum wave-packet trajectories. On the contrary, since for v 2 different angles contribute with different weights/signs, it is naturally more sensitive to small differences in the survival probability and hence is more sensitive to these oscillations.
From figure 5, we also observe that the elliptic flow for all states goes to zero for central collisions (centrality = 0%). This is expected using the smooth (optical) Glauber model used herein for the hydrodynamic initial condition. In the future it will be interesting to assess the degree of Υ elliptic flow in the most central classes taking into account also JHEP03(2021)235 geometrical fluctuations in the initial condition (energy density, flow velocities, etc.). We note, however, that since we use smooth hydrodynamic initial conditions v 2 must go to zero for central collisions. This can be used as a check on our calculation of v 2 [Υ] for all three states. We also note that in the opposite limit, that of ultra-peripheral collisions (centrality = 100%), one expects the v 2 for all states to go to zero since the QGP lifetime goes to zero in this limit. This again provides a non-trivial test of our computation of v 2 for bottomonium states.
As alluded to above, from figure 5 one notices that the dependence of v 2 on centrality can contain oscillations (non-monotonic behavior). These oscillations arise in HQQD due to quantum oscillations in the overlaps of the in-medium bottomonium states with their corresponding vacuum states. Wave-packets which experience different effective path-lengths can exit the QGP at different points during this quantum oscillation. Since the characteristic period for these oscillations is on the order of 1 fm/c, this can have an observable effect on the HQQD predictions for v 2 as a function of centrality.
A comparison of HQQD predictions for v 2 [Υ(1s)] with experimental data collected by the ALICE [82] and CMS [83] collaborations in three different transverse momentum bins, 0-4, 4-6, and 6-15 GeV, is shown in figure 6. The HQQD predictions and experimental results are integrated over centrality in the range 5-60%. One sees from this figure that HQQD predicts results consistent with v 2 being zero, slightly negative, and small but positive value in the lowest p T bin, central p T bin, and highest p T bin, respectively. A similar trend has also been predicted earlier by an adiabatic approximation based model in ref. [53]. According to ref. [53], the negative v 2 observed in the central p T bin is related to the transverse expansion of the QGP overtaking bottomonia states which have escaped from near the surface of the QGP. To demonstrate that this was indeed the cause for the negative v 2 , the authors ref. [53] considered the case in which the transverse expansion of the QGP was turned off, finding that in this case v 2 was always positive. Returning to the theory to data comparison, figure 6 demonstrates that HQQD has a reasonable agreement with the available experimental data given current experimental uncertainties. Finally, we note that even given the experimental uncertainties, we see a similar trend in the experimental results as a function of p T as predicted by HQQD.
In figure 7, we present the centrality dependence of v 2 [Υ(1s)] in 10-30%, 30-50%, 50-90%, and 10-90% centrality bins and compare our HQQD predictions with experimental data from the CMS collaboration [83]. We show a comparison of the HQQD prediction with the experimental result integrated over 10-90% centrality in the right panel. As can be seen from this figure, in the case of the integrated v 2 [Υ(1s)] in the 10-90% bin (right panel), there is quite reasonable agreement, within experimental uncertainties, between the HQQD prediction and the data reported by the CMS collaboration. In the left panel of the figure 7, good agreement between HQQD predictions and CMS data can be seen in the 10-30% bin, however, there are larger differences in the other two centrality bins but still within 2σ of the experimental data. Hopefully, higher experimental statistics from future runs will reduce the experimental uncertainties in the near future.
In  [82] and CMS [83] collaborations. The shaded red boxes associated with each HQQD prediction indicate the systematic variation observed when changing the model parameters in the range specified in eq. (6.1). Note that, for the highest p T bin, the systematic error estimate is so small that it is indistinguishable from the horizontal error bar. currently available from the CMS collaboration and has been shown as a green triangle in the 10-90% panel (right). As can be seen from the right panel of the figure, the integrated result for v 2 [Υ(2s)] is within the experimental uncertainties, however, at the very top end of them. We again hope that more statistics will help in making more constraining comparisons in the near future. Note that for v 2 [Υ(2s)] and v 2 [Υ(3s)] one sees a larger systematic uncertainty than for v 2 [Υ(1s)].
In figure 9, we present HQQD predictions for the elliptic flow of Υ(2s) and Υ(3s) states as a function of transverse momentum in the 5-60% centrality bin. As can be seen from this figure, a sizable v 2 is predicted by HQQD for the Υ(3s) because of the strong path length dependence of R AA [Υ(3s)] between the short and long sides of the QGP fireball. Similar to the Υ(1s), for the central values of our model parameters, HQQD predicts a negative v 2 for the Υ(2s) in the lowest two p T -bins. The reason behind this negative v 2 is again that the transverse expansion of the QGP occurs more rapidly along the short side than the long side which has the affect of overtaking bottomonium states which had previously escaped the QGP with φ ∼ 0. One can see that HQQD makes a prediction that v 2 is positive in the highest p T -bin for both states. Turning to the model variation shown as red and blue box in figure 9, one sees that in the 4-6 GeV p T bin the variation of the model parameters can result in a positive v 2 [Υ(2s)]. This large systematic HQQD uncertainty is related to the quantum-mechanical oscillations discussed previously in the context of figure 5. Given these larger uncertainties, however, one still observes an ordering of the elliptic flow of the excited states with v 2 [Υ(3s)] > v 2 [Υ(2s)] in the highest two p T bins.
Finally, in table 4 we compare the predictions of HQQD for various observables with the corresponding experimental results from the ALICE, ATLAS, and CMS collaborations. In each row of the table, the results are integrated over centrality and transverse momentum in the ranges shown in the middle column, the first line of the third column is the experimental data (if available), and the second line of the third column shows the HQQD prediction.  with available experimental data. In each row, the first column indicates the observable, the second column indicates the source of the experimental result and relevant cuts, and the third column shows the experimental result on the first line and the HQQD prediction on the second line. For both HQQD predictions and the experimental results, the first uncertainty reported is statistical uncertainty and the second is systematic uncertainty. In the case of HQQD, the systematic uncertainities were estimated by varying the HQQD model parameters in the range specified in eq. (6.1).

Conclusions and outlook
In this paper, we have extended our previous work [44] by allowing for variation in the underlying HQQD model parameters in order to assess the impact such variation has on HQQD predictions for R AA [Υ] and v 2 [Υ]. In addition, herein, we provided more details concerning the numerical method used in HQQD and the method by which final state feed down effects were included. Similar to our previous paper we used real-time quantum evolution averaged over a large set of bottomonium trajectories (2 million). After averaging over all wave-packet trajectories, for each set of model parameters we were able to obtain precise estimates for R AA . For the central values of the model parameters, we found quite reasonable agreement with available experimental data for R AA [Υ(1s)], however, we found that HQQD over-predicts the suppression of Υ(2s) compared to experimental observations. For the case in which we decreased the initialization time and increased the effective Debye mass, we found that the model predicts perhaps too much Υ(1s) suppression as can be seen from the lower limits indicated in figure 1. For the opposite case in which we increased the initialization time and decreased the effective Debye mass, we found that the model
Using this same set of model parameters, we then looked at the impact of our model assumptions on v 2 [Υ]. We found that the results for v 2 [Υ(1s)] were relatively insensitive to the choice of model parameters, with the maximum v 2 [Υ(1s)] as a function of centrality being on the order of 1% in all three cases considered. We found stronger variations in the HQQD predictions for v 2 [Υ(2s)] and v 2 [Υ(3s)], however, we still observed a maximal v 2 [Υ(2s)] on the order of 3% and v 2 [Υ(2s)] on the order of 5%. For each of the sets, we found that the ordering v 2 [Υ(1s)] v 2 [Υ(2s)] v 2 [Υ(3s)] remains valid. This ordering is expected due to the fact that excited states experience more suppression and, hence, more path-length dependence of their propagation than the ground state. All of our HQQD predictions for integrated suppression and elliptic flow of bottomonium states were finally summarized in table 4 which includes estimates of both the systematic and statistical model uncertainties.
With respect to the HQQD over prediction for Υ(2s) suppression, we note that HQQD includes in-medium thermal noise only through the imaginary part of the potential included in the real-time dynamics. This allowed us to extract the survival probability associated with the thermal-noise averaged quantum wave-function, however, does not fully take into account in-medium transitions between different singlet/octet and angular momentum states. There is currently ongoing work to include true stochastic noise at the level of the real-time Schrödinger equation solution, either through explicit noisy potentials [26,27,31,[84][85][86][87][88] or through solution of the Lindblad equation for the reduced density matrix [28,30,[32][33][34]. Preliminary results obtained using the quantum trajectories method to solve the Lindblad equation indicate that more accurate inclusion of noise effects including, e.g., in-medium singlet-octet transitions, will result in only small changes in R AA [Υ(1s)] [89]. Conversely, these studies indicate that the excited states can be less suppressed when including in-medium stochastic singlet-octet transitions. This is expected due to the possibility of "quantum regeneration" in which a state transitions to a finite octet state and then back to a singlet = 0 state, either through single or multiple in-medium quantum state transitions [89]. It is possible that inclusion of the effect of in-medium quantum jumps will help to reduce the tension between HQQD predictions for R AA [Υ(2s)] and current experimental data.
In closing, we note that herein we made use of a potential whose real part was obtained from the internal energy. There were three reasons for choosing the internal energy over the free energy: (1) it was found in prior studies using the adiabatic approximation that using the free energy resulted in too much suppression for all states (see ref. [46]), (2) there is some theoretical support for this choice when the system is in equilibrium due to entropy considerations, see e.g. refs. [5,6,64], and (3) this same potential has been used in many prior works that utilized the adiabatic approximation. Another constraint on the potential was that it possesses the correct low temperature limit, i.e. that it should reduce to a realvalued Cornell potential tuned to the known spectrum of bottomonium bound states. This was done because we wanted to (a) have a more realistic phenomenological description of bottomonium excited states and (b) match smoothly onto the known vacuum eigenstates JHEP03(2021)235 which are used, in practice, to project out the final survival probabilities. Of course, it would also be nice to add some non-perturbative contributions to the imaginary part of the potential, however, on this front things are much less constrained. One possibility would be to use lattice extractions of the imaginary part of the potential similar to what was done in ref. [50]. Such investigations are planned for the future. An additional possibility is that the choice of the model potential used herein can be improved to include, for example, the effects of non-equilibrium momentum-space anisotropies in the QGP. The inclusion of momentum-space anisotropy has been shown to reduce quarkonium suppression for all states due to a reduction in the effective in-medium Debye mass [20]. We, again, leave this for future work.