aainstitutetext: Physik-Department, Technische Universität München, James-
Franck-Str. 1, 85748 Garching, Germanybbinstitutetext: Institute for Advanced
Study, Technische Universität München, Lichtenbergstrasse 2 a, 85748 Garching,
Germanyccinstitutetext: Munich Data Science Institute, Technische Universität
München, Walther-von-Dyck-Strasse 10, 85748 Garching, Germanyddinstitutetext:
Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago
de Compostela. E-15782, Galicia, Spaineeinstitutetext: Department of Physics,
Kent State University, Kent, OH 44242, United States

# Heavy quarkonium dynamics at next-to-leading order in the binding energy
over temperature

Nora Brambilla d Miguel Ángel Escobedo e Ajaharul Islam e Michael Strickland e
Anurag Tiwari a Antonio Vairo a Peter Vander Griend<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Using the potential non-relativistic quantum chromodynamics (pNRQCD) effective
field theory, we derive a Lindblad equation for the evolution of the heavy-
quarkonium reduced density matrix that is accurate to next-to-leading order
(NLO) in the ratio of the binding energy of the state to the temperature of
the medium. The resulting NLO Lindblad equation can be used to more reliably
describe heavy-quarkonium evolution in the quark-gluon plasma at low
temperatures compared to the leading-order truncation. For phenomenological
application, we numerically solve the resulting NLO Lindblad equation using
the quantum trajectories algorithm. To achieve this, we map the solution of
the three-dimensional Lindblad equation to the solution of an ensemble of one-
dimensional Schrödinger evolutions with Monte-Carlo sampled quantum jumps.
Averaging over the Monte-Carlo sampled quantum jumps, we obtain the solution
to the NLO Lindblad equation without truncation in the angular momentum
quantum number of the states considered. We also consider the evolution of the
system using only the complex effective Hamiltonian without stochastic jumps
and find that this provides a reliable approximation for the ground state
survival probability at LO and NLO. Finally, we make comparisons with our
prior leading-order pNRQCD results and experimental data available from the
ATLAS, ALICE, and CMS collaborations.

###### Keywords:

heavy quarkonium suppression, NLO Lindblad equation, heavy-ion collision,
quantum trajectories method, open quantum systems

††preprint: TUM-EFT 169/22

##  1 Introduction

In the pioneering work of Matsui and Satz Matsui:1986dk , heavy quarkonium
suppression was proposed as a signal of the formation of a quark gluon plasma
(QGP). Since then, quarkonium suppression studies have been an important part
of the heavy-ion program of several experimental facilities. For recent
studies at LHC and RHIC experiments, see refs. Acharya:2020kls ; ATLAS5TeV ;
Sirunyan:2018nsz ; CMS:2017ycw ; CMS:2020efs ; ALICE:2019pox ; STAR:2013kwk ;
PHENIX:2014tbe ; STAR:2016pof ; CMSupsilonQM2022 . In ref. Matsui:1986dk ,
Matsui and Satz conjectured that the screening of chromoelectric fields
generated by the medium at distances proportional to the inverse of the Debye
mass induces quarkonium dissociation and consequently quarkonium suppression
in medium. In other words, the heavy quark-antiquark potential acquires a
screening factor $\exp({-m_{D}r})$, where $r$ is the quark-antiquark distance
and $m_{D}$ the Debye mass which is proportional to the temperature. Within
this screening picture, quarkonia states of different radii dissociate at
different temperatures providing a QGP thermometer.

The last decades have seen a change in this screening driven dissociation
paradigm. A seminal perturbative calculation of the heavy-quark potential by
Laine, et. al., found that in the screening regime $r\sim m_{D}^{-1}$, in
addition to Debye screening of the real part of the potential, a large
imaginary contribution to the potential also arises Laine:2006ns .
Perturbative calculations based on nonrelativistic effective theories of QCD
at finite temperature Brambilla:2008cx ; Beraudo:2007ky ; Escobedo:2008sy ;
Brambilla:2010vq ; Brambilla:2011sg ; Brambilla:2013dpa  and perturbative QCD
resummation Dumitru:2009fy  confirmed this result. Different temperature
regimes have been treated using pNRQCD; outside the screening regime,
different patterns of finite temperature corrections to the potential emerge,
with all possibilities possessing imaginary contributions.

At temperatures relevant in current heavy ion collision experiments, the
imaginary part of the potential dominates the screening contributions, and it
is thus the imaginary contributions to the potential, rather than screening
effects, which trigger in-medium quarkonium dissolution and suppression. The
real part of the potential may be approximated by the quark-antiquark gauge
invariant free energy in some cases Berwein:2017thy ; Bazavov:2018wmo . The
imaginary part of the potential has its origin in two effects: Landau damping
(related to parton dissociation Brambilla:2013dpa ) and singlet to octet
transitions (an effect particular to QCD Brambilla:2008cx  related to the
gluodissociation process Brambilla:2011sg ). Nonperturbative lattice QCD and
classical real-time lattice measurements of the imaginary part of the
potential were recently performed in refs. Rothkopf:2011db ; Rothkopf:2019ipj
; Petreczky:2018xuh ; Bala:2019cqu  and Laine:2007qy ; Lehmann:2020fjt ;
Boguslavski:2020bxt , respectively. The imaginary part of the potential
indicates that quarkonium has a medium-induced decay width and, therefore, a
finite lifetime. There are situations in which the imaginary part of the
potential is as large as the real part and cannot be treated as a
perturbation, and, as mentioned, in the small coupling limit, this happens at
a temperature smaller than the screening temperature. Extensive
phenomenological studies on the impact of the imaginary part of the static
potential were performed in refs. Strickland:2011mw ; Strickland:2011aa ;
Krouppa:2015yoa ; Islam:2020bnp ; Islam:2020gdv .

In addition to the thermal width and screening effects, an accurate
description of in-medium quarkonium evolution must account for the
recombination process Braun-Munzinger:2000csl ; Thews:2000rj  in which an
unbound heavy quark-antiquark pair recombines inside the medium to form a new
bound state. One may distinguish between correlated recombination and
uncorrelated recombination. The former describes the recombination of heavy
quarks which initially belonged to a bound state which dissociated; the latter
increases with the number of heavy quarks created in the collision and is,
therefore, believed to have a significant role in charmonium production at the
LHC.

It is challenging to find a theoretical description that can consistently
describe all these effects. On the one hand, when considering screening, we
aim to understand if bound state formation is possible. This requires a
quantum description. On the other hand, decays and recombination require
consideration of transitions in which the heavy quark pairs change their color
state due to the interaction with the medium. Therefore, we need a formalism
that is quantum and describes exchanges of momentum, color, energy, etc., with
the medium. The Open Quantum Systems (OQS) formalism Breuer:2002pc  offers a
natural possibility to achieve such a description. Within this framework, we
regard the heavy quark-antiquark state as an open quantum system interacting
with an environment (the medium). The key object here is the reduced density
matrix, obtained from the full density matrix by performing a trace over the
degrees of freedom of the environment. The resulting equation describing the
evolution of the reduced density matrix is called the master equation. In
recent years, many works have addressed quarkonium suppression using this
formalism Akamatsu:2011se ; Akamatsu:2020ypb ; Akamatsu:2014qsa ;
Miura:2019ssi ; Brambilla:2016wgg ; Brambilla:2017zei ; Brambilla:2019tpt ;
Sharma:2019xum ; Blaizot:2015hya ; Blaizot:2017ypk ; Blaizot:2018oev ;
Blaizot:2021xqa ; Yao:2018nmy ; Yao:2020xzw ; Yao:2020eqy ; Katz:2015qja . In
this work, we focus on an approach combining the OQS framework with the use of
effective field theories (EFTs) to exploit the non-relativistic nature of
heavy quarkonium.

Due to its nonrelativistic nature, i.e., $v\ll 1$, where $v$ is the velocity
of the heavy quarks around the center of mass, heavy quarkonium possesses a
number of widely separated energy scales. At zero temperature, we identify:
the heavy quark mass $M$, the inverse of the typical radius $1/r\sim Mv$, and
the binding energy $E\sim Mv^{2}$. An additional relevant scale is
$\Lambda_{\rm QCD}$ at which nonperturbative effects become dominant.
Nonperturbative effects and the breaking of the weak coupling expansion when
$\alpha_{\rm s}/v\sim 1$ challenge the use of perturbation theory to describe
quarkonium physics, whereas at the same time lattice QCD computations struggle
to simultaneously accomodate widely separated scales on sufficiently large and
fine lattices. Effective field theories address and manage these problems.

An EFT is a quantum field theory that gives the same results as the more
general theory at any order of the expansion in a small parameter but is
restricted to a reduced kinematical region. The Lagrangian of an EFT has an
infinite number of terms. However, given a desired precision, only a finite
number of them are needed. Nonrelativistic QCD (NRQCD) Caswell:1985ui ;
Bodwin:1994jh  is an EFT that is equivalent to QCD for energies much smaller
than $M$. It is useful to go a step further and use potential NRQCD (pNRQCD)
Pineda:1997bj ; Brambilla:1999xf ; Brambilla:2004jw , an EFT equivalent to QCD
for energies much smaller than $Mv$. Such a theory makes it explicit that the
appropriate zeroth order problem for the quarkonium bound state dynamics
amounts to solving a Schrödinger equation and gives a field theoretical
definition of the potentials as Wilson coefficients. In pNRQCD, the degrees of
freedom are not heavy quarks and antiquarks but the color singlet and octet
states that one can form with pairs of them. At leading order (LO) in the
power counting, the evolution of the singlet field is given by a Schrödinger
equation. However, at higher orders, interactions with gauge fields with
energy smaller than $Mv$ become important and contain retardation effects.
Systematic calculations of energy levels, decays, and transitions are made
possible by the EFT and the power counting allows one to estimate the
uncertainty in the prediction of any observable.

Non-relativistic EFTs can also be used to study quarkonium at finite
temperature Brambilla:2008cx ; Escobedo:2008sy ; Brambilla:2010vq ;
Brambilla:2011sg ; Brambilla:2013dpa . In this case, we must consider also the
energy scales induced by the medium, e.g., the temperature $T$. In particular,
pNRQCD can be realized in medium and combined with the OQS framework to derive
a master equation governing the evolution of in-medium heavy quarkonium.
References Brambilla:2016wgg ; Brambilla:2017zei ; Brambilla:2019tpt  analyzed
the regime $Mv\gg T$. In this case, the matching between NRQCD and pNRQCD is
unaffected by the medium, and one can perform the calculation of the relevant
diagrams contributing to the in-medium evolution using the $T=0$ pNRQCD
Lagrangian. Using the closed path real-time formalism, the authors of refs.
Brambilla:2016wgg ; Brambilla:2017zei  derived a master equation for the in-
medium heavy quarkonium color singlet and color octet density matrices. We
emphasize that the theory, methods, and results are fully quantum, nonabelian,
and heavy quark number conserving. With the further assumption $T\gg E$ the
master equation could be cast in Lindblad form Gorini:1975nb ; Lindblad:1975ef
at leading order in the $E/T$ expansion. In a series of works
Brambilla:2016wgg ; Brambilla:2017zei ; Brambilla:2020qwo ; Brambilla:2021wkt
, a number of phenomenological predictions were made for bottomonium
suppression in the temperature regime $Mv\gg T\gg E$; physically, this amounts
to using bottomonium states of small radius as probes of a strongly coupled
QGP. In the most recent works Brambilla:2020qwo ; Brambilla:2021wkt , the
QTraj code of ref. Omar:2021kra , which implements the Monte Carlo quantum
trajectories algorithm Dalibard:1992zz , was used to solve the Lindblad
equation at substantially reduced computational cost compared to previous
works, while simultaneously implementing a realistic medium evolution by
coupling to a state of the art 3D viscous hydrodynamics code. In this way, one
obtains a good description of existing experimental results on bottomonium
suppression and in-medium observables including the nuclear modification
factor $R_{AA}$ and the elliptic flow $v_{2}$. We note that these results do
not require, in principle, the fixing of free parameters to experimental data:
the strong coupling $\alpha_{\rm s}$ we take from the Particle Data Group
Zyla:2020zbs , for the bottom mass, we work in the $1S$ scheme where
$m_{b}=m_{\Upsilon(1S)}/2$ with the $\Upsilon(1S)$ mass also from the Particle
Data Group, while the nonperturbative transport coefficients $\kappa$ and
$\gamma$ describing the strongly coupled QGP may be taken from direct or
indirect lattice measurements, as they are expressed in terms of gauge
invariant correlators of chromoelectric fields.

In the previous phenomenological studies Brambilla:2020qwo ; Brambilla:2021wkt
, the requirement that the inequality $T\gg E$ be fulfilled over the course of
the evolution necessitated the termination of the coupling to the medium at
$T_{f}=250$ MeV. In the present paper, we aim to improve on this by
considering the next-to-leading order (NLO) corrections in the $E/T$
expansion. Previous computations in other models Blaizot:2018oev ;
Miura:2019ssi ; DeBoni:2017ocl  indicate that some of the NLO corrections are
related to the drag force on a particle moving in the medium. In this way,
heavy quarkonium loses energy to the medium thus bringing it closer to a
thermal distribution. We note that in the Abelian case considered in refs.
Blaizot:2018oev ; Miura:2019ssi ; DeBoni:2017ocl  the NLO corrections in the
$E/T$ expansion bring the system to thermal equilibrium.

In this paper, we compute the first corrections in the $E/T$ expansion to the
master equation of pNRQCD in the region $Mv\gg T\sim gT\gg E$ and explore
their phenomenological impact. As explained, we work in QCD, in a fully
quantum and nonabelian setup, considering the case of quarkonium evolution in
a strongly coupled medium. Moreover, we implement a realistic hydrodynamic
evolution. Phenomenological results at NLO in such a case have never been
obtained before. As such, our results may be relevant not only to the study of
the nonequilibrium evolution of quarkonium in medium but also to the study of
the approach to equilibrium. We plan to address the second point in a future
publication and here focus on the first point, the phenomenological impact of
NLO corrections.

This paper is organized as follows. In sec. 2, we extend our open quantum
system (OQS) treatment of in-medium heavy quarkonium to include corrections at
NLO in the $E/T$ expansion and construct an NLO evolution equation of Lindblad
form for the heavy quarkonium reduced density matrix. Taking advantage of the
isotropy of the problem, we expand the Lindblad equation in spherical
harmonics to rewrite the 3-dimensional evolution in terms of a 1-dimensional
evolution; we discuss the form of the collapse operators and in-medium widths
in this expansion. In sec. 3, we discuss the implementation of the quantum
trajectories method at this order. In sec. 4, we present our results. On the
theoretical side, we discuss the effect of the NLO terms on the in-medium
widths and survival probabilities and the effect of quantum jumps. On the
phenomenological side, we present our results for the nuclear modification
factor $R_{AA}$ of the $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$ as a
function of number of participating nucleons $N_{\text{part}}$ and transverse
momentum $p_{T}$. We extend the evolution of the Linblad equation down to
$T_{f}=190$ MeV and still obtain a satisfactory description of the
experimental data. We comment on the dependence and features of our results in
relation to the present uncertainty of the lattice determinations of the
transport coefficients $\kappa$ and $\gamma$. We conclude and provide an
outlook on future prospects in sec. 5.

##  2 Derivation of the NLO Lindblad equation

In this section, we expand the evolution equation obeyed by the heavy
quarkonium reduced density matrix to include terms of order $E/T$ and derive a
Lindblad equation accurate at this order. In this way, we go beyond the
analysis at leading order in $E/T$ presented in refs. Brambilla:2016wgg ;
Brambilla:2017zei ; Brambilla:2020qwo ; Brambilla:2021wkt .

###  2.1 The pNRQCD master equation

The evolution of the reduced density matrix of heavy quarkonium in the regime
$Mv\gg T$ was derived in refs. Brambilla:2016wgg ; Brambilla:2017zei . If one
does not make any assumption on the relation between $T$ and $E$, then the
evolution equations that one obtains are

$\displaystyle\frac{d\rho_{s}(t)}{dt}$
$\displaystyle=-i[h_{s},\rho_{s}(t)]-\Sigma_{s}\rho_{s}(t)-\rho_{s}(t)\Sigma^{\dagger}_{s}+\Xi_{so}(\rho_{o}(t)),$
(1) $\displaystyle\frac{d\rho_{o}(t)}{dt}$
$\displaystyle=-i[h_{o},\rho_{o}(t)]-\Sigma_{o}\rho_{o}(t)-\rho_{o}(t)\Sigma_{o}^{\dagger}+\Xi_{os}(\rho_{s}(t))+\Xi_{oo}(\rho_{o}(t)),$
(2)

where, assuming the in medium correlators to be isotropic and time-translation
invariant and observing the evolution for a time $t$ larger than any other
length scale in the problem,111 These assumptions account for the difference
between eq. (8) and eq. (D.2) of ref. Brambilla:2017zei .

