# Heavy flavour hadron production in relativistic heavy ion collisions at RHIC and LHC in EPOS4HQ

Jiaxing Zhao, Joerg Aichelin, Pol-Bernard Gossiaux, Vitalii Ozvenchuk, Klaus Werner

► To cite this version:

Jiaxing Zhao, Joerg Aichelin, Pol-Bernard Gossiaux, Vitalii Ozvenchuk, Klaus Werner. Heavy flavour hadron production in relativistic heavy ion collisions at RHIC and LHC in EPOS4HQ. Physical Review C, 2024, 110 (2), pp.024909. 10.1103/PhysRevC.110.024909. hal-04438004

HAL Id: hal-04438004

https://hal.science/hal-04438004v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Heavy-flavor hadron production in relativistic heavy ion collisions at energies available at BNL
RHIC and at the CERN LHC in the EPOS4HQ framework

Jiaxing Zhao ●,¹ Joerg Aichelin,¹ Pol Bernard Gossiaux ●,¹ Vitalii Ozvenchuk,² and Klaus Werner ●¹

¹SUBATECH, Nantes University, IMT Atlantique, IN2P3/CNRS, 4 rue Alfred Kastler, 44307 Nantes cedex 3, France

²Institute of Nuclear Physics, Polish Academy of Sciences, PL-02-668 Cracow, Poland

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

(Received 4 March 2024; accepted 22 July 2024; published 19 August 2024)

Employing the recently developed EPOS4HQ event generator, we study the production of different heavy-flavor mesons in relativistic heavy-ion collisions at BNL Relativistic Heavy Ion Collider and CERN Large Hadron Collider energies. The transverse momentum spectra, yield ratio, nuclear modification factor, and elliptic flow can be well described in the EPOS4HQ framework. We furthermore analyze the processes which modify these observables as compared to pp collisions and are at the origin of the experimentally determined nuclear modification factor $R_{AA}$.

DOI: 10.1103/PhysRevC.110.024909

I. INTRODUCTION

Theory predicts that matter at high temperature and/or high density forms a plasma of quarks and gluons (QGP) in which partons are not confined anymore. One of the main goals of heavy-ion reactions at ultrarelativistic energies is to study this new form of deconfined matter. In the last years ample evidence that this state exists has been gathered. The observation of strangeness enhancement [1], flow harmonics compatible with hadronization after an hydrodynamic expansions [2], jet quenching [3], and quarkonium suppression [4] points in this direction.

Presently the main objective is to get a more quantitative understanding and investigate the properties of the QGP matter.

The agreement of observables like spectra and flow harmonics with viscous hydrodynamical calculations reveals the strongly coupled nature of the QGP. The very small viscosity $\eta/s$ [5], which is necessary to bring the calculations in agreement with data, is close to the lowest possible limit given by anti-de Sitter/conformal field theory correspondence (AdS/CFT) [6]. So we have strong evidence that a system has been created that behaves as an almost ideal fluid.

Despite the fact that we have a good qualitative understanding about the space-time evolution of the system, there are many open questions when it comes to a quantitative analysis. The above-mentioned probes (yields, spectra, flow harmonics) give only an indirect access to the properties of an expanding system, one has to rely on model calculations or simulations—which have uncertainties. So we need additional

probes to get complementary information about the system, done in the same model and/or simulation. Heavy flavor hadrons have turned out to be an almost ideal probe to study the time evolution of the QGP due to the following reasons: (a) The heavy quark mass ($m_c$ = 1.5 GeV and $m_b$ = 4.5 GeV) is much larger than the quantum chromodynamics (QCD) cutoff, $\Lambda_{QCD}$ ≈ 200 MeV. Therefore their production can be well described by perturbative QCD (pQCD). (b) Heavy quarks are produced at the early stage of heavy ion collisions and witness all later stages of the collision. (c) The heavy quark mass is much larger than the typical temperature of the QGP medium, which is about a couple of hundred of MeV, estimated by the spectrum of the directly produced photons [7]. Consequently, the mass of heavy quark changes little in the hot medium and their number is conserved during the evolution.

With these advantages, heavy flavor physics has attracted a lot of attention from both the theoretical and experimental communities and several models have been advanced to describe the heavy flavor observables. The Parton-Hadron-String Dynamics (PHSD) [8,9] approach, in which the heavy quark physics is embedded in an approach which described the light hadrons as well, is based in the dynamical quasi-particle model (DQPM), which respects the equation of state of strongly interacting matter. The Catania model [10–12] is also based on a DQPM approach and describes the expanding medium by a Boltzmann equation. Based on the EPOS3 event generator we have also studied heavy quark production [13] using an elastic and inelastic pQCD cross sections. The Linear Boltzmann Transport (LBT) model [14] solves as well a Boltzmann equation for the heavy quarks including elastic and inelastic collisions but the medium is described by viscous hydrodynamics. Other models like TAMU [15,16], Duke [17], and Torino [18,19] use a Fokker-Planck equation to describe the dynamics of heavy quarks and model the expanding QGP by ideal or viscous hydrodynamics. Apart from PHSD and EPOS, all these models have in common that they concentrate on heavy quark physics only and do not take advantage of

<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>
the fact that the many light hadron observables allow for assessing the expansion of the QGP despite different expansion properties of the QGP possibly a strong influence on the heavy quark observables [20].

Thanks to the high statistic data, which are now or will soon become available, the error of the key observables, such as the enhancement of the baryon to meson ratio and the elliptic flow, is strongly reduced. This allows for a detailed quantitative comparison between theory and experiment. In this paper we compare the heavy hadron observables with the results of EPOS4HQ. EPOS4HQ is the heavy hadron extension of the recently advanced EPOS4 approach [21–24], which has been successfully used to study the light hadron observables. In the EPOS4HQ approach the heavy flavor production has been substantially improved in comparison with the former EPOSLHC approach. Heavy flavor quarks can now be produced in a hard process, as well as by gluon splitting and flavor excitation. The interaction of heavy partons with the QGP includes elastic and gluon emission reactions as in Refs. [25,26]. After the hadronization the heavy hadrons still have final state interactions, modeled by ultrarelativistic quantum molecular dynamics (UrQMD). The calculations are based on version EPOS4.0.1.s9.

In this paper we will systematically investigate all published heavy flavor observables. It is organized as follows. The heavy quark initial production and the medium evolution in EPOS4HQ are discussed in Sec. II. Next, the heavy quark initial spectra and energy loss mechanism is the subject of Secs. III and IV. In Section V, we will present the hadronization in EPOS4HQ, which differs substantially from former heavy quark calculations with EPOS2/EPOS3/EPOSLHC by including baryons and excited mesons states. The results and the comparison with experimental data are shown in Sec. VI. A conclusion will be given in Sec. VII.

## II. EPOS4

### A. EPOS4 primary interactions

A fundamental ingredient of the EPOS4 approach [21–24] is the observation that multiple partonic scatterings must strictly happen in parallel, and not sequentially, based on very elementary considerations concerning time-scales. To take this into account, EPOS4 brings together ancient knowledge about S-matrix theory (to deal with parallel scatterings) and modern concepts of perturbative QCD and saturation, going much beyond the usual factorization approach. The parallel scattering principle requires sophisticated Monte Carlo techniques, inspired by those used in statistical physics to investigate the Ising model.

In the EPOS4 approach, we distinguish “primary scatterings” and “secondary scatterings”. The former refer to the above-mentioned parallel scatterings with the initial nucleons (and their partonic constituents) being involved, happening at very high energies instantaneously. The theoretical tool is S-matrix theory, using a particular form of the proton-proton scattering S matrix (“classical” Gribov-Regge approach [27–30]). Within such an approach, one can deduce the very important AGK theorem [30], which leads to factorization and

<figure><img src="image_3.png" /><figcaption>FIG. 1. Rigorous parallel scattering scenario, for $n = 3$ parallel scatterings, including nonlinear effects via saturation scales. The red symbols should remind us that the parts of the diagram representing nonlinear effects are replaced by simply using saturation scales.</figcaption></figure>

