Hadronization of Heavy Quarks

Heavy-flavor hadrons produced in ultra-relativistic heavy-ion collisions are a sensitive probe for studying hadronization mechanisms of the quark-gluon-plasma. In this work, we survey how different transport models for the simulation of heavy-quark diffusion through a quark-gluon plasma in heavy-ion collisions implement hadronization and how this affects final-state observables. Utilizing the same input charm-quark distribution in all models at the hadronization transition, we find that the transverse-momentum dependence of the nuclear modification factor of various charm hadron species has significant sensitivity to the hadronization scheme. In addition, the charm-hadron elliptic flow exhibits a nontrivial dependence on the elliptic flow of the hadronizing partonic medium.


I. INTRODUCTION
At high temperatures and vanishing net-baryon density, lattice-discretized simulations of Quantum Chromodynamics (QCD) at finite temperature [1,2] predict a transition of hadronic matter to a new phase of matter, called quarkgluon plasma (QGP), which is a state of deconfined quarks and gluons.This transition has been realized experimentally via heavy-ion collisions at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC).It is widely accepted by now that, upon initial impact of the incoming nuclei, the system briefly evolves through a pre-equilibrium stage before a locally thermalized QGP is formed.Large pressure gradients in the fireball drive the QGP to expand and cool rapidly.Once the system reaches the crossover region with temperatures around the pseudo-critical temperature, T pc , the partons constituting the QGP undergo colorneutralization into hadronic bound states -a process generically denoted as hadronization.
Quantitative studies of the properties of hot QCD matter are a central objective of ultrarelativistic heavy-ion collisions, which, however, require suitable probes that can be observed in the final state while still carrying information about the hot and dense phase phases.A promising probe in this regard is heavy-flavor (HF) hadrons which have several advantageous features [3]: • The heavy-quark (HQ) mass is much larger than the nonperturbative QCD scale, m c , m b ≫ Λ QCD , and thus the production in nuclear collisions mainly occurs in initial hard scatterings, which can be reasonably well described by perturbative QCD; • The HQ mass is much larger than the typical temperature of the hot medium; this implies that their number is expected to be effectively conserved during the evolution due to high production thresholds and that their diffusion is akin to a Brownian motion which enables rather direct access to pertinent transport parameters.
• The determination of HQ fragmentation functions can be carried out at next-to-leading order in the production process within an HQ mass expansion using the methods of HQ effective theory (HQET) [4].
Several research groups have developed transport approaches that describe the dynamics of heavy quarks from their creation at the onset of a heavy-ion collision through their evolution in the QCD medium until their detection as heavy hadrons (or their decay products) in experiment [5][6][7][8][9][10][11][12][13][14][15][16][17].After about two decades of developments, these models are largely able to describe key experimental data, such as the scaled ratio between the transverse momentum (p T ) spectra in heavy-ion collisions and proton-proton collisions, R AA (p T ), and the azimuthal asymmetry in the p T spectra quantified by the elliptic-flow coefficient, v 2 (p T ), which have become increasingly precise in recent years [18][19][20][21][22][23].These transport approaches have become rather complex, aimed at a comprehensive description of the initial production of the heavy quarks (including nuclear effects in the initial state), their interactions with the QGP, hadronization, and the interactions of heavy hadrons in the hadronic phase.Most of these processes cannot be rigorously described using first-principles QCD as they occur in the non-perturbative regime, especially at low transverse momentum, p T .In addition, the transport calculations need to be embedded into realistic simulations of the rapidly expanding QCD medium, the dynamics of which vary in the different approaches.
The fair success of these models in describing HF measurements, given their differing implementation of the HQ dynamics and medium evolution, translates into significant uncertainties regarding the scientific conclusions that can be obtained from them.In order to gain a better understanding of the relevant commonalities and differences between these models a series of comparative studies have been performed in recent years: For example, in Ref. [24] benchmark HQ interactions within different medium evolutions and different hadronization schemes have been investigated, in Ref. [10] different descriptions of the interaction of heavy quarks with the QGP have been confronted, and in Ref. [25] the influence of different initial-production mechanisms of heavy quarks and expansion scenarios of the QGP have been analyzed.These works also contain discussions of the theoretical approaches that are being employed in the description of the HQ transport through the QGP medium.Significant commonalities were found between different approaches, thereby improving our understanding of HQ dynamics in a QGP medium.Initial studies of the different hadronization models were also carried out [24].With this article, we provide a dedicated study on the hadronization process in HF transport models by conducting in-depth comparative studies that have not been done before.
The paper is organized as follows.In Sec.II, we provide a general introduction to hadronization mechanisms for heavy quarks in relativistic heavy-ion collisions, followed by a detailed description of the implementation of these hadronization mechanisms into different dynamical models of HQ transport in Sec.III.In Sec.IV, we define a common environment of an expanding fireball that all models for HQ hadronization are subjected to.Section V provides the results and discussion of systematic comparisons of these models, facilitated by identical HQ distributions at the hadronization hyper-surface as a common input.Conclusions are given in Sec.VI.

II. HADRONIZATION MECHANISMS
In this section, we provide a brief overview of the hadronization schemes that will be encountered in the model implementations, specifically quark fragmentation typically encoded in universal fragmentation functions in Sec.II A and quark recombination in Sec.II B, specifically instantaneous coalescence models (ICMs) based on phase-space Wigner density (Sec.II B 1), and the resonance recombination model (RRM) which is carried out in momentum space (Sec.II B 2).

A. Fragmentation
Independent fragmentation of partons into hadrons is the standard way to describe hadronization in elementary collision systems, such as pp and e + e − .It is based on the assumption that the differential cross section for a hadron H with momentum P H factorizes into a hard production cross section for a high-energy parton (i) and the fragmentation fragmentation function, D i→H , as [26] E dσ H where E i dσi d 3 pi is the differential perturbative cross section for a parton with momentum p i , and z = P H /p i is the momentum fraction carried by the hadron.The symbol ⊗ denotes a generic convolution, including the relevant fragmentation fraction z Jacobian.The fragmentation function, D i→H , is a non-perturbative quantity but it is considered to be universal and usually extracted from experiments such as e + e − collisions.There are various choices of fragmentation functions see, e.g., the review in Ref. [27].The most frequently used fragmentation function for heavy quarks is from Peterson [28] in which the HF meson carries most of the momentum of the heavy quark Q, where ϵ = m 2 q /m 2 Q , the ratio of the light-and heavy-quark mass squared, is an adjustable parameter.It is commonly fixed by experimental data, e.g., in e + e − collisions.
Equation ( 1) is considered to be valid if p H ≫ Λ QCD , i.e., for high-p T D-mesons.Another popular choice is the HQET-based fragmentation function.For the production of a pseudoscalar meson, H P , or a vector meson, H V from a heavy quark Q, the fragmentation function can be expressed as [4,29] Here, the adjustable parameter is r = (m H − m Q )/m H , which can be interpreted as the ratio of the constituent mass of the light quark to the meson mass.r can also be fixed by experimental data.
Finally, some groups directly use the fragmentation function provided in PYTHIA with a modified Lund string fragmentation function [30].It can be expressed as where r, a α , a β , and b are free parameters, and m T = T is the transverse energy.In PYTHIA, r = 1.32, a α = a β = 0.68, and b = 0.98 are the default values.This fragmentation function depends separately on the HQ transverse momentum.

