The multiple-charm hierarchy in the statistical hadronization model

In relativistic nuclear collisions the production of hadrons with light (u,d,s) quarks is quantitatively described in the framework of the Statistical Hadronization Model (SHM). Charm quarks are dominantly produced in initial hard collisions but interact strongly in the hot fireball and thermalize. Therefore charmed hadrons can be incorporated into the SHM by treating charm quarks as 'impurities' with thermal distributions, while the total charm content of the fireball is fixed by the measured open charm cross section. We call this model SHMc and demonstrate that with SHMc the measured multiplicities of single charm hadrons in lead-lead collisions at LHC energies can be well described with the same thermal parameters as for (u,d,s) hadrons. Furthermore, transverse momentum distributions are computed in a blast-wave model, which includes the resonance decay kinematics. SHMc is extended to lighter collision systems down to oxygen-oxygen and includes doubly- and triply-charmed hadrons. We show predictions for production probabilities of such states exhibiting a characteristic and quite spectacular enhancement hierarchy.


Introduction
The statistical hadronization model (SHM) is the standard tool to predict and describe hadron abundances produced in relativistic nuclear collisions [1]. The main physics assumption underlying the SHM is that, near the phase boundary between the quark-gluon plasma (QGP) at high temperature and confined hadronic matter at lower temperature, the fireball formed in such collisions is close to thermal equilibrium. In the large volume limit applicable for Pb-Pb collisions at LHC energies or Au-Au collisions at RHIC energies the produced hadrons can then be precisely described by using a grand canonical partition function based on the hadron-resonance gas (HRG) and with residual interactions deduced using the S-matrix approach of [2]. We note that this HRG statistical operator provides an equation of state that is very close to that emerging from lattice QCD (lQCD) studies in the hadronic phase [3]. Furthermore, the pseudo-critical temperature T pc at µ B = 0, which is now determined in lQCD calculations [4,5] with great precision: T pc = 156.5±1.5 MeV [4], agrees within (small) uncertainties with the chemical freeze-out temperature obtained from the SHM analysis of light-flavour hadron production data [1,2].
How to extend the SHM to the charm sector, i.e., to SHMc, was outlined more than 20 years ago [6] and further developed in [7][8][9][10]. The main idea behind this development is as follows: The charm quark mass m c is much larger than T pc and hence thermal production of charm quarks or hadrons is strongly Boltzmann suppressed. However, with increasing center-of-mass energy the total charm production cross section which results from initial hard collisions increases strongly. If the so produced charm quarks thermalize in the hot fireball they participate in the thermal evolution as 'impurities', their total yield being determined by the charm cross section, not by the fireball temperature. Quantitatively, this is described by the charm balance equation [6,9] leading to a fugacity g c . Roughly from √ s N N > 15 GeV on this will lead to an enhancement of hadrons with charm compared to a purely thermal description, see, e.g., Fig. 1 in [9] and the discussion below. Apart from canonical corrections [7,9] the enhancement scales ∝ (g c ) α where α is the number of charm quarks in a given hadron. Evidence for the thermalization of charm quarks in the fireball is discussed in [2]. Charm quarks are deconfined inside the QGP, thermalize within the QGP and hadronize at the QCD phase boundary into open and hidden charm hadrons. This SHMc was used to predict [7,11] charmonium yields in Pb-Pb collisions at LHC energies long before the LHC turned on. It provides an excellent description of charmonium production [9,[12][13][14] without any new parameters and this success represents compelling evidence for this new production mechanism on the hadronizing QGP phase boundary.
In the present paper we explore the predictions of the SHMc for the production of open charm mesons and baryons. Early predictions for open charm hadrons were made already in [7], and in [8] for baryons with α > 1, but in the absence of experimental data in the relevant low transverse momentum region these early investigations were not pursued further. The situation changed recently when the STAR collaboration at RHIC [15] as well as the ALICE [16][17][18][19] and CMS [20] collaborations at the LHC published first results with Au and Pb beams. It is therefore timely to provide a concise description of the SHMc in the charm sector, to compare results based on this approach to the newly available data and to extend the predictions to the multi-charm sector. We note that the only additional information needed for SHMc predictions are the total open charm cross section and as complete as possible information on the mass spectrum of states in the charm sector. Apart from those there are no free parameters in our approach.
In Section 2 we discuss the SHMc formalism including the charm-balance equation and fugacities, the information on the total open charm cross section, and the hadron mass spectrum in the charm sector. In addition, we will lay out the framework for extending our results to lighter colliding systems of Xe-Xe, Kr-Kr, Ar-Ar and O-O, which could be studied in future runs of LHC. For the study of system size dependence of D meson R AA in a dynamical heavy flavour framework see ref. [21]. For these systems and, in particular for the evaluation of production yields of multi-charm hadrons, a detailed description in terms of canonical thermodynamics is required and outlined. This leads to thermal predictions for rapidity densities of all charmed hadrons in all colliding systems investigated here.
In section 3 we discuss the most up-to date information of the hadron mass spectrum in the charm sector. In particular we review the theoretical and experimental motivation of additional yet-undiscovered charmed hardon states.
In section 4 we present the description of transverse momentum spectra for charmed hadrons using a blast-wave approach. This includes a comparison of results for different freeze-out surfaces. An integral part of this approach is the incorporation of resonance decays into the calculation of spectra. In this section we also outline the 'core-corona' picture which is important to describe the high transverse momentum and centrality dependence of charm hadron production.
Results and comparisons to data are discussed in section 5. In this section we first compare SHMc predictions to data of D-mesons and make a prediction for Λ c baryons. With the same approach and no new inputs aside from masses and quantum numbers of charm hadrons we show how a whole hierarchy of predictions emerges depending on whether we deal with single, double, or triple charm hadrons. Because of the above discussed enhancement of production yields for states with multiple charm these predictions will be tested in the upcoming LHC Run3 and Run4 at least for a selected number of states with α ≤ 2. With a new ALICE3 experiment [22] a large part of the whole mass spectrum of charmed mesons and baryons should be within reach. These experiments can therefore bring completely new information on the degree of deconfinement and mechanism of hadronization of charm quarks in the hot fireball. We conclude this paper with a brief summary and outlook.

Heavy quarks in the statistical hadronization model
Here we recapitulate the physics ideas and formalism behind the SHMc with special focus on the multi-charm sector. For more detail on the original development see [1,7,9]. Our main emphasis will be on the description of yields and transverse momentum spectra for open charm hadrons with α ≤ 3, produced in Pb-Pb collisions at LHC energy. We will also provide expressions to describe the change of yields when going to lighter collision systems including Ar-Ar and O-O and discuss briefly what can be expected. The production of charmonia or charmonium-like states has recently been investigated, see [1,14] and will not be discussed here. Our approach can also be used to make predictions for open charm hadron production at lower energies such as at the RHIC, SPS and FAIR facilities and for higher energies expected at a possible Future Circular Collider [23]. The model can be straightforwardly extended to the beauty sector without conceptual changes or new parameters except for the total open beauty cross section and the corresponding hadronic mass spectrum. However, SHM might need to be modified for beauty hadrons, if future data reveal only partial thermalization of beauty quarks in the QCD medium.

Multi-charm hadrons, charm balance equation and the charm fugacity factor
Our starting point is the charm balance equation [6] where N cc ≡ dN cc /dy denotes the rapidity density of charm quark pairs produced in early, hard collisions and the (grand-canonical) thermal densities for open and hidden charm hadrons are given by n th i,j,k . The index i runs over all open charm states h i oc,1 = D, D s , Λ c , Ξ c , · · · ,Ω c with one valence charm or anti-charm quark, the index j over all hidden charm states h j hc = J/ψ, χ c , ψ , · · · , and the index k over open charm states h k oc,2 = Ξ cc · · · ,Ω cc with two charm or anti-charm quarks. We leave out here states with 3 charm or anti-charm quarks as their contribution to the sum is negligible for realistic masses and values of g c and they have yet to be discovered. These thermal densities are computed using the latest version of the SHMc [1,14] with the chemical freeze-out temperature T cf = 156.5 MeV and the fireball volume per unit rapidity at mid-rapidity V = 4997 ± 455 fm 3 as appropriate for the most central 10% Pb-Pb collisions at LHC energy √ s N N = 5.02 TeV. In the appendix we also give results for the 30-50% centrality interval and at mid-rapidity. Scaling with the measured charged particle pseudo-rapidity density the corresponding volume in this centrality bin is V = 1238 ± 113 fm 3 . For the results shown below, the uncertainties in volume were not propagated, because they are sub-leading compared to the uncertainty in g c discussed below. The total number of charm quark pairs N cc produced in a Pb-Pb collision is a quantity that should be determined by measurement of all hadrons with open or hidden charm. Following this prescription, the only (additional) input parameter of the SHMc, N cc , is determined by experiment. In particular, we note that N cc already includes all nuclear effects in charm production as compared to pp collisions, takes into account potential additions to the charm yield from thermal production in the QGP as well as potential losses due to charm quark annihilation. In practice, using this prescription is, however, difficult since the measurement of all open and hidden charm hadrons needs to be performed without cuts in transverse momentum. Achieving a precision measurement of N cc is one of the priorities for the upgraded ALICE experiment in LHC Run3 and Run4.
In the absence of a measured charm production cross section in Pb-Pb collisions we obtain N cc at mid-rapidity from the measured charm cross section dσ cc /dy in pp collisions by multiplication with the appropriate nuclear thickness function for Pb-Pb collisions and taking into account nuclear modifications. The procedure is described in detail below.
The pp data were measured at √ s = 5.02 and 7 TeV at mid-rapidity [24][25][26][27]. To apply to Pb-Pb collisions, the cross sections are multiplied with the nuclear thickness function and folded with a factor accounting for nuclear modification effects such as shadowing, energy loss or saturation effects. The estimate of this factor is based on the analysis of prompt D 0 and J/ψ production in p-Pb collisions at 5.02 and 8.16 TeV. We used the data from the LHCb collaboration [28][29][30] at forward rapidity, and of J/ψ production at mid-rapidity measured by the ALICE collaboration in pp and p-Pb collisions at 5.02 TeV [26,27]. The √ s = 8.16 and 7.0 TeV data are interpolated to 5.02 TeV using the measured data at other center-of-mass energies and the functional form obtained from perturbative QCD (FONLL) [31]. For mid-rapidity, we obtain a reduction factor of 0.65 ± 0.12, resulting in a value of dσ cc /dy = 0.532±0.096 mb. The corresponding factor for y = 2.0-4.5 is 0.70±0.08 leading to a differential charm production cross section of dσ cc /dy = 0.334 ± 0.053 mb. To obtain the charm quark rapidity density for Pb-Pb collisions of a given centrality, the pp cross section is then multiplied with the mean nuclear thickness function T AA as described in [32]. We neglect in the procedure based on results from pp and p-Pb collisions potential contributions to the differential charm cross section in Pb-Pb collisions from thermal charm production as well as reductions from charm quark annihilation. For LHC both contributions were estimated to be very small and negligible for lower energies [9,33].
We note here that the charm balance equation should contain canonical corrections for more peripheral collisions or for lighter collision systems, i.e., whenever the number of charm pairs is not large compared to 1 [34,35]. The charm balance Eq. 2.1 needs then to be modified accordingly. To that end we define where N oc,1 is the rapidity density of charm quarks bound in hadrons h i oc,1 with one valence charm quark, N oc,2 is the rapidity density of charm quarks bound in hadrons h k oc,2 with two valence charm quarks, and N hc is the rapidity density of charm-(anti-charm) quark pairs bound in hidden charm hadrons h j hc . This defines the total rapidity density of charm quarks, neglecting triply charmed states, as N tot c = N oc,1 + N oc,2 + N hc . Note that the value of N tot c itself depends on the charm fugacity g c . Then the modified charm balance equation using the canonical corrections reads: Here, the I α are modified Bessel functions. For hadrons with 2 or 3 charm quarks there are generally additional terms which are, however, very small because of the small charm densities, and are neglected here (see, e.g. sect. 3.2 in [35]). Solving Eq. 2.3 for g c then determines the charm fugacity factor at 5.02 TeV. For central (0-10%) Pb-Pb collisions and the above discussed differential charm cross section at mid-rapidity (implying dN cc /dy=12.95±2.27) this leads to g c = 29.6 ± 5.2, with the uncertainty determined by the uncertainty in the open charm cross section for Pb-Pb collisions. The rapidity density of open charm hadrons of type h i oc,α with α = 1, 2 charm quarks can then be obtained from the computed thermal densities n th i as : . (2.4) The large value of g c = 29.6 ± 5.2 for central Pb-Pb collisions for charm production at midrapidity (see Fig. 1 in the following section) implies very large enhancements for charmed hadrons compared to what is obtained in the purely thermal case. In the absence of canonical corrections the enhancement factor is (nearly) 900 for doubly charmed, and 2.6 · 10 4 for triply charmed hadrons. For central Pb-Pb collisions at 5.02 TeV the canonical correction factors are in fact close to 1: 0.98, 0.92, and 0.84 for α = 1, 2, 3 charm quarks respectively, for the central value of the differential charm cross section at mid-rapidity, see Fig. 2 below. If these enhancement factors are realized in nature then even very massive triply charmed hadrons may come into reach experimentally. For hidden charm states Eq. 2.4 reduces to The enhancement factors expressed in Eqs. 2.4 and 2.5 come about because of the assumption that all charm quark reach thermal equilibrium at least for temperatures close to T cf . In that case the heavy quarks are completely uncorrelated and the resulting statistical weight is just g α c . We note that this implies deconfinement of the heavy quarks over the volume V , as discussed below.
We also stress that all hadron rapidity densities discussed above are computed as rapidity densities for a volume and hence rapidity window of width of ∆y = 1. The rationale behind this is that one cannot combine charm quarks into hadrons over large rapidity distances as they are causally disconnected because hadrons have a finite formation time τ f ≈ 1 fm and large rapidity correlations can only be established at very early times τ 1 fm [36,37]. The value of ∆y is somewhat arbitrary and a range of ∆y = 1 − 3 was explored in the past and for colliders a weak dependence was found [7]. We finally note the asymptotic form of the modified Bessel functions I α (x). For small argument x and order α this reads: where Γ is the Euler Gamma function. For large x the modified Bessel functions approach This implies that the canonical suppression disappears for large arguments x, i.e., the system has reached the grand-canonical limit. For small x, I 0 ≈ 1 and the canonical suppression factor approaches 1 Γ(α+1) (x/2) α .

Dependence on mass number of the colliding nuclei
In the following we provide information on how to also compute the yields for (multi-)charm hadrons produced in lighter collision systems such as Xe-Xe, Kr-Kr, Ar-Ar and O-O. Of course, these calculations are valid as long as the charm quarks produced in initial hard collisions reach or closely approach kinetic equilibrium in the hot fireball formed in the collision. This has to be carefully checked when one plans to study the production of charm hadrons in such small systems. In addition, we have not included in these exploratory calculations any contributions due to corona effects. Their importance will increase as the colliding systems become smaller. For the system O-O where the nuclear densities never reach a central plateau we expect very substantial corrections which need to be studied carefully if one wants to look for QGP effects in such very light systems. For more discussion on the corona effect see section 4 below.
To understand the charm hadron yield dependence on mass number A of the colliding nuclei we first determine the A dependence of g c . From the charm balance Eqs. 2.1 and 2.3 we note that N cc ∝ A 4/3 since charm is produced in hard collisions and we are interested in central nuclear collisions [38]. Noting further that the volume V ∝ A we immediately obtain that g c ∝ A 1/3 in the grand-canonical limit. In the canonical limit, i.e., for small charm densities, one obtains g c ∝ A −1/3 using the properties of the modified Bessel functions near the origin (see Eqs. 2.6 and 2.7). However, at LHC energies charm densities are not so small and the grand-canonical approximation is a good approximation for the heavier systems Xe-Xe and Kr-Kr and leads to a 20% correction for Ar-Ar. The correction becomes large for the O-O system. In Fig. 1 we show the result of the A dependence of g c as obtained by numerical solution of Eq. 2.3.
The rather strong deviation from the A 1/3 dependence observed for the O-O system is caused by the changes in the canonical correction factor due to the transition from grand-canonical to canonical thermodynamics where the A dependence of g c is expected to approach the A −1/3 scaling as discussed above. For the rapidity range 2.5-4 the nonmonotonic feature of the curves is more pronounced, as the system is deeper into the canonical regime, see Fig. 2.
In Fig. 2 we present the dependence on mass number A of the canonical correction factors f can for the production of charm hadron h i in A-A collisions. They are defined as: . (2.8) The curves on the left and right side are again obtained at rapidity |y| < 0.5 and rapidity 2.5-4, respectively. They are evaluated for charm hadrons with the expression given in equation 2.3. The A dependence of g c needs to be obtained numerically and is displayed in Fig. 1 above.
With the A-dependence of g c and of the canonical corrections factors at hand we can now compute the yield of any charmed hadron in the SHMc as function of mass number A. In section 5 below we will present our results on yields and transverse momentum distributions.
To get a more intuitive understanding of these results we assume, in the following, that the A dependence of g c can be described by the above grand-canonical relation g c ∝ A 1/3 . As can be seen from Fig. 1, this is well fulfilled, at the better than 10% (1%) level, for A ≥ 40 (80). Keeping these small deviations in mind, we can provide a good estimate of the A dependence of charm hadron yields provided we stay with A ≥ 40 , i.e., Ar-Ar collisions, by making use of Eq. 2.4 and the above defined canonical suppression factors f can . This leads to the scaling relation for the production of hadron h i with α charm quarks in collision systems of A-A. Using this relation and the yields for charm hadrons produced in Pb-Pb collisions as displayed in Table 1, see section 5 below, the yields can be computed for charm hadrons yields in lighter systems from Ar-Ar to Xe-Xe. For very light systems such as O-O the full approach as discussed above should always be used. In Fig. 3 the system size dependence of selected hadron yields is displayed for midrapidity (left panel) and forward rapidity (right panel). The band for each hadron species correspond to different charm production cross sections as indicated in the figure. Note the  change in A dependence for open and hidden charm states as a consequence of the absence of the canonical suppression for the latter (compare Eq. 2.5 and 2.4 above).

The canonical volume
The volume V appearing in Eq. 2.1 is usually set equal to the fireball volume at chemical freeze-out V determined by the requirement that the measured rapidity density of charged particles divided by V equals the thermal density of charged particles after strong decays at chemical freeze-out [1]. Employing a connection between momentum rapidity and spacetime rapidity, this volume, corresponding to one unit of rapidity, is a fraction of the entire fireball. To consider such a sub-volume is meaningful since, at high collision energies, equilibration is achieved only locally and not globally. This leads to the picture at freezeout of a string of fireballs lined up in rapidity and filling the entire gap between the rapidities of the two beams (or between beam and target in fixed target mode). The thermal parameters of these fireballs could differ, albeit at LHC we expect a slow variation with rapidity. Only at low collisions energies (AGS energy and below) one should think of one global thermalized system. We note in this context that in [8] it was assumed that the fireball volume comprises all rapidities up to but excluding beam and target rapidities, hence is significantly larger than what is discussed here. When computing the canonical suppression factor f can defined in Eq. 2.8, a new scale enters the problem. To obtain the argument of the Bessel functions, the differential cross section or multiplicity needs to be multiplied with the width of a rapidity interval ∆y which then can be associated with a canonical volume V can over which the relevant quantum number is conserved. For the conservation of baryon number we have recently learned, in the context of net-proton fluctuations, that this volume V can may be significantly larger, not smaller than V [36,39]. Very recent results concerning canonical strangeness suppression [40] at the LHC point also in that direction. Since charm quarks are all produced in the very early phase of the collision we could expect that the canonical volume for charm V can is similarly large, implying a reduced role of canonical suppression and yields larger than computed with V = V can . This would affect in particular predicted yields for multi-charm hadrons from lighter collision systems such as Ar-Ar or O-O. In the numbers given below for (multiple) charm production yields canonical suppression is included. To stay on the conservative side and in the absence of measurements of V can for charm we have, in the following employed only one volume setting V can = V , implying that the canonical corrections for the smallest collision systems could be less severe when more information on V can becomes available.

Charm hadron production and deconfinement of charm quarks
Early on it was realized [7,12,41] that a successful description of the measured yields of charmonia in the SHMc would imply deconfinement for charm quarks. The measurements at RHIC and, in particular, LHC energy lend support to this interpretation [1]. Here we briefly discuss what could be learned on deconfinement from analysis of multi-charm meson and, in particular, baryon production data.
In the SHMc the production of hadrons with α charm quarks is enhanced by a factor (g c ) α compared to what is expected in a purely thermal approach, see Eq. 2.4. Since g c ≈ 30 for central Pb-Pb collisions, the expected enhancements for multi-charm hadron production are very substantial and produce a distinctive hierarchy in their yield pattern, as shown below. That pattern results only if the charm quarks making up the final hadron are uncorrelated prior to hadronization as is expected for fully deconfined ('no strings attached') charm quarks. We note that even the residual correlation imposed by overall baryon number and charm conservation will be very small if the measurement window is of order one unit in rapidity [36].
Production of multi-charm hadrons in the (confined) hadronic phase would also be very small as it would necessarily have to involve exotic multi-particle collisions. To illustrate this point, the following estimates are based on energy conservation and on masses of 4.8 GeV for Ω ccc [42] and 3.62 GeV for Ξ cc [43]. For the most exotic case of Ω ccc production a possible production path is via collisions such as 3D + mπ →p + Ω ccc with m = 3. For the Ξ cc baryon the analogous rate equation reads 2D + mπ →p + Ξ cc with m = 7. But many other processes such as Λ c + D → Ξ cc + π or Λ c + 2D → Ω ccc + π are imaginable. While the rates for all these processes will be enhanced compared to purely thermal estimates by a fugacity factors (g c ) α , they will, nevertheless, be very small because of the low D meson and Λ c density of 1.2 · 10 −3 fm −3 (for D 0 , the highest for D mesons) and 2.6 · 10 −4 fm −3 for g c = 29.6 at chemical freeze-out entering at the same power of α. These rates will fall very rapidly with temperature during the hadronic expansion [44]. Also the phase after chemical freeze-out is by construction not in equilibrium. How to constrain the rate for such multi-particle collisions is totally unclear due to the unknown amplitudes for these different possible many-body collision processes. Similar arguments apply for charmonia, where the dominant channel would be D +D → J/ψ + π. Here, even the extension to ψ involves at least one more unknown parameter. This is to be contrasted with the SHMc approach where there is no free parameters. The experimental observation of a significant number of hadrons with multiple charm in relativistic nuclear collisions hence provides a unique opportunity to test the 'deconfinement' prediction and get quantitative information on the degree of deconfinement achieved in the hot fireball.
The full predictions of the model, including the contribution from the low density corona, are presented for a selection of species in Table 1 for Pb-Pb collisions at 5.02 TeV, for the 0-10% and 30-50% centralities (mid-rapidity values). For these hadrons, the production cross sections in pp collisions have recently been measured by ALICE at midrapidity [18,26,45,46] and those are employed for the calculation of the corona component (we have employed the ratio ψ(2S)/(J/ψ)=0.15 [1]). The model predictions for the core part for all systems for the two rapidity ranges are available in numerical form as auxiliary file with the arXiv version of the publication.

Charm hadron spectrum and SHMc
The spectrum of open charm hadrons incorporated in the SHMc includes all mesons and baryons established experimentally as given by the PDG [43]. This includes 27 D mesons and their anti-particles with angular momenta from 0 to 3 and masses up to 3 GeV. There are 36 established singly-charmed baryons and as many anti-baryons in the mass range up to 3.12 GeV. The known angular momenta are low, mostly 1/2 and 3/2 with one established 5/2 state. The thermal population of the charmed hadrons is strong enough so that the density of the ground state D 0 is quadrupled due to feeding from strong decays, the Λ c density is increased by a factor 5 due to feeding. There has been discussion recently that the number of charmed baryons, in particular, could be significantly larger. Fourth order susceptibilities were constructed and evaluated in lQCD calculations [47] and compared to results from HRG calculations of the same quantities in the temperature range up to the pseudo-critical temperature. The ratios were chosen such that they are particularly sensitive to contributions from the charmed baryon sector in the HRG. It was found that the lQCD results are significantly (at least 40%) above the HRG calculation based on the states established by the PDG in 2012, while adding to the HRG charmed baryon states obtained from a lQCD calculation [48], resulted in good agreement up to the pseudo-critical temperature. The authors of [47] view this as evidence for so far unobserved charmed hadrons contributing to the thermodynamics in the cross over region. Indeed, while the spectrum of [48] is consistent with the number of known states in the mass range above the respective ground state, about 200 additional baryons with total angular momenta up to 7/2 are predicted. Most of these states are significantly higher in mass. For the positive parity states there is a mass gap of about 500-600 MeV, the gap is only of the order of 400 MeV for the negative parity states (that are generally about 300 MeV higher in mass). The situation is only different for the negative parity Ξ c states, where the new states start right at the mass of the highest experimentally established state at 3123 MeV. Accordingly, at a freeze-out temperature T cf = 156.5 MeV the thermal weights are significantly lower. Still, due to their large number and in part also higher degeneracy factors the feeding of ground state charmed baryons could be significantly affected. In this context it is interesting to note that a wealth of new XYZ states were found at the LHC while only 1 additional Λ c , 2 Ξ c and 5 Ω c states were newly discovered (compare e.g. the PDG2012 and PDG2020 compilations).
Triggered by the surprizingly large fragmentation of charm into Λ c measured in pp collisions at 7 and 5.02 TeV by the ALICE collaboration [18,19,49], He and Rapp [50] incorporated into a SHM calculation a hadron spectrum resulting from a relativistic quark model calculation [51] exhibiting a very large number of additional charmed baryons with angular momenta up to 11/2 and both parities. The additional charmed baryons from the RQM calculation have by and large smaller masses than resulting from lQCD [48], falling in part even into the mass range of the known states. Using this charmed baryon spectrum and a temperature of 170 MeV, the authors of [50] find a doubling of the Λ c ground state population as compared to the PDG spectrum and predict a yield in line with the ALICE experimental data.
It should be noted that this poses a conceptual problem because it implies that charmed baryons exist at a temperature significantly above the pseudo-critical temperature for the chiral phase transition, while this is explicitly not supported by lQCD calculations. In [47] it is argued that cumulants on net charm fluctuations indicate that above T pc the charm degrees of freedom are no longer described by an uncorrelated gas of charmed hadrons but that rather the emergence of deconfined charm states sets in just near the chiral cross over transition. On the other hand, Petreczky [52] notes that while the ratio of fourth order baryon-charm susceptibilities around and above the pseudo-critical temperature of the chiral transition is much above the values for the HRG but still below the free quark gas value, that fact could be understood if charm hadron like excitations would still exist above T pc possibly up to 200 MeV. This is not the baseline of the predictions of this publication where deconfinement of all flavors at T pc is assumed. The predictions presented below will provide a stringent test of charm deconfinement and settle this discussion once a large enough dynamic range in mass and charm quantum number is covered by experimental data. Finally we quote recent lQCD results [53] where comparisons of Euclidean correlators to perturbative spectral functions were found to be indicative of charmonium melting in lQCD very close to T pc .
While the questions raised here are debated in the community, we want to give an indication in this publication how the SHMc predictions given below would be affected by a large number of yet undiscovered charmed baryons behaving like simple resonances. To this extent we have performed also calculations where the statistical weight of all excited charmed baryons was tripled and the corresponding change in the predictions by the SHMc is given in section 5 where hadron yields are presented. Finally it should be noted that, even if the above plethora of charmed baryons exists, a treatment as simple resonances in the SHMc could be too naive and a situation could arise similar to the light quark sector. In a recent study [54], the SHM was augmented by 180 nonstrange and 300 strange baryons predicted by lQCD. When they were treated a simple additional resonances, their presence showed a significant impact on particularly the proton yield, strongly deteriorating agreement with experimental data. Proper treatment of the pion-nucleon interaction by the S-matrix approach and using all measured phase shifts [2] completely cancelled out the effect of these additional states. This strong effect of the S-matrix approach could be traced [55] to non-resonant and repulsive components in the pion-nucleon interaction for some partial waves. Whether such a situation could arise in the charm baryon sector depends, among other things, on the widths of the additional states, and is currently completely unexplored. We have assumed that all additional resonances are narrow Breit-Wigner-type resonances.

Transverse momentum spectra of charm hadrons
In the SHM fitted to integrated particle yields no assumption is made about the form of the momentum spectra of produced particles. Therefore the transverse momentum dependence must be supplied by additional modelling of the particle freeze-out.
In hydrodynamical modelling of heavy ion collisions the soft momentum part of particle spectra is obtained by the Cooper-Frye [56] integral over the freeze-out surface and subsequently passing to the hadronic afterburner to perform resonance decays and possible hadronic rescattering. The blast-wave model [57,58] is motivated by the same physics picture, but realized in simpler but approximate way to generate the p T spectra. The thermal particle spectra are obtained from a simple freeze-out surface with a given freeze-out temperature and with parametrized radial velocity profile. This thermal blast-wave model has been used extensively in the past to fit and characterize the experimentally measured identified particle spectra [59][60][61][62].
For boost-invariant and azimuthally symmetric freeze-out surfaces dσ µ , the Cooper-Frye integral can be reduced to a one-dimensional integral along the freeze-out contour in the τ -r plane [57,58]: where 2J + 1 accounts for spin-degeneracy. Here we consider a freeze-out surface defined by a single-valued function τ (r) in the range 0 < r < r max . The freeze-out kernels K eq 1,2 (p T , u r ) can be calculated analytically for the Boltzmann distribution f (p) = exp(− m 2 + p 2 /T ) of initial particles on the freeze-out surface and takes the well-known form in terms of modified Bessel functions [57,58] K eq 1 (p T , u r ) = 4πm T I 0 where m T = m 2 + p 2 T and T is the (constant) freeze-out temperature. The 4-velocity u r = β/ 1 − β 2 is given in terms of radial velocity β(r), which is commonly parametrized by a power function with two parameters β max and n β(r) = β max r n r n max . In this paper the spectra of charmed hadrons formed in the core, i.e. by hadronization of the hot QGP fireball, are evaluated by using the velocity profile from a (3+1)D viscous hydrodynamics code MUSIC with IP-Glasma initial conditions tuned to the light flavor hadron observables [63,64]. The velocity profile and best fit with β max = 0.62 and n = 0.85 for 0-10% centrality bin is shown in Fig. 4 (we use β max = 0.60 and n = 0.85 for 30-50% centrality bin). The fit uncertainties of the parameters β max and n are 0.005 and 0.05, respectively.
Different types of freeze-out surfaces have been used in the past, for example, the constant Bjorken time freeze-out surface introduced in ref. [57] τ (r) = τ fo (4.4) or constant proper time time surface of [65] τ (r) = τ 2 fo + r 2 . (4.5) In ref. [65] the velocity flow was restricted to be a Hubble-like u µ = x µ /τ fo and parallel to the norm of the surface. For parametrized velocity in Eq. 4.3, u µ is no longer proportional to dσ µ . However, one can consider a third type of the surface for which this condition is still true: τ (r) = τ fo + r 0 dr β(r ) and using Eq. 4.3 we get The three freeze-out surfaces are depicted in Fig. 5 (left). Without loss of generality, the freeze-out time is taken to be equal to τ fo = r max and r max itself can be determined by requiring the freeze-out volume per unit rapidity to be equal to a given value, e.g. V = 4997 fm 3 in central Pb-Pb collisions. Note, however, that the integration variable r can be rescaled to x = r/r max with the result that r appears as normalization in front of the integral. Since we replace the overall normalization by that obtained from the SHMc, knowledge of r max is not required, and the only parameters left are the dimensionless parameters β max and n, as discussed above.
As we did in a previous publication for the J/ψ spectrum [14], the spectra for various charmed hadrons are computed using this velocity profile as input for a blast-wave parameterization in terms of temperature, flow velocity profile and mass of the hadron. The temperature we use is the chemical freeze-out temperature T cf = 156.5 MeV obtained from fitting the yields of light flavor hadrons and nuclei as measured by ALICE for Pb-Pb collisions at √ s = 2.76 TeV [1,2]. We studied the effects of the uncertainties of the blast wave parameters β max and n on the hadron spectra. The resulting variations in the spectra are less than 10% and in the ratios to D 0 less than 3%.
In Fig. 5 (right) we show the D 0 spectra for the three freeze-out surfaces. We see that the difference in the absolute spectra is small and lies within the uncertainty band, which is mostly due to uncertainty in g c at these low momenta. In addition, given the still large experimental uncertainties we do not expect the precise form of the freeze-out surface to be the most important factor and we will use a constant freeze-out time surface as the default choice. We emphasize here that for particle ratios, e.g. Λ c /D 0 , even this small difference mostly cancels.
One of the limitations of the standard blast-wave model is that it does not include the momentum modification of particle spectra due to the feed-down caused by resonance decays. Recently, a very efficient way of computing such modifications was derived [67] and applied in blast-wave fits with resonance decay feed-down [68] and hydrodynamic simulations [69]. Here we compute the momentum resolved decay feed-down to long lived charmed mesons and baryons using the FastReso computer code [70]. In total we perform calculations for 76 2-body and 10 3-body decays of charmed mesons and baryons. In practise, this procedure replaces thermal Boltzmann freeze-out kernels in Eq. 4.1 with numerically computed total final particle kernels. We use the same temperature and radial velocity profiles as in a standard blast-wave model. In Fig. 6 (left) we show the full decay spectra of charmed hadrons over their initial thermal spectra. In addition, in Fig. 6 (right)  Figure 6. Left: ratios of different particle spectra with feed-down contribution to thermal spectra (note that the corona contribution is not included here). Dashed lines correspond to the ratio of integrated yields (these ratios were previously used to scale thermal spectra in SHMc). Right: feed-down contribution to Λ + c from different decay channels. For details see text.
we show the selected decay-channel contributions to Λ c spectra. The feed-down contributions preferentially accumulate at low momentum and can be as large as 5 times that of thermal spectra for Λ c . The dotted lines in Fig. 6 (left) show the ratio of full over thermal p T -integrated yields in SHMc. These feed-down factors were used previously to scale the thermal spectra without accounting for p T dependence of the feed-down. One can see rather good agreement between the naive and exact scaling of the spectra for p T 3 GeV, where most of the particles are. As low momentum is the only region where core charmed hadron production is dominant, we find in practice very small differences between full decay spectra and scaled thermal spectra in this momentum range. Nevertheless, in the plots below we will use the spectra obtained with decay kernels from FastReso. Finally, the high momentum power-law tail actually observed in experimental particle spectra is not described by hydrodynamics. Instead it can be modelled using a corecorona picture [14]. Even in nucleus-nucleus collisions at small impact parameter, a number of nucleon-nucleon collisions take place in the so-called corona-region where the overlap density is a small fraction of the maximum density achieved in the collision. In this overlap volume where nucleons undergo on average one or less collisions, we assume that no QGP is formed and, hence, treat the collisions as pp-like. On the contrary, in the core part, we assume full thermalization of produced charm quarks. We define the corona region as having 10% of the central density ρ 0 . In a heavy nucleus at rest the central nucleon number density is ρ 0 = 0.16 fm −3 . The p T shape of the cross section measured in pp collision is parametrized by where the coefficients C, p 0 and n are obtained from a fit to experimental distributions for from [66]. The pp data needed to compute the corona part are taken from [18,45]. The model band width at low and high p T are driven by the uncertainties of g c and pp spectra fits, respectively, as described in the text. each particle species [18,26,45] and the total integral of the function is set to experimentally measured integrated cross section dσ/dy. The fit is found to describe the measured cross sections well within the uncertainties in the whole p T range considered. We then scale the pp differential cross section by the overlap function T corona AA to account for the number of binary nucleon-nucleon collisions in the corona.
In summary, for each of the charmed hadrons under consideration the p T spectra are obtained by summing the soft momentum spectrum from the the blast-wave model with resonance decays and the high momentum tail from the corona part. The uncertainty bands are obtained by varying g c . In addition, the uncertainty on the corona part also includes the uncertainty of the fit to the pp data [18,26,45]. This uncertainty is assumed to be uncorrelated for different particle species and is the dominant source of uncertainties for particle spectra and their ratios at high p T , although it cancels for R AA .

Results for Pb-Pb and lighter collision systems
In the following we will describe predictions from the SHMc as well as the comparison of results from SHMc with the currently available data. For simplicity, we will only consider Pb-Pb collisions at √ s NN = 5.02 TeV and 0-10% centrality, and predictions for 30-50% centrality will be given in Appendix A. The model predictions for all particle species and the two centrality bins are available in numerical form as auxiliary file with the arXiv data [66]. The pp data needed to compute the corona part are taken from [18,26,45]. The model band width at low and high p T are driven by the uncertainties of g c and pp spectra fits, respectively, as described in the text.
version of the publication. By far the best series of experiments exists for D mesons produced in Pb-Pb collisions, see [66].

Transverse momentum distributions
In Fig. 7 we show the comparison between the SHMc predictions and data for spectra and nuclear modification factor R AA as a function of transverse momentum p T . The transverse momentum dependence is obtained as explained in detail in section 4 above. Note that there are no new parameters used here apart from the hydrodynamics input discussed in section 4. The transverse momentum spectrum for D 0 mesons is very well described in particular in the purely thermal ("core") region for p T ≤ 4 GeV. In the transition region between core and corona as well for the high momentum tail we notice that the data are under-predicted for both the p T spectrum and the R AA . This suggests that the corona description is somewhat schematic and could be further optimized. The corresponding distribution for the Λ c baryon are displayed in the lower panels of Fig. 7. We note that these spectra and distributions are obtained with the unmodified charm resonance spectrum discussed below.
In Fig. 8 we show the corresponding distributions for D + , D * + , D + s and Λ c , plotted as a ratio to the D 0 spectrum. In this normalized plot, the charm cross section which determines the charm fugacity parameter g c , is eliminated. For the three D-mesons we observe very good agreement with the experimental observations. For the Λ c baryon the structure of the distribution changes quite strongly: a clear maximum appears near p T = 4.5 GeV. Within  Figure 9. Mass dependence of yields dN/dy for various hadron species for Pb-Pb collisions at midrapidity. The left panel is for absolute yields and the right panel is for yields per degree of freedom (2J + 1). In this plot also the primordial (prior to decays) values are shown as lines, corresponding to hadrons with charm-quark or anti-quark content of 0, 1, 2, and 3 (respective powers of g c ).
the framework of the SHMc this maximum appears as a consequence of a superposition of collective flow (hydrodynamic expansion) and change of hadronization regime from bulk (statistical hadronization) to jets, much as it is observed also for the Λ/K ratio in the (u,d,s) sector [71].

Integrated yields
In this section we discuss results for momentum integrated particle yields, which for constant temperature freeze-out assumed in the SHMc, do not depend on the details of the freeze-out surface and velocity prametrizations discussed in section 4 In Fig. 9 we show the mass dependence of rapidity distributions dN/dy for selected charm hadrons at mid-rapidity. The selection includes D 0 mesons at the lower masses and includes many multi-charm states including the hypothetical Ω ccc baryon at the high mass end of the plot. All are stable against decays via strong interactions. Already the left plot exhibits clear structures whose origin becomes clear with the plot at the right hand side, where the yields are divided by the angular momentum degeneracy. Since we are in the 'Boltzmann' regime where all masses M are much larger then the temperature T cf = 156.5 MeV, the degeneracy-normalized particle yields scale in the SHMc as ∝ M 3/2 exp(−M/T cf ). In a log plot over 7 decades this function looks essentially like a straight line for fixed charm quark number. The color code separates particles with α = 1, 2, 3 charm quarks. The line at the far left corresponds to α = 0 and coincides with that determined for (u,d,s) hadrons in [1]. The deviation clearly visible for α = 1 is due to feeding from hadronically unstable resonances. The grouping into three distinct regions is what is called in the introduction 'the charm hadron hierarchy'.
In Fig. 10 we show the total yields, the sum of core and corona components, for selected hadron species for which the data in pp collisions, used for the calculations of the corona   component, are available. We include in the plot a scenario of charm baryon enhancement, implemented via tripled statistical weights for excited charmed baryons, which leads to an increase of the total thermal charm densities by 18%. Note that the additional charmed baryon resonances are all assumed to be narrow Breit-Wigner-type resonances, as discussed in section 3. We demonstrate that the equivalent increase in the input charm cross section (from 0.53 to 0.63 mb) leads to a significant increase in the predicted yield for the charmed baryons, while the yields of all the rest of the species remain unchanged 1 . The numerical values for the case of the PDG hadron spectrum are shown in Table 1. One notices that some of the uncertainties are asymmetric and this originates either from SHMc, as the g c values are characterized by (slightly) asymmetric uncertainties and from the corona component via the experimental production cross section for pp collisions.
In Table 2 we have compiled the expected luminosity, rapidity density for Ω ccc production, inelastic cross section corresponding to the 10% most central collisions, and expected yields for Ω ccc production in 5 different collision systems at top LHC energy and for a run time of 10 6 s. The beam parameters are from [73], the rapidity densities and yields for Ω ccc production are our predictions. The predictions are per unit rapidity for the 10% most central collisions but contain no efficiency and acceptance corrections. Nevertheless, substantial yields can be expected. Even though the expected luminosity increases by 4 orders of magnitude when moving from Pb-Pb to O-O, the yield in O-O is comparable to that for Pb-Pb, and that at a price of about 10 collisions per bunch crossing for O-O [73]. Furthermore, corona effects will be much increased when going to such a small system. Which of the systems is optimal for QGP-related research will have to be carefully optimized.

Conclusions and Outlook
In the present paper we have explored a range of predictions made within the framework of the SHMc with focus on hadrons with open charm. Most important is the comparison to recent ALICE measurements on D mesons [66] and predictions for Λ c baryons. As baseline for SHMc predictions we kept the chemical freeze-out temperature T cf = 156.5 MeV determined from the analysis of (u,d,s) hadrons. As only additional input we used the open charm cross section based on pp measurements from the ALICE and LHCb collaborations and extrapolated to the Pb-Pb system using hard collision scaling and a correction for nuclear modifications obtained from an analysis of recently measured p-Pb open and hidden charm data. The transverse momentum distributions were obtained in a novel, hydro-inspired approach including resonance decays. Without any further assumptions and parameters all D meson yields and low transverse momentum distributions in Pb-Pb collisions are well described. The situation is less well settled in the Λ c baryon sector. Recent ALICE measurements in pp and p-Pb collisions [19] indicate enhanced production of Λ c baryons compared to what was expected based on e + e − and ep data on fragmentation into charmed baryons. For an account of ALICE preliminary data including those from Pb-Pb collisions see Fig. 4 in [74]. These preliminary data have led to new charm baryon production models including "missing" charm baryons [75]. We have therefore provided predictions for Λ c production in Pb-Pb collisions using the current experimental information on the charm baryon resonance spectrum [43] as well as with an increased number of charm baryons. New data on this puzzling situation are expected soon from both the CERN ALICE and LHCb collaborations.
The success of the description of yields and low transverse momentum spectra of open charm hadrons by the SHMc also demonstrates that the hadronization of open and hidden charm takes place at or close to the QCD phase boundary. It further demonstrates that open and hidden charm data can be reproduced with one common hadronization mechanism.
Our predictions for Pb-Pb collisions imply very large enhancements for hadrons with 2 or 3 charm quarks compared to pure thermal production with charm fugacity g c = 1. The enhancement will be predominantly visible at low transverse momentum p T , see, e.g., Fig. 7. For multi-charmed baryons these enhancements lead to an impressive and quite spectacular hierarchy, see Fig. 9. To test these predictions is a challenge for future charm production experiments in LHC Run3 and Run4 and ultimately one of the important goals for the ALICE3 'all Silicon' experiment [22]. Fundamental new information on the hadronization and deconfinement of charm quarks should be the rewards for the efforts to build such a detector.
Polish Ministry of Science and Higher Education. V.V. is supported by a research grant (Grant No. 00025462) from VILLUM FONDEN, the Danish National Research Foundation (Danmarks Grundforskningsfond), and the Carlsberg Foundation (Carlsbergfondet).  Figure 11. Spectra (left) and R AA (right) of D 0 (top) and Λ c in Pb-Pb collisions at √ s NN = 5.02 TeV and 30-50% centrality. Pb-Pb data for D-meson distributions taken from [66]. The pp data needed to compute the corona part are taken from [18,45]. The model band width at low and high p T are driven by the uncertainties of g c and pp spectra fits, respectively, as described in the text.  to ALICE data [66]. The pp data needed to compute the corona part are taken from [18,26,45]. The model band width at low and high p T are driven by the uncertainties of g c and pp spectra fits, respectively, as described in the text.