Transverse momentum distributions of charmonium states with the statistical hadronization model

Calculations and predictions are presented within the framework of the statistical hadronization model for transverse momentum spectra of the charmonium states J/$\psi$, $\psi(2S)$ and $X(3872)$ produced in nucleus-nucleus collisions at LHC energies. The results are confronted with available data and exhibit very good agreement by using particle flow profiles from state-of-the-art hydrodynamic calculations. For $X(3872)$ production in Pb-Pb collisions we predict a transverse momentum distribution similar in shape to that for J/$\psi$ with a strong enhancement at low transverse momenta and a production yield of about 1\% relative to that for J/$\psi$.


Introduction
Ultra-relativistic heavy-ion collisions are widely used to investigate the evolution from the hot and dense phase of QCD with quarks and gluons as degrees of freedom, the quark-gluon plasma (QGP), towards hadronization into color-neutral objects as degrees of freedom. Analysis of the abundances of the resulting hadrons provides a reliable tool to characterize the phase boundary in the phase diagram of the strong interaction [1]. The statistical hadronization model (SHM) is successfully used to predict and describe hadron abundances produced in relativistic nuclear collisions [2]. The underlying assumption is that at the latest at hadronization the fireball formed in such collisions is close to thermal equilibrium such that hadron yields can be characterized by a grand canonical partition function where baryon number, the 3-component of the isospin, and strangeness are conserved on average and a rapid hadrochemical freeze-out takes place at the phase boundary. These assumptions connect the microscopic content of a heavy-ion collision with macroscopic at-m T CF , the primordial or thermal yields follow the relation m 3/2 exp(−m/T CF ). The total yields include contributions from resonance decays and reproduce particle yields from data over nine orders of magnitude from the lightest mesons up to nuclei with mass number A = 4.
A striking feature in the figure is the difference between the light-flavored loosely bound state hyper-triton 3 Λ H and the charmonium state J/ψ. Although their masses are very similar, the yield of J/ψ is about three orders of magnitude larger. The origin of this enhancement is due to the statistical hadronization of charm quarks, formed in hard initial processes, where the number of charm quarks is conserved throughout the evolution of the collision. Technically, the conservation 1 It should be noted, that throughout this letter, we use c = k B = = 1.  of charm quarks leads to a fugacity in the SHM for charmed hadrons [4] which is, however, not a free parameter but determined by the measured charm cross section. Charm quarks are not confined inside the QGP, thermalize within the QGP and hadronize at the QCD phase boundary into open and hidden charm hadrons. This SHM for charmed hadrons (SHMC) provides an excellent description of charmonium production [5][6][7] without any new parameters and represents compelling evidence, as demonstrated in Figure 1, for this new production mechanism. A more detailed account of the SHMC is given below. Furthermore, a large degree of thermalization is observed in the spectra and the elliptic flow of D-mesons and their decay electrons [8,9]. A number of recent measurements have established the SHMC process (sometimes dubbed '(re)generation') as the dominant production mechanism of J/ψ in heavy-ion collisions at LHC energies [10][11][12][13]. It is therefore appealing and important to extend the intriguing results of J/ψ production beyond yields to particle spectra and to more complex charmonium as well as open charm states to further investigate the SHMC mechanism. In the present publication we focus on charmonium states. Predictions for the open charm sector will be the subject of a future publication.
Loosely bound states such as ψ(2S) and, more dramatically, the potential tetra-quark charmonium state X(3872) are of particular interest. Its observation by the Belle collaboration [14] and the subsequent confirmation by the CDF [15], D0 [16], and BaBar [17] col-laborations showed that it is a narrow charmonium-like resonance and the close vicinity of the particle mass to the D 0D * 0 production threshold suggests that the particle could be a charm meson molecule with a very small binding energy [18]. At the LHC, the state X(3872) was first observed by the LHCb collaboration [19], which later also determined [20,21] its quantum numbers, In addition to the X(3872), the ψ(2S) is a natural choice when expanding the studies on the SHMC mechanism beyond the J/ψ.
In this letter, we present calculations for the yields and transverse momentum spectra of the charmonium states J/ψ, ψ(2S), and X(3872) for heavy-ion collisions at √ s NN = 5 TeV. Results will be presented for the current collision system Pb-Pb and, in the case of X(3872), also for Kr-Kr collisions where much larger luminosities are possible at the LHC.

