Bottomonium transport in a strongly coupled quark-gluon plasma

TL;DR

This paper develops a semiclassical transport model combining nonperturbative reaction rates and hydrodynamic evolution to study bottomonium suppression and regeneration in quark-gluon plasma, aligning with LHC data.

Contribution

It introduces a nonperturbative, lattice-constrained transport approach that integrates reaction rates with hydrodynamics for bottomonium in heavy-ion collisions.

Findings

Enhanced dissociation and regeneration rates compared to previous models.

The approach describes centrality dependence of bottomonium yields at LHC energies.

Discrepancies observed at large transverse momenta.

Abstract

Quarkonium production in high-energy heavy-ion collisions remains a key probe of the quark-gluon plasma formed in these reactions, but the development of a fully integrated nonperturbative approach remains a challenge. Toward this end, we set up a semiclassical transport approach that combines nonperturbative reaction rates rooted in lattice-constrained $T$-matrix interactions with a viscous hydrodynamic medium evolution. Bottomonium suppression is computed along trajectories in the hydrodynamic evolution while regeneration is evaluated via a rate equation extended to a medium with spatial gradients. The much larger reaction rates compared to previous calculations markedly enhance both dissociation and regeneration processes. This, in particular, requires a reliable assessment of bottomonium equilibrium limits and of the non-thermal distributions of the bottom quarks transported through the expanding medium. Within current uncertainties our approach can describe the centrality dependence of bottomonium yields measured in Pb-Pb ($\sqrt{s_{_{\rm NN}}}$=5.02\,TeV) collisions at the LHC, while discrepancies are found at large transverse momenta.

Equations

$$\frac{\mathrm{d}{N}_{{\rm Y}}(\tau ;p)}{\mathrm{d}\tau }=-{\Gamma }_{{\rm Y}}(T(\tau ),p)N_{{\rm Y}}(\tau ;p),$$

$$\frac{\mathrm{d}{N}_Y(\tau )}{\mathrm{d}\tau }=-{\Gamma }_Y\left(T(\tau )\right)\left[N_Y(\tau )-N_Y^{\mathrm{eq}}\left(T(\tau )\right)\right],$$

$$N_Y^{\mathrm{eq}}={\gamma }_b^2{V}_{\mathrm{FB}}\frac{d_Y}{2{\pi }^2}T{m}_Y^2{K}_2\left(\frac{m_Y}{T}\right),$$

$$N_{b\bar{b}}=\frac{1}{2}{\gamma }_b{n}_{\mathrm{op}}{V}_{\mathrm{FB}}\frac{I_1\left({\gamma }_b{n}_{\mathrm{op}}{V}_{\text{corr }}\right)}{I_0\left({\gamma }_b{n}_{\mathrm{op}}{V}_{\text{corr }}\right)}+{\gamma }_b^2{n}_{\mathrm{hid}}{V}_{\mathrm{FB}},$$

$$V_{\mathrm{corr}}=\frac{4}{3}\pi (r_0+⟨v_b⟩\tau {)}^3,$$

$$\mathcal{R}(\tau )=1-\exp \left[-{\int }_{{\tau }_0}^{\tau }\frac{\mathrm{d}{\tau }^{\prime }}{{\tau }_b}\right],$$

$$V_{\mathrm{act}}(\tau )=\frac{S(\tau )}{s[T_{\mathrm{avg}}(\tau )]},$$

$$S(\tau )={\int }_{T(\tau ,\mathbf{x})>T_{\mathrm{f}}}{\mathrm{d}}^{3}xs[T(\tau ,\mathbf{x})];$$

$$N_{b\bar{b}}=\frac{\mathrm{d}{\sigma }_{b\bar{b}}}{\mathrm{d}y}\frac{N_{\mathrm{coll}}}{{\sigma }_{pp}^{\mathrm{inel}}}\Delta y,$$

