Regeneration of bottomonia in an open quantum systems approach

We demonstrate the importance of quantum jumps in the nonequilibrium evolution of bottomonium states in the quark-gluon plasma. Based on nonrelativistic effective field theory and the open quantum system framework, we evolve the density matrix of color singlet and octet pairs. We show that quantum regeneration of singlet states from octet configurations is necessary to understand experimental results for the suppression of both bottomonium ground and excited states. The values of the heavy-quarkonium transport coefficients used are consistent with recent lattice QCD determinations.

The suppression of bottomonium production in heavyion collisions relative to proton-proton collisions provides strong evidence of the formation of a hot, deconfined quark-gluon plasma (QGP) [1][2][3][4][5][6][7][8][9][10][11]. Studies of quarkonium suppression were triggered in 1986 by the idea of Matsui and Satz relating suppression to the Debye screening of interactions in a color-ionized QGP [12]. In recent years, detailed quantum chromodynamics (QCD) calculations have demonstrated that in-medium heavyquarkonium suppression is due to two effects: screening and the generation of a large imaginary part of the in-medium heavy-quarkonium potential, with the latter effect dominating over the former in the temperature regime of interest [13][14][15][16][17][18][19][20]. Although this new understanding represents a fundamental paradigm shift, it has further reinforced the importance of heavy-quarkonium suppression as an ideal probe of the QGP.
Prior studies [21,22] have shown that the in-medium quantum evolution of heavy quarkonium depends on the heavy-quarkonium momentum diffusion coefficientκ and its dispersive counterpartγ, with both defined in QCD in terms of nonperturbative correlators of chromoelectric fields. In the context of the in-medium evolution equations,κ andγ appear as in-medium corrections to the imaginary and real parts of the heavy quarkonium potential, respectively. In this work, we present results obtained using state of the art determinations ofκ [23,24] andγ [24] from unquenched lattice QCD measurements. In contrast with earlier determinations [25][26][27], these recent measurements tend to favor somewhat larger values ofκ and a thermal part ofγ ≈ 0 [28]. A combination of a large in-medium width, related to a largeκ, and a small thermal mass shift, related to a small thermal part ofγ, provide further support for changing the original screening picture of heavy quarkonium suppression, as proposed by Matsui and Satz, into a suppression mechanism dominated by in-medium dissociation due to the thermal width.
Effective field theories (EFTs), specifically potential nonrelativistic quantum chromodynamics (pNRQCD) [29][30][31], allow us to systematically exploit the hierarchies of energy scales of the quarkonium system and the thermal medium. Based on these hierarchies, one can factorize the contributions from the different scales and thus significantly simplify the problem, while still performing the calculations in QCD in a fully quantum setting. Furthermore, the formalism of open quantum systems (OQS) [32] allows us to rigorously treat a quantum system that evolves coupled to and out of equilibrium with a thermal bath. Making use of these two theoretical tools, we are able to describe the evolution of the system, fully taking into account its quantum non-Abelian nature and the out-of-equilibrium evolution [33].
In the remainder of this Letter, we present our method, discussing pNRQCD and the heavy-quarkonium evolution equations, and results including the nuclear modification factor R AA of the Υ(1S), Υ(2S), and Υ(3S) (and double ratios thereof) with respect to the number of participating nucleons N part and transverse momentum p T . We compare our results to measurements of the ALICE [1], ATLAS [2], and CMS [3,10] Collaborations and observe good agreement with the experimental data for all observables. Our main findings are twofold: namely that the experimental data cannot be reproduced in the absence of quantum regeneration and, in accordance with recent lattice measurements, a larger in-medium width and smaller mass shift, i.e., largerκ and smaller |γ|, than previously measured are favored by current bottomonium suppression results.
Methodology -Heavy quarkonium states are characterized by separated energy scales making them ideal for description using EFTs. Heavy-heavy bound states are characterized by the heavy quark mass M and the nonrelativistic relative velocity v ≪ 1. Within this context, M , M v, and M v 2 are referred to as the hard, soft, and ultrasoft scales, respectively. Integrating out the hard scale from QCD gives rise to nonrelativistic QCD (NRQCD) [34,35]. Use of NRQCD to describe heavy-heavy bound states is problematic as both the soft scale characterizing the momentum transfer and the ultrasoft scale corresponding to the binding energy E, remain dynamical, and no unambiguous power counting can be assigned to the operators appearing in the NRQCD Lagrangian without additional assumptions. This issue is remedied by integrating out the soft scale giving rise to potential NRQCD (pNRQCD) [29][30][31].
pNRQCD implements an expansion in the inverse of the heavy quark mass M −1 and the bound-state radius r ∼ (M v) −1 at the Lagrangian level and is thus ideally suited to describe small, heavy-heavy bound states. After integrating out the soft scale, the individual heavy quarks are no longer resolved, and the degrees of freedom in the pNRQCD Lagrangian are heavy-heavy bound states in color singlet and color octet configurations and gluons and light quarks at the ultrasoft scale. Matching pNRQCD onto NRQCD gives rise to non-local potentials at lowest order in the power counting. If we assume that the radius of the system is smaller than 1/Λ QCD , then the potentials are the attractive and repulsive Coulombic potentials in the singlet and octet sectors, respectively. This assumption may apply to the lowest bottomonium states.
In Refs. [21,22] the authors used pNRQCD and the OQS formalism to derive a master equation for a heavy, Coulombic quarkonium realizing the following hierarchy of scales in a strongly coupled (T ∼ gT ) thermal QCD medium (the QGP) where a 0 is the Bohr radius of the bound state. They derived, furthermore, a Lindblad equation [36,37] valid in the additional limit πT ≫ E. References [38,39] solved this Lindblad equation by developing an open source code, called QTraj [40], which implements the quantum trajectories algorithm [41]. Reference [42] used the QTraj code to solve a Lindblad equation accurate up to and including terms of order E/(πT ) allowing for an extension of the region of validity to a lower final temperature T F nearer the pseudocritical temperature T pc of the QGP phase transition. The next-to-leading order (NLO) Lindblad equation is where ρ(t) and H represent the quarkonium density matrix and Hamiltonian including in-medium corrections and the C n i are collapse operators encoding the in-medium width ρ s,o (t) and h s,o (t) represent the heavy quarkonium singlet, octet density matrix and vacuum Hamiltonian; r i and p i are the position and momentum operators, respectively, associated with the radial coordinate; and ∆V so = −∆V os = V s − V o is the difference of the singlet and octet potentials. For further details on the above Lindblad equation and its derivation we refer the reader to Refs. [42,43].
The heavy quarkonium momentum diffusion coefficient κ and its dispersive counterpartγ are defined aŝ whereẼ a i is a chromoelectric field with Wilson lines attached (see Eqs. (2.9) and (2.10) of Ref. [42]) and N c = 3 is the number of colors. We note that the above transport coefficients are in the adjoint representation as obtained in Refs. [21,22,44]. They are given in terms of nonperturbative correlators that should be evaluated in QCD with nonperturbative methods and contain the information about heavy-quarkonium coupling to the QGP.
For the simulation details used to solve the Lindblad equation including temporal and spatial discretizations, initial state, etc., see the introduction of Sec. 4 of Ref. [42]; results in the present work differ only in the values of the transport coefficientsκ andγ and the inclusion of quantum jumps, i.e., quantum regeneration.
In order to accurately describe the in-medium, nonequilibrium evolution of heavy quarkonium, one must include dissociation and recombination. The Lindblad equation accounts for these processes via the collapse operators given in Eqs. (4) and (5), which implement transitions among the different color and angular momentum states. The Hamiltonian and the anticommutator terms in Eq. (2) preserve the quantum numbers of the state while reducing its norm; the first term in parentheses in Eq. (2) changes the quantum numbers of the system and ensures that the overall evolution is trace preserving. The former terms implement dissociation while the latter term implements quantum jumps which are responsible for quantum recombination. Due to computational costs, in our most recent work [42], we implemented only dissociation and observed reasonable agreement with experimental measurements of the nuclear modification factor for the ground state. In this work, we include the effect of quantum jumps and observe this effect to be of critical importance for simultaneously describing the suppression of both the ground and excited bottomonium states.
In all cases discussed herein, we generated approximately 100,000 -200,0000 physical trajectories. When including quantum jumps, we generated approximately 30 quantum trajectories per physical trajectory [48]. We initialized the quantum wave function as a highly-localized Gaussian delta function [38,39,42,49] and terminated the evolution along each physical trajectory when the temperature dropped below T F = 190 MeV [50]. We computed the survival probabilities of the Υ(1S, 2S, 3S) and χ b (1P, 2P ) states by taking the ratio of the final and initial overlap probabilities obtained using vacuum Coulomb eigenstates. We then took into account final state feed down by applying a feed down matrix constructed from data available from the Particle Data Group [51]. The feed down matrix used was the same as in our prior works [38,39,42,49]. We analyzed the agreement of the resulting NLO QTraj predictions for R AA (N part ) with data available from the ALICE, ATLAS, and CMS experiments and found that, without including quantum jumps, it was impossible to simultaneously describe the suppression of the Υ(1S) and Υ(2S) states. Contrarily, when quantum jumps were included, we found that it was possible to quantitatively describe the suppression of both states. In addition, we found that this allowed for a good description of the suppression of the Υ(3S). For details concerning the measure we used to assess agreement with experimental data, we refer the reader to the supplemental material associated with this Letter. Therein, we present plots of the agreement measure as a function ofκ andγ. From our analysis, we concluded that the best description of the data was obtained using the valuesκ ≈ 4 and γ ≈ 0 from our original set of simulations. In all figures, we present results obtained with these values ofκ andγ and assess the effect of including quantum regeneration. In Figs. 1 and 2, we present a comparison of the NLO QTraj results and experimental data as functions of N part and transverse momentum, respectively. The vertical axis in both plots is logarithmic in order to better resolve the suppression of the excited states. In these figures and all subsequent figures, the solid and dashed lines correspond to the results obtained with and without quantum jumps, respectively. The shaded bands indicate the statistical uncertainty associated with the average over quantum and physical trajectories. As can be seen from Figs. 1 and 2, without quantum jumps included, one under predicts R AA for excited states. This is not particular to the case ofκ = 4 andγ = 0 and was observed in our prior simulations in which a temperature-dependentκ was used [38,39,42,49]. When quantum jumps are included, we find excellent agreement with the experimentally measured suppression of both the ground and excited states. For the fully integrated R AA [1S, 2S, 3S], we obtain 0.3598 ± 0.0028, 0.1118 ± 0.0019, and 0.0775 ± 0.0007, respectively.
In Figs. 3 -5, we compare our predictions for the 2S to 1S, 3S to 1S, and 3S to 2S double ratios defined via (n AA [nS]/n AA [n ′ S])/(n pp [nS]/n pp [n ′ S]) as a function of N part . In the case of the CMS data, for the 2S to 1S double ratio, we inferred this observable by taking the ratio of the reported R AA [nS] to R AA [n ′ S] data as no update of the double ratio was reported in Ref. [10]. As a result, for this data set, we do not report the statistical and systematic uncertainties separately. For the ATLAS 2S to 1S double ratio, the black and red error bars correspond to statistical and systematic uncertainties, respectively. In the case of the CMS 3S to 1S double ratio, we inferred this observable from the 3S to 2S double ratio reported by CMS (shown in Fig. 5) and the inferred 2S to 1S double ratio presented in Fig. 3. For the ATLAS 3S to 1S double ratio, we use their reported combined 2S+3S to 1S double ratio [3]. As Figs. 3 -5 demonstrate, this observable is only well explained when quantum regeneration is included. Without quantum regeneration, the 2S to 1S and 3S to 1S double ratios are under predicted and the 3S to 2S double ratio is slightly over predicted. Finally, in Figs. 6 and 7, we present predictions for the 2S to 1S and 3S to 2S double ratios as a function of p T . The data for the 2S to 1S double ratio labeled 'CMS (2023)' were inferred from the reported p T dependence of R AA [1S] and R AA [2S]. For the ATLAS and 'CMS (2019)' data, the Collaborations reported their computed 2S to 1S double ratio directly. As can be seen from these two figures, without jumps, these double ratios are in poorer agreement with the data, again indicating the importance of quantum regeneration for describing the suppression of excited states in heavy-ion collisions. We note, however, that both with and without jumps our NLO QTraj simulations predict a weak dependence of these double ratios on p T . Conclusions -In this Letter we have shown that the inclusion of quantum regeneration allows for a quite satisfactory description of the world's collected data on suppression of bottomonium ground and excited states in a QGP. Leaving out quantum regeneration considerably drifts predictions away from the data. These results confirm in a fully quantum and rigorous QCD setting semi-classical studies of singlet-octet transitions [52][53][54][55][56]. We can thus evolve the singlet and octet density matrix considering the effect of medium-induced color and angular momentum transitions between the states. Looking to the future, it will be interesting to examine the effect of temperature dependence in the heavy-quarkonium transport coefficients, the effect of fluctuating hydrodynamical backgrounds [49], and the possibility of further improvements to the underlying pNRQCD treatment.