Heavy quarks in the statistical hadronization model
In the SHMC it is assumed that charm quarks 2 are produced in initial hard scatterings and that during the QGP phase the number of (anti-)charm quarks is conserved, i.e. the thermal production or annihilation are negligible at LHC energies [5]. Charm quarks, together with all other quarks, are screened for T T CF and do not form color-less bound states in the fireball volume V, where µ B is consistent with zero for LHC energies. The quarks thermalize in the QGP before the hadronization and rapid freeze-out at the phase boundary.
The (anti-)charm hadron densities computed in canonical statistical mechanics, n th X , are anchored to the number of produced cc pairs, N cc , by a balance equation where the quantity N cc is interpolated via FONLL [22,23] from charm cross-section measurements [24][25][26][27] in the corresponding rapidity region. Shadowing is taken into account when calculating N cc for nucleusnucleus collisions using rapidity-dependent measurements of the nuclear modification factor in protonnucleus collisions of D-mesons and J/ψ [28,29], where interpolations, if necessary, are done via model calculations [30][31][32]. A canonical suppression factor, I 1 (g c n th oc V)/I 0 (g c n th oc V), is applied to the open charm densities computed in grand canonical statistical mechanics to obtain n th X in Eq. 1. Here, I i with i = {0, 1} are modified Bessel functions and n th oc are the thermal open charm meson densities. This correction factor gains importance for decreasing charm quark densities, i.e. N cc 1. In order to fulfill the balance equation (Eq. 1), the charm quark fugacity g c is needed. As explained before,it is not a free parameter but entirely defined by N cc and n th X . The density profile of the colliding nuclear species is taken into account by a core-corona approach. The 'core' part represents the central region of the colliding nuclei where nucleons undergo many scatterings and are assumed to produce a QGP. The core fraction is normalized to the thermal yields. The 'corona' area includes the nucleons in the outer region of the colliding nuclei. In the overlap zone with on average one or less collisions, no QGP formation is assumed. Rather the yields from the corona fraction are modeled by protonproton (pp) differential distributions scaled by the number of binary nucleon-nucleon collisions in the corona. The radius at which the nucleon density is not sufficient anymore to produce more than one inelastic scattering is taken from the nuclear charge density distributions. This density is found to be approximately at 10 % of the central nuclear density. To get a feeling for the sensitivity to this estimate, we also give the result for a value of 20 %. The fraction of the core and the corona part is estimated by a Glauber simulation [33].
The calculation of the transverse momentum spectra is based on the following consideration: the charm quarks are assumed in local thermal equilibrium in the fireball formed in the collision. At hadronization, i.e. at T CF , the charmonia states then 'inherit' the random thermal motion of the charm quarks superimposed with the collective velocity of the expanding QGP fluid. This velocity is, hence, extracted from state-of-the-art viscous hydrodynamic modelling of light flavor observables.
Specifically, this modelling is based on the (3+1)D viscous hydrodynamic simulation framework MU-SIC [34] with IP-Glasma as initial condition [35]. This framework does not assume boost-invariance. Rather the relevant parameters are adjusted to reproduce the charged particle multiplicity distributions measured by the ALICE collaboration as function of the pseudorapidity η [36,37]. The equation of state and the shear viscosity are parameterized from lattice QCD calculations [38]. The freeze-out hyper-surface is then taken at the chemical freeze-out temperature T = T CF = 156.5 MeV, as determined by thermal fits, see [7].
The results from these hydrodynamic simulations for the velocity profile n and the transverse flow velocity β s T are fed into a blast-wave function formulated for a Hubble-type expansion of the hyper-surface in a boostinvariant scenario as relevant for mid-rapidity data at LHC energy. The result, inspired by earlier calculations reported in [39,40], is: where ρ = atanh(β s T (r/R) n ) and I i and K i are modified Bessel functions with i = {0, 1}. The transverse mass m T = m 2 + p 2 T is obtained using the corresponding particle mass and the temperature is that of chemical freeze-out T = T CF . This functional form (Eq. 2) is then used to model the thermal part of the spectrum and is normalized to the core fraction. To implement the boost-noninvariance discussed above in an approximate way for the description of data at forward rapidity, we adjust the flow velocity according to the (3+1)D hydrodynamic calculations, keeping otherwise the same functional form of the blast-wave description. We have tested numerically that using eq. 2 for large masses such as that of J/ψ yields results very close to those obtained with the standard blast wave description of [39]. The overall normalization factor is obtained from the SHMC as used in Figure 1.
The shape of the corona fraction is modeled by a fit to the transverse momentum spectra measured in pp collision, given by where C, p 0 and n are fit parameters. Since no published data are yet available for the p T spectrum of the J/ψ in pp collisions for mid-rapidity at √ s NN = 5 TeV, an interpolation/extrapolation procedure is applied to model the spectral shape of the corona following [28].
In case of the ψ(2S), the pp cross section is taken from measurements at the corresponding collision energy and rapidity [41]. The transverse momentum spectrum of ψ(2S) in pp collisions at low p T is not measured at mid-or forward rapidity. The shape is approximated by the J/ψ transverse momentum spectrum from [41], are based on a charm cross section at mid-rapidity |y| < 0.9 including shadowing as discussed above. In addition to the full spectrum calculation, the contributions for the thermal core part and the corona are shown. While at low p T the uncertainties are due to the charm cross section, at high p T the uncertainties come from the uncertainty of the corona thickness.
which is consistent with the finding that the production ratio of ψ(2S) and J/ψ are constant over a wide range in transverse momentum [42]. The X(3872) p T -differential cross section is measured at mid-rapidity for 10 < p T /GeV < 30 in pp collisions at √ s = 7 TeV [43]. A universal collision energy scaling is assumed for the transverse momentum spectrum to extrapolate the spectrum from √ s = 7 TeV to 5 TeV for an estimate of the corona shape of the X(3872) spectrum. 'Hybrid' J/ψ spectra in pp collisions are created from the ALICE and CMS results at 5 and 7 TeV to cover the low-and high-p T region simultaneously [44][45][46].These spectra are fitted by the functional form eq. 3 and the ratio of the transverse momentum spectra parameterizations is used to weight the X(3872) p T spectrum at √ s = 7 TeV to estimate the shape of the spectrum at √ s = 5 TeV. Finally, the X(3872) cross section in pp collisions used for the determination of the corona fraction is estimated by an extrapolation from the p T range measured by CMS down to zero p T using Eq. 3 and the corresponding p T spectrum at √ s = 5 TeV.

