Charmed hadron chemistry in relativistic heavy-ion collisions

A solid hadronization model is essential for understanding hadronic observables in high-energy nuclear collisions, while still remains a challenge due to its non-perturbative nature. We develop a comprehensive coalescence model that includes a complete set of both $s$ and $p$-wave hadronic states. With a strict energy-momentum conservation implemented, the boost invariance of the coalescence probability and the thermal limit of the produced hadron spectrum are guaranteed. By combining our hadronization scheme with an advanced Langevin-hydrodynamics model that incorporates both elastic and inelastic energy loss of heavy quarks inside the realistic QGP medium, we provide a satisfactory description of the nuclear modification factors of $D$ mesons and $D$-decayed electrons, together with the $p_\mathrm{T}$-integrated and differential $\Lambda_c/D^0$ and $D_s/D^0$ ratios. Our study indicates that the in-medium size of charmed hadrons should be larger than the size in vacuum. It is also found that the inclusion of the $p$-wave states and the radial flow of the QGP is essential to understand the chemical composition of the charmed hadrons observed in relativistic heavy-ion collisions.

A solid hadronization model is essential for understanding hadronic observables in high-energy nuclear collisions, while still remains a challenge due to its non-perturbative nature. We develop a comprehensive coalescence model that includes a complete set of both s and p-wave hadronic states. With a strict energy-momentum conservation implemented, the boost invariance of the coalescence probability and the thermal limit of the produced hadron spectrum are guaranteed. By combining our hadronization scheme with an advanced Langevin-hydrodynamics model that incorporates both elastic and inelastic energy loss of heavy quarks inside the realistic QGP medium, we provide a satisfactory description of the nuclear modification factors of D mesons and D-decayed electrons, together with the pT-integrated and differential Λc/D 0 and Ds/D 0 ratios. Our study indicates that the in-medium size of charmed hadrons should be larger than the size in vacuum. It is also found that the inclusion of the p-wave states and the radial flow of the QGP is essential to understand the chemical composition of the charmed hadrons observed in relativistic heavy-ion collisions.
Introduction. -Relativistic heavy-ion collisions provide a unique opportunity to study the color deconfined state of nuclear matter, known as the Quark-Gluon Plasma (QGP). Hadronization is a crucial bridge between many theoretical calculations at the partonic level and the final hadronic observables in order to understand the hadron chemistry at low transverse momenta (p T ) [1][2][3][4] and the high-p T hadron and jet spectra [5][6][7]. Heavy quarks are clean probes of the QGP since they are produced in the early stage of nuclear collisions and their flavor numbers are conserved during interactions with the QGP. Extensive studies have been devoted to the production, evolution and nuclear modification of heavy quarks and heavy flavor hadrons in relativistic nuclear collisions [8][9][10][11]. Recent experiments have observed a non-monotonic p T dependence of the D meson nuclear modification factor (R AA ) at low p T [12] and a large baryon-to-meson ratio [13,14] at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC). They call for a solid hadronization model for understanding the detailed production mechanisms of the heavy flavor hadrons in heavy-ion collisions.
Various coalescence approaches have been developed for the hadronization process in high-energy nuclear collisions. Simplified assumptions, such as the equal-velocity coalescence [15] or the coalescence between neighboring quarks [16], have been applied to grasp some of the key features of the charmed hadron chemistry. However, these models lack sophisticated calculations of microscopic coalescence probabilities. References [17,18] have developed a resonant recombination model in which the coalescence probability is connected with the resonant scattering rate of heavy quarks inside the QGP, but the information about the detailed hadron structures has not been included. A third category of coalescence is based on the sudden approximation of wave function projection, which was first developed for light nuclei production from nucleons [19][20][21][22] and then extensively applied to study charmed hadron and jet hadron production from constituent quarks [23][24][25][26][27][28][29][30][31]. In this work, we further develop this approach by including a complete set of s and p-wave hadronic states for a proper normalization of the total coalescence probability of zero-momentum charm quarks with a reasonable in-medium size of hadrons, which turns out to be important for understanding the observed Λ c /D 0 ratio. A rigorous 4-momentum conservation is introduced to our coalescence model which guarantees the boost invariance of the coalescence probability and the thermal limit of the produced hadron spectra. By combining our new hadronization approach with the advanced Langevin-hydrodynamics model that simulates elastic and inelastic energy loss of heavy quarks inside a realistic QGP medium, we obtain a satisfactory description of both the heavy flavor hadron R AA and the chemical compositions of charmed hadrons in relativistic heavy-ion collisions.
The coalescence model. -The momentum distribution of hadrons produced from quark coalescence is given by where p ′ h is the 3-momentum of the hadron, and p i is that of each constituent quark with i running from 1 to 2 (3) for meson (baryon) formation. The Wigner function W represents the probability for quarks to coalesce, and is calculated from the wavefunction overlap between the free quark states and the hadron bound state. Assuming a harmonic oscillator potential for a twobody system, we have the following Wigner functions for s and p-wave mesons [32,33]: Here, the position space component of the Wigner function has been averaged over the volume V , which will be cancelled by the volume factor associated with f i of light quarks; g h is the statistical factor taking into account of the spin-color degree of freedom; k is the relative momentum between the two constituent quarks defined in their center-of-momentum frame (the meson rest frame); σ = 1/ √ µω where µ is the reduced mass of the two-body system and ω is the oscillator frequency which can be directly related to the meson radius [23]. The thermal masses of light quarks at the hadronization temperature T c are taken as m u,d = 0.3 GeV and m s = 0.4 GeV, and the charm quark mass is taken as m c = 1.8 GeV, as suggested by the T -matrix approach that takes into account the mass contribution from the confining interaction [34,35]. This coalescence model can be extended to a three-body system by first combining two quarks and then combining their center-of-momentum with the third quark. We symmetrize all possible inner configurations for the formation of baryons, e.g., no matter which two of the three quarks are combined first, and no matter whether the first or the second step of combination carries the non-zero orbital angular momentum for p-wave baryon formation.
In this work, we incorporate all s and p-wave hadron states allowed by the spin-orbit coupling, which covers nearly all charmed hadron species reported in the PDG [36]. To compare to experimental data, we decay all D * 0 's to D 0 's, 68% D * + 's to D 0 's and 32% D * + 's to D + 's. For other excited states of D 0 and D + whose decay branching ratios are not clearly measured, we assume 50% probability to D 0 and 50% to D + . All excited states of D s are decayed to D s . All excited states of Λ c , together with all states of Σ c are decayed to Λ c .
In Fig. 1, we present the coalescence probabilities of charm quarks to different hadron species inside a static medium with temperature T = 160 MeV, as a function of the charm quark momentum p c . These probabilities are obtained via integrating the hadron spectrum as given by Eq. (1), where the charm quark distribution is taken as δ(p i − p c ), and the light quark distribution is taken as g i V /(e Ei/T + 1)/(2π) 3 , with E i and g i being the light quark energy and degree of freedom respectively. The contribution from thermal gluons is taken into account by converting a pair of thermal gluons (with mass m g = 0.3 GeV) into a quark-anti-quark pair gg → qq, equally distributed between u, d and s flavors [27,30]. The oscillator frequency is adjusted as ω = 0.24 GeV such that the total coalescence probability for a zero-momentum charm quark via all possible channels is summed up to unity when both s and p wave contributions are included, since charm quarks with zero momentum do not fragment into hadrons. This value of ω corresponds to an in-medium wavefunction radius of r = 3/(2µω) = 0.97 fm for D 0 , which is larger than its size in vacuum r = 0.83 fm when ω = 0.33 GeV is taken [23,37]. In principle, different oscillator frequencies can be applied to different hadron species. However, as a minimal model, we here apply a common parameter value to all charmed hadrons before theoretical calculations of their in-medium sizes become available. As shown in Fig. 1, the inclusion of the p-wave contributions can significantly enhance the total coalescence probability of charm quarks. Such enhancement is stronger for Λ c baryons than D 0 and D s mesons, thus the inclusion of p-wave contributions can also increase the baryon-to-meson (Λ c /D 0 ) ratio from our coalescence model. We would like to point out that our approach by tuning the ω value to normalize the coalescence probability is more sensible than just applying an overall normalization factor for the Wigner function [18,27] whose physical meaning seems vague. We also note that Refs. [23,27,30] only take the s-wave contribution and perform the calculations on a 2-dimensional transverse plane (for cylindrical medium).
To simulate the charm quark hadronization at the boundary of the QGP (T c = 160 MeV), we first boost each charm quark into the rest frame of the local fluid cell, and then determine the coalescence channel based on the momentum-dependent probabilities as shown in Fig. 1(b). For a selected channel, the charmed hadron momentum is then determined by its differential spectrum [Eq. (1)] before being boosted back to the global center-of-mass frame of heavy-ion collisions. Charm quarks that do not coalesce are converted to hadrons via Pythia using the Peterson fragmentation function [38]. A long-standing defect of the coalescence model is the lack of the energy conservation which breaks the boost invariance of the coalescence probability. In many earlier studies, the produced hadron with momentum p ′ h in Eq. (1) is forced to be on-shell. In our current model, we first assign the total energy of the combined quarks to the produced hadron (E ′ h = i E i ) and then decay it into an on-shell charmed hadron and a pion. The final energy and momentum (E h , p h ) of each charmed hadron is determined by the corresponding 2-body decay process based on its initial energy and momentum (E ′ h , p ′ h ) from Eq. (1). This setup guarantees the 4-momentum conservation and preserves the boost invariance. Figure 2 shows that our improved coalescence approach respects the thermal limit of hadronization. Here we use the Maxwell-Boltzmann distribution at T = 160 MeV for both charm and light quarks, and then compute the momentum distribution of D 0 via coalescence. As shown in Fig. 2, the D 0 spectrum from our coalescence model follows a Maxwell-Boltzmann distribution at T = 160 MeV. In contrast, although the old coalescence model without energy conservation seems to follow the thermal limit in the high energy region (above 6 GeV) where the hadron binding energy can be neglected, deviation appears at low energy. We have verified that the same conclusion also holds for D s and Λ c .
Nuclear modification of charmed hadrons. -To study the charmed hadron production and their nuclear modification in relativistic heavy-ion collisions, we employ the Langevin-hydrodynamics model developed in our previ- ous work [41,42] for heavy quark evolution inside the QGP prior to their hadronization via our improved hybrid fragmentation-coalescence approach.
The momentum distribution of heavy quarks produced in the initial hard scatterings is calculated using the leading order pQCD cross sections convoluted with the CTEQ5 parton distribution functions [43], which provides a reliable p T spectrum of D 0 as observed in protonproton collisions when combined with the Pythia fragmentation. For heavy-ion collisions, we use the EPS09 parameterization [44] to account for the cold nuclear matter effects. The spatial distribution of heavy quarks is initialized with the Monte-Carlo-Glauber model.
The space-time evolution profiles of the QGP medium, such as the local temperature and flow velocity, are simulated with a (2+1)-dimensional viscous hydrodynamic model [45][46][47]. The specific shear viscosity η/s = 0.08 is constrained by the soft hadron spectra observed at RHIC.
After the QGP evolution commences (τ 0 = 0.6 fm), the time evolution of charm quark momentum inside the thermal medium is described using the following modified Langevin equation [41,42]: The first two terms on the right represent the drag force and the thermal random force, inherited from the conventional Langevin equation, for the elastic energy loss of heavy quarks. The third term is introduced to describe the recoil force exerted on heavy quarks when they emit medium-induced gluons. It can be simulated using the differential gluon radiation spectrum which is adopted from the higher-twist energy loss calculation [48][49][50] in this work. As detailed in Refs. [41,42], this transport model has only one free parameter. By convention, we quote the position space diffusion coefficient D s . For each time step, we update the momentum of each heavy quark based on Eq. (4) in the local rest frame of the expanding medium, until the charm quark crosses the T c = 160 MeV hypersurface. Then heavy quarks are converted into hadrons using our new fragmentationcoalescence approach. In Fig. 3, we present the nuclear modification factor R AA 's of D 0 mesons and D-decayed electrons. The diffusion coefficient is set as D s (2πT ) = 3 to describe the high-p T data where R AA is dominated by the in-medium energy loss of heavy quarks and the uncertainties from the hadronization process is negligible. Our result shows that fragmentation alone is not sufficient to describe the non-monotonic p T dependence of the D-meson R AA at low p T . The inclusion of coalescence mechanism is necessary to yield the bump structure of R AA between 1-3 GeV observed in the experimental data. Note that the magnitude of the coalescence enhancement depends on the strength of the QGP radial flow, as we will see later.
Charmed hadron chemistry. -The hadron chemistry from the coalescence model is sensitive to the collective flow of the QGP medium. Also baryons experience stronger boost by the medium flow than mesons do. These can be clearly seen from Fig. 4. In the lower panels, we show D 0 and Λ c spectra from our full simulations of central Au-Au collisions at RHIC with the presence of the QGP flow. The upper panels show the simulation results when we switch off the fluid velocity during the hadronization process. One can see that the radial flow of the QGP can push the dominance of coalescence over fragmentation to higher p T region. With the presence of medium flow, the coalescence mechanism can dominate the D 0 formation up to 5 GeV, and the Λ c up to 8 GeV. This will directly affect the Λ c /D 0 ratio up to 8 GeV. Our result is consistent with the findings The pT-integrated Λc/D 0 ratio (a), the pT-differential Λc/D 0 ratio (b) and the pT-differential Ds/D 0 ratio (c), compared to the STAR data [51,52].
presented in Refs. [27,30] where a simplified parametrization of medium flow was implemented. We also note that the so-called space-momentum correlation suggested in Ref. [18] is automatically included in our current and previous studies, because the local temperature and flow information of the surrounding fluid naturally enters the hadronization process of heavy quarks in our model.
In Fig. 5, we present our result for Λ c /D 0 and D s /D 0 ratios. As shown in Fig. 5(a), the p T integrated ratio of Λ c /D 0 (from 3 to 6 GeV) increases with the participant number of heavy-ion collisions due to the stronger radial flow of the QGP in more central collisions. This ratio is significantly reduced if the medium flow is switched off during the hadronization process. The p T differential ratios of Λ c /D 0 is shown in Fig. 5(b). Again, the coalescence mechanism greatly enhances the Λ c /D 0 ratio as compared to applying fragmentation alone. We note that many theoretical calculations in literature compare their central collision results to the 10-80% centrality data. Our calculation shows a sizable difference between results in the two centrality bins. For a wide centrality region (e.g. 10-80%), we compute the hadron spectra via dN h /dp T = Σ c P (c) dN bin is the probability of finding heavy quark events in a given smaller centrality bin "c", with N (c) bin being the binary collision number within this bin. Note that if one uses only one average medium profile for the large centrality bin, much smaller nuclear modification would be obtained, i.e., a smaller Λ c /D 0 ratio here and also a larger R AA for the D-decayed electrons in Fig. 3(b). Figure 5(c) shows the D s /D 0 ratio; we see that this ra-tio is also significantly enhanced by the coalescence process, due to the combinational effects of the larger D s mass compared to D 0 and the enhanced strangeness production in a thermal medium compared to that in the vacuum fragmentation.
It is worth noting that Pythia fragmentation deviates from the measured Λ c /D 0 ratio in proton-proton collisions [53]. Several studies are devoted to understanding such discrepancy [15,54]. However, as shown in Fig. 4, the charmed hadron production within the p T regime measured at RHIC is dominated by the coalescence of charm quarks and the QGP partons, thus the inaccuracy of Pythia should have a minor impact on our charmed hadron chemistry presented in this work.
Summary. -We have developed a comprehensive fragmentation-coalescence approach for studying heavy quark hadronization in heavy-ion collisions. Both s and p-wave states have been included in our coalescence model; they are sufficient to cover nearly all charmed hadron states reported in PDG. We find that including p-waves increases the Λ c /D 0 ratio. It is also necessary for a full 3-dimensional calculation to normalize the charm-QGP coalescence probability at zero momentum using a reasonable in-medium size of hadrons. In our improved coalescence model, a strict 4-momentum conservation is implemented by first forming the off-shell excited hadron states and then decaying them into the ground states; this guarantees the boost invariance of the coalescence probability and the thermal limit of the produced hadrons. By combining our new hadronization approach with the advanced Langevin-hydrodynamics model, we have obtained satisfactory descriptions of the D meson R AA , D-decayed electron R AA , as well as the p T -integrated and differential Λ c /D 0 and D s /D 0 ratios. Our study shows that the coalescence mechanism and the QGP collective flow are both important for understanding the charmed hadron chemistry observed in relativistic heavy-ion collisions. It also suggests that the in-medium size of charmed hadrons should be larger than the size in vacuum, which may be tested by the hadronic model calculations in the future.