# Bottomonium suppression from the three-loop QCD potential

Nora Brambilla<EMAIL_ADDRESS>Technical University of Munich, TUM
School of Natural Sciences, Physics Department, James-Franck-Strasse 1, 85748
Garching, Germany  Institute for Advanced Study, Technische Universität
München, Lichtenbergstrasse 2 a, 85748 Garching, Germany  Munich Data Science
Institute, Technische Universität München, Walther-von-Dyck-Strasse 10, 85748
Garching, Germany  Tom Magorsch<EMAIL_ADDRESS>Technical University of
Munich, TUM School of Natural Sciences, Physics Department, James-Franck-
Strasse 1, 85748 Garching, Germany  Michael Strickland<EMAIL_ADDRESS>Department of Physics, Kent State University, Kent, Ohio 44242, United States
Antonio Vairo<EMAIL_ADDRESS>Technical University of Munich, TUM School
of Natural Sciences, Physics Department, James-Franck-Strasse 1, 85748
Garching, Germany  Peter Vander Griend<EMAIL_ADDRESS>Department
of Physics and Astronomy, University of Kentucky, Lexington, Kentucky 40506,
United States  Theoretical Physics Department, Fermilab, P.O. Box 500,
Batavia, Illinois 60510, United States

###### Abstract

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.

††preprint: TUM-EFT 186/23, FERMILAB-PUB-24-0029-V