binary scaling in nuclear scatterings, which is not trivial in a multiple scattering scheme. However, introducing energy-momentum sharing [31] (which is absolutely crucial for realistic event-by-event simulations, AGK is violated (and so is factorization and binary scaling). The main new development in EPOS4 [21–24] is a way to accommodate simultaneously: (1) rigorous parallel scattering, (2) energy-momentum sharing, (3) AGK theorem and factorization for hard processes, by introducing (in a very particular way) saturation, compatible with recent “low-x-physics” considerations [32–43].

Validity of AGK means that we can do the same as models based on factorization (defining $n$ and using parton distribution functions) to study very hard processes, but this represents only a very small fraction of all possible applications, and there are very interesting cases outside the applicability of that approach. A prominent example, one of the highlights of the past decade in our domain, concerns collective phenomena in small systems. It has been shown that high-multiplicity $pp$ events show very similar collective features as earlier observed in heavy ion collisions [44]. High multiplicity means automatically “multiple parton scattering”. As discussed earlier, this means that we have to employ the full parallel scattering machinery developed earlier, based on S-matrix theory. We cannot use the usual parton distribution functions (representing the partonic structure of a fast nucleon), we have to treat the different scatterings (happening in parallel) individually, for each one we have a parton evolution according to some evolution function $E$ (representing the partonic structure of a fast parton), as sketched in Fig. 1.

We still have Dokshitzer–Gribov–Lipatov–Altarelli–Parisii (DGLAP) evolution, for each of the scatterings [we only show the spacelike (SL) cascade], but we introduce saturation scales. But, most importantly, these scales are not constants, they depend on the number of scatterings, and they depend as well on $x^+$ and $x^-$. An example of a multiple scattering AA configuration is shown in Fig. 2.

### B. EPOS4 heavy quark issues

At each step in the SL cascade, there is the possibility of quark-antiquark production, and in the Born process as well. In the following, we discuss in particular the case of heavy
<figure><img src="image_4.png" /><figcaption>FIG. 2. Rigorous parallel scattering scenario, for $n = 3$ parallel scatterings for a collision of a nucleus A with a nucleus B, including nonlinear effects via saturation scales.</figcaption></figure>

flavor quarks, with the general notation $Q$ for quarks and $\bar{Q}$ for antiquarks (for details see Ref. [22]). Heavy flavor may be produced in different ways, as shown in Fig. 3. Starting from a gluon, a $Q-\bar{Q}$ pair may be produced in the SL cascade, as shown in Fig. 3(a), provided the virtuality is large enough. The number of allowed flavors is considered to be depending on the virtuality (variable flavor number scheme). It is also possible to create a $Q-\bar{Q}$ in the Born process, via $g+g \to Q+\bar{Q}$ or $q+\bar{q} \to Q+\bar{Q}$ (for light flavor quarks $q$), as shown in Fig. 3(b), and finally $Q-\bar{Q}$ may be produced in the time-like cascade (TLC), via $g \to Q+\bar{Q}$, as shown in Fig. 3(c). Let us first consider the $Q-\bar{Q}$ production in the space-like cascade (SLC). We may have the situation as shown in Fig. 4(a), where a heavy flavor parton (here a $\bar{Q}$) is emitted, and the corresponding antiparticle (here a $Q$) continues the SLC. But before reaching the Born process, it is emitted, and a gluon continues the SLC. The two heavy flavor partons have in general low transverse momenta. Another possibility is shown in Fig. 4(b), where a heavy flavor parton produced in the SLC “survives” till the Born process, and the latter has most likely the form $Q+l \to Q+l$, with $l$ being a light flavor parton. Other than the production during the SLC, heavy flavor may be produced in the Born process, via $g+g \to Q+\bar{Q}$ or $q+\bar{q} \to Q+\bar{Q}$ (for light flavor quarks $q$), as shown in Fig. 4(c). Finally, heavy flavor may be produced during the time-like cascade, as shown in Fig. 5, either initiated from a TLC parton in the SLC [Fig. 5(a)], or initiated from an outgoing parton of the Born process [Fig. 5(b)]. In the first case, the transverse momenta are in general small.

The next step will be, for a given Feynman diagram, to construct the color flow diagram. Let us take the graph of

<figure><img src="image_5.png" /><figcaption>FIG. 3. Different possibilities to create heavy flavor, (a) in the space-like cascade (SLC), (b) in the Born process, (c) in the time-like cascade (TLC).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4. Heavy flavor production (a,b) in the SLC and (c) in the Born process. The magenta point indicates the Born process.</figcaption></figure>

Fig. 5(b), i.e., heavy flavor production during the TLC of an outgoing Born parton. As usual, the gluons are emitted to either side with equal probability, so a possible color flow diagram is the one shown in Fig. 6. We identify three chains of partons: 1-2-3-4-5, 6-7-8, and 9-10-11. The initial TL partons (the horizontal blue lines with arrows) or most likely quarks and antiquarks (in any case 3 and $\bar{3}$ color representations). Let us assume that 3 is a quark, and 6 an antiquark (light flavor, both), then the two chains containing heavy flavor are of the form $\bar{Q}-g-q$ and $\bar{q}-g-Q$, in both cases, the heavy flavor partons are “end partons” in the chains.

These chains of partons are finally mapped (in a unique fashion) to kinky strings, where each parton corresponds to a kink, as shown in Fig. 7. The general mapping procedure (chains of partons to kinky strings) as well as the string decay procedures (producing so-called “prehadrons”) are described in detail in [31].

### C. EPOS4 core-corona method and fluid evolution

The above-mentioned parallel scattering happens at zero time. After that, we obtain a more or less important number of

<figure><img src="image_7.png" /><figcaption>FIG. 5. Heavy flavor production in the TLC.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. A possible color flow diagram corresponding to the graph of Fig. 5(b).</figcaption></figure>

prehadrons. We employ a core-corona procedure [24,45–47], where the prehadrons, considered at a given proper time $\tau_0$, are separated into “core” and “corona” prehadrons, depending on the energy loss of each prehadron when traversing the “matter” composed of all the others. Corona prehadrons (per definition) can escape, whereas core prehadrons lose all their energy and constitute what we call “core”, which acts as an initial condition for a hydrodynamic evolution [47,48].

The evolution of the core ends whenever the energy density falls below some critical value $\epsilon_{FO}$, which marks the point where the fluid “decays” into hadrons. It is not a switch from fluid to particles, it is a sudden decay, called “hadronization”. Let us consider a (randomly chosen, but typical) 5.02 ATeV lead-lead scattering. In Fig. 8, we plot the energy density in the transverse plane (x, y). We consider two snapshots, namely at the start time of the hydroevolution $\tau_0$ (upper plot) and a later time $\tau_1$ close to final freeze-out (lower plot). The initial distribution has an elongated shape (just by accident, due to the random positions of interacting partons). One can clearly see that the final distributions are as well elongated, but perpendicular to the initial ones, as expected in a hydro-dynamical expansion due to the larger pressure perpendicular

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

<figure><img src="image_10.png" /><figcaption>FIG. 8. Energy density in the transverse plane (x, y) for a Pb-Pb collision at $\sqrt{s_{NN}}$ = 5.02 TeV with an impact parameter of 10.4 fm. The upper plot represents the start time $\tau_0$ (of the hydroevolution), and the lower plot a later time $\tau_1$, close to final freeze-out.</figcaption></figure>

to the elongation, which creates a strong collective (elliptic) flow. More examples can be found in [24].

In EPOS4, as discussed in detail in [24], we developed a new procedure of energy-momentum flow through the “freeze-out (FO) hypersurface” defined by $\epsilon_{FO}$, which allows defining an effective invariant mass, decaying according to microcanonical phase space into hadrons, which are then Lorentz boosted according to the flow velocities computed at the FO hypersurface. We also developed new and very efficient methods for the microcanonical procedure [24]. Also in the full scheme, including primary and secondary interactions, energy-momentum and flavors are conserved.

### III. INITIAL CHARM QUARK MOMENTUM DISTRIBUTION

EPOS2 and EPOS3 provided only the interaction points at which charm and anticharm quarks have been produced.

<figure><img src="image_11.png" /><figcaption>FIG. 7. The chains 1-2-3-4-5, 6-7-8, and 9-10-11 are mapped to kinky strings (red lines). The black points indicate the kinks, which carry the parton momenta.</figcaption></figure>
Their transverse momentum had been chosen to reproduce the fixed order next to leading log (FONLL) calculations [49,50] and energy and momentum has not been conserved at the interaction points. Correlations between the heavy quark and antiquark, imposed by the pQCD processes, have not been conserved in the standard version.

EPOS4, on the contrary, provides the correlation, imposed by the pQCD processes, between the momenta of the charm and the anticharm quark which are created at the same vertex. These correlations are interesting as they can be used to study the thermalization of heavy quarks in a QGP [51,52]. These pQCD correlations are difficult to compare with other approaches. Only PYTHIA provides this correlation as well. In the standard version of PYTHIA, however, only the leading order cross section is calculated and therefore the opening angle between $c$ and $\bar{c}$ is 180°. Including initial and final state interactions this distribution is modified but a one-to-one correspondence with the diagrams included in EPOS4 is not possible. The initial single (anti)charm quark $p_T$ distribution, given by EPOS4HQ for $pp$ collisions can be compared with the FONLL calculations [49,50], the most advanced pQCD approach for describing single charm quark distributions. In Fig. 9 we display the FONLL and EPOS4 single charm and bottom quark distributions for 5.02 TeV $pp$ collisions. We see that the EPOS4 distribution is at the upper limit of the error bars of the FONLL calculation for low $p_T$ and consistent with FONLL for high $p_T$. This means that EPOS4 gives a larger charm production cross section, which is consistent with the experimental data obtained recently [53]. When produced in heavy-ion collisions, shadowing, saturation effects, and the Cronin effect modify the initial (anti)charm quark distribution. In Fig. 10 we compare for PbPb at $\sqrt{s_{NN}} = 5.02$ TeV the extrapolated $pp~p_T$ distribution (black dashed line) with the PbPb $p_T$ distribution, which takes the above-mentioned cold nuclear matter effects into account (blue line). At large $p_T$ both distributions are almost identical but at low $p_T$ one sees differences due to the different saturation scales in $pp$ and PbPb collisions, which affect the $p_T$ spectra close to this scale.

#### IV. EPOS4HQ: HEAVY QUARK ENERGY LOSS IN HOT MEDIUM

In contradistinction to EPOS4, in EPOS4HQ heavy quarks interact with the partons of the QGP, formed by the light partons and gluons from the core. We include in this study both elastic [25] and radiative [26] collisions. To determine where the interaction takes place, we calculate the interaction rate and move the heavy quark to the interaction point, select whether the QGP parton is a gluon or a quark, and draw the momentum of the scattering QGP parton from its corresponding thermal distribution. The thermal distribution is determined by the local temperature and mean velocity at the position at which the collision takes place. The scattering cross sections of the heavy quark with gluons and light quarks are calculated by pQCD matrix elements with a running coupling constant.

The pQCD elastic scattering cross section diverges for a small momentum transfer in the $t$ and $u$ channels. These

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

<figure><img src="image_13.png" /><figcaption>FIG. 9. The initial transverse momentum distribution of charm (upper) and bottom (lower) quarks for $pp$ collisions at 5.02 TeV. The solid black line is the EPOS4HQ calculation, while the orange (blue) dashed line presents the FONLL [49] result. The uncertainty of the FONLL result is given by the orange (blue) shaded area.</figcaption></figure>

infrared divergences are healed by the Debye screening mass $m_D(T)$ of gluons in the hot medium, which is calculated in the hard thermal loop (HTL) approach. It serves as a regulator of the propagator of the exchanged gluon. Scattering at high momentum transfer is, on the contrary, described by a free gluon propagator for massless gluons. A smooth transition of the energy loss between both regimes can be assured by an effective Debye mass $m_{\text{eff}} = \kappa m_D(T)$ in the gluon propagato, with $\kappa = 0.2$ [25].

The pQCD inelastic scattering cross section has been calculated in [26]. This cross section contains five matrix elements for gluon emission from the heavy quark and the light quark and gluon, respectively. Also for the inelastic cross section the momentum of the plasma particle is chosen by a
<figure><img src="image_14.png" /><figcaption>FIG. 10. The transverse momentum distribution of charm quarks at creation and before hadronization for central PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. The transverse momentum distribution observed in pp collisions, multiplied by the number of initial binary collisions, is shown as a black dashed line.</figcaption></figure>

Monte Carlo approach from the local thermal distribution. As in the elastic cross section, the gluon propagator is regulated by $m_{eff} = \kappa m_D(T)$. For the gluon emission vertex a constant $\alpha_s = 0.3$ is used. The emitted gluon is considered as massless. The different limits of the pQCD cross section calculations as well as more details of the approach have been discussed in Ref. [26].

The validity of this approach can be confirmed by comparing its spatial diffusion coefficient $D_s$ [54] with lattice calculations. Figure 11 shows $2\pi T D_s$ of our model, in com-

<figure><img src="image_15.png" /><figcaption>FIG. 11. $2\pi T D_s$ obtained in our approach in comparison with lattice data from [55]. The full red line corresponds to our collisional+radiative model, while the orange long-dashed line represents the contribution from the sole collisional part; the dashed curves illustrate the equivalent calculation if the IR regulator is taken as the Debye mass ($\kappa = 1$).</figcaption></figure>

parison with recent unquenched lattice data [55]. The former lattice data for $2\pi T D_s$ were calculated in a quenched approximation, so we do not include them here. We see that both coefficients agree well for $\kappa = 0.2$ whereas there are large differences for the choice $\kappa = 1$. Inelastic collisions have only a negligible influence on $D_s$.

Both the elastic as well as the inelastic collisions have been already employed in the EPOS2 and EPOS3 frameworks to describe heavy meson data in heavy-ion collisions at CERN Large Hadron Collider (LHC) energies [13,52]. We use this theoretical framework without modification also in this new EPOS4HQ version. In this paper the K factor for elastic as well as for inelastic collisions, which has been varied in the past Refs. [13,52] in calculations in which the HQ part was coupled to EPOS2 or EPOS3, is equal one, so the calculated pQCD cross sections are not modified by an overall factor.

Figure 10 displays the change of the momentum distribution of heavy quarks caused by the interactions with the QGP for central PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. As said, the blue line is the initial transverse momentum distribution. The $p_T$ distribution of the heavy quark just before they hadronize is represented by the red line. At low $p_T$ there is little difference between the curves because the initial $p_T$ of the heavy quarks is close to the value one expects if the heavy quarks come to a thermal equilibrium with the light quarks, so the momentum change is small. For large momenta, on the contrary, we observe a $p_T$ shift to lower $p_T$ values of the order of 5 GeV.

## V. HADRONIZATION SCHEME

When the QGP is expanded to the critical value of the local energy density $\epsilon_{FO} = 0.57 \text{GeV/fm}^3$), the system hadronizes and the heavy quark converts into a heavy flavor hadron. In the EPOS4HQ framework, there are two ways in which the heavy quark can hadronize, either by fragmentation or by coalescence. In the coalescence process, the heavy quarks coalesce with light quarks from the hypersurface at which they are localized. The fluid at the hypersurface has an average velocity. The momentum of the light quark, with which the parton coalesces, is selected randomly from the thermal distribution with the fluid temperature $T_{FO}$, which is related to the critical energy density $\epsilon_{FO}$ via the equation of state, and the average fluid velocity at the hypersurface.