$\displaystyle\Sigma_{s}$ $\displaystyle=r_{i}A_{i}^{so\dagger},$ (3)
$\displaystyle\Sigma_{o}$
$\displaystyle=\frac{1}{N_{c}^{2}-1}r_{i}A_{i}^{os\dagger}+\frac{N_{c}^{2}-4}{2(N_{c}^{2}-1)}r_{i}A_{i}^{oo\dagger},$
(4) $\displaystyle\Xi_{so}(\rho_{o}(t))$
$\displaystyle=\frac{1}{N_{c}^{2}-1}\left(A_{i}^{os\dagger}\rho_{o}(t)r_{i}+r_{i}\rho_{o}(t)A_{i}^{os}\right),$
(5) $\displaystyle\Xi_{os}(\rho_{s}(t))$
$\displaystyle=A_{i}^{so\dagger}\rho_{s}(t)r_{i}+r_{i}\rho_{s}(t)A_{i}^{so},$
(6) $\displaystyle\Xi_{oo}(\rho_{o}(t))$
$\displaystyle=\frac{N_{c}^{2}-4}{2(N_{c}^{2}-1)}\left(A_{i}^{oo\dagger}\rho_{o}(t)r_{i}+r_{i}\rho_{o}(t)A_{i}^{oo}\right),$
(7)

with

$\displaystyle A_{i}^{uv}$
$\displaystyle=\frac{g^{2}}{6N_{c}}\int^{\infty}_{0}\text{d}s\,e^{-ih_{u}s}r_{i}e^{ih_{v}s}\langle\tilde{E}^{a}_{j}(0,\vec{0})\tilde{E}^{a}_{j}(s,\vec{0})\rangle,$
(8) $\displaystyle\tilde{E}^{a}_{i}(s,\vec{0})$
$\displaystyle=\Omega(s)^{\dagger}E^{a}_{i}(s,\vec{0})\Omega(s),$ (9)
$\displaystyle\Omega(s)$
$\displaystyle=\text{exp}\left[-ig\int_{-\infty}^{s}\text{d}s^{\prime}A_{0}(s^{\prime},\vec{0})\right].$
(10)

$h_{u,v}$ is the singlet or octet Hamiltonian: $h_{s,o}=\vec{p}^{2}/M+V_{s,o}$
where $V_{s}=-C_{f}\alpha_{\rm s}(1/a_{0})/r$ and $V_{o}=\alpha_{\rm
s}(1/a_{0})/2N_{c}r$, with $C_{f}=(N_{c}^{2}-1)/2N_{c}=4/3$ the Casimir of the
fundamental representation, $N_{c}=3$ the number of colors, and $\alpha_{\rm
s}(1/a_{0})$ the strong coupling evaluated at the inverse of the Bohr radius
$a_{0}$. Details about the derivation of these results can be found in refs.
Brambilla:2017zei ; Brambilla:2019tpt .

These equations can be written as the master equation

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

where

$\rho(t)=\begin{pmatrix}\rho_{s}(t)&0\\\
0&\rho_{o}(t)\end{pmatrix},\text{\quad}H=\begin{pmatrix}h_{s}+\text{Im}(\Sigma_{s})&0\\\
0&h_{o}+\text{Im}(\Sigma_{o})\end{pmatrix},$ (12)
$L_{i}^{0}=\begin{pmatrix}0&0\\\ 0&1\end{pmatrix}r_{i}\text{,
\quad}L_{i}^{1}=\begin{pmatrix}0&0\\\
0&\frac{N_{c}^{2}-4}{2(N_{c}^{2}-1)}A_{i}^{oo\dagger}\end{pmatrix},$ (13)
$L_{i}^{2}=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}r_{i}\text{,
\quad}L_{i}^{3}=\begin{pmatrix}0&\frac{1}{N_{c}^{2}-1}A_{i}^{os\dagger}\\\
A_{i}^{so\dagger}&0\end{pmatrix},$ (14)

and

$h=\begin{pmatrix}0&1&0&0\\\ 1&0&0&0\\\ 0&0&0&1\\\ 0&0&1&0\end{pmatrix}.$ (15)

If the matrix $h$ were a completely positive matrix, its eigenvalues would be
positive and we could make a change of variables to bring eq. (11) into
Lindblad form. However, as the matrix $h$ is not completely positive definite,
this manipulation is not possible. Any Markovian, trace preserving evolution
that preserves complete positivity of the density matrix can be written as a
Lindblad equation. As eq. (11) is Markovian and trace preserving but does not
preserve complete positivity of the density matrix, it follows that eq. (11)
allows for negative probabilities, a feature similar to the widely used
Caldeira–Leggett master equation Caldeira:1982iu . We note that complete
positivity is necessary for use of the quantum trajectories method.

###  2.2 Lindblad equation to leading-order in $E/T$

Although, as explained above, the evolution specified in eq. (11) does not, in
general, take the form of a Lindblad equation, in specific temperature regimes
the operators characterizing the medium interaction simplify, and the
resulting evolution equation can be rewrtitten in Lindblad form. We consider
the case $T\gg E$. In eq. (8), the thermal chromoelectric correlator decays
rapidly at times $s>1/T$, while the quarkonium Hamiltonian $h_{s,o}$ scales as
the binding energy $E$. The arguments of the exponentials in $A_{i}^{uv}$ thus
scale like $E/T$. In the limit $T\gg E$, the exponentials can be set to $1$,
and $A_{i}^{uv}$ can be written in terms of the transport coefficients
$\kappa$ and $\gamma$

$A_{i}^{uv}=\frac{r_{i}}{2}\left(\kappa-i\gamma\right)+\cdots,$ (16)

where the ellipsis indicates terms of order $E/T$ and higher; $\kappa$ is the
heavy quarkonium momentum diffusion coefficient, and $\gamma$ is its
dispersive counterpart

$\displaystyle\kappa$
$\displaystyle=\frac{g^{2}}{6N_{c}}\int_{0}^{\infty}\text{d}s\,\left\langle\left\\{\tilde{E}^{a}_{i}(s,\vec{0}),\tilde{E}^{a}_{i}(0,\vec{0})\right\\}\right\rangle,$
(17) $\displaystyle\gamma$
$\displaystyle=-i\frac{g^{2}}{6N_{c}}\int_{0}^{\infty}\text{d}s\,\left\langle\left[\tilde{E}^{a}_{i}(s,\vec{0}),\tilde{E}^{a}_{i}(0,\vec{0})\right]\right\rangle.$
(18)

At leading order in $E/T$, the $L^{n}$ are not linearly independent; as a
consequence, there exists a rotation of the vector $(L^{0},L^{1},L^{2},L^{3})$
that makes it orthogonal to the eigenspaces of the negative eigenvalues of the
matrix $h_{nm}$ (see appendix B for details). This allows to write the master
equation in Lindblad form

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

with Hamiltonian

$H=\begin{pmatrix}h_{s}+\text{Im}(\Sigma_{s})&0\\\
0&h_{o}+\text{Im}(\Sigma_{o})\end{pmatrix},$ (20)

where

$\text{Im}\left(\Sigma_{s}\right)=\frac{r^{2}}{2}\gamma,\quad\text{Im}\left(\Sigma_{o}\right)=\frac{N_{c}^{2}-2}{2(N_{c}^{2}-1)}\frac{r^{2}}{2}\gamma,$
(21)

and collapse operators

$\displaystyle C_{i}^{0}=$
$\displaystyle\sqrt{\frac{\kappa}{N_{c}^{2}-1}}r_{i}\begin{pmatrix}0&1\\\
\sqrt{N_{c}^{2}-1}&0\end{pmatrix},$ (22) $\displaystyle C_{i}^{1}=$
$\displaystyle\sqrt{\frac{\kappa(N_{c}^{2}-4)}{2(N_{c}^{2}-1)}}r_{i}\begin{pmatrix}0&0\\\
0&1\end{pmatrix}.$ (23)

This equation was first derived and solved in refs. Brambilla:2016wgg ;
Brambilla:2017zei  and further studied in refs. Brambilla:2020qwo ;
Brambilla:2021wkt  using the QTraj code of ref. Omar:2021kra .

###  2.3 Lindblad equation to next-to-leading order in $E/T$

In the previous subsections, we presented the master equation governing the
in-medium evolution of Coulombic quarkonium and showed that in the $T\gg E$
limit it takes the form of a Lindblad equation. In this subsection, we
consider contributions to $A_{i}^{uv}$ linear in $E/T$ and investigate the
conditions under which a Lindblad equation can still be derived. Making use of
the relation

$i\frac{g^{2}}{6N_{c}}\int_{0}^{\infty}\text{d}t\,t\,\Big{\langle}\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}=\frac{\kappa}{4T},$
(24)

which follows from the fluctuation-dissipation theorem (see appendix A for a
detailed derivation), we obtain corrections to eqs. (21)-(23) of order $E/T$
by expanding the exponentials of eq. (8) and retaining terms up to order
$h_{u,v}$. This results in $A_{i}^{uv}$ taking the form

$\displaystyle A_{i}^{uv}=$
$\displaystyle\frac{r_{i}}{2}(\kappa-i\gamma)+\kappa\left(-\frac{ip_{i}}{2MT}+\frac{\Delta
V_{uv}}{4T}r_{i}\right)+\cdots,$ (25)

where $\Delta V_{uv}=V_{u}-V_{v}$ is the difference of the $u$ and $v$
(singlet or octet) potentials and the ellipsis indicates terms of order
$(E/T)^{2}$ and higher. With these contributions to $A_{i}^{uv}$, the
operators $L_{i}^{1}$ and $L_{i}^{3}$ in the master equation take the form

$\displaystyle
L_{i}^{1}=\frac{N_{c}^{2}-4}{2(N_{c}^{2}-1)}\begin{pmatrix}0&0\\\
0&1\end{pmatrix}\left[\frac{r_{i}}{2}(\kappa+i\gamma)+\kappa\frac{ip_{i}}{2MT}\right],$
(26) $\displaystyle\begin{aligned}
L_{i}^{3}=&\begin{pmatrix}0&\frac{1}{N_{c}^{2}-1}\\\
0&0\end{pmatrix}\left[\frac{r_{i}}{2}(\kappa+i\gamma)+\kappa\left(\frac{ip_{i}}{2MT}+\frac{\Delta
V_{os}}{4T}r_{i}\right)\right]\\\ &+\begin{pmatrix}0&0\\\
1&0\end{pmatrix}\left[\frac{r_{i}}{2}(\kappa+i\gamma)+\kappa\left(\frac{ip_{i}}{2MT}+\frac{\Delta
V_{so}}{4T}r_{i}\right)\right],\end{aligned}$ (27)

while $L_{i}^{0}$ and $L_{i}^{2}$ keep the form given in eqs. (13) and (14),
respectively.

In contrast to the case at order $(E/T)^{0}$ considered in subsection 2.2, the
above master equation cannot be written in Lindblad form as the $L^{n}_{i}$
are linearly independent. However, the vector $(L^{0},L^{1},L^{2},L^{3})$ can
be rotated in such a way that only components of order $E/T$ project on the
eigenspaces of the negative eigenvalues of the matrix $h_{nm}$. Discarding
these components amounts to neglecting terms in the evolution equation of
order $(E/T)^{2}$. Hence, we can still write a Lindblad equation that is
accurate at order $E/T$

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

with Hamiltonian

$H=\begin{pmatrix}h_{s}+\text{Im}(\Sigma_{s})&0\\\
0&h_{o}+\text{Im}(\Sigma_{o})\end{pmatrix},$ (29)

where

$\text{Im}\left(\Sigma_{s}\right)=\frac{r^{2}}{2}\gamma+\frac{\kappa}{4MT}\\{r_{i},p_{i}\\}\text{,
\quad}\text{Im}\left(\Sigma_{o}\right)=\frac{N_{c}^{2}-2}{2(N_{c}^{2}-1)}\left(\frac{r^{2}}{2}\gamma+\frac{\kappa}{4MT}\\{r_{i},p_{i}\\}\right),$
(30)

and collapse operators

$\displaystyle\begin{aligned}
C_{i}^{0}=&\sqrt{\frac{\kappa}{N_{c}^{2}-1}}\begin{pmatrix}0&1\\\
0&0\end{pmatrix}\left(r_{i}+\frac{ip_{i}}{2MT}+\frac{\Delta
V_{os}}{4T}r_{i}\right)\\\ &+\sqrt{\kappa}\begin{pmatrix}0&0\\\
1&0\end{pmatrix}\left(r_{i}+\frac{ip_{i}}{2MT}+\frac{\Delta
V_{so}}{4T}r_{i}\right),\end{aligned}$ (31) $\displaystyle
C_{i}^{1}=\sqrt{\frac{\kappa(N_{c}^{2}-4)}{2(N_{c}^{2}-1)}}\begin{pmatrix}0&0\\\
0&1\end{pmatrix}\left(r_{i}+\frac{ip_{i}}{2MT}\right).$ (32)

For details, see appendix B.

At the price of introducing some spurious terms at order $(E/T)^{2}$, we have
derived a Lindblad equation accurate up to order $E/T$ and can thus make use
of the quantum trajectories algorithm to solve it. Note that similar
strategies have been used to improve the Caldeira–Leggett model diosi1993high
; diosi1993calderia ; gao1997dissipative ; vacchini2000completely . We can
understand the anti-commutator term $\\{r_{i},p_{i}\\}$ multiplying $\kappa$
in $\text{Im}\left(\Sigma_{s,o}\right)$ as a drag force. This identification
is clear if one uses the procedure described in ref. Blaizot:2017ypk  to
obtain a Langevin equation from eq. (28).

###  2.4 Lindblad equation to NLO in $E/T$ in the spherical basis

The six collapse operators in Eqs. (31) and (32) (along with the Hamiltonian
$H$) encode the full 3-dimensional evolution of Coulombic quarkonium of
binding energy $E$ propagating in a thermal medium of temperature $T$ in the
regime $T\gg E$.222 In practice, adding $E/T$ corrections serve to relax the
parametric condition $T\gg E$ to a milder $T\gtrsim E$ in the phenomenological
applications of sec. 4.  The collapse operators implement transitions between
quarkonium states of different color and angular momentum. This can be made
manifest in the angular momentum sector by projecting the density matrix onto
the spherical harmonics $Y_{lm}$

$\rho^{lm;l^{\prime}m^{\prime}}=\int
d\Omega\,d\Omega^{\prime}\,Y_{lm}\,\rho\,Y_{l^{\prime}m^{\prime}}.$ (33)

As the plasma is isotropic, i.e., there is no preferred direction in space,
only diagonal elements $l=l^{\prime}$ and $m=m^{\prime}$ are nonzero;
furthermore, all information can be encoded in

$\rho^{l}=\sum_{m}\rho^{lm;lm}.$ (34)

We thus project the Lindblad equation onto the spherical harmonics and sum
over the magnetic quantum number $m$. After this projection, the Lindblad
equation is an infinite dimensional matrix equation, which is the tensor
product of the $2\times 2$ color space, the infinite dimensional angular
momentum space, and the radial wave function. The density matrix is diagonal
and takes the form

$\rho=\begin{pmatrix}\rho^{0}_{s}&0&\ldots&0&0&\cdots\\\
0&\rho^{1}_{s}&\ldots&0&0&\ldots\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\ldots&\rho^{0}_{o}&0&\ldots\\\ 0&0&\ldots&0&\rho_{o}^{1}&\ldots\\\
\vdots&\vdots&&\vdots&\vdots&\ddots\end{pmatrix},$ (35)

with the Hamiltonian being an analogous diagonal matrix, the elements of which
are given by

$h^{l}_{s,o}=-\frac{1}{M}\left(\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial
r}\right)+V_{s,o}+\frac{l(l+1)}{Mr^{2}}.$ (36)

The collapse operators take the form

$\displaystyle\begin{aligned}
C^{0}=&\sqrt{\frac{\kappa}{N_{c}^{2}-1}}\begin{pmatrix}0&1\\\
0&0\end{pmatrix}\otimes\left(O^{-}\sqrt{\frac{l}{2l+1}}C^{\downarrow}_{o\to
s}+O^{+}\sqrt{\frac{l+1}{2l+1}}C^{\uparrow}_{o\to s}\right)\\\
&+\sqrt{\kappa}\begin{pmatrix}0&0\\\
1&0\end{pmatrix}\otimes\left(O^{-}\sqrt{\frac{l}{2l+1}}C^{\downarrow}_{s\to
o}+O^{+}\sqrt{\frac{l+1}{2l+1}}C^{\uparrow}_{s\to o}\right),\end{aligned}$
(37) $\displaystyle
C^{1}=\sqrt{\frac{\kappa(N_{c}^{2}-4)}{2(N_{c}^{2}-1)}}\begin{pmatrix}0&0\\\
0&1\end{pmatrix}\otimes\left(O^{-}\sqrt{\frac{l}{2l+1}}C^{\downarrow}_{o\to
o}+O^{+}\sqrt{\frac{l+1}{2l+1}}C^{\uparrow}_{o\to o}\right),$ (38)

where the $2\times 2$ matrices are responsible for color transitions and the
infinite dimensional matrices $O^{\pm}_{l^{\prime},l}=\delta_{l^{\prime},l\pm
1}$ for angular momentum transitions. The operators $C_{u\to
v}^{\uparrow,\downarrow}$ act on $\rho^{l}_{u}$ contributing to its partial
width to $\rho^{l\pm 1}_{v}$; they are given explicitly by