Results
In Figure 2, the calculations are shown for the transverse momentum spectrum of the J/ψ at mid-rapidity for Pb-Pb collisions at √ s NN = 5 TeV. In addition to the full spectrum calculations the two components, thermal core fraction and corona, are plotted to visualize their contributions. One can see how dramatically different in shape the predicted p T spectrum in the hadronizing QGP is due to the large fugacity, see Eq. 1 and the strong collective transverse expansion. In minimum bias pp collisions such effects are expected to be very small.
In the upper panel of Figure 3, calculations for the transverse momentum spectrum of J/ψ in Pb-Pb collisions at a collision energy of √ s NN = 5 TeV are compared to the measurement by the ALICE collaboration at forward rapidity [11]. The model describes the data for low p T values but falls below the data for p T 4.5 GeV. This suggests additional production mechanisms such as J/ψ production from gluon fragmentation in jets which could contribute towards increasing transverse momentum. In addition, predictions for the ψ(2S) are shown for the same collision system and rapidity window for central collisions. The lower panel shows the production ratio of ψ(2S)/J/ψ.
In Figure 4, predictions for the X(3872) transverse momentum spectrum at √ s NN = 5 TeV are shown for Pb-Pb and Kr-Kr collisions. The branching ratio BR(X(3872) → J/ψπ + π − ) is only known approxi- mately, but can be constrained to a range between 3.2% and about 20% [47]. All results are therefore shown as BR(X(3872) → J/ψπ + π − ) × d 2 σ/dp T dy, where the branching ratio is chosen to 10%. Heavy-ion data taking in LHC Run4 or beyond the long shutdown 4 of the LHC, i.e. after the year 2030, might enable these measurements. We also studied X(3872) production in collisions between lighter beams such as 84 Kr nuclei to assess, in future experiments, the impact of potential increased luminosity vs reduced cross section. In this case, the yields in the Kr-Kr collisions are found to be a factor 4 − 5 lower than for Pb-Pb collisions for p T 5 GeV.

Summary and conclusions
Calculations and predictions using the SHMC for the transverse momentum spectra of different charmonium states are presented and, whenever measurements are available, confronted with data. The spectra calculations are based on the yields obtained from the SHMC as well as a modified blast wave function with input from state-of-the-art hydrodynamical calculations for the flow profiles. Without any further modifications or parameters a very good description of the measured J/ψ transverse momentum spectra in the region p T 5 GeV is obtained. These results further strengthen the evidence that charmonia are produced from deconfined charm quarks at the QCD phase boundary as their yields and spectra are established at the LQCD chiral crossover temperature.
The comparison of the calculations with data for J/ψ in Pb-Pb collisions at √ s NN = 5 TeV at forward ra-pidity shows a good agreement for low p T but falls below the data for increasing momenta (p T 5 GeV). This suggests additional production mechanisms beyond (re)generation at high charmonium momenta.
Predictions are presented for ψ(2S) production in Pb-Pb collisions at √ s NN = 5 TeV at forward rapidity. Statistics collected in LHC heavy-ion data taking in LHC Run 3 might be sufficient for a comparison with data. These results will be important to understand open issues such as the possible presence, in the QGP, of colorless bound states. We further presented calculations for the transverse momentum spectrum of J/ψ mesons produced in Pb-Pb collisions at mid-rapidity which can be tested soon.
Predictions for the production of X(3872) for different centrality regimes of nuclear collisions at √ s NN = 5 TeV are presented for Pb-Pb and Kr-Kr systems. The measurements might become available after the long shutdown 3 of the LHC and will provide a deeper understanding of the mechanism of charmonium and exotica production at the QGP phase boundary and consequently about the QCD dynamics of deconfinement and the hadronization process.