The melting and abundance of open charm hadrons

Ratios of cumulants of conserved net charge fluctuations are sensitive to the degrees of freedom that are carriers of the corresponding quantum numbers in different phases of strong interaction matter. Using lattice QCD with 2+1 dynamical flavors and quenched charm quarks we calculate second and fourth order cumulants of net charm fluctuations and their correlations with other conserved charges such as net baryon number, electric charge and strangeness. Analyzing appropriate ratios of these cumulants we probe the nature of charmed degrees of freedom in the vicinity of the QCD chiral crossover region. We show that for temperatures above the chiral crossover transition temperature, charmed degrees of freedom can no longer be described by an uncorrelated gas of hadrons. This suggests that the dissociation of open charm hadrons and the emergence of deconfined charm states sets in just near the chiral crossover transition. Till the crossover region we compare these lattice QCD results with two hadron resonance gas models --including only the experimentally established charmed resonances and also including additional states predicted by quark model and lattice QCD calculations. This comparison provides evidence for so far unobserved charmed hadrons that contribute to the thermodynamics in the crossover region.

Ratios of cumulants of conserved net charge fluctuations are sensitive to the degrees of freedom that are carriers of the corresponding quantum numbers in different phases of strong interaction matter. Using lattice QCD with 2+1 dynamical flavors and quenched charm quarks we calculate second and fourth order cumulants of net charm fluctuations and their correlations with other conserved charges such as net baryon number, electric charge and strangeness. Analyzing appropriate ratios of these cumulants we probe the nature of charmed degrees of freedom in the vicinity of the QCD chiral crossover region. We show that for temperatures above the chiral crossover transition temperature, charmed degrees of freedom can no longer be described by an uncorrelated gas of hadrons. This suggests that the dissociation of open charm hadrons and the emergence of deconfined charm states sets in just near the chiral crossover transition. Till the crossover region we compare these lattice QCD results with two hadron resonance gas models -including only the experimentally established charmed resonances and also including additional states predicted by quark model and lattice QCD calculations. This comparison provides evidence for so far unobserved charmed hadrons that contribute to the thermodynamics in the crossover region.

Introduction
Bound states of heavy quarks, in particular the charmonium state J/ψ and its excitations as well as the heavier bottomonium states, are sensitive probes for deconfining features of the quark-gluon plasma (QGP) [1]. Different excitations of these states are expected to dissolve at different temperatures in the QGP, giving rise to a characteristic sequential melting pattern [2]. Recent lattice QCD calculations of thermal hadron correlation functions suggest that certain quarkonium states survive as bound states in the QGP well beyond the pseudo-critical temperature of the chiral crossover transition T c = (154 ± 9) MeV [3]; the J/ψ and its pseudo-scalar partner η c disappear at about 1.5T c [4], while the heavier bottomonium ground states can survive even up to 2T c [5,6].
Light quark bound states, on the other hand, dissolve already at or close to the pseudo-critical temperature, T c , reflecting the close relation between the chiral crossover and deconfinement of light quark degrees of freedom. This leads to a sudden change in the bulk thermodynamic observables and is even more apparent in the behavior of fluctuations of conserved charges, i.e. baryon number, electric charge or strangeness [7,8]. The sudden change of ratios of different moments (cumulants) of net-charge fluctuations and their correlations in the transition region directly reflects the change of degrees of freedom that carry the relevant conserved charges. The total number of hadronic degrees of freedom, i.e. the detailed hadronic mass spectrum also influences bulk thermodynamics. For instance, the strong rise of the trace anomaly (ǫ − 3P )/T 4 , found in lattice QCD calculations may be indicative for contributions of yet unobserved hadron resonances [9].
Recently it has been shown that the large set of fourth order cumulants of charge fluctuations and cross-correlations among fluctuations of conserved charges allows for a detailed analysis of the change from hadronic to partonic degrees of freedom in different charge sectors [10]. For instance, changes of degrees of freedom in the strange meson and baryon sectors of hadronic matter can be analyzed separately by choosing appropriate combinations of charge fluctuation observables. This led to the conclusion that a description of strong interaction matter in terms of uncorrelated hadronic degrees of freedom breaks down for all strange hadrons in the chiral crossover region, i.e. at T 160 MeV [10], which suggests that strangeness gets dissolved at or close to T c . This finding has been confirmed with the analysis presented in [11].
A more intriguing question is what happens to the charmed sector of the hadronic medium at the QCD transition temperature. While it seems to be established that charmonium states, i.e. bound states with hidden charm, still exist in the QGP at temperatures well above T c , this may not be the case for heavy-light mesons or baryons, i.e. open charm mesons (D, D s ) [12,13] or charmed baryons (Λ c , Σ c , Ξ c , Ω c ). To address this question we calculate cumulants of net-charm fluctuations as well as correlations between moments of net-charm fluctuations and moments of net baryon number, electric charge or strangeness fluctuations. Motivated by the approach outlined in Ref. [10] we analyze ratios of observables that may, at low temperature, be interpreted as contributions of open charm hadrons to the partial mesonic or baryonic pressure of strong interaction matter. We show that a description of net charm fluctuations in terms of models of uncorrelated hadrons breaks down at temperatures close to the chiral crossover temperature. We furthermore show that at low temperatures the partial pressure calculated in the open charm sector is larger than expected from hadron resonance gas (HRG) model calculations based on all experimentally measured charmed resonances as listed in the particle data tables [14]. It, however, agrees well with an HRG based on charm resonances from quark model [15,16,17,18] and lattice QCD calculations [19,20,21]. This points at the existence and thermodynamic importance of additional, experimentally so far not established, open charm hadrons.