$\displaystyle C^{\uparrow}_{u\to v}=$ $\displaystyle r\left(1+\frac{\Delta
V_{uv}}{4T}\right)+\frac{1}{2MT}\left(\frac{\partial}{\partial
r}-\frac{l}{r}\right),$ (39) $\displaystyle C^{\downarrow}_{u\to v}=$
$\displaystyle r\left(1+\frac{\Delta
V_{uv}}{4T}\right)+\frac{1}{2MT}\left(\frac{\partial}{\partial
r}+\frac{l+1}{r}\right).$ (40)

In this way, the 3-dimensional problem has been reduced to a 1-dimensional
problem and the 6 collapse operators to 2. For details, see appendix C.

This same procedure was carried out in ref. Brambilla:2017zei  truncating the
angular momentum at $l=1$. Since then, works based on the QTraj code
Omar:2021kra  have, however, made it possible to account for states of
arbitrarily high angular momentum Brambilla:2020qwo ; Brambilla:2021wkt .

The decay width of the state $\rho^{l}_{u}(r)$ is given by the trace of the
anticommutator term of the Lindblad equation

$\Gamma\left[\rho_{u}^{l}(r)\right]\equiv\sum_{n=0}^{1}\text{Tr}\left[C^{n\,\dagger}C^{n}\rho^{l}_{u}(r)\right],$
(41)

where on the left $\Gamma[\rho_{u}^{l}(r)]$ represents a functional returning
the width and on the right the trace is over color, angular momentum, and
position. Partial widths with respect to each of the quantum numbers are given
by

$\displaystyle\Gamma^{\uparrow,\downarrow}\rho_{u}(r)=$
$\displaystyle\sum_{n=0}^{1}\text{Tr}_{l}\left[C^{n\,\dagger}|l\pm
1\rangle\langle l\pm 1|C^{n}\rho^{l}_{u}(r)\right],$ (42)
$\displaystyle\Gamma_{u\to v}\rho^{l}(r)=$
$\displaystyle\sum_{n=0}^{1}\text{Tr}_{c}\left[C^{n\,\dagger}|v\rangle\langle
v|C^{n}\rho^{l}_{u}(r)\right],$ (43)

where the traces are only over angular momentum and color, respectively, and
$|l\pm 1\rangle\langle l\pm 1|$ and $|v\rangle\langle v|$ are states of
angular momentum $l\pm 1$ and color $v$, respectively. Computing the operators
$\Gamma^{\uparrow,\downarrow}$ and $\Gamma_{u\to v}$ explicitly, we find

$\displaystyle\Gamma^{\uparrow}$ $\displaystyle=\frac{l+1}{2l+1},$ (44)
$\displaystyle\Gamma^{\downarrow}$ $\displaystyle=\frac{l}{2l+1},$ (45)
$\displaystyle\Gamma_{s\to o}$
$\displaystyle=\kappa\left\\{r^{2}\left(1+\frac{\Delta
V_{so}}{4T}\right)^{2}-\frac{3}{2MT}-\frac{\Delta
V_{so}}{4MT^{2}}+\frac{\mathcal{D}^{2}}{4M^{2}T^{2}}\right\\},$ (46)
$\displaystyle\Gamma_{o\to s}$
$\displaystyle=\frac{\kappa}{N_{c}^{2}-1}\left\\{r^{2}\left(1+\frac{\Delta
V_{os}}{4T}\right)^{2}-\frac{3}{2MT}-\frac{\Delta
V_{os}}{4MT^{2}}+\frac{\mathcal{D}^{2}}{4M^{2}T^{2}}\right\\},$ (47)
$\displaystyle\Gamma_{o\to o}$
$\displaystyle=\frac{\kappa(N_{c}^{2}-4)}{2(N_{c}^{2}-1)}\left\\{r^{2}-\frac{3}{2MT}+\frac{\mathcal{D}^{2}}{4M^{2}T^{2}}\right\\},$
(48)

where

$\mathcal{D}^{2}=-\left(\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}\right)+\frac{l(l+1)}{r^{2}}\,.$
(49)

The partial width of the state $\rho^{l}_{u}(r)$ to a state of angular
momentum $l\pm 1$ and color $v$ is given by

$\Gamma^{\uparrow,\downarrow}_{u\to
v}\left[\rho^{l}_{u}(r)\right]\equiv\text{Tr}_{r}\left[\Gamma^{\uparrow,\downarrow}\Gamma_{u\to
v}\rho^{l}_{u}(r)\right],$ (50)

where on the left $\Gamma^{\uparrow,\downarrow}_{u\to v}[\rho^{l}_{u}(r)]$
represents a functional returning the partial width and the trace is over the
radial coordinate $r$.

##  3 The quantum trajectories algorithm at NLO

The Lindblad equation given in eq. (28) with the collapse operators of eqs.
(31) and (32) describes the evolution of the heavy quarkonium reduced density
matrix at next-to-leading-order accuracy in the $E/T$ expansion. To solve it
numerically, we make use of the quantum trajectories algorithm used to solve
the LO evolution equations in refs. Brambilla:2021wkt ; Brambilla:2020qwo ;
Omar:2021kra . In the quantum trajectories approach, one evolves a large set
of independently sampled quantum evolutions of the wave function. Observables,
e.g., the overlap with vacuum eigenstates, are computed along each sampled
quantum trajectory and averaged to obtain the final predictions for the
observable in question.

Here we make use of the “waiting time approach” in which one evolves the wave
function using the complex effective Hamiltonian until the norm squared of the
wave function falls below a random number $p_{1}$ uniformly sampled in
$[0,1]$, at which point additional random numbers are generated to determine
the outgoing quantum numbers and corresponding collapse operator to apply to
the wave function Omar:2021kra . After application of the selected collapse
operator, the wave function is normalized, and the process is repeated until
the simulation end time is reached. This method has the benefit that between
quantum jumps the evolution proceeds with a fixed color state labeled by $c$
and fixed angular momentum labeled by $l$. For the description below, we note
that as the system is isotropic we can write the Hamiltonian and collapse
operators more compactly by using the reduced wave function $u(t,r)\equiv
rR(t,r)$, where the three-dimensional wave function for fixed $l$ and $m$ is
given by $\psi(t,\vec{r})=R(t,r)Y_{lm}(\theta,\phi)$. In subsections 3.1, 3.2,
and 3.3, we give the forms of the collapse operators, widths, and in-medium
Hamiltonian acting on the reduced wave function; we refer to quantities acting
on the reduced wave function as being in the reduced spherical representation
and denote them with an overbar. More details on the reduced spherical
representation are given in appendix C. In subsection 3.4, we give the quantum
trajectories algorithm as implemented to obtain the phenomenological results
in sec. 4. In subsection 3.5, we describe the procedure used to evolve the
wave function between quantum jumps.

###  3.1 Jump operators in the reduced spherical representation

There are six jump operators present in the Lindblad equation. They represent
the six possible physical transitions for the quarkonium state inside the QGP.
A color singlet state can transition into an octet and vice versa; in
addition, the octet state can also transition to another octet state. For each
color transition, we have the two additional possibilities of an angular
momentum jump up $(l\rightarrow l+1$) or down ($l\rightarrow l-1$), thus
totalling six. When expressed as operators acting on the reduced wave
function, the six jump operators are

$\displaystyle\overline{C}^{\uparrow}_{s\to o}$ $\displaystyle=$
$\displaystyle r-\frac{N_{c}\alpha_{\rm
s}}{8T}+\frac{1}{2MT}\left(\frac{\partial}{\partial r}-\frac{l+1}{r}\right),$
(51) $\displaystyle\overline{C}^{\downarrow}_{s\to o}$ $\displaystyle=$
$\displaystyle r-\frac{N_{c}\alpha_{\rm
s}}{8T}+\frac{1}{2MT}\left(\frac{\partial}{\partial r}+\frac{l}{r}\right),$
(52) $\displaystyle\overline{C}^{\uparrow}_{o\to s}$ $\displaystyle=$
$\displaystyle r+\frac{N_{c}\alpha_{\rm
s}}{8T}+\frac{1}{2MT}\left(\frac{\partial}{\partial r}-\frac{l+1}{r}\right),$
(53) $\displaystyle\overline{C}^{\downarrow}_{o\to s}$ $\displaystyle=$
$\displaystyle r+\frac{N_{c}\alpha_{\rm
s}}{8T}+\frac{1}{2MT}\left(\frac{\partial}{\partial r}+\frac{l}{r}\right),$
(54) $\displaystyle\overline{C}^{\uparrow}_{o\to o}$ $\displaystyle=$
$\displaystyle r+\frac{1}{2MT}\left(\frac{\partial}{\partial
r}-\frac{l+1}{r}\right),$ (55) $\displaystyle\overline{C}^{\downarrow}_{o\to
o}$ $\displaystyle=$ $\displaystyle
r+\frac{1}{2MT}\left(\frac{\partial}{\partial r}+\frac{l}{r}\right).$ (56)

Since the wave function is normalized after application, the normalization of
these operators is arbitrary. In the code, we multiply all of these operators
by a factor of $T$ so that they are regular when $T=0$.

###  3.2 Width operators in the reduced spherical representation

When expressed as operators acting on the reduced wave function, the octet-
singlet and octet-octet width operators take the form

$\displaystyle\overline{\Gamma}_{o\to s}^{\uparrow}$
$\displaystyle=\frac{\hat{\kappa}T^{3}}{N_{c}^{2}-1}\frac{l+1}{2l+1}\left[\left(r+\frac{N_{c}\alpha_{\rm
s}}{8T}\right)^{2}-\frac{3}{2MT}+\frac{\overline{{\cal
D}}^{2}}{(2MT)^{2}}-\frac{1}{2MT}\left(\frac{N_{c}\alpha_{\rm
s}}{4T}\right)\frac{1}{r}\right],$ (57) $\displaystyle\overline{\Gamma}_{o\to
s}^{\downarrow}$ $\displaystyle=\frac{l}{l+1}\overline{\Gamma}_{o\to
s}^{\uparrow}\,,$ (58) $\displaystyle\overline{\Gamma}_{o\to o}^{\uparrow}$
$\displaystyle=\hat{\kappa}T^{3}\frac{N_{c}^{2}-4}{2(N_{c}^{2}-1)}\frac{l+1}{2l+1}\left[r^{2}-\frac{3}{2MT}+\frac{\overline{{\cal
D}}^{2}}{(2MT)^{2}}\right]\,,$ (59) $\displaystyle\overline{\Gamma}_{o\to
o}^{\downarrow}$ $\displaystyle=\frac{l}{l+1}\overline{\Gamma}_{o\to
o}^{\uparrow}\,,$ (60)

with

$\displaystyle\overline{{\cal D}}^{2}=-\frac{\partial^{2}}{\partial
r^{2}}+\frac{l(l+1)}{r^{2}}\,.$ (61)

We do not list the singlet-octet width operator because singlet states must
always transition to an octet state and the outgoing angular momentum quantum
number can be chosen based on the fact that, similar to eq. (60), the up and
down transitions are related by $\overline{\Gamma}^{\downarrow}_{s\to
o}/\overline{\Gamma}^{\uparrow}_{s\to o}=l/(l+1)$.

###  3.3 Effective Hamiltonian in the reduced spherical representation

The effective Hamiltonian for singlet and octet evolution is defined by
$H^{\rm eff}_{s,o}=h_{s,o}+\text{Im}(\Sigma_{s,o})-i\Gamma_{s,o}/2$ with
$\Gamma_{s}=\sum_{i\in\\{\uparrow,\downarrow\\}}\Gamma_{s\rightarrow o}^{i}$
and $\Gamma_{o}=\sum_{i\in\\{\uparrow,\downarrow\\}}(\Gamma_{o\rightarrow
s}^{i}+\Gamma_{o\rightarrow o}^{i}$). When expressed as operators acting on
the reduced wave function, the singlet effective Hamiltonian
$\overline{H}^{\rm eff}_{s}$ is given by

$\displaystyle\text{Re}[\overline{H}^{\rm eff}_{s}]$
$\displaystyle=\frac{\overline{\cal D}^{2}}{M}-\frac{C_{f}\,\alpha_{\rm
s}}{r}+\frac{\hat{\gamma}T^{3}}{2}r^{2}+\frac{\hat{\kappa}T^{2}}{4M}\\{r,p_{r}\\}\,,$
(62) $\displaystyle\text{Im}[\overline{H}^{\rm eff}_{s}]$
$\displaystyle=-\frac{\hat{\kappa}T^{3}}{2}\left[\left(r-\frac{N_{c}\alpha_{\rm
s}}{8T}\right)^{2}-\frac{3}{2MT}+\frac{{\overline{\cal
D}}^{2}}{(2MT)^{2}}+\frac{1}{2MT}\left(\frac{N_{c}\alpha_{\rm
s}}{4T}\right)\frac{1}{r}\right],$ (63)

where $p_{r}=-i\partial_{r}$. Similarly, the octet effective Hamiltonian
$\overline{H}^{\rm eff}_{o}$ is given by

$\displaystyle\text{Re}[\overline{H}^{\rm eff}_{o}]$
$\displaystyle=\frac{\overline{\mathcal{D}}^{2}}{M}+\frac{1}{2N_{c}}\frac{\alpha_{\rm
s}}{r}+\frac{N_{c}^{2}-2}{2(N_{c}^{2}-1)}\left[\frac{\hat{\gamma}T^{3}}{2}r^{2}+\frac{\hat{\kappa}T^{2}}{4M}\\{r,p_{r}\\}\right],$
(64) $\displaystyle\text{Im}[\overline{H}^{\rm eff}_{o}]$
$\displaystyle=-\frac{\hat{\kappa}T^{3}}{2(N_{c}^{2}-1)}\left[\left(r+\frac{N_{c}\alpha_{\rm
s}}{8T}\right)^{2}-\frac{3}{2MT}+\frac{{\overline{\cal
D}}^{2}}{(2MT)^{2}}-\frac{1}{2MT}\left(\frac{N_{c}\alpha_{\rm
s}}{4T}\right)\frac{1}{r}\right]$ $\displaystyle~{}\hskip
11.38092pt{-}\frac{\hat{\kappa}T^{3}}{4(N_{c}^{2}-1)}\left[r^{2}-\frac{3}{2MT}+\frac{{\overline{\cal
D}}^{2}}{(2MT)^{2}}\right],$ (65)

where $\hat{\kappa}=\kappa/T^{3}$ and $\hat{\gamma}=\gamma/T^{3}$.

###  3.4 Description of the algorithm

The NLO quantum trajectories algorithm proceeds as follows:

  1. 1.

Initialize the reduced wave function $u(t_{0},r)=u_{0}(r)$ at $N$ points which
are uniformly distributed on $(0,L)$ with $L$ being the box size. The wave
function is assumed to obey Dirichlet boundary conditions at the boundaries
with $u(t,0)=u(t,L)=0$ at all times. Additionally, we specify the initial
color state $c\in\\{0,1\\}$, with 0 being a color singlet and 1 being a color
octet state, and the initial orbital angular momentum quantum number $l\geq
0$.

  2. 2.

Normalize the wave function such that $|u|^{2}=\langle u|u\rangle$ = 1.

  3. 3.

Generate a random number $p_{1}$ uniformly distributed in $[0,1]$.

  4. 4.

Using $\overline{H}^{\rm eff}_{s,o}$ specified in eqs. (62)-(65), evolve the
wave function in time until the jump time $t_{j}$. A quantum jump is triggered
by the norm squared of the wave function falling below $p_{1}$ at time
$t_{j}$, i.e., $|u(t_{j},r)|^{2}<p_{1}$. If the maximum evolution time is
reached before a jump is triggered, i.e., $t_{j}>t_{f}$, terminate execution
at time $t_{f}$, otherwise, perform a quantum jump at time $t_{j}$ by
proceeding to step 5.

  5. 5.

Initiate a quantum jump by generating two additional random numbers $p_{2}$
and $p_{3}$ uniformly distributed in $[0,1]$. These numbers will be used to
select new angular momentum and color states.

  6. 6.

If the color state is singlet ($c=0$):

    1. (a)

Set the color state to octet, i.e. $c\rightarrow 1$.

    2. (b)

If $p_{2}<{l}/{(2l+1)}$, $l\rightarrow l-1$, else $l\rightarrow l+1$.

Else ($c=1$):

    1. (a)

Evaluate the transition amplitudes
$\overline{\Gamma}_{1}=\langle\overline{\Gamma}_{o\to o}^{\downarrow}\rangle$,
$\overline{\Gamma}_{2}=\langle\overline{\Gamma}_{o\to o}^{\uparrow}\rangle$,
$\overline{\Gamma}_{3}=\langle\overline{\Gamma}_{o\to s}^{\downarrow}\rangle$,
and $\overline{\Gamma}_{4}=\langle\overline{\Gamma}_{o\to
s}^{\uparrow}\rangle$ using the operators listed in eqs. (57)-(60) above. From
these, obtain the normalized amplitudes
$\tilde{\Gamma}_{i}=\overline{\Gamma}_{i}/\sum_{j=1}^{4}\overline{\Gamma}_{j}$.
The total probability for a downward angular momentum transition is given by
$p_{\downarrow}=\tilde{\Gamma}_{1}+\tilde{\Gamma}_{3}$, and the total
probability for an octet-singlet transition is given by $p_{\rm o\to
s}=\tilde{\Gamma}_{3}+\tilde{\Gamma}_{4}$.

    2. (b)