$$f_{b\bar{b}}(\tau )=\frac{N_{b\bar{b}}^{\text{in}}(\tau )}{N_{b\bar{b}}^{\text{total}}},$$

$$N_{b\bar{b}}^{\text{eff}}(\tau )=f_{b\bar{b}}(\tau )\cdot N_{b\bar{b}}(\tau ).$$

$$R_{\mathrm{AA}}^Y(N_{\mathrm{part}})=\frac{N_Y^{\mathrm{AA}}(N_{\mathrm{part}};p_{\mathrm{T}})}{N_Y^{pp}(p_{\mathrm{T}})N_{\mathrm{coll}}(N_{\mathrm{part}})},$$

$$\frac{{\mathrm{d}}^3{N}_Y^{\mathrm{reg}}\left(\mathbf{p}\right)}{{\mathrm{d}}^3\mathbf{p}}=C_{\mathrm{reg}}{d}_Y\int {\mathrm{d}}^3{\mathbf{p}}_b{\mathrm{d}}^3{\mathbf{p}}_{\bar{b}}\frac{{\mathrm{d}}^3{N}_b}{{\mathrm{d}}^3{\mathbf{p}}_b}\frac{{\mathrm{d}}^3{N}_{\bar{b}}}{{\mathrm{d}}^3{\mathbf{p}}_{\bar{b}}}$$

$$\text{hfill}\times {\delta }^{(3)}\left(\mathbf{p}-{\mathbf{p}}_b-{\mathbf{p}}_{\bar{b}}\right)w\left(\mathbf{k}\right).$$

Full text

††journal: Physics Letters B

\affiliation

[first]organization= Cyclotron Institute and Department of Physics and Astronomy, Texas A&M University, addressline=, city=College Station, postcode=77843-3366, state=Texas, country=USA

1 Introduction

The fundamental force between a heavy quark and its antiquark is well established by now, essentially following the Cornell potential Eichten et al. [1975] that was originally introduced on phenomenological grounds to describe charmonium and bottomonium spectra. The long-distance linear part of this potential can be considered as one of the most direct manifestations of confinement in the spectroscopy of hadrons. As such, its medium modifications in a quark-gluon plasma (QGP) play a pivotal role in understanding the fate of the confining force in QCD matter, and how this force affects the QGP’s properties. Heavy-quarkonium production in ultra-relativistic heavy-ion collisions (URHICs) is therefore a prime observable to address these questions. Modern transport approaches have reached a fair description of pertinent observables from SPS via RHIC to the LHC to date Rapp et al. [2010], Braun-Munzinger and Stachel [2010], Liu et al. [2011], Mocsy et al. [2013], Ferreiro [2014], Liu et al. [2015], Chen and Zhao [2017], Wolschin [2020], Brambilla et al. [2023], Wu and Rapp [2024], Song et al. [2025], Daddi-Hammou et al. [2025]. However, the microscopic ingredients to these calculations, as well as the medium evolution model required to carry out the transport simulations, vary considerably Andronic and others [2024]. In particular, the connection to the strong-coupling nature of the QGP, as has been deduced from transport calculations of individual charm and bottom quarks (and their hadronization) Beraudo and others [2018], Acharya and others [2022], has not yet been established. Furthermore, constraints on the quarkonium transport parameters that can be obtained from first-principles lattice-QCD computations have not been systematically implemented either.