In this study, we neglect the coordinate information of quarks and require only that the hadronization hypersurface is the same for the heavy and light quarks. In the numerical simulation the coalescence formula is applied in the center-of-mass frame of heavy and light quarks.

The momentum distribution of Heavy Flavor (HF) hadrons, produced via the coalescence process, can be calculated by

$$ \frac{dN}{d^3\mathbf{P}} = g_H \sum_{N_Q} \int \prod_{i=1}^{k} \frac{d^3\mathbf{p}_i}{(2\pi)^3} f(\mathbf{p}_i) \times W_H(\mathbf{p}_1, \dots, \mathbf{p}_i) \delta^{(3)}\left(\mathbf{P} - \sum_{i=1}^{N} \mathbf{p}_i\right), \quad (1) $$

where $g_H$ is the degeneracy factor of color and spin. **P** and **p_i** are the momenta of heavy flavor hadrons and the constituent
quarks, respectively. The $\delta$ function conserves momentum. The integration is carried at the point where the heavy quark crosses the hadronization hypersurface, which is given by EPOS4. The summation is performed over all heavy quarks in the system. $f(\mathbf{p}_i)$ is the momentum space distribution of the constituent in the heavy hadron at the moment of hadronization with $k = 2$ for mesons and 3 for baryons. $W_H(\mathbf{p}_1, \dots, \mathbf{p}_k)$ is the Wigner density of a given heavy hadron $H$, which can be constructed from the heavy hadron wave function, the solution of the Schrödinger equation, and which is approximated here by a three-dimensional harmonic oscillator state with the same root mean square radius. After integration over the coordinate space, the Wigner density of the mesonic ground state can be expressed as

