Bottomonium suppression from the three-loop QCD potential

We compute the suppression of bottomonium in the quark-gluon plasma using the three-loop QCD static potential. The potential describes the spin-averaged bottomonium spectrum below threshold with a less than 1% error. Within potential nonrelativistic quantum chromodynamics and an open quantum systems framework, we compute the evolution of the bottomonium density matrix. The values of the quarkonium transport coefficients are obtained from lattice QCD measurements of the bottomonium in-medium width and thermal mass shift; we additionally include for the first time a vacuum contribution to the dispersive coefficient ${\gamma}$. Using the three-loop potential and the values of the heavy quarkonium transport coefficients, we find that the resulting bottomonium nuclear modification factor is consistent with experimental observations, while at the same time reproducing the lattice measurements of the in-medium width.


I. INTRODUCTION
The study of the deconfined quark-gluon plasma (QGP) is the main goal of heavy-ion collision experiments.It provides unique insight into some sectors of the phase diagram of Quantum Chromodynamics (QCD) and a window to conditions that existed in the early universe.Heavy quarkonium suppression in the QGP is one of the core processes to probe the properties of the QGP, as heavy quarkonium is amenable to a clean experimental reconstruction and a robust theoretical description.The idea of heavy quarkonium suppression traces its origin to Matsui and Satz, who proposed that a screening of the quark-antiquark potential for J/ψ mesons traveling through the QGP produced in heavy-ion collisions should lead to a lower number of particles measured compared to proton-proton collisions [1].In the past, this effect has been studied extensively for bottomonium [2], where due to the heavy bottom quark mass, various calculational simplifications are possible, including the use of effective field theories (EFTs) [3] and open quantum systems [4][5][6][7][8][9][10][11].
Ever more precise measurements of the heavy quarkonium nuclear suppression factor, the observable of interest in quantifying heavy quarkonium suppression, by the ALICE [12], ATLAS [13], and CMS [14,15] collaborations necessitate matchingly precise theoretical calculations.
To describe the evolution of heavy quarkonium in medium, it is advantageous to use EFTs that exploit the hierarchy of energy scales exhibited by nonrelativistic bound states.Integrating out modes associated with the heavy quark mass M from QCD leads to nonrelativistic QCD (NRQCD) [16,17].Further, integrating out soft modes associated with the momentum transfer scale given by M times the relative velocity v, which is also the scale of the inverse of the heavy quark-antiquark distance in the bound state, leads to potential NRQCD (pNRQCD) [18,19] and, at the Lagrangian level, to the emergence of quarkonium potentials.For weakly coupled quarkonia, the potentials may be computed in perturbation theory; at leading order in the nonrelativistic and coupling expansion, these are an attractive color-singlet and repulsive color-octet static potential.
Combining nonrelativistic EFTs with the open quantum system (OQS) framework, it is possible to derive a master equation for the evolution of the quarkonium density matrix [4,5,10].
A weak coupling treatment of the quarkonium soft modes is justified for quarkonia, whose typical radius is much smaller than the inverse of the hadronic scale Λ QCD .This is certainly the case for the lowest-lying bottomonium states and may be extended, although more marginally, to all bottomonium states below the open flavor threshold.At leading order, the weak-coupling static potential coincides with the Coulomb potential.It is well known that the Coulomb potential describes poorly the quarkonium spectrum [20].However, it has been shown that including higherorder corrections to the static potential, which is known up to three loops [21][22][23], and subtracting from the perturbative series the leading renormalon contribution by reabsorbing it into a renormalon subtracted mass [24] leads to a well-behaved perturbative series and a fair reproduction of the bottomonium spectrum below threshold (see, for instance, Refs.[25][26][27]).In this work, we use the three-loop renormalon subtracted version of the QCD static potential put forward in Refs.[27][28][29][30].We treat this potential, together with the kinetic energy, as the leading order contribution to the pNRQCD Hamiltonian and, therefore, to the Schrödinger equation describing the bound state.For an early application based on the same idea, see Ref. [31].We further supplement the perturbative static potential with its leading nonperturbative correction.This correction has been derived rigorously in pNRQCD and has the form of a harmonic potential, whose strength is given by the time integral of a suitable nonlocal chromoelectric correlator [19,32].The Lindblad equation describing the time evolution of the bottomonium density matrix in the QGP formed in heavy-ion collisions has been derived in pNRQCD in Refs.[4,5].The Lindblad equation depends on the static potential defined above and on heavy quarkonium transport co-efficients.The transport coefficients may be extracted from fitting available lattice data of the bottomonium thermal mass shifts and decay widths [33] with the corresponding quantities computed in pNRQCD.Finally, we solve the Lindblad equation in the computational setup of Ref. [34].
We work in the strict quantum Brownian limit, i.e., we assume the system correlation time to be much larger than the environment correlation time, which is proportional to the inverse of the temperature, and we neglect higher-order corrections.We find an accurate description of the LHC data on bottomonium suppression.
The remainder of the paper is structured as follows: Section II A briefly presents the theoretical background for describing heavy quarkonium suppression using EFTs and open quantum systems.
In section II B, we introduce the three-loop renormalon subtracted static potential and its leading nonperturbative correction, and in section II C, we discuss the quarkonium transport coefficients used in the simulation.We present our results for the nuclear modification factor in section III and conclude in section IV.