In the present work we progress in this direction by implementing our recently developed nonperturbative reaction rates Wu et al. [2025], Tang et al. [2024] into a state-of-the-art viscous hydrodynamic simulation of Pb-Pb collisions at the LHC Strickland [2015]. Both components represent a notable advance over our previous work where the quarkonium reaction rates, albeit based on nonperturbative binding energies, were calculated with a perturbative coupling to the medium and applied in a schematic fireball expansion. The reaction rates are based on nonperturbative $T$-matrix interactions that have been thoroughly constrained by recent lattice-QCD data in both heavy- and light-parton sectors Tang et al. [2024]. Pertinent transport coefficients for heavy-quark (HQ) diffusion have been found to be in the range required by open heavy-flavor (HF) phenomenology at the LHC Acharya and others [2022] and agree fairly well with the HQ diffusion coefficient from lattice-QCD Altenkort et al. [2024]. In this way, the nonperturbative rates and the viscous-hydro evolution used in the present work, constitute a largely parameter-free combination of the two main building blocks of a heavy-quarkonium transport model in URHICs, and both components embody the notion of a strongly coupled QGP. That being said, there remain significant uncertainties in our current calculations of observables that we will elaborate on. Specifically, regeneration processes, which couple to the open HF sector (both its kinetics and chemistry), contribute to these uncertainties in several respects. In our applications to experiment, we therefore focus on bottomonium ($Y$) observables, where regeneration yields are believed to be smaller than for charmonia and thus are expected to provide better sensitivity to the inelastic reaction rates.

The remainder of this paper is organized as follows: In Sec. 2 we lay out our implementation of dissociation reactions into trajectories obtained from relativistic hydrodynamics including pre-equilibrium effects associated with nuclear shadowing and bottomonium formation times. In Sec. 3 we set up the rate equation framework for use in hydrodynamics that enables the computation of regeneration reactions in the presence of a spatially non-uniform medium evolution, including escape effects of bottom quarks and their thermal relaxation. In Sec. 4 we compare the results of our updated framework to available bottomonium data in Pb-Pb collisions at the LHC, and in Sec. 5 we summarize and conclude.

2 $Y$ Suppression along trajectories from hydrodynamics

For the bulk medium evolution we employ 3+1-dimensional anisotropic hydrodynamics simulations based on a lattice-QCD equation of state Borsanyi et al. [2010], Bazavov and others [2014] and parameters (e.g., viscosities) tuned to reproduce experimental data for soft-hadron production Strickland [2015], Alqahtani and Strickland [2021]. Focusing on the Pb-Pb (5.02 TeV) system, we initialize $Y$ trajectories in space by an optical Glauber model using a binary-collision profile of hard production, restricted to the rapidity interval of interest. The initial $Y$ momenta are sampled from an azimuthally symmetric distribution in transverse momentum ($p_{\rm T}$) that reproduces measured spectra in proton-proton ($pp$) collisions. Due to their large mass and relatively small binding energies we assume straight-line motion, i.e., neglect elastic rescattering. Each state, $Y$=$\Upsilon(1S,2S,3S)$, $\chi_{b}(1P,2P)$, is then dissociated along its trajectory according to its individual momentum- and temperature-dependent rates Wu et al. [2025] according to the loss term of the rate equation,

[TABLE]

which we integrate down to a final temperature of $T_{\rm f}$=170 MeV (lower values do not significantly affect our final results). The trajectories exhibit a considerable range in temperature and lifetime (cf. the upper panel of Fig. 1 for central Pb-Pb collisions), while the trajectory average is rather similar to previous fireball calculations Du et al. [2017], albeit with a somewhat smaller lifetime in the hydrodynamic environment (cf. the lower panel of Fig. 1 for different collision centralities).

Cold-nuclear-matter effects at the LHC are chiefly expected from nuclear shadowing of the parton distribution functions (PDFs) in the incoming nuclei, which affects the initial hard production of bottomonia. We implement this as in our previous work Wu and Rapp [2024], by adopting a suppression of 20-30% on the integrated $Y$ production at mid-rapidity and 25-35% at forward rapidity, also compatible with the studies of Ref. Du et al. [2017] where 25% and 30% was used, respectively. The $p_{\rm T}$ dependence is taken from nuclear PDFs within the EPPS21 fit Eskola et al. [2022] at 5.02 TeV.