$$W(p_r) = (2\sqrt{\pi}\sigma)^3 e^{-\sigma^2 p_r^2}. \qquad (2)$$

$p_r$ is the relative momentum between the two constituent quarks in the center-of-mass (c.m.) frame, $p_r = |E_2\mathbf{p}_1 - E_1\mathbf{p}_2|/(E_1 + E_2)$, $E_1(\mathbf{p}_1)$ and $E_2(\mathbf{p}_2)$ are the energies (momenta) of the quark and antiquark in the heavy hadron c.m. frame, respectively. Baryons are treated as two two-body systems (baryons are produced by recombining two particles first and then by using their center of mass to recombine this diquark with the third quark). The light quarks (antiquarks) are assumed to be thermalized. In the rest system of the hypersurface their distribution is given by

$$f_q(\mathbf{p}_q) = \frac{1}{(2\pi)^3} \frac{g}{e^{E_q/T_{\text{FO}}} + 1}, \qquad (3)$$

where $g = 6$ is the statistical factor. $E_q = \sqrt{\mathbf{p}_q^2 + m_q^2}$ is the energy of the light quark. The quark masses are $m_{u/d} = 0.1$ GeV, $m_s = 0.3$ GeV.

The heavy quark coalescence probability $P_{\text{coal}}$ in a static hot medium with an energy density $\epsilon_{\text{FO}} = 0.57$ GeV/fm$^3$ for the ground states of the open heavy flavor states, e.g., $D^0, D^+, D_s, \Lambda_c, \Xi_c$, and $\Omega_c$ for charm, $B^0, B^-, B_s, \Lambda_b, \Xi_b$, and $\Omega_b$ for bottom, is obtained by integrating Eq. (1) with a $\delta$ distribution $\delta^3(\mathbf{p}_1 - \mathbf{p}_Q)$ for the momentum of the heavy quark and with a $\delta$ distribution $\delta^3(\mathbf{p}_2 - \mathbf{p}_q)$ for the momentum of the quark. $p_q$ is sampled from the thermal distribution, Eq. (3). We neglect here the rarely produced heavy flavor hadrons, such as multi-charmed (bottomed) baryons, $B_c$, and quarkonium. For heavy flavor mesons, the averaged square radius can be expressed as $\langle r^2 \rangle = \frac{m_Q^2+m_q^2}{2(m_Q+m_q)^2}\sigma^2$ with charm quark mass $m_c = 1.5$ GeV and $m_b = 4.5$ GeV for bottom quarks. The root-mean-square radius of the ground state charmed meson has been calculated by the two-body Dirac equation [56]. It gives $\sqrt{\langle r^2 \rangle} = 0.85$ fm, for in-medium $D^0$. So, the corresponding width $\sigma = 3.725 \text{ GeV}^{-1}$ for $D^0$. In the absence of theoretical studies of the in-medium radius of the charmed baryon, we take the same width $\sigma$ for any two-quark systems in charmed baryons. Because the reduced masses are very similar, the averaged radius of bottom mesons is comparable to that of charmed mesons. So, we take the same width for bottom mesons and any two-quark systems in bottom baryons.

The probability that excited states are created, which can strongly decay into the ground states, is estimated via the statistic model. There the hadron density at the temperature

<table><caption>TABLE I. Number of excited states considered in EPOS4HQ ($N_{\text{excited}}$) and their relative contributions to the ground states, $R_i$.</caption><thead><tr><th></th><th>D<sup>0</sup></th><th>D<sup>+</sup></th><th>D<sub>s</sub></th><th>&Lambda;<sub>c</sub></th><th>&Xi;<sub>c</sub></th><th>&Omega;<sub>c</sub></th></tr></thead><tbody><tr><td>N<sub>excited</sub></td><td>9</td><td>8</td><td>5</td><td>91</td><td>91</td><td>53</td></tr><tr><td>R</td><td>2.91</td><td>0.96</td><td>2.06</td><td>7.63</td><td>3.97</td><td>4.22</td></tr><tr><td></td><td>B<sup>0</sup></td><td>B<sup>-</sup></td><td>B<sub>s</sub></td><td>&Lambda;<sub>b</sub></td><td>&Xi;<sub>b</sub></td><td>&Omega;<sub>b</sub></td></tr><tr><td>N<sub>excited</sub></td><td>6</td><td>6</td><td>2</td><td>91</td><td>91</td><td>53</td></tr><tr><td>R</td><td>2.95</td><td>2.95</td><td>2.26</td><td>7.67</td><td>3.92</td><td>5.99</td></tr></tbody></table>

$T_{\text{FO}}$ is given by [57]

$$n_i = \frac{g_i}{2\pi^2} T_{\text{FO}} m_i^2 K_2\left(\frac{m_i}{T_{\text{FO}}}\right), \qquad (4)$$

where $g_i$ is the spin isospin degeneracy. $m_i$ is the mass of the hadron. $K_2$ is the second-order Bessel function. In our study, we consider almost all possible excited states, also the missing baryons, which are predicted by the quark model [58] and lattice QCD [59,60], as shown in Table I. For each ground state hadron $D, D_s, \Lambda_c, \Xi_c$, and $\Omega_c$ we calculate the density of their excited states $m$ and the ratio $R^m = n_{\text{excited}}^m/n_{\text{ground}}$. This ratio is momentum-independent. Finally we sum up $R = \sum R^m$ and multiply the ground state momentum-dependent coalescence probability, which is calculated via Eq. (1) with the Wigner function as discussed above, by $1 + R$ to obtain its effective momentum distribution with all excited state contributions. The sum of the effective momentum distributions for all hadrons gives the total coalescence probability $P_{\text{coal}}(p)$. Same as for the bottom sector.

The probability that a charm (bottom) quark with a given momentum $\mathbf{p}_c(\mathbf{p}_b)$ forms a specific hadron by coalescence is shown in Fig. 12. Heavy quarks, which do not hadronize via coalescence, will fragment into a heavy-flavor hadron. The fragmentation probability is therefore $1 - P_{\text{coal}}$. In EPOS4HQ, we use the heavy quark effective theory-based fragmentation function [50,61]. The fragmentation ratios to various charmed hadrons are taken as the $e^+e^-$ collisions [62] and are shown in Table II. After the hadronization, all charmed hadrons and light hadrons evolve together in the hadronic phase, which is controlled by the UrQMD [63]. Bottom hadrons do not interact in the hadronic phase.

## VI. RESULTS

In this section, we present the EPOS4HQ results for open heavy flavor hadrons. We are as exhaustive as possible and