II. THEORY
A. Open quantum system formulation of quarkonium suppression Heavy quarkonium exhibits distinct hierarchically ordered energy scales, making it amenable to an EFT description.Using pNRQCD, in Refs.[4,5], a master equation was derived for the density matrix of a quarkonium state traveling through a strongly coupled, thermal QGP and satisfying the conditions where a 0 is the quarkonium Bohr radius, T the QGP temperature, and E the typical quarkonium energy.Because 1/E is the system correlation time, the strict Brownian limit corresponds to expanding in E/(πT ) and keeping only the leading term; beyond the strict Brownian limit, the Lindblad equation has been studied in Ref. [10] at order E/(πT ).In the strict Brownian limit, the master equation can be cast into a Lindblad form [35,36] with Hamiltonian and Lindblad operators The operators h s and h o encode the contributions from the energy modes that have been integrated out in constructing pNRQCD.These are, according to the hierarchy of Eq. ( 1), all the modes associated with energy scales larger than the thermal scales and Λ QCD .By construction, therefore, h s and h o are the in-vacuum, perturbative Hamiltonians of a heavy quark-antiquark pair in a color singlet and octet configuration, respectively, where r is the heavy quark-antiquark distance, p = −i∇ r the relative momentum, V pert The low-energy and thermal dynamics in the evolution equation are encoded in transport coefficients.In deriving Eqs. ( 3)-( 5), we have considered only the heavy quarkonium momentum diffusion coefficient κ and its dispersive counterpart, γ.The transport coefficients κ and γ are defined in terms of chromoelectric correlators as where Ω(t 2 , t 1 ) ab = P exp −ig ) is an adjoint Wilson line.The Lindblad equation provides the full three-dimensional evolution of the quarkonium density matrix.Utilizing the spherical symmetry of the system, we follow the procedure given in detail in Appendix C of Ref. [10] to bring the evolution equation into the form of a one-dimensional Lindblad equation describing the evolution of the radial wave function; this greatly reduces the computational cost of solving the Lindblad equation.The solution of the resulting equation is found by using the open-source code QTraj [34], which employs the Monte Carlo wave function method [37] to solve the Lindblad equation in a computationally efficient and embarrassingly parallel manner.