Since typical thermalization times used at the LHC ($\sim$0.2 fm/$c$) are comparable to, or even smaller than, the inverse $Y$ binding energies, bound-state formation times are likely to play a role. Within our semiclassical approach, we follow the standard procedure Gerland et al. [1998] of mimicking the quantum evolution of the forming wave packet as a reduction in the inelastic reaction rate, growing linearly in time from a pointlike $b\bar{b}$ to its asymptotic value as $\tau/\tau_{\rm form}$ for $\tau\leq\tau_{\rm form}$. The quantum formation times are estimated based on the vacuum binding energies as $\tau_{\rm form}=C/E_{B}$, where the coefficient reflects the uncertainty within this approximation; e.g., we find that a variation of $C$ between 1 and 2 causes a $\sim$10% difference in the final results for the $3S$ suppression in peripheral collisions where the formation time is most relevant.

Finally, after evaluating the number of states that survive traversing the QGP we account for late-time feed-down from excited $Y$ states, following Ref. Islam and Strickland [2020].

3 Regeneration in hydrodynamics

To evaluate regeneration within a hydrodynamic background with space-time-dependent temperature and flow, we employ a thermodynamically-motivated approach, building on our previous applications in an isotropic medium Du et al. [2017]. Alternatively, one could pursue a coupled transport approach of diffusing $b$ and $\bar{b}$ quarks with the kinetics of $Y$ suppression and regeneration. Efforts in this direction have been undertaken Yao and Mehen [2019], Yao and Müller [2019], Du and Rapp [2022], but a nonperturbative implementation remains to be developed.

The challenge is to evaluate the additional transport parameter, i.e., the equilibrium limit, $N_{Y}^{\rm eq}\left(T(\tau)\right)$, figuring in the gain term of the kinetic rate equation,

[TABLE]

for each $Y$ state in a spatially non-uniform flowing medium. We neglect (“off-diagonal”) transitions between bound states as they require color-singlet exchanges which are subleading in deconfined matter; in addition, our reaction rates in the sQGP are large and rapidly equilibrate the excited states while the in-medium feeddown to the ground state is thermally suppressed. For a uniform fireball at temperature $T$ and volume $V_{\rm FB}$ one has

[TABLE]

with spin-degeneracy $d_{Y}$, mass $m_{Y}$ and modified Bessel function of the second kind, $K_{2}$. The fugacity $\gamma_{b}=\gamma_{\bar{b}}$ is fixed by the conservation of $b\bar{b}$ pairs in the fireball,

[TABLE]

with open- and hidden-bottom densities $n_{\rm op}$ and $n_{\rm hid}$, respectively, and modified Bessel functions of the first kind, $I_{0,1}$. In the QGP, we include $b$-quarks and the $S$-wave meson resonances in $n_{\rm op}$ as in Ref. Du et al. [2017], but further add $S$-wave baryon states Kaczmarek et al. [2025] (the latter reduce $N_{Y}^{\rm eq}$ equilibrium limit by about $\sim$10%).

In the canonical ensemble, the correlation volume, $V_{\rm corr}$, in $I_{0}/I_{1}$ ensures exact $b\bar{b}$ conservation starting from an initial production volume of the pair. Its time evolution is taken as Grandchamp et al. [2004]

[TABLE]

where $r_{0}=1$ fm characterizes an initial strong-interaction range, and $\langle v_{b}\rangle=0.65\,c$ the average $b$-quark recoil velocity estimated from $B$-meson $p_{\rm T}$ spectra. The sensitivity of the regeneration yield to the parameterization of this volume turns out to be rather weak Du et al. [2017].

Incomplete $b$-quark thermalization is accounted for in a relaxation time approximation Grandchamp and Rapp [2002], via a correction factor,

[TABLE]

which has been shown to produce fairly reliable results Song et al. [2012]. Non-thermal $b$-quark spectra are less favorable for bound-state formation and thus reduce regeneration. The pertinent relaxation time, $\tau_{b}\simeq 7.5$ fm/$c$, is calculated from the same interactions as the $Y$ reaction rates Tang et al. [2024].