B. Parton Recombination
Experimental results indicate that in the hot and dense medium, created in relativistic heavy-ion collisions, an additional hadronization mechanism is at work.This was motivated by the enhancement of baryon-to-meson ratios, such as p/π and Λ/K ratio, observed in heavy-ion collisions at both RHIC [31][32][33] and the LHC [34,35], and an approximate quark number scaling of the elliptic flow [18,[36][37][38].The parton recombination model was proposed and successfully used to describe these phenomena [39][40][41][42], primarily in the regime of intermediate p T ≲ 5 GeV (the high-p T region is dominated by fragmentation).The most popular phenomenological model for the recombination process is the coalescence model, which assumes that all quarks hadronize at an equaltemperature hypersurface and that the coalescence probability is given by a Wigner function [41,43].Given that HF observables exhibit similar features as the aforementioned lightquark observables (e.g., the Λ c /D 0 ratio exhibits an enhancement [35] akin to the p/π or Λ/K ratio), it is natural to presume that heavy quarks also hadronize via parton recombination in the low-and intermediate-p T domains.
In the following, we briefly recollect the main features of ICMs (Sec.II B 1) and RRM (Sec.II B 2).

Instantaneous Coalescence (ICM)
In the sudden coalescence models, the momentum spectrum of the produced hadron H can be written as where g H is a degeneracy factor of color and spin; P and p i are the momenta of HF hadron and constituents quarks, respectively.The δ-function ensures the conservation of 3-momentum, dσ i is a hypersurface element, and F (x 1 ...x N , p 1 ...p N ) represents the constituent-anti-/quark phase-space distribution functions.It amounts to a product of single-particle distribution functions when neglecting correlations between the constituents.The Wigner density of the hadron, W (x...x N , p 1 ...p N ), is used as the recombination probability and can be constructed from the hadron wave function.For the two-body case, it can be defined as where r and p r are the relative distance and momentum in the center-of-mass (CM) frame of the hadron, respectively.The wave function ψ can, e.g., be obtained by solving the twobody Schrödinger equation.In practice, the Wigner density is not used directly.Rather it is assumed that for a hadron in an S-wave state, it can be approximated by a Gaussian Wigner density (which can be obtained analytically when using, e.g., the spherical harmonic oscillator) which gives the same rootmean-square (rms) radius as the original Wigner density, The parameter σ in the Wigner function is related to the rms radius of the hadron, σ 2 = 2⟨r 2 ⟩/3.Some authors relate the width to the rms charge radius [44,45] instead of mass radius (see, e.g., in Eq. ( 19) below).For baryons, the three-body problem can be reduced into two two-body problems by separating the CM motion with the help of Jacobian transformations, for details see Ref. [46].There are two relative coordinates, ρ and λ as shown in Fig. 1.The relative wave function of the three-body system in the ground state is assumed to be factorized as Ψ = ψ(ρ)ψ(λ).Consequently, the Wigner density can be expressed as a product of two Gaussians, where p ρ and p λ are the relative momentum corresponding to the relative coordinates, ρ and λ.The widths σ ρ and σ λ are related to the inner radii via The rms radius can be obtained by solving the three-body Schrödinger equation.We note that for comparisons to experimental data, one needs to consider excited states as they contribute through decay feeddown to the ground states.For radial and angular excited states, the Wigner function will become complicated but is still manageable by utilizing the wave function of the 3-D isotropic harmonic oscillator.
The probability density, dN H (P, p Q ), to produce a heavy hadron, H, with momentum P from a heavy quark, Q, with momentum p Q is defined by introducing a delta distribution centered at p Q into F in Eq. ( 5).One then defines the total recombination probability for this heavy quark to hadronize through coalescence as While this quantity is boost invariant, it is not the case for the spectra defined at Eq. ( 5), due to the violation of the energymomentum conservation inherent to the coalescence process.
A particular frame thus needs to be chosen, usually taken as the local fluid rest frame.

Resonance Recombination (RRM)
The resonance recombination model (RRM) has been derived from a Boltzmann equation where a recombination rate is obtained by utilizing a resonant quark anti-quark cross section (or pertinent scattering matrix element) for hadron fomration [47], instead of the Wigner density.In the long time limit, one can formulate this as a coalescence probability while still using off-equilibrium anti-/quark distributions akin to a linearresponse theory.In a similar spirit, a different recombination criterion has recently been proposed [17] for HF hadron production in heavy-ion collisions.The main difference between these models is whether the recombination is carried out in momentum space or utilizes phase-space criteria.In particular, the RRM obeys 4-momentum conservation and satisfies the equilibrium limit, i.e., it produces equilibrium hadron distributions if equilibrated quark distributions are used as input (including collective-flow effects).It has been implemented on a hydrodynamic hypersurface in Ref. [6] and extended to HF baryons in Ref. [15].
Let us briefly recall the main steps in deriving the RRM.For mesons, one starts from a distribution function, f M , that satisfies the Boltzmann equation, where R and P are the position and momentum of the meson, γ = E M /m is a Lorentz factor with meson mass m and energy E M , and Γ is the width of the meson.The first and second terms on the right-hand side are loss and gain terms.
In equilibrium, the collisional term disappears, yielding The integrated meson yield can then be expressed as where f q and f q are the phase space distributions of quark and antiquark, and p r = p 1 − p 2 and P = p 1 + p 2 are their relative and total momentum, and v rel is the relative velocity between quark and antiquark.In its current implementation, the production cross section is assumed to be of a Breit-Wigner form, where g σ = g M /(g q g q ) is the spin-color degeneracy, and k denotes the 3-momentum of the quarks (heavy or light) in the CM frame.s is the invariant mass.For baryons, a two-step process has been employed [15], where the first two quarks form a bound diquark, followed by its recombination with another quark into a baryon.

III. MODEL DESCRIPTIONS
While the basic features of hadronization via fragmentation or recombination are well-defined and understood, different models vary in their implementation of the respective hadronization schemes.In this section, we will provide detailed descriptions of the hadronization schemes used in each model setup.For HQ initial production and energy loss in these models, we refer to the original literature [5, 6, 8, 11-17, 24, 48].Here we limit ourselves to the comparison of those quantities that directly influence the hadronization scheme.These may include: 1.The functional form of the recombination probability; 2. Parameters such as the relative distance and relative momentum in the recombination probability; 3. The distribution of light quarks which coalesce with the heavy quark; 4. The quark masses and widths in the Wigner function; 5. The number of resonance states accounted for; 6.The choice of fragmentation function.
We need to define and distinguish different Lorentz frames used in the following discussion: a) The CM frame of two colliding nuclei, which is also called the laboratory (Lab) frame; b) The fluid rest frame, is usually based on a hydrodynamic description of the expansion of the hot medium; different regions of the medium have different fluid velocities.The fluid rest frame is the stationary frame inside the flowing fluid cell; c) The CM frame of the light and heavy quark or the diquark and the heavy quark (in the case of baryons).
As decay feeddowns from excited charm hadrons play an important role in the abundance (and spectra) of the various ground states (in particular, for D 0 , D s , and Λ c ) in the following analyses, we have included a table summarizing the excited states accounted for by each model in Tab.I.
In the remainder of this section, we discuss the hadronization models by the following groups: Catania (Sec.III A), Duke (Sec.III B), LBT (Sec.III C), Nantes (Sec.III D), PHSD (Sec.III E), TAMU (Sec.III F), Torino (Sec.III G), LANL (Sec.III H), as well as a summary table of the number of excited states involved in each model, and then a discussion of the various fragmentation schemes employed (Sec.III I).