B. Heavy quarkonium potential
In previous works [6,10,11], the perturbative color-singlet and color-octet potentials were taken at leading order in α s , i.e., V pert α s , the strong coupling, evaluated at the inverse of the Bohr radius a 0 determined from the solution of the self-consistency equation Taking M b = 4.850 GeV in order to reproduce the Υ(1S) mass from the solution of the Schrödinger equation with V pert s ≈ V c s and solving Eq. ( 9) by running α s at one-loop accuracy using RunDec [38] from α s (M Z = (91.1876± 0.0021) GeV) = 0.1181 ± 0.0011, we find 1/a 0 = 1.11GeV.
To improve on this description, we consider here the case where V pert s and V pert o incorporate the three-loop expression.The three-loop color singlet potential has been computed in Refs.[22,23] and the three-loop color octet potential in Ref. [39].At three loops (3L), both the singlet and octet potential are infrared divergent with the divergence reabsorbed into nonperturbative corrections to the energy [19,21]; we call ν us the corresponding renormalization scale.Higher-loop corrections to the static potential are poorly convergent, a behavior that may be traced back to the order Λ QCD renormalon affecting the potential.The renormalon may be subtracted and reabsorbed into a redefinition of the heavy quark mass [40,41].We adopt the renormalon subtraction scheme RS ′ introduced in Ref. [24]; we call ν f the renormalon factorization scale.At very short distances, potentially large logarithms of the form log(νr), where ν is the renormalization scale of the strong coupling, are resummed in the running of α s by setting ν = 1/r.At larger distances, keeping the renormalization scale constant ensures the convergence of the perturbative expansion.We call ν r the scale separating the two regions, while still keeping both in the perturbative regime.The resulting potential has then the form [27-30] The explicit expression of the coefficients s,RS ′ can be found in Appendix A. The potential depends on the four scales ν, ν f , ν r , and ν us .We set ν f = ν us = 1 GeV and ν r = ν = 1/a 0 = 1.12 GeV, where 1/a 0 is calculated using Eq. ( 9) for M b = 4.921 GeV (for the choice of M b , see the discussion at the end of this subsection).We evaluate α s using RunDec [38] from α s (M Z ) and evaluate α k+1 s to (4 − k)-loop accuracy in Eq. (10).
In the scale setting of Eq. ( 1), the quarkonium potential also gets nonperturbative and thermal corrections.The leading one is encoded in the term proportional to γr 2 /2 in the Hamiltonian H given in Eq. ( 3). 1 The in-vacuum, nonperturbative (np) correction to the color singlet static potential, reads [19,32] The quantity γ(T = 0) stands for the in-vacuum, zero-temperature part of the chromoelectric correlator defined in Eq. ( 8).The fact that γ(T = 0) is different from zero should be contrasted with the fact that the quarkonium momentum diffusion coefficient κ vanishes in vacuum as a consequence of quarkonium decay into an unbound heavy quark-antiquark pair being kinematically forbidden in vacuum.An estimate of γ(T = 0) is provided in Eq. ( 18) of section II C. The full in vacuum color-singlet potential that we consider here is then Its perturbative part is the three-loop, renormalon subtracted, renormalization group improved color-singlet potential defined in Eq. (10).Its nonperturbative part is given in Eq. ( 11).
Similarly we define a three-loop, renormalon subtracted, renormalization group improved coloroctet potential [39,42] The coefficients o and the Coulomb potentials V c s , V c o is shown in Fig. 1.Since the expressions of the potentials in Eqs. ( 12) and ( 13) rely on assuming the soft modes to be perturbative, the curves are theoretically justified only at distances shorter than 1/Λ QCD , which we may indicatively take to be about 0.3 fm.The curves for V 3L o and V c o exhibit no qualitative difference, while V 3L+np s shows a much steeper slope than V c s at distances between 0.15 fm and 0.3 fm (and larger).The large, positive slope of the potential V 3L+np s at those distances is suggestive of an onset of confinement.
The bottomonium spectrum follows from solving the radial Schrödinger equation where u nl (r) is the reduced radial wave function of the state, E nl the corresponding binding energy  I. Since we neglect relativistic corrections to the potential, which are 1/M b suppressed, we do not expect to reproduce the spectrum exactly; nevertheless, we note that in contrast to the Coulomb potential, the potential V 3L+np s yields a spectrum accurate to within 1%.Moreover, the potential V 3L s has been instrumental in reproducing the fine structure of the bottomonium spectrum [27], electromagnetic decays [28] and electromagnetic transition widths [29,30].
The singlet and octet Coulomb potentials V c s and V c o (dashed lines) compared to the improved potentials V 3L+np s and V 3L o (solid lines).For the improved potentials, there is a discontinuity in the first derivative at 1/ν r = 0.197 fm.This can be seen by the kink in the solid blue curve, while it is barely noticeable in the solid orange curve that describes the octet potential.
In contrast to κ, γ contains vacuum contributions.To isolate the vacuum part, we write where γ(T = 0) denotes the vacuum contibution to γ while κ and γ are dimensionless coefficients of the thermal contributions.Up to now, there are no direct lattice QCD determinations of the quantities κ and γ.We note that lattice measurements of a related quantity, the heavy quark momentum diffusion coefficient, have existed for a number of years (for a recent determination see Ref. [46]) and anticipate lattice measurements of κ and γ in the near future.
To obtain an estimate for γ(T = 0), we parametrize the nonperturbative part of the chromoelectric correlator as where |0⟩ is the vacuum state, ⟨0| E 2 (0) |0⟩ the chromoelectric condensate, and Λ E the gluelump mass [19,32,47].This parameterization interpolates in the simplest possible way between the condensate at t = 0 and a large time behavior where the spectral decomposition of the correlator is dominated by the lowest-lying gluelump with the quantum numbers of the chromoelectric field.
Likewise, we extract values for κ by calculating the width Γ and fitting it to the lattice measurements of the in-medium width reported in Ref. [33].In our setting, the quarkonium width, which is entirely thermal, is given by [44] Since different b-masses were used in the lattice calculation [33], we rescale the reported lattice widths by a factor M lattice