The charmed hadron resonance gas
While light quark fluctuations can be quite well described by a hadron resonance gas [22] built up from experimentally measured resonances that are listed in the particle data tables [14] it is not at all obvious that this suffices in the case of the heavy open charm resonances. The particle data tables only list a few measured open charm resonances. Many more are predicted in the relativistic quark model [15,16,17,18] and lattice QCD [20,21] calculations. In fact, the large set of excited charmed mesons and baryons found in lattice QCD calculations closely resembles the excitation spectrum predicted in quark model calculations. It is expected that many new open flavor states will be detected in upcoming experiments at Jefferson Laboratory, FAIR and the LHC [16,23,24,25]. If these resonances are indeed part of the charmed hadron spectrum of QCD, they become excited thermally and contribute to the thermodynamics of the charmed sector of a hadron resonance gas. They will show up as intermediate states in the hadronization process of a quark-gluon plasma formed in heavy ion collisions and influence the abundances of various particle species [26]. Heavy-light bound states also play an important role in the break-up of quarkonium bound states. In lattice QCD calculations their contribution becomes visible in the analysis of the heavy quark potential where they can help to explain the non-vanishing expectation value of the Polyakov loop at low temperatures [27,28].
In order to explore the significance of a potentially large additional set of open charm resonances in thermodynamic calculations at low temperature we have constructed HRG models based on different sets of open charm resonances. In addition to the HRG model that is based on all experimentally observed charmed hadrons (PDG-HRG), we also construct an HRG model based on a set of charmed hadrons calculated in a quark model (QM-HRG) where we used the charmed meson [17] and charmed baryon [18] spectrum calculated by Ebert et al. 1 .
One may wonder whether all the resonances calculated in a quark model exist or are stable and long-lived enough to contribute to e.g. the pressure of charmed hadrons. However, as highly excited states with masses much larger than the ground state energy in a given quark flavor channel are strongly Boltzmann suppressed, they play no significant role in thermodynamics. For this reason we also need not consider multiple charmed baryons or open charm hybrid states that have been identified in lattice QCD calculations [20,21] but generally have masses more than (0.8-1) GeV above those of the ground state resonances. We explore the impact of such heavy states by introducing different cut-offs to the maximum mass up to which open charm resonances are taken into account in the HRG model. For instance, QM-HRG-3 includes all charmed hadron resonances determined in quark model calculations that have masses less than 3 GeV.
We calculate the open charm meson (M C (T, µ)) and baryon (B C (T, µ)) pressure in units of T 4 , such that the total charm contribution to the pressure is written as P C (T, µ)/T 4 = M C (T, µ) + B C (T, µ). As the charmed states are all heavy compared to the scale of the temperatures relevant for the discussion of the thermodynamics in the vicinity of the QCD crossover transition, a Boltzmann approximation is appropriate for all charmed hadrons, Here, µ = (µ B , µ Q , µ S , µ C ),μ ≡ µ/T and g i are the degeneracy factors for the different states with electric charge Q i , strangeness S i and charm C i . Results from calculations of open charm meson and baryon pressures using different HRG models are shown in Fig. 1. The influence of additional states predicted by the quark model is clearly visible already in the QCD crossover transition region. At T c , differences between PDG-HRG (dashed lines) and QM-HRG (solid lines) in the baryon sector are as large as 40% while they are negligible in the meson sector. This reflects that the experimentally known meson spectrum is more complete than the baryon spectrum.
In the open charm meson sector, the well established excitations cover a mass range of about 700 MeV above the ground state D, D s -mesons. In the charmed baryon sector much less is known, for instance, experimentally well known excitations of Ξ c range up to 350 MeV above the ground state and in the doubly strange charmed baryon sector only two Ω c states separated by 100 MeV are well established.
As a consequence of the limited knowledge of the charmed baryon spectrum compared to the open charm meson spectrum, the ratio of partial pressures in the baryon and meson sectors differs strongly between the PDG-HRG and the QM-HRG. This is shown in Fig. 1 (top). Significant differences between the QM-HRG-3 and PDG-HRG results also indicate that almost half of the enhanced contributions actually comes from additional charmed baryons that are lighter than the heaviest PDG state. Similar conclusions can be drawn when analyzing partial pressures in the strange-charmed hadron sector or the electrically charged charmed hadron sectors.