A. CATANIA
The Catania model adopts a coalescence plus fragmentation scheme for HF hadronization.For the coalescence part, a Gaussian shape in space and momentum for the Wigner function is used (including only S-wave states), where N q is the number of constituent quarks and A W is a normalization constant that has been fixed to guarantee that in the limit of vanishing charm-quark momentum in the fluid rest frame, p c → 0, all charm quarks hadronize by coalescence into a heavy hadron.This is imposed by requiring that the total coalescence probability gives lim p→0 P tot = 1.Going into the CM frame of the particles involved in the process the relative coordinates are defined as follows.For mesons, they are given by while for baryons they are defined as for the diquark, and r 2 , p r2 for the quark-diquark are given by as shown in Fig. 1.
The width in the Wigner function is linked to the size of the hadron and, in particular, to the rms charge radius of the hadron, 3 for mesons and baryons, respectively; R cm is the CM coordinate.For mesons one has and for baryons where Q i is the charge of the i-th quark.The charge radius of the hadrons is taken according to the constituent-quark model.In the Catania model one has 475GeV −1 GeV for D s mesons.For Λ c , the rms charge radius is taken as ⟨r 2 ⟩ ch = 0.15fm 2 .Actually, for baryons, there is only one free parameter, because the two widths are related by the oscillator frequency ω through the reduced masses, µ i , via where µ 1 is the reduced mass of two-body system and µ 2 is the reduced mass of the two-body system with the third particle.The corresponding widths are σ 1 = 5.556GeV −1 and σ 2 = 2.924GeV −1 .
In the Catania model, all first excited states for D mesons and Λ c are included.For the resonances, the coalescence probability is multiplied by a suppression factor that takes into account the Boltzmann probability to populate an excited state of energy E + ∆E, at temperature T and a statistical factor, , where m R is the resonance mass and m G is the mass of the ground state.
If a heavy quark hadronizes via the coalescence process, one or two light quarks will be sampled on the hadronization hypersurface.The light quarks are assumed to be thermally distributed.Their momentum distribution in the lab frame satisfies where m T = p 2 T + m 2 q , and β T is determined by the velocity field given by the bulk.The factor g = 6 is the spin-color degeneracy.The presence of gluons in the QGP is taken into account by converting them to quarks and anti-quark pairs according to the flavor compositions; τ is the proper time of the hadronization hypersurface.A uniform distribution in coordinate space is assumed.In this model comparison, the quark masses are taken as m u,d = 0.3 GeV, m s = 0.38 GeV, and m c = 1.5 GeV.
Heavy quarks that do not hadronize via coalescence are converted to hadrons by fragmentation.Catania uses the Peterson fragmentation with ϵ = 0.1 for charmed mesons and ϵ = 0.02 for charmed baryons.
For more details, we refer the reader to Refs.[7,48].

B. Duke
The hadronization model used in the Duke model is also a combination of coalescence and fragmentation.The coalescence probability is given by the momentum-space Wigner function, and for the results shown here only S-wave mesons are taken into account.The Wigner function is taken as a Gaussian, where p r is the relative momentum defined in the CM frame, The width in the Gaussian distribution can be obtained by σ = 1/ √ µω with µ being the reduced mass of the light and heavy quark.Duke takes ω = 0.106 GeV for all charmed hadrons.So, the related widths are, σ(D) = 6.23 GeV −1 and σ(D s ) = 5.22 GeV −1 .Both ground states and first excited states of D mesons are considered.However, the hadronization implementation utilized for this work [49] does not include baryonic states.Later implementations of the Duke hadronization model [8] do include baryons such as the Λ c , Σ c , and Ω c .The light-quark distribution in the fluid rest frame satisfies where g i = 6 is the statistical factor.The quark masses are taken as m u,d = 0.3 GeV, m s = 0.475 GeV, and m c = 1.27GeV.For fragmentation, Duke uses the string fragmentation in PYTHIA 6.4.The main references for Duke's recombination plus fragmentation approach to HQ hadronization are [8,49].

C. LBT
Hadronization applied in the LBT model is also described by coalescence plus fragmentation.The coalescence probability is given by the momentum space Wigner function, which can be expressed for S-wave and P -wave as where the spatial part of the Wigner function has been integrated over coordinate space.The integral will be canceled by the volume factor associated with the momentum distribution functions of light quarks; g h is the statistical factor; p r is the relative momentum between the two constituent quarks in the CM frame, where E 1 and E 2 are the energies of the quark and antiquark.Baryons are treated as two two-body systems (produced by recombining two quarks first and then using their CM to recombine with the third one).All S and P -wave hadron states, allowed by the spin-orbit coupling in a full 3-dimensional calculation, are included, which covers nearly all charmed hadron species (Ξ c and Ω c are also included) as reported by the particle data group (PDG).
The width of the Wigner function is determined by σ = 1/ √ µω with µ the reduced mass and ω the oscillator frequency, ω = 0.24 GeV, for all charmed hadrons (both S and P -wave states), which is tuned so that the total coalescence probability for a zero-momentum charm quark (which cannot fragment) is equal to one when all Sand P -wave charmed mesons and baryons are included in the calculation.The quark masses are taken as m u,d = 0.3 GeV, m s = 0.4 GeV, m c = 1.8 GeV, and the related widths are σ(D) = 4.02 GeV −1 and σ(D s ) = 3.57GeV −1 .For Λ c , there are two widths, σ 1 = 4.93 GeV −1 and σ 2 = 2.87 GeV −1 , corresponding to the diquark system and quark-diquark system, respectively.Light quarks in the hot medium are thermalized.Their momentum distribution in the fluid rest frame satisfies where g i = 6 is the statistical factor.
The traditional instantaneous coalescence model only respects the 3-momentum conservation but violates the energy conservation.A solution of this problem is proposed in Ref. [11] by allowing quarks to combine into an off-shell state of a charmed hadron first and then decay into its on-shell state by emitting a pion.The pion and the final state charmed hadron are generated back-to-back in the rest frame of the offshell charmed hadron, and then boosted back into the global frame using the velocity of this off-shell state.Since sufficiently large masses are applied for thermal quarks, the kinematics of such decay is usually allowed.In rare cases when the decay is forbidden by kinematics, a photon is emitted instead of a pion.
Heavy quarks that do not hadronize via coalescence are converted to hadrons by fragmentation.The fragmentation process used in the LBT model is PYTHIA 6.4 and the default Peterson fragmentation is used with ϵ = 0.05 for charmed hadrons.
For more details on the hadronization model applied to the LBT calculation, please see Ref. [11].

D. Nantes
In the Nantes model, heavy quarks hadronize either by coalescence or by fragmentation.The recombination probability is evaluated in the thermal restframe and normalized to one for vanishing HQ momentum in that frame.Thus, in the lab frame, it depends on the fluid cell velocity and the orientation of the hadronization hypersurface.The coalescence mechanism is based on the model of Dover [50].
The coalescence probability is given by where u Q is the four-velocity of the heavy quark and with α d = 0.51; R c is fixed for a given a d from the normalisation condition In coordinate space, one can obtain σ r = √ 2R c , but in momentum space one has σ p = 1/( √ 2α d µ) by definition of the F Φ , with µ as the reduced mass.The widths, σ, in coordinate and momentum space are therefore not related by the uncertainty relation.The Nantes model is restricted to the hadronization of c (resp.b) into D (resp.B) mesons and only the ground states are accounted for.
Light quarks are assumed to be thermalized.Their distribution is given by the Boltzmann-Jüttner distribution in the fluid rest frame, GeV.Heavy quarks that do not coalesce form mesons by fragmentation.The fragmentation function used in the Nantes model is the HQET-based fragmentation function [4,29].
For more details, we refer the reader to Ref. [12].