III. RESULTS
We solve Eq. ( 2) using the open-source code QTraj [34] choosing for V pert s and V pert o the threeloop potentials V 3L s and V 3L o .For the medium evolution, we include input from 3+1D quasiparticle anisotropic hydrodynamics simulations [50].We initiate the coupling to the hydrodynamic evolution at τ med = 0.6 fm and calculate the evolution until a final temperature of T f = 250 MeV is reached, which we take as the lowest temperature for which the strict Brownian limit is valid [6].For more details on the numerical setup, we refer the reader to Refs.[10,11].Using QTraj, we sample approximately 140 000 physical trajectories from the hydrodynamics simulation and 30 quantum trajectories per physical trajectory.At the end of the evolution, we project the wave function of the final state onto the bottomonium eigenstates for Υ(1S), Υ(2S), Υ(3S) and χ b (1P ), χ b (2P ), which are relevant through feed-down, obtained from solving the Schrödinger equation for the improved potential V 3L+np s .By taking the ratio of the final state overlap to the initial state overlap, we obtain the survival probability of the corresponding state.Including the effect of feed-down as described in Refs.[6,10], we finally obtain predictions for the bottomonium nuclear modification factor R AA , which is the bottomonium yield in heavy-ion collisions normalized by the bottomonium yield in proton-proton collisions.The results for R AA in terms of the number of participating nucleons, N part , obtained from using the three-loop potentials V 3L s and V 3L o in the Lindblad equation, are shown in Figure 3.The colored bands correspond to the statistical uncertainties and the uncertainties from κ.The dash-dotted and dashed lines indicate the curves of κ = 1.72 and κ = 2.04, respectively.The data are taken from the ALICE [12], ATLAS [13], and CMS [14,15] experiments.
For the numerical computation, we use the transport coefficient values obtained in section II C, i.e., γ = 0, γ(T = 0) = 0.017 GeV 3 and κ = 1.88 ± 0.16.The bands indicate the uncertainty in κ and the statistical uncertainty from the quantum trajectory algorithm.The dash-dotted and dashed lines are the curves for κ = 1.72 and κ = 2.04, respectively.To assess the compatibility with the data, we calculate the weighted relative mean square error (WMSE) where the sum runs over the discrete values of N part corresponding to the different measurements by the experiments, R exp AA,i denotes the ith measurement, R QTraj AA,i the corresponding prediction from QTraj, and σ exp i and σ QTraj i the relative uncertainties.The data are measurements by the AL-ICE [12], ATLAS [13], and CMS [14,15] collaborations.The QTraj predictions and corresponding uncertainties at N part given by the experiments are obtained by linear extrapolation between QTraj prediction points.The WMSE for the three-loop potentials with the above parameters is ∆ 2 = 2.1 × 10 −3 .In contrast, results obtained by solving the Lindblad equation with V pert s and V pert o set to the Coulomb potentials V c s and V c o with γ = 0, γ(T = 0) = 0 and κ = 0.33 and projecting onto bottomonium eigenstates obtained by solving the Schrödinger equation for the Coulomb potential V c s give ∆ 2 = 0.22, showing that the three-loop potential is better suited to describe the data than the simple Coulomb potential.We note that the wave functions calculated with V c s and V 3L+np s , especially the wave functions of the 2S and 3S states, differ considerably due to the steeper slope of the V 3L+np s potential at large distances.
Finally, we remark that the results with the Coulomb potential of Refs.[6,51] were obtained using different parameters for the b-mass and α s than here and, importantly, different values for the transport coefficients.In this work, however, we constrain the transport coefficients to reproduce the thermal widths of the S-and P -wave states recently computed in lattice QCD.Under this constraint, the three-loop potential provides a better description of the bottomonium nuclear modification factor data from the LHC experiments.

