Bottomonium suppression and elliptic flow from real-time quantum evolution

We compute the suppression and elliptic flow of bottomonium using real-time solutions to the Schr\"{o}dinger equation with a realistic in-medium complex-valued potential. To model the initial production, we assume that, in the limit of heavy quark masses, the wave-function can be described by a lattice-smeared (Gaussian) Dirac delta wave-function. The resulting final-state quantum-mechanical overlaps provide the survival probability of all bottomonium eigenstates. Our results are in good agreement with available data for $R_{AA}$ as a function of $N_{\rm part}$ and $p_T$ collected at $\sqrt{s_{\rm NN}} =$ 5.02 TeV. In the case of $v_2$ for the various states, we find that the path-length dependence of $\Upsilon(1s)$ suppression results in quite small $v_2$ for $\Upsilon(1s)$. Our prediction for the integrated elliptic flow for $\Upsilon(1s)$ in the $10{-}90$% centrality class is $v_2[\Upsilon(1s)] = 0.0026 \pm 0.0007$. We additionally find that, due to their increased suppression, excited bottomonium states have a larger elliptic flow and we make predictions for $v_2[\Upsilon(2s)]$ and $v_2[\Upsilon(3s)]$ as a function of centrality and transverse momentum. Similar to prior studies, we find that it is possible for bottomonium states to have negative $v_2$ at low transverse momentum.