Calculation of charm fluctuations in (2+1)-flavor lattice QCD
In order to detect changes in the relevant degrees of freedom that are the carriers of charm quantum numbers at low and high temperatures as well as to study their properties we calculate dimensionless generalized susceptibilities of conserved charges, Here P denotes the total pressure of the system. In the following we also use the convention to drop a superscript in χ BQSC klmn when the corresponding subscript is zero.
For our analysis of net charm fluctuations we use gauge field configurations generated with the highly improved staggered quark (HISQ) action [29]. Use of the HISQ action in the charm sectors includes the so-called ǫ-term and thus makes our calculations free of tree-level order (am c ) 4 discretization errors [29], where m c is the bare charm quark mass in units of the lattice spacing. These dynamical (2+1)-flavor QCD calculations have been carried out with a strange quark mass (m s ) that has been tuned to its physical value and light (u, d) quarks with mass m l /m s = 1/20. In the continuum limit, the latter corresponds to a light pseudo-scalar mass of about 160 MeV. The charm quark sector is treated within the quenched approximation, neglecting the effects of charm quark loops. Within the temperature range relevant for the present study, the quenched approximation for the charm quarks is very well justified. Various lattice QCD calculations using dynamical charm have confirmed that contributions of dynamical charm quarks to bulk thermodynamic quantities, including the gluonic part of the trace anomaly as well as the susceptibilities of light, strange and charm quarks, remain negligible even up to temperatures as high as 300 MeV [30,31]. We note that these quantities directly probe the influence of virtual quark pairs on observables calculated at a fixed value of the temperature. Unlike in these cases there is no simple observable known that would allow us to directly calculate the pressure at fixed temperature. This may be the reason for differences seen in current calculations of the pressure [30,31] using quenched or dynamical charm. In this work, we only use observables that are of the former type and also do not require any multiplicative or additive renormalization.
The line of constant physics for the charm quark has been determined at zero temperature by calculating the spin-averaged charmonium mass [32], 1 4 (m ηc + 3m J/ψ ). For this purpose we used gauge field configurations generated by hotQCD on lattices of size 32 4 and 32 3 · 48 in the range of gauge couplings, 6.39 ≤ β = 10/g 2 ≤ 7.15 [3,22]. On finite temperature lattices with temporal extent N τ = 8, this covers the temperature range 2 156.8 MeV ≤ T ≤ 330.2 MeV. On these lattices and for the slightly largerthan-physical light quark mass value used in our calculations the transition temperature is 158(3) MeV, i.e. about 4 MeV larger than the continuum extrapolated results at the physical values of the light and strange quark masses [3]. We consider this difference of about 3% as the typical systematic error for all temperature values quoted for our analysis, which is not extrapolated to the physical point in the continuum limit.
The line of constant physics for the charm quark sector is well parametrized by with R(β) denoting the two-loop β-function of massless 3-flavor QCD and c 0 = 56.0, c 2 = 1.16 · 10 6 , d 2 = 8.67 · 10 3 . On this line the charm quark mass varies by less than 5%. The ratio of charm and strange quark masses, m c /m s , varies by about 10%, with m c /m s = 12.42 at β = 6.39 and m c /m s = 11.28 at β = 7.15. For most of our calculations we use data sets on lattices of size 32 3 · 8. A subset of these configurations has already been used for the analysis of strangeness fluctuations [10]. These data sets have been enlarged and now contain up to 16700 configurations at the lowest temperature, separated by 10 time units in rational hybrid Monte Carlo updates. Some additional calculations have been performed on coarser 24 3 ·6 lattices, with fixed m c /m s = 12, in order to check cut-off effects also in the charm quark sector. We summarize the statistics exploited in this calculation in Table. 1. We calculate all the moments of net charm fluctuations needed to construct up to fourth order cumulants that correlate net-charm fluctuations with net baryon number, electric charge and strangeness fluctuations. As the calculation of charm fluctuations is fast we can afford to use on each gauge field configuration up to 6000 Gaussian distributed random source vectors for the inversion of the charmed fermion matrix. This leaves us with statistical errors that mainly arise from fluctuations in the light and strange quark sectors where we have used 1500 random source vectors for the inversion of the corresponding fermion matrices.

