# Hadronization of heavy quarks

Jiaxing Zhao<sup>1</sup>, Jörg Aichelin,<sup>1</sup> Pol Bernard Gossiaux,<sup>1</sup> Andrea Beraudo<sup>2</sup>, Shanshan Cao,<sup>3</sup> Wenkai Fan,<sup>4</sup> Min He,<sup>5</sup>

Vincenzo Minissale<sup>6,7</sup> Taesoo Song<sup>8</sup>, Ivan Vitev<sup>9</sup>, Ralf Rapp,<sup>10</sup> Steffen Bass<sup>4</sup>, Elena Bratkovskaya,<sup>8,11,12</sup>

Vincenzo Greco<sup>6,7</sup> and Salvatore Plumari<sup>6,7</sup>

<sup>1</sup>SUBATECH, Nantes University, IMT Atlantique, Centre National de la Recherche Scientifique IN2P3,
4 rue Alfred Kastler, 44307 Nantes cedex 3, France

<sup>2</sup>Sezione di Torino, Istituto Nazionale di Fisica Nucleare, Via Pietro Giuria 1, I-10125 Torino, Italy

<sup>3</sup>Institute of Frontier and Interdisciplinary Science, Shandong University, Qingdao, Shandong 266237, China

<sup>4</sup>Department of Physics, Duke University, Durham, North Carolina 27708, USA

<sup>5</sup>Department of Applied Physics, Nanjing University of Science and Technology, Nanjing 210094, China

<sup>6</sup>Dipartimento di Fisica e Astronomia “Ettore Majorana,” Università degli Studi di Catania, Catania, Italy

<sup>7</sup>Laboratori Nazionali del Sud, Istituto Nazionale di Fisica Nucleare, Catania, Italy

<sup>8</sup>GSI Helmholtzzentrum für Schwerionenforschung GmbH, Planckstrasse 1, 64291 Darmstadt, Germany

<sup>9</sup>Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

<sup>10</sup>Cyclotron Institute and Department of Physics and Astronomy, Texas A&M University, College Station, Texas 77843, USA

<sup>11</sup>Institute for Theoretical Physics, Johann Wolfgang Goethe Universität, Frankfurt am Main, Germany

<sup>12</sup>Helmholtz Research Academy Hessen for FAIR, GSI Helmholtz Center for Heavy Ion Research, Campus Frankfurt,
60438 Frankfurt, Germany

<figure><img src="image_1.png" /></figure>

(Received 24 November 2023; accepted 29 April 2024; published 24 May 2024)

Heavy-flavor hadrons produced in ultrarelativistic heavy-ion collisions are a sensitive probe for studying hadronization mechanisms of the quark-gluon-plasma. In this paper, 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.

DOI: 10.1103/PhysRevC.109.054912

## 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 quark-gluon 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 preequilibrium 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 pseudocritical temperature, $T_{PC}$, the partons constituting the QGP undergo color neutralization 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 phases. A promising probe in this regard is heavy-flavor (HF) hadrons which have several advantageous features [3].

(1) The heavy-quark (HQ) mass is much larger than the nonperturbative QCD scale, $m_c$, $m_b \gg \Lambda_{QCD}$, and thus the production in nuclear collisions mainly occurs in initial hard scatterings, which can be reasonably well described by perturbative QCD.

(2) 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

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.</footnote>
diffusion is akin to a Brownian motion which enables
rather direct access to pertinent transport parameters.

(3) 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].

Over the last years major experimental efforts have been
underway to measure HF observables with increasing preci-
sion and breadth, in particular the production of various open
heavy-flavor hadrons, referred to as HF hadrochemistry (see,
e.g., Refs. [5–7]) in the charm sector. A large emphasis on
HF observables is also planned for upcoming heavy-ion runs
at RHIC with sPHENIX and the LHC where also the bottom
hadrochemistry will come into focus. Current results include
the observation of a marked increase of heavy-baryon yields
in hadronic collision systems, over those observed in the e⁻e⁺
collisions, which imply production mechanisms beyond uni-
versal fragmentation processes. It is the task for theoretical
approaches to cope with this development and advance to
more sophisticated approaches to make an interpretation if
these expected experimental results are possible.

Several research groups have developed transport ap-
proaches 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 [8–20].
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<sub>T</sub>) spectra
in heavy-ion collisions and proton-proton collisions, R<sub>AA</sub>(p<sub>T</sub>),
and the azimuthal asymmetry in the p<sub>T</sub> spectra quantified
by the elliptic-flow coefficient, v<sub>2</sub>(p<sub>T</sub>), which have become
increasingly precise in recent years [6,21–25]. These transport
approaches have become rather complex, aimed at a com-
prehensive 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 pro-
cesses cannot be rigorously described using first-principles
QCD as they occur in the nonperturbative regime, especially
at low transverse momentum, p<sub>T</sub>. 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 mea-
surements, 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. [26] benchmark HQ
interactions within different medium evolutions and different
hadronization schemes have been investigated, in Ref. [13]
different descriptions of the interaction of heavy quarks with
the QGP have been confronted, and in Ref. [27] the influence
of different initial-production mechanisms of heavy quarks
and expansion scenarios of the QGP has 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 [26]. With this paper, 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.
As such our paper does not contain genuinely new develop-
ments but is rather an attempt to highlight underlying model
assumptions and their manifestations in observable quanti-
ties under controlled conditions. Whereas the hadronization
of high-momentum light and heavy quarks is described by
fragmentation functions, the description of the hadronization
of heavy quarks with low momenta, which hadronize when
exiting a QGP, is a challenging task as heavy quarks are not
in equilibrium with the QGP partons. Therefore, the standard
procedure to hadronize the QGP partons via the Cooper-Frye
prescription [28] is not applicable and one has to resort to
more involved approaches, which invoke a combination of
coalescence and fragmentation mechanisms.

The current paper follows up on a sequence of previ-
ous papers [14,26,27], to which many groups which work
on transport approaches for heavy flavor dynamics in high-
energy heavy ion collisions have volunteered to contribute,
representing a large part of the international community.

The paper is organized as follows. In Sec. II, we pro-
vide 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
hypersurface 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, instantaneous coales-
cence 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 col-
lision 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*<sub>H</sub> factorizes into a hard production cross
section for a high-energy parton (*i*) and the fragmentation
fragmentation function, $D_{i \to H}$, as [29]

$$E \frac{d\sigma_H}{d^3P_H} = E_p \frac{d\sigma_i}{d^3p_i} \otimes D_{i \to H}(z), \quad (1)$$

where $E_i \frac{d\sigma_i}{d^3p_i}$ 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 $\otimes$ denotes a generic convolution, including the relevant fragmentation fraction $z$ Jacobian. The fragmentation function, $D_{i \to H}$, is a nonperturbative 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. [30]). The most frequently used fragmentation function for heavy quarks is from Peterson [31] in which the HF meson carries most of the momentum of the heavy quark $Q$:

$$D_{Q \to H}(z) \propto \frac{1}{z\left[1 - \frac{1}{z} - \frac{\epsilon}{1-z}\right]^2}, \quad (2)$$

where $\epsilon = m_Q^2/m_Q^2$, 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 \gg \Lambda_{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,32]

$$\begin{aligned} D_{Q \to H_p} & \propto \frac{r z (1-z)^2}{[1 - (1-r)z]^6} [6 - 18(1-2r)z \\ & \quad + (21 - 74r + 68r^2)z^2 \\ & \quad - 2(1-r)(6 - 19r + 18r^2)z^3 \\ & \quad + 3(1-r)^2(1 - 2r + 2r^2)z^4], \\ D_{Q \to H_V} & \propto \frac{r z (1-z)^2}{[1 - (1-r)z]^6} [2 - 2(3-2r)z \\ & \quad + 3(3-2r+4r^2)z^2 - 2(1-r)(4-r+2r^2)z^3 \\ & \quad + 3(1-r)^2(3-2r+2r^2)z^4]. \end{aligned} \quad (3)$$

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 [33]. It can be expressed as

$$D_{Q \to H} \propto \frac{1}{z^{1+r b m_Q^2}} z^{a_\alpha} \left(\frac{1-z}{z}\right)^{a_\beta} \exp\left(-\frac{b m_T^2}{z}\right), \quad (4)$$