E. PHSD
PHSD is a kinetic approach and therefore not based on a hydrodynamical expansion of the QGP.To make it comparable to the other models it has to be modified.For the comparison, the calculation is carried out in a box of size 10 × 10 × 10 fm 3 .The heavy quark is placed at the center of the box and antiquarks are uniformly distributed in the box with a random momentum given by the thermal distribution at T = 0.18 GeV and with a random mass from the spectral function of antiquark within the dynamical quasi-particle model (DQPM) [51,52], that is, in the fluid rest frame, with g = 6 the color-spin degeneracy factor; M and γ are the pole mass and spectral width of antiquark, respectively.The coalescence probabilities for the S-wave and P -wave states are given by W P (r, p r ) = 2s + 1 36 16 3 where s is the spin of meson.The relative distance and momentum in the CM frame are, respectively, The width σ is related to the rms radius of the D meson through for the S-wave state and for the P -wave state, and ⟨r 2 M ⟩ = 0.9 fm for both Swave and P -wave states.Different from the charm quark with m c = 1.5 GeV, the light quarks have no fixed mass but are off-shell and sampled with the given distribution of Eq. ( 33) [51,52].The relative momentum coordinates, Eq. ( 35), are defined as in the Catania model, Eq. ( 16), while for the Duke and LBT, Eq.s (24) and ( 27) it appears as a different definition in terms of the energy and not the masses of the recombining quarks; however, we notice that in the center of mass of the pair, the two definitions correspond to the same quantity, being p 1 = −p 2 .
Only Sand P -wave mesons are taken into account for coalescence (not taking into account the principle quantum number), while charm baryons such as Λ c are produced only through fragmentation.After production, a P -wave state immediately decays into an S-wave state by emitting a pion.The mass of the P state is taken as 2.46 GeV for the u and d sectors and as 2.57 GeV for the s sector.Since the probability is low, the coalescence process is repeated many times until the probability becomes large enough at low p T .It turns out that about ten trials make the coalescence probability similar to those in heavy-ion collisions at RHIC and LHC energies [13,14,53].
If the heavy quarks do not hadronize via the coalescence process, they will convert into charmed mesons or Λ c via fragmentation.The Peterson fragmentation function with ϵ = 0.01 for charmed mesons is used in the PHSD model.
For more details, we refer the reader to Refs.[13,14].

F. TAMU
The TAMU model applies a combined recombination and fragmentation approach as well.The former is realized via the resonance recombination model (RRM) [47], where the recombination probability for the two-body case is controlled by resonance amplitudes and is currently expressed as a relativistic Breit-Wigner cross section, where g σ = g M /(g q g q ) is a statistical weight, k denotes the quark 3-momentum in the CM frame, and m and Γ are the mass and resonance width of the meson, respectively.The recombination probability for a charm quark with momentum p * c in the fluid rest frame to form a meson can be expressed as and for baryons as where f q is the thermal light-quark distribution, with a degeneracy factor g i = 6; v 12 rel and v dqc rel are the relative velocities between the two light quarks, which form a diquark, and between the diquark and the charm quark, respectively (recall that ).The light-and heavy-quark masses are taken as m u,d = 0.3 GeV, m s = 0.4 GeV and m c = 1.5 GeV, and the diquark mass as m ud = 0.7 GeV.The resonance widths for meson, diquark, and baryon are taken as Γ M = 0.1 GeV, Γ dq = 0.2 GeV, and Γ B = 0.3 GeV; γ M , γ dq , and γ B denoted the Lorentz-gamma factors of the meson, diquark, and baryon states, respectively.An overall parameter P 0 is introduced to adjust the recombination probability equal to 1 for a charm quark at rest in the thermal frame.
The RRM in the TAMU model includes a large set of excited states of charmed hadrons that include all states listed by the PDG [54], but go well beyond that in terms of "missing" charm-baryon states as predicted by the relativistic-quark model (RQM) [55] and lattice QCD [56,57].The excited Λ c states and Σ c states are assumed to feed the ground state Λ c with a 100% branching fraction.
Different from the most recent version of RRM [15], where space-momentum correlations (SMCs) between the diffusing quarks and the phase space distributions of the thermal quarks in the hydrodynamically expanding medium are incorporated, the SMCs are neglected for the present study.This means that the RRM is carried out in momentum space only.
At the end of the diffusion through the QGP, when the charm quark with momentum p c T in the lab frame has entered a fluid cell that hadronizes, its momentum is boosted into the fluid rest frame, where it has a momentum p * c .The recombination probability, P rec (p * c ), is then calculated via Eqs.( 39) and (40), i.e., for a p c T we obtain the corresponding recombination probability P rec (p * c ).Finally, the recombination probability P rec (p c T ) in the lab frame can be obtained by dividing the sum of probabilities P rec (p c T ) by the number of incoming particles with p c T .Heavy quarks that do not hadronize via RRM are fragmented into HF hadrons using HQET-based fragmentation functions.The charm-quark fragmentation probability is identified as 1−P rec (p c T ), which is partitioned over all D, D s , and Λ c baryons according to their respective fractions determined in pp collisions [54].

G. Torino
The hadronization model developed by the Torino group is based on a local color neutralization mechanism: once a heavy quark reaches the hadronization hypersurface it undergoes recombination with the closest opposite color charge, forming an excited color-singlet object which will eventually decay into a final charm hadron plus a soft particle allowing for exact four-momentum conservation.The mechanism works as follows.First, one defines a hadronization hypersurface, identified by the temperature T H = 155 MeV.Once a charm quark reaches a fluid cell belonging to such a hypersurface element it is recombined with a thermal particle from the same cell.This thermal particle can be either a (anti-)quark or a (anti-)diquark.Light quarks and diquarks are assumed to be thermally distributed in the local rest frame of the fluid cell.The algorithm consists of the following steps: 1. Randomly extract the medium particle involved in the recombination process, with a statistical weight given by its corresponding thermal density at the temperature T H , namely where the spin degeneracy g s and isospin degeneracy g I are factorized; M is the quark (diquark) mass taken from the default values of PYTHIA 6.4 (m u/d = 0.33 GeV, m s = 0.5 GeV, m (ud)0 = 0.579 GeV, m (sl)0 = 0.804 GeV, m (ss)1 = 1.0936GeV,...) and the −/+ sign refers to fermions/bosons, respectively.2. Once the particle species is selected, determine its threemomentum in the local rest frame of the fluid cell from a thermal distribution.
3. Boost the four-momentum of the medium particle to the laboratory frame and recombine it with the heavy quark, constructing the cluster C; 4. Calculate the invariant mass M C of the cluster.If the latter is smaller than the mass of the lightest charmed hadron in that given flavor channel, re-sample the medium particle.If not, after defining an intermediate cutoff M max set to M max = 3.8 GeV, there are two possibilities: (1) If M C < M max the cluster decays isotropically in its own rest frame into two hadrons, which are then boosted back to the laboratory frame.The daughter-charmed hadron will carry the baryon number and strangeness of the parent cluster.(2) If M C > M max the cluster hadronizes via string-fragmentation simulated through PYTHIA 6.4.
Only ground-state hadrons, like D 0 , D + , D s , Λ c , Ξ c and Ω c , are populated in the model.The resonance decay is roughly encoded in the decay of the parent clusters, following the PDG branching ratios for resonances of mass around M C whenever possible.
For more details, we refer the reader to Refs.[17,58], where the model is also applied to charmed-hadron production in pp collisions.

