Bottomonium suppression using a lattice QCD vetted potential

We estimate bottomonium yields in relativistic heavy-ion collisions using a lattice QCD vetted, complex-valued, heavy-quark potential embedded in a realistic, hydrodynamically evolving medium background. We find that the lattice-vetted functional form and temperature dependence of the proper heavy-quark potential dramatically reduces the dependence of the yields on parameters other than the temperature evolution, strengthening the picture of bottomonium as QGP thermometer. Our results also show improved agreement between computed yields and experimental data produced in RHIC 200 GeV/nucleon collisions. For LHC 2.76 TeV/nucleon collisions, the excited states, whose suppression has been used as a vital sign for quark-gluon-plasma production in a heavy-ion collision, are reproduced better than previous perturbatively-motivated potential models; however, at the highest LHC energies our estimates for bottomonium suppression begin to underestimate the data. Possible paths to remedy this situation are discussed.


Ultrarelativistic heavy-ion collision (URHIC) experiments being performed at the Large Hadron Collider (LHC) at CERN and the Relativistic Heavy Ion Collider (RHIC) at
Brookhaven National Laboratory aim to recreate a primordial state of nuclear matter known as the quark-gluon plasma (QGP). Comparisons between hydrodynamic simulations that describe the evolution of the bulk matter and experimental data suggest that LHC URHICs produce a QGP with an initial temperature on the order of T 0 = 600−700 MeV at τ 0 = 0.25 fm/c [1][2][3][4][5], for beam energies in the range of 2. 76 -5.02 TeV/nucleon. In addition, analysis of the collective flow of the matter produced in URHICs indicates that the QGP behaves like a nearly inviscid relativistic fluid [1][2][3][4][5]. Although the light hadronic states, such as pions, are disassociated at temperatures around the pseudocritical temperature T c 155 MeV, it was predicted that, due to their large binding energies, bound states of heavy quarks could survive up to temperatures on the order of 600 MeV. Due to their short formation times, such heavy-quark bound states probe the history of the QGP and their suppression relative to production in pp collisions was proposed to be a clear signal of the creation of the QGP stemming from temperature-dependent screening of the color force [6,7]. In practice, experimentalists do see a reduced yield of heavy quarkonium compared to elementary collisions both at RHIC and LHC [8,9].
As we have learned over the course of the last years, however, the suppression of heavy quarkonium may be accompanied by the regeneration of the observed quarkonium states (a) at the phase boundary between the QGP and the hadronic phase [10][11][12] or (b) dynamically due to in-medium recombination of QQ pairs [13,14]. If a sizable number of QQ pairs is created in the initial stages of a collision, the probability for a free quark antiquark pair to coalesce into a bound state at some point during the QGP lifetime can become significant. Measurements of charmonium yields e.g. have shown that for this lighter flavor indeed regeneration becomes quite important at LHC energies [15]. For bottomonium, only recently have there been studies carried out that include a regeneration component, see e.g. [16]. In the case of bottomonium, one expects regeneration to be less important due to its much heavier rest mass and large vacuum binding energy ( 1 GeV) and this is borne out by detailed calculations. That said, it is desirable to have a unified framework that includes the effects of both in-medium suppression and regeneration.
Recent measurements at LHC have also shown an unambiguous signal for the elliptic flow of the J/Ψ particle, which implies that charm quarks are at least partially kinetically equilibrated with the surrounding medium [17][18][19]. A similar observation for bottomonium has not been made and it is expected that bottomonium does not yet participate in the collective motion of the bulk at 5.02 TeV. Since equilibration is intimately related to a loss of memory, charmonium at LHC appears to provide us with a window into the late stages of the collision, while bottomonium is considered to act as probe of the full evolution of the QGP. As a result, a detailed understanding of heavy quarkonium in-medium will open the possibility to infer the time-dependent properties of the nuclear matter produced in a heavy-ion collision.
It is an open question which properties the different species are most sensitive to, as candidates of course temperature or the shear viscosity to entropy ratio come to mind.
In this paper, we focus solely on bottomonium and investigate its suppression at both RHIC and LHC energies. To this end we combine a non-relativistic description of the quarkantiquark bound state in terms of an in-medium potential with a realistic dynamical model of the bulk matter created in the collision. The potential model used is vetted by comparing it to lattice QCD calculations of the real and imaginary parts of the heavy-quark potential, providing the first model to attempt to constrain the full complex potential using lattice input. In addition, since bottomonium states are believed to be formed early on in a collision (τ < 1 fm/c), they can be sensitive to the early-time non-equilibrium dynamics of the QGP.
The structure of our paper is as follows. In section II we review how the complex-valued heavy-quark potential used for the description of bottomonium is defined and extracted using lattice QCD and appropriately parametrized for use in phenomenological applictions.
Section III collects the relevant details on the dynamical evolution model of anisotropic hydrodynamics and how it is connected to the lattice-vetted potential via a spacetime-dependent anisotropic Debye mass. We present the results of our computation in section IV and conclude with a discussion and outlook in section V.