If $p_{2}<p_{\downarrow}$, $l\rightarrow l-1$, else $l\rightarrow l+1$.

    3. (c)

If $p_{3}<p_{\rm o\to s}$, $c\rightarrow 0$.

  7. 7.

Based on the choices made above, apply the appropriate collapse/jump operator
to the wavefunction $u(r)\rightarrow\overline{C}u(r)$ with
$\overline{C}\in\\{\overline{C}_{s\to o}^{\downarrow},\overline{C}_{s\to
o}^{\uparrow},\overline{C}_{o\to s}^{\downarrow},\overline{C}_{o\to
s}^{\uparrow},\overline{C}_{o\to o}^{\downarrow},\overline{C}_{o\to
o}^{\uparrow}\\}$ given in eqs. (51)-(56) above.

  8. 8.

Repeat starting from step 2.

###  3.5 Evolution between jumps

In step 4 above, we must evolve the wave function forward in time using the
complex effective Hamiltonian $\overline{H}^{\rm eff}_{s,o}$. Because of the
appearance of the anticommutator term $\\{r,p_{r}\\}$, it is not easy to
straightforwardly apply the same split-step pseudospectral method
(Suzuki–Trotter) as was used in ref. Omar:2021kra . Instead, here we use the
Crank–Nicolson method to update the wave function between quantum jumps. With
this method, one approximates the infinitesimal time evolution operator as

$e^{-i\overline{H}^{\rm eff}_{s,o}\Delta t}\simeq\frac{1-i\Delta
t\overline{H}^{\rm eff}_{s,o}/2}{1+i\Delta t\overline{H}^{\rm
eff}_{s,o}/2}\,,$ (66)

where the Hamiltonian has implicit time dependence. To proceed, one expresses
this as

$\left(1+\frac{i}{2}\Delta t\overline{H}^{\rm eff}_{s,o}\right)\psi(t+\Delta
t)=\left(1-\frac{i}{2}\Delta t\overline{H}^{\rm eff}_{s,o}\right)\psi(t)\,.$
(67)

To evolve the wave function one step forward in time, one applies the terms on
the right to obtain the “half updated” wave function $X\equiv\left(1-i\Delta
t\overline{H}^{\rm eff}_{s,o}/2\right)\psi(t)$. The final step in the
Crank–Nicolson algorithm is to solve

$A_{ij}\psi_{j}(t+\Delta t)=X_{i}\,,$ (68)

for $\psi_{j}(t+\Delta t)$ where $i,j$ are spatial coordinates and
$A_{ij}=\left(1+i\Delta t\overline{H}^{\rm eff}_{s,o}/2\right)_{ij}$.
Typically $A_{ij}$ is tridiagonal, in which case there are many optimized
solvers. For the results below, we make use of the optimized sparse matrix
solver spsolve provided by the open-source Armadillo package Armadillo . The
Armadillo package is also used to handle all matrix multiplications required
to compute the action of the effective Hamiltonian, width, and jump operators.

##  4 Results

We consider $5.02$ TeV Pb-Pb collisions with the background temperature
evolution given by 3+1D quasiparticle anisotropic hydrodynamics which was
tuned to reproduce experimentally observed soft hadron spectra, elliptic flow,
and HBT radii Alqahtani:2020paa ; Alalawi:2021jwn . For this purpose, smooth
optical Glauber initial conditions were used, and the resulting initial
central temperature was $T_{0}=$ 630 MeV at $\tau_{0}=0.25$ fm and the shear
viscosity to entropy density ratio was $\eta/s=0.159$. All other transport
coefficients, such as the bulk viscosity, were self-consistently determined in
terms of $\eta/s$ and the lattice-based equation of state used in the
evolution. The hydrodynamic evolution used is the same as that used in the
prior papers Brambilla:2020qwo ; Brambilla:2021wkt .

For all results presented herein, we used a one-dimensional lattice with
$\texttt{NUM}=2048$ points and $\texttt{L}=40\,\mathrm{GeV}^{-1}$ with a
corresponding lattice spacing of $a\simeq 0.0195\,\mathrm{GeV}^{-1}$. The
Crank–Nicolson evolution time step was taken to be
$\texttt{dt}=0.001\,\mathrm{GeV}^{-1}$. This spatiotemporal discretization was
shown in prior work to have only small lattice spacing and finite size effects
Omar:2021kra . Finally, we note that the code used to generate all results in
this paper is available publicly in the NLO branch of the QTraj code
repository qtraj-download .

As emphasized in the introduction, the parameters entering the evolution
equations do not need to be fit to experimental data. These are the heavy
quark mass $M$, the strong coupling $\alpha_{\rm s}$, and the transport
coefficients $\kappa$ and $\gamma$. We work in the bottom sector and, as such,
take $M=m_{b}=m_{\Upsilon(1S)}/2=4.73$ GeV with $m_{\Upsilon(1S)}$ from ref.
Zyla:2020zbs . We compute the Bohr radius $a_{0}$ from the 1-loop relation

$a_{0}=\frac{2}{C_{f}\,\alpha_{\rm s}(1/a_{0})\,m_{b}},$ (69)

and evaluate $\alpha_{\rm s}$ at the inverse of the Bohr radius using the
1-loop running with $N_{f}=3$ flavors and $\Lambda_{\overline{\rm
MS}}^{N_{f}=3}=332$ MeV Zyla:2020zbs  obtaining $\alpha_{\rm
s}(1/a_{0})=0.468$.

The transport coefficients $\kappa$ and $\gamma$ are fixed from direct and
indirect lattice measurements, respectively. In ref. Brambilla:2020siz , the
heavy quark momentum diffusion coefficient CasalderreySolana:2006rq ;
CaronHuot:2007gq  was measured directly in a quenched lattice simulation over
a wide range of temperatures. Here we identify $\kappa(T)$ with this lattice
determination of the quark momentum diffusion coefficient, while we postpone a
discussion on this identification and its associated uncertainty to appendix
D. We perform QTraj simulations using three parametrizations of the
dimensionless quantity $\hat{\kappa}(T)=\kappa/T^{3}$ which we denote
$\hat{\kappa}_{L}(T)$, $\hat{\kappa}_{C}(T)$, and $\hat{\kappa}_{U}(T)$ and
which correspond to the lower, central, and upper “fit” curves, respectively,
of fig. 13 of ref. Brambilla:2020siz ; we use this variation as an estimate of
our systematic uncertainty due to $\kappa$. As no direct lattice measurements
of $\gamma$ currently exist, we make use of an indirect estimate from
unquenched lattice simulations. At zeroth order in the $E/T$ expansion,
projection of the in-medium heavy quarkonium self energy onto the $1S$ state
yields the relation $\delta M(1S)=(3/2)a_{0}^{2}\gamma$, where $\delta M(1S)$
is the in-medium mass shift of the $1S$ state, cf. eq. (78) of ref.
Brambilla:2017zei  and ref. Brambilla:2019tpt . From this relation, $\gamma$
is accessible from unquenched lattice measurements of $\delta M(1S)$. In ref.
Brambilla:2019tpt , unquenched lattice measurements of $\delta
M[\Upsilon(1S)]$ from refs. Kim:2018yhk ; Aarts:2011sm  were used to place
bounds on $\hat{\gamma}(T)=\gamma/T^{3}$ of approximately
$-3.5\leq\hat{\gamma}\leq 0$. We perform simulations using $\hat{\gamma}$ in
this range and use this as an estimation of our systematic uncertainty due to
$\gamma$. We note that more recent lattice studies Larsen:2019bwy ;
Shi:2021qri  seem to favor $\delta M(\Upsilon(1S))\simeq 0$ and thus
$\hat{\gamma}\simeq 0$.

To initialize the simulation, we assume that at $\tau=0$ fm the bottomonium
wave function is in a singlet state with a smeared delta function initial
condition.333One can also consider octet initialization; however, we find that
the off-diagonal octet-singlet overlap is negligible, similar to the finding
at leading order Brambilla:2020qwo . In practice, the initial reduced wave
function is given by a Gaussian delta function multiplied by a power of $r$
appropriate for the initial angular momentum state $l$, i.e.,

$u_{\ell}(t_{0})\propto r^{l+1}e^{-r^{2}/(ca_{0})^{2}},$ (70)

with $u$ normalized to one when summed over the entire (one-dimensional)
lattice volume. We set $c=0.2$ noting that observables do not appear to show
significant dependence on this parameter for $c$ below this value (cf. fig. 6
of ref. Omar:2021kra ) while the computational cost increases dramatically as
$c$ is decreased. For the results presented in the main body of the text, we
evolve the initial wave function using the vacuum potential
($\hat{\kappa}=\hat{\gamma}=0$) from $\tau=0$ fm to $\tau_{\rm med}$ = 0.6 fm
at which time we turn on the medium interactions. This is the same medium
initialization time scale as used in refs. Brambilla:2020qwo ;
Brambilla:2021wkt . To assess the dependence on this assumption, in appendix
E, we present results obtained with $\tau_{\rm med}$ = 0.25 fm, which
corresponds to the earliest time for which the anisotropic hydrodynamic
simulation results are available.

###  4.1 Comparison of LO and NLO singlet-octet widths

In this subsection, we analyze the effect of the inclusion of higher order
terms in the $E/T$ expansion. Projecting the singlet to octet width
$\Gamma_{s\to o}$ given in eq. (46) onto $1S$, $2S$, and $3S$ states, we find

$\displaystyle\langle 1S|\Gamma_{s\to o}|1S\rangle$
$\displaystyle=3a_{0}^{2}\kappa\left\\{1-\frac{2N_{c}^{2}-1}{2(N_{c}^{2}-1)}\frac{E}{T}+\frac{(2N_{c}^{2}-1)^{2}}{12(N_{c}^{2}-1)^{2}}\left(\frac{E}{T}\right)^{2}\right\\},$
(71) $\displaystyle\langle 2S|\Gamma_{s\to o}|2S\rangle$
$\displaystyle=42a_{0}^{2}\kappa\left\\{1-\frac{5N_{c}^{2}-1}{28(N_{c}^{2}-1)}\frac{E}{T}+\frac{7N_{c}^{4}-4N_{c}^{2}+1}{672(N_{c}^{2}-1)^{2}}\left(\frac{E}{T}\right)^{2}\right\\},$
(72) $\displaystyle\langle 3S|\Gamma_{s\to o}|3S\rangle$
$\displaystyle=207a_{0}^{2}\kappa\left\\{1-\frac{10N_{c}^{2}-1}{138(N_{c}^{2}-1)}\frac{E}{T}+\frac{12N_{c}^{4}-4N_{c}^{2}+1}{7452(N_{c}^{2}-1)^{2}}\left(\frac{E}{T}\right)^{2}\right\\},$
(73)

where $E=1/(Ma_{0}^{2})$ is the magnitude of the Coulombic binding energy; the
above expressions make explicit the $E/T$ expansion. We plot the various
contributions to the bottomonium widths in fig. 1. Reading the plot from left
to right, we observe the expansion to converge for high $T$ and, comparing the
solid to dotted curves, upon inclusion of additional terms in the $E/T$
expansion. Furthermore, comparing the blue to orange to green curves, we
observe better convergence properties for states of higher principal quantum
due to their lower Coulombic binding energy.

Figure 1:  The singlet to octet widths of the 1S, 2S, and 3S states including
$E/T$ corrections. $\Gamma_{m}(nS)$ represents the singlet to octet width of
the state $nS$ including terms up to order $(E/T)^{m}$. The black line at
unity represents the limit of perfect convergence of the $E/T$ expansion,
i.e., $E/T\to 0$. The plot spans from $T=158$ MeV to $T=600$ MeV; the gray
shaded area on the left represents the temperature region $T<T_{f}=190$ MeV
not included in the phenomenological results presented in this work.

###  4.2 Comparisons of LO and NLO survival probabilities with complex $H_{\rm
eff}$

In this subsection, we investigate the effect of the NLO evolution terms on
the survival probabilities of the $\Upsilon(1S)$ state. To model the medium
evolution, we follow the procedure of ref. Brambilla:2020qwo  and sample a
large number of physical bottomonium trajectories through the QGP recording
the temperature along each trajectory. We bin the events in centrality
(0-100%) and average over the physical trajectories in a given centrality bin
to obtain a trajectory-averaged temperature evolution in each bin; plots of
the temperature evolution in each centrality bin considered can be found in
fig. 4 of ref. Brambilla:2020qwo .

Figure 2:  The survival probability as a function of $N_{\text{part}}$ of the
$1S$ state calculated using $H_{\rm eff}$ evolution without jumps. In the top
row, we compare results obtained using LO and NLO evolution both with
$T_{f}=190$ MeV to highlight the effect of the inclusion of higher order terms
in the $E/T$ expansion. In the bottom row, we compare our current state of the
art NLO evolution with $T_{f}=190$ MeV to LO evolution with $T_{f}=250$ MeV as
used in previous works. The bands indicate variation with respect to
$\hat{\kappa}(T)$ (left) and $\hat{\gamma}$ (right). The central curves
represent the central values of $\hat{\kappa}(T)$ and $\hat{\gamma}$, and the
dashed and dot-dashed lines represent the lower and upper values,
respectively, of $\hat{\kappa}(T)$ and $\hat{\gamma}$.

Using this trajectory averaged temperature evolution, we perform QTraj
simulations including only the LO and the NLO terms of the evolution. For the
LO evolution, we run simulations down to $T_{f}=190$ MeV and $T_{f}=250$ MeV;
for the NLO evolution, we use $T_{f}=190$ MeV. We plot these results in fig.
2. Comparing the LO and NLO results both with $T_{f}=190$ MeV, we observe the
former to fall below the latter for all but the most peripheral events where
very little hydrodynamic evolution takes place. Comparing the LO results with
$T_{f}=250$ MeV to the NLO results with $T_{f}=190$ MeV, we observe the NLO
results to fall above the LO results in central collisions and below them in
peripheral collisions. There are two competing trends at work here: that of
the NLO terms to reduce suppression and that of lower $T_{f}$ to increase it.
The former dominates in central collisions with greater total time spent
traversing the QGP while the latter dominates in peripheral collisions for
which a high value of $T_{f}$ can lead to no hydrodynamic evolution due to low
initial temperatures.

Figure 3: Jumps vs. no jumps. Top row is $\hat{\kappa}$ variation; bottom row
is $\hat{\gamma}$ variation. Solid, dashed, and dot-dashed curves represent
the central, lower, and upper values, respectively, of $\hat{\kappa}$ and
$\hat{\gamma}$.

###  4.3 Assessing the effects of quantum jumps on the survival probabilities

In this subsection, we assess the effects of quantum jumps on our NLO results.
We make use of the same trajectory-averaged temperature profile as in the
previous subsection. In fig. 3, we present results for the survival
probability of the $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$ states
as functions of centrality. For the results including full evolution with
stochastically sampled jumps, we report only the central (mean) values and
note that the statistical uncertainties are sub-leading when compared to the
uncertainty resulting from $\kappa$ and $\gamma$ variation. For identical
parameters, we observe the jumps to increase the yield of the states. The
relative effect appears to become more pronounced for the excited states. We
note, however, that the low yields of the excited states make this difference
small in absolute terms. As in the leading order case explored in refs.
Brambilla:2020qwo ; Brambilla:2021wkt , we find that at NLO the effect of the
stochastic jumps is much smaller than the theoretical uncertainties stemming
from $\kappa$ and $\gamma$ variation. For $R_{AA}$, the evolution with
$H_{\text{eff}}$ provides a reasonable estimation at a significantly reduced
computational cost.

###  4.4 NLO $H_{\rm eff}$ results for $R_{AA}$ with sampled physical
trajectories

In this subsection we present results of $H_{\rm eff}$ evolution using an
ensemble of 80,000 sampled physical trajectories for each set of values for
$\hat{\kappa}(T)$ and $\hat{\gamma}$. Unlike the results presented in the
previous two subsections, here we compute the survival probability on a
trajectory-by-trajectory basis and average the resulting survival
probabilities rather than using a single trajectory-averaged temperature
evolution in each centrality bin. Using individually sampled physical
trajectories allows us to more faithfully include the physics of states
produced at a large distance from the center of the overlap region (the corona
effect). This is particularly important for quantitatively understanding
excited state survival probabilities, since the corona region contributes a
large fraction of surviving states (see, e.g., fig. 3 of ref. Islam:2020bnp ).

#### Excited state feed down

The Lindblad equation, which the QTraj code solves, gives the probability of a
heavy quark-antiquark pair exiting the QGP with specific quantum numbers. In
order to compare to experimental measurements of $R_{AA}$, this raw survival
probability data must be post processed to account for the probability of a
state decaying after exiting the QGP but before reaching the detector. This
effect is denoted feed down and can be accounted for by using the feed down
matrix $F$ defined as relating the experimental and direct production cross
sections, $\vec{\sigma}_{\text{exp}}=F\vec{\sigma}_{\text{direct}}$. The cross
section vectors correspond to the states considered, while $F$ is a matrix the
values of which are fixed by the branching fractions of the excited states. In
our analysis, the states considered are $\vec{\sigma}=\\{\Upsilon(1S),\,$
$\Upsilon(2S),\,$ $\chi_{b0}(1P),\,$ $\chi_{b1}(1P),\,$ $\chi_{b2}(1P),\,$
$\Upsilon(3S),\,$ $\chi_{b0}(2P),\,$ $\chi_{b1}(2P),\,$ $\chi_{b2}(2P)\\}$.
The entries of $F$ are