To define the “active” fireball volume, $V_{\rm FB}=V_{\rm act}$, of the hydrodynamic medium in which regeneration can take place, we approximate the local environment by the trajectory-averaged temperature, $T_{\rm avg}(\tau)$, sampled along the $Y$ paths as shown in Fig. 1. The active volume is then evaluated from the total entropy, $S(\tau)$, and local entropy density, $s(T)$, in the hydro cell,

[TABLE]

with

[TABLE]

the integration is over the hydrodynamic eigenvolume for temperatures larger than the freeze-out temperature $T_{\rm f}$. In Fig. 2 we display $S(\tau)$ and the resulting $V_{\rm act}(\tau)$ in Pb-Pb (5 TeV) collisions at mid-rapidity for various centralities with $T_{\rm f}$=170 MeV (as in our treatment of suppression). After an initial increase in the build-up phase, the system cools and $S(\tau)$ decreases more rapidly while $V_{\rm act}(\tau)$ grows, reflecting the fireball’s expansion.

The total $b\bar{b}$-pair number is set by the production cross section and the inelastic cross section in $pp$ collisions, multiplied by the number of binary collisions $N_{\rm coll}$ for a given centrality and rapidity bin, $\Delta y=1.8$,

[TABLE]

and corrected for shadowing as described above.

Much like for the bottomonia, $b$ quarks may exit the active fireball volume before the temperature has dropped below $T_{f}$. We estimate the fraction of $b\bar{b}$ pairs that remain within the fireball at each time step,

[TABLE]

by sampling their initial positions according to a Glauber-model distribution and their $p_{\rm T}$ from FONLL distributions, assuming back-to-back kinematics. Initially, all $b\bar{b}$ pairs are inside the fireball (except for very peripheral collisions where they have a finite probability to be produced outside the region with $T\geq T_{f}$). By applying this fraction to the total number of $b\bar{b}$ pairs, we obtain the effective number of $b\bar{b}$ pairs present in the fireball at each time step,

[TABLE]

Regeneration processes become operative once the medium has cooled to the relevant melting temperature, $T_{\rm m}^{Y}$, for each state. In our earlier works, we assumed this temperature to be defined by the vanishing of the pertinent binding energies, $E_{B}^{Y}$ Du et al. [2017]. However, in recent work Tang et al. [2025] it was discovered (by using a complex-pole analysis of the bottomonium $T$-matrices) that resonance-like correlations can persist even beyond the temperature where the nominal binding energy vanishes. Therefore, we have investigated two different scenarios for the onset of regeneration (defined by the lower limit of the integral, $\tau_{\rm m}$, in the gain term of Eq. (2)), namely $\tau_{\rm m}=\tau(E_{B}=0)$ at which $E_{B}^{Y}$ vanishes and $\tau_{\rm m}=\tau(T=T_{\rm m})$ when pole of the $T$-matrix disappears. It turns out that the resulting differences in the final regeneration yields are very small.

4 Predictions for LHC $Y$ data

In our comparisons to experimental data we include error bands where the theoretical uncertainties due to nuclear shadowing and $Y$ formation times, $\tau_{\rm form}^{Y}$ (as discussed above) are added in quadrature. We focus on Pb-Pb collisions at $\sqrt{s_{{}_{\rm NN}}}=5.02\,$TeV, specifically the nuclear modification factor

[TABLE]

where $N_{Y}^{\rm AA}(N_{\rm part})$ is the inclusive $Y$ yield measured in nucleus-nucleus collisions for a given number, $N_{\rm part}$ of participating nucleons, $N_{Y}^{pp}(p_{\rm T})$ is the corresponding yield in $pp$ collisions at the same collision energy, and $N_{\rm coll}(N_{\rm part})$ is the number of primordial binary nucleon–nucleon collisions for a given centrality selection. We discuss the centrality dependence in Sec. 4.1, followed by the $p_{\rm T}$ spectra in Sec. 4.2.