H. LANL
The LANL model only includes fragmentation processes.In pp collisions, heavy quarks fragment into HF hadrons which is described by the HQET-based fragmentation function [4,29].
However, in high-energy heavy-ion collisions, such as the ones at the LHC, the gluon contribution to heavy-meson production can be as large as ∼ 50% [59].Contributions to parton energy loss and in-medium parton showers in general, to hadron production, can be accounted for in different but related ways [60].In the LANL model an effective modification of the fragmentation function [61,62] is employed as where P (ϵ) is the probability that the heavy quark loses a fraction of its energy, ϵ, due to multiple gluon emissions, and dN g /dϵ is the distribution of the average gluons as a function of ϵ = ω/E; D H/c is the fragmentation function of the charm quark into charmed hadrons, whereas D H/g is the fragmentation function of gluons into charm hadrons.For more details, please see Refs.[16,59].

I. Fragmentation Function Comparison
Different charm-quark fragmentation functions as used in the approaches described above are compared in Fig. 2. The Catania (light blue), LBT (green), and PHSD (dark blue) groups use the Peterson fragmentation function [28] with the parameter ϵ = 0.1, 0.05, 0.01, respectively.The Nantes (orange), TAMU (orange), and LANL (purple) groups use HQET-based fragmentation functions.The TAMU model uses one parameter for the pseudoscalar meson channel, r = 0.1 for the D 0 meson, while resonance states are obtained by the scaling law r M /r D 0 = ((m M − m c )/m M )/((m D 0 − m c )/m D 0 ).In the Nantes model, the charm quark is randomly chosen to fragment via the pseudoscalar or the vector channel with a branching ratio of 0.418.The r parameter is set to r = 0.1 as well.The LANL model uses fragmentation into scalar and vector charm mesons with the r parameter set to r = 0.2.The Duke and Torino groups directly utilize PYTHIA 6.4 for fragmentation.
We observe that the fragmentation functions differ considerably in the different approaches.Four of them (Nantes, PHSD, TAMU, and Torino) are rather narrowly peaked around 0.9 and therefore transfer a large momentum fraction to the D mesons, whereas three others (LANL, LBT, and Catania) show a rather broad distribution with a significant probability to produce D mesons which carry half or less of the momentum of the c quark.The ordering is not so much governed by the formalism but rather by the parameters used in each formalism.

Catania
Duke LANL we set pT = 3 GeV and mc = 1.5 GeV for the plot.Here the solid lines in the LANL and Nantes models are pseudoscalar channels, and dashed lines are vector channels.

IV. SETUP OF THE MODEL COMPARISON
To scrutinize the hadronization mechanisms realized in different evolution frameworks, we embed the different model implementations into an otherwise identical environment.In the present study, the expanding QGP is modeled by a simple fireball model [5].The fireball is cylindrically oriented along the z-direction with an elliptic cross sectional area.The long (a) and short (b) semi-axes of the ellipse are oriented, respec- tively, along yand x-direction and are given by The velocity of the fireball surface along the long and short axis is given by v a (t) = ȧ(t) and v b (t) = ḃ(t).Hadronization is assumed to occur on a hypersurface R µ = (t, R), which is a 4-dimensional surface with equal temperature T (R µ ) = T H .This corresponds to a hadronization time τ = τ (T c , R) in the above set-up, which amounts to t =4.86 fm/c at a final temperature of T H = 180 MeV.This defines the coordinates and velocities of the fluid cell at the hadronization hypersurface.The chosen temperature value, T H , is somewhat larger than the pseudo-critical temperature of the chiral transition (around T pc ≃ 155 MeV) [56], but the latter does not have to coincide with the (onset) of hadronization, which is expected to be a continuous process over a finite temperature interval.For the comparative study carried out here, this is of minor importance, where we rather need a fireball with a simple and well-defined flow field (and do not aim at quantitative comparisons to data).
The light quarks are assumed to be thermalized in almost all models.By taking a thermal distribution for u, d, and s quarks, their three momenta can be sampled in the fluid rest frame.Due to the different velocities in the xand y-directions, the light quarks boosted to the lab frame will carry a non-zero elliptic flow, which is defined as the second-order coefficient of the Fourier expansion of the azimuthal distribution.The resulting transverse-momentum distribution and the elliptic flow of light quarks are displayed in Fig. 3. Charm quarks are also assumed to be uniformly distributed on the hadronization hypersurface, and their momentum dis-tribution in the lab frame is taken from the EMMI rapid reaction task framework [24], displayed in the top panel of Fig. 4. The distribution can be parameterized as dN/dp T = 53.15(pT /GeV)/(4.423+ (p T /GeV) 2 ) 2.808 .The corresponding charm-quark elliptic flow, which is shown in the bottom panel of Fig. 4, will be specified below.
In realistic simulations of heavy-ion reactions, charm quarks move outward through a hypersurface element and typically possess a velocity that is correlated with the one of the local fluid environment, resulting in space-momentum correlations.In this study, we neglect these correlations: charm quarks are sampled isotropically in each fluid cell at the hypersurface.(c) the full framework of fragmentation plus coalescence, as realized in the respective model.
For each case, provide a "hadronization ratio" [24] H where D includes D 0 and D + .D s and Λ c are the prompt yields (direct + strong decay).
2. The elliptic flow v 2 of all charmed hadrons with an anisotropic charm-quark spectrum according to with v 2 (p T ) = 0.0577p T /(1 + 0.4212p 1.6367 T ), which is shown as the black dashed line in Fig. 4 (p T in units of GeV).
The calculation should consider two cases: (a) only-fragmentation (assuming all c quarks to hadronize through fragmentation); (b) the full hadronization scheme used by the respective model (recombination + fragmentation).
3. The H AA and v 2 of directly produced D 0 , D s , and Λ c (without feeddown contributions).
4. dN (D 0 )/dp T of direct D 0 mesons (without feeddown contributions) produced by a c-quark with p c T = 3 and 10 GeV.
For the subsequent discussion we define the yield for a given charm hadron with transverse momentum, P T , originating from a c-quark with momentum p c T as where dN H encodes the physics of the recombination.This suggests to organize a model comparison in terms of the 3 ingredients P rec , dN H and D c→H , which is part of the motivation for the tasks to provide not only the "mixed" (or total) hadronization results but also "only-fragmentation" and "only-recombination" cases.

V. RESULTS OF THE MODEL COMPARISONS
The comparisons laid out in the previous section are carried out in this section and organized as follows.In Sec.V A we first collect the inclusive production fractions of the groundstate hadrons D, D s , and Λ c , for the cases where only fragmentation or only recombination is accounted for, and for the combined scenario in each approach, followed by a study of the p T -differential recombination probabilities.In Sec.V B  we scrutinize the p T dependence of charm-hadron production using the hadronization ratio, H AA [24] and variants thereof, as well as ratios of different hadron species (hadro-chemistry).In Sec.V C we analyze how the p T -dependent elliptic flow of the input charm-quark distribution is converted into that of various hadrons.