II. POTENTIAL DESCRIPTION OF HEAVY-QUARKONIUM
Herein, we use a lattice QCD vetted, non-relativistic potential based description of heavyquarkonium in a thermal medium to compute its real-time evolution in a heavy-ion collision.
The potential originates in a systematic treatment of heavy quarkonium in QCD, based on the effective field theories non-relativistic QCD (NRQCD) and potential NRQCD [39,40].
These frameworks exploit the inherent separation of scales between the heavy quark rest mass m c,b , the surrounding medium temperature, as well as the characteristic scale of QCD Λ QCD to dramatically simplify the description of heavy quarkonium. Instead of having to consider a full quantum field-theoretical boundary value problem for Dirac spinor fields, one may go over in a first step to an initial-value problem for two-component Pauli spinors (NRQCD). In a second step this non-relativistic theory may be matched to a further simplified description in terms of coupled color singlet ψ S (r, t) and color octet wavefunctions ψ O (r, t) (pNRQCD).
In the latter, the interaction among the heavy-quarks, as well as their interaction with the medium is captured in both potential and non-potential contributions. Depending on the concrete scale hierarchy of the problem at hand, the potential contributions may dominate and the relevant real-time evolution of heavy quarkonium reduces to a Schrödinger equation.
For realistic settings, such as those encountered in heavy-ion collisions at RHIC and LHC, the matching coefficients of the effective theory, i.e. the potential cannot be reliably determined using perturbation theory. Nevertheless vital insight had been gained by evaluating pNRQCD using the hard-loop approximation [40,41]. In particular it was found that the proper in-medium potential must assume complex values at high temperatures, as was first discussed in [41] and subsequently extended to a momentum-space anisotropic medium in [22][23][24]32]. This fact, in particular, implies that purely real model potentials, such as the popular color singlet free energies or the internal energies do not constitute valid descriptions of the relevant in-medium quarkonium physics.
The proper heavy-quark potential is related to a real-time QCD quantity, the rectangular Wilson loop, via the process of matching, i.e. one selects a correlation function in the effective theory pNRQCD and in the underlying microscopic theory QCD which carry the same physics content and identifies them at an appropriate scale. In our case the unequal time correlation function of a heavy quarkonium singlet state may be identified with the Wilson loop in the static limit Since the Wilson loop obeys a simple equation of motion [41] i∂ t W (r, t) = Φ(r, t)W (r, t), with a time-and space-dependent complex function Φ(r, t), the potential picture is applicable as long as Φ asymptotes towards a time independent value at late times. This value in general is complex and the corresponding potential may be formally defined as As such, this real-time definition is not yet amenable to an evaluation in non-perturbative lattice QCD, which is simulated in unphysical Euclidean time. Instead, one has to take a detour via the spectral decomposition of the Wilson loop to relate the Euclidean and Minkowski domain [42,43] W (τ, r) = dωe −ωτ ρ (ω, r) ↔ dωe −iωt ρ (ω, r) = W (t, r).
Inserting (3) into (2) tells us that the real and imaginary part of the potential are related to the position and width of the lowest lying peak structure within the Wilson loop spectrum.
It has been shown [44] that if a potential picture is applicable, the Wilson loop spectrum actually contains a well defined lowest lying peak of skewed Lorentzian form from which the values of the potential can be straightforwardly extracted via a χ 2 -fit.
Note however that the extraction of spectral functions from Euclidean lattice data is an illposed inverse problem, which has only recently been successfully tackled in the context of the heavy-quark potential. With the help of a novel Bayesian approach [45] [46], as well as for full QCD with N f = 2 + 1 light medium quark flavors based on simulations by the HotQCD collaboration. In both cases, the applicability of the potential picture was confirmed, at all temperatures considered, as the Wilson spectral functions showed a well defined peak of Lorentzian shape.
To utilize the discrete values of the potential obtained from lattice QCD two further steps need to be taken, for which we follow the same strategy as laid out in [47,48].
Three parameters enter this expression, which characterize the non-perturbative vacuum physics of the quarkonium bound state: the strong coupling α s , the string tension σ, and a constant shift c. Note that we absorb the factor C F into our definition of the strong coupling α s = g 2 C F 4π . In order to introduce the effects of a thermal medium, we adopt a prescription well known in classical electrodynamics, i.e. in the case of the Coulombic contribution, one Fourier transforms Gauss' law and subsequently modifies the right hand side by dividing it with an in-medium permittivity . We use here the permittivity of the QCD medium computed in hard-thermal loop perturbation theory The idea is that the non-perturbative physics of the bound state is encoded in the Cornell form of the T = 0 potential, whose modification is driven by a weakly-coupled gas of quarks and gluons. Combining (5) and (4) thus leads to integro-differential equations for the inmedium modified Coulombic and string-part of the T = 0 potential as discussed in detail in [49]. Since the permittivity is complex, the in-medium potential also contains an imaginary part. In contrast to purely perturbative computations, which capture only the Coulombic contribution to the potential, the in-medium potential here receives additional contributions to its real and imaginary part from the string-like part of the T = 0 Cornell potential. The explicit expressions for the Coulombic part are with which coincide with the results of Ref. [41]. The additional and novel string-like contribution on the other hand reads for the real part, where the strength of the in-medium modification is characterized by the parameter µ 4 = m 2 D σ/α s . For its imaginary part we have where ψ corresponds to the following Wronskian An important aspect of these expressions is that, once the vacuum parameters of the Cornell potential are fixed, only a single temperature-dependent parameter remains, the Debye mass While the analytic parametrization was derived in a straightforward fashion, it relies on several assumptions and thus needs to be validated on real lattice QCD data. Using both quenched [46] and full QCD simulations with N f = 2 + 1 light flavors [47] it was shown that the lattice values of the potential were indeed excellently reproduced by the generalized Gauss law parametrization (see Fig.1). After fixing α s , σ, and c using low temperature ensembles, the real-part of the potential was fitted by tuning m D . Once m D is fixed, the parametrization makes a prediction for Im [V], which for quenched QCD simulations showed quantitative agreement at high temperatures and, as expected, became less accurate around the phase transition. In full QCD, no robust determination of the imaginary part has been achieved so far, however, the tentative values extracted, again, showed very good agreement with the Gauss-law parametrization at high and intermediate temperatures [47].
The values for the Debye mass related to the full QCD in-medium potential also showed clear deviations from the perturbative predictions in the phenomenologically relevant regime between T C < T < 3T C .
Even though it has been established that the Gauss-law parametrization is capable of reproducing the lattice QCD in-medium potential, its values found on the lattice may not be applied directly to phenomenological computations due to the presence of lattice artifacts.
If a continuum extrapolation of the T = 0 potential in the thermodynamic limit were available, we could directly determine the parameters α s , σ and c from first principles.
Here instead we select these parameters in a phenomenological fashion, i.e. we tune them such that the vacuum bottomonium spectrum below the B-meson threshold is reproduced.
The correct quark mass to be used in such a computation is the renormalon subtracted In this study we will compute the Debye mass self-consistently from the dynamical evolution of the bulk and use its value to implement the in-medium modification of the Cornell potential with the above parameters. In Fig. 2 we show the real (left) and imaginary part (right) of the actual potential used for different values of the Debye mass. In order understand better which values of m D play a role in the evolution of heavy quarkonium we note that lattice QCD studies showed that in a thermal QCD medium close to the crossover transition the ratio of m D /T ≈ 1 and grows to m D /T ≈ 2 as temperature is increased to Real-part of the heavy-quark potential from the Gauss-Law α_s=0.504 (29) σ =0.415 (15)GeV Imaginary-part of the heavy-quark potential from the Gauss-Law With the real and imaginary parts of the potential specified, we can now turn to the method used to fold the dynamical evolution of the full three-dimensional QGP evolution together with information about the real and imaginary parts of the resulting binding energies. In order to do so, however, we will must first make an extension of the result presented in the previous section to the case of a QGP with momentum-space anisotropies. This is necessary for a realistic model of QGP evolution and quarkonia suppression. The simplest form for the soft-particle distribution function that can be used to take into account QGP momentum-space anisotropies is a generalization of an isotropic phase-space distribution which is squeezed or stretched along one direction in momentum space, defined byn, with a parameter −1 < ξ < ∞. In a heavy-ion collision the directionn can be identified with the beam line direction,n =ẑ. The resulting one-particle distribution function is given by the following spheroidal "Romatschke-Strickland" form [33,50] f (p, ξ, Λ) ≡ f iso ( p 2 + ξ(p ·n) 2 /Λ) , where Λ is the transverse temperature scale and f iso is the isotropic thermal distribution function associated with the soft degrees of freedom in the QGP.