$F_{ij}=\left\\{\begin{matrix}\text{branching fraction $j$ to $i$},&\text{for
}i<j,\\\ 1,&\text{for }i=j,\\\ 0,&\text{for }i>j,\end{matrix}\right.$ (74)

where the branching fractions are taken from the Particle Data Group
Zyla:2020zbs  (see also eq. (6.4) of ref. Brambilla:2020qwo ).

Finally, the nuclear modification factor $R_{AA}^{i}$ for state $i$ is given
by

$R^{i}_{AA}(c,p_{T},\phi)=\frac{\left(F\cdot
S(c,p_{T},\phi)\cdot\vec{\sigma}_{\text{direct}}\right)^{i}}{\vec{\sigma}_{\text{exp}}^{i}}\,,$
(75)

where $S(c,p_{T},\phi)$ represents the survival probabilities computed from
the QTraj code in the form of a diagonal matrix; $c$ labels the centrality
class, $p_{T}$ the transverse momentum, and $\phi$ the azimuthal angle. We use
the integrated experimental cross sections
$\vec{\sigma}_{\text{exp}}=\\{57.6$, 19, 3.72, 13.69, 16.1, 6.8, 3.27, 12.0,
$14.15\\}$ nb, computed from the measurements of refs. Sirunyan:2018nsz ;
Aaij:2014caa  as explained in sec. 6.4 of ref. Brambilla:2020qwo .

Figure 4:  $R_{AA}$ for the $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$
as a function of $N_{\rm part}$. The left panel shows variation of
$\hat{\kappa}\in\\{\kappa_{L}(T),\kappa_{C}(T),\kappa_{U}(T)\\}$ and the right
panel shows variation of $\hat{\gamma}$ in the range $-3.5\leq\hat{\gamma}\leq
0$. In both panels, the solid line corresponds to
$\hat{\kappa}=\hat{\kappa}_{C}(T)$ and the best fit value of
$\hat{\gamma}=-2.6$. The experimental measurements shown are from the ALICE
Acharya:2020kls , ATLAS ATLAS5TeV , and CMS Sirunyan:2018nsz ;
CMSupsilonQM2022  collaborations.

Figure 5:  $R_{AA}$ for the $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$
as a function of $p_{T}$. The bands and experimental data sources are the same
as fig. 4.

#### NLO predictions for $R_{AA}$ using $H_{\rm eff}$ evolution

We now turn to our results for $R_{AA}$ of the 1S, 2S, and 3S states. In figs.
4 and 5, we present our NLO predictions for $R_{AA}$ as a function of $N_{\rm
part}$ and $p_{T}$, respectively. The results shown in both figures include
the effect of excited state feed down using the method described in the
previous subsection, and, as in previous sections, we take the decoupling
temperature to be $T_{f}$ = 190 MeV and $\tau_{\rm med{}}=0.6$ fm.444We
consider the case of $T_{f}$ = 190 MeV and $\tau_{\rm med}$ = 0.25 fm in
appendix E. For the results presented in figs. 4 and 5, we do not include the
effect of quantum jumps due to the high numerical demand when including them;
quantum jumps require approximately 50-100x more computational resources than
using only the complex effective Hamiltonian evolution. Based on the results
presented and discussed in sec. 4.3, we expect that using $H_{\rm eff}$
without quantum jumps is a reliable approximation for the suppression of the
$\Upsilon(1S)$; however, we expect this approximation to over-predict the
amount of excited state suppression compared to the complete solution of the
Lindblad equation.

In the left panel of fig. 4, we show the variation of $\hat{\kappa}$ in the
set
$\hat{\kappa}\in\\{\hat{\kappa}_{L}(T),\hat{\kappa}_{C}(T),\hat{\kappa}_{U}(T)\\}$
while holding $\hat{\gamma}=-2.6$. This value of $\hat{\gamma}$ was chosen as
to best reproduce the collected $\Upsilon(1S)$ suppression data. In the right
panel of fig. 4, we show the variation of $\hat{\gamma}$ in the range
$-3.5\leq\hat{\gamma}\leq 0$. As in the left panel, the solid line corresponds
to $\hat{\gamma}=-2.6$. In both the left and right panels, the experimental
data presented are from the ALICE Acharya:2020kls , ATLAS ATLAS5TeV , and CMS
Sirunyan:2018nsz ; CMSupsilonQM2022  collaborations.

Compared to the LO results presented in ref. Brambilla:2021wkt , we find that
the variation of $R_{AA}$ with $\hat{\kappa}$ is reduced, while the variation
with $\hat{\gamma}$ is increased. The former is due to the fact that, for
fixed $\hat{\kappa}$, going from LO to NLO reduces the singlet to octet decay
width. As a result, varying $\hat{\kappa}$ over a fixed range results in a
smaller variation of the dimensionful singlet to octet decay width. Regarding
the $\hat{\gamma}$ variation, it is not entirely clear to us why the $R_{AA}$
variation becomes larger when going from LO to NLO, however, this indicates
that there is a higher degree of dynamical quantum mixing between states at
NLO during the real-time evolution. Apart from these changes in the
sensitivity to $\hat{\kappa}$ and $\hat{\gamma}$, we find that our NLO
predictions with $T_{f}=$190 MeV and our LO predictions with $T_{f}=250$ MeV
presented in ref. Brambilla:2021wkt  are in good agreement, with the NLO
result being in slightly better agreement with data, particularly at low
$N_{\rm part}$. In both figs. 4 and 5, we find excellent agreement between our
NLO $H_{\rm eff}$ predictions and the experimental data for the $R_{AA}[1S]$,
however, for the 2S and 3S excited states, the $H_{\rm eff}$ predictions are
somewhat lower than the experimental results particularly for the most central
collisions. As demonstrated in fig. 5, when integrated over centrality, the
$p_{T}$-dependence of $R_{AA}$ is better reproduced.

Looking forward, the underpredictions of $R_{AA}[2S]$ and $R_{AA}[3S]$
indicate that it is perhaps necessary to include the effect of quantum jumps
on the excited bottomonium states in order to draw firm conclusions. As
discussed in sec. 4.3, the inclusion of quantum jumps does not strongly affect
the $1S$ survival probability, whereas it increases both the $2S$ and $3S$
survival probabilities relative to the $H_{\rm eff}$ (no jump) evolution. We
plan to present results including the effects of quantum jumps in a followup
paper.

##  5 Conclusions and outlook

In the paper, we have studied bottomonium suppression by combining pNRQCD with
an open quantum system framework. Working in the regime $Mv\gg T\gg E$, we
have gone beyond prior studies Brambilla:2020qwo ; Brambilla:2021wkt  by
deriving and numerically solving a Lindblad-type evolution equation that is
accurate to NLO in $E/T$. The goal of extending prior calculations to NLO in
$E/T$ was (1) to allow the pNRQCD+OQS framework to be applied at lower
temperatures than is permitted by a LO truncation and (2) to include terms
which are necessary to describe the approach of the system to equilibrium at
late times.

On the first point, we have shown that when going from LO to NLO there is a
sizable correction to the singlet decay width for $1S$ bottomonium and smaller
corrections for the $2S$ and $3S$ bottomonium states. The size of the
corrections for the $1S$, in particular, stems from the fact that for the $1S$
state $E/T$ is not small at phenomenologically relevant temperatures. Based on
our findings in sec. 4.1, the inclusion of NLO corrections has allowed us to
extend the pNRQCD+OQS treatment down to temperatures in the vicinity of the
QCD phase transition; in practice, we have lowered the decoupling temperature
from $T_{f}$ = 250 MeV at LO to $T_{f}$ = 190 MeV at NLO. We have found that
including the NLO corrections improves the description of $R_{AA}[1S]$ when
compared to experimental data, particularly at low values of $N_{\rm part}$,
which correspond to peripheral and semi-peripheral collisions where the plasma
temperature generated is low. In terms of sensitivity to the parameters, at
NLO we have found that the pNRQCD+OQS predictions for $R_{AA}[1S]$ are less
sensitive to $\kappa$ since the magnitude of the decay width decreases when
going from LO to NLO; however, the sensitivity to $\gamma$ is slightly
increased even though the NLO corrections do not affect the mass shift.555 We
note that
$\left\langle\psi\left|\left\\{r_{i},p_{i}\right\\}\right|\psi\right\rangle=0$
for $\psi^{*}(r)=\psi(r)$ as is the case for eigenstates of the vacuum
Hamiltonian in the reduced spherical representation.  As a result, the
increased sensitivity must be due to increased quantum state mixing during the
real-time evolution.

Regarding the question of thermalization, qualitatively Abelian computations
lead us to expect that the inclusion of NLO corrections in $E/T$ improves the
approach to thermal equilibrium of the master equation Blaizot:2017ypk .
Checking that the inclusion of the NLO corrections leads to a thermal ensemble
of bottomonium states at late times is beyond the scope of this work since
this requires the solution of the Lindblad equation including jumps to be
carried out on a much longer time scale than presented here. This is due to
the rather slow rate of equilibration for bottomonium states as seen in prior
studies Akamatsu:2018xim ; Miura:2019ssi . However, we do observe that NLO
corrections make the $\Upsilon(1S)$ population more stable. This is compatible
with the fact that finite energy effects included at NLO make transitions that
liberate energy more likely than those that absorb energy.

###### Acknowledgements.

N.B., P.V.G. and A.V. acknowledge support by the DFG cluster of excellence
ORIGINS funded by the Deutsche Forschungsgemeinschaft under Germany’s
Excellence Strategy - EXC-2094-390783311. This work has also been supported by
European Research Council project ERC-2018-ADG-835105 YoctoLHC; by Maria de
Maetzu excellence program under project CEX2020-001035-M; by Spanish Research
State Agency under project PID2020-119632GB-I00; and by Xunta de Galicia
(Centro singular de investigación de Galicia accreditation 2019-2022), by
European Union ERDF. M.S., A.I., and A.T. were supported by the U.S.
Department of Energy, Office of Science, Office of Nuclear Physics Award No.
DE-SC0013470. M.S. also thanks the Ohio Supercomputer Center for support under
the auspices of Project No. PGS0253.

##  Appendix A Correlator identities and transport coefficients

In this appendix, we derive the correlator identities necessary to write the
evolution equations at NLO in $E/T$ in terms of the transport coefficients
$\kappa$ and $\gamma$. The transport coefficient $\kappa$ is defined in terms
of the chromoelectric correlator as

$\kappa=\frac{g^{2}}{6N_{c}}\text{Re}\int^{+\infty}_{-\infty}\text{d}t\,\Big{\langle}T\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}=\frac{g^{2}}{6N_{c}}\int^{+\infty}_{0}\text{d}t\,\Big{\langle}\left\\{\tilde{E}^{a}_{i}(t,\vec{0}),\tilde{E}^{a}_{i}(0,\vec{0})\right\\}\Big{\rangle}.$
(76)

In ref. Brambilla:2017zei , the dispersive counterpart of $\kappa$, denoted
$\gamma$, was first defined as

$\gamma=\frac{g^{2}}{6N_{c}}\text{Im}\int^{+\infty}_{-\infty}\text{d}t\,\Big{\langle}T\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}=-i\frac{g^{2}}{6N_{c}}\int^{+\infty}_{0}\text{d}t\,\Big{\langle}\left[\tilde{E}^{a}_{i}(t,\vec{0}),\tilde{E}^{a}_{i}(0,\vec{0})\right]\Big{\rangle},$
(77)

where the $T$ appearing in the correlators expressions represents the time
ordering operator.

Integrals over in-medium correlators at LO in $E/T$ can be written in terms of
these two transport coefficients, since

$\frac{g^{2}}{6N_{c}}\int_{0}^{\infty}\text{d}t\,\Big{\langle}\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}=\frac{1}{2}\left(\kappa+i\gamma\right).$
(78)

For the NLO terms, we have to calculate also integrals with an additional
factor of $t$ in the integrand. In order to compute them, we start by
approximating Akamatsu:2012vt ; Blaizot:2017ypk 666 This amounts at expanding
the chromoelectric correlator around its instantaneous limit.

$i\int_{0}^{\infty}\text{d}t\,t\Big{\langle}\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}\simeq\frac{1}{2}\left.\frac{\text{d}D^{>}(\omega)}{\text{d}\omega}\right|_{\omega=0},$
(79)

where $D^{>}(\omega)$ is the Fourier transform of the chromoelectric
correlator,

$\int_{-\infty}^{\infty}\text{d}t\,e^{i\omega
t}\Big{\langle}\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}=D^{>}(\omega).$
(80)

Then we define the finite temperature correlators

$\displaystyle D^{>}(t,t^{\prime})$
$\displaystyle=\Big{\langle}\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(t^{\prime},\vec{0})\Big{\rangle}$
$\displaystyle=\frac{1}{Z(T)}\text{Tr}\left(\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(t^{\prime},\vec{0})e^{-H/T}\right),$
(81) $\displaystyle D^{<}(t,t^{\prime})$
$\displaystyle=\Big{\langle}\tilde{E}^{a}_{i}(t^{\prime},\vec{0})\tilde{E}^{a}_{i}(t,\vec{0})\Big{\rangle}$
$\displaystyle=\frac{1}{Z(T)}\text{Tr}\left(\tilde{E}^{a}_{i}(t^{\prime},\vec{0})\tilde{E}^{a}_{i}(t,\vec{0})e^{-H/T}\right),$
(82)

where $H$ is the Hamiltonian, and $Z(T)$ is the partition function defined by

$Z(T)=\text{Tr}\left(e^{-H/T}\right).$ (83)

That $e^{-H/T}$ is the time evolution operator in imaginary time, i.e.,
$e^{-H/T}E^{a}_{i}(t,\vec{0})e^{H/T}=E^{a}_{i}(t+i/T,\vec{0})$, implies the
Kubo–Martin–Schwinger relation

$D^{>}(t,t^{\prime})=D^{<}(t+i/T,t^{\prime}).$ (84)

Assuming $D^{>}(t,t^{\prime})=D^{>}(t-t^{\prime})$,
$D^{<}(t,t^{\prime})=D^{<}(t-t^{\prime})$, and setting $t^{\prime}=0$, we find

$\displaystyle D^{>}(t)=D^{<}(-t)=D^{<}(t+i/T)$ $\displaystyle\Rightarrow$
$\displaystyle\int_{-\infty}^{\infty}\text{d}t\,e^{i\omega
t}D^{>}(t)=\int_{-\infty}^{\infty}\text{d}t\,e^{i\omega
t}D^{<}(-t)=\int_{-\infty}^{\infty}\text{d}t\,e^{i\omega t}D^{<}(t+i/T)$
$\displaystyle\Rightarrow$ $\displaystyle
D^{>}(\omega)=D^{<}(-\omega)=e^{\omega/T}D^{<}(\omega)$
$\displaystyle\Rightarrow$
$\displaystyle\left.\frac{\text{d}D^{>}(\omega)}{\text{d}\omega}\right|_{\omega=0}=\left.-\frac{\text{d}D^{<}(\omega)}{\text{d}\omega}\right|_{\omega=0}=\frac{1}{T}D^{<}(\omega=0)+\left.\frac{\text{d}D^{<}(\omega)}{\text{d}\omega}\right|_{\omega=0}$
$\displaystyle\Rightarrow$
$\displaystyle\left.\frac{\text{d}D^{>}(\omega)}{\text{d}\omega}\right|_{\omega=0}=\left.-\frac{\text{d}D^{<}(\omega)}{\text{d}\omega}\right|_{\omega=0}=\frac{1}{2T}D^{<}(\omega=0).$
(85)

Inserting the last equality into eq. (79), we have

$i\int_{0}^{\infty}\text{d}t\,t\Big{\langle}\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}=\frac{1}{4T}D^{<}(\omega=0).$
(86)

Taking $\kappa$ as written in eq. (76), we can bring it to the form

$\displaystyle\kappa$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{6N_{c}}\int^{+\infty}_{0}\text{d}t\,\Big{\langle}\left\\{\tilde{E}^{a}_{i}(t,\vec{0}),\tilde{E}^{a}_{i}(0,\vec{0})\right\\}\Big{\rangle}=\frac{g^{2}}{6N_{c}}\int^{+\infty}_{0}\text{d}t\,\left(D^{>}(t)+D^{<}(t)\right)$
(87) $\displaystyle=$
$\displaystyle\frac{g^{2}}{6N_{c}}\int^{0}_{-\infty}\text{d}t\,D^{<}(t)+\frac{g^{2}}{6N_{c}}\int^{+\infty}_{0}\text{d}t\,D^{<}(t)=\frac{g^{2}}{6N_{c}}\int^{\infty}_{-\infty}\text{d}t\,D^{<}(t)$
$\displaystyle=$ $\displaystyle\frac{g^{2}}{6N_{c}}D^{<}(\omega=0).$

Finally, from eq. (86), we obtain an expression for the integral of the
chromoelectric correlator times one power of $t$ in terms of the transport
coefficient $\kappa$

