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

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\"odinger 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.


Introduction
In the pioneering work of Matsui and Satz [1], 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. [2][3][4][5][6][7][8][9][10][11]. In ref. [1], 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 ∼ m −1 D , in addition to Debye screening of the real part of the potential, a large imaginary contribution to the potential also arises [12]. Perturbative calculations based on nonrelativistic effective theories of QCD at finite temperature [13][14][15][16][17][18] and perturbative QCD resummation [19] 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 [20,21]. The imaginary part of the potential has its origin in two effects: Landau damping (related to parton dissociation [18]) and singlet to octet transitions (an effect particular to QCD [13] related to the gluodissociation process [17]). Nonperturbative lattice QCD and classical real-time lattice measurements of the imaginary part of the potential were recently performed in refs. [22][23][24][25] and [26][27][28], 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. [29][30][31][32][33].
In addition to the thermal width and screening effects, an accurate description of inmedium quarkonium evolution must account for the recombination process [34,35] 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 [36] 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 [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52]. 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 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 ∼ M v, and the binding energy E ∼ M v 2 . An additional relevant scale is Λ QCD at which nonperturbative effects become dominant. Nonperturbative effects and the breaking of the weak coupling expansion when α s /v ∼ 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) [53,54] 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) [55][56][57], an EFT equivalent to QCD for energies much smaller than M v. 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 M v 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 [13,[15][16][17][18]. 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 [41][42][43] analyzed the regime M v 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. [41,42] 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 E the master equation could be cast in Lindblad form [58,59] at leading order in the E/T expansion. In a series of works [41,42,60,61], a number of phenomenological predictions were made for bottomonium suppression in the temperature regime M v T E; physically, this amounts to using bottomonium states of small radius as probes of a strongly coupled QGP. In the most recent works [60,61], the QTraj code of ref. [62], which implements the Monte Carlo quantum trajectories algorithm [63], 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 α s we take from the Particle Data Group [64], for the bottom mass, we work in the 1S scheme where m b = m Υ(1S) /2 with the Υ(1S) mass also from the Particle Data Group, while the nonperturbative transport coefficients κ and γ 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 [60,61], the requirement that the inequality T 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 [40,47,65] 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. [40,47,65] 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 M v T ∼ gT 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 inmedium 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 Υ(1S), Υ(2S), and Υ(3S) as a function of number of participating nucleons N 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 κ and γ. We conclude and provide an outlook on future prospects in sec. 5.

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. [41,42,60,61].

The pNRQCD master equation
The evolution of the reduced density matrix of heavy quarkonium in the regime M v T was derived in refs. [41,42]. If one does not make any assumption on the relation between T and E, then the evolution equations that one obtains are 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, 1 ) the Casimir of the fundamental representation, N c = 3 the number of colors, and α 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. [42,43].
These equations can be written as the master equation (2.13) 1 These assumptions account for the difference between eq. (2.8) and eq. (D.2) of ref. [42].
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. (2.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. (2.11) is Markovian and trace preserving but does not preserve complete positivity of the density matrix, it follows that eq. (2.11) allows for negative probabilities, a feature similar to the widely used Caldeira-Leggett master equation [66]. We note that complete positivity is necessary for use of the quantum trajectories method.

Lindblad equation to leading-order in E/T
Although, as explained above, the evolution specified in eq. (2.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 E. In eq. (2.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 uv i thus scale like E/T . In the limit T E, the exponentials can be set to 1, and A uv i can be written in terms of the transport coefficients κ and γ where the ellipsis indicates terms of order E/T and higher; κ is the heavy quarkonium momentum diffusion coefficient, and γ is its dispersive counterpart 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 with Hamiltonian where and collapse operators This equation was first derived and solved in refs. [41,42] and further studied in refs. [60,61] using the QTraj code of ref. [62].

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 E limit it takes the form of a Lindblad equation. In this subsection, we consider contributions to A uv i linear in E/T and investigate the conditions under which a Lindblad equation can still be derived. Making use of the relation where ∆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 uv i , the operators L 1 i and L 3 i in the master equation take the form while L 0 i and L 2 i keep the form given in eqs. (2.13) and (2.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 with Hamiltonian and collapse operators 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 [67][68][69][70]. We can understand the anti-commutator term {r i , p i } multiplying κ in Im (Σ s,o ) as a drag force. This identification is clear if one uses the procedure described in ref. [46] to obtain a Langevin equation from eq. (2.28).

Lindblad equation to NLO in E/T in the spherical basis
The six collapse operators in Eqs. (2.31) and (2.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 E. 2 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 As the plasma is isotropic, i.e., there is no preferred direction in space, only diagonal elements l = l and m = m are nonzero; furthermore, all information can be encoded in 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 × 2 color space, the infinite dimensional angular momentum space, and the radial wave function. The density matrix is diagonal and takes the form with the Hamiltonian being an analogous diagonal matrix, the elements of which are given by The collapse operators take the form where the 2 × 2 matrices are responsible for color transitions and the infinite dimensional matrices O ± l ,l = δ l ,l±1 for angular momentum transitions. The operators C ↑,↓ u→v act on ρ l u contributing to its partial width to ρ l±1 v ; they are given explicitly by 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. [42] truncating the angular momentum at l = 1. Since then, works based on the QTraj code [62] have, however, made it possible to account for states of arbitrarily high angular momentum [60,61].
The decay width of the state ρ l u (r) is given by the trace of the anticommutator term of the Lindblad equation where on the left Γ[ρ l u (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 where the traces are only over angular momentum and color, respectively, and |l ± 1 l ± 1| and |v v| are states of angular momentum l ± 1 and color v, respectively. Computing the operators Γ ↑,↓ and Γ u→v explicitly, we find The partial width of the state ρ l u (r) to a state of angular momentum l ± 1 and color v is given by where on the left Γ ↑,↓ u→v [ρ l u (r)] represents a functional returning the partial width and the trace is over the radial coordinate r.

The quantum trajectories algorithm at NLO
The Lindblad equation given in eq. (2.28) with the collapse operators of eqs. (2.31) and (2.32) describes the evolution of the heavy quarkonium reduced density matrix at next-toleading-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. [60][61][62]. 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 [62]. 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) ≡ rR(t, r), where the three-dimensional wave function for fixed l and m is given by ψ(t, r) = R(t, r)Y lm (θ, φ). 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.

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 → l + 1) or down (l → l − 1), thus totalling six. When expressed as operators acting on the reduced wave function, the six jump operators are 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.

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 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. (3.10), the up and down transitions are related by Γ ↓ s→o /Γ ↑ s→o = l/(l + 1).

Effective Hamiltonian in the reduced spherical representation
The effective Hamiltonian for singlet and octet evolution is defined by H eff . When expressed as operators acting on the reduced wave function, the singlet effective Hamiltonian H eff s is given by

Description of the algorithm
The NLO quantum trajectories algorithm proceeds as follows: 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 ∈ {0, 1}, with 0 being a color singlet and 1 being a color octet state, and the initial orbital angular momentum quantum number l ≥ 0.
3. Generate a random number p 1 uniformly distributed in [0, 1].  .15), 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.

Using
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. If the color state is singlet (c = 0): (a) Set the color state to octet, i.e. c → 1.
7. Based on the choices made above, apply the appropriate collapse/jump operator to the wavefunction u 8. Repeat starting from step 2.

Evolution between jumps
In step 4 above, we must evolve the wave function forward in time using the complex effective Hamiltonian H 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. [62]. 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 where the Hamiltonian has implicit time dependence. To proceed, one expresses this as 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 ≡ 1 − i∆tH eff s,o /2 ψ(t). The final step in the Crank-Nicolson algorithm is to solve for ψ j (t + ∆t) where i, j are spatial coordinates and

. 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 [71]. The Armadillo package is also used to handle all matrix multiplications required to compute the action of the effective Hamiltonian, width, and jump operators.

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 [72,73]. For this purpose, smooth optical Glauber initial conditions were used, and the resulting initial central temperature was T 0 = 630 MeV at τ 0 = 0.25 fm and the shear viscosity to entropy density ratio was η/s = 0.159. All other transport coefficients, such as the bulk viscosity, were self-consistently determined in terms of η/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 [60,61]. For all results presented herein, we used a one-dimensional lattice with NUM = 2048 points and L = 40 GeV −1 with a corresponding lattice spacing of a 0.0195 GeV −1 . The Crank-Nicolson evolution time step was taken to be dt = 0.001 GeV −1 . This spatiotemporal discretization was shown in prior work to have only small lattice spacing and finite size effects [62]. 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 [74].
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 α s , and the transport coefficients κ and γ. We work in the bottom sector and, as such, take M = m b = m Υ(1S) /2 = 4.73 GeV with m Υ(1S) from ref. [64]. We compute the Bohr radius a 0 from the 1-loop relation and evaluate α s at the inverse of the Bohr radius using the 1-loop running with N f = 3 flavors and Λ obtaining α s (1/a 0 ) = 0.468. The transport coefficients κ and γ are fixed from direct and indirect lattice measurements, respectively. In ref. [75], the heavy quark momentum diffusion coefficient [76,77] was measured directly in a quenched lattice simulation over a wide range of temperatures.
Here we identify κ(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κ(T ) = κ/T 3 which we denoteκ L (T ),κ C (T ), andκ U (T ) and which correspond to the lower, central, and upper "fit" curves, respectively, of fig. 13 of ref. [75]; we use this variation as an estimate of our systematic uncertainty due to κ. As no direct lattice measurements of γ 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 δM (1S) = (3/2)a 2 0 γ, where δM (1S) is the in-medium mass shift of the 1S state, cf. eq. (78) of ref. [42] and ref. [43]. From this relation, γ is accessible from unquenched lattice measurements of δM (1S). In ref. [43], unquenched lattice measurements of δM [Υ(1S)] from refs. [78,79] were used to place bounds onγ(T ) = γ/T 3 of approximately −3.5 ≤γ ≤ 0. We perform simulations usingγ in this range and use this as an estimation of our systematic uncertainty due to γ. We note that more recent lattice studies [80,81] seem to favor δM (Υ(1S)) 0 and thuŝ γ 0.
To initialize the simulation, we assume that at τ = 0 fm the bottomonium wave function is in a singlet state with a smeared delta function initial condition. 3 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., 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. [62]) 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 (κ =γ = 0) from τ = 0 fm to τ 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. [60,61]. To assess the dependence on this assumption, in appendix E, we present results obtained with τ med = 0.25 fm, which corresponds to the earliest time for which the anisotropic hydrodynamic simulation results are available.

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 Γ s→o given in eq. (2.46) onto 1S, 2S, and 3S states, we find where E = 1/(M a 2 0 ) 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.

Comparisons of LO and NLO survival probabilities with complex H eff
In this subsection, we investigate the effect of the NLO evolution terms on the survival probabilities of the Υ(1S) state. To model the medium evolution, we follow the procedure of ref. [60] 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. [60].
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.  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.

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 Υ(1S), Υ(2S), and Υ(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 κ and γ 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. [60,61], we find that at NLO the effect of the stochastic jumps is much smaller than the theoretical uncertainties stemming from κ and γ variation. For R AA , the evolution with H eff provides a reasonable estimation at a significantly reduced computational cost.

NLO H eff results for R AA with sampled physical trajectories
In this subsection we present results of H eff evolution using an ensemble of 80,000 sampled physical trajectories for each set of values forκ(T ) andγ. Unlike the results presented in the previous two subsections, here we compute the survival probability on a trajectoryby-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. [32]).

Excited state feed down
The Lindblad equation, which the QTraj code solves, gives the probability of a heavy quarkantiquark 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, σ exp = F σ 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 for i > j, (4.6) where the branching fractions are taken from the Particle Data Group [64] (see also eq. (6.4) of ref. [60]).
Finally, the nuclear modification factor R i AA for state i is given by where S(c, p T , φ) 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 φ the azimuthal angle. We use the integrated experimental cross sections σ 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. [4,82] as explained in sec. 6.4 of ref. [60]. QTraj -Υ(1S)  , ATLAS [3], and CMS [4,11] collaborations.
NLO predictions for R AA using H 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 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 τ med = 0.6 fm. 4 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 eff without quantum jumps is a reliable approximation for the suppression of the Υ(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κ in the setκ ∈ {κ L (T ),κ C (T ),κ U (T )} while holdingγ = −2.6. This value ofγ was chosen as to best reproduce the collected Υ(1S) suppression data. In the right panel of fig. 4, we show the variation ofγ in the range −3.5 ≤γ ≤ 0. As in the left panel, the solid line corresponds toγ = −2.6. In both the left and right panels, the experimental data presented are from the ALICE [2], ATLAS [3], and CMS [4,11] collaborations.
Compared to the LO results presented in ref. [61], we find that the variation of R AA withκ is reduced, while the variation withγ is increased. The former is due to the fact that, for fixedκ, going from LO to NLO reduces the singlet to octet decay width. As a result, varyingκ over a fixed range results in a smaller variation of the dimensionful singlet to octet decay width. Regarding theγ 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κ andγ, we find that our NLO predictions with T f =190 MeV and our LO predictions with T f = 250 MeV presented in ref. [61] are in good agreement, with the NLO result being in slightly better agreement with data, particularly at low N part . In both figs. 4 and 5, we find excellent agreement between our NLO H eff predictions and the experimental data for the R AA [1S], however, for the 2S and 3S excited states, the H 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 eff (no jump) evolution. We plan to present results including the effects of quantum jumps in a followup paper.

Conclusions and outlook
In the paper, we have studied bottomonium suppression by combining pNRQCD with an open quantum system framework. Working in the regime M v T E, we have gone beyond prior studies [60,61] 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 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 κ since the magnitude of the decay width decreases when going from LO to NLO; however, the sensitivity to γ is slightly increased even though the NLO corrections do not affect the mass shift. 5 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 [46]. 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 [40,83]. However, we do observe that NLO corrections make the Υ(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.

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 κ and γ. The transport coefficient κ is defined in terms of the chromoelectric correlator as In ref. [42], the dispersive counterpart of κ, denoted γ, was first defined as 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 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 [46,84] where D > (ω) is the Fourier transform of the chromoelectric correlator, Then we define the finite temperature correlators where H is the Hamiltonian, and Z(T ) is the partition function defined by That e −H/T is the time evolution operator in imaginary time, i.e., e −H/T E a i (t, 0)e H/T = E a i (t + i/T, 0), implies the Kubo-Martin-Schwinger relation , and setting t = 0, we find Inserting the last equality into eq. (A.4), we have Taking κ as written in eq. (A.1), we can bring it to the form Finally, from eq. (A.11), we obtain an expression for the integral of the chromoelectric correlator times one power of t in terms of the transport coefficient κ

B Conditions for obtaining a Lindblad equation
We aim to write a Lindblad equation equivalent to the master equation given in eq. (2.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 n i ) 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. (2.15)), the problem of reducing the master equation given in eq. (2.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 1 n,m=0 It therefore exists if the equation 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. (B.2) has a solution. First, we examine the case when L 0 and L 1 are linearly dependent, i.e., 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., where 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 = 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. (2.26)). In this case, the left-hand side of eq. (B.2) can be written as .
Equivalently, one can solve eq. (B.2) for L 1 as given in eq. (B.7), i.e., whose solution is eq. (B.10) up to corrections of order 2 to C † C.

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 Since the density matrix is radially symmetric, i.e., it possesses no preferred direction, it is diagonal in l and m, i.e., ρ l m ;lm = ρ lm;lm δ ll δ mm .
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 † i C i is written the same in the Cartesian and spherical bases, i.e., where it should be noted that C where in the second line we have used ρ lm;l m = ρ lm;lm δ ll δ mm and in the last line we have used lm| S |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 The only non-trivial calculation is that of the term C i ρC † 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 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 l ||C (1) ||l 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 l ||r (1) ||l
To calculate l ||r (1) ||l , we make use of the Wigner-Eckart relation. To evaluate the reduced matrix element l ||r (1) ||l , we consider the case of q = m = m = 0, i.e., from eq. (C.8), we have l ||r (1) ||l = l 0|z|l0 l0; 10|l 0 , (C.10) where z = r cos θ can be expressed in terms of spherical harmonics using cos θ = (4π)/3 Y 10 . By using the following identities for the spherical harmonic functions and Clebsch-Gordon coefficients, and l0; 10|l 0 = l + 1 2l + 1 δ l ,l+1 − l 2l + 1 δ l ,l−1 , (C.12) we obtain l ||r (1) ||l = l 0|r cos θ|l0 l0; 10|l 0 To evaluate the reduced momentum operator l ||p (1) ||l we need to evaluate l ||p (1) ||l = l 0|p z |l0 l0; 10|l 0 . (C.14) The momentum operator p z = −i∂ z in spherical coordinates is given by The calculation of the term cos θ ∂ r follows that of r cos θ and reduces to l 0| cos θ∂ r |l0 l0; 10|l 0 Evaluating the second term using the relation where x = cos θ 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 to write (C.22) Using eq. (C.12), this becomes l 0| sin θ ∂ θ |l0 = l l + 1 and, furthermore, Combining eqs. (C.14), (C. 15), (C. 16), and (C.24), we obtain We note that the reduced matrix elements l ± 1||p (1) ||l are not Hermitian. To derive the Hermitian conjugate of the p ↑/↓ reduced matrix elements, one can use the following result for the Hermitian conjugate of the partial derivative operator −i∂ r in the radial basis We use the following identity for simplifying sums over q and m (which can be shown for any general vector operator C (1) ) m,q lm| C (1) q ρ( r, r , t)( where C ↑,↓ = l ± 1||C (1) ||l . Finally, using eqs. (C.13), (C. 25) and (C.27), we obtain (C. 28) In summary, we have shown that starting from the three-dimensional Lindblad equation given in eq. (2.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. (C.26) 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 ψ(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., dr r 2 ψ * (r)Âφ(r) ≡ dr u * ψ (r)Â u u φ (r) . (C.29) For example, the radial derivative operator −i∂ r , which is not Hermitian in the radial wave function basis, is Hermitian in the reduced wave function basis (u-space) which gives [−i(∂ r + 1/r)] † = −i(∂ r − 1/r) ⇒ (−i∂ r ) † = −i∂ 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 (ψ-space) (C.34)

D Transport coefficient κ
As has been pointed out in refs. [85,86], the transport coefficient κ 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 κ, as defined in eq. (A.1), has not been calculated directly. Indirect determinations based on the pNRQCD definition of the quarkonium thermal width, which is proportional to κ, 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 [43]. 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. [60,61] and parameterize the temperature dependence of κ based on the lattice determination of the heavy quark momentum diffusion coefficient done in ref. [75]. Using only indirect lattice determinations of κ would not significantly change the central value of our results but somewhat enlarges the uncertainties due to the variation of κ. In fig. 6, we estimate these errors by running simulations withκ = 0.24 and 4.2 (withγ = 2.6) as taken from eq. (36) of ref. [43]. We observe larger uncertainty compared to theκ variation presented in figs. 4 and 5, but, even with this largerκ variation, we find smaller uncertainty compared to theγ variation.

E Dependence on the medium initialization time
In this appendix, we consider the dependence of our results on the assumed medium initialization time τ med . In the body of the paper, we assumed τ med = 0.6 fm. Here we consider the effect of an early medium initialization time by choosing instead τ 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 part and p T , respectively, including the effect of excited state feed down. In these figures, the central solid line corresponds tô the central value ofγ, both choices for τ med provide a good description of the data. From the comparison with the case of τ med = 0.6 fm presented in the main body of the paper, we see that the uncertainty of our inferred central value ofγ is ∆γ ∼ 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 [87,88]. In the future, it may be possible to use methods such as those presented in ref. [89] to include the effect of early-time momentum-space anisotropies in a systematic fashion.