<table><caption>TABLE II. Fragmentation ratio of charm (bottom) quark to charmed (bottomed) hadrons in percent.</caption><thead><tr><th>D<sup>0</sup></th><th>D<sup>+</sup></th><th>D<sub>s</sub></th><th>&Lambda;<sub>c</sub></th><th>&Xi;<sub>c</sub></th><th>&Omega;<sub>c</sub></th></tr></thead><tbody><tr><td>60.8%</td><td>24.0%</td><td>8.0%</td><td>6.0%</td><td>1.0%</td><td>0.2%</td></tr><tr><td>B<sup>0</sup></td><td>B<sup>-</sup></td><td>B<sub>s</sub></td><td>&Lambda;<sub>b</sub></td><td>&Xi;<sub>b</sub></td><td>&Omega;<sub>b</sub></td></tr><tr><td>42.4%</td><td>42.4%</td><td>8.0%</td><td>6.0%</td><td>1.0%</td><td>0.2%</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>FIG. 12. The momentum dependent coalescence probabilities of charm (upper) and bottom (lower) quark in a hot and static medium with an energy density of 0.57 GeV/fm$^3$.</figcaption></figure>

show the results for all published experimental data on transverse momentum and elliptic flow. We present as well yield ratios like $R_{AA}$ and interpret these results with the help of additional information, which the EPOS4HQ approach provides.

A. Transverse momentum spectra

1. Charm hadrons

The transverse momentum spectra of $D^0$, $D^+$, $D_s$, $\Lambda_c$, $\Xi_c$, and $\Omega_c$ in PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV and for 0–10%, 30–50%, and 60–80% centrality are shown in Figs. 13 and 14. The EPOS4HQ calculations (thick black line) are compared with the ALICE data: $D^0$ and $D^+$ from [64], $D_s$ form [65,66], and $\Lambda_c$ from [67]. For $D^0$ ($\Lambda_c$) we show as well the initial $p_T$ distribution of those charm quarks which are finally entrained in $D^0$ ($\Lambda_c$) hadrons (dotted blue line)

<figure><img src="image_17.png" /><figcaption>FIG. 13. $p_T$ spectra of charmed hadrons $D^0$, $D^+$, $D_s$, $\Lambda_c$, $\Xi_c^0$, and $\Omega_c$ in the 0–10% centrality class of PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. The experimental data of $D$ and $D^+$ [64], $D_s$ [65,66], $\Lambda_c$ [67] are from the ALICE Collaboration. The thick black line is the EPOS4HQ result, the dotted blue and dashed green lines are the $p_T$ distribution of charm quarks at production and before hadronization, respectively, which finally are part of the $D^0$ ($\Lambda_c$).</figcaption></figure>

as well as that before hadronization (dashed green line) to demonstrate how the passage through the QGP and the subsequent hadronization modifies the initial distribution. At high transverse momenta the $p_T$ distribution before the $c$ quarks hadronizes is strongly suppressed as compared to the initial distribution, testifying to the energy loss of the $c$ quark in the QGP. At low momentum, when the momentum of the $c$ quark is of the order of the averaged momentum of the QGP partons, collisions have the consequence that the $c$ quarks approach an equilibrium with the QGP. Comparing the dashed green line to the black line, we can see clearly the influence of the hadronization on the transverse momentum change. High $p_T$ charm hadronizes almost exclusively via fragmentation, which shifts the spectra to lower $p_T$. Low $p_T$ charm hadronizes by combining with another light antiquark (or two light quarks for baryons), which shifts the momentum to a higher region. This momentum shift is more pronounced for $\Lambda_c$ due to two light quarks entrained. Where data are available, the transverse momentum distributions of EPOS4HQ are for all centrality bins close to the measured ones.

The $p_T$ spectra of different charmed hadrons, produced in AuAu collisions at RHIC at $\sqrt{s_{NN}} = 200$ GeV and in the 0–10% and 10–40% centrality bins, are shown in Fig. 15. The EPOS4HQ results (thick black line) are compared with the STAR data [68] for $D$ mesons and [69] for $D_s$ mesons.
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 14. $p_T$ spectra of charmed hadrons $D^0$, $D^+$, $D_s$, $\Lambda_c$, $\Xi_c^0$, and $\Omega_c$ in the 30–50 % and 60–80 % centrality classes of PbPb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV. The experimental data of D and $D^+$ [65], $D_s$ [65], $\Lambda_c$ [67] are from ALICE.</figcaption></figure>

For the $D^0$ data we show as well the initial $p_T$ distribution of those charm quarks which are finally entrained in a $D^0$ meson (dotted blue line) as well as that before hadronization (dashed

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

<figure><img src="image_21.png" /><figcaption>FIG. 15. $p_T$ spectra of charmed hadrons $D^0$, $D^+$, $D_s$, $\Lambda_c$, $\Xi_c^0$, and $\Omega_c$ in the 0–10 % and 10–40 % centrality classes of AuAu collisions at $\sqrt{s_{NN}}$ = 200 GeV. The experimental data of D [68], $D_s$ [69] are from STAR. The thick black line is the EPOS4HQ result, the dotted blue and dashed green lines are the $p_T$ distribution of charm quarks at production and before hadronization, respectively, which finally entrained in $D^0$ mesons.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>FIG. 16. $p_T$ spectra of charmed hadrons $B^0$, $B^-$, $B_s$, $\Lambda_b$, $\Xi_b$, and $\Omega_b$ in the 0–100 % centrality class of PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. The experimental data of $B^0$ [70] and $B_s$ [71] are from the CMS collaboration. The thick black line is the EPOS4HQ result, the dotted blue and dashed green lines are the $p_T$ distribution of bottom quarks at production and before hadronization, respectively, which finally are entrained in $B^0$ mesons.</figcaption></figure>

green line) to demonstrate how the passage through the QGP and the subsequent hadronization modifies the initial distribution. As compared to the LHC results the energy loss while traversing the QGP shifts the momentum less and therefore the spectrum at high $p_T$ is less suppressed. Also at RHIC we reproduce nicely the available experimental data.

## 2. Bottom hadrons

The transverse momentum spectra of bottom hadrons $B^0$, $B^-$, $B_s$, $\Lambda_b$, $\Xi_b$, and $\Omega_b$ in the 0–100 % centrality class of PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV are shown in Fig. 16. The experimental data of $B^0$ [70] and $B_s$ [71] are from the CMS collaboration. For the $B^0$ data we show as well the initial $p_T$ distribution of the charm quarks (dotted blue line) as well as that before hadronization (dashed green line), which are finally part of $B^0$ mesons. The momentum loss of a bottom quark while passing through the plasma is less than that of a charm quark. This is due to the kinematics of the interaction of the large mass b quarks with the QGP partons, which suppresses the momentum transfer as compared to the interaction of c quarks.

Another source of information about B mesons are the decay products, especially the nonprompt $J/\psi$ s. Their $p_T$

<figure><img src="image_23.png" /><figcaption>FIG. 17. $p_T$ spectrum of $B^0$ and nonprompt $J/\psi$ in the 0–10 % PbPb collisions at 5.02 TeV. The experimental data are from ATLAS [72] and ALICE [73].</figcaption></figure>

spectrum has been measured for central events in PbPb at $\sqrt{s_{NN}} = 5.02$ GeV by the ATLAS [72] and the ALICE [73] collaborations. Their results are displayed in Fig. 17 and compared to that of EPOS4HQ calculations. We display in this figure as well the $p_T$ distribution of the $B^0$ mesons in this centrality bin. In our calculation we assumed that nonprompt $J/\psi$ come exclusively from $B^0$ decays and that the experimental $p_T$ differential branching ratio $B(p_T)$ is the same in $pp$ as in PbPb. Then we can calculate the nonprompt $J/\psi$ $p_T$ spectrum in PbPb by multiplying the calculated $B^0(p_T)$ spectrum in PbPb with $B(p_T)$. Doing this we assume that shadowing does not influence the branching ratio, an assumption which creates an additional uncertainty.

## B. Yield ratio