##  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/\psi$ mesons traveling through the QGP
produced in heavy-ion collisions should lead to a lower number of particles
measured compared to proton-proton collisions Matsui and Satz (1986). In the
past, this effect has been studied extensively for bottomonium Andronic  _et
al._ (2024), where due to the heavy bottom quark mass, various calculational
simplifications are possible, including the use of effective field theories
(EFTs) Brambilla  _et al._ (2005) and open quantum systems Brambilla  _et al._
(2017, 2018, 2021a); Akamatsu (2022); Yao (2021); Miura  _et al._ (2022);
Brambilla  _et al._ (2022, 2023). Ever more precise measurements of the heavy
quarkonium nuclear suppression factor, the observable of interest in
quantifying heavy quarkonium suppression, by the ALICE Acharya  _et al._
(2021), ATLAS Aad  _et al._ (2023), and CMS Sirunyan  _et al._ (2019);
Tumasyan  _et al._ (2023) 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) Caswell and
Lepage (1986); Bodwin  _et al._ (1995). 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)
Pineda and Soto (1998); Brambilla  _et al._ (2000) 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 Brambilla  _et al._ (2017, 2018, 2022).

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 $\Lambda_{\textrm{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
Brambilla and Vairo (1999). However, it has been shown that including higher-
order corrections to the static potential, which is known up to three loops
Brambilla  _et al._ (1999); Anzai  _et al._ (2010); Smirnov  _et al._ (2010),
and subtracting from the perturbative series the leading renormalon
contribution by reabsorbing it into a renormalon subtracted mass Pineda (2001)
leads to a well-behaved perturbative series and a fair reproduction of the
bottomonium spectrum below threshold (see, for instance, Refs. Brambilla  _et
al._ (2001); Sumino (2016); Peset  _et al._ (2018)). In this work, we use the
three-loop renormalon subtracted version of the QCD static potential put
forward in Refs. Kiyo  _et al._ (2010); Pineda and Segovia (2013); Peset  _et
al._ (2018); Segovia  _et al._ (2019). 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. Titard and
Yndurain (1994). 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
Brambilla  _et al._ (2000); Brambilla and Vairo (2000).

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. Brambilla  _et al._ (2017, 2018). The Lindblad equation depends on the
static potential defined above and on heavy quarkonium transport coefficients.
The transport coefficients may be extracted from fitting available lattice
data of the bottomonium thermal mass shifts and decay widths Larsen  _et al._
(2019) with the corresponding quantities computed in pNRQCD. Finally, we solve
the Lindblad equation in the computational setup of Ref. Omar  _et al._
(2022). 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.1 briefly
presents the theoretical background for describing heavy quarkonium
suppression using EFTs and open quantum systems. In section II.2, we introduce
the three-loop renormalon subtracted static potential and its leading
nonperturbative correction, and in section II.3, 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

###  II.1 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. Brambilla
_et al._ (2017, 2018), a master equation was derived for the density matrix of
a quarkonium state traveling through a strongly coupled, thermal QGP and
satisfying the conditions

$M\gg 1/a_{0}\gg\pi T\sim gT,\Lambda_{\textrm{QCD}}\gg E\,,$ (1)

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/(\pi T)$ and keeping
only the leading term; beyond the strict Brownian limit, the Lindblad equation
has been studied in Ref. Brambilla  _et al._ (2022) at order $E/(\pi T)$. In
the strict Brownian limit, the master equation can be cast into a Lindblad
form Lindblad (1976); Gorini  _et al._ (1976)

$\frac{d\rho(t)}{dt}=-i[H,\rho(t)]+\sum_{n}\left[C^{n}_{i}\rho(t)C^{n\dagger}_{i}-\frac{1}{2}\left\\{C^{n\dagger}_{i}C^{n}_{i},\rho(t)\right\\}\right],$
(2)

with Hamiltonian

$H=\begin{pmatrix}h_{s}+\displaystyle\frac{r^{2}}{2}\gamma&0\\\
0&\displaystyle
h_{o}+\frac{N^{2}_{c}-2}{2(N^{2}_{c}-1)}\frac{r^{2}}{2}\gamma\end{pmatrix},$
(3)

and Lindblad operators

$\displaystyle C^{0}_{i}$
$\displaystyle=\sqrt{\frac{\kappa}{N^{2}_{c}-1}}r_{i}\begin{pmatrix}0&1\\\
\sqrt{N^{2}_{c}-1}&0\end{pmatrix},$ (4) $\displaystyle C^{1}_{i}$
$\displaystyle=\sqrt{\frac{\kappa(N^{2}_{c}-4)}{2(N^{2}_{c}-1)}}r_{i}\begin{pmatrix}0&0\\\
0&1\end{pmatrix}.$ (5)

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 $\Lambda_{\textrm{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,

$h_{s}=\frac{\bm{p}^{2}}{M}+V_{s}^{\textrm{pert}}(r),\qquad\qquad
h_{o}=\frac{\bm{p}^{2}}{M}+V_{o}^{\textrm{pert}}(r),$ (6)

where $\bm{r}$ is the heavy quark-antiquark distance,
$\bm{p}=-i\bm{\nabla}_{r}$ the relative momentum, $V_{s}^{\textrm{pert}}$ the
color singlet and $V_{o}^{\textrm{pert}}$ 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 $T_{F}=1/2$.

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 $\kappa$ and its dispersive
counterpart, $\gamma$. The transport coefficients $\kappa$ and $\gamma$ are
defined in terms of chromoelectric correlators as

$\displaystyle\kappa$
$\displaystyle=\frac{g^{2}}{3N_{c}}\text{Re}\int^{\infty}_{0}dt\,\braket{E^{a\,i}(t)\Omega(t,0)^{ab}E^{b\,i}(0)},$
(7) $\displaystyle\gamma$
$\displaystyle=\frac{g^{2}}{3N_{c}}\text{Im}\int^{\infty}_{0}dt\,\braket{E^{a\,i}(t)\Omega(t,0)^{ab}E^{b\,i}(0)},$
(8)

where
$\displaystyle\Omega(t_{2},t_{1})^{ab}=P\exp\left(-ig\int^{t_{2}}_{t_{1}}dtA^{c}_{0}(t)(-if^{abc})\right)$
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. Brambilla  _et al._
(2022) 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 Omar  _et al._ (2022), which employs the Monte Carlo wave function
method Mølmer  _et al._ (1993) to solve the Lindblad equation in a
computationally efficient and embarrassingly parallel manner.

###  II.2 Heavy quarkonium potential

In previous works Brambilla  _et al._ (2021a, 2022, 2023), the perturbative
color-singlet and color-octet potentials were taken at leading order in
$\alpha_{s}$, i.e., $V_{s}^{\textrm{pert}}(r)\approx
V_{s}^{c}(r)\equiv-4\,\alpha_{s}/(3r)$ and $V_{o}^{\textrm{pert}}(r)\approx
V_{o}^{c}(r)\equiv\alpha_{s}/(6r)$, with $\alpha_{s}$, the strong coupling,
evaluated at the inverse of the Bohr radius $a_{0}$ determined from the
solution of the self-consistency equation

$a_{0}=\frac{3}{2\,M_{b}\,\alpha_{s}(1/a_{0})}\,.$ (9)

Taking $M_{b}=$4.850\text{\,}\mathrm{GeV}$$ in order to reproduce the
$\Upsilon(1S)$ mass from the solution of the Schrödinger equation with
$V_{s}^{\textrm{pert}}\approx V_{s}^{c}$ and solving Eq. (9) by running
$\alpha_{s}$ at one-loop accuracy using RunDec Herren and Steinhauser (2018)
from $\alpha_{s}(M_{Z}=$91.1876\pm 0.0021\text{\,}\mathrm{GeV}$)=$0.1181\pm
0.0011$$, we find $1/a_{0}=$1.11\text{\,}\mathrm{GeV}$$.

To improve on this description, we consider here the case where
$V_{s}^{\textrm{pert}}$ and $V_{o}^{\textrm{pert}}$ incorporate the three-loop
expression. The three-loop color singlet potential has been computed in Refs.
Anzai  _et al._ (2010); Smirnov  _et al._ (2010) and the three-loop color
octet potential in Ref. Anzai  _et al._ (2013). At three loops (3L), both the
singlet and octet potential are infrared divergent with the divergence
reabsorbed into nonperturbative corrections to the energy Brambilla  _et al._
(1999, 2000); we call $\nu_{\textrm{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 $\Lambda_{\textrm{QCD}}$
renormalon affecting the potential. The renormalon may be subtracted and
reabsorbed into a redefinition of the heavy quark mass Pineda (1998); Hoang
_et al._ (1999). We adopt the renormalon subtraction scheme RS′ introduced in
Ref. Pineda (2001); we call $\nu_{f}$ the renormalon factorization scale. At
very short distances, potentially large logarithms of the form $\log(\nu r)$,
where $\nu$ is the renormalization scale of the strong coupling, are resummed
in the running of $\alpha_{s}$ by setting $\nu=1/r$. At larger distances,
keeping the renormalization scale constant ensures the convergence of the
perturbative expansion. We call $\nu_{r}$ the scale separating the two
regions, while still keeping both in the perturbative regime. The resulting
potential has then the form Kiyo  _et al._ (2010); Pineda and Segovia (2013);
Peset  _et al._ (2018); Segovia  _et al._ (2019)

$V_{s}^{\textrm{pert}}(\nu,\nu_{f},\nu_{r},\nu_{\textrm{us}},r)\approx
V_{s}^{\textrm{3L}}(\nu,\nu_{f},\nu_{r},\nu_{\textrm{us}},r)=\begin{cases}\sum^{3}_{k=0}V^{(k)}_{s,\text{RS}^{\prime}}\alpha^{k+1}_{s}(1/r)\quad\text{if}\;\,r<\nu^{-1}_{r}\\\
\sum^{3}_{k=0}V^{(k)}_{s,\text{RS}^{\prime}}\alpha^{k+1}_{s}(\nu)\qquad\text{if}\;r>\nu^{-1}_{r}\end{cases}.$
(10)

The explicit expression of the coefficients $V^{(k)}_{s,\text{RS}^{\prime}}$
can be found in Appendix A. The potential depends on the four scales $\nu$,
$\nu_{f}$, $\nu_{r}$, and $\nu_{\textrm{us}}$. We set
$\nu_{f}=\nu_{\textrm{us}}=$1\text{\,}\mathrm{GeV}$$ and
$\nu_{r}=\nu=1/a_{0}=$1.12\text{\,}\mathrm{GeV}$$, where $1/a_{0}$ is
calculated using Eq. (9) for $M_{b}=$4.921\text{\,}\mathrm{GeV}$$ (for the
choice of $M_{b}$, see the discussion at the end of this subsection). We
evaluate $\alpha_{s}$ using RunDec Herren and Steinhauser (2018) from
$\alpha_{s}(M_{Z})$ and evaluate $\alpha^{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 $\gamma r^{2}/2$ in the Hamiltonian $H$ given in Eq.
(3).111 Higher order terms not considered here ensure the cancellation of the
$\nu_{\textrm{us}}$ dependence of the perturbative potential Brambilla  _et
al._ (2000). The in-vacuum, nonperturbative (np) correction to the color
singlet static potential, reads Brambilla  _et al._ (2000); Brambilla and
Vairo (2000)

$V^{\text{np}}_{s}(r)=\frac{r^{2}}{2}\,\gamma(T=0)\,.$ (11)

The quantity $\gamma(T=0)$ stands for the in-vacuum, zero-temperature part of
the chromoelectric correlator defined in Eq. (8). The fact that $\gamma(T=0)$
is different from zero should be contrasted with the fact that the quarkonium
momentum diffusion coefficient $\kappa$ vanishes in vacuum as a consequence of
quarkonium decay into an unbound heavy quark-antiquark pair being
kinematically forbidden in vacuum. An estimate of $\gamma(T=0)$ is provided in
Eq. (18) of section II.3. The full in vacuum color-singlet potential that we
consider here is then

$V^{\text{3L+np}}_{s}(r)=V^{\text{3L}}_{s}(\nu,\nu_{f},\nu_{r},\nu_{\textrm{us}},r)+V^{\text{np}}_{s}(r)\,.$
(12)

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 color-octet potential Bali and Pineda (2004); Anzai  _et al._ (2013)

$V_{o}^{\textrm{pert}}(\nu,\nu_{f},\nu_{r},\nu_{\textrm{us}},r)\approx
V_{o}^{\textrm{3L}}(\nu,\nu_{f},\nu_{r},\nu_{\textrm{us}},r)=\begin{cases}\sum^{3}_{k=0}V^{(k)}_{o,\text{RS}^{\prime}}\alpha^{k+1}_{s}(1/r)\quad\text{if}\;\,r<\nu^{-1}_{r}\\\
\sum^{3}_{k=0}V^{(k)}_{o,\text{RS}^{\prime}}\alpha^{k+1}_{s}(\nu)\qquad\text{if}\;r>\nu^{-1}_{r}\end{cases}.$
(13)

The coefficients $V^{(k)}_{o,\text{RS}^{\prime}}$ are listed in Appendix A.

A comparison between the potentials $V^{\text{3L+np}}_{s}$,
$V^{\text{3L}}_{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/\Lambda_{\textrm{QCD}}$, which we may indicatively take to be about 0.3 fm.
The curves for $V^{\text{3L}}_{o}$ and $V^{c}_{o}$ exhibit no qualitative
difference, while $V^{\text{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^{\text{3L+np}}_{s}$ at those distances is
suggestive of an onset of confinement.

Figure 1: The singlet and octet Coulomb potentials $V^{c}_{s}$ and $V^{c}_{o}$
(dashed lines) compared to the improved potentials $V^{\text{3L+np}}_{s}$ and
$V^{\text{3L}}_{o}$ (solid lines). For the improved potentials, there is a
discontinuity in the first derivative at
$1/\nu_{r}=$0.197\text{\,}\mathrm{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.

The bottomonium spectrum follows from solving the radial Schrödinger equation

$-\frac{1}{M_{b}}\frac{d^{2}u_{nl}(r)}{dr^{2}}+\left(V_{s}(r)+\frac{l(l+1)}{M_{b}r^{2}}u_{nl}(r)\right)=E_{nl}u_{nl}(r)\,,$
(14)

where $u_{nl}(r)$ is the reduced radial wave function of the state, $E_{nl}$
the corresponding binding energy (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\text{\,}\mathrm{GeV}$$. The bottomonium masses are understood as
spin averaged. The $b$-mass that follows from solving the Schrödinger equation
for $V_{s}=V_{s}^{c}$ is $M_{b}=$4.850\text{\,}\mathrm{GeV}$$, whereas the one
that follows from solving the Schrödinger equation for
$V_{s}=V^{\text{3L+np}}_{s}$, with $\gamma(T=0)$ given in Eq. (18), is
$M_{b}=$4.921\text{\,}\mathrm{GeV}$$. The results for the spectrum are shown
in Table 1. 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^{\text{3L+np}}_{s}$ yields a spectrum accurate to within 1%. Moreover, the
potential $V^{\text{3L}}_{s}$ has been instrumental in reproducing the fine
structure of the bottomonium spectrum Peset  _et al._ (2018), electromagnetic
decays Kiyo  _et al._ (2010) and electromagnetic transition widths Pineda and
Segovia (2013); Segovia  _et al._ (2019).

| PDG |  $V^{c}_{s}$ | $V^{\text{3L+np}}_{s}$
---|---|---|---
$M(1S)/\,$\mathrm{GeV}$$ |  $9.445$ |  $9.445$ |  $9.445$
$M(2S)/\,$\mathrm{GeV}$$ |  $10.017$ |  $9.635$ |  $10.042$
$M(3S)/\,$\mathrm{GeV}$$ |  $10.355$ |  $9.670$ |  $10.395$
$M(1P)/\,$\mathrm{GeV}$$ |  $9.888$ |  $9.635$ |  $9.887$
$M(2P)/\,$\mathrm{GeV}$$ |  $10.251$ |  $9.670$ |  $10.279$
Table 1: The spin-averaged bottomonium spectrum from PDG data Zyla  _et al._
(2020), from solving the Schrödinger equation with the Coulomb potential
$V^{\text{c}}_{s}$, and from solving the Schrödinger equation with the
potential $V^{\text{3L+np}}_{s}$ and $\gamma(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\text{\,}\mathrm{GeV}$$ for $V^{\text{c}}_{s}$ and
$M_{b}=$4.921\text{\,}\mathrm{GeV}$$ for $V^{\text{3L+np}}_{s}$.

###  II.3 Transport coefficients

According to Eqs. (2)-(5), the interaction between heavy quark-antiquark pairs
and the QGP is encoded in the transport coefficients $\kappa$ and $\gamma$,
defined in Eq. (7) and (8), respectively Brambilla  _et al._ (2017, 2018,
2019); Eller  _et al._ (2019). In contrast to $\kappa$, $\gamma$ contains
vacuum contributions. To isolate the vacuum part, we write

$\displaystyle\kappa$ $\displaystyle=\hat{\kappa}T^{3}\,,$ (15)
$\displaystyle\gamma$ $\displaystyle=\gamma(T=0)+\hat{\gamma}T^{3}\,,$ (16)

where $\gamma(T=0)$ denotes the vacuum contibution to $\gamma$ while
$\hat{\kappa}$ and $\hat{\gamma}$ are dimensionless coefficients of the
thermal contributions. Up to now, there are no direct lattice QCD
determinations of the quantities $\kappa$ and $\gamma$. 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. Brambilla  _et al._ (2020)) and anticipate lattice measurements of
$\kappa$ and $\gamma$ in the near future.

To obtain an estimate for $\gamma(T=0)$, we parametrize the nonperturbative
part of the chromoelectric correlator as

$\bra{0}E^{a\,i}(t)\Omega(t,0)^{ab}E^{b\,i}(0)\ket{0}=\bra{0}E^{2}(0)\ket{0}\mathrm{e}^{-i\Lambda_{E}t}\,,$
(17)

where $\ket{0}$ is the vacuum state, $\bra{0}E^{2}(0)\ket{0}$ the
chromoelectric condensate, and $\Lambda_{E}$ the gluelump mass Foster and
Michael (1999); Brambilla  _et al._ (2000); Brambilla and Vairo (2000). 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. The perturbative part of the correlator
is scaleless and proportional to $1/t^{4}$ and, therefore, vanishes in
$\gamma(T=0)$ when regularized in dimensional regularization. We take the
gluon condensate to be $\bra{0}\alpha_{s}G^{2}\ket{0}=$(6.4\pm
0.3)\text{\times}{10}^{-2}\text{\,}{\mathrm{GeV}}^{4}$$ Narison (2023), which
leads to $\bra{0}g^{2}E^{2}(0)\ket{0}=$-0.2\text{\,}{\mathrm{GeV}}^{4}$$, and
the chromoelectric gluelump mass $\Lambda_{E}=$1.25\text{\,}\mathrm{GeV}$$
Bali and Pineda (2004); Herr  _et al._ (2024). Inserting Eq. (17) into Eq.
(8), performing the integral, and regularizing at infinity, we obtain

$\gamma(T=0)=-\frac{\bra{0}g^{2}E^{2}(0)\ket{0}}{3N_{c}\Lambda_{E}}=$0.017\text{\,}{\mathrm{GeV}}^{3}$\,.$
(18)

The thermal part $\hat{\gamma}$ can be extracted from lattice measurements of
the in-medium heavy quarkonium mass shift Larsen  _et al._ (2019). In our
setting, the thermal mass shift is given by Brambilla  _et al._ (2019)

$\delta
M=\frac{\hat{\gamma}\,T^{3}}{2}\,\int_{0}^{\infty}dr\,u_{nl}^{*}r^{2}\,u_{nl}\,,$
(19)

where $u_{nl}$ is the solution of the Schrödinger equation (14). Since the
results reported in Ref. Larsen  _et al._ (2019) suggest a vanishing thermal
mass shift, we set $\hat{\gamma}=0$. Note that the vacuum part $\gamma(T=0)$
is not accessible through the thermal mass shift.

Figure 2: The calculated widths $\Gamma$ for the Coulomb potential $V_{s}^{c}$
and the improved potential $V^{\text{3L+np}}_{s}$ compared to the lattice data
from Ref. Larsen  _et al._ (2019). The left panel shows calculations for $l=0$
states compared to $\Upsilon$ data, and the right panel shows calculations for
$l=1$ states compared to $\chi_{b}$ lattice data. In both cases, we use the
extracted and averaged values $\hat{\kappa}=$0.33$\pm 0.04$ for the Coulomb
potential and $\hat{\kappa}=$1.88$\pm 0.16$ for the potential
$V^{\text{3L+np}}_{s}$. The lattice widths displayed are rescaled by
$\left({M^{\text{lattice}}_{b}}/{M_{b}}\right)^{2}$ with
$M_{b}=$4.921\text{\,}\mathrm{GeV}$$.

Likewise, we extract values for $\hat{\kappa}$ by calculating the width
$\Gamma$ and fitting it to the lattice measurements of the in-medium width
reported in Ref. Larsen  _et al._ (2019). In our setting, the quarkonium
width, which is entirely thermal, is given by Brambilla  _et al._ (2019)

$\Gamma=\hat{\kappa}\,T^{3}\,\int_{0}^{\infty}dr\,u_{nl}^{*}r^{2}\,u_{nl}\,.$
(20)

Since different $b$-masses were used in the lattice calculation Larsen  _et
al._ (2019), we rescale the reported lattice widths by a factor
$\left({M^{\text{lattice}}_{b}}/{M_{b}}\right)^{2}$, where
$M^{\text{lattice}}_{b}$ denotes the $b$-mass used in the lattice computation
(see table II of Ref. Larsen  _et al._ (2019)), and $M_{b}$ is the $b$-mass
used in our computation. In practice, we rescale using
$M_{b}=$4.921\text{\,}\mathrm{GeV}$$, since rescaling with the $b$-mass
determined for the Coulomb potential, $M_{b}=$4.850\text{\,}\mathrm{GeV}$$,
leads to compatible widths within the error bars. Due to the Brownian limit,
the width in Eq. (20) may not be accurate for the lowest temperatures; we have
checked, however, that excluding the lattice data for
$T<$200\text{\,}\mathrm{MeV}$$ in the fit does not change qualitatively our
results. In the case of the potential $V^{\text{3L+np}}_{s}$, a weighted least
squares fit of $\hat{\kappa}$ from the $\Upsilon$ data gives
$\hat{\kappa}_{\Upsilon}=$2.08\pm 0.20$$, and, from the $\chi_{b1}$ data,
$\hat{\kappa}_{\chi}=$1.52\pm 0.26$$. The weighted average of these two
determinations is $\hat{\kappa}=$1.88\pm 0.16$$. Applying the same procedure
for the Coulomb potential, we find $\hat{\kappa}_{\Upsilon}=$1.16\pm 0.11$$
and $\hat{\kappa}_{\chi}=$0.23\pm 0.04$$ leading to the average
$\hat{\kappa}=$0.33\pm 0.04$$. We plot the temperature dependence of the width
for the lowest-lying states with $l=0$ and $l=1$ in the left and right panels
of Fig. 2, respectively. The plots show with an orange line the width for the
final average $\hat{\kappa}=$0.33$\pm 0.04$, obtained from the Coulomb
potential, and with a blue line the width for the final average
$\hat{\kappa}=$1.88\pm 0.16$$, obtained from the improved potential
$V^{\text{3L+np}}_{s}$. We see that the Coulomb fit is unable to describe both
the $\Upsilon$ and $\chi_{b}$ data at the same time, while the fit based on
the improved potential $V^{\text{3L+np}}_{s}$ provides a more reasonable
description of the lattice data in both cases.

##  III Results

We solve Eq. (2) using the open-source code QTraj Omar  _et al._ (2022)
choosing for $V_{s}^{\textrm{pert}}$ and $V_{o}^{\textrm{pert}}$ the three-
loop potentials $V^{\text{3L}}_{s}$ and $V^{\text{3L}}_{o}$. For the medium
evolution, we include input from 3+1D quasiparticle anisotropic hydrodynamics
simulations Alqahtani and Strickland (2021). We initiate the coupling to the
hydrodynamic evolution at $\tau_{\text{med}}=$0.6\text{\,}\mathrm{fm}$$ and
calculate the evolution until a final temperature of
$T_{f}=$250\text{\,}\mathrm{MeV}$$ is reached, which we take as the lowest
temperature for which the strict Brownian limit is valid Brambilla  _et al._
(2021a). For more details on the numerical setup, we refer the reader to Refs.
Brambilla  _et al._ (2022, 2023). 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 $\Upsilon(1S)$, $\Upsilon(2S)$, $\Upsilon(3S)$ and $\chi_{b}(1P)$,
$\chi_{b}(2P)$, which are relevant through feed-down, obtained from solving
the Schrödinger equation for the improved potential $V^{\text{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. Brambilla  _et al._ (2021a, 2022),
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_{\text{part}}$, obtained
from using the three-loop potentials $V^{\text{3L}}_{s}$ and
$V^{\text{3L}}_{o}$ in the Lindblad equation, are shown in Figure 3.

Figure 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_{\text{part}}$. The colored solid lines show results obtained with
$\hat{\kappa}=$1.88$$, $\gamma(T=0)=$0.017\text{\,}{\mathrm{GeV}}^{3}$$ and
$\hat{\gamma}=0$ using the three-loop potentials $V^{\text{3L}}_{\text{s}}$
and $V^{\text{3L}}_{\text{o}}$ in the evolution equation. The colored bands
correspond to the statistical uncertainties and the uncertainties from
$\hat{\kappa}$. The dash-dotted and dashed lines indicate the curves of
$\hat{\kappa}=$1.72$$ and $\hat{\kappa}=$2.04$$, respectively. The data are
taken from the ALICE Acharya  _et al._ (2021), ATLAS Aad  _et al._ (2023), and
CMS Sirunyan  _et al._ (2019); Tumasyan  _et al._ (2023) experiments.

For the numerical computation, we use the transport coefficient values
obtained in section II.3, i.e., $\hat{\gamma}=0$,
$\gamma(T=0)=$0.017\text{\,}{\mathrm{GeV}}^{3}$$ and $\hat{\kappa}=$1.88\pm
0.16$$. The bands indicate the uncertainty in $\hat{\kappa}$ and the
statistical uncertainty from the quantum trajectory algorithm. The dash-dotted
and dashed lines are the curves for $\hat{\kappa}=$1.72$$ and
$\hat{\kappa}=$2.04$$, respectively. To assess the compatibility with the
data, we calculate the weighted relative mean square error (WMSE)

$\Delta^{2}=\frac{1}{N_{\text{data}}}\sum_{i}\frac{\left(\frac{R^{\text{exp}}_{AA,i}-R^{\text{QTraj}}_{AA,i}}{R^{\text{exp}}_{AA,i}}\right)^{2}}{\left(\sigma^{\text{exp}}_{i}\big{/}R^{\text{exp}}_{AA,i}\right)^{2}+\left(\sigma^{\text{QTraj}}_{i}\big{/}R^{\text{QTraj}}_{AA,i}\right)^{2}}\,,$
(21)

where the sum runs over the discrete values of $N_{\rm part}$ corresponding to
the different measurements by the experiments, $R^{\text{exp}}_{AA,i}$ denotes
the $i$th measurement, $R^{\text{QTraj}}_{AA,i}$ the corresponding prediction
from QTraj, and $\sigma^{\text{exp}}_{i}$ and $\sigma^{\text{QTraj}}_{i}$ the
relative uncertainties. The data are measurements by the ALICE Acharya  _et
al._ (2021), ATLAS Aad  _et al._ (2023), and CMS Sirunyan  _et al._ (2019);
Tumasyan  _et al._ (2023) collaborations. The QTraj predictions and
corresponding uncertainties at $N_{\text{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
$\Delta^{2}=$2.1\text{\times}{10}^{-3}$$. In contrast, results obtained by
solving the Lindblad equation with $V_{s}^{\textrm{pert}}$ and
$V_{o}^{\textrm{pert}}$ set to the Coulomb potentials $V^{c}_{s}$ and
$V^{c}_{o}$ with $\hat{\gamma}=0$, $\gamma(T=0)=0$ and $\hat{\kappa}=$0.33$$
and projecting onto bottomonium eigenstates obtained by solving the
Schrödinger equation for the Coulomb potential $V^{c}_{s}$ give
$\Delta^{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_{s}^{c}$ and $V^{\text{3L+np}}_{s}$, especially
the wave functions of the $2S$ and $3S$ states, differ considerably due to the
steeper slope of the $V^{\text{3L+np}}_{s}$ potential at large distances.

Finally, we remark that the results with the Coulomb potential of Refs.
Brambilla  _et al._ (2021a, b) were obtained using different parameters for
the $b$-mass and $\alpha_{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 $\gamma(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 1. 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/(\pi T)$ corrections, as done for the Coulomb
potential in Ref. Brambilla  _et al._ (2022).

## Acknowledgments

N.B., T.M., and A.V. acknowledge support by the DFG cluster of excellence
ORIGINS funded by the Deutsche Forschungsgemeinschaft (DFG) under Germany’s
Excellence Strategy - EXC-2094-390783311. The work of N.B. is supported by the
DFG Grant No. BR 4058/5-1 “Open Quantum Sytems and Effective Field Theories
for hard probes of hot and/or dense medium”. N.B., T.M., and A.V. acknowledge
support from STRONG-2020, European Union’s Horizon 2020 research and
innovation program under grant agreement No. 824093. M.S. was supported by the
U.S. Department of Energy, Office of Science, Office of Nuclear Physics
(Nuclear Theory) under contract number DE-SC0013470. P.V.G. was supported by
the U.S. Department of Energy Award No. DE-SC0019095. P.V.G. is grateful for
the support and hospitality of the Fermilab theory group. Fermilab is operated
by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the
United States Department of Energy.

##  Appendix A Potentials

For convenience, we provide the explicit expressions of the coefficients
$V^{(k)}_{s,\text{RS}^{\prime}}$ and $V^{(k)}_{o,\text{RS}^{\prime}}$ used in
the potentials (10) and (13). In the singlet case, we have

$V^{(0)}_{s,\text{RS}^{\prime}}=-\frac{C_{F}}{r}\,,$ (22)
$V^{(1)}_{s,\text{RS}^{\prime}}=-\frac{C_{F}}{4\pi r}\left[a_{1}+2\log(r\nu
e^{\gamma_{E}})\beta_{0}\right]+2\delta m^{(1)}_{\text{RS}^{\prime}}\,,$ (23)
$V^{(2)}_{s,\text{RS}^{\prime}}=-\frac{C_{F}}{(4\pi)^{2}r}\left[(a_{2}+\log(r\nu
e^{\gamma_{E}})(4a_{1}\beta_{0}+2\beta_{1}))+\left(\frac{\pi^{2}}{3}+4\log^{2}(r\nu
e^{\gamma_{E}})\right)\beta^{2}_{0}\right]+2\delta
m^{(2)}_{\text{RS}^{\prime}}\,,$ (24)
$\displaystyle\begin{split}V^{(3)}_{s,\text{RS}^{\prime}}=&-\frac{C_{F}}{(4\pi)^{3}r}\bigg{[}(a_{3}+a_{1}\beta^{2}_{0}\pi^{2}+\frac{5\pi^{2}}{6}\beta_{0}\beta_{1}+16\zeta^{3}(3)\beta^{3}_{0}\\\
&+(2\pi^{2}\beta^{3}_{0}+6a_{2}\beta_{0}+4a_{1}\beta_{1}+2\beta_{2}+144\pi^{2})\log(r\nu
e^{\gamma_{E}})\\\ &+(12a_{1}\beta^{2}_{0}+10\beta_{0}\beta_{1})\log^{2}(r\nu
e^{\gamma_{E}})+8\beta^{3}_{0}\log^{3}(r\nu
e^{\gamma_{E}})+144\pi^{2}\log\frac{\nu_{\textrm{us}}}{\nu}\bigg{]}+2\delta
m^{(3)}_{\text{RS}^{\prime}}\,,\end{split}$ (25)

where $\beta_{i}$ are the coefficients of the $\beta$-function. We work in the
convention where $\beta_{0}=11C_{A}/3-4T_{F}n_{f}/3$, with $n_{f}$ the number
of massless flavors, which we set to $n_{f}=3$ in our bottomonium calculations
Brambilla  _et al._ (2002). The scale $\nu$ is the renormalization scale and
$\nu_{\text{us}}$ the factorization scale separating the soft from the
ultrasoft region (see Sec. II.2 for discussion of the scales). The SU($N$)
Casimirs are $C_{F}=T_{F}(N_{c}^{2}-1)/N_{c}$, $C_{A}=2T_{F}N_{c}$ with
$N_{c}$ the number of colors and $T_{F}=1/2$. The expressions for
$a_{1},a_{2}$, and $a_{3}$ were computed in Refs. Fischler (1977); Billoire
(1980), Refs. Peter (1997); Schroder (1999) and Refs. Anzai  _et al._ (2010);
Smirnov  _et al._ (2010), respectively. The ultrasoft correction in
$V^{(3)}_{s,\text{RS}^{\prime}}$ was computed in Ref. Brambilla  _et al._
(1999). The renormalon contributions $\delta m_{\text{RS}^{\prime}}$ are given
by Pineda (2001)

$\delta m_{\text{RS}^{\prime}}=\nu_{f}\sum^{3}_{k=0}\delta
m^{(k)}_{\text{RS}^{\prime}}\left(\frac{\nu_{f}}{\nu}\right)\alpha^{k+1}_{s}(\nu)\,,$
(26)

with

$\displaystyle\delta m_{\mathrm{RS}^{\prime}}^{(0)}=0,\qquad\delta
m_{\mathrm{RS}^{\prime}}^{(1)}\left(\frac{\nu_{f}}{\nu}\right)=N_{m}\frac{\beta_{0}}{2\pi}S(1,b),$
(27) $\displaystyle\delta
m_{\mathrm{RS}^{\prime}}^{(2)}\left(\frac{\nu_{f}}{\nu}\right)=N_{m}\left(\frac{\beta_{0}}{2\pi}\right)\left[S(1,b)\frac{2d_{0}\left(\nu,\nu_{f}\right)}{\pi}+\left(\frac{\beta_{0}}{2\pi}\right)S(2,b)\right],$
(28) $\displaystyle\begin{split}&\delta
m_{\mathrm{RS}^{\prime}}^{(3)}\left(\frac{\nu_{f}}{\nu}\right)=N_{m}\left(\frac{\beta_{0}}{2\pi}\right)\bigg{[}S(1,b)\frac{3d_{0}^{2}\left(\nu,\nu_{f}\right)+2d_{1}\left(\nu,\nu_{f}\right)}{\pi^{2}}+\left(\frac{\beta_{0}}{2\pi}\right)S(2,b)\frac{3d_{0}\left(\nu,\nu_{f}\right)}{\pi}\\\
&\qquad\qquad\qquad\qquad\qquad\qquad+\left(\frac{\beta_{0}}{2\pi}\right)^{2}S(3,b)\bigg{]},\end{split}$
(29)

where

$S(n,b)=\sum_{k=0}^{2}c_{k}\frac{\Gamma(n+1+b-k)}{\Gamma(1+b-k)},\quad
d_{k}\left(\nu,\nu_{f}\right)=\frac{\beta_{k}}{2^{1+2k}}\ln\left(\frac{\nu}{\nu_{f}}\right),$
(30)

and

$c_{0}=1,\quad
c_{1}=\frac{\beta_{1}^{2}-\beta_{0}\beta_{2}}{4\beta_{0}^{4}b},\quad
c_{2}=\frac{\beta_{1}^{4}+4\beta_{0}^{3}\beta_{1}\beta_{2}-2\beta_{0}\beta_{1}^{2}\beta_{2}+\beta_{0}^{2}\left(-2\beta_{1}^{3}+\beta_{2}^{2}\right)-2\beta_{0}^{4}\beta_{3}}{32\beta_{0}^{8}b(b-1)}\,.$
(31)

The scale $\nu_{f}$ is the renormalon factorization scale and
$b=\beta_{1}/(2\beta^{2}_{0})$. For $N_{m}$ we use $N_{m}=0.563126$ from Ref.
Ayala  _et al._ (2014). There are also other determinations, e.g.,
$N_{m}=0.574974$ from Ref. Pineda (2001), and $N_{m}=0.535$ from Ref. Komijani
(2017), which only slightly alter the spectrum and the bottom mass.

For the octet potential, we use the three-loop one from Ref. Anzai  _et al._
(2013). The coefficients $V^{(k)}_{o,\text{RS}}$ read

$V^{(0)}_{o,\text{RS}^{\prime}}=-\frac{C^{[8]}}{r}\,,$ (32)
$V^{(1)}_{o,\text{RS}^{\prime}}=-\frac{C^{[8]}}{4\pi
r}\left[a^{[8]}_{1}+2\beta_{0}\log(r\nu e^{\gamma_{E}})\right]-\delta
V^{(1)}_{o,\text{RS}^{\prime}}\,,$ (33)
$V^{(2)}_{o,\text{RS}^{\prime}}=-\frac{C^{[8]}}{16\pi^{2}r}\left[a^{[8]}_{2}+\beta^{2}_{0}\left(\frac{\pi^{2}}{3}+4\log^{2}(r\nu
e^{\gamma_{E}})\right)+\log(r\nu
e^{\gamma_{E}})(4a^{[8]}_{1}\beta_{0}+2\beta_{1})\right]-\delta
V^{(2)}_{o,\text{RS}^{\prime}}\,,$ (34)
$\displaystyle\begin{split}V^{(3)}_{o,\text{RS}^{\prime}}=-&\frac{C^{[8]}}{64\pi^{2}r}\bigg{[}a^{[8]}_{3}+8\beta^{3}_{0}\big{(}\log^{3}(r\nu
e^{\gamma_{E}})+\frac{\pi^{2}}{4}\log(r\nu e^{\gamma_{E}})+2\zeta(3)\big{)}\\\
&+(10\beta_{0}\beta_{1}+12\beta^{2}_{0}a^{[8]}_{1})\bigg{(}\log^{2}(r\nu
e^{\gamma_{E}})+\frac{\pi^{2}}{12}\bigg{)}\\\
&+(2\beta_{2}+4\beta_{1}a^{[8]}_{1}+6\beta_{0}a^{[8]}_{2})\log(r\nu
e^{\gamma_{E}})+\frac{64\pi^{2}}{3}N^{3}_{c}\log(r\nu)\bigg{]}-\delta
V^{(3)}_{o,\text{RS}^{\prime}}\,,\end{split}$ (35)

with $C^{[8]}=C_{F}-C_{A}/2$. The renormalon contributions $\delta
V^{(i)}_{o,\text{RS}^{\prime}}$ have the same form as $\delta
m^{(i)}_{\text{RS}^{\prime}}$ but with normalization $N_{V_{o}}$ instead of
$N_{m}$, which is set to $N_{V_{o}}=$0.114$$ Bali and Pineda (2004). The
coefficients $a^{[8]}_{1}$ and $a^{[8]}_{2}$ are defined as

$a^{[8]}_{1}=a_{1}\,,$ (36)
$a^{[8]}_{2}=a_{2}+N^{2}_{c}\pi^{2}(\pi^{2}-12)\,,$ (37)

where $a_{1}$ and $a_{2}$ are from the singlet potential. The coefficient
$a^{[8]}_{3}$ has the form

$a^{[8]}_{3}=a^{[8],(3)}_{3}n^{3}_{f}+a^{[8],(2)}_{3}n^{2}_{f}+a^{[8],(1)}_{3}n_{f}+a^{[8],(0)}_{3}\,,$
(38)

with

$a^{[8],(3)}_{3}=-\left(\frac{20}{9}\right)^{3}T^{3}_{F}\,,$ (39)
$a^{[8],(2)}_{3}=\bigg{(}\frac{12541}{243}+\frac{368\zeta(3)}{3}+\frac{64\pi^{4}}{135}\bigg{)}C_{A}T^{2}_{F}+\bigg{(}\frac{14002}{81}-\frac{416\zeta(3)}{3}\bigg{)}C_{F}T^{2}_{F}\,,$
(40)
$a^{[8],(1)}_{3}=-327.2N^{2}_{c}+\frac{66133}{648}-\frac{112\pi^{2}}{9}-\frac{272\zeta(3)}{3}+\frac{8\pi^{4}}{3}-\frac{32\pi^{2}\zeta(3)}{3}+20\zeta(5)-12.6\frac{1}{N^{2}_{c}}\,,$
(41) $a^{[8],(0)}_{3}=-17.0N_{c}+499.4N^{3}_{c}-97.6N^{3}_{c}\,.$ (42)

## References

  * Matsui and Satz (1986) T. Matsui and H. Satz, Phys. Lett. B 178, 416 (1986).
  * Andronic  _et al._ (2024) A. Andronic _et al._ , (2024), arXiv:2402.04366 [nucl-th] .
  * Brambilla  _et al._ (2005) N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Rev. Mod. Phys. 77, 1423 (2005), arXiv:hep-ph/0410047 .
  * Brambilla  _et al._ (2017) N. Brambilla, M. A. Escobedo, J. Soto, and A. Vairo, Phys. Rev. D 96, 034021 (2017), arXiv:1612.07248 .
  * Brambilla  _et al._ (2018) N. Brambilla, M. A. Escobedo, J. Soto, and A. Vairo, Phys. Rev. D 97, 074009 (2018), arXiv:1711.04515 .
  * Brambilla  _et al._ (2021a) N. Brambilla, M. A. Escobedo, M. Strickland, A. Vairo, P. Vander Griend, and J. H. Weber, JHEP 05, 136 (2021a), arXiv:2012.01240 .
  * Akamatsu (2022) Y. Akamatsu, Prog. Part. Nucl. Phys. 123, 103932 (2022), arXiv:2009.10559 [nucl-th] .
  * Yao (2021) X. Yao, Int. J. Mod. Phys. A 36, 2130010 (2021), arXiv:2102.01736 [hep-ph] .
  * Miura  _et al._ (2022) T. Miura, Y. Akamatsu, M. Asakawa, and Y. Kaida, Phys. Rev. D 106, 074001 (2022), arXiv:2205.15551 [nucl-th] .
  * Brambilla  _et al._ (2022) N. Brambilla, M. A. Escobedo, A. Islam, M. Strickland, A. Tiwari, A. Vairo, and P. Vander Griend, JHEP 08, 303 (2022), arXiv:2205.10289 .
  * Brambilla  _et al._ (2023) N. Brambilla, M. A. Escobedo, A. Islam, M. Strickland, A. Tiwari, A. Vairo, and P. Vander Griend, Phys. Rev. D 108, L011502 (2023), arXiv:2302.11826 .
  * Acharya  _et al._ (2021) S. Acharya _et al._ (ALICE), Phys. Lett. B 822, 136579 (2021), arXiv:2011.05758 .
  * Aad  _et al._ (2023) G. Aad _et al._ (ATLAS), Phys. Rev. C 107, 054912 (2023), arXiv:2205.03042 .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), Phys. Lett. B 790, 270 (2019), arXiv:1805.09215 .
  * Tumasyan  _et al._ (2023) A. Tumasyan _et al._ (CMS), (2023), arXiv:2303.17026 .
  * Caswell and Lepage (1986) W. E. Caswell and G. P. Lepage, Phys. Lett. B 167, 437 (1986).
  * Bodwin  _et al._ (1995) G. T. Bodwin, E. Braaten, and G. P. Lepage, Phys. Rev. D 51, 1125 (1995), [Erratum: Phys.Rev.D 55, 5853 (1997)], arXiv:hep-ph/9407339 .
  * Pineda and Soto (1998) A. Pineda and J. Soto, Nucl. Phys. B Proc. Suppl. 64, 428 (1998), arXiv:hep-ph/9707481 .
  * Brambilla  _et al._ (2000) N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Nucl. Phys. B 566, 275 (2000), arXiv:hep-ph/9907240 .
  * Brambilla and Vairo (1999) N. Brambilla and A. Vairo, in  _13th Annual HUGS AT CEBAF_ (1999) pp. 151–220, arXiv:hep-ph/9904330 .
  * Brambilla  _et al._ (1999) N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Phys. Rev. D 60, 091502 (1999), arXiv:hep-ph/9903355 .
  * Anzai  _et al._ (2010) C. Anzai, Y. Kiyo, and Y. Sumino, Phys. Rev. Lett. 104, 112003 (2010), arXiv:0911.4335 [hep-ph] .
  * Smirnov  _et al._ (2010) A. V. Smirnov, V. A. Smirnov, and M. Steinhauser, Phys. Rev. Lett. 104, 112002 (2010), arXiv:0911.4742 [hep-ph] .
  * Pineda (2001) A. Pineda, JHEP 06, 022 (2001), arXiv:hep-ph/0105008 .
  * Brambilla  _et al._ (2001) N. Brambilla, Y. Sumino, and A. Vairo, Phys. Lett. B 513, 381 (2001), arXiv:hep-ph/0101305 .
  * Sumino (2016) Y. Sumino, PoS LL2016, 011 (2016), arXiv:1607.03469 [hep-ph] .
  * Peset  _et al._ (2018) C. Peset, A. Pineda, and J. Segovia, Phys. Rev. D 98, 094003 (2018), arXiv:1809.09124 [hep-ph] .
  * Kiyo  _et al._ (2010) Y. Kiyo, A. Pineda, and A. Signer, Nucl. Phys. B 841, 231 (2010), arXiv:1006.2685 [hep-ph] .
  * Pineda and Segovia (2013) A. Pineda and J. Segovia, Phys. Rev. D 87, 074024 (2013), arXiv:1302.3528 [hep-ph] .
  * Segovia  _et al._ (2019) J. Segovia, S. Steinbeißer, and A. Vairo, Phys. Rev. D 99, 074011 (2019), arXiv:1811.07590 [hep-ph] .
  * Titard and Yndurain (1994) S. Titard and F. J. Yndurain, Phys. Rev. D 49, 6007 (1994), arXiv:hep-ph/9310236 .
  * Brambilla and Vairo (2000) N. Brambilla and A. Vairo, in  _5th Workshop on QCD (QCD 2000)_ (2000) pp. 8–16, arXiv:hep-ph/0004192 .
  * Larsen  _et al._ (2019) R. Larsen, S. Meinel, S. Mukherjee, and P. Petreczky, Phys. Rev. D 100, 074506 (2019), arXiv:1908.08437 [hep-lat] .
  * Omar  _et al._ (2022) H. B. Omar, M. A. Escobedo, A. Islam, M. Strickland, S. Thapa, P. Vander Griend, and J. H. Weber, Comput. Phys. Commun. 273, 108266 (2022), arXiv:2107.06147 .
  * Lindblad (1976) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * Mølmer  _et al._ (1993) K. Mølmer, Y. Castin, and J. Dalibard, J. Opt. Soc. Am. B 10, 524 (1993).
  * Herren and Steinhauser (2018) F. Herren and M. Steinhauser, Comput. Phys. Commun. 224, 333 (2018), arXiv:1703.03751 .
  * Anzai  _et al._ (2013) C. Anzai, M. Prausa, A. V. Smirnov, V. A. Smirnov, and M. Steinhauser, Phys. Rev. D 88, 054030 (2013), arXiv:1308.1202 [hep-ph] .
  * Pineda (1998) A. Pineda,  _Heavy Quarkonium and Nonrelativistic Effective Field Theories_ , Ph.D. thesis, Unversitat de Barcelona (1998).
  * Hoang  _et al._ (1999) A. H. Hoang, M. C. Smith, T. Stelzer, and S. Willenbrock, Phys. Rev. D 59, 114014 (1999), arXiv:hep-ph/9804227 .
  * Bali and Pineda (2004) G. S. Bali and A. Pineda, Phys. Rev. D 69, 094001 (2004), arXiv:hep-ph/0310130 .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Brambilla  _et al._ (2019) N. Brambilla, M. A. Escobedo, A. Vairo, and P. Vander Griend, Phys. Rev. D 100, 054025 (2019), arXiv:1903.08063 [hep-ph] .
  * Eller  _et al._ (2019) A. M. Eller, J. Ghiglieri, and G. D. Moore, Phys. Rev. D 99, 094042 (2019), [Erratum: Phys.Rev.D 102, 039901 (2020)], arXiv:1903.08064 [hep-ph] .
  * Brambilla  _et al._ (2020) N. Brambilla, V. Leino, P. Petreczky, and A. Vairo, Phys. Rev. D 102, 074503 (2020), arXiv:2007.10078 [hep-lat] .
  * Foster and Michael (1999) M. Foster and C. Michael (UKQCD), Phys. Rev. D 59, 094509 (1999), arXiv:hep-lat/9811010 .
  * Narison (2023) S. Narison, in  _26th High-Energy Physics International Conference in QCD_ (2023) arXiv:2309.05342 [hep-ph] .
  * Herr  _et al._ (2024) J. Herr, C. Schlosser, and M. Wagner, Phys. Rev. D 109, 034516 (2024), arXiv:2306.09902 [hep-lat] .
  * Alqahtani and Strickland (2021) M. Alqahtani and M. Strickland, Eur. Phys. J. C 81, 1022 (2021), arXiv:2008.07657 [nucl-th] .
  * Brambilla  _et al._ (2021b) N. Brambilla, M. A. Escobedo, M. Strickland, A. Vairo, P. Vander Griend, and J. H. Weber, Phys. Rev. D 104, 094049 (2021b), arXiv:2107.06222 .
  * Brambilla  _et al._ (2002) N. Brambilla, Y. Sumino, and A. Vairo, Phys. Rev. D 65, 034001 (2002), arXiv:hep-ph/0108084 .
  * Fischler (1977) W. Fischler, Nucl. Phys. B 129, 157 (1977).
  * Billoire (1980) A. Billoire, Phys. Lett. B 92, 343 (1980).
  * Peter (1997) M. Peter, Nucl. Phys. B 501, 471 (1997), arXiv:hep-ph/9702245 .
  * Schroder (1999) Y. Schroder, Phys. Lett. B 447, 321 (1999), arXiv:hep-ph/9812205 .
  * Ayala  _et al._ (2014) C. Ayala, G. Cvetič, and A. Pineda, JHEP 09, 045 (2014), arXiv:1407.2128 [hep-ph] .
  * Komijani (2017) J. Komijani, JHEP 08, 062 (2017), arXiv:1701.00347 [hep-ph] .