where $r$, $a_\alpha$, $a_\beta$, and $b$ are free parameters, and $m_T = \sqrt{m_Q^2 + p_T^2}$ is the transverse energy. In PYTHIA, $r = 1.32$, $a_\alpha = a_\beta = 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/\pi$ and $\Lambda/K$ ratios, observed in heavy-ion collisions at both RHIC [5,34,35] and the LHC [7,36], and an approximate quark number scaling of the elliptic flow [21,37–39]. The parton recombination model was proposed and successfully used to describe these phenomena [40–43], primarily in the regime of intermediate $p_T \lesssim 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 equal-temperature hypersurface and that the coalescence probability is given by a Wigner function [42,44]. Given that HF observables exhibit similar features as the aforementioned light-quark observables (e.g., the $\Lambda_c/D^0$ ratio exhibits an enhancement [7] akin to the $p/\pi$ or $\Lambda/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 the RRM (Sec. II B 2).

### 1. Instantaneous coalescence (ICM)

In the sudden coalescence models, the momentum spectrum of the produced hadron $H$ can be written as

$$\frac{dN_H}{d^3\mathbf{P}} = g_H \int \prod_{i=1}^{N} \frac{d^3 p_i}{(2\pi)^3 E_i} p_i \cdot d\sigma_i F(\mathbf{x}_1 \dots \mathbf{x}_N, \mathbf{p}_1 \dots \mathbf{p}_N) \\ \times W(\mathbf{x}_1 \dots \mathbf{x}_N, \mathbf{p}_1 \dots \mathbf{p}_N) \delta^{(3)}\left(\mathbf{P} - \sum_{i=1}^{N} \mathbf{p}_i\right), \quad (5)$$

where $g_H$ is a degeneracy factor of color and spin; $\mathbf{P}$ and $\mathbf{p}_i$ are the momenta of HF hadron and constituents quarks, respectively. The $\delta$ function ensures the conservation of three-momentum, $d\sigma_i$ is a hypersurface element, and $F(\mathbf{x}_1 \dots \mathbf{x}_N, \mathbf{p}_1 \dots \mathbf{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(\mathbf{x}_1 \dots \mathbf{x}_N, \mathbf{p}_1 \dots \mathbf{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

$$W(\mathbf{r}, \mathbf{p}_r) = \int d^3\mathbf{y} e^{-i\mathbf{p}_r \cdot \mathbf{y}} \psi\left(\mathbf{r} + \frac{\mathbf{y}}{2}\right) \psi^*(\mathbf{r} - \frac{\mathbf{y}}{2}), \quad (6)$$

where $\mathbf{r}$ and $\mathbf{p}_r$ are the relative distance and momentum in the center-of-mass (CM) frame of the hadron, respectively. The wave function $\psi$ can, e.g., be obtained by solving the two-body 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
<figure><img src="image_2.png" /><figcaption>FIG. 1. Illustration of the relative constituent-quark coordinates for mesons and baryons.</figcaption></figure>

root-mean-square (rms) radius as the original Wigner density:

$$W(\mathbf{r}, \mathbf{p}_r) = 8 \exp \left( -\frac{r^2}{\sigma^2} - p_r^2 \sigma^2 \right). \quad (7)$$

The parameter $\sigma$ in the Wigner function is related to the rms radius of the hadron, $\sigma^2 = 2\langle r^2 \rangle/3$. Some authors relate the width to the rms charge radius [45,46] instead of mass radius [see, e.g., 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. [47]). There are two relative coordinates, $\rho$ and $\lambda$, as shown in Fig. 1. The relative wave function of the three-body system in the ground state is assumed to be factorized as $\Psi = \psi(\rho)\psi(\lambda)$. Consequently, the Wigner density can be expressed as a product of two Gaussians:

$$W(\rho, \lambda, \mathbf{p}_{\rho}, \mathbf{p}_{\lambda}) = 8^2 \exp \left( -\frac{\rho^2}{\sigma_{\rho}^2} - p_{\rho}^2 \sigma_{\rho}^2 \right) \times \exp \left( -\frac{\lambda^2}{\sigma_{\lambda}^2} - p_{\lambda}^2 \sigma_{\lambda}^2 \right), \quad (8)$$

where $\mathbf{p}_{\rho}$ and $\mathbf{p}_{\lambda}$ are the relative momenta corresponding to the relative coordinates, $\rho$ and $\lambda$. The widths $\sigma_{\rho}$ and $\sigma_{\lambda}$ are related to the inner radii via

$$\sigma_{\rho}^{2} = \frac{2}{3} \langle \rho^{2} \rangle, \quad \sigma_{\lambda}^{2} = \frac{2}{3} \langle \lambda^{2} \rangle. \quad (9)$$

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 feed-down 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 three-dimensional isotropic harmonic oscillator.

The probability density, $dN_H(\mathbf{P}, \mathbf{p}_Q)$, to produce a heavy hadron, $H$, with momentum $\mathbf{P}$ from a heavy quark, $Q$, with momentum $\mathbf{p}_Q$ is defined by introducing a delta distribution centered at $\mathbf{p}_Q$ into $F$ in Eq. (5). One then defines the total recombination probability for this heavy quark to hadronize through coalescence as

$$P_{\text{rec}}(\mathbf{p}_Q) = \sum_H \int d\mathbf{P} \, dN_H(\mathbf{P}, \mathbf{p}_Q). \quad (10)$$

While this quantity is boost invariant, it is not the case for the spectra defined at Eq. (5), due to the violation of the energy-momentum conservation inherent to the coalescence process. A particular frame thus needs to be chosen, usually taken as the local fluid rest frame.

## 2. Resonance recombination (RRM)

The RRM has been derived from a Boltzmann equation where a recombination rate is obtained by utilizing a resonant quark antiquark cross section (or pertinent scattering matrix element) for hadron formation [48], 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 linear-response theory. In a similar spirit, a different recombination criterion has recently been proposed [20] 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 four-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. [9] and extended to HF baryons in Ref. [18].

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,

$$\left(\frac{\partial}{\partial t} + \mathbf{v} \cdot \nabla\right) f_M(\mathbf{R}, \mathbf{P}) = -\frac{\Gamma}{\gamma} f_M(\mathbf{R}, \mathbf{P}) + \beta(\mathbf{R}, \mathbf{P}), \quad (11)$$

where $\mathbf{R}$ and $\mathbf{P}$ are the position and momentum of the meson, $\gamma = E_M/m$ is a Lorentz factor with meson mass $m$ and energy $E_M$, and $\Gamma$ 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

$$-\int \frac{d^3 \mathbf{R} d^3 \mathbf{P}}{(2\pi)^3} \frac{\Gamma}{\gamma} f_M^{\text{eq}} + \int \frac{d^3 \mathbf{R} d^3 \mathbf{P}}{(2\pi)^3} \beta = 0. \quad (12)$$

The integrated meson yield can then be expressed as

$$\begin{aligned} N_M &= \int \frac{d^3 \mathbf{R} d^3 \mathbf{P}}{(2\pi)^3} f_M^{\text{eq}} \\ &= \frac{\gamma}{\Gamma} \int \frac{d^3 \mathbf{R} d^3 \mathbf{P}}{(2\pi)^3} \int \frac{d^3 \mathbf{p}_r}{8(2\pi)^3} f_q(\mathbf{R}, \mathbf{P}, \mathbf{p}_r) f_{\bar{q}}(\mathbf{R}, \mathbf{P}, \mathbf{p}_r) \\ &\quad \times \sigma(s) v_{\text{rel}}, \end{aligned} \quad (13)$$

where $f_q$ and $f_{\bar{q}}$ are the phase-space distributions of quark and antiquark, $\mathbf{p}_r = \mathbf{p}_1 - \mathbf{p}_2$ and $\mathbf{P} = \mathbf{p}_1 + \mathbf{p}_2$ are their relative and total momentum, and $v_{\text{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:

$$\sigma(s) = g_{\sigma} \frac{4\pi}{k^2} \frac{(\Gamma m)^2}{(s - m^2)^2 + (\Gamma m)^2}, \quad (14)$$

where $g_{\sigma} = g_M/(g_q g_{\bar{q}})$ is the spin-color degeneracy, and $k$ denotes the three-momentum of the quarks (heavy or light) in the CM frame. $s$ is the invariant mass. For baryons, a two-step
<table><caption>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.</caption><thead><tr><th></th><th>D</th><th>D<sub>s</sub></th><th>Λ<sub>c</sub></th></tr></thead><tbody><tr><td>Catania</td><td>D<sup>0</sup>, D<sup>+</sup>, D<sup>*0</sup>, D<sup>*+</sup></td><td>D<sub>s</sub>, D<sub>s</sub><sup>*+</sup></td><td>Λ<sub>c</sub>, Λ<sub>c</sub>(2595), Λ<sub>c</sub>(2625), Σ<sub>c</sub>(2455), Σ<sub>c</sub>(2520)</td></tr><tr><td>Duke</td><td>D<sup>0</sup>, D<sup>+</sup>, D<sup>*0</sup>, D<sup>*+</sup></td><td>All S- and P-wave D<sub>s</sub></td><td>All S- and P-wave Λ<sub>c</sub> and Σ<sub>c</sub></td></tr><tr><td>LBT</td><td>All S- and P-wave D<sup>0</sup> and D<sup>+</sup></td><td>All S- and P-wave D<sub>s</sub></td><td>All S- and P-wave Λ<sub>c</sub> and Σ<sub>c</sub></td></tr><tr><td>Nantes</td><td>D<sup>0</sup></td><td>Most S- and P-wave D<sup>0</sup></td><td>Most S- and P-wave D<sub>s</sub></td></tr><tr><td>PHSD</td><td>Most S- and P-wave D<sup>0</sup></td><td>Most S- and P-wave D<sub>s</sub></td><td>Most S- and P-wave D<sub>s</sub></td></tr><tr><td>TAMU</td><td>PDG</td><td>PDG</td><td>RQM</td></tr><tr><td>Torino</td><td>D<sup>0</sup></td><td>D<sub>s</sub></td><td>Λ<sub>c</sub></td></tr><tr><td>LANL</td><td>D<sup>0</sup>, D<sup>+</sup>, D<sup>*0</sup>, D<sup>*+</sup></td><td>D<sub>s</sub></td><td>Λ<sub>c</sub></td></tr></tbody></table>

process has been employed [18], 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 mostly rather well defined and understood, different models vary in their implementation of the respective hadronization schemes. This fact limits the predictive power of transport approaches and is a consequence of the lack of fundamental QCD-based methods for evaluating nonperturbative process. The main objective of studying the hadronization models in the present paper, relying on a combination of fragmentation and coalescence processes, is essentially restricted to testing how the different hadronization schemes respond to a common well-defined configuration of a QGP akin to what one encounters in heavy-ion collisions (same freeze-out hypersurface, same distribution of heavy quarks).

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 [8,9,11,14–20,26,49]. We emphasize that all these quantities have not been modified with respect to the original publications with a few exceptions which will be explicitly mentioned below. Here we limit ourselves to the comparison of those quantities that directly influence the hadronization scheme. Key features 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; and (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 frame; (b) the fluid rest frame, usually based on a hydrodynamic description of the expansion of the hot medium (different regions of the medium have different fluid velocities), which is the stationary frame inside the flowing fluid cell; and (c) the CM frame of the light and heavy quark or the diquark and the heavy quark (in the case of baryons).

As decay feed-downs 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 $\Lambda_c$) in the following analyses, we have included a table summarizing the excited states accounted for by each model in Table 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):

$$W(r_i, p_{ri}) = \prod_{i=1}^{N_q-1} A_W \exp \left( -\frac{r_i^2}{\sigma_i^2} - p_{ri}^2 \sigma_i^2 \right), \quad (15)$$

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 \to 0$, all charm quarks hadronize by coalescence into a heavy hadron. This is imposed by requiring that the total coalescence probability gives $\lim_{p \to 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

$$r_1 = |\mathbf{x}_1 - \mathbf{x}_2|, \quad p_{r1} = \frac{|m_2 \mathbf{p}_1 - m_1 \mathbf{p}_2|}{m_1 + m_2}, \quad (16)$$

while for baryons they are defined as

$$r_1 = \frac{|\mathbf{x}_1 - \mathbf{x}_2|}{\sqrt{2}}, \quad p_{r1} = \sqrt{2} \frac{|m_2 \mathbf{p}_1 - m_1 \mathbf{p}_2|}{m_1 + m_2} \quad (17)$$

for the diquark, and $r_2$ and $p_{r2}$ for the quark-diquark are given by

$$r_2 = \sqrt{\frac{2}{3}} \left| \frac{m_1 \mathbf{x}_1 + m_2 \mathbf{x}_2}{m_1 + m_2} - \mathbf{x}_3 \right|, \\ p_{r2} = \sqrt{\frac{3}{2}} \frac{|m_3(\mathbf{p}_1 + \mathbf{p}_2) - (m_1 + m_2)\mathbf{p}_3|}{m_1 + m_2 + m_3}, \quad (18)$$

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, $\langle r^2 \rangle_{\text{ch}} = \sum_{i=1}^{N} ((x_i - R_{\text{CM}})^2)$ with N = 2 and 3 for mesons and baryons, respectively; $R_{\text{CM}}$ is the CM coordinate. For mesons one has

$$ \langle r^2 \rangle_{\text{ch}} = \frac{3}{2} \frac{Q_1 m_2^2 + Q_2 m_1^2}{(m_1 + m_2)^2} \sigma_1^2, \quad (19) $$

and for baryons

$$ \langle r^2 \rangle_{\text{ch}} = \frac{3}{2} \frac{Q_1 m_2^2 + Q_2 m_1^2}{(m_1 + m_2)^2} \sigma_1^2 + \frac{3}{2} \frac{(Q_1 + Q_2)m_3^2 + Q_3(m_1 + m_2)^2}{(m_1 + m_2 + m_3)^2} \sigma_2^2, \quad (20) $$

where $Q_i$ is the charge of the ith quark. The charge radius of the hadrons is taken according to the constituent-quark model. In the Catania model one has $\langle r^2 \rangle_{\text{ch}} = 0.184 \text{ fm}^2$ for $D^+$, meaning $\sigma = \sigma_1 = 3.546 \text{ GeV}^{-1}$ for D mesons, and $\langle r^2 \rangle_{\text{ch}} = 0.083 \text{ fm}^2$ for $D_s^+$, meaning $\sigma_1 = 2.475 \text{ GeV}^{-1} \text{ GeV}$ for $D_s$ mesons. For $\Lambda_c$, the rms charge radius is taken as $\langle r^2 \rangle_{\text{ch}} = 0.15 \text{ fm}^2$. Actually, for baryons, there is only one free parameter, because the two widths are related by the oscillator frequency $\omega$ through the reduced masses, $\mu_i$, via $\sigma_i = 1/\sqrt{\mu_i\omega}$:

$$ \begin{aligned} \mu_1 &= \frac{m_1 m_2}{m_1 + m_2}, \\ \mu_2 &= \frac{(m_1 + m_2) m_3}{m_1 + m_2 + m_3}, \end{aligned} \quad (21) $$

where $\mu_1$ is the reduced mass of the two-body system and $\mu_2$ is the reduced mass of the two-body system with the third particle. The corresponding widths are $\sigma_1 = 5.556 \text{ GeV}^{-1}$ and $\sigma_2 = 2.924 \text{ GeV}^{-1}$.

In the Catania model, all first excited states for D mesons and $\Lambda_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 + \Delta E$, at temperature T and a statistical factor, $(m_R/m_G)^{3/2} \times \exp(-(m_R - m_G)/T)$, 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 laboratory frame satisfies

$$ f_q = \frac{g \tau m_T}{(2\pi)^3} \exp \left( - \frac{\gamma_T (m_T - p_T \cdot \beta_T)}{T} \right), \quad (22) $$

where $m_T = \sqrt{p_T^2 + m_q^2}$, and $\beta_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 antiquark pairs according to the flavor compositions; $\tau$ 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 \text{ GeV}$, $m_s = 0.38 \text{ GeV}$, and $m_c = 1.5 \text{ GeV}$.

Heavy quarks that do not hadronize via coalescence are converted to hadrons by fragmentation. Catania uses the Perterson fragmentation with $\epsilon = 0.1$ for charmed mesons and $\epsilon = 0.02$ for charmed baryons.

For more details, we refer the reader to Refs. [10,49].

## 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,

$$ W(p_r) = g_h \frac{(2\sqrt{\pi}\sigma)^3}{V} e^{-\sigma^2 p_r^2}, \quad (23) $$

where $p_r$ is the relative momentum defined in the CM frame:

$$ p_r = \frac{|E_2 \mathbf{p}_1 - E_1 \mathbf{p}_2|}{E_1 + E_2}. \quad (24) $$

The width in the Gaussian distribution can be obtained by $\sigma = 1/\sqrt{\mu\omega}$ with $\mu$ being the reduced mass of the light and heavy quark. Duke takes $\omega = 0.106 \text{ GeV}$ for all charmed hadrons. So, the related widths are $\sigma(D) = 6.23 \text{ GeV}^{-1}$ and $\sigma(D_s) = 5.22 \text{ GeV}^{-1}$. Both ground states and first excited states of D mesons are considered. However, the hadronization implementation utilized for this paper [50] does not include baryonic states. Later implementations of the Duke hadronization model [11] do include baryons such as $\Lambda_c$, $\Sigma_c$, and $\Omega_c$.

The light-quark distribution in the fluid rest frame satisfies

$$ f_q = \frac{V}{(2\pi)^3} \frac{g_i}{e^{E_i/T} + 1}, \quad (25) $$

where $g_i = 6$ is the statistical factor. The quark masses are taken as $m_{u,d} = 0.3 \text{ GeV}$, $m_s = 0.475 \text{ GeV}$, and $m_c = 1.27 \text{ GeV}$. 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 Refs. [11,50].

## 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

$$ \begin{aligned} W_S(p_r) &= g_h \frac{(2\sqrt{\pi}\sigma)^3}{V} e^{-\sigma^2 p_r^2}, \\ W_P(p_r) &= g_h \frac{(2\sqrt{\pi}\sigma)^3}{V} \frac{2}{3} \sigma^2 p_r^2 e^{-\sigma^2 p_r^2}, \end{aligned} \quad (26) $$

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:

$$p_r = \frac{|E_2\mathbf{p}_1 - E_1\mathbf{p}_2|}{E_1 + E_2}, \qquad (27)$$

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 three-dimensional calculation, are included, which covers nearly all charmed hadron species ($\Xi_c$ and $\Omega_c$ are also included) as reported by the particle data group (PDG).

The width of the Wigner function is determined by $\sigma = 1/\sqrt{\mu\omega}$ with $\mu$ the reduced mass and $\omega$ the oscillator frequency, $\omega = 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 1 when all S- and 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, and $m_c = 1.8$ GeV, and the related widths are $\sigma(D) = 4.02 \text{ GeV}^{-1}$ and $\sigma(D_s) = 3.57 \text{ GeV}^{-1}$. For $\Lambda_c$, there are two widths, $\sigma_1 = 4.93 \text{ GeV}^{-1}$ and $\sigma_2 = 2.87 \text{ 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

$$f_q = \frac{V}{(2\pi)^3} \frac{g_i}{e^{E_i/T} + 1}, \qquad (28)$$

where $g_i = 6$ is the statistical factor.

The traditional instantaneous coalescence model only respects the three-momentum conservation but violates the energy conservation. A solution of this problem is proposed in Ref. [14] 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 off-shell 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 $\epsilon = 0.05$ for charmed hadrons.

For more details on the hadronization model applied to the LBT calculation, please see Ref. [14].

### D. Nantes

In the Nantes model, heavy quarks hadronize either by coalescence or by fragmentation. The recombination probability is evaluated in the thermal rest frame and normalized to 1 for vanishing HQ momentum in that frame. Thus, in the laboratory 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 [51].

The coalescence probability is given by

$$W(x_Q, x_q, p_Q, p_q) = \exp\left(\frac{(x_q - x_Q)^2 + [(x_q - x_Q) \cdot u_Q]^2}{2R_c^2}\right) F_{\Phi}(p_Q, p_q), \qquad (29)$$

where $u_Q$ is the four-velocity of the heavy quark and

$$F_{\Phi}(p_Q, p_q) = \exp\left(\frac{(p_Q/m_Q - p_q/m_q)^2}{2\alpha_d^2}\right) \qquad (30)$$

with $\alpha_d = 0.51$; $R_c$ is fixed for a given $a_d$ from the normalization condition

$$\frac{R_c}{\sqrt{2\pi}\hbar} = \frac{m_Q + m_q}{m_Q m_q} \times \left( 4\pi\alpha_d^2 K_2\left(\frac{1}{\alpha_d^2}\right) e^{\frac{1}{\alpha_d^2}} \right)^{-1/3} \qquad (31)$$

In coordinate space, one can obtain $\sigma_r = \sqrt{2}R_c$, but in momentum space one has $\sigma_p = 1/(\sqrt{2}\alpha_d\mu)$ by definition of the $F_{\Phi}$, with $\mu$ as the reduced mass. The widths, $\sigma$, 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:

$$f_q = g_i \exp(-\sqrt{m^2 + p^2}/T) \qquad (32)$$

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 and $m_c = 1.5$ GeV, while in the new version adopted since 2018, $m_{u,d} = 0.3$ GeV and $m_c = 1.5$ 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,32].

For more details, we refer the reader to Ref. [15].

### 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 \times 10 \times 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 the antiquark within the dynamical quasiparticle model [52,53], that is, in the fluid rest frame:

$$\frac{df_{\bar{q}}}{dm} = \frac{2m^2\gamma}{\pi(m^2 - M^2)^2 + (2m\gamma)^2} \frac{g}{e^{E_i/T} + 1} \qquad (33)$$

with $g = 6$ the color-spin degeneracy factor; M and $\gamma$ are the pole mass and spectral width of the antiquark, respectively.
The coalescence probabilities for the S-wave and P-wave states are given by

$$ 
\begin{aligned}
W_S(r, p_r) &= \frac{8(2s + 1)}{36} e^{-\frac{r^2}{\sigma^2} - \sigma^2 p_r^2}, \\
W_P(r, p_r) &= \frac{2s + 1}{36} \left( \frac{16}{3} \frac{r^2}{\sigma^2} + \frac{16}{3} \sigma^2 p_r^2 - 8 \right) e^{-\frac{r^2}{\sigma^2} - \sigma^2 p_r^2},
\end{aligned}
\quad (34) $$

where s is the spin of the meson. The relative distance and momentum in the CM frame are, respectively,

$$ 
\begin{aligned}
r &= |\mathbf{r}_1 - \mathbf{r}_2|, \\
p_r &= \frac{|m_2\mathbf{p}_1 - m_1\mathbf{p}_2|}{m_1 + m_2}.
\end{aligned}
\quad (35) $$

The width σ is related to the rms radius of the D meson through

$$ \sigma_s^2 = \frac{4(m_1 + m_2)^2}{3(m_1^2 + m_2^2)} \langle r_M^2 \rangle \quad (36) $$

for the S-wave state and

$$ \sigma_p^2 = \frac{3}{5} \sigma_s^2 \quad (37) $$

for the P-wave state, and $\sqrt{\langle r_M^2 \rangle} = 0.9$ fm for both S-wave 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) [52,53]. The relative momentum coordinates, Eq. (35), are defined as in the Catania model, Eq. (16), while for the Duke and LBT, Eqs. (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 $\mathbf{p}_1 = -\mathbf{p}_2$.

Only S- and P-wave mesons are taken into account for coalescence (not taking into account the principle quantum number), while charm baryons such as $\Lambda_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 [16,17,54].

If the heavy quarks do not hadronize via the coalescence process, they will convert into charmed mesons or $\Lambda_c$ via fragmentation. The Peterson fragmentation function with $\epsilon = 0.01$ for charmed mesons is used in the PHSD model.

For more details, we refer the reader to Refs. [16,17].

## F. TAMU

The TAMU model applies a combined recombination and fragmentation approach as well. The former is realized via the RRM [48], where the recombination probability for the two-body case is controlled by resonance amplitudes and is currently expressed as a relativistic Breit-Wigner cross sec-

tion:

$$ \sigma(s) = g_{\sigma} \frac{4\pi}{k^2} \frac{(\Gamma m)^2}{(s - m^2)^2 + (\Gamma m)^2}, \quad (38) $$

where $g_σ = g_M/(g_q g_{\bar{q}})$ is a statistical weight, k denotes the quark three-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

$$ P_{\text{rec}}(p_c^*) = P_0 \int \frac{d^3 p_q}{(2\pi)^3} f_q(p_q) \frac{\gamma_M}{\Gamma_M} \sigma(s) v_{\text{rel}}, \quad (39) $$

and for baryons as

$$ P_{\text{rec}}(p_c^*) = P_0 \int \frac{d^3 p_{q1} d^3 p_{q2}}{(2\pi)^6} f_{q1}(p_{q1}) f_{q2}(p_{q2}) \frac{\gamma_B \gamma_{[qq]}}{\Gamma_B \Gamma_{[qq]}} \times \sigma(s_{12}) v_{\text{rel}}^{12} \sigma(s_{[qq]}) v_{\text{rel}}^{[qq]c}, \quad (40) $$

where $f_q$ is the thermal light-quark distribution

$$ f_q = g_i e^{-E/T} \quad (41) $$

with a degeneracy factor $g_i = 6$; $v_{\text{rel}}^{12}$ and $v_{\text{rel}}^{[qq]c}$ are the relative velocities between the two light quarks, which form a diquark, and between the diquark and the charm quark, respectively [recall that $v_{\text{rel}}^{12} = \sqrt{(p_1^\mu p_{2\mu})^2 - m_1^2 m_2^2 / (E_1 E_2)}$]. 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 $\Gamma_M = 0.1$ GeV, $\Gamma_{[qq]} = 0.2$ GeV, and $\Gamma_B = 0.3$ GeV; $\gamma_M$, $\gamma_{[qq]}$, and $\gamma_B$ denote 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 includes all states listed by the PDG [55], but goes well beyond that in terms of “missing” charm-baryon states as predicted by the relativistic-quark model (RQM) [56] and lattice QCD [57,58]. The excited $\Lambda_c$ states and $\Sigma_c$ states are assumed to feed the ground state $\Lambda_c$ with a 100% branching fraction.

Different from the most recent version of RRM [18], 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 paper. 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_T^c$ in the laboratory 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_{\text{rec}}(p_c^*)$, is then calculated via Eqs. (39) and (40), i.e., for a $p_T^c$ we obtain the corresponding recombination probability $P_{\text{rec}}(p_c^*)$. Finally, the recombination probability $P_{\text{rec}}(p_T^c)$ in the laboratory frame can be obtained by dividing the sum of probabilities $P_{\text{rec}}(p_T^c)$ by the number of incoming particles with $p_T^c$.
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_{\text{rec}}(p_T^c)$, which is partitioned over all $D, D_s$, and $\Lambda_c$ baryons according to their respective fractions determined in pp collisions [55].

For more details, we refer the reader to Refs. [9,18,48,55].

## 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

$$n = g_s g_I \frac{T_H M^2}{2\pi^2} \sum_{n=1}^{\infty} \frac{(\pm 1)^{n+1}}{n} K_2\left(\frac{nM}{T_H}\right), \quad (42)$$

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.579$ GeV, $m_{sl} = 0.804$ GeV, $m_{(ss)} = 1.0936$ GeV,...) and the $-/+$ sign refers to fermions/bosons, respectively.

(2) Once the particle species is selected, determine its three-momentum 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, resample the medium particle. If not, after defining an intermediate cutoff $M_{\text{max}}$ set to $M_{\text{max}} = 3.8$ GeV, there are two possibilities.

(a) If $M_C < M_{\text{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.

(b) If $M_C > M_{\text{max}}$ the cluster hadronizes via string fragmentation simulated through PYTHIA 6.4.

Only ground-state hadrons, like $D^0, D^+, D_s, \Lambda_c, \Xi_c$, and $\Omega_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. [20,59], 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,32].

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 $\approx 50\%$ [60]. Contributions to parton energy loss and in-medium parton showers in general, to hadron production, can be accounted for in different but related ways [61]. In the LANL model an effective modification of the fragmentation function [62,63] is employed as

$$\mathcal{D}_{H/c}(z) \Rightarrow \int_0^{1-z} d\epsilon P(\epsilon) \frac{1}{1-\epsilon} \mathcal{D}_{H/c}\left(\frac{z}{1-\epsilon}\right) + \int_z^1 d\epsilon \frac{dN_g}{d\epsilon}(\epsilon) \frac{1}{\epsilon} \mathcal{D}_{H/g}\left(\frac{z}{\epsilon}\right), \quad (43)$$

where $P(\epsilon)$ is the probability that the heavy quark loses a fraction of its energy, $\epsilon$, due to multiple gluon emissions, and $dN_g/d\epsilon$ is the distribution of the average gluons as a function of $\epsilon = \omega/E$; $\mathcal{D}_{H/c}$ is the fragmentation function of the charm quark into charmed hadrons, whereas $\mathcal{D}_{H/g}$ is the fragmentation function of gluons into charm hadrons. For more details, please see Refs. [19,60].

### 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 [31] with the parameter $\epsilon = 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
<figure><img src="image_3.png" /><figcaption>FIG. 2. The normalized fragmentation function ($c \to D$) used by the different groups. The Duke and Torino groups use the fragmentation function of PYTHIA 6.4, which is $p_T$ dependent, as shown in Eq. (4); we set $p_T = 3$ GeV and $m_c = 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.</figcaption></figure>

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.

## 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 paper, the expanding QGP is modeled by a simple fireball model [8]. The fireball is cylindrically oriented along the z direction with an elliptic cross-sectional area. The long (a) and short (b) semiaxes of the ellipse are oriented, respectively, along y and x direction and are given by

$$ 
\begin{aligned}
a(t) = a_0 &+ v_\infty \left[ t - \frac{1 - \exp(-At)}{A} \right] \\ 
& - \Delta v \left[ t - \frac{1 - \exp(-Bt)}{B} \right], \\ 
b(t) = b_0 &+ v_\infty \left[ t - \frac{1 - \exp(-At)}{A} \right] \\ 
& + \Delta v \left[ t - \frac{1 - \exp(-Bt)}{B} \right].
\end{aligned} 
\quad (44) 
$$

The velocity of the fireball surface along the long and short axis is given by $v_a(t) = \dot{a}(t)$ and $v_b(t) = \dot{b}(t)$. The parameters are taken as $a_0 = 5.562$ fm, $b_0 = 4.450$ fm, $v_\infty = 0.52$, $\Delta v = 0.122$, $A = 0.55/fm$, and $B = 1.3/fm$ for Pb-Pb collisions at

<figure><img src="image_4.png" /><figcaption>FIG. 3. The transverse-momentum spectrum (upper panel) and elliptic flow, $v_2$ (lower panel), of light quarks, u/d and s, at the hadronization hypersurface in the laboratory frame.</figcaption></figure>

$\sqrt{s_{NN}} = 2.76$ TeV and impact parameter $b = 7$ fm, resulting in $v_b(t) > v_a(t)$.

Hadronization is assumed to occur on a hypersurface $R_\mu = (t, \mathbf{R})$, which is a four-dimensional surface with equal temperature $T(R_\mu) = T_H$. This corresponds to a hadronization time $\tau = \tau(T_c, \mathbf{R})$ in the above setup, 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 pseudocritical temperature of the chiral transition (around $T_{PC} \simeq 155$ MeV) [57], 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 x and y directions, the light quarks boosted to the laboratory frame will carry a nonzero elliptic flow, which is defined as the second-order coefficient of the Fourier expansion of the azimuthal distribution. The resulting transverse-momentum
<figure><img src="image_5.png" /><figcaption>FIG. 4. The transverse-momentum spectrum in the laboratory frame (upper panel) and elliptic flow, $v_2$ (lower panel), of charm quarks on the hadronization hypersurface. The black dashed line shows the $v_2$ from the parametrized formula.</figcaption></figure>

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 distribution in the laboratory frame is taken from the EMMI rapid reaction task framework [26], displayed in the top panel of Fig. 4. The distribution can be parametrized as $dN/dp_T = 53.15(p_T/\text{GeV})/[4.423 + (p_T/\text{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 paper, we neglect these correlations: charm quarks are sampled isotropically in each fluid cell at the hypersurface.

In such prepared environment, the different groups carried out four different calculations.

(1) The yield of different charm hadrons is calculated considering (a) only fragmentation (assuming that all heavy hadrons are created by fragmentation), (b) only coalescence/recombination (all hadrons are produced by coalescence/recombination), and (c) the full framework of fragmentation plus coalescence, as realized

in the respective model. For each case, provide a “hadronization ratio” [26]:

$$ 
\begin{align} 
H_{AA}(D) &= \frac{dN_D/dp_T}{dN_c/dp_T}, \\ 
H_{AA}(D_s) &= \frac{dN_{D_s}/dp_T}{dN_c/dp_T}, \\ 
H_{AA}(\Lambda_c) &= \frac{dN_{\Lambda_c}/dp_T}{dN_c/dp_T}, 
\end{align} 
\quad (45) $$

where D includes $D^0$ and $D^+$. $D_s$ and $\Lambda_c$ are the prompt yields (direct and strong decay).

(2) The elliptic flow $v_2$ of all charmed hadrons with an anisotropic charm-quark spectrum is calculated according to

$$ \frac{dN_c}{d\bar{p}_T} \Rightarrow \frac{dN_c}{dp_T} \times [1 + 2v_2(p_T) \cos(2\phi)], \quad (46) $$

with $v_2(p_T) = 0.0577p_T/(1 + 0.4212p_T^{1.6367})$, 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) and (b) the full hadronization scheme used by the respective model (recombination and fragmentation).

(3) The $H_{AA}$ and $v_2$ of directly produced $D^0$, $D_s$, and $\Lambda_c$ (without feed-down contributions) are calculated.

(4) $dN(D^0)/dp_T$ of direct $D^0$ mesons (without feed-down contributions) produced by a c quark with $p_T^c = 3$ and 10 GeV is calculated.

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_T^c$ as

$$ dN_H^{\text{mixed}}(\mathbf{p}_T, \mathbf{p}_T^c) = dN_H(\mathbf{p}_T, \mathbf{p}_T^c) + [1 - P_{\text{rec}}(\mathbf{p}_T^c)] \mathcal{D}_{c \to H}, \quad (47) $$

where $dN_H$ encodes the physics of the recombination. This suggests to organize a model comparison in terms of the three ingredients $P_{\text{rec}}$, $dN_H$, and $\mathcal{D}_{c \to 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. VA we first collect the inclusive production fractions of the ground-state hadrons $D$, $D_s$, and $\Lambda_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. VB we scrutinize the $p_T$ dependence of charm-hadron production using the hadronization ratio, $H_{AA}$ [26] and variants thereof, as well as ratios of different hadron species (hadrochemistry). In Sec. VC we analyze how the $p_T$-dependent elliptic flow of the input charm-quark distribution is converted into that of various hadrons.
<table><caption>TABLE II. Percentage distributions of three different ground-state hadron species to which a c quark hadronizes. Presented are the results for only fragmentation (left three columns), only recombination (middle three columns), and mixed hadronization (right three columns).</caption><thead><tr><th rowspan="2"></th><th colspan="3">Only fragmentation</th><th colspan="3">Only recombination</th><th colspan="3">Mixed hadronization</th></tr><tr><th>D (%)</th><th>D<sub>s</sub> (%)</th><th>Λ<sub>c</sub> (%)</th><th>D (%)</th><th>D<sub>s</sub> (%)</th><th>Λ<sub>c</sub> (%)</th><th>D (%)</th><th>D<sub>s</sub> (%)</th><th>Λ<sub>c</sub> (%)</th></tr></thead><tbody><tr><td>Catania</td><td>78.3</td><td>8.0</td><td>13.7</td><td></td><td></td><td></td><td>48.8</td><td>6.8</td><td>24.3</td></tr><tr><td>Duke</td><td>100</td><td></td><td></td><td>100</td><td></td><td></td><td>100</td><td></td><td></td></tr><tr><td>LBT</td><td>37.8</td><td>5.4</td><td>3</td><td>50.3</td><td>14.6</td><td>20.8</td><td>54.7</td><td>12.1</td><td>15.3</td></tr><tr><td>Nantes</td><td>100</td><td></td><td></td><td>100</td><td></td><td></td><td>100</td><td></td><td></td></tr><tr><td>PHSD</td><td>81</td><td>10</td><td></td><td>67</td><td>33</td><td></td><td>75</td><td>20</td><td></td></tr><tr><td>TAMU</td><td>60.7</td><td>11.5</td><td>24.1</td><td></td><td></td><td></td><td>50.2</td><td>16.2</td><td>22.8</td></tr><tr><td>Torino</td><td></td><td></td><td></td><td></td><td></td><td></td><td>50.6</td><td>17.9</td><td>20.4</td></tr><tr><td>LANL</td><td>77.8</td><td>10</td><td>11.9</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

A. Charm-hadron production fractions

We start our analysis by listing the $p_T$-integrated production fractions of the ground-state hadrons in Table II, i.e., $D$ mesons as well as those carrying an additional conserved quantum number, i.e., the charm-strange $D_s$ meson and the $\Lambda_c$ baryon. We do that for three 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 only-recombination 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 1, indicating that other charm hadrons with additional conserved quantum numbers contribute (such as charm-strange or doubly charm baryons).

Next, we turn to the recombination probability, $P_{\text{rec}}(p_T^c)$, as a function of charm-quark transverse momentum, $p_T^c$, 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-laboratory momentum, $p_T^c$), that a charm quark hadronizes via recombination. In Fig. 5 (top), $P_{\text{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 $\Lambda_c$, $\Xi_c$, $\Omega_c$, and all excited states). The color coding is the same as that in Fig. 2. The fragmentation probability $P_{\text{frag}}(p_T^c)$ is given by $1 - P_{\text{rec}}(p_T^c)$.

There is no quantitative guidance from QCD on how $P_{\text{rec}}$ should behave as a function of $p_T^c$. Only at large $p_T^c$, the domain of Dokshitzer-Gribov-Lipatov-Altarelli-Parisi jet evolution, hadronization is dominated by universal fragmentation functions and therefore $P_{\text{rec}} \to 0$. This lack of guidance is a main reason why the $p_T^c$ dependence of $P_{\text{rec}}$ is rather different in the different approaches. In some approaches fragmentation also contributes at very small $p_T^c$, whereas in others hadronization at small $p_T^c$ is exclusively due to recombination. The latter is based on the assumption that the probability goes to 1 when $p_T^c \to 0$ (usually applied in the thermal rest frame), whereas in the former case the analytical recombination equations are simply extended to $p_T^c = 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 laboratory frame which can cause issues with Lorentz invariance. In fact, also fragmentation functions, usually evaluated in the laboratory frame, are not guaranteed to be Lorentz invariant. While the $p_T^c$ dependence for relative contributions from recombination and fragmentation covers a large variation in the models, we nevertheless observe

<figure><img src="image_6.png" /><figcaption>FIG. 5. The recombination probability of a charm quark of transverse momentum $p_T^c$ in the laboratory frame, into all charmed mesons (upper panel) and into all charmed hadrons (lower panel).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 6. The normalized spectrum, $dN_D/dp_T$, of direct $D^0$'s resulting from c quarks with fixed initial transverse momentum $p_T^c=3.0$ and 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.</figcaption></figure>

three 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_T^c$ (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_T^c$ dependence and show rather good agreement.

## B. Transverse-momentum distributions of charm hadrons

We first carry out task 4 by studying what the transverse-momentum ($p_T^D$) distribution of D mesons looks like in the different models if one hadronizes a c quark of a given transverse momentum, $p_T^c$. If created by a fragmentation process, the D meson carries only a fraction of the $p_T^c$ 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_T^D$, especially for ground-state hadrons. The distributions of $p_T^D$ for a charm quark with $p_T^c = 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, 12—requiring the charm- and light-quark masses in the calculations to be $m_c = 1.5$ GeV, $m_{u/d} = 0.3$ GeV, and $m_s = 0.4$ GeV. In addition, we also impose $\sigma = 0.5$ fm (charmed mesons) and $\sigma_\rho = \sigma_\lambda = 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_T^D$, as compared to $p_T^c$, with a width of the order of 1 GeV for $p_T^c = 3$ GeV, which increases with $p_T^c$ (see 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_T^D < p_T^c + p_T^g$. If we weight fragmentation and recombination according to Fig. 5, the form and the maximum of the distributions become rather different, even for $p_T^c$ as low as 3 GeV (see the right panel of Fig. 6). For $p_T^c = 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_T^c = 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 $\Lambda_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 feed-down contributions. As expected, fragmentation, where the heavy hadron has a smaller $p_T$ than the heavy parton, yields $H_{AA}$ values well below 1 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 to the heavy quark, is smaller and therefore the $H_{AA}$ is larger. The uncertainties in the fragmentation functions cause an uncertainty in the high-$p_T$ $H_{AA}$ of D mesons of almost a factor of 2, and similar for the $D_s$. For the $\Lambda_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 $\approx 2-3$, primarily due to the large number of excited resonances which have been constrained by a fit to pp data [55]. 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$.
<figure><img src="image_8.png" /><figcaption>FIG. 7. The $H_{AA}$ of $D, D_s$, and $\Lambda_c$. For each charm hadron, the only-fragmentation, only-coalescence, and mixed-hadronization processes are plotted.</figcaption></figure>

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 feed-down “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$ [14]. A similar, albeit stronger, trend is observed for $\Lambda_c$ production. However, above $p_T^c = 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 recom-

bination 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 leveling off in the total $H_{AA}$, or even a slight maximum structure in the Nantes model as a consequence of the stronger effect in the recombination 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_T^c = 0$ hadronizes exclusively via recombination, $D$ mesons with $p_T^D = 0$ can also originate from fragmentation of $c$ quarks with a finite $p_T^c$. Therefore, at $p_T^D = 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
<figure><img src="image_9.png" /><figcaption>FIG. 8. The rescaled $H_{AA}^r$ of $D$, $D_s$, and $\Lambda_c$. For each charmed hadron, the only-fragmentation, only-coalescence, and mixed-hadronization processes are plotted.</figcaption></figure>

for the Catania model which falls short of the other models by a factor of 2–3, as was borne out already in the inclusive fractions listed in Table II. For $\Lambda_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%, see Table II). At higher $p_T$ the TAMU model sticks out due to the large feed-down contribution from excited baryons.

The model differences in the charm hadrochemistry 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 a rescaled ratio, $H_{AA}^r \equiv H_{AA}/R$ with

$$R = \frac{\int(dN_c/dp_T)H_{AA}dp_T}{\int(dN_c/dp_T)dp_T} \qquad (48)$$

for coalescence-only and mixed-hadronization processes, and

$$R = \frac{\int_{10}^{20} (dN_c/dp_T) H_{AA} dp_T}{\int_{10}^{20} (dN_c/dp_T) dp_T} \qquad (49)$$

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_{AA}^r$'s (left panels) at high $p_T$ now rather closely agree for $p_T \gtrsim 8$ GeV among all models for all three hadron types. Even at low $p_T \lesssim 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_{AA}^r$, 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 $\Lambda_c$, where two classes seem to emerge, i.e., Catania and Torino vs PHSD, TAMU, and LBT. Somewhat surprisingly, the high-$p_T$ $H_{AA}^r$ for the mixed case shows a similar convergence between the models for all hadron species (with the exception of the $\Lambda_c$ from TAMU) as seen in the fragmentation-only case (left column), despite a rather different normalization.
<figure><img src="image_10.png" /><figcaption>FIG. 9. The yield ratio of $D_s/D$ and $\Lambda_c/D$. For each charm hadron, the only-fragmentation, only-coalescence, and mixed-hadronization processes are plotted.</figcaption></figure>

Next, we study the recombination and fragmentation prob-
abilities for different charmed hadrons, by plotting in Fig. 9
the yield ratios $D_s/D$ (top row) and $\Lambda_c/D$ (bottom row) as
a function of $p_T$. Here $D$ means $D^0$ plus $D^+$. For fragmen-
tation only (left panels) we see an almost constant ratio for
all approaches, indicating that the forms of the fragmentation
functions are very similar. The absolute values differ, how-
ever, slightly for $D_s$, while for $\Lambda_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 feed-down, whereas
for PHSD a falling trend is found characteristic of recombi-
nation with little or no feed-down. For the mixed scenario
 fragmentation plus recombination), with the weights plot-
ted in Fig. 5, all $D_s/D$ results show an increase over the
only-fragmentation 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 neutral-
ization mechanism (Torino), followed by TAMU (including
resonance feed-down), and less for Catania, LBT, and PHSD
which use Wigner functions. For the $\Lambda_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 trans-
lates 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
feed-down of heavy charm-baryon resonances present in their
fragmentation function.

To eliminate the feed-down 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 $\Lambda_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 feed-down contribu-
tions 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
<figure><img src="image_11.png" /><figcaption>FIG. 10. The $H_{AA}$ of the direct $D^0$, $D_s$, and $\Lambda_c$. For each charmed hadron, the only-coalescence and mixed-hadronization processes are plotted.</figcaption></figure>

hard process and therefore commonly assumed to carry van-
ishing $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 two contributions are
not necessarily independent [64]). It is therefore desirable to
assess these two contributions individually and to understand

how the hadronization process modifies the $v_2$. Here we fo-
cus on how the $v_2$ of the heavy meson is generated during
hadronization. However, since no straightforward decompo-
sition 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 $\Lambda_c$ (bottom row) including feed-down contributions, for the only-fragmentation scenario (left column) and the mixed-hadronization scenario
<figure><img src="image_12.png" /><figcaption>FIG. 11. The elliptic flow, $v_2$, of D, $D_s$, and $\Lambda_c$ (including feed-down contributions) for a c quark distribution given by Eq. (46). For each charm hadron, the only-fragmentation (left panels) and mixed-hadronization (right panels) scenarios are displayed.</figcaption></figure>

(right column). The color coding is the same as in Fig. 2, supplemented with a dashed line to benchmark the $v_2$ of the input charm-quark distribution. In almost all models, the only-fragmentation $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(\Lambda_c) > v_2(D_s) > v_2(D)$ for basically all models. This hierarchy is mainly caused by two factors.

(a) The $v_2$ of $\Lambda_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 [65].

Passing from $D \to D_s \to \Lambda_c$, the agreement between the models becomes worse. For the $\Lambda_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
<figure><img src="image_13.png" /><figcaption>FIG. 12. The elliptic flow, $v_2$, of directly produced $D^0$ (left panel), $D_s$ (right panel), and $\Lambda_c$ (right panel) using the c quark distributions defined in Eq. (46). For each charmed hadron, the results of the mixed hadronization process are plotted.</figcaption></figure>

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 feed-down 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 feed-down 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 [18,20]. As model studies show, these correlations may enhance $v_2^D$ 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 [18,20].

## 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 momentum 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.

(1) The recombination probability, $P_{\text{rec}}(p_T^c)$, and hence the fragmentation probability $1 - P_{\text{rec}}(p_T^c)$, 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_T^c$ reach of recombination substantially, which is especially important for the intermediate $p_T^c$ region where recombination and fragmentation mechanisms compete. Finally, different normalization conditions are applied for the low-momentum limit of recombination.

(2) The shift in $p_T$ during the hadronization process, defined as $p_T^D - p_T^c$, 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)].

(3) The functional form of the hadronization ratio, $H_{AA}(i) = \frac{dN_i/dp_T}{dN_c/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_{\text{rec}}(p_T^c)$ 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 focusing on direct production (i.e., no feed-down), where different shapes of the $H_{AA}$ at low $p_T$ exhibit the varying kinematics in the recombination models.

(4) 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 $\Lambda_c$, although the relative enhancement differs somewhat, especially for the $\Lambda_c$, with the fragmentation fraction (and its feed-down contributions) playing an important role. For this observable the decay feed-down 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 feed-down. 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 process.

## ACKNOWLEDGMENTS

We acknowledge support by Deutsche Forschungsgemeinschaft Grant No. CRC-TR 211, “Strong-interaction matter

under extreme conditions,” Project No. 315477589-TRR 211. This work is supported by the European Union’s Horizon 2020 research and innovation program under Grant No. 824093 (STRONG-2020), and by the U.S. Department of Energy (DOE), Office of Science, Office of Nuclear Physics through the Topical Collaboration in Nuclear Theory on Heavy-Flavor Theory for QCD Matter under Grant No. DE-SC0023547. R.R. has been supported by NSF Grants No. PHY-1913286 and No. PHY-2209335. S.A.B. and W.F. acknowledge support by NSF within the framework of the JETSCAPE collaboration, under Grants No. OAC-2004571 (CSSI:X-SCAPE) and No. ACI-1550225, as well as by DOE Grant No. DE-FG02-05ER41367. V.G., V.M., and S.P. acknowledge the support of a linea 2 UniCT Grant for the HQCDyn project. S.C. is supported by National Natural Science Foundation of China Grants No. 12175122 and No. 2021-867. I.V. is supported in part by the Laboratory Directed Research and Development program at LANL

[1] C. Bernard, T. Burch, C. DeTar, J. Osborn, S. Gottlieb, E. B. Gregory, D. Toussaint, U. M. Heller, and R. Sugar, Phys. Rev. D **71**, 034504 (2005).

[2] M. Cheng et al., Phys. Rev. D **74**, 054507 (2006).

[3] R. Rapp and H. van Hees, in *Heavy Quarks in the Quark-Gluon Plasma*, Quark-Gluon Plasma 4 (World Scientific, Singapore, 2010), pp. 111–206.

[4] E. Braaten, K. Cheung, S. Fleming, and T. C. Yuan, Phys. Rev. D **51**, 4819 (1995).

[5] J. Adam et al. (STAR Collaboration), Phys. Rev. Lett. **124**, 172301 (2020).

[6] S. Acharya et al. (ALICE Collaboration), J. High Energy Phys. **01** (2022) 174.

[7] S. Acharya et al. (ALICE Collaboration), Phys. Lett. B **839**, 137796 (2023).

[8] H. van Hees, V. Greco, and R. Rapp, Phys. Rev. C **73**, 034913 (2006).

[9] M. He, R. J. Fries, and R. Rapp, Phys. Rev. C **86**, 014903 (2012).

[10] V. Minissale, S. Plumari, and V. Greco, Phys. Lett. B **821**, 136622 (2021).

[11] S. Cao, G.-Y. Qin, and S. A. Bass, Phys. Rev. C **92**, 024907 (2015).

[12] S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C **94**, 014909 (2016).

[13] S. Cao et al., Phys. Rev. C **99**, 054907 (2019).

[14] S. Cao, K.-J. Sun, S.-Q. Li, S. Y. F. Liu, W.-J. Xing, G.-Y. Qin, and C. M. Ko, Phys. Lett. B **807**, 135561 (2020).

[15] P. B. Gossiaux, R. Bierkandt, and J. Aichelin, Phys. Rev. C **79**, 044906 (2009).

[16] T. Song, H. Berrehrah, D. Cabrera, J. M. Torres-Rincon, L. Tolos, W. Cassing, and E. Bratkovskaya, Phys. Rev. C **92**, 014910 (2015).

[17] T. Song, H. Berrehrah, D. Cabrera, W. Cassing, and E. Bratkovskaya, Phys. Rev. C **93**, 034906 (2016).

[18] M. He and R. Rapp, Phys. Rev. Lett. **124**, 042301 (2020).

[19] H. T. Li, Z. L. Liu, and I. Vitev, Phys. Lett. B **816**, 136261 (2021).

[20] A. Beraudo, A. De Pace, M. Monteno, M. Nardi, and F. Prino, Eur. Phys. J. C **82**, 607 (2022).

[21] L. Adamczyk et al. (STAR Collaboration), Phys. Rev. Lett. **118**, 212301 (2017).

[22] J. Adam et al. (STAR Collaboration), Phys. Rev. C **99**, 034908 (2019).

[23] S. Acharya et al. (ALICE Collaboration), Phys. Lett. B **813**, 136054 (2021).

[24] A. M. Sirunyan et al. (CMS Collaboration), Phys. Lett. B **816**, 136253 (2021).

[25] A. M. Sirunyan et al. (CMS Collaboration), Phys. Lett. B **782**, 474 (2018).

[26] R. Rapp et al., Nucl. Phys. A **979**, 21 (2018).

[27] Y. Xu et al., Phys. Rev. C **99**, 014902 (2019).

[28] F. Cooper and G. Frye, Phys. Rev. D **10**, 186 (1974).

[29] J. F. Owens, Rev. Mod. Phys. **59**, 465 (1987).

[30] A. Andronic et al., Eur. Phys. J. C **76**, 107 (2016).

[31] C. Peterson, D. Schlatter, I. Schmitt, and P. M. Zerwas, Phys. Rev. D **27**, 105 (1983).

[32] M. Cacciari, P. Nason, and R. Vogt, Phys. Rev. Lett. **95**, 122001 (2005).

[33] T. Sjöstrand, S. Mrenna, and P. Z. Skands, J. High Energy Phys. **05** (2006) 026.

[34] S. S. Adler et al. (PHENIX Collaboration), Phys. Rev. Lett. **91**, 172301 (2003).

[35] J. Adams et al. (STAR Collaboration), Phys. Rev. Lett. **98**, 062301 (2007).

[36] J. Adam et al. (ALICE Collaboration), Phys. Rev. C **93**, 034913 (2016).

[37] J. Adams et al. (STAR Collaboration), Phys. Rev. Lett. **92**, 052302 (2004).

[38] A. Adare et al. (PHENIX Collaboration), Phys. Rev. Lett. **98**, 162301 (2007).

[39] B. B. Abelev et al. (ALICE Collaboration), J. High Energy Phys. **06** (2015) 190.

[40] R. C. Hwa and C. B. Yang, Phys. Rev. C **67**, 034902 (2003).

[41] V. Greco, C. M. Ko, and P. Levai, Phys. Rev. Lett. **90**, 202302 (2003).
[42] R. J. Fries, B. Muller, C. Nonaka, and S. A. Bass, Phys. Rev. C 68, 044902 (2003).

[43] D. Molnar and S. A. Voloshin, Phys. Rev. Lett. 91, 092301 (2003).

[44] R. J. Fries, V. Greco, and P. Sorensen, Annu. Rev. Nucl. Part. Sci. 58, 177 (2008).

[45] C.-W. Hwang, Eur. Phys. J. C 23, 585 (2002).

[46] C. Albertus, J. E. Amaro, E. Hernandez, and J. Nieves, Nucl. Phys. A 740, 333 (2004).

[47] J. Zhao, K. Zhou, S. Chen, and P. Zhuang, Prog. Part. Nucl. Phys. 114, 103801 (2020).

[48] L. Ravagli and R. Rapp, Phys. Lett. B 655, 126 (2007).

[49] S. Plumari, V. Minissale, S. K. Das, G. Coci, and V. Greco, Eur. Phys. J. C 78, 348 (2018).

[50] S. Cao, G.-Y. Qin, and S. A. Bass, Phys. Rev. C 88, 044907 (2013).

[51] C. B. Dover, U. W. Heinz, E. Schnedermann, and J. Zimanyi, Phys. Rev. C 44, 1636 (1991).

[52] H. Berrehrah, E. Bratkovskaya, W. Cassing, P. B. Gossiaux, J. Aichelin, and M. Bleicher, Phys. Rev. C 89, 054901 (2014).

[53] P. Moreau, O. Soloveva, L. Oliva, T. Song, W. Cassing, and E. Bratkovskaya, Phys. Rev. C 100, 014911 (2019).

[54] T. Song and G. Coci, Nucl. Phys. A 1028, 122539 (2022).

[55] M. He and R. Rapp, Phys. Lett. B 795, 117 (2019).

[56] D. Ebert, R. N. Faustov, and V. O. Galkin, Phys. Rev. D 84, 014025 (2011).

[57] A. Bazavov et al., Phys. Lett. B 737, 210 (2014).

[58] M. Padmanath, R. G. Edwards, N. Mathur, and M. J. Peardon, PoS LATTICE2014, 084 (2015).

[59] A. Beraudo, A. De Pace, D. Pablos, F. Prino, M. Monteno, and M. Nardi, Phys. Rev. D 109, L011501 (2024).

[60] D. P. Anderle, T. Kaufmann, M. Stratmann, F. Ringer, and I. Vitev, Phys. Rev. D 96, 034028 (2017).

[61] Z.-B. Kang, F. Ringer, and I. Vitev, J. High Energy Phys. 03 (2017) 146.

[62] I. Vitev and M. Gyulassy, Phys. Rev. Lett. 89, 252301 (2002).

[63] I. Vitev, Phys. Lett. B 639, 38 (2006).

[64] M. He, H. van Hees, and R. Rapp, Prog. Part. Nucl. Phys. 130, 104020 (2023).

[65] M. He, R. J. Fries, and R. Rapp, Phys. Rev. Lett. 110, 112301 (2013).