Next we present the yield ratio between different heavy flavor hadrons as a function of $p_T$. Yield ratios of charmed hadrons are important because they present an experimental probe to study fragmentation functions into the different hadron species. This is especially true at high $p_T$ where on the one side hadrons are almost exclusively produced by fragmentation and on the other side hardly any theoretical calculation of these fragmentation functions is available. Yield ratios can also eliminate systematic uncertainties and have therefore a smaller error than absolute yields. They may also contribute to solve physics questions. For example, the yield ratio of $D_s/D^0$ can reveal the strangeness enhancement in the QGP phase of the heavy-ion collisions. The enhancement of the baryon to meson ratio, as compared to the $e^+e^-$ collisions, e.g., $p/\pi$, $\Lambda/K$, shows the importance of hadronization by recombination which is absent in $e^+e^-$ collision where hadronization is exclusively described by fragmentation functions. This enhancement has also been observed in the heavy flavor sector,
<figure><img src="image_24.png" /><figcaption>FIG. 18. Yield ratio of $\Lambda_c/D^0$, $\Xi_c^0/D^0$, and $D_s/D^0$ in the 0–10 % (left) and 30–50 % (right) centrality class of PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. The experimental data are from ALICE, $\Lambda_c/D^0$ [67], $D_s/D^0$ [65,66].</figcaption></figure>

like in the $\Lambda_c/D^0$ ratio, and is even present in high energy pp collisions [67,74].

In Fig. 18 we present the yield ratios of $\Lambda_c/D^0$, $\Xi_c^0/D^0$, and $D_s/D^0$ in the 0–10 % (left) and 30–50 % (right) centrality classes of PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. The experimental data are from the ALICE collaboration, $\Lambda_c/D^0$ from [67], $D_s/D^0$ from [65,66]. We observe that in both centrality classes the experimental $D_s/D^0$ ratio is well reproduced by the EPOS4HQ calculations. This is a strong hint that in the QGP, from where the light quarks originate, the strange quark multiplicity corresponds to the equilibrium value for a system close to $T_{FO}$, as assumed in the EPOS4HQ approach. Also the $\Lambda_c/D^0$ ratio of EPOS4HQ is for both centrality classes and in the whole $p_T$ interval in agreement with experiment [67], including the peak structure of the distribution at around $p_T = 3$ GeV. This ratio is well above the ratio measured in $e^+e^-$ collisions. The enhancement of the baryon to meson ratio in heavy ion collisions can be well understood within the quark recombination/coalescence model [10,15,56,75,76] and therefore the agreement between theory and data gives evidence that at low $p_T$ the majority of $\Lambda_c$ is created in coalescence processes. At higher $p_T$, where fragmentation dominates, the agreement shows that the employed fragmentation functions [50,61] are realistic. We also display the prediction for the $\Xi_c/D^0$ ratio. The shape is similar to the $\Lambda_c/D^0$, but the value is reduced to one half of the former due to the strange quark.

<figure><img src="image_25.png" /><figcaption>FIG. 19. Nuclear modification factor $R_{AA}$ of $D^0$ in the 0–10 %, 10–30 %, and 30–50 % PbPb collisions at 5.02 TeV. The experimental data are from ALICE [64,65]. The black thick (orange lines) lines are the ratio of $D^0$ in PbPb and pp after (before) UrQMD. The blue dotted and green dashed lines are the ratio of charm quarks at production and before hadronization, respectively, which finally are entrained in $D^0$ mesons.</figcaption></figure>

### C. Nuclear modification factor $R_{AA}$

Heavy quarks are produced in primary hard scattering processes, which occur in the early stage of the heavy-ion collisions. Therefore it is useful to define the nuclear modification factor

$$R_{AA} = \frac{1}{\langle N_{\text{coll}} \rangle} \frac{dN_{AA}}{dp_T}, \qquad (5)$$

which compares the $p_T$ distribution in AA collisions with that in pp collisions multiplied by $\langle N_{\text{coll}} \rangle$, the average number of initial hard binary collisions. If the $p_T$ distribution, observed in heavy-ion collisions, is just a superposition of that observed in pp collisions we expect $R_{AA} = 1$. In Fig. 19 the full black
lines show the $R_{AA}$ for the finally observed hadrons, calculated in the EPOS4HQ approach, and the experimental data [64,65] are marked as red points. We see that $R_{AA}$ is smaller than 1. This has two reasons:

(i) Nuclear shadowing modifies the parton distribution function of a nucleus in comparison to that of a proton. In [64] the measured integrated $R_{AA}$ for $D^0$ mesons (0.689 for the centrality [0–10 %] and for $|y| < 0.5$, what EPOS4HQ reproduces within the error bars) is compared with perturbative QCD calculations of $D^0$-meson production including only initial-state effects modeled using two different sets of nuclear parton distribution functions (PDF), namely, nCTEQ15 and EPPS16. These calculations show, however with a large error bar, that the major fraction of the suppression comes from the PDF's, i.e., from the fact that in a nucleus the parton distribution function is different than in a proton. This is confirmed by EPOS4HQ. The $R_{AA}$ for the initially produced charm quarks as around 1 at large $p_T$ but is almost constant ~0.7 at low $p_T$ values, as seen in Fig. 19, top, as a blue dotted line. There we display $R_{AA}$ as a function of $p_T$ for [0–10 %] central PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. After production, heavy quarks propagate in the QGP and interact with the thermal partons by exchanging energy and momentum via elastic and inelastic collisions, as discussed in Sec. III. Due to the energy loss, at high $p_T$ the initial momentum of a c quark will shift to a lower momentum, the spectrum becomes softer and $R_{AA}$ falls well below one (green dashed line). At low $p_T$ the momentum of the c quark increases due to these collisions and the c-quarks approaches thermal equilibrium. This leads to a slight increase of $R_{AA}$.

(ii) Second, due to hadronization by coalescence, which is more important in AA than in pp collisions, we observe at low $p_T$ an increased production of heavy baryons and therefore a lower production of heavy mesons. This lowers $R_{AA}$ at low $p_T$, as seen, when comparing the green dashed and the orange line.

Both processes lead to a complex structure of $R_{AA}$ with a maximum at intermediate $p_T$. For less central collisions, shown in Fig. 19 in the middle for [0–30 %] and in the bottom row for [30–50 %] centrality, the form of $R_{AA}$ as a function of $p_T$ is similar but $R_{AA}$ increases because the energy loss in the plasma gets smaller. The $R_{AA}$ of $D_s$ and $\Lambda_c$ are also investigated and shown in Figs. 20 and 21, which show a good agreement with the experimental data.

Figure 22 shows the centrality dependence of $R_{AA}$ for $D^0$ mesons in AuAu collisions at $\sqrt{s_{NN}} = 200$ GeV. We display $R_{AA}(p_T)$ for two different centrality classes [0–10 %] and [10–40 %] and for the $p_T$ interval for which experimental data are available and not only extrapolations. The STAR data are from [68]. In this figure we display as well the $R_{AA}$ of the those c quarks which are later part of $D^0$ mesons. The blue dotted line shows $R_{AA}$ at the moment when the c quarks are produced. Also at RHIC energies for large $p_T$ c quarks $R_{AA}$ is equal one because shadowing modifies only the c-quark

<figure><img src="image_26.png" /><figcaption>FIG. 20. Nuclear modification factor $R_{AA}$ of $D_s$ in the 0–10 %, 10–30 %, and 30–50 % PbPb collisions at 5.02 TeV. The experimental data are from ALICE [66].</figcaption></figure>

distribution at low $p_T$. There the suppression of c quarks in AuAu can reach 50% as compared to pp. The interaction of heavy quarks with the QGP (difference between green long dashed line and the blue dotted line) modifies strongly the $R_{AA}$ at high $p_T$. Coalescence dominates for low $p_T$ heavy quarks and therefore the momentum of the D meson is larger than that of the heavy c quark embedded, leading to an enhancement of $R_{AA}$ at small $p_T$ (difference between dashed orange and long dashed green line). The final hadronic interactions (difference between black and long dashed green line) have little influence of the final form of $R_{AA}(p_T)$.

Finally we compare in Fig. 23 $R_{AA}(p_T)$ of $B^+$ mesons, measured by the CMS collaboration for PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV with the EPOS4HQ calculations. B mesons are too heavy to approach equilibrium with the QGP. For this reason we observe only a momentum shift of the b quarks
<figure><img src="image_27.png" /><figcaption>FIG. 21. Nuclear modification factor $R_{AA}$ of $\Lambda_c$ in the 0–10%, 10–30%, and 30–50% PbPb collisions at 5.02 TeV. The experimental data are from ALICE [67] and CMS [77].</figcaption></figure>