4.1 Centrality dependence

Figure 3 displays our results for $R_{\rm AA}(N_{\rm part})$ for $\Upsilon(1S)$, $\Upsilon(2S)$, and $\Upsilon(3S)$ at both mid- and forward rapidity, compared to CMS Sirunyan and others [2019], Tumasyan and others [2024], ATLAS Aad and others [2023], and ALICE Acharya and others [2019, 2021] data.

At mid-rapidity, displayed in the upper panel, our predictions for $\Upsilon(1S)$, $\Upsilon(2S)$ and $\Upsilon(3S)$ (left, middle and right panel, respectively) agree fairly well with CMS and ATLAS data. The quality of the data description is comparable to our previous calculations based on a perturbative coupling to a quasiparticle QGP Du et al. [2017]. However, the present description is conceptually superior since the calculated rates are rooted in constraints from lattice-QCD and selfconsistently coupled to a strongly interacting medium with broad parton spectral functions representing large collision rates in the QGP liquid. The composition of the yields is also different: the larger $Y$ reaction rates imply significantly more suppression but also more regeneration. The latter becomes the leading contribution for $\Upsilon(1S)$ in central collisions, while it is the dominant source of $\Upsilon(2S)$ and $\Upsilon(3S)$ starting from rather peripheral collisions, for $N_{\rm part}\gtrsim 50$.

This picture essentially persists at forward rapidity, where the partitioning between suppressed primordial and regenerated production is very similar to mid-rapidity, yielding fair agreement with ALICE data except for a tension in peripheral collisions where $\Upsilon(1S)$ production tends to be over-predicted.

4.2 Transverse-momentum spectra

While the predictions for the inclusive yields with the nonperturbative rates seem to work out reasonably well, they do not render a decisive signature compared to other scenarios. But a larger regeneration yield may leave more distinct features in the $p_{\rm T}$ spectra. A reliable calculation of the latter requires to account for the non-equilibrium spectra of the $b$ quarks as they diffuse through the fireball (as mentioned above), which we have thus far not explicitly accounted for. Here, we approximate the $Y$ $p_{\rm T}$ spectra by an instantaneous coalescence model (ICM) Du et al. [2017]

[TABLE]

We employ transported $b$-quark spectra, $\mathrm{d}^{3}N_{b}/\mathrm{d}^{3}\mathbf{p}_{b}$, from HQ diffusion simulations He et al. [2014] at an average regeneration temperature for each $Y$ state, i.e., T=280, 180, 170 MeV for $\Upsilon(1S)$, $\Upsilon(2S)$ and $\Upsilon(3S)$, respectively. In the Wigner function, $w(k)$, we insert the same radii as used in the interference factor of the reaction rates from Ref. Wu et al. [2025]. The overall normalization, $C_{\rm reg}$, is fixed to the values obtained from the rate equation at the end time, $\tau(T_{\rm f})$. The ICM results are combined with the $p_{\rm T}$-dependent suppressed primordial component into the $R_{\rm AA}(p_{\rm T})$ for minimum bias Pb-Pb collisions for inclusive $1S$, $2S$, and $3S$ bottomonia at mid-rapidity and compared to CMS Sirunyan and others [2019], Tumasyan and others [2024], and ATLAS Aad and others [2023] data in Fig. 4. For all 3 states, a maximum structure at relatively low $p_{\rm T}\lesssim 10$ GeV develops which is most pronounced for the $\Upsilon(1S)$, less so for the $\Upsilon(2S)$ and essentially only a shoulder for the $\Upsilon(3S)$. The data are not inconsistent with these structures but do not show compelling evidence in favor of them either. At high $p_{\rm T}$ our calculations tend to underestimate the data for all 3 states, even though the primordial component increases at higher $p_{\rm T}$ due to smaller reaction rates and time-dilated formation time effects. There could be several reasons for that, e.g., longer thermalization times for the onset of the hydro evolution especially in peripheral collisions (currently assumed to be constant at $\tau_{0}$=0.25 fm) or the onset of different production mechanisms Aronson et al. [2018] (e.g., gluon splitting Bedjidian and others [2004]) in connection with an explicit quantum evolution Blaizot and Ollitrault [1987], which might reduce the suppression of the primordial component. Also, space-momentum correlations between the $b$ and $\bar{b}$ quarks, not included here, are known to extend the relevance of regeneration processes He et al. [2022].