Introduction
Relativistic heavy-ion collision experiments at Brookhaven National Laboratory's (BNL) Relativistic Heavy Ion Collider (RHIC) and the European Organization for Nuclear Research's (CERN) Large Hadron Collider (LHC) study the behavior of matter subject to extreme conditions. The goal of these experiments is to create and study the properties of the Quark-Gluon Plasma (QGP) which is expected to be produced when the energy density of matter exceeds approximately 1 GeV/fm 3 . Detailed lattice studies have demonstrated that quantum chromodynamics (QCD) has a pseudo-critical temperature of approximately T pc 155 MeV [1,2]. Since the QGP is a color-ionized phase of matter, it is expected to strongly affect the propagation of both individual partons and hadronic bound states.
Hadrons composed of light quarks are expected to disassociate at temperatures around, or just above, T pc . For heavy-quarkonium bound states, such as the J/ψ and Υ, however, it was predicted in the late 1980s that such states could survive into the QGP phase due to their large binding energy [3][4][5][6]. QCD-based model and lattice gauge theory calculations have found that the J/ψ and Υ disassociation temperatures are approximately 250-400 MeV and 450-700 MeV, respectively [7][8][9][10]. Before one reaches these high temperatures, however, one expects there to be partial suppression of heavy-quarkonium bound states due to in-medium breakup processes related to, for example, Landau damping and singlet-octet transitions. Because of this, one can use heavy quarkonium bound states as an internally generated probe of the QGP, with their survival probabilities depending on QGP properties such its initial temperature and the size of expected non-equilibrium deviations.
In the pioneering papers of Karsch, Matsui, and Satz (KMS) they made the first predictions that heavy quarkonia would "melt" in the QGP [3,4]. These studies were based on a non-relativistic potential-based model and the disassociation temperature of states were obtained by finding when the binding energy of the state goes to zero or r → ∞. 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 and, for sufficiently heavy quarks, e.g. bottom quarks, one can construct a non-relativistic effective field theory and then solve the non-relativistic Schrödinger equation with the resulting in-medium potential.This process can be made more formal using effective field theory methods to integrate out different energy/momentum scales, resulting in potential-based non-relativistic QCD (pNRQCD) [11][12][13][14][15].
Based on high-temperature quantum field theory calculations and complementary effective field theory calculations, it is now known 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 [16][17][18][19][20][21][22][23][24][25]. This imaginary part has been shown to be related to gluon dissociation or parton free dissociation of the states in Refs. [26,27]. The imaginary part of the potential makes the quantum evolution non-unitary (non-Hermitian Hamiltonian), which can be understood in the context of open quantum systems in which there is a heavy-quark bound state coupled to a thermal heat bath [28][29][30][31][32][33][34][35][36][37]. We note that, in the context of transport models, in-medium breakup is also included in studies of bottomonium and charmonium suppression [38][39][40][41][42][43][44].
In this paper, we 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 wavepacket trajectories. We consider only √ s N N = 5.02 TeV Pb-Pb collisions herein. For each trajectory, the states are in a quantum linear superposition and we 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. We do not include explicit time-dependent noise contributions in the potential and, as a result, for singlet evolution the obtained survival probabilities correspond to those associated with the average wave-function (averaged over thermal fluctuations) [33]. Although this is an approximation, it is a very reasonable starting point for updated phenomenological studies. Many phenomenological studies presented in the past have used this approximation to solve for the evolution of the average wave-function, however, they additionally made use of the adiabatic approximation which allows one to compute the instantaneous breakup rate for a given state from time-independent solutions to the Schrödinger equation [22,[45][46][47][48][49][50][51][52][53]. This approximation throws out potentially important physics such as quantum state mixing due to the time-dependent inmedium potential.
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.

Methodology
We solve the real-time Schrödinger equation with a complex-valued potential of the form V (r) = V R (r)+V I (r). We assume that, in vacuum, the heavy-quarkonium potential is given by a Cornell potential with finite stringbreaking distance where a = 0.409 is the effective coupling, σ = 0.21 GeV 2 is the string tension, and r SB = 1.25 fm is the string breaking distance. With this tuning of the vacuum potential, and assuming M b = 4.7 GeV, we obtain vacuum masses of The real-part of the finite-temperature single quarkantiquark potential is taken to be given by the internalenergy associated with the Karsch-Mehr-Satz (KMS) potential [45,46] where m 2 D = 4πN c (1 + N f /6)α s T 2 /3 is the in-medium gluonic Debye mass. Although there are arguments to support the use of the internal energy in thermally equilibrated systems [5,6,64], it is unclear what the correct prescription is in the non-equilibrium case. For this reason, one can consider the chosen real-part of the potential as a model choice. To match smoothly onto the zero temperature limit we use In the limit that T → 0, Eq. (3) reduces to Eq. (1). The imaginary part of the potential is taken from a leading-order resummed perturbative QCD calculation of Laine et al [16]. We evaluate the strong coupling α s at the scale µ = 2πT and use threeloop running [63] with Λ M S = 344 MeV, which reproduces the lattice result for the running coupling α s (5 GeV) = 0.2034 [65].
Using this complex potential, we then numerically solve the time-dependent Schrödinger equation on a discrete lattice. The method used is manifestly unitary for realvalued potentials and is based on a split-step pseudospectral method [66,67]. This algorithm allows for higher code accuracy and speed compared to traditional finite-size difference methods such as the Crank-Nicolson method and can be easily implemented on massively parallel architectures such as graphics cards [54]. Due to the central nature of the potential, for a fixed orbital angular momentum , we can reduce the problem to solving a one-dimensional Schrödinger equation for the scaled wavefunction u (r) = rψ (r). For the results reported herein, we used N = 4096 points with L = r max = 19.7 fm, resulting in a lattice spacing of a 0.0048 fm. We compute the in-medium suppression for = 0 and = 1 states, separately.
Due to the local nature of heavy quarkonium production, one can assume that the initial quantum mechanical wave-function is given by a Dirac delta function. Since herein, we discretize space on a finite lattice, one must regulate the delta function. 1 For this purpose, we choose a Gaussian initial wave-function with ∆ = 0.04 fm. For a given , such an initial state is a quantum superposition of many eigenstates of the Schrödinger equation. After evolving the wave-function forward in time, the probability to find a given vacuum state can be obtained by computing the overlap of the in-medium quantum wave-function with the vacuum basis states. In this way, one can obtain the survival probability of each state. Due to the fact that the Hamiltonian for this system is non-Hermitian, one finds that these overlaps decay in time, which physically reflects the in-medium breakup of bottomonium states. Note that this is different than what has been done in prior works which compute bottomonium suppression using real-time solutions to the Schrödinger-Langevin equation [31,68,69], since here the noise is encoded in the imaginary part of the potential, we use a realistic 3+1D hydrodynamics background tuned to data, and we solve the 3+1D Schrödinger equation for ensembles of trajectories.
Since each wave-packet propagating through the QGP experiences a different temperature along its trajectory, we numerically solve the time-dependent Schrödinger equation for a large set of bottomonium trajectories (1.2 million). Initial bottomonium production is Monte-Carlo sampled, assuming that the initial transverse spatial distribution is proportional to the binary overlap profile of the two colliding nuclei, N bin AA (x, y). We exploit the approximate boost-invariance of the QGP at mid-centrality and assume that all bottomonia have zero rapidity, y = 0. For the transverse momentum distribution, we assume that all states have a p T -distribution proportional to p T /(p 2 T + M 2 ) 2 , where M is the average mass of all states being considered. We assume that the initial azimuthal angle φ is distributed uniformly between 0 and 2π. Once the initial position, momentum, and azimuthal angle are sampled, we then record the QGP temperature along the trajectory followed by the quantum wave-packets. Herein, we assume that each quantum wave-packet's velocity is constant and, hence, they propagate along a straight line trajectory.
For the background temperature evolution, we use the output of a 3+1D quasiparticle anisotropic hydrodynamics (aHydroQP) code which has been tuned to reproduce soft hadron multiplicities, elliptic flow, etc at √ s N N = 5.02 TeV [60,62]. We use smooth optical Glauber initial conditions and the parameters used for the aHydroQP runs correspond to an initial central temperature of T 0 = 630 MeV at τ 0 = 0.25 fm/c, with a constant specific shear viscosity of 4πη/s = 2 [62]. For each trajectory sampled, we  evolve the quantum state using the in-medium complex potential until the local temperature is below the QGP transition temperature, T QGP = 155 MeV. We evolve the quantum wave-packets using the vacuum potential starting at τ = 0 fm/c and turn on the in-medium potential at τ = τ med = 0.4 fm/c. Whenever the temperature on a given trajectory drops below T QGP , we use the vacuum potential for its evolution.
After each quantum state is propagated along its trajectory, we convert the survival probabilities into particle number by multiplying by (1)    [63]. 3 For example, the final number of 1s states produced can be computed as N QGP . Note that each column of F must sum to unity in order to preserve bottom number.
To compute R AA , we divide the final number of bottomonium states produced by the number of binary collisions in the sampled centrality class times the post feeddown pp production cross-section for each state. Since we know the reaction plane (provided by aHydroQP) one has Ψ RP = 0 and, as a result, one can compute v n by simply averaging cos(nφ) over all particles, v n ≡ cos(nφ) , in a given p T and centrality bin. For both R AA and v 2 , we report the statistical uncertainty associated with the sum over the sampled quantum wave-packet trajectories. The resulting model will be referred to as Heavy Quarkonium Quantum Dynamics (HQQD) in what follows.

Results
In Fig. 1 we present HQQD predictions for the suppression of Υ(1s), Υ(2s), and Υ(3s) states as a function of N part . For this Figure, in HQQD we applied a transverse momentum cut of p T < 30 GeV. We compare with results obtained by the ALICE [70], ATLAS [71], and CMS [72] collaborations, shown as circles, squares, and triangles, respectively. From this Figure, we see that HQQD does a quite reasonable job in describing the N part dependence of R AA [Υ(1s)], however, HQQD predicts a somewhat smaller R AA [Υ(2s)] than the experimental results. Similar conclusions can be obtained from    [76]. From the results shown in Fig. 2, we see that HQQD predicts a very weak dependence of R AA [Υ] on p T , with only a small decrease at momentum less than the mass scale of the bottomonium states. The increased suppression at low-p T can be attributed to such wave-packets having, on average, a longer effective lifetime inside the QGP fireball (due to their lower velocities).
In Fig. 3, we present our results for the elliptic flow of Υ(1s), Υ(2s), and Υ(3s) states as a function of centrality. For this Figure, we impose p T < 50 GeV and compute v 2 in 10 equally spaced centrality bins from 0-100%. The bands in this Figure show the statistical uncertainty associated with the mean values extracted in each bin. As can be seen from this Figure, there is a clear ordering of the elliptic flow, with the Υ(3s) state having the largest flow and the Υ(1s) the smallest. This is in agreement with expectations, since the source of the elliptic flow in all cases is the suppression of the states and, hence those with stronger suppression will have a larger elliptic flow. One other thing that is evident from Fig. 3 is that the elliptic flow of all states goes to zero for central collisions (left hand side of the plot). This, of course, is a consequence of our choice of non-fluctuating optical Glauber initial conditions and provides a non-trivial test of the HQQD calculation of v 2 . If one includes geometric fluctuations in the initial hydrodynamic variables (energy density, etc.), one would expect to see small, but finite, values for the elliptic flow of all states in central collisions. On the right hand side of Fig. 3 one sees that the elliptic flow for all states goes to zero. This, again, agrees with expectations since the QGP lifetime in such events is zero.
One other feature visible in Fig. 3 is the non-monotonic We turn next to Fig. 4 in which we present a comparison of HQQD predictions for v 2 [Υ(1s)] with experimental data collected by the ALICE [77] and CMS [78] collaborations in three different transverse momentum bins: 0-4, 4-6, and 6-15 GeV. For both HQQD and the experiments, the results are integrated over centrality in the range 5-60%. As can be seen from this Figure, HQQD predicts a result consistent with zero in the lowest p T bin, a slightly negative result in the central bin, and a small but positive value in the highest momentum bin. This trend (positive near zero, then negative, and then positive again) and the overall magnitude of v 2 predicted by HQQD is similar to what has been predicted previously using a model which relies on the adiabatic approximation [53]. In Ref. [53] it was posited that the explanation for this negative v 2 is related to the transverse expansion of the QGP overtaking bottomonia states which have escaped from near the surface of the QGP. 4 With respect to the comparisons with experimental data, we find reasonable agreement with available data, given current experimental uncertainties, and one sees a similar trend in the three centrality classes as predicted by HQQD.
In Fig. 5, we present a comparison of HQQD with experimental data from the CMS collaboration for the centrality dependence of v 2 [Υ(1s)]. All results are binned into three centrality bins: 10-30%, 30-50%, and 50-90%. In the rightmost panel of Fig. 5, we show the experimental result integrated over 10-90% centrality compared to the HQQD prediction in the same centrality interval. From  this Figure we see that the integrated v 2 [Υ(1s)] in the 10-90% class is in agreement, within uncertainties, with the experimental data provided by CMS. In the separate bins (left panel ), we see good agreement in the 10-30% bin, however, in the other two bins we larger differences, albeit still within 2σ of the HQQD predictions. In the future, hopefully higher statistics will allow for more constraining comparisons between HQQD and experiment.
In Fig. 6, we present HQQD predictions for v 2 [Υ(2s)] and v 2 [Υ(3s)] in the same centrality bins as Fig. 5. For v 2 [Υ(2s)], there is currently only one integrated data point available from the CMS collaboration, which is shown as a green triangle in the 10-90% panel (right). Comparing the integrated results, we see that v 2 [Υ(2s)] is currently within the reported experimental uncertainties, however, at the very top end of them. Again, increased statistics will allow for more accurate comparisons in the future. In the left panel of Fig. 6 we see that the flow of v 2 [Υ(3s)] can be on the same order of magnitude as the experimentally observed v 2 [J/ψ] [77,78].
In Fig. 7, we present HQQD predictions for v 2 [Υ(2s)] and v 2 [Υ(3s)] as a function of transverse momentum using the same p T -bins as Fig. 4 in order to allow for easier comparison with experimental data in the future. From this Figure we see that the Υ(3s) can develop a sizable v 2 solely due to path length differences between the short and long sides of the QGP fireball. Turning to the Υ(2s) we see that, similar to the Υ(1s), HQQD predicts a negative v 2 in the lowest two p T -bins. This once again is related to the fact that the QGP expands more rapidly along the short side than the long side, which can have the affect of overtaking bottomonium states which had previously escaped the QGP with φ ∼ 0. In the highest p T -bin shown, we see that HQQD predicts positive v 2 for both states.
Finally, in Table 1 we present comparisons between HQQD predictions for various observables and the corre-  sponding experimental results from the ALICE, ATLAS, and CMS collaborations. In this Table, the results are integrated over centrality and transverse momentum in the ranges shown in the middle column and the last row shows the HQQD prediction for v 2 [Υ(3s)]. We do not indicate the rapidity cuts used by each experimental collaboration, which correspond to 2.5 < y < 4.0, |y| < 1.5, and |y| < 2.4 for the ALICE, ATLAS, and CMS collaborations, respectively. From this Table we see that all HQQD predictions are within the combined statistical and systematic uncertainties reported for each measurement. We once again note that, for the Υ(2s), HQQD seems to predict slightly too much suppression, however, the HQQD predictions are still compatible with experimental results within uncertainties.

Conclusions and outlook
In this paper we used real-time quantum evolution to compute the suppression and elliptic flow of bottomonium states and presented the details of the resulting HQQD model. Using HQQD, we sampled a large set of bottomonium trajectories (1.2 million). For the HQQD hydrodynamic background, we used anisotropic hydrodynamics to provide the 3+1D temperature field through which the states were propagated. Given this background, we then solved the time-dependent Schrödinger equation with a complex potential and obtained the survival probability for each bottomonium state by computing the quantum mechanical overlap of the in-medium evolved wave-packet with the vacuum eigenstate of the state of interest.
After averaging over all wave-packet trajectories, we were able to obtain precise estimates for R AA which are in quite reasonable agreement with available experimental data. We then presented predictions of HQQD for the elliptic flow of Υ(1s), Υ(2s), and Υ(3s) and compared those to available experimental data. In the case of v 2 we, once again, found reasonable agreement between theory and experiment, with the integrated v 2 being consistent with the experimental data within uncertainties (Table 1). With respect to v 2 , we emphasized two model observations: (1) that v 2 for the various states can be non-monotonic (oscillating) due to quantum mechanical oscillations in the timeevolved overlaps and (2) that v 2 for Υ(1s) and Υ(2s) can be negative in intermediate transverse momentum bins, e.g. 4 < p T < 6 GeV. The first observation is novel, however, the second has been observed previously in calculations of v 2 using the adiabatic approximation [53]. In order to make the HQQD predictions presented herein, we had to make a set of model choices corresponding to, for example, the choice of the real and imaginary parts of the quarkantiquark potential and the precise initialization time for medium interaction τ med . We did not attempt to estimate the systematic uncertainties associated with these choices, but plan to in a forthcoming paper.
Our use of real-time solutions allowed us to go beyond the adiabatic approximation. Overall, we found our HQQD results to be qualitatively consistent with previous adiabatic approximation results, however, with HQQD one has a more complete description of the quantum dynamics. The use of real-time solutions allowed us, for example, to include the effect of quantum-mechanical state mixing due to the time-dependent in-medium quark-antiquark potential. Looking to the future, we plan to more fully include the effect of thermal noise in the underlying evolution. In this work, we evolved the states with a complex Hamiltonian which is appropriate for describing the evolution of the average wave-function of the system. As a result, the system remains in the singlet configuration and there can be no transitions between different angular momentum states.
The description in terms of the average wave-function is not yet a complete description, however, in practice one finds that the in-medium wave-function evolution for the ground state is well-approximated by the evolution of the average wave-function subject to a complex Hamiltonian when including noisy potentials [79,80]. For excited state suppression, it may be important to go beyond the complex Hamiltonian approach used herein. For this purpose, one can either introduce a noisy potential, with the noise spectrum set by the imaginary part of the quark-antiquark potential [28,29,33,[81][82][83][84][85] or instead solve the resulting Lindblad equation including both singlet and octet states in order to describe the evolution of the full density ma-trix [30,32,[34][35][36]. Preliminary results obtained using the quantum trajectories method to solve the Lindblad equation indicate that the more accurate inclusion of noise effects and singlet-octet transitions will result in only small changes in R AA [Υ(1s)], however, these studies indicate that the excited states are less suppressed when including initial state octet production and in-medium stochastic singlet-octet transitions [79]. This will hopefully result in better agreement between HQQD and experimental results for R AA [Υ(2s)]. Finally, we mention again that it would also be interesting to study the effect of geometric fluctuations in the initial state on v 2 [Υ].