towards lower $p_T$ values. The data agree quite nicely with the EPOS4HQ calculations.

From the confrontation of the EPOS4HQ results with calculations employing the generalized parton distribution function we can conclude that the shadowing in EPOS4HQ is well described. From the agreement of the EPOS4HQ results with data we can furthermore conclude that also the other processes, which influence the final $p_T$ distributions of open heavy flavour mesons, are quite reasonable described.

### D. Elliptic flow $v_2$

The elliptic flow $v_2$ is another key observable to study the interaction of heavy quarks with the QGP. In noncentral heavy ion collisions, the initial spatial anisotropy of the overlap region is converted into an anisotropic azimuthal distribution in momentum space of the final particles at low $p_T$. This anisotropy can be characterized in terms of Fourier

<figure><img src="image_28.png" /><figcaption>FIG. 22. Nuclear modification factor $R_{AA}$ of $D^0$ in the 0–10% and 10–40% AuAu collisions at 200 GeV. The experimental data are from STAR [68]. The black thick (orange lines) lines are the ratio of $D^0$ in AuAu and $pp$ after (before) UrQMD. The blue dotted and green dashed lines are the ratio of charm quarks at production and before hadronization, respectively, which finally are entrained in $D^0$ mesons.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>FIG. 23. Nuclear modification factor $R_{AA}$ of $B^+$ in the 0–100% PbPb collisions at 5.02 TeV. The experimental data are from the CMS collaboration [70].</figcaption></figure>
coefficients

$$v_n \equiv \langle \cos[n(\Phi - \Psi_n)] \rangle, \qquad (6)$$

where $\Phi$ is the azimuthal angle of the particle and $\Psi_n$ is the azimuthal angle of the event plane for the nth-order harmonics. Elliptic flow, $v_2$, is the second-order coefficient, whose value is in hydrodynamical calculations proportional to the initial spatial eccentricity of the overlap region. Being produced in a hard process, at creation heavy quarks have $v_2 = 0$. So any observed $v_2$ of heavy mesons is due to the interaction of heavy quarks with QGP partons or of heavy mesons with other hadrons. The final elliptic flow of heavy flavor hadrons comes from the heavy quark itself and also from the light antiquark, which merges with the heavy quark in the hadronization process. Also the hadronization process itself may contribute slightly. The measurement of $v_2$ of heavy-flavor hadrons at low $p_T$ can therefore help to quantify the coupling between heavy quarks and the hot medium as well as to understand the hadronization mechanism. At high $p_T$, heavy quarks hadronize via fragmentation. In this case, the elliptic flow $v_2$ together with $R_{AA}$ can be used to explore the path-length dependence of the in-medium energy loss as in the almond shaped overlap region the path length of the heavy quark in the QGP medium depends on its azimuthal direction.

In EPOS4 the event plane $\Psi_n$ can be determined, which is not possible in experiments with a limited acceptance. Therefore a multiparticle correlation or cumulants method has been developed to determine $v_2$. Depending on the number $n$ of particles, whose correlations are calculated, the result is named $v_2\{n\}$. The higher the order, the more nonflow effects, such as resonance decays and jets, are eliminated. This technique has been widely applied in the light-flavor sector [82]. Recently it has also been used for the heavy flavor flow analysis [80].

In Fig. 24 we show the $p_T$ dependence of $v_2\{2\}$ and $v_2\{4\}$ of $D^0$ mesons measured by the CMS collaboration [80] at midrapidity for PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV in different centrality bins. From top to bottom we display the results for the centrality bins [0–10 %], [10–30 %], and [30–50 %]. We see that indeed the experimental $v_2\{4\}$ is slightly smaller than $v_2\{2\}$. In older publications usually $v_2\{2\}$ or the equivalent scalar product (SP) approach has been used. The EPOS4HQ results, obtained with Eq. (6) and shown as well, should be compared with $v_2\{4\}$.

We display in this figure as well the EPOS4HQ results for different times during the evolution of the heavy-ion reaction. Initially the flow of heavy quarks (dotted blue line) is compatible with zero. After passing through the QGP (dashed green line) the c quarks have already acquired $v_2$ by collisions with the QGP partons. Hadronization does not change $v_2$ substantially as can be concluded by comparing the orange short-dashed line ($v_2$ of $D^0$ mesons after hadronization) with the green dashed line but shifts the maximum to larger $p_T$ values, mainly because the $p_T$ of the D meson is different from that of the c quark. Hadronic interactions (the difference between the orange and the black line) still change slightly the value of $v_2(p_T)$. The final distribution is given by the black line which agrees quite nicely with the experimental data. For the [30–50 %] centrality bin we show in Fig. 24 as well $v_2$ of the $B^0$ mesons, as a magenta dashed-dotted line, for which no

<figure><img src="image_30.png" /><figcaption>FIG. 24. Elliptic flow $v_2$ of $D^0$ in the 0–10 %, 10–30 %, and 30–50 % centrality classes of PbPb collisions at $\sqrt{s_{NN}} = 5.02$ TeV. The experimental data are from the ALICE [78] and CMS [79,80] collaborations. The black thick (orange lines) lines are the $v_2$ of $D^0$ after (before) UrQMD. The blue dotted and green dashed lines are $v_2$ of charm quarks at production and before hadronization, respectively. The magenta dashed-dotted line in the middle panel is the $v_2$ of $B^0$ after UrQMD.</figcaption></figure>

data exists. At low $p_T$ it is considerably smaller than that of the $D^0$ mesons due to the large mass of the transfer of $v_2$ to b quarks during the interactions with the QGP partons being less efficient. At higher $p_T$, where the $v_2$ is created by the path length difference, $v_2$ of $B^0$ mesons approaches that of $D^0$ mesons because the mass difference gets increasingly less important.

At RHIC $v_2(p_T)$ of $D^0$ mesons has been measured by the STAR collaboration [81]. The data for semicentral AuAu collisions at $\sqrt{s_{NN}} = 200$ GeV are displayed in Fig. 25 and compared with EPOS4HQ results. We display in this figure also our results for central collisions for which no data are available yet.
<figure><img src="image_31.png" /><figcaption>FIG. 25. Elliptic flow $v_2$ of $D^0$ in the 0–10 % and 10–40 % centrality classes of AuAu collisions at $\sqrt{s_{NN}}$ = 200 GeV. The experimental data are from STAR [81].</figcaption></figure>

## VII. SUMMARY

We used the newly developed EPOS4 approach, which describes quite well a multitude of light hadron observables, to study the physics of heavy quarks. For this we implemented the same description of the dynamics of heavy quarks in the QGP as in the former EPOS3/EPOS2 calculations and included a new hadronization approach, which allows to create in the hadronization process all known heavy hadron species. This approach is named EPOS4HQ. In this approach we investigated the production of heavy hadrons in ultrarelativistic heavy-ion collisions at RHIC and LHC energies. We compared transverse momentum spectra, yield ratio, nuclear modification, and elliptic flow of different charmed hadrons for both PbPb collisions at $\sqrt{s_{NN}}$ = 5.02 TeV and AuAu collisions at $\sqrt{s_{NN}}$ = 200 GeV. As already in pp collisions,

we find a quite good agreement of all observables with the available experimental data. This agreement is considerable better than for the results we obtained using the same description of the heavy quark/heavy meson dynamics in EPOS2 and EPOS3. We can also extend this approach to RHIC energies which was not possible with the EPOS3/EPOS2 approach.