IV. CONCLUSIONS
In the paper, we compute bottomonium suppression in heavy-ion collisions from an evolution equation that encompasses the three-loop quarkonium potential and the vacuum part of the dispersive transport coefficient γ(T = 0).There are a number of advantages to this approach.First, the three-loop potential plus the addition of a small nonperturbative correction derived from the operator product expansion (with the bottom mass fixed on the 1S state) provides a realistic description of the spin-averaged bottomonium spectrum; see Table I.Second, the thermal decay width that follows from the bottomonium wave functions is consistent with that computed in lattice QCD; this allows to fit reasonably well the quarkonium momentum diffusion coefficient on the lattice data; see Fig. 2. Finally, the obtained bottomonium nuclear modification factor shown in Fig. 3 provides a more accurate description of the LHC data in the strict Brownian limit than previous analyses with the Coulomb potential.In the future, this study could be extended beyond the Brownian limit by including E/(πT ) corrections, as done for the Coulomb potential in Ref. [10].
and V pert o the color octet static potential computed in perturbative QCD.The number of colors is N c = 3, and, for further use, the Casimir of the adjoint representation is C A = 3, the Casimir of the fundamental representation is C F = 4/3 and the color matrix normalization is

(
hence the mass of the state is 2M b + E nl ), and V s the color singlet static potential.The b-mass M b is fixed by fitting the 1S mass to the PDG value M (1S) = 9.445 GeV.The bottomonium masses are understood as spin averaged.The b-mass that follows from solving the Schrödinger equation for V s = V c s is M b = 4.850 GeV, whereas the one that follows from solving the Schrödinger equation for V s = V 3L+np s , with γ(T = 0) given in Eq. (18), is M b = 4.921 GeV.The results for the spectrum are shown in Table 1S)/ GeV 9.445 9.445 9.445 M (2S)/ GeV 10.017 9.635 10.042 M (3S)/ GeV 10.355 9.670 10.395 M (1P )/ GeV 9.888 9.635 9.887 M (2P )/ GeV 10.251 9.670 10.279 TABLE I: The spin-averaged bottomonium spectrum from PDG data [43], from solving the Schrödinger equation with the Coulomb potential V c s , and from solving the Schrödinger equation with the potential V 3L+np s and γ(T = 0) as in Eq. (18).The b-mass is fixed by fitting the 1S prediction to the PDG value yielding M b = 4.850 GeV for V c s and M b = 4.921 GeV for V 3L+np s .

b /M b 2 ,FIG. 2 :
FIG.2:The calculated widths Γ for the Coulomb potential V c s and the improved potential V 3L+np s

FIG. 3 :
FIG. 3: The nuclear modification factor R AA for the 1S, 2S, and 3S bottomonium states as a function of the number of participating nucleons N part .The colored solid lines show results obtained with κ = 1.88, γ(T = 0) = 0.017 GeV 3 and γ = 0 using the three-loop potentials V 3L s and V 3L o in the evolution equation.