$i\frac{g^{2}}{6N_{c}}\int_{0}^{\infty}\text{d}t\,t\Big{\langle}\tilde{E}^{a}_{i}(t,\vec{0})\tilde{E}^{a}_{i}(0,\vec{0})\Big{\rangle}=\frac{\kappa}{4T}.$
(88)

##  Appendix B Conditions for obtaining a Lindblad equation

We aim to write a Lindblad equation equivalent to the master equation given in
eq. (11). There are three distinct vector spaces labeling color, spatial
direction, and the transition operator. The color vector space labels the
entries of the transition operator matrices, i.e., $(L_{i}^{n})_{ab}$ where
$a$ and $b$ take the values $0$ and $1$. The spatial direction is labeled by
the index $i$ taking the values $1$, $2$, and $3$. The transition operators
are also labeled by the indices $m$ or $n$ taking the values from $0$ to $3$.

Because the matrix $h$ is block diagonal (see eq. (15)), the problem of
reducing the master equation given in eq. (11) into Lindblad form becomes the
problem of finding a collapse operator for each of the blocks. Let us
consider, for instance, the block elements $h_{00}=h_{11}=0$ and
$h_{10}=h_{01}=1$. The collapse operator $C$ we seek must satisfy

$\sum_{n,m=0}^{1}h_{nm}\left(L^{n}\rho(t)L^{m\,\dagger}-\frac{1}{2}\left\\{L^{m\,\dagger}L^{n},\rho(t)\right\\}\right)=C\rho(t)C^{\dagger}-\frac{1}{2}\left\\{C^{\dagger}C,\rho(t)\right\\}.$
(89)

It therefore exists if the equation

$\begin{pmatrix}L^{0\,\dagger}&L^{1\,\dagger}\end{pmatrix}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\begin{pmatrix}L^{0}\\\
L^{1}\end{pmatrix}=L^{0\,\dagger}L^{1}+L^{1\,\dagger}L^{0}=C^{\dagger}C$ (90)

has a solution. An equivalent condition defines a second collapse operator
from $L^{2}$ and $L^{3}$, when considering the block elements
$h_{22}=h_{33}=0$ and $h_{23}=h_{32}=1$.

We investigate under which circumstances eq. (90) has a solution. First, we
examine the case when $L^{0}$ and $L^{1}$ are linearly dependent, i.e.,

$L^{1}=(a+ib)L^{0},$ (91)

with $a$ and $b$ real numbers. This is the case at LO in the $E/T$ expansion,
where both $L^{0}$ and $L^{1}$ are proportional to $r_{i}$. Inserting eq. (91)
into eq. (90), we obtain

$2aL^{0\,\dagger}L^{0}=C^{\dagger}C,$ (92)

the solution of which is

$C=\sqrt{2a}L^{0}.$ (93)

When rewriting this result in matrix form starting from the left-hand side of
eq. (90),

$\begin{pmatrix}L^{0\,\dagger}&(a-ib)L^{0\,\dagger}\end{pmatrix}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\begin{pmatrix}L^{0}\\\
(a+ib)L^{0}\end{pmatrix}=\begin{pmatrix}0&\sqrt{2a}L^{0\,\dagger}\end{pmatrix}\begin{pmatrix}-1&0\\\
0&1\end{pmatrix}\begin{pmatrix}0\\\ \sqrt{2a}L^{0}\end{pmatrix},$ (94)

one may rotate the vector $(L^{0},L^{1})$ in such a way that it becomes
orthogonal to the eigenspace of the matrix $h_{nm}$ with negative eigenvalue.

We now examine the case when $L^{0}$ and $L^{1}$ are neither linearly
dependent nor orthogonal, i.e.,

$L^{1}=(a+ib)L^{0}+\epsilon L^{1}_{(1)},$ (95)

where $\epsilon$ is a complex scalar and $L^{1}_{(1)}$ is nonzero and linearly
independent of $L^{0}$, i.e., there exists no scalar $c$ for which
$L^{0}=cL^{1}_{(1)}$. This is the case at NLO in the $\epsilon=E/T$ expansion,
where $L^{1}$ contains both a term proportional to $r_{i}$, and hence to
$L_{0}$, and a term proportional to $p_{i}$ (see eq. (26)). In this case, the
left-hand side of eq. (90) can be written as

$\displaystyle\begin{pmatrix}L^{0\,\dagger}&\left[(a-ib)L^{0\,\dagger}+\epsilon^{*}L^{1\,\dagger}_{(1)}\right]\end{pmatrix}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\begin{pmatrix}L^{0}\\\
\left[(a+ib)L^{0}+\epsilon^{*}L^{1}_{(1)}\right]\end{pmatrix}=$
$\displaystyle\hskip
56.9055pt\begin{pmatrix}\frac{\epsilon^{*}}{\sqrt{2a}}L^{1\,\dagger}_{(1)}&\left[\sqrt{2a}L^{0\,\dagger}+\frac{\epsilon^{*}}{\sqrt{2a}}L^{1\,\dagger}_{(1)}\right]\end{pmatrix}\begin{pmatrix}-1&0\\\
0&1\end{pmatrix}\begin{pmatrix}\frac{\epsilon}{\sqrt{2a}}L^{1}_{(1)}\\\
\left[\sqrt{2a}L^{0}+\frac{\epsilon}{\sqrt{2a}}L^{1}_{(1)}\right]\end{pmatrix}.$
(96)

This is not of Lindblad form as the negative eigenvalue of the matrix $h_{nm}$
contributes through the vector component $\epsilon L^{1}_{(1)}/\sqrt{2a}$.
However, in the context of an expansion in $\epsilon$, discounting this
component affects the master equation only at order $\epsilon^{2}$. Hence, if
we aim to obtain a master equation that is accurate to order $\epsilon$, it
holds that

$\begin{pmatrix}L^{0\,\dagger}&\left[(a-ib)L^{0\,\dagger}+\epsilon^{*}L^{1\,\dagger}_{(1)}\right]\end{pmatrix}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\begin{pmatrix}L^{0}\\\
\left[(a+ib)L^{0}+\epsilon^{*}L^{1}_{(1)}\right]\end{pmatrix}=C^{\dagger}C+O(\epsilon^{2}),$
(97)

with collapse operator

$C=\sqrt{2a}L^{0}+\frac{\epsilon}{\sqrt{2a}}L^{1}_{(1)}.$ (98)

Equivalently, one can solve eq. (90) for $L^{1}$ as given in eq. (95), i.e.,

$\displaystyle L^{0\,\dagger}\left[(a+ib)L^{0}+\epsilon
L^{1}_{(1)}\right]+\left[(a-ib)L^{0\,\dagger}+\epsilon^{*}L^{1\,\dagger}_{(1)}\right]L^{0}$
$\displaystyle=2aL^{0\,\dagger}L^{0}+\epsilon
L^{0\,\dagger}L^{1}_{(1)}+\epsilon^{*}L^{1\,\dagger}_{(1)}L^{0}$
$\displaystyle=C^{\dagger}C,$ (99)

whose solution is eq. (98) up to corrections of order $\epsilon^{2}$ to
$C^{\dagger}C$.

##  Appendix C Lindblad equation in the spherical basis

In this section, we examine the angular momentum structure of the Lindblad
equation. We define the density matrix projected onto states of definite
azimuthal quantum number $l$ and magnetic quantum number $m$ as

$\rho^{l^{\prime}m^{\prime};lm}=\langle
l^{\prime}m^{\prime}|\rho|lm\rangle=\int\text{d}\Omega\,Y_{l^{\prime}m^{\prime}*}(\theta,\phi)\,\rho\,Y_{lm}(\theta,\phi)\,.$
(100)

Since the density matrix is radially symmetric, i.e., it possesses no
preferred direction, it is diagonal in $l$ and $m$, i.e.,

$\rho^{l^{\prime}m^{\prime};lm}=\rho^{lm;lm}\delta_{ll^{\prime}}\delta_{mm^{\prime}}\,.$
(101)

Radial symmetry, furthermore, implies equal probability for all polarizations
within an orbital. All information can, therefore, be encoded in

$\rho^{l}=\sum_{m}\rho^{lm;lm},\quad\rho^{lm;lm}=\frac{1}{2l+1}\rho^{l}\,.$
(102)

We note that for our calculation (both at LO and NLO), the Lindblad equation
contains only vector operators. Given a Lindblad equation with a vector jump
operator ($\vec{C}=\\{C_{x},C_{y},C_{z}\\}$), one can first go to the
spherical basis given by the operators
$\\{C^{(1)}_{+1},C^{(1)}_{0},C^{(1)}_{-1}\\}$, defined as

$C^{(1)}_{\pm 1}=\mp\frac{1}{\sqrt{2}}\left(C_{x}\pm iC_{y}\right),\quad
C^{(1)}_{0}=C_{z}\,,$ (103)

where, anticipating the use of the Wigner–Eckart theorem, we explicitly label
the rank of the spherical tensor operator with the superscript $(1)$. An inner
product of the form $C^{\dagger}_{i}C_{i}$ is written the same in the
Cartesian and spherical bases, i.e.,

$\sum_{i}C_{i}^{\dagger}C_{i}=\sum_{q}C^{(1)\dagger}_{q}C^{(1)}_{q}\quad,\quad\sum_{i}C_{i}\rho\,C^{\dagger}_{i}=\sum_{q}C^{(1)}_{q}\rho\,C^{(1)\dagger}_{q}\,,$
(104)

where it should be noted that $C^{(1)}_{+1}=-(C^{(1)}_{-1})^{\dagger}$ and
$C^{(1)}_{-1}=-(C^{(1)}_{+1})^{\dagger}$.

We start by projecting the Lindblad equation onto the spherical basis. The
terms containing scalar operators (such as $\text{Re}[H^{\mathrm{eff}}_{s/o}]$
and $\text{Im}[H^{\mathrm{eff}}_{s/o}]$) are easy to project

$\displaystyle\sum_{m}\bra{lm}S\rho\ket{lm}$
$\displaystyle=\sum_{m,m^{\prime},l^{\prime}}\bra{lm}S\ket{l^{\prime}m^{\prime}}\bra{l^{\prime}m^{\prime}}\rho\ket{lm}$
$\displaystyle=\sum_{m}\bra{lm}S\ket{lm}\rho^{lm;lm}$
$\displaystyle=S\rho_{l}\,,$ (105)

where in the second line we have used
$\rho^{lm;l^{\prime}m^{\prime}}=\rho^{lm;lm}\delta_{ll^{\prime}}\delta_{mm^{\prime}}$
and in the last line we have used $\bra{lm}S\ket{lm}=S$, which is true for any
scalar operator $S$ that is a function of radial operators only. As a result,
we obtain for the scalar terms

$\displaystyle\sum_{m}\bra{lm}\frac{\mathrm{d}\rho(\vec{r},\vec{r}^{\,\prime},t)}{\mathrm{d}t}\ket{lm}$
$\displaystyle=\frac{\mathrm{d}\rho_{l}(r,r^{\prime},t)}{\mathrm{d}t},$
$\displaystyle\sum_{m}\bra{lm}H^{\mathrm{eff}}_{s/o},\rho(\vec{r},\vec{r}^{\prime},t)\ket{lm}$
$\displaystyle=H^{\mathrm{eff}}_{s/o}\rho_{l}(r,r^{\prime},t)\,.$ (106)

The only non-trivial calculation is that of the term $C_{i}\rho
C^{\dagger}_{i}$ as the operators $C_{i}$ are vector operators and give rise
to transitions between different states in color and angular momentum. The
evaluation of this term can be done using the Wigner–Eckart theorem. The
Wigner–Eckart theorem states that we can write the angular momentum matrix
element of a vector operator $C^{(1)}_{q}$ as

$\langle l^{\prime}m^{\prime}|C^{(1)}_{q}|lm\rangle=\langle
lm;1q|l^{\prime}m^{\prime}\rangle\langle l^{\prime}||C^{(1)}||l\rangle\,.$
(107)

The advantage of using the spherical basis is clear now. The angular matrix
element of the operators $C^{(1)}_{0,\pm 1}$ can be expressed as
Clebsch–Gordon coefficients $\langle lm;1q|l^{\prime}m^{\prime}\rangle$
multiplied by the reduced matrix element $\langle
l^{\prime}||C^{(1)}||l\rangle$.

###  C.1 Projecting the jump term $C_{i}\rho C^{\dagger}_{i}$ on the spherical
basis

For the calculation of the jump term, we begin by projecting onto the
eigenstate $\ket{lm}$

$\displaystyle\sum_{m,q}\bra{lm}C^{(1)}_{q}\rho C^{(1)\,\dagger}_{q}\ket{lm}$
$\displaystyle=\sum_{\begin{subarray}{c}q,m,l_{1},m_{1},\\\
l_{2},m_{2}\end{subarray}}\bra{lm}C^{(1)}_{q}\ket{l_{1}m_{1}}\bra{l_{1}m_{1}}\rho\ket{l_{2}m_{2}}\bra{l_{2}m_{2}}C^{(1)\,\dagger}_{q}\ket{lm}$
$\displaystyle=\sum_{q,m,l_{1},m_{1}}\bra{lm}C^{(1)}_{q}\ket{l_{1}m_{1}}\rho^{l_{1}m_{1};l_{1}m_{1}}\bra{l_{1}m_{1}}C^{(1)\,\dagger}_{q}\ket{lm}$
$\displaystyle=\sum_{q,m,l_{1},m_{1}}\bra{l_{1}m_{1};1q}lm\rangle\langle
l||C^{(1)}||l_{1}\rangle\rho^{l_{1}m_{1};l_{1}m_{1}}\braket{lm;1q}{l_{1}m_{1}}\langle
l_{1}||C^{(1)\,\dagger}||l\rangle.$ (108)

The Lindblad operators at NLO can be parametrized as $C_{i}=ar_{i}+bp_{i}$
where $a$ and $b$ are radially symmetric complex matrix operators. To proceed
further, we need to evaluate the reduced matrix elements $\langle
l^{\prime}||C^{(1)}||l\rangle$ for the two cases separately,
$C_{i}=r_{i},p_{i}$. We present the details for each of these calculations in
the next subsections.

###  C.2 Calculation of $\langle l^{\prime}||r^{(1)}||l\rangle$

To calculate $\langle l^{\prime}||r^{(1)}||l\rangle$, we make use of the
Wigner–Eckart relation. To evaluate the reduced matrix element $\langle
l^{\prime}||r^{(1)}||l\rangle$, we consider the case of $q=m=m^{\prime}=0$,
i.e., from eq. (107), we have

$\langle l^{\prime}||r^{(1)}||l\rangle=\frac{\langle
l^{\prime}0|z|l0\rangle}{\langle l0;10|l^{\prime}0\rangle},$ (109)

where $z=r\cos{\theta}$ can be expressed in terms of spherical harmonics using
$\cos\theta=\sqrt{(4\pi)/3}\,Y_{10}$. By using the following identities for
the spherical harmonic functions and Clebsch–Gordon coefficients,

$\int\text{d}\,\Omega\,Y_{l_{1}m_{1}}\,Y_{l_{2}m_{2}}\,Y^{*}_{l_{3}m_{3}}=\sqrt{\frac{(2l_{1}+1)(2l_{2}+1)}{4\pi(2l_{3}+1)}}\langle
l_{1}m_{1};\,l_{2}m_{2}|l_{3}m_{3}\rangle\langle
l_{1}0;\,l_{2}0|l_{3}0\rangle\,,$ (110)

and

$\langle
l0;\,10|l^{\prime}0\rangle=\sqrt{\frac{l+1}{2l+1}}\delta_{l^{\prime},l+1}-\sqrt{\frac{l}{2l+1}}\delta_{l^{\prime},l-1}\,,$
(111)

we obtain

$\langle l^{\prime}||r^{(1)}||l\rangle=\frac{\langle
l^{\prime}0|r\cos\theta|l0\rangle}{\langle
l0;10|l^{\prime}0\rangle}=\sqrt{\frac{l+1}{2l+3}}\delta_{l^{\prime},l+1}r-\sqrt{\frac{l}{2l-1}}\delta_{l^{\prime},l-1}r\,.$
(112)

###  C.3 Calculation of $\langle l^{\prime}||p^{(1)}||l\rangle$

To evaluate the reduced momentum operator $\langle
l^{\prime}||p^{(1)}||l\rangle$ we need to evaluate

$\langle l^{\prime}||p^{(1)}||l\rangle=\frac{\langle
l^{\prime}0|p_{z}|l0\rangle}{\langle l0;10|l^{\prime}0\rangle}\,.$ (113)

The momentum operator $p_{z}=-i\partial_{z}$ in spherical coordinates is given
by

$p_{z}=-i\left[\cos\theta\,\partial_{r}-\frac{\sin\theta}{r}\partial_{\theta}\right].$
(114)

The calculation of the term $\cos{\theta}\,\partial_{r}$ follows that of
$r\cos\theta$ and reduces to

$\frac{\langle l^{\prime}0|\cos\theta\partial_{r}|l0\rangle}{\langle
l0;10|l^{\prime}0\rangle}=\sqrt{\frac{l+1}{2l+3}}\delta_{l^{\prime},l+1}\,\partial_{r}-\sqrt{\frac{l}{2l-1}}\delta_{l^{\prime},l-1}\,\partial_{r}\,.$
(115)

Evaluating the second term using the relation