5 Summary

We have presented a new calculation of quarkonium transport in heavy-ion collisions within a fully integrated nonperturbative framework. The two main building blocks are recently obtained microscopic reaction rates based on lattice-QCD constraints and a viscous hydrodynamic medium evolution tuned to light-hadron data. This merges the nonperturbative micro-physics of bound-state structure and medium coupling with the macro-physics of a hydrodynamic evolution through the same equation of state in a selfconsistent way, while also establishing consistency with the nonperturbative interactions required to describe the individual diffusion of heavy quarks in the sQGP. While both building blocks did not involve tunable parameters, their current implementation still involves a number of modeling components, e.g., the quarkonium equilibrium limit in evaluating regeneration processes in a spatially non-uniform medium, and the concepts of bound-state formation times and melting temperatures in the early evolution to mimic quantum effects in our semiclassical treatment. We have elaborated them via parameter variations and included the most relevant ones in our final results. We have estimated uncertainties from variations in our input parameters, most notably nuclear shadowing in the primordial production of $Y$’s and $b\bar{b}$ pairs affecting initial conditions and the magnitude of regeneration. Both turned out to be rather modest as part of our theoretical error estimates.

In our applications to experiment we have focused on bottomonium production in Pb-Pb collisions at the LHC. The agreement with the centrality dependence of various $Y$ states turned out to be comparable to previous fireball calculations with perturbative coupling to a quasiparticle medium. However,the much larger nonperturbative reaction rates predict a marked change in the composition of the $Y$ yields, with stronger suppression and a marked increase of regeneration yields. Consequently, the latter become the leading source of $\Upsilon(1S)$ states in semi-/central collisions, while it dominates $\Upsilon(2S)$ and $\Upsilon(3S)$ production starting from peripheral collisions on. This leaves significant footprints in the transverse-momentum spectra, with a maximum in the nuclear modification factor at low $p_{\rm T}$ while toward high $p_{\rm T}$, our results are somewhat below the data. An even stronger footprint might come from the $v_{2}$ spectra. We still expect the $v_{2}$ of the $\Upsilon(1S)$ to be fairly small since its kinetics mostly occurs early in the fireball evolution. However, for all excited states the $v_{2}$ should be quite substantial, since their production is dominated by regeneration relatively late in the evolution where the $b$ quarks are likely to have developed appreciable collectivity.

We envisage a number of future directions and improvements. Clearly, the discrepancy at high momentum needs to be addressed. It will also be critical to test how the nonperturbative micro-/macro-approach fares for charmonia and $B_{c}$ mesons. A quantum transport treatment for the earlier phases is certainly in order; our results suggest that it is vital to accomplish that with a realistic regeneration mechanism encompassing all $b\bar{b}$ pairs in the system, which thus far is not available. Furthermore, a more microscopic treatment of regeneration is in order, to improve the control over the $p_{\rm T}$ dependence in the recombination of the diffusing heavy quarks, in particular also for the elliptic flow. Work in some of these directions is in progress.

Acknowledgments

We thank Jacob Boyd, Sabin Thapa and Ramona Vogt for valuable discussions. This work is supported by the U.S. National Science Foundation under grant no. PHY-2209335 and the Department of Energy via the Topical Collaboration in Nuclear Theory on Heavy-Flavor Theory (HEFTY) for QCD Matter under award no. DE-SC0023547.