Partial pressure of open charm hadrons from fluctuations and correlations
Our analysis of higher order cumulants of net charm fluctuations and their correlations with net baryon number, electric charge and strangeness, closely follows the concepts developed for our analysis of strangeness fluctuations [10]. The large charm quark mass, m c ≫ T , however leads to some simplifications. First of all, for temperatures a few times the QCD transition temperature, Boltzmann statistics is still a good approximation for a free charm quark gas. In the high temperature phase we can thus compare our results with cumulants derived from a free massive quark-antiquark gas in the Boltzmann approximation, where we used explicitly the quantum numbers of charm quarks. Another simplification occurs at low temperatures, where we expect a hadron resonance gas to provide a good description of cumulants of net charge fluctuations. At these temperatures, the pressure of the hadronic medium receives contributions from different open charm mesons and baryons. Using the fact that these hadrons carry integer conserved charges for baryon number (|B| ≤ 1), electric charge (|Q| ≤ 2), strangeness (|S| ≤ 2) and charm (|C| ≤ 3), we can separate the total open charm contribution to the pressure in terms of different mesonic (M C ) and baryonic (B C,i with i ≡ |C| = 1, 2, 3) sectors, In this work, we further motivate the decomposition of the open charm pressure in terms of partial pressures in different electric charge and strangeness sectors. In such cases, we decompose the corresponding partial pressures as, (6) Due to the large charm quark mass, the masses of charmed baryons with |C| = 2 or 3 are substantially larger than those of the |C| = 1 hadrons; e.g. ∆ = m C=2 − m C=1 ≃ 1.2 GeV. Even at T ≃ 200 MeV, i.e. well beyond the validity range of any HRG model, the contribution of a |C| = 2 hadron to P C (T, µ)/T 4 thus is suppressed by a factor exp(−∆/T ) ≃ 10 −3 relative to that of a corresponding |C| = 1 hadron. The latter thus will dominate the total partial charm pressure, P C (T, µ)/T 4 ≃ M C (T, µ)+B C,1 (T, µ). Similarly the baryon contributions to the charged and strange partial charm pressures will be dominated by |C| = 1 baryons only.
The dominance of the |C| = 1 sector in all fluctuation observables involving open charm hadrons is immediately apparent from the temperature dependence of second and fourth order cumulants of net-charm fluctuations, χ C 2 and χ C 4 , as well as the correlations between moments of net baryon number and charm fluctuations (BC-correlations). As long as the strong interaction medium can be described by a gas of uncorrelated hadrons these observables have simple interpretations in terms of partial pressure contributions M C and B C,i evaluated at µ = 0, where n, m > 0 and n or n + m are even, respectively. Here and in the following we often omit the arguments of the functions M C (T, 0), B C,i (T, 0). The quantity (χ C 4 − χ C 2 )/12 is an upper bound for the contribution to the pressure from the |C| > 1 channels in the open charm sector. For all temperature values analyzed by us, we find that this quantity is less than 0.2% of χ C 2 . In fact, for temperatures T ≤ 200 MeV the difference vanishes within errors. This may easily be understood as this difference is only sensitive to contributions of baryons with charm |C| = 2, 3; i.e. χ C 4 − χ C 2 = 12B C,2 + 72B C,3 in a gas of uncorrelated hadrons. We thus conclude that up to negligible corrections all cumulants of net-charm fluctuations, χ C n , with n > 0 and even, directly give the total open charm contribution to the pressure in an HRG, P C ≡ P C (T, 0) ≃ χ C 2 . Moreover, each of the off-diagonal BC-correlations, χ BC nm , with n + m > 0 and even, approximates well the partial pressure of charmed baryons, B C ≡ B C (T, 0) ≃ χ BC mn . In Fig. 2 (right) we show lattice QCD data for χ C 4 /χ C 2 . In the crossover region this ratio is close to unity. This confirms that at low temperature the charm fluctuations χ C 2 and χ C 4 indeed are equally good representatives for the open charm partial pressure.