A. Anisotropic Debye mass
The Debye mass, which we use to evaluate the heavy-quark potential in our approach, is self consistently computed from the dynamical evolution of the soft bulk matter. The conventional isotropic m D is defined via an integral over the isotropic distribution function where p 2 ≡ p 2 = p 2 ⊥ + p 2 z . Since lattice studies are restricted to a thermal and isotropic state, we have to use additional input to incorporate the effects of a QGP momentum-space anisotropy on the potential itself. For this purpose, we use the results contained in Sec. 2.2 of Ref. [26] where it was shown that, empirically, one can incorporate the momentum-space anisotropy parameter ξ into the potential by making the Debye mass depend on the angle with respect to the beam line direction. The result obtained was with a = 16/π 2 , b = 1/2, d = 3/2, e = 1/3, and With these values in hand, we simply tabulate the real and imaginary parts of the potential functions with respect to m D and r, and then replace m D → µ in the isotropic potential to obtain the corresponding anisotropically modified potential. The parameter ξ sets the degree of momentum-space anisotropy, with −1 < ξ < 0 corresponding to a prolate distribution, with more momentum along the beam line direction than transverse to it and ξ > 0 corresponding to an oblate distribution. The anisotropy parameter is dynamically tied to the bulk QGP evolution using anisotropic hydrodynamics (aHydro) which we will discuss shortly [5,21,[34][35][36][37][38][51][52][53][54][55].
Since the potential no longer has full rotational symmetry one has to go beyond a single (radial) dimension when solving the Schrödinger equation. With the assumed form for the one-particle distribution function, only one symmetry direction is broken (spheroidal symmetry), and one can, in principle, simply use a two-dimensional solver. Here we have chosen to make use of a full three-dimensional solution since codes for this are already available even for complex-valued potentials [56,57]. For this purpose, the real-time Schrödinger equation is solved in imaginary time. The grid spacing for Υ(nS) states is chosen as a = 0.15 fm on N = 256 3 regularly-spaced grid points. Due to the larger size of p-wave bottomonia, our grid spacing for the χ b (mP ) states is set to a = 0.175 fm. Starting from a randomized three-dimensional wavefunction, we evolve in imaginary time until the ground state wavefunction converges to within a given tolerance. Using "snapshots" of the wavefunction stored during the imaginary-time evolution, we can project out the low-lying excited states [56].
Additionally, by fixing the symmetry (symmetric vs anti-symmetric) of the initial random wavefunction we can select the s-wave or p-wave states independently [56]. In this way we can obtain the wavefunctions of the 1s, 2s, 3s, 1p, and 2p Upsilon states and in turn we can compute the real and imaginary parts of their respective energies. For more details concerning the numerical method, we refer the reader to Refs. [56,57].

C. Anisotropic hydrodynamics equations and initial conditions
To proceed, we assume that the underlying bulk one-particle distribution function is well approximated by Eq. (11) at all points in spacetime. In addition, we assume that the bulk evolution is well described using hydrodynamical degrees of freedom, such as energy density, pressures, and viscous corrections. This assumption has been tested by comparing the predictions of anisotropic hydrodynamics to exact solutions of the Boltzmann equation [58][59][60][61][62][63][64][65][66][67] (see also [68] for more comparisons to kinetic theory) where it has been shown to reproduce the exact kinetic evolution even far from equilibrium. The equations used herein are obtained using the zeroth and first moments of the Boltzmann equation in the relaxation-time approximation. For details about the dynamical equations used and their physical content, we refer the reader to Sec. IV of Ref. [69].
In order to solve the aHydro dynamical equations, one has to make a reasonable assumption about the initial conditions at the initial longitudinal proper-time for the hydrodynamic evolution, τ = τ 0 . For our initial conditions, we take the system to be isotropic in momentum space (ξ = 0) with zero transverse flow and Bjorken flow in the longitudinal direction. In the transverse plane, the initial energy density is computed from a linear combination of smooth Glauber wounded-nucleon and binary-collision profiles with a binary mixing factor of α = 0.15. In the longitudinal direction, we used a "tilted" profile with a central plateau and Gaussian tails resulting in a profile function of the form , with ς = arctanh(z/t) being spatial rapidity. The parameters entering the longitudinal profile function were fitted to the pseudorapidity distribution of charged hadrons with the results being ∆ς = 2.3 and σ ς = 1.6. The first quantity sets the width of the central plateau and the second sets the width of the Gaussian "wings". The resulting initial energy density at a given transverse position x ⊥ and spatial rapidity ς was computed using where W A,B (x ⊥ ) is the wounded nucleon density for nuclei A and B, C(x ⊥ ) is the binary collision density, and g(ς) is the "tilt function". The tilt function g(ς) = 0 if ς < −y N , g(ς) = (ς + y N )/(2y N ) if −y N ≤ ς ≤ y N , and g(ς) = 1 if ς > y N where y N = log(2 √ s N N /(m p + m n )) is the nucleon momentum rapidity [70].

D. Computing the suppression
The solution of the 3+1d aHydro dynamical equations gives us hard momentum scale Λ, the momentum-anisotropy parameter ξ, and, consequently, the anisotropic Debye mass µ as a function of proper time, transverse coordinate x ⊥ , and spatial rapidity ς. Solving the Schrödinger equation in addition provides the real and imaginary parts of the binding energy of a given state as a function of Λ and ξ. Folding these together gives us the real and imaginary parts of the binding energy as a function of proper time, transverse coordinate x ⊥ , and spatial rapidity ς: Re[E bind (τ, x ⊥ , ς)] and Im[E bind (τ, x ⊥ , ς)], respectively. Full details of the method used can be found in Ref. [26]. Here we summarize the important points.
We use the imaginary part of the binding energy to provide information about the decay rate of a given state and the real part of the binding energy to tell us when the state becomes completely unbound. The imaginary part of the binding energy can be related to the decay rate, Γ, using ∂ τ n i = −Γ i n i , where i indexes the state in question, giving Γ i = −2Im[E i ] [26]. Finally, using the fact Im[E bind ] = −Im[E] we obtain The value of 10 GeV in the second case was chosen in order to quickly suppress states which are fully unbound and, in practice, the results do not depend significantly on this value as long as it is large enough to quickly disassociate the state under consideration.
We then integrate the instantaneous decay rate, Γ, obtained in this manner over propertime to extract the dimensionless logarithmic suppression factor where τ form (p T ) is the lab-frame formation time of the state in question. The formation time of a state in its local rest frame can be estimated by the inverse of its vacuum binding energy [71].  factor For implementing cuts in centrality we compute R AA for finite impact parameter b and map centrality to impact parameter in the standard manner. Finally, to compare with experimental observations we average R AA (x ⊥ , ς) over x ⊥ . For this operation, we use a production probability distribution which proportional to the overlap density [26]

IV. RESULTS
In this section, we present our numerical results and compare them to data from RHIC and LHC. The background dynamics discussed in the previous section have been determined for three different values of the shear viscosity to entropy ratio varied around the proposed quantum bound η/s = 1/4π, with the initial temperature tuned so that the soft particle multiplicity is held fixed. The corresponding initial temperature values for the appropriate beam energies are listed in Table I and lead to hydrodynamic evolution, which reproduces light particle spectra and azimuthal flow well.  there are always states at the edge of the plasma where the temperature/density is low where the states are largely unsuppressed. This results in a kind of universal "halo" survival probability which is related to the geometry of the fireball and its temperature profile.

Feed down fractions
Since it is not the primordial bottomonium states themselves that are measured in the detector but instead the decay di-leptons, feed down must be factored to calculate the inclusive R AA for each state observed in experiment as described in the previous section.
The relevant feed down fractions are obtained from averaging over p T from experimental feed down yields [72] and are listed in Table II. The inclusive Υ(1S) R AA is presented in Fig. 4  (Color online) bottomonium suppression for the LHC run2 parameter set based on the perturbative model of [56] (top three curves) and the lattice-vetted heavy-quark potential (bottom three curves). Note that due to the large imaginary part in the lattice-vetted potential the suppression is consistently larger in than in the perturbative model potential. Interestingly we find that virtually no dependence on the chosen values for the shear viscosity to entropy ratio is observed with the lattice-vetted potential. This behavior is found both at RHIC and LHC energies.
from the perturbative potential. The reason is that the former features a stronger imaginary part and thus the bottomonium states are more easily dissociated. Secondly we find that that R AA computed with the lattice-vetted heavy-quark potential is virtually independent of the η/s parameter of the aHydro background evolution. This behavior is consistently observed both at RHIC and LHC energies and has important consequences for the role bottomonium can play as a probe of the QGP. The less the suppression depends on parameters other than the temperature, the more bottomonium can be used as a as dynamical thermometer of nuclear matter under extreme conditions.
In order to meaningfully compare to experimental results, we need to quantify the uncertainty in the used potential. When vetting the potential with lattice QCD simulations it was observed that the lattice potential values could be fitted with a Debye mass with around 10-20% error. We, therefore, have repeated our calculations including a modest ±15% variation of m D which leads to an error estimate on the R AA theory curve quantifying the potential uncertainty. Since our results are essentially independent of the shear viscosity parameter, we set 4πη/s = 2 consistent with recent particle spectra fits [5].  estimates agree with the data within the still relatively large error bars but are slightly lower than the experimental data. On the other hand, the trend in the excited state data points is excellently reproduced, touching also the point at the lowest centrality bin, providing a better description overall than the perturbative model results (see Ref. [31] for a compilation of the prior results).
That being said, when moving to the higher energy of run2 at the LHC, we find that the trend of stronger suppression continues in our estimates of bottomonium suppression. At this energy the lattice-vetted model overpredicts the amount of suppression of both the Υ(1S) and Υ(2S). This means that now our R AA systematically overestimates the suppression for both states. The discrepancy is larger for more central collisions while for smaller N part we still find reasonable agreement with the data.
In Fig. 7 we plot the nuclear modification factor as function of p T integrated over all centrality classes at (left) 2.76 TeV and (right) 5.02 TeV. Similarly to our findings in terms of centrality at LHC run1 energies, the agreement here is best for Υ(2S) and acceptable for Υ(1S). There is a tendency visible to slightly overestimate the suppression when using the lattice-vetted potential in the computation.
Finally, in Fig. 8   We have found that our estimates for the suppression of primordial bottomonium show the best agreement with experimental results at RHIC energies, where they reproduce the measured data within errors. The larger imaginary part in the lattice-vetted potential compared to previously used model potentials leads to a stronger suppression, which in the case of RHIC energies is the reason for the excellent reproduction of the experimental data.
At LHC run1 we find, on the other hand, hints of a tendency to overestimate the suppression compared to the data, while at 5.02 TeV we clearly see too much suppression when using the lattice-vetted potential. Two reasons for this behavior immediately come to mind. The first is related to the physics mechanism we assume underlies the imaginary part of the potential, the second one is related to the fact that no regeneration component has been included in our current study.
Our calculation uses a simple implementation of bottomonium suppression as discussed in Sec. III. By using the imaginary part of the potential directly in the Schrödinger equation we have assumed that it arises solely from gluo-dissociation of the heavy-quarkonium. Studies of bottomonium in perturbative pNRQCD have shown however that the imaginary part contains both contributions from gluo-dissociation and Landau damping, which leads to the excitation of bottomonium states without decay of the QQ pair. Eventually we will need to understand how to disentangle these two contributions, since the latter, as shown in exploratory studies in the context of open-quantum systems [73,74] leads to a weaker suppression. I.e. our inability to disentangle the underlying mechanisms of the imaginary part have led us to assign it fully to a loss channel which may overestimate the actual suppression.
One possible way for future improvement of the proper potential based real-time description of the heavy quarkonium evolution, would be to utilize the framework of open-quantum systems. As was proposed e.g. in [75], the imaginary part of the heavy-quark potential, which arises in the description of the unequal time correlation functions of Eq. (1) is unraveled into either a stochastic dynamics on the level of the Schrödinger equation or equivalently into a master equation for the density matrix of quarkonium states. Recent developments in this field include a formulation of a Lindblad type coupled dynamics for color singlet and octet degrees of freedom [76][77][78], which has been evaluated in a perturbative setting in [79]. A first attempt to describe bottomonium dynamics was presented in [73], however, the medium evolution used in this paper was not yet realistic.
In addition, we have estimated the suppression of bottomonium states without including a possibility for regeneration. While for charmonium there are clear indications for the presence of a regeneration components in the observed yields [15], no such signals have been firmly established for bottomonium yet. One hint towards the onset of regeneration could be the small change in the bottomonium suppression when going from 2.76 TeV to 5.02 TeV.
While the lattice-vetted potential predicts suppression that increases with beam energy, the data appears to not decrease as rapidly as suggested by our calculations. This would allow for the addition of a regeneration component to approach the data from below, however, we are not yet in a position to quantitatively estimate the magnitude of this effect.
In the end it will be necessary to explore both paths to come to a robust phenomenological understanding of the observed yields. The model of bottomonium suppression needs to be made more flexible in order to accommodate the different underlying physics mechanisms inherent in the imaginary part, e.g. via a stochastic Schödinger equation and the inclusion of regeneration via a rate equation framework is desirable.