SUPPLEMENTAL MATERIAL
In this supplement, we provide details concerning our parameter scan. We considered sixteen combinations, varyinĝ κ ∈ {2, 3, 4, 5} andγ ∈ {−3.5, −2.6, 0, 1}. For each of the sixteen combinations, we computed the quadratic differences between our next-to-leading order result for R AA [1S, 2S] as a function of N part and the central values of the corresponding experimental data. We did not include R AA [3S] in the analysis since only one experiment has reported data for R AA [3S]. For each experiment i ∈ {ALICE, ATLAS, CMS} and state j ∈ {1S, 2S}, we computed the distance measure where k indexes the experimental data points in N part . At each point, we took the central value of the experimental data points and QTraj predictions, with the QTraj predictions being linearly interpolated in N part to match the precise values of N part reported by each experiment. We then constructed an average over experiments for each state ∆ 2 (j) ≡ 1 3 i ∆ 2 (i, j) and a combined measure ∆ 2 ≡ j ∆ 2 (j). These measures are shown in Fig. 8. The top and bottom rows show the results obtained with and without jumps (quantum regeneration), respectively. To generate these figures and the conclusions reported next, the values of ∆ 2 (j) and ∆ 2 were interpolated using cubic and quadratic interpolations inκ andγ, respectively. When including jumps, we found that there was a reasonable overlap in the best descriptions for the 1S and 2S states and a minimum in the combined ∆ 2 of ∆ 2 min = 0.065 at κ = 3.9 andγ = 0. Without jumps, we found that the 1S and 2S best description regions do not overlap, with the 2S description being best whenκ ∼ 2.8 andγ = 1, which results in a poor description of R AA [1S] as a function of N part . Without jumps, the minimal value of the combined ∆ 2 was ∆ 2 min = 0.22 atκ = 3.3 andγ = 1. The resulting best description without jumps does not describe well the suppression of either the 1S or 2S states. Based on these findings, in the main body of the text we compare results obtained from the caseκ = 4 andγ = 0, which are the values closest to the interpolated minimum for ∆ 2 when including quantum jumps.