$Y^{*}_{l0}(\theta,\phi)=Y_{l0}(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}}P_{l}(\cos\theta),$
(116)

one finds

$\displaystyle\langle l^{\prime}0|\sin\theta\partial_{\theta}|l0\rangle$
$\displaystyle=\int\text{d}\Omega\,Y_{l0}^{*}(\theta,\phi)\sin\theta\,\partial_{\theta}\,Y_{l0}(\theta,\phi)$
$\displaystyle=\frac{1}{2}\sqrt{(2l^{\prime}+1)(2l+1)}\int_{-1}^{-1}\text{d}x\,P_{l^{\prime}}(x)(x^{2}-1)P^{\prime}_{l}(x)\,,$
(117)

where $x=\cos\theta$ and the prime indicates a derivative with respect to the
argument of the Legendre polynomial. We make use of a number of orthogonality
and recurrence properties of the Legendre polynomials

$\displaystyle(x^{2}-1)P^{\prime}_{l}$
$\displaystyle=l\left[xP_{l}(x)-P_{l-1}(x)\right]$
$\displaystyle=l\left[P_{1}(x)P_{l}(x)-P_{l-1}(x)\right],$ (118)
$\displaystyle\int_{-1}^{1}\text{d}x\,P_{l_{1}}(x)P_{l_{2}}(x)P_{l_{3}}(x)$
$\displaystyle=\frac{2}{2l_{3}+1}\langle l_{1}0;\,l_{2}0|l_{3}0\rangle^{2},$
(119) $\displaystyle\int_{-1}^{1}\text{d}x\,P_{l_{1}}(x)P_{l_{2}}(x)$
$\displaystyle=\frac{2}{2l_{1}+1}\delta_{l_{1}l_{2}}\,,$ (120)

to write

$\frac{1}{2}\sqrt{(2l^{\prime}+1)(2l+1)}\int_{-1}^{-1}\text{d}x\,P_{l^{\prime}}(x)(x^{2}-1)P^{\prime}_{l}(x)=l\sqrt{\frac{2l+1}{2l^{\prime}+1}}\left(\langle
l0;\,10|l^{\prime}0\rangle^{2}-\delta_{l^{\prime},l-1}\right).$ (121)

Using eq. (111), this becomes

$\langle
l^{\prime}0|\sin\theta\,\partial_{\theta}|l0\rangle=l\frac{l+1}{\sqrt{(2l+3)(2l+1)}}\delta_{l^{\prime},l+1}-l\frac{l+1}{\sqrt{(2l-1)(2l+1)}}\delta_{l^{\prime},l-1}\,,$
(122)

and, furthermore,

$\frac{\langle l^{\prime}0|\sin\theta\,\partial_{\theta}|l0\rangle}{\langle
l0;10|l^{\prime}0\rangle}=l\sqrt{\frac{l+1}{2l+3}}\delta_{l^{\prime},l+1}+(l+1)\sqrt{\frac{l}{2l-1}}\delta_{l^{\prime},l-1}\,.$
(123)

Combining eqs. (113), (114), (115), and (123), we obtain

$\displaystyle\langle l^{\prime}||p^{(1)}||l\rangle=$
$\displaystyle\frac{\displaystyle\langle
l^{\prime}0|\left[-i\left(\cos\theta\,\partial_{r}-\frac{\sin\theta}{r}\,\partial_{\theta}\right)\right]|l0\rangle}{\langle
l0;10|l^{\prime}0\rangle}$ $\displaystyle=$
$\displaystyle\hphantom{{}-{}}\sqrt{\frac{l+1}{2l+3}}\left[-i\left(\partial_{r}-\frac{l}{r}\right)\right]\delta_{l^{\prime},l+1}-\sqrt{\frac{l}{2l-1}}\left[-i\left(\partial_{r}+\frac{l+1}{r}\right)\right]\delta_{l^{\prime},l-1}\,$
(124)

We note that the reduced matrix elements $\langle l\pm 1||p^{(1)}||l\rangle$
are not Hermitian. To derive the Hermitian conjugate of the
$p^{\uparrow/\downarrow}$ reduced matrix elements, one can use the following
result for the Hermitian conjugate of the partial derivative operator
$-i\partial_{r}$ in the radial basis

$\langle\psi|-i\partial_{r}|\phi\rangle=\int\mathrm{d}r\
r^{2}\psi^{*}(r)[-i\partial_{r}\phi(r)]=\int\mathrm{d}r\
r^{2}\left[-i\left(\partial_{r}+\frac{2}{r}\right)\psi(r)\right]^{*}\phi(r)\,,$
(125)

which gives $(-i\partial_{r})^{\dagger}=-i\left(\partial_{r}+2/r\right)$.

We use the following identity for simplifying sums over $q$ and $m$ (which can
be shown for any general vector operator $C^{(1)}$)

$\displaystyle\sum_{m,q}\bra{lm}C^{(1)}_{q}\rho(\vec{r},\vec{r}^{\,\prime},t)(C^{(1)}_{q})^{\dagger}$
$\displaystyle\ket{lm}=\frac{2l-1}{2l+1}C^{\downarrow}\rho^{l-1}(r,r^{\prime},t)C^{\downarrow\,\dagger}+\frac{2l+3}{2l+1}C^{\uparrow}\rho^{l+1}(r,r^{\prime},t)C^{\uparrow\,\dagger}$
(126)

where $C^{\uparrow,\downarrow}=\langle l\pm 1||C^{(1)}||l\rangle$. Finally,
using eqs. (112), (C.3) and (126), we obtain

$\displaystyle\sum_{m,i}\bra{lm}(ar_{i}+bp_{i})\rho(\vec{r},\vec{r}^{\,\prime},t)(ar_{i}+bp_{i})^{\dagger}\ket{lm}=$
$\displaystyle\hskip 85.35826pt\frac{l+1}{2l+1}\left[ar-
bi\left(\partial_{r}-\frac{l}{r}\right)\right]\rho^{l-1}\left[a^{\dagger}r-b^{\dagger}i\left(\partial_{r}+\frac{l+2}{r}\right)\right]$
$\displaystyle\hskip 85.35826pt+\frac{l}{2l+1}\left[ar-
bi\left(\partial_{r}+\frac{l+1}{r}\right)\right]\rho^{l+1}\left[a^{\dagger}r-b^{\dagger}i\left(\partial_{r}-\frac{l-1}{r}\right)\right].$
(127)

In summary, we have shown that starting from the three-dimensional Lindblad
equation given in eq. (28) and using the rotational symmetry present in the
problem one can reduce it to a one-dimensional Lindblad equation with
operators acting only on the radial part of the wave function.

###  C.4 Projection into the reduced radial space

A further simplification of both the Lindblad operators and required radial
differential operators is possible when one projects them into the reduced
wave function basis ($u$-space). We have already noted that in eq. (125) the
derivative operator is not Hermitian. In numerical implementations this
creates a problem since it is not straightforward to implement such a non-
Hermitian derivative on a one-dimensional lattice. To fix this problem, as is
familiar in the construction of a one-dimensional effective Hamiltonian for
the hydrogen atom, we redefine the radial wave function $\psi(r)=u(r)/r$,
where $u(r)$ is the reduced wave function. We then project all operators into
the reduced wave function basis ($u$-space), i.e.,

$\int\mathrm{d}r\ r^{2}\psi^{*}(r)\hat{A}\phi(r)\equiv\int\mathrm{d}r\
u_{\psi}^{*}(r)\hat{A}_{u}u_{\phi}(r)\,.$ (128)

For example, the radial derivative operator $-i\partial_{r}$, which is not
Hermitian in the radial wave function basis, is Hermitian in the reduced wave
function basis ($u$-space)

$\displaystyle\int\mathrm{d}r\
r^{2}\psi^{*}(r)\left[-i\partial_{r}\phi(r)\right]=\int\mathrm{d}r\
r^{2}\left[-i\left(\partial_{r}+\frac{2}{r}\right)\psi(r)\right]^{*}\phi(r)\,,$
$\displaystyle\Rightarrow$ $\displaystyle\int\mathrm{d}r\
r^{2}\frac{u_{\psi}^{*}(r)}{r}\left[-i\partial_{r}\frac{u_{\phi}(r)}{r}\right]=\int\mathrm{d}r\
r^{2}\left[-i\left(\partial_{r}+\frac{2}{r}\right)\frac{u_{\psi}(r)}{r}\right]^{*}\frac{u_{\phi}(r)}{r}\,,$
$\displaystyle\Rightarrow$ $\displaystyle\int\mathrm{d}r\
u_{\psi}^{*}(r)\left[-i\left(\partial_{r}-\frac{1}{r}\right)u_{\phi}(r)\right]=\int\mathrm{d}r\
\left[-i\left(\partial_{r}+\frac{1}{r}\right)u_{\psi}(r)\right]^{*}u_{\phi}(r)\,,$
(129)

which gives
$[-i(\partial_{r}+{1}/{r})]^{\dagger}=-i(\partial_{r}-{1}/{r})\Rightarrow(-i\partial_{r})^{\dagger}=-i\partial_{r}$
in $u$-space. We denote the spatial derivative in $u$-space as $p_{r}$ to not
cause any confusion. The following operators have different forms in the
reduced wave function basis ($u$-space) compared to their counterparts in the
radial basis ($\psi$-space)

$\displaystyle\frac{\vec{p}^{\,2}}{M}$
$\displaystyle\rightarrow-\frac{1}{M}\partial_{r}^{2}+\frac{l(l+1)}{Mr^{2}}=\frac{\mathcal{\overline{D}}^{2}}{M}\,,$
(130) $\displaystyle\sum_{i}\\{p_{i},r_{i}\\}$ $\displaystyle\rightarrow
p_{r}r+rp_{r}\,,$ (131)
$\displaystyle-i\sqrt{\frac{l+1}{2l+1}}\left(\partial_{r}-\frac{l}{r}\right)$
$\displaystyle\rightarrow-i\sqrt{\frac{l+1}{2l+1}}\left(\partial_{r}-\frac{l+1}{r}\right),$
(132)
$\displaystyle-i\sqrt{\frac{l}{2l+1}}\left(\partial_{r}+\frac{l+1}{r}\right)$
$\displaystyle\rightarrow-i\sqrt{\frac{l}{2l+1}}\left(\partial_{r}+\frac{l}{r}\right).$
(133)

##  Appendix D Transport coefficient $\kappa$

As has been pointed out in refs. Eller:2019spw ; Binder:2021otw , the
transport coefficient $\kappa$ that enters into the quarkonium evolution
equation and the heavy quark momentum diffusion coefficient computed on the
lattice differ in the time ordering of the fields in the correlator. They both
vanish in vacuum at all orders in perturbation theory and dimensional
regularization, but they could differ due to finite temperature effects.
Investigations are under way, but no quantitative conclusion has been reached
yet, neither in weakly coupled thermal field theory, nor in lattice thermal
QCD. In lattice thermal QCD, the coefficient $\kappa$, as defined in eq. (76),
has not been calculated directly. Indirect determinations based on the pNRQCD
definition of the quarkonium thermal width, which is proportional to $\kappa$,
and unquenched lattice measurements of the quarkonium thermal width provide a
range that fully overlaps with the direct lattice determinations of the heavy
quark momentum diffusion coefficient Brambilla:2019tpt . In the body of the
paper, in the absence of a direct lattice determination of the correlator
defining the quarkonium momentum diffusion coefficient, we follow the approach
of refs. Brambilla:2020qwo ; Brambilla:2021wkt  and parameterize the
temperature dependence of $\kappa$ based on the lattice determination of the
heavy quark momentum diffusion coefficient done in ref. Brambilla:2020siz .

Figure 6:  Predictions for $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$
$R_{AA}$ as a function of $N_{\rm part}$ (left) and $p_{T}$ (right) with
$\tau_{\text{med}}=0.6$ fm, NLO evolution, and $T_{f}=190$ MeV. In both
panels, the lower, central, and upper curves correspond to
$\hat{\kappa}=\\{0.24,\kappa_{C}(T),4.2\\}$ with $\hat{\gamma}=-2.6$.

Using only indirect lattice determinations of $\kappa$ would not significantly
change the central value of our results but somewhat enlarges the
uncertainties due to the variation of $\kappa$. In fig. 6, we estimate these
errors by running simulations with $\hat{\kappa}=0.24$ and $4.2$ (with
$\hat{\gamma}=2.6$) as taken from eq. (36) of ref. Brambilla:2019tpt . We
observe larger uncertainty compared to the $\hat{\kappa}$ variation presented
in figs. 4 and 5, but, even with this larger $\hat{\kappa}$ variation, we find
smaller uncertainty compared to the $\hat{\gamma}$ variation.

##  Appendix E Dependence on the medium initialization time

In this appendix, we consider the dependence of our results on the assumed
medium initialization time $\tau_{\rm med}$. In the body of the paper, we
assumed $\tau_{\rm med}$ = 0.6 fm. Here we consider the effect of an early
medium initialization time by choosing instead $\tau_{\rm med}$ = 0.25 fm,
which is the earliest proper time available in the hydrodynamic simulations.
In figs. 7 and 8, we present our NLO predictions for $R_{AA}$ as a function of
$N_{\rm part}$ and $p_{T}$, respectively, including the effect of excited
state feed down. In these figures, the central solid line corresponds to
$\hat{\gamma}=-1.75$ as opposed to the value of $\hat{\gamma}=-2.6$ used in
figs. 4 and 5. Due to the shift in the central value of $\hat{\gamma}$, both
choices for $\tau_{\rm med}$ provide a good description of the data. From the
comparison with the case of $\tau_{\rm med}$= 0.6 fm presented in the main
body of the paper, we see that the uncertainty of our inferred central value
of $\hat{\gamma}$ is $\Delta\hat{\gamma}\sim 1$. We note that starting the
Lindblad evolution at 0.25 fm introduces additional uncertainties related to
the fact that at early times in the QGP evolution the system is far from
isotropic thermal evolution Almaalol:2020rnu ; Strickland:2013uga . In the
future, it may be possible to use methods such as those presented in ref.
Dong:2021gnb  to include the effect of early-time momentum-space anisotropies
in a systematic fashion.

Figure 7: $R_{AA}$ predictions for the $\Upsilon(1S)$, $\Upsilon(2S)$, and
$\Upsilon(3S)$ as a function of $N_{\rm part}$ with $\tau_{\text{med}}=0.25$
fm. The left panel shows variation of
$\hat{\kappa}\in\\{\kappa_{L}(T),\kappa_{C}(T),\kappa_{U}(T)\\}$ and the right
panel shows variation of $\hat{\gamma}$ in the range $-3.5\leq\hat{\gamma}\leq
0$. In both panels, the solid line corresponds to
$\hat{\kappa}=\hat{\kappa}_{C}(T)$ and $\hat{\gamma}=-1.75$. The experimental
measurements shown are from the ALICE Acharya:2020kls , ATLAS ATLAS5TeV , and
CMS Sirunyan:2018nsz ; CMSupsilonQM2022  collaborations.

Figure 8:  $R_{AA}$ for the $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$
as a function of $p_{T}$. The value of $\tau_{\rm med}$, the meaning of the
bands, and experimental data sources are the same as in fig. 7.