A. Charm-Hadron Production Fractions
We start our analysis by listing the p T -integrated production fractions of the ground-state hadrons in Tab.II, i.e., Dmesons as well as those carrying an additional conserved quantum number, i.e., the charm-strange D s meson and the Λ c baryon.We do that for 3 scenarios where all charm quarks are forced to hadronize via either fragmentation or recombination, or in the "realistic" scenario (i.e., as used in practice) of both mechanisms contributing.Not all entries are filled, especially for the "recombination-only" case, as it is not straightforward to enforce 100% recombination in several models.Also note that even in the mixed scenario the fractions do not necessarily add up to one, indicating that other charm hadrons with additional conserved quantum numbers contribute (such as charmstrange or doubly-charm baryons).
Next, we turn to the recombination probability, P rec (p c T ), as a function of charm-quark transverse momentum, p c T , shown in Fig. 5.This represents the probability, calculated in the local thermal rest frame for most models (but displayed as a function of the charm-lab momentum, p c T ) that a charm quark hadronizes via recombination.In Fig. 5, top, P rec includes hadronization to all charmed mesons (such as D 0 , D + , D s , and all excited states), while in the bottom panel we also include charmed baryons (such as Λ c , Ξ c , Ω c , and all excited states).The color coding is the same as that in Fig. 2. The fragmentation probability P f rag (p c T ) is given by 1−P rec (p c T ).There is no quantitative guidance from QCD on how P rec should behave as a function of p c T .Only at large p c T , the domain of DGLAP jet evolution, hadronization is dominated by universal fragmentation functions and therefore P rec → 0. This lack of guidance is a main reason why the p c T dependence of P rec is rather different in the different approaches.In some approaches fragmentation also contributes at very small p c T , whereas in others hadronization at small p c T is exclusively due to recombination.The latter is based on the assumption that the probability goes to one when p c T → 0 (usually applied in the thermal rest frame), whereas in the former case the analytical recombination equations are simply extended to p c T = 0, which implies a significant dependence on the parameters in the respective coalescence model.One should also recall that most coalescence models are performed in the thermal rest frame, while others are directly applied in the lab frame which can cause issues with Lorentz invariance.In fact, also fragmentation functions, usually evaluated in the lab frame, are not guaranteed to be Lorentz invariant.While the p c T dependence for relative contributions from recombination and fragmentation covers a large variation in the models, we nevertheless observe 3 groups.The hardest dependence is found for the TAMU model, likely due to the large set of excited resonances that go out to rather large masses, thus facilitating charm-quark recombination at relatively large momenta and further augmented by the Lorentz boost due to the collective fireball flow.The next hardest distributions are from LBT and Catania; for the former, off-shell production (with subsequent pion decay) is presumably the dominant mechanism at high p c T (akin to the excited states in the TAMU model), while the Catania model features rather small wave function radii compared to the other ICMs (Duke, Nantes, and PHSD), which form the class with the softest p c T dependence and show rather good agreement.