We studied as well in detail how the different observables are modified during the heavy ion collisions. We displayed the initial distribution of c quarks, the modification of the c-quark distribution due to heavy quark-QGP interactions, the influence of hadronization and of the final hadronic rescattering. This allowed for an understanding of these physical processes which influence the heavy quark observables. We observed that the initial parton distribution function in large nuclei shows a strong suppression of low $p_T$ heavy quarks as compared to the scaled pp parton distribution function. This agrees quantitatively with calculations, which use explicitly the generalized parton distribution function. At high $p_T$, where the parton distribution functions in nucleons and heavy nuclei are similar, the interaction of the heavy quarks with the QGP medium are responsible for the low $R_{AA}$ values. The other processes are subdominant but influence quantitatively the results. The elliptic flow at low $p_T$ is as well created by the interaction of the heavy quarks with the QGP medium but its maximum is shifted in the hadronization process. At high $p_T$ it is the path length difference which creates the finite $v_2$ values. Despite EPOS4HQ being based on the fundamental Gribov-Regge approach and describing all presently available heavy hadron data in heavy ion collisions it should be noted that EPOS4HQ needs parametrization to cope with kinematic ranges, which are not accessible by more fundamental approaches. Therefore it cannot be proven whether all the subprocesses are described correctly. The agreement of the EPOS4 results with the light hadron section reduced the uncertainty considerably. An experimental study of the system size dependence would further help to reduce the uncertainties.

## ACKNOWLEDGMENTS

We thank E. Bratkovskaya and T. Song for fruitful discussions. This work is supported by the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 824093 (STRONG-2020).

[1] G. Agakishiev et al. (STAR Collaboration), Phys. Rev. Lett. 108, 072301 (2012).

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

[3] V. Khachatryan et al. (CMS Collaboration), J. High Energy Phys. 04 (2017) 039.

[4] T. Matsui and H. Satz, Phys. Lett. B 178, 416 (1986).

[5] D. Teaney, Phys. Rev. C 68, 034913 (2003).

[6] P. K. Kovtun, D. T. Son, and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005).

[7] S. Acharya et al. (ALICE Collaboration), arXiv:2211.04384 [nucl-ex].

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

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

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

[11] F. Scardina, S. K. Das, V. Minissale, S. Plumari, and V. Greco, Phys. Rev. C 96, 044905 (2017).
[12] V. Minissale, S. Plumari, and V. Greco, Phys. Lett. B **821**, 136622 (2021).

[13] M. Nahrgang, J. Aichelin, S. Bass, P. B. Gossiaux, and K. Werner, Phys. Rev. C **91**, 014904 (2015).

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

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

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

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

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

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

[20] P. B. Gossiaux, S. Vogel, H. van Hees, J. Aichelin, R. Rapp, M. He, and M. Bluhm, arXiv:1102.1114 [hep-ph].

[21] K. Werner, Phys. Rev. C **108**, 064903 (2023).

[22] K. Werner and B. Guiot, Phys. Rev. C **108**, 034904 (2023).

[23] K. Werner, arXiv:2310.09380 [hep-ph].

[24] K. Werner, Phys. Rev. C **109**, 014910 (2024).

[25] P. B. Gossiaux and J. Aichelin, Phys. Rev. C **78**, 014904 (2008).

[26] J. Aichelin, P. B. Gossiaux, and T. Gousset, Phys. Rev. D **89**, 074018 (2014).

[27] V. N. Gribov, Zh. Eksp. Teor. Fiz. **53**, 654 (1967).

[28] V. N. Gribov, Sov. Phys. JETP **29**, 483 (1969).

[29] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. **15**, 438 (1972).

[30] V. A. Abramovsky, V. N. Gribov, and O. V. Kancheli, Yad. Fiz. **18**, 595 (1973).

[31] H. J. Drescher, M. Hladik, S. Ostapchenko, T. Pierog, and K. Werner, Phys. Rep. **350**, 93 (2001).

[32] L. V. Gribov, E. M. Levin, and M. G. Ryskin, Phys. Rep. **100**, 1 (1983).

[33] L. D. McLerran and R. Venugopalan, Phys. Rev. D **49**, 2233 (1994).

[34] L. D. McLerran and R. Venugopalan, Phys. Rev. D **49**, 3352 (1994).

[35] A. Kovner, L. McLerran, and H. Weigert, Phys. Rev. D **52**, 3809 (1995).

[36] Y. V. Kovchegov, Phys. Rev. D **54**, 5463 (1996).

[37] Y. V. Kovchegov and D. H. Rischke, Phys. Rev. C **56**, 1084 (1997).

[38] J. Jalilian-Marian, A. Kovner, L. McLerran, and H. Weigert, Phys. Rev. D **55**, 5414 (1997).

[39] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Nucl. Phys. B **504**, 415 (1997).

[40] Y. V. Kovchegov and A. H. Mueller, Nucl. Phys. B **529**, 451 (1998).

[41] A. Krasnitz and R. Venugopalan, Nucl. Phys. B **557**, 237 (1999).

[42] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert, Phys. Rev. D **59**, 034007 (1999).

[43] J. Jalilian-Marian, A. Kovner, and H. Weigert, Phys. Rev. D **59**, 014015 (1998).

[44] V. Khachatryan et al. (CMS Collaboration), J. High Energy Phys. 09 (2010) 091.

[45] K. Werner, Phys. Rev. Lett. **98**, 152301 (2007).

[46] K. Werner, I. Karpenko, T. Pierog, M. Bleicher, and K. Mikhailov, Phys. Rev. C **82**, 044904 (2010).

[47] K. Werner, B. Guiot, I. Karpenko, and T. Pierog, Phys. Rev. C **89**, 064903 (2014).

[48] I. Karpenko, P. Huovinen, and M. Bleicher, Comput. Phys. Commun. **185**, 3016 (2014).

[49] M. Cacciari, M. Greco, and P. Nason, J. High Energy Phys. 05 (1998) 007.

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

[51] X. Zhu, M. Bleicher, S. L. Huang, K. Schweda, H. Stoecker, N. Xu, and P. Zhuang, Phys. Lett. B **647**, 366 (2007).

[52] M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, Phys. Rev. C **90**, 024907 (2014).

[53] S. Acharya et al. (ALICE Collaboration), J. High Energy Phys. 12 (2023) 086.

[54] H. Berrehrah, P. B. Gossiaux, J. Aichelin, W. Cassing, J. M. Torres-Rincon, and E. Bratkovskaya, Phys. Rev. C **90**, 051901(R) (2014).

[55] L. Altenkort, D. de la Cruz, O. Kaczmarek, R. Larsen, G. D. Moore, S. Mukherjee, P. Petreczky, H.-T. Shu, and S. Stendebach (HotQCD Collaboration), Phys. Rev. Lett. **132**, 051902 (2024).

[56] J. Zhao, S. Shi, N. Xu, and P. Zhuang, arXiv:1805.10858 [hep-ph].

[57] A. Andronic, P. Braun-Munzinger, K. Redlich, and J. Stachel, Phys. Lett. B **659**, 149 (2008).

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

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

[60] P. Madanagopalan, R. G. Edwards, N. Mathur, and M. J. Peardon, PoS (LATTICE2014), 084 (2015).

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

[62] M. Lisovyi, A. Verbytskyi, and O. Zenaiev, Eur. Phys. J. C **76**, 397 (2016).

[63] S. A. Bass et al., Prog. Part. Nucl. Phys. **41**, 255 (1998).

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

[65] S. Acharya et al. (ALICE Collaboration), J. High Energy Phys. 10 (2018) 174.

[66] S. Acharya et al. (ALICE Collaboration), Phys. Lett. B **827**, 136986 (2022).

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

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

[69] J. Adam et al. (STAR Collaboration), Phys. Rev. Lett. **127**, 092301 (2021).

[70] A. M. Sirunyan et al. (CMS Collaboration), Phys. Rev. Lett. **119**, 152301 (2017).

[71] A. M. Sirunyan et al. (CMS Collaboration), Phys. Lett. B **796**, 168 (2019).

[72] M. Aaboud et al. (ATLAS Collaboration), Eur. Phys. J. C **78**, 762 (2018).

[73] S. Acharya et al. (ALICE Collaboration), J. High Energy Phys. 02 (2024) 066.

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

[75] Y. Oh, C. M. Ko, S. H. Lee, and S. Yasui, Phys. Rev. C **79**, 044905 (2009).

[76] 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).
[77] A. Tumasyan et al. (CMS Collaboration), J. High Energy Phys. 01 (2024) 128.

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

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

[80] A. Tumasyan et al. (CMS Collaboration), Phys. Rev. Lett. 129, 022001 (2022).

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

[82] S. Chatrchyan et al. (CMS Collaboration), Phys. Rev. C 87, 014902 (2013).