## References

  * (1) T. Matsui and H. Satz, _$J/\psi$ suppression by quark-gluon plasma formation_, _Phys. Lett._ B178 (1986) 416.
  * (2) ALICE collaboration, S. Acharya et al., _$\Upsilon$ production and nuclear modification at forward rapidity in Pb-Pb collisions at $\mathbf{\sqrt{\textit{s}_{\textbf{NN}}}=5.02}$ TeV_, 2011.05758.
  * (3) Songkyo Lee (ATLAS Collaboration), “Quarkonium production in Pb+Pb collisions with ATLAS.” Quark Matter 2020 https://indico.cern.ch/event/792436/contributions/3535775/, 2017.
  * (4) CMS collaboration, A. M. Sirunyan et al., _Measurement of nuclear modification factors of $\Upsilon$(1S), $\Upsilon$(2S), and $\Upsilon$(3S) mesons in PbPb collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}=$ 5.02 TeV_, _Phys. Lett. B_ 790 (2019) 270 [1805.09215].
  * (5) CMS collaboration, A. M. Sirunyan et al., _Suppression of Excited $\Upsilon$ States Relative to the Ground State in Pb-Pb Collisions at $\sqrt{s_{\mathrm{NN}}}$=5.02 TeV_, _Phys. Rev. Lett._ 120 (2018) 142301 [1706.05984].
  * (6) CMS collaboration, A. M. Sirunyan et al., _Measurement of the azimuthal anisotropy of Image 1 and Image 2 mesons in PbPb collisions at sNN=5.02TeV_ , _Phys. Lett. B_ 819 (2021) 136385 [2006.07707].
  * (7) ALICE collaboration, S. Acharya et al., _Measurement of $\Upsilon(1{\rm S})$ elliptic flow at forward rapidity in Pb-Pb collisions at $\sqrt{s_{\rm{NN}}}=5.02$ TeV_, _Phys. Rev. Lett._ 123 (2019) 192301 [1907.03169].
  * (8) STAR collaboration, L. Adamczyk et al., _Suppression of $\Upsilon$ production in d+Au and Au+Au collisions at $\sqrt{s_{NN}}$=200 GeV_, _Phys. Lett. B_ 735 (2014) 127 [1312.3675].
  * (9) PHENIX collaboration, A. Adare et al., _Measurement of $\Upsilon(1S+2S+3S)$ production in $p+p$ and Au$+$Au collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV_, _Phys. Rev. C_ 91 (2015) 024913 [1404.2246].
  * (10) STAR collaboration, L. Adamczyk et al., _$\Upsilon$ production in U + U collisions at $\sqrt{{s}_{NN}}=$ 193 GeV measured with the STAR experiment_, _Phys. Rev. C_ 94 (2016) 064904 [1608.06487].
  * (11) Soohwan Lee (CMS Collaboration), “Observation of the $\Upsilon(3S)$ meson and sequential suppression of states in PbPb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV.” Quark Matter 2022, Krakow, Poland, https://indico.cern.ch/event/895086/contributions/4716202/, 2022.
  * (12) M. Laine, O. Philipsen, P. Romatschke and M. Tassler, _Real-time static potential in hot QCD_ , _JHEP_ 03 (2007) 054 [hep-ph/0611300].
  * (13) N. Brambilla, J. Ghiglieri, A. Vairo and P. Petreczky, _Static quark-antiquark pairs at finite temperature_ , _Phys. Rev._ D78 (2008) 014017 [0804.0993].
  * (14) A. Beraudo, J.-P. Blaizot and C. Ratti, _Real and imaginary-time Q anti-Q correlators in a thermal medium_ , _Nucl. Phys. A_ 806 (2008) 312 [0712.4394].
  * (15) M. A. Escobedo and J. Soto, _Non-relativistic bound states at finite temperature (I): The Hydrogen atom_ , _Phys. Rev. A_ 78 (2008) 032520 [0804.0691].
  * (16) N. Brambilla, M. A. Escobedo, J. Ghiglieri, J. Soto and A. Vairo, _Heavy Quarkonium in a weakly-coupled quark-gluon plasma below the melting temperature_ , _JHEP_ 09 (2010) 038 [1007.4156].
  * (17) N. Brambilla, M. A. Escobedo, J. Ghiglieri and A. Vairo, _Thermal width and gluo-dissociation of quarkonium in pNRQCD_ , _JHEP_ 12 (2011) 116 [1109.5826].
  * (18) N. Brambilla, M. A. Escobedo, J. Ghiglieri and A. Vairo, _Thermal width and quarkonium dissociation by inelastic parton scattering_ , _JHEP_ 05 (2013) 130 [1303.6097].
  * (19) A. Dumitru, Y. Guo and M. Strickland, _The Imaginary part of the static gluon propagator in an anisotropic (viscous) QCD plasma_ , _Phys.Rev._ D79 (2009) 114003 [0903.4703].
  * (20) M. Berwein, N. Brambilla, P. Petreczky and A. Vairo, _Polyakov loop correlator in perturbation theory_ , _Phys. Rev. D_ 96 (2017) 014025 [1704.07266].
  * (21) TUMQCD collaboration, A. Bazavov, N. Brambilla, P. Petreczky, A. Vairo and J. H. Weber, _Color screening in (2+1)-flavor QCD_ , _Phys. Rev. D_ 98 (2018) 054511 [1804.10600].
  * (22) A. Rothkopf, T. Hatsuda and S. Sasaki, _Complex Heavy-Quark Potential at Finite Temperature from Lattice QCD_ , _Phys. Rev. Lett._ 108 (2012) 162001 [1108.1579].
  * (23) A. Rothkopf, _Heavy quarkonium in extreme conditions_ , _Phys. Rept._ 858 (2020) 1 [1912.02253].
  * (24) P. Petreczky, A. Rothkopf and J. Weber, _Realistic in-medium heavy-quark potential from high statistics lattice QCD simulations_ , _Nucl. Phys. A_ 982 (2019) 735 [1810.02230].
  * (25) D. Bala and S. Datta, _Nonperturbative potential for the study of quarkonia in QGP_ , _Phys. Rev. D_ 101 (2020) 034507 [1909.10548].
  * (26) M. Laine, O. Philipsen and M. Tassler, _Thermal imaginary part of a real-time static potential from classical lattice gauge theory simulations_ , _JHEP_ 09 (2007) 066 [0707.2458].
  * (27) A. Lehmann and A. Rothkopf, _Proper static potential in classical lattice gauge theory at finite T_ , _JHEP_ 07 (2021) 067 [2012.10089].
  * (28) K. Boguslavski, B. S. Kasmaei and M. Strickland, _The imaginary part of the heavy-quark potential from real-time Yang-Mills dynamics_ , _JHEP_ 21 (2020) 083 [2102.12587].
  * (29) M. Strickland, _Thermal $\Upsilon(1s)$ and $\chi_{b1}$ suppression in $\sqrt{s_{NN}}=2.76$ TeV Pb-Pb collisions at the LHC_, _Phys.Rev.Lett._ 107 (2011) 132301 [1106.2571].
  * (30) M. Strickland and D. Bazow, _Thermal Bottomonium Suppression at RHIC and LHC_ , _Nucl. Phys._ A879 (2012) 25 [1112.2761].
  * (31) B. Krouppa, R. Ryblewski and M. Strickland, _Bottomonia suppression in 2.76 TeV Pb-Pb collisions_ , _Phys. Rev._ C92 (2015) 061901 [1507.03951].
  * (32) A. Islam and M. Strickland, _Bottomonium suppression and elliptic flow using Heavy Quarkonium Quantum Dynamics_ , _JHEP_ 21 (2020) 235 [2010.05457].
  * (33) A. Islam and M. Strickland, _Bottomonium suppression and elliptic flow from real-time quantum evolution_ , _Phys. Lett. B_ 811 (2020) 135949 [2007.10211].
  * (34) P. Braun-Munzinger and J. Stachel, _(Non)thermal aspects of charmonium production and a new look at J / psi suppression_ , _Phys. Lett. B_ 490 (2000) 196 [nucl-th/0007059].
  * (35) R. L. Thews, M. Schroedter and J. Rafelski, _Enhanced $J/\psi$ production in deconfined quark matter_, _Phys. Rev._ C63 (2001) 054905 [hep-ph/0007323].
  * (36) H. Breuer and F. Petruccione, _The theory of open quantum systems_. Oxford University Press, Oxford, 2002.
  * (37) Y. Akamatsu and A. Rothkopf, _Stochastic potential and quantum decoherence of heavy quarkonium in the quark-gluon plasma_ , _Phys. Rev._ D85 (2012) 105011 [1110.1203].
  * (38) Y. Akamatsu, _Quarkonium in quark-gluon plasma: open quantum system approaches re-examined_ , 2009.10559.
  * (39) Y. Akamatsu, _Heavy quark master equations in the Lindblad form at high temperatures_ , _Phys. Rev._ D91 (2015) 056002 [1403.5783].
  * (40) T. Miura, Y. Akamatsu, M. Asakawa and A. Rothkopf, _Quantum Brownian motion of a heavy quark pair in the quark-gluon plasma_ , _Phys. Rev. D_ 101 (2020) 034011 [1908.06293].
  * (41) N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, _Quarkonium suppression in heavy-ion collisions: an open quantum system approach_ , _Phys. Rev._ D96 (2017) 034021 [1612.07248].
  * (42) N. Brambilla, M. A. Escobedo, J. Soto and A. Vairo, _Heavy quarkonium suppression in a fireball_ , _Phys. Rev._ D97 (2018) 074009 [1711.04515].
  * (43) N. Brambilla, M. A. Escobedo, A. Vairo and P. Vander Griend, _Transport coefficients from in medium quarkonium dynamics_ , _Phys. Rev. D_ 100 (2019) 054025 [1903.08063].
  * (44) R. Sharma and A. Tiwari, _Quantum evolution of quarkonia with correlated and uncorrelated noise_ , _Phys. Rev. D_ 101 (2020) 074004 [1912.07036].
  * (45) J.-P. Blaizot, D. De Boni, P. Faccioli and G. Garberoglio, _Heavy quark bound states in a quark–gluon plasma: Dissociation and recombination_ , _Nucl. Phys. A_ 946 (2016) 49 [1503.03857].
  * (46) J.-P. Blaizot and M. A. Escobedo, _Quantum and classical dynamics of heavy quarks in a quark-gluon plasma_ , _JHEP_ 06 (2018) 034 [1711.10812].
  * (47) J.-P. Blaizot and M. A. Escobedo, _Approach to equilibrium of a quarkonium in a quark-gluon plasma_ , _Phys. Rev._ D98 (2018) 074007 [1803.07996].
  * (48) J.-P. Blaizot and M. A. Escobedo, _Phenomenological study of quarkonium suppression and the impact of the energy gap between singlets and octets_ , _Phys. Rev. D_ 104 (2021) 054034 [2106.15371].
  * (49) X. Yao and T. Mehen, _Quarkonium in-medium transport equation derived from first principles_ , _Phys. Rev. D_ 99 (2019) 096028 [1811.07027].
  * (50) X. Yao, W. Ke, Y. Xu, S. A. Bass and B. Müller, _Coupled Boltzmann Transport Equations of Heavy Quarks and Quarkonia in Quark-Gluon Plasma_ , _JHEP_ 21 (2020) 046 [2004.06746].
  * (51) X. Yao and T. Mehen, _Quarkonium Semiclassical Transport in Quark-Gluon Plasma: Factorization and Quantum Correction_ , _JHEP_ 21 (2020) 062 [2009.02408].
  * (52) R. Katz and P. B. Gossiaux, _The Schrödinger–Langevin equation with and without thermal fluctuations_ , _Annals Phys._ 368 (2016) 267 [1504.08087].
  * (53) W. Caswell and G. Lepage, _Effective Lagrangians for Bound State Problems in QED, QCD, and Other Field Theories_ , _Phys. Lett. B_ 167 (1986) 437.
  * (54) G. T. Bodwin, E. Braaten and G. Lepage, _Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium_ , _Phys. Rev. D_ 51 (1995) 1125 [hep-ph/9407339].
  * (55) A. Pineda and J. Soto, _Effective field theory for ultrasoft momenta in NRQCD and NRQED_ , _Nucl. Phys. B Proc. Suppl._ 64 (1998) 428 [hep-ph/9707481].
  * (56) N. Brambilla, A. Pineda, J. Soto and A. Vairo, _Potential NRQCD: An Effective theory for heavy quarkonium_ , _Nucl. Phys. B_ 566 (2000) 275 [hep-ph/9907240].
  * (57) N. Brambilla, A. Pineda, J. Soto and A. Vairo, _Effective field theories for heavy quarkonium_ , _Rev. Mod. Phys._ 77 (2005) 1423 [hep-ph/0410047].
  * (58) V. Gorini, A. Kossakowski and E. Sudarshan, _Completely positive dynamical semigroups of N level systems_ , _J. Math. Phys._ 17 (1976) 821.
  * (59) G. Lindblad, _On the Generators of Quantum Dynamical Semigroups_ , _Commun. Math. Phys._ 48 (1976) 119.
  * (60) N. Brambilla, M. A. Escobedo, M. Strickland, A. Vairo, P. Vander Griend and J. H. Weber, _Bottomonium suppression in an open quantum system using the quantum trajectories method_ , _JHEP_ 05 (2021) 136 [2012.01240].
  * (61) N. Brambilla, M. A. Escobedo, M. Strickland, A. Vairo, P. Vander Griend and J. H. Weber, _Bottomonium production in heavy-ion collisions using quantum trajectories: Differential observables and momentum anisotropy_ , _Phys. Rev. D_ 104 (2021) 094049 [2107.06222].
  * (62) H. B. Omar, M. A. Escobedo, A. Islam, M. Strickland, S. Thapa, P. Vander Griend et al., _QTRAJ 1.0: A Lindblad equation solver for heavy-quarkonium dynamics_ , _Comput. Phys. Commun._ 273 (2022) 108266 [2107.06147].
  * (63) J. Dalibard, Y. Castin and K. Molmer, _Wave-function approach to dissipative processes in quantum optics_ , _Phys. Rev. Lett._ 68 (1992) 580.
  * (64) Particle Data Group collaboration, P. A. Zyla et al., _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (65) D. De Boni, _Fate of in-medium heavy quarks via a Lindblad equation_ , _JHEP_ 08 (2017) 064 [1705.03567].
  * (66) A. O. Caldeira and A. J. Leggett, _Path integral approach to quantum Brownian motion_ , _Physica A_ 121 (1983) 587.
  * (67) L. Diósi, _On high-temperature markovian equation for quantum brownian motion_ , _EPL (Europhysics Letters)_ 22 (1993) 1.
  * (68) L. Diósi, _Calderia-leggett master equation and medium temperatures_ , _Physica A: Statistical Mechanics and its Applications_ 199 (1993) 517.
  * (69) S. Gao, _Dissipative quantum dynamics with a lindblad functional_ , _Physical review letters_ 79 (1997) 3101.
  * (70) B. Vacchini, _Completely positive quantum dissipation_ , _Physical review letters_ 84 (2000) 1374.
  * (71) Armadillo contributors, “Armadillo: C++ library for linear algebra and scientific computing.” http://arma.sourceforge.net/download.html, 2022\.
  * (72) M. Alqahtani and M. Strickland, _Bulk observables at 5.02 TeV using quasiparticle anisotropic hydrodynamics_ , 2008.07657.
  * (73) H. Alalawi, M. Alqahtani and M. Strickland, _Resummed Relativistic Dissipative Hydrodynamics_ , _Symmetry_ 14 (2022) 329 [2112.14597].
  * (74) Strickland et al, “QTraj 1.0 - Public Repository.” https://bitbucket.org/michael_strickland/qtraj-fftw, 2021.
  * (75) N. Brambilla, V. Leino, P. Petreczky and A. Vairo, _Lattice QCD constraints on the heavy quark diffusion coefficient_ , _Phys. Rev. D_ 102 (2020) 074503 [2007.10078].
  * (76) J. Casalderrey-Solana and D. Teaney, _Heavy quark diffusion in strongly coupled N=4 Yang-Mills_ , _Phys. Rev. D_ 74 (2006) 085012 [hep-ph/0605199].
  * (77) S. Caron-Huot and G. D. Moore, _Heavy quark diffusion in perturbative QCD at next-to- leading order_ , _Phys. Rev. Lett._ 100 (2008) 052301 [0708.4232].
  * (78) S. Kim, P. Petreczky and A. Rothkopf, _Quarkonium in-medium properties from realistic lattice NRQCD_ , _JHEP_ 11 (2018) 088 [1808.08781].
  * (79) G. Aarts, C. Allton, S. Kim, M. P. Lombardo, M. B. Oktay, S. M. Ryan et al., _What happens to the $\Upsilon$ and $\eta_{b}$ in the quark-gluon plasma? Bottomonium spectral functions from lattice QCD_, _JHEP_ 11 (2011) 103 [1109.4496].
  * (80) R. Larsen, S. Meinel, S. Mukherjee and P. Petreczky, _Thermal broadening of bottomonia: Lattice nonrelativistic QCD with extended operators_ , _Phys. Rev. D_ 100 (2019) 074506 [1908.08437].
  * (81) S. Shi, K. Zhou, J. Zhao, S. Mukherjee and P. Zhuang, _Heavy Quark Potential in QGP: DNN meets LQCD_ , 2105.07862.
  * (82) LHCb collaboration, R. Aaij et al., _Study of $\chi_{{\mathrm{b}}}$ meson production in $\mathrm{p}$ $\mathrm{p}$ collisions at $\sqrt{s}=7$ and $8{\mathrm{\,TeV}}$ and observation of the decay $\chi_{{\mathrm{b}}}\mathrm{(3P)}\rightarrow\Upsilon\mathrm{(3S)}{\gamma}$_, _Eur. Phys. J._ C74 (2014) 3092 [1407.7734].
  * (83) Y. Akamatsu, M. Asakawa, S. Kajimoto and A. Rothkopf, _Quantum dissipation of a heavy quark from a nonlinear stochastic Schrödinger equation_ , _JHEP_ 07 (2018) 029 [1805.00167].
  * (84) Y. Akamatsu, _Real-time quantum dynamics of heavy quark systems at high temperature_ , _Phys. Rev. D_ 87 (2013) 045016 [1209.5068].
  * (85) A. M. Eller, J. Ghiglieri and G. D. Moore, _Thermal Heavy Quark Self-Energy from Euclidean Correlators_ , _Phys. Rev. D_ 99 (2019) 094042 [1903.08064].
  * (86) T. Binder, K. Mukaida, B. Scheihing-Hitschfeld and X. Yao, _Non-Abelian electric field correlator at NLO for dark matter relic abundance and quarkonium transport_ , _JHEP_ 01 (2022) 137 [2107.03945].
  * (87) D. Almaalol, A. Kurkela and M. Strickland, _Non-equilibrium attractor in high-temperature QCD plasmas_ , _Phys. Rev. Lett._ 125 (2020) 122302 [2004.05195].
  * (88) M. Strickland, _Thermalization and isotropization in heavy-ion collisions_ , 1312.2285.
  * (89) L. Dong, Y. Guo, A. Islam and M. Strickland, _Effective Debye screening mass in an anisotropic quark gluon plasma_ , _Phys. Rev. D_ 104 (2021) 096017 [2109.01284].