B. Transverse-Momentum Distributions of Charm Hadrons
We first carry out task no. 4 by studying how the transversemomentum (p D T ) distribution of D-mesons looks like in the different models if one hadronizes a c-quark of a given transverse momentum, p c T .If created by a fragmentation process, the D-meson carries only a fraction of the p c T of the c-quark, as shown in Fig. 2. For recombination one generally expects the opposite trend, as the momentum of the D-meson is the (vector) sum of the HQ and light-quark momenta, and most algorithms favor an alignment of the light-quark p T with p c T , especially for ground-state hadrons.The distributions of p D T for a charm quark with p c T = 3 and 10 GeV are displayed in Fig. 6 left for the case that all D-mesons are created via recombination.To minimize the influence of parameter choices, we have decided to evaluate quantities pertaining to the directly produced hadrons -hence Figs. 6, 10, and 12 -requiring the charm-and light-quark masses in the calculations to be m c = 1.5 GeV, m u/d = 0.3 GeV, m s = 0.4 GeV.In addition, we also impose σ = 0.5 fm (charmed mesons) and σ ρ = σ λ = 0.5 fm (charmed baryons) for all models that use the Wigner function as recombination probability.
For recombination, almost all models show a slight shift of the maximum towards a larger D-meson momentum, p D T , as compared to p c T , with a width of the order of 1 GeV for p c T = 3 GeV, which increases with p c T , cf. right panel of Fig. 6.Only in the LBT model, the produced charm hadron, which is usually off-shell at first and subsequently decays into a pion and an on-shell charm hadron, one has p D T < p c T + p q T .If we weight fragmentation and recombination according to Fig. 5, the form and the maximum of the distributions become rather different, even for p c T as low as 3 GeV, see the right panel of Fig. 6.For p c T = 10 GeV most mesons are produced by fragmentation and therefore all models besides the TAMU model show a shift exclusively to lower D-meson p T (while the recombination portion in the TAMU model still generates a component above the c-quark momentum).For p c T = 3 GeV we see that contributions from both recombination and fragmentation are operative.Next, we turn to the H AA of D, D s , and Λ c production, displayed in the top, middle, and bottom panels, respectively, of Fig. 7. To exhibit the influence of the hadronization mechanisms we display the H AA for the case when all hadrons are produced exclusively via fragmentation (left), exclusively via recombination (middle), and in a combination of both (right); all results include feeddown contributions.As expected, fragmentation, where the heavy hadron has a smaller p T than the heavy parton, yields H AA values well below one for large p T , which are directly related to the pertinent fragmentation functions, recall Fig. 2: if the fragmentation function peaks at large values of z (like that of the Nantes group), the momentum change of the heavy hadron, relative the heavy quark, is smaller and therefore the H AA is larger.The uncertainties in the fragmentation functions cause an uncertainty in the highp T H AA of D mesons of almost a factor of 2, and similar for the D s .For the Λ c the models that use fragmentation functions from e + e − collisions are in fairly good agreement, while in the TAMU model the H AA is larger over the entire p T range by roughly a factor of ∼2-3, primarily due to the large number of excited resonances which have been constrained by a fit to pp data [54].At low p T fragmentation functions are not reliable, and therefore the extrapolation is rather uncontrolled unless constrained by pp data.Usually, this extrapolation does not affect the predicted spectra in heavy-ion collisions, since recombination dominates at low p T .
The middle panels of Fig. 7 show the results for models when enforcing hadronization by recombination only.In the Nantes model, which only considers D mesons, one finds the expected coalescence effect where the addition of a light quark pushes the hadron distribution out to higher p T , entailing a deficit at low p T .For a model that has more than one heavy flavor hadron, the trend of H AA for a specific state depends not only on the coalescence probability but also on the decay contributions.In the PHSD model resonance feeddown "repopulates" the lower p T region at the expense of higher p T where, consequently, the H AA of D decreases.This effect is much less pronounced in the LBT model, despite its implementation of off-shell production with subsequent pion decays.For D s , both the PHSD and LBT model show a decreasing H AA at large p T , due to the smaller coalescence probability for D s compared to D [11].A similar, albeit stronger, trend is observed for Λ c production.However, above p c T = 5 GeV, P rec is (well) below 0.1 for all models (except TAMU), and therefore the detailed form of H AA from recombination-only at high p T does not play a role for the total H AA , which is the combination of recombination and fragmentation weighted with P rec and 1 − P rec , respectively.
The total (or "mixed") H AA , which is a combination of the fragmentation and recombination (see Eq. ( 47)), is shown in the right panel of Fig. 7. Its high-p T behavior is almost identical to the fragmentation-only H AA .However, towards lower p T the increase of H AA is now due to the recombination contribution.At still lower p T , most models develop a levelingoff in the total H AA , or even a slight maximum structure in the Nantes model as a consequence of the stronger effect in the re-combination part (see middle panel).In the TAMU, LBT, and PHSD models, a "flow bump" is not generated for the given fireball configuration, pointing at differences in how the parton collectivity (both the charm quark and the collective flow of the medium) is converted into the one of hadrons.
Even if in some models a heavy quark with p c T = 0 hadronizes exclusively via recombination, D-mesons with p D T = 0 can also originate from fragmentation of c-quarks with a finite p c T .Therefore, at p D T = 0 we see a difference between the pure recombination and the mixed scenario.For D s mesons, a fair agreement of the H AA 's is found, except for the Catania model which falls short of the other models by a factor of 2-3, as was borne out already of the inclusive fractions listed in Tab.II.For Λ c baryons, one finds an appreciable sensitivity to the recombination kinematics at low p T , even though the inclusive fractions are not very different (between 15% and 24%, cf.Tab.II).At higher p T the TAMU model sticks out due to the large feeddown contribution from excited baryons.
The model differences in the charm hadro-chemistry are a substantial part of the discrepancies found in Fig. 7.In an attempt to eliminate this effect and focus on the kinematics, we display in Fig. 8 for coalescence-only and mixed-hadronization processes, and R = for fragmentation-only processes, where the integration limits are chosen under the premise that at high p T fragmentation functions can be reasonably well adjusted to experimental data.We observe that the fragmentation-only H r AA 's (left panels) at high p T now rather closely agree for p T ≳ 8 GeV among all models for all 3 hadron types.Even for at low p T ≲ 5 GeV, the agreement is not bad (with an outlier from LBT), which, however, may not be very significant in practice.For the recombination-only scenario (middle panels in Fig. 8), the model curves are now closer together, but the features in the p T dependence remain as discussed in the context of Fig. 7.
For the "mixed" H r AA , the D-meson results are now rather similar (in contrast to the H AA in Fig. 7), with values around 1 at low p T , where most of the normalization is operative.The shape differences are more pronounced for D s and Λ c , where two classes seem to emerge, i.e., Catania and Torino vs. PHSD, TAMU and LBT.Somewhat surprisingly, the high-p T H r AA for the mixed case shows a similar convergence between the models for all hadron species (with the exception of the Λ c from TAMU) as seen in the fragmentation only case (left column), despite a rather different normalization.
Next, we study the recombination and fragmentation probabilities for different charmed hadrons, by plotting in fragmentation-only (left panels) we see an almost constant ratio for all approaches, indicating that the form of the fragmentation functions are very similar.The absolute values differ, however, slightly for D s , while for Λ c the TAMU model gives a twice larger multiplicity than in all other models (again due to charm-baryon resonance feed-down adjusted to reproduce to pp data).For the recombination-only scenario, we only have two results for D s , which differ quite a bit.The LBT results show an increased ratio towards low p T , presumably a consequence again of the pion decay feeddown, whereas for PHSD a falling trend is found characteristic of recombination with little or no feeddown.For the mixed scenario (fragmentation plus recombination), with the weights plotted in Fig. 5, all D s /D results show an increase over the onlyfragmentation scenario, especially at low p T , and in all approaches this is likely due to the enhanced s-quark content in the QGP, most pronounced for the local color neutralization mechanism (Torino), followed by TAMU (including reso-nance feeddown), and less for Catania, LBT, and PHSD which use Wigner functions.For the Λ c /D ratio, the Torino and TAMU results are rather large again, while for Catania the tuning of the fragmentation function to pp data also translates into a large value in the heavy-ion environment, while the LBT result, with a fragmentation function tuned to e + e − collisions, remains rather small, similar to the e + e − results.The TAMU result becomes the largest at high p T due to the feeddown of heavy charm-baryon resonances present in their fragmentation function.
To eliminate the feeddown from excited states, which vary largely throughout the models, we compare in Fig. 10 the H AA 's for directly produced hadrons, to be confronted with the inclusive results in Fig. 7.The left column shows the results for recombination only whereas the right column contains the combined results for recombination plus fragmentation.From top to bottom, we plot the results for D, D s , and Λ c .For the two recombination-only results, LBT and PHSD, the H AA value for both D 0 and D s drops markedly compared to the prompt ones plotted in Fig. 7.This points to the importance of the excited states in both models.The decay of excited states shifts the heavy-meson momentum.From the near-constant H AA as a function of p T for directly produced heavy hadrons observed in PHSD, one may postulate that the decrease of H AA towards large p T observed for prompt hadrons is due to resonance decay.However, such an effect is not observed in the LBT approach, calling for a more detailed investigation of the role of excited states.For the mixed-hadronization, a similar pattern is observed, but also here the degree of feed-down varies substantially so that the previously decent agreement for the prompt H AA 's fans out substantially (where TAMU drops the most, as expected due to its larger extensive feed-down).Thus, settling the issue of feeddown contributions will play a critical role in establishing better consensus in the model approaches.

C. Elliptic Flow of Charm Hadrons
Another quantity, which is critical in the interpretation of experimental results, is the elliptic flow, v 2 , the second-order coefficient of the Fourier analysis of the azimuthal distribution of produced particles.Heavy quarks, initially produced in a hard process and therefore commonly assumed to carry vanishing v 2 (although that may not necessarily be correct in the presence of initial state correlations), are expected to acquire a finite v 2 by interactions with light partons from the medium, as the latter develop a finite v 2 during the expansion of the QGP due to an initial spatial eccentricity of the interaction zone.In addition, when the heavy quark recombines with a light quark to form a heavy meson, the underlying kinematics implies that the heavy mesons inherit further v 2 contributions from the light quarks (strictly speaking, these 2 contributions are not necessarily independent [63]).It is therefore desirable to assess these two contributions individually and to under- stand how the hadronization process modifies the v 2 .Here we focus on how the v 2 of the heavy meson is generated during hadronization.However, since no straightforward decomposition such as Eq. ( 47) can be given for the v 2 , considering the v 2 associated with an only-coalescence setting turns out to be problematic for the present analysis.We therefore focus on the well-defined cases of only-fragmentation and mixed hadronization.
We start by displaying in Fig. 11 the v 2 as a function of p T for D 0 (top row), D s (middle row), and Λ c (bottom row) including feeddown contributions, for the only-fragmentation scenario (left column) and the mixed-hadronization (right column).The color coding is the same as in Fig. 2, supplemented with a dashed line to the v 2 of the input charm-quark distribution.In almost all models, the onlyfragmentation v 2 of the heavy hadrons is very similar to that of the heavy quark.This is expected because the fragmentation into hadrons is collinear with the HQ direction.In the mixed scenario, we observe v 2 (Λ c ) > v 2 (D s ) > v 2 (D) for basically all models.This hierarchy is mainly caused by two reasons: (a) The v 2 of Λ c is larger than that of D s and D due to an additional light quark picked up via recombination; (b) The v 2 of D s is larger than that of D due to a smaller fraction of D s produced by fragmentation (relative to D 0 ), or, in other words, due to the larger recombination fraction facilitated by a higher concentration of strange quarks in the QGP relative to pp collisions [64].Passing from D → D s → Λ c , the agreement between the models becomes worse.For the Λ c , the extent to which the v 2 is increased over the meson-v 2 is close to a factor of 2 for Catania, LBT, and Torino, while it is less for TAMU due to its relatively large fragmentation fraction, recall Fig. 7. Thus, the relative increase of the baryon v 2 over the meson v 2 could potentially provide a measure of the fragmentation fraction of charm baryons in Pb-Pb collisions.
Finally, to have a more direct comparison of basic kinematics in v 2 production within the various models, we eliminate the effect of feeddown contributions and collect the results for the v 2 of the directly produced hadrons in Fig. 12. First, we recall that v 2 from the fragmentation process is essentially the same as that of the charm quark.Compared to Fig. 11, we see that directly produced hadrons have a considerably lower v 2 than for the case including feeddown from excited states.
In the real transport calculations that include the full QGP evolution, strong momentum-coordinate space correlations of the heavy quarks with the surrounding medium should be taken into account [15,17].As model studies show, these correlations may enhance v D 2 of the order of 50% in certain p T regions, especially at intermediate p T (more precisely in the transition region from recombination to fragmentation processes) where the mere increase of the recombination yield is the leading effect [15,17].