Melting of open charm hadrons
In order to determine the validity range of an uncorrelated hadron resonance gas model description of the open charm sector of QCD, without using details of the open charm hadron spectrum, we analyze ratios of cumulants of correlations between net charm fluctuations and net-baryon number fluctuations (BC-correlations) as well as cumulants of net charm fluctuations (χ C n ). As motivated in the previous section, a consequence of the dominance of the |C| = 1 charmed baryon sector in thermodynamic considerations is that, to a good approximation, BC-correlations in the hadronic phase obey simple relations as, χ BC nm ≃ χ BC
The ratio of any two of these susceptibilities, i.e. χ BC nm /χ BC kl thus will be unity in a hadron resonance gas irrespective of its composition and the details of the baryon resonance spectrum. In Fig. 2 (left) we show the ratio χ BC 13 /χ BC 22 . It clearly suggests that above the crossover region, an uncorrelated gas of charmed baryons does no longer provide an appropriate description of the BC-correlations. Also shown in this figure is the ratio χ BC 11 /χ BC 13 . It is consistent with unity for all temperatures because the relation χ BC 1n = χ BC 11 not only holds in a non-interacting charmed hadron gas (Eq. 8), but also is valid in an uncorrelated charmed quark gas, as is easily seen from Eq. 4. Higher order derivatives with respect to baryon chemical potentials, on the other hand, distinguish between the hadronic and partonic phases. E.g., one finds that for n being odd, χ BC n1 /χ BC 11 = 1 in a hadron gas and 3 1−n in an uncorrelated charm quark gas.
Subtracting any of the BC-correlations from the quadratic or quartic charm fluctuations provides an approximation for the open charm meson pressure in a gas of uncorrelated hadrons. We thus expect for instance, the relation to hold at low temperatures. Their ratio thus should be unity at low temperatures as long as the HRG description is valid. Fig. 2 (right) shows the ratio of the two observables introduced in Eq. 9. It is obvious from the figure that also in the meson sector, an HRG model description breaks down in the crossover region at or close to T c .
The behavior seen in Fig. 2 for correlations between net charm fluctuations and net baryon number fluctuations, in fact, is quite similar to the behavior seen in the strangeness sector (BS-correlations) [10] as well as in the light quark sector which dominates the correlations between net electric charge and net baryon number (BQ-correlations) [22]. In Fig. 3 we show a comparison of ratios of cumulants of such correlations. For the BS and BQ correlations with the lighter quarks we have two additional data points below 156 MeV. In the charm sector we choose a ratio of cumulants involving higher order derivatives in the charm sector as correlations involving only first order derivatives have large statistical errors. These ratios all should be unity in a gas of uncorrelated hadrons. It is apparent from Fig. 3 that such a description breaks down for charge correlations involving light, strange, or charm quarks in or just above the chiral crossover region. These data are taken from Ref. [10,33]. The shaded region shows the chiral crossover region as in Fig. 2. Horizontal lines on the right side show corresponding results for an uncorrelated quark gas. It should be noted that this limiting value is not defined for χ BQ 31 /χ BQ 11 since the denominator as well as the numerator vanishes in perturbation theory up to O(g 4 ).

Abundance of open charm hadrons
We now turn to the analysis of ratios of charge correlations and fluctuations that are, in contrast to the ratios shown in Fig. 2, sensitive to some details of the open charm hadron spectrum. We construct partial pressure components for the electrically charged charmed mesons and the strangecharm mesons, M QC ≃ χ QC 13 −χ BQC 112 and M SC ≃ χ SC 13 −χ BSC 112 , respectively. We also consider the partial pressure of all open charm mesons M C = χ C 4 − χ BC 13 as motivated in Eq. 9. Using these observables we construct ratios with cumulants, which in an HRG receive contributions only from different charmed baryon sectors in the numerator, In an HRG, the first ratio just gives the ratio of charmed baryon and meson pressure, R BC 13 HRG = B C /M C . In the two other cases, the numerator is a weighted sum of partial charmed baryon pressures in charge sectors |X| = 1 and |X| = 2 with X = Q and S, respectively. These ratios are shown in Fig. 4.  (10)). The dashed lines (PDG-HRG) are predictions for an uncorrelated hadron gas using only the PDG states. The solid lines (QM-HRG) are similar HRG predictions including also the states predicted by the quark model of Ref. [17,18]. The dotted lines (QM-HRG-3) are the same QM predictions, but only including states having masses < 3 GeV. The shaded region shows the QCD crossover region as in Fig. 2. The horizontal lines on the right hand side denote the infinite temperature non-interacting charm quark gas limits for the respective quantities. The lattice QCD data have been obtained on lattices of size 32 3 · 8 (filled symbols) and 24 3 · 6 (open symbols).
HRG model predictions for these ratios strongly depend on the relative abundance of the charmed baryons over open charm mesons. Shown in Fig. 4 are results obtained from the PDG-HRG calculation (dashed lines) and the QM-HRG (solid lines). Clearly in the temperature range of the QCD crossover transition, the lattice QCD data for these ratios are much above the PDG-HRG model results. In all the cases, the deviation from the PDG-HRG at T = 160 MeV is 40% or larger. As discussed in Sec. 2, this may not be too surprising as only a few charmed baryons have so far been listed in the particle data tables. The lattice QCD results instead show good agreement with an HRG constructed from open charm meson and baryon spectra calculated in a relativistic quark model [17,18]. The difference in PDG-HRG and QM-HRG model calculations mainly arises from the baryon sector (see Fig. 1). The observables shown in Fig. 4 thus provide first-principles evidence for a substantial contribution of experimentally so far unobserved charmed baryons to the pressure of a hadron resonance gas 3 . This is also consistent with a large set of additional charmed baryon resonances that are predicted in lattice QCD calculations [21].

Conclusions
We have calculated second and fourth order cumulants of net charm fluctuations and their correlations with fluctuations of other conserved charges, i.e. baryon number, electric charge and strangeness. Ratios of such cumulants indicate that a description of the thermodynamics of open charm degrees of freedom in terms of an uncorrelated charmed hadron gas is valid only up to temperatures close to the chiral crossover transition temperature.
This suggests that open charm hadrons start to dissolve already close to the chiral crossover. Moreover, observables that are sensitive to the ratio of the partial open charm meson and baryon pressures as well as their counterparts in the electrically charged charm sector and the strange-charm sector suggest that a large number of so far experimentally not measured open charm hadrons will contribute to bulk thermodynamics close to the melting temperature. This should be taken into account when analyzing the hadronization of charmed hadrons in heavy ion collision experiments.
So far our analysis has been performed by treating the charm quark sector in quenched approximation using fully dynamical (2+1)-flavor gauge field configurations as thermal heat bath. This, in fact, seems to be appropriate for the situation met in heavy ion collisions, where charm quarks are not generated thermally but are embedded into the thermal heat bath of light and strange quarks through hard collisions at early stages of the collision.
We also do not expect that the cumulant ratios analyzed here will change significantly by treating also the charm sector dynamically. This, however, should be verified in future calculations.