VI. CONCLUSIONS
In this paper, we have conducted specific comparisons between existing models of the hadronization of charm quarks.Hadronization is required to connect the computations of their transport through the QGP to experimental observables.Specifically, we have investigated how the yields and momen- tum spectra of produced HF hadrons emerge from the same underlying charm-quark distribution that includes both radial and elliptic flow.Most models employ a combination of recombination and fragmentation mechanisms, but the concrete implementations differ substantially.From the comparisons of the resulting hadron spectra, the following observations can be made: • The recombination probability, P rec (p c T ), and hence the fragmentation probability 1 − P rec (p c T ), varies strongly in the different approaches.An important cause of this difference is likely the role of excited states.In some models only ground-state hadrons are included, while others include a rich spectrum of excited states (in particular for baryons), predicted by relativistic quark models.This also extends the p c T reach of recombination substantially, which is especially important for the intermediate p c T region where recombination and fragmentation mechanisms compete.Finally, different normalization conditions are applied for the low-momentum limit of recombination.
• The shift in p T during the hadronization process, defined as p D T − p c T , is generally positive for recombination and negative for fragmentation processes.Surprisingly, the model variation in the former turns out to be smaller than in the latter (recall the comparison of Fig. 2 vs. Fig.6 (left panel)).
• The functional form of the hadronization ratio, H AA (i) = dNi/dp T dNc/dp T , at high p T is mostly controlled by the excited states which are included in the fragmentation process, but exhibits the largest model variations at small p T , especially when split up into different hadron species.Also here the excited states play a key factor for the difference, together with the model assumptions for P rec (p c T ) in the small-momentum limit.The difference between the various recombination models in how the collectivity of the c-quarks is converted into that of charm hadrons becomes particularly apparent when fo-  (46).For each charmed hadron, the results of the mixed hadronization process are plotted.cusing on direct production (i.e., no feeddown), where different shapes of the H AA at low p T exhibit the varying kinematics in the recombination models.
• The elliptic flow of the fragmentation component of charm hadrons agrees fairly well among the models, closely coinciding with that of the parent c-quark.The elliptic flow from the full hadronization models also agrees on a qualitative mass hierarchy for the maximal v 2 of D, D s , and Λ c , although the relative enhancement differs somewhat, especially for the Λ c , with the fragmentation fraction (and its feeddown contributions) playing an important role.For this observable the decay feeddown also appears to play a crucial role.
Future progress on these issues will likely require a more systematic resorting to basic principles (e.g., by developing closer connections to the properties of the QGP) and constraints, as well as the identification of suitable benchmarks in model comparisons (e.g., by scrutinizing the direct production of charm hadrons).Constraints from pp collisions, including multiplicity dependencies, for which high-precision data will become increasingly abundant, may help disentangle contributions from fragmentation and recombination processes.This could be very effective to constrain the hadrochemistry, both inclusive and as a function of p T and thus to reduce uncertainties associated with feeddown.The latter turned out to be a major source of the model variations in hadron observables, largely governed by the p T -dependent recombination fraction which controls the interplay of recombination and fragmentation.More ambitious constraints can be envisaged from lattice QCD, such as in-medium effects on the partonic and hadronic properties in the recombination pro-cess.

FIG. 1 .
FIG. 1. Illustration of the relative constituent-quark coordinates for mesons and baryons.
with g i = 6 the degeneracy factor.The essential change in the new version of the Nantes model is the light-quark mass.In the old version, m u,d = 0.1 GeV, m c = 1.5 GeV, while in the new version adopted since 2018, m u,d = 0.3 GeV, m c = 1.5

FIG. 2 .
FIG. 2. The normalized fragmentation function (c → D) used by the different groups.The Duke and Torino group use the fragmentation function of PYTHIA 6.4, which is pT dependent, as shown in Eq. 4;we set pT = 3 GeV and mc = 1.5 GeV for the plot.Here the solid lines in the LANL and Nantes models are pseudoscalar channels, and dashed lines are vector channels.
D + , D * 0 , D * + Ds, D * + s Λc, Λc(2595), Λc(2625), Σc(2455), Σc(2520) Duke D 0 , D + , D * 0 , D * + --LBT All S and P -wave D 0 and D + All S and P -wave Ds All S and P -wave Λc and Σc Nantes D 0 --PHSD Most S and P -wave D 0 Most S and P -wave Ds -D + , D * 0 , D * + --TABLE I. Number of excited states involved in each model.For the Torino model, the contribution of excited states is encoded in the decay of the clusters.

2 FIG. 3 .
FIG. 3. The transverse-momentum spectrum (upper panel) and elliptic flow, v2 (lower panel), of light quarks, u/d and s, at the hadronization hypersurface in the lab frame.

2 ] 2 FIG. 4 .
FIG.4.The transverse-momentum spectrum in the Lab frame (upper panel) and elliptic flow, v2, (lower panel) of charm quarks on the hadronization hypersurface.The black dashed line shows the v2 from the parameterized formula.

FIG. 5 .
FIG. 5.The recombination probability of a charm quark of transverse momentum p cT in the lab frame, into all charmed mesons (upper panel) and into all charmed hadrons (lower panel).

FIG. 6 .
FIG.6.The normalized spectrum, dND/dpT , of direct D 0 's resulting from c-quarks with fixed initial transverse momentum p c T =3.0, 10.0 GeV (indicated by the gray-dashed vertical lines).The only-coalescence and mixed-hadronization processes are plotted in the left and right panel, respectively.

FIG. 9 .
FIG.9.The yield ratio of Ds/D and Λc/D.For each charm hadron, the only-fragmentation, only-coalescence, and mixed-hadronization processes are plotted.

FIG. 10 .
FIG.10.The HAA of the direct D 0 , Ds, and Λc.For each charmed hadron, the only-coalescence and mixed-hadronization processes are plotted.

TABLE II .
Percentage distributions of 3 different ground-state hadron species to which a c-quark hadronizes.Presented are the results for only fragmentation (left 3 columns), only recombination (middle 3 columns), and mixed hadronization (right 3 columns).
a rescaled ratio, H r AA ≡ H AA /R with R = (dN c /dp T )H AA dp T (dN c /dp T )dp T