2021 213 13 October

Prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ $R_{\rm AA}$ in Pb–Pb

ALICE Collaboration††thanks: See Appendix A for the list of collaboration
members ALICE Collaboration

The production of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$
mesons was measured at midrapidity ($|y|<\leavevmode\nobreak\ $0.5) in Pb–Pb
collisions at the centre-of-mass energy per nucleon–nucleon pair
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$
with the ALICE detector at the LHC. The D mesons were reconstructed via their
hadronic decay channels and their production yields were measured in central
(0–10%) and semicentral (30–50%) collisions. The measurement was performed up
to a transverse momentum ($p_{\rm T}$) of 36 or 50 $\mathrm{GeV}/c$ depending
on the D meson species and the centrality interval. For the first time in
Pb–Pb collisions at the LHC, the yield of ${\rm D^{0}}$ mesons was measured
down to $p_{\rm T}=$ 0, which allowed a model-independent determination of the
$p_{\rm T}$-integrated yield per unit of rapidity (${\rm d}N/{\rm d}y$). A
maximum suppression by a factor 5 and 2.5 was observed with the nuclear
modification factor ($R_{\rm AA}$) of prompt D mesons at $p_{\rm T}=$
6–8$\leavevmode\nobreak\ \mathrm{GeV}/c$ for the 0–10% and 30–50% centrality
classes, respectively. The D-meson $R_{\rm AA}$ is compared with that of
charged pions, charged hadrons, and ${\rm J}/\psi$ mesons as well as with
theoretical predictions. The analysis of the agreement between the measured
$R_{\rm AA}$, elliptic ($v_{2}$) and triangular ($v_{3}$) flow, and the model
predictions allowed us to constrain the charm spatial diffusion coefficient
$D_{s}$. Furthermore the comparison of $R_{\rm AA}$ and $v_{2}$ with different
implementations of the same models provides an important insight into the role
of radiative energy loss as well as charm quark recombination in the
hadronisation mechanisms.

##  1 Introduction

Ultra-relativistic heavy-ion collisions allow the study of quantum
chromodynamics (QCD) at high energy density and temperature. According to
lattice QCD (lQCD) calculations, these extreme conditions lead to a transition
of hadronic matter into a strongly-interacting medium, called quark–gluon
plasma (QGP), in which quarks and gluons are deconfined and chiral symmetry is
partially restored [1, 2, 3, 4, 5]. A QGP is formed and studied in high-energy
heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and at the
CERN Large Hadron Collider (LHC). The existing measurements indicate that the
QGP behaves as a strongly-coupled low-viscosity liquid-like system [6] and its
lifetime at the energy densities reached at the LHC is of the order of 10
$\mathrm{fm}/c$ [7].

Heavy quarks, like charm and beauty, are mostly produced in primary hard
scattering processes between the partons of the incoming nuclei, which occur
in the early stages of the collisions. The time scales for heavy-quark
production ($\leq$ 0.07 $\mathrm{fm}/c$ for ${\rm c}\overline{\rm c}$ and
$\leq$ 0.02 $\mathrm{fm}/c$ for ${\rm b}\overline{\rm b}$ pairs [8]) are
shorter than the QGP formation time, which is about 0.3–1.5 $\mathrm{fm}/c$ at
LHC energies [9]. Therefore, heavy quarks experience the full space–time
evolution of the hot and dense QCD medium.

During their propagation in the QGP, heavy quarks interact with the medium
constituents, by exchanging energy and momentum via elastic [10, 11, 12] and
inelastic [13, 14] processes. Low-momentum heavy quarks (i.e. $\lesssim$ 3–4
$\mathrm{GeV}/c$) mainly interact via elastic scatterings and are expected to
acquire some collective behaviour of the system, like radial and anisotropic
azimuthal flow, as a consequence of multiple interactions with the medium [15,
16]. The typical momentum exchange in the interactions of heavy quarks with
the medium is small compared to the charm-quark mass. Therefore, heavy quarks
undergo Brownian motion in the medium characterised by many small-momentum
kicks. The degree of heavy-quark thermalisation provides unique insight into
the thermalisation process of the system. Substantial theoretical efforts were
made to describe the charm-quark transport in a hydrodynamically expanding
medium [17]. The relevant transport parameter is the spatial diffusion
coefficient $D_{s}$, which does not depend on the precise value of the heavy-
quark mass and, hence, is a medium property. It can be estimated via the
comparison of the measured differential yields as a function of transverse
momentum ($p_{\rm T}$) and azimuthal distributions of D mesons with model
predictions.

For high-$p_{\rm T}$ quarks (i.e. $\gtrsim$ 8–10 $\mathrm{GeV}/c$) the main
manifestation of the interactions with the medium is energy loss, which can
occur via collisional processes and gluon radiation [18], with the latter
expected to be the dominant mechanism. In particular, the amount of radiative
energy loss is influenced by the colour charge of the parton interacting with
the medium. According to QCD calculations, gluons are expected to lose more
energy than quarks, due to their stronger coupling to the medium [14, 13]. At
LHC energies, light-flavour particles in the momentum interval 5 $<p_{\rm T}<$
20 $\mathrm{GeV}/c$ are expected to originate mostly from the fragmentation of
gluons produced in hard scattering processes, while at higher $p_{\rm T}$ the
contribution of light-quark jets becomes relevant [19]. On the other hand,
charm mesons provide an experimental tag for a quark parent at all momenta.
Thus, the comparison of heavy-flavour hadron production with that of light-
flavour particles at high $p_{\rm T}$ can provide insight into this aspect. In
addition, the energy loss of quarks can be affected by several mass-dependent
effects like the dead-cone effect, which reduces the small-angle gluon
radiation for quarks with moderate energy-over-mass values [20, 21, 22, 23].
Similarly, collisional energy loss is predicted to depend on the quark mass
and to be smaller for heavy quarks [24].

Heavy-flavour hadron yields and momentum distributions can be influenced by
the modification of the parton distribution function (PDF) due to initial
state effects, like nuclear shadowing [18]. It was also suggested that low-
momentum heavy quarks could hadronise via recombination with other quarks from
the medium, in addition to the fragmentation in the vacuum [16, 25].

The aforementioned effects can be investigated with the measurement of the
nuclear modification factor $R_{\rm AA}$ of heavy-flavour hadrons. It is
defined as the ratio of the $p_{\rm T}$-differential production yield in
nucleus–nucleus collisions (${\rm d}N_{\rm AA}/{\rm d}p_{\rm T}$) and the
production cross section in proton–proton collisions (${\rm d}\sigma_{\rm
pp}/{\rm d}p_{\rm T}$) scaled by the average nuclear overlap function $\langle
T_{\rm AA}\rangle$

$R_{\rm AA}(p_{\rm T})={1\over\left\langle T_{\rm AA}\right\rangle}{{\rm
d}N_{\rm AA}/{\rm d}p_{\rm T}\over{\rm d}\sigma_{\rm pp}/{\rm d}p_{\rm T}},$
(1)

where $\langle T_{\rm AA}\rangle$ is defined as the average number of binary
nucleon–nucleon collisions ($\langle N_{\rm coll}\rangle$), which can be
estimated via Glauber model calculations [26, 27, 28], divided by the
inelastic nucleon–nucleon cross section.

Measurements of prompt D-meson production were performed by the ALICE
Collaboration in Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=$
2.76 TeV [29] and 5.02 $\mathrm{TeV}$ [30] in the 0–10% and 30–50% centrality
classes for both collision energies and in the 60–80% for the highest one. The
D-meson yields show a suppression in central (0–10$\%$) and semi-central
(30–50$\%$) collisions that reaches up to a factor of 5 and 2.5, respectively,
at $p_{\rm T}$ of 6–10 GeV/c. It then decreases with increasing $p_{\rm T}$,
while the suppression factor is of 1.25 without a pronounced $p_{\rm T}$
dependence in peripheral (60–80$\%$) collisions. The prompt ${\rm D^{0}}$
production was also measured by the CMS Collaboration in the most central
Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$ [31] in the $p_{\rm T}$ range
2–100 $\mathrm{GeV}/c$ and the result is in agreement with ALICE measurements.
Furthermore, also the STAR Collaboration measured the production of the prompt
${\rm D^{0}}$ mesons in Au–Au collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=$ 200 GeV [32] and the comparison with this result can provide insight
into the $\sqrt{s}$ dependence of energy-loss effects. Conversely, the D-meson
nuclear modification factor measured by the ALICE Collaboration in p–Pb
collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\
\mathrm{TeV}$ [33], where an extended QGP is not expected to form, is
compatible with unity for the whole measured $p_{\rm T}$ range.

Complementary information on the interactions of heavy quarks with the medium
constituents and their possible thermalisation can be obtained from the
measurements of the anisotropies in the D-meson azimuthal distributions, which
are characterised in terms of the Fourier coefficients [34, 35]. In
particular, the second order coefficient $v_{2}$, called elliptic flow, is
defined as $v_{2}=\langle\cos(\mathrm{2}\varphi)\rangle$ where $\varphi$ is
the azimuthal angle. Recently, the ALICE Collaboration reported the latest
measurements of prompt D-meson azimuthal anisotropic flow [36] in Pb–Pb
collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\
\mathrm{TeV}$.

In this paper, the measurements of the $p_{\rm T}$-differential yields and
nuclear modification factors of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm
D^{*+}}$ mesons (i.e. produced via the hadronisation of charm quarks or from
the decays of excited charmed hadron states) and their charge conjugates, in
central (0–10$\%$) and semicentral (30–50$\%$) Pb–Pb collisions at
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$
are reported. The statistical and systematic uncertainties have been reduced
due to the Pb–Pb data sample collected with the ALICE detector at the end of
2018, that is larger by a factor eight (four) for central (semicentral)
collisions with respect to the data sample collected in 2015 used for the
previous publication [30], and to the production cross section measured in pp
collisions at the same centre-of-mass energy [37]. The larger data sample
allowed the reduction of the statistical uncertainties and the extension of
the $p_{\rm T}$ reach of the measurements down to $p_{\rm T}=$ 0 for ${\rm
D^{0}}$ mesons. Moreover, the previous measurements [30] were affected by the
uncertainty on the $\sqrt{s}$ scaling of the pp production cross section. The
new measurements do not include such uncertainty since they use the measured
pp production cross section at the same collision energy [37].

The paper is structured as follows. The data sample and the experimental
apparatus are briefly introduced in Section 2. The D-meson reconstruction
procedure and the corrections applied to the raw yields are presented in
Section 3, while the estimation of the systematic uncertainties is described
in Section 4. The results are presented in Section 5, together with a
comparison with the charged-pion, proton, and prompt and non-prompt ${\rm
J}/\psi\leavevmode\nobreak\ R_{\rm AA}$, along with model calculations. In
addition, a simultaneous comparison of the measured $R_{\rm AA}$ and elliptic
flow $v_{2}$ [36] at low and intermediate $p_{\rm T}$ with transport-model
calculations is discussed, together with an estimate of the spatial diffusion
coefficient. This comparison will provide insights into the participation of
the charm quark in the collective expansion of the medium [38], as well as on
the role of the different hadronisation mechanisms. Conclusions and
perspectives are drawn in Section 6.

##  2 Detector and data sample

The ALICE experimental apparatus includes several detectors for particle
tracking and identification as well as electromagnetic calorimeters at
midrapidity ($|\eta|<0.9$), a forward muon spectrometer ($-4<\eta<-2.5$), and
a set of forward and backward detectors used for triggering, background
rejection, and event characterisation. A detailed description of the apparatus
and its performance can be found in Refs. [39, 40].

The main detectors used for the track reconstruction and particle
identification in the analysis are the Inner Tracking System (ITS) [41], the
Time Projection Chamber (TPC) [42], and the Time-Of-Flight (TOF) detector
[43]. They are located inside a large solenoidal magnet, providing a uniform
magnetic field of 0.5 T parallel to the LHC beam direction. Charged-particle
trajectories are reconstructed from their hits in the ITS and the TPC. The ITS
consists of a six-layer silicon detector, which is used for tracking and for
the reconstruction of primary and secondary vertices. The TPC is used for
track reconstruction and particle identification (PID) via the measurement of
the specific energy loss d$E/$d$x$. The TOF detector extends at intermediate
momenta the PID capabilities of the TPC via the measurement of the time-of-
flight of charged particles from the interaction point to the detector.

Event selection and characterisation are performed with the V0 detector [44]
and the Zero Degree Calorimeter (ZDC) [45]. The V0 detector consists of two
scintillator arrays, which cover the full azimuth in the pseudorapidity
intervals $-3.7<\eta<-1.7$ (V0C) and $2.8<\eta<5.1$ (V0A), and is used for
triggering, event selection, and centrality determination. A minimum-bias
interaction trigger was provided by the coincidence of signals in the V0A and
the V0C. Furthermore, two separate trigger classes, which consist of an online
event selection based on the V0 signal amplitude, were used during the 2018
Pb–Pb data taking period in order to enrich the sample of central and
semicentral collisions. Events due to the interaction of the beams with
residual gas in the vacuum pipe were rejected offline using the V0 and the ZDC
timing information. Only events with a primary vertex reconstructed within
$\pm 10$ cm from the centre of the detector along the beam axis were
considered in the analysis.

Collisions were divided into centrality classes, defined in terms of
percentiles of the hadronic Pb–Pb cross section, using the sum of the V0
signal amplitudes, as described in details in Ref. [30]. The centrality
classes used in the present analysis, together with the corresponding average
nuclear overlap function $\left\langle T_{\rm AA}\right\rangle$ [46, 47], the
number of events ($N_{\rm events}$), and the recorded integrated luminosity
$L_{\rm int}$ [48] in each class, are summarised in Table 1.

Table 1: Average nuclear overlap function, number of events, and recorded integrated luminosities for the two centrality classes used in the analysis. Centrality class |  $\langle T_{\rm AA}\rangle$ (mb-1) | $N_{\rm events}$ | $L_{\rm int}(\mu\mathrm{b}^{-1})$
---|---|---|---
00–10% | $23.26\pm 0.17$ | $100\times 10^{6}$ | 130.5 $\pm$ 0.5
30–50% | $3.92\pm 0.06$ | $85\times 10^{6}$ | 55.5 $\pm$ 0.2

##  3 Analysis techniques

###  3.1 D-meson raw yields

The D mesons and their charge conjugates were reconstructed via the following
hadronic decay channels: ${\rm D}^{0}\to{\rm K}^{-}\pi^{+}$ (with branching
ratio, $\mathrm{BR}=3.950\pm 0.031\%$), ${\rm D}^{+}\to{\rm
K}^{-}\pi^{+}\pi^{+}$ ($\mathrm{BR}=9.38\pm 0.16\%$), and ${\rm D}^{*+}\to{\rm
D}^{0}\pi^{+}$ ($\mathrm{BR}=67.7\pm 0.5\%$) [49]. The analysis was based on
the reconstruction of decay vertices displaced from the primary vertex,
exploiting the separation of few hundreds $\mu$m due to the weak decays of
${\rm D^{0}}$ and ${\rm D^{+}}$ (c$\tau$ of $\sim$ 123 and $\sim$ 312 $\mu$m,
respectively [49]). In the case of the strong decay of the ${\rm D^{*+}}$
meson, the decay topology of the produced ${\rm D^{0}}$ was exploited. This
method allowed the reconstruction of ${\rm D^{0}}$ candidates for $p_{\rm
T}>1\leavevmode\nobreak\ \mathrm{GeV}/c$ in both centrality classes and of
${\rm D^{+}}$ (${\rm D^{*+}}$) candidates for $p_{\rm T}>$ 2.5
(3)$\leavevmode\nobreak\ \mathrm{GeV}/c$ and for $p_{\rm T}>$ 2
(2)$\leavevmode\nobreak\ \mathrm{GeV}/c$ in the 0–10% and 30–50% centrality
classes, respectively.

${\rm D^{0}}$ and ${\rm D^{+}}$ candidates were built using pairs and triplets
of tracks with proper charge-sign combination, $|\eta|<0.8$, $p_{\rm T}>$ 0.4
GeV/$c$, a minimum number of two hits (out of six) in the ITS, with at least
one in the two innermost layers, at least 70 out of 159 crossed TPC pad rows,
and a fit quality $\chi^{2}/$ndf $<$ 1.25 in the TPC (where ndf is the number
of degrees of freedom involved in the track fit procedure). ${\rm D^{*+}}$
candidates were reconstructed by combining ${\rm D^{0}}$ candidates with
tracks, which were required to have $|\eta|<$ 0.8, $p_{\rm T}>$ 0.1
$\mathrm{GeV}/c$, and at least three associated hits in the ITS.

These track selection criteria reduce the D-meson acceptance in rapidity,
which, depending on $p_{\rm T}$, varies from $|y|>0.5$ at low $p_{\rm T}$ and
$|y|>$ 0.8 for $p_{\rm T}>$ 5 GeV/$c$. Consequently, a $p_{\rm T}$-dependent
fiducial acceptance selection, $|y|<y_{\mathrm{fid}}(p_{\rm T})$, was applied,
with $y_{\mathrm{fid}}(p_{\rm T})$ defined according to a second-order
polynomial function increasing from 0.5 to 0.8 in the range $0<p_{\rm
T}<5\leavevmode\nobreak\ \mathrm{GeV}/c$, and fixed to 0.8 for $p_{\rm
T}>5\leavevmode\nobreak\ \mathrm{GeV}/c$.

The D-meson selection strategy is similar to the one used in previous analyses
and the variables used to distinguish between signal and background candidates
are based on the impact parameters (i.e. the distance of closest approach
between the track and the primary vertex in the plane transverse to the beam
direction) of the decay particles, the separation of the primary and secondary
vertices, and the pointing of the reconstructed D-meson momentum to the
primary vertex, as described in Refs. [50, 30]. These selection criteria tend
to enhance reconstructed feed-down D mesons, originating from beauty-hadron
decays, over those promptly produced. An additional selection was therefore
applied on the normalised difference between the measured and expected impact
parameters of each of the decay particles, which allows for a significant
rejection of background candidates and feed-down D mesons, while keeping a
large fraction of prompt D mesons as reported in Ref. [51]. Further reduction
of the combinatorial background was obtained by applying PID for charged pions
and kaons with the TPC and TOF detectors. A $\pm$3$\sigma$ window around the
expected mean values of d$E/$d$x$ and time-of-flight was used for the
identification, where $\sigma$ is the resolution on these two quantities. A
tighter PID selection criterion ($\pm 2\sigma$) was used for ${\rm D^{+}}$
(for $p_{\rm T}<$ 3 GeV/$c$) and ${\rm D^{*+}}$ candidate daughters in central
collisions because of the large background of track triplets.

An alternative analysis method, not based on the reconstruction and selection
of the displaced decay-vertex topology, was previously developed and applied
for the two-body decay of ${\rm D^{0}}$ mesons in pp and p–Pb collisions, in
order to extend the measurement down to $p_{\rm T}=0$ [52, 37]. This
technique, which is used for the first time in the analysis of Pb–Pb
collisions, does not use geometrical selections on the displaced decay vertex
(which are not effective at very low $p_{\rm T}$), but relies mainly on
particle identification and on the estimation and subtraction of the
combinatorial background. The ${\rm D^{0}}$ candidates were formed combining
pairs of tracks with opposite charge sign (unlike sign, ULS), with $|\eta|<$
0.8 and $p_{\rm T}>$ 0.4 $\mathrm{GeV}/c$. Single-track and particle
identification selections as well as the fiducial acceptance selection of the
reconstructed candidates were performed using the same strategy adopted in the
analysis with decay-vertex reconstruction previously described. The
contribution to the invariant-mass distribution of ULS K$\pi$ pairs due to
combinatorial background is estimated with the event-mixing technique. The
events are grouped in pools, before the event-mixing, based on the primary
vertex position along z and on the event multiplicity, quantified by the
number of track segments reconstructed in the two innermost layers of the ITS.
Then, the kaon tracks of a given event are mixed with the pion tracks of other
events. The top left panel of Fig. 1 shows the invariant-mass distribution of
ULS K$\pi$ pairs together with that of the background estimated with the
event-mixing technique in the interval $0<p_{\rm T}<1\leavevmode\nobreak\
\mathrm{GeV}/c$. The latter distribution was normalised to match the yield of
ULS pairs at one edge of the invariant-mass interval considered for the
extraction of the ${\rm D^{0}}$ raw yield.

The D-meson raw yields, including particles and antiparticles, were obtained
from binned maximum-likelihood fits to the ${\rm D^{0}}$ and ${\rm D^{+}}$
candidate invariant mass ($M$) distribution and to the mass difference $\Delta
M=M(\mathrm{K}\pi\pi)-M(\mathrm{K}\pi)$ distribution for ${\rm D^{*+}}$
candidates. Figure 1 shows examples of fits for these distributions in the
0–10$\%$ centrality class and for four $p_{\rm T}$ intervals. In the analysis
of the ${\rm D^{0}}$ meson without decay-vertex reconstruction, the fit is
performed after subtracting the estimated background from the ULS K$\pi$
invariant-mass distribution as shown in the top right panel.

Figure 1: Invariant-mass (mass-difference) distribution of ${\rm D^{0}}$,
${\rm D^{+}}$, and ${\rm D^{*+}}$ meson candidates in different $p_{\rm T}$
intervals for the centrality class 0–10$\%$. The fit functions are composed of
a Gaussian function for the signal and an additional term for the background,
as described in detail in the text. The values for the Gaussian mean $\mu$,
width $\sigma$, and raw yield S are also reported. Top row: ${\rm
D^{0}}$-meson candidates with $0<p_{\rm T}<1\leavevmode\nobreak\
\mathrm{GeV}/c$ without reconstructing the decay vertex. The invariant-mass
distributions are shown before (left) and after (right) the subtraction of the
combinatorial background estimated from event-mixing. Middle row: ${\rm
D^{0}}$-meson candidates with $1<p_{\rm T}<1.5\leavevmode\nobreak\
\mathrm{GeV}/c$ with reconstruction of the decay vertex before (left) and
after (right) background subtraction. The width of the Gaussian in this and in
the previous $p_{\rm T}$ interval is fixed to the value obtained from
simulations. Bottom row: ${\rm D^{+}}$-meson candidates with $5<p_{\rm
T}<5.5\leavevmode\nobreak\ \mathrm{GeV}/c$ (left) and ${\rm D^{*+}}$-meson
candidates with $16<p_{\rm T}<24\leavevmode\nobreak\ \mathrm{GeV}/c$ (right).

The fit function was composed of a Gaussian term to describe the signal and an
exponential function for the background for ${\rm D^{0}}$ and ${\rm D^{+}}$
candidates at intermediate and high $p_{\rm T}$ in the analysis with vertex
reconstruction. At low $p_{\rm T}$ ($p_{\rm T}<3\leavevmode\nobreak\
\mathrm{GeV}/c$ for ${\rm D^{0}}$ and $p_{\rm T}<5\leavevmode\nobreak\
\mathrm{GeV}/c$ for ${\rm D^{+}}$), a second-order polynomial function is used
to model the background invariant-mass shape in both centrality classes. In
the ${\rm D^{0}}$ analysis without vertex reconstruction, instead, the
background was parametrised for $p_{\rm T}<$ ($>$) 2$\leavevmode\nobreak\
\mathrm{GeV}/c$ with a fourth-order (third-order) polynomial function. The
$\Delta M$ distribution of ${\rm D^{*+}}$ candidates was fitted with a
Gaussian function for the signal and the following function for the
background: $a\sqrt{\Delta M-m_{\pi}}{\rm e}^{b(\Delta M-m_{\pi})}$, where
$m_{\pi}$ is the pion mass and $a$ and $b$ are free parameters.

In some cases a ${\rm D^{0}}$-meson candidate can be also identified as a
${\overline{{\rm D}}\,^{0}}$ meson when the two daughter tracks are compatible
with both the kaon and pion hypothesis, leading to an irreducible correlated
background. Hence, in the ${\rm D^{0}}$ meson analyses, an additional term was
included in the fit function to take into account this contribution to the
invariant mass of signal candidates, which is called reflections. It was
parametrised with a double-Gaussian distribution based on detailed Monte Carlo
simulations. This contribution ranges between 2% and 3% of the raw signal
depending on $p_{\rm T}$ [30]. Given the critical signal extraction for the
${\rm D^{0}}$ meson at low $p_{\rm T}$ in the analysis with vertex
reconstruction, due to the small signal-to-background ratio, the width of the
Gaussian for the signal was fixed to the value obtained from the simulation
for $p_{\rm T}<1.5\leavevmode\nobreak\ \mathrm{GeV}/c$. We verified that the
widths from the simulation were consistent with those extracted from the data.
For the same reason, in the analysis without vertexing the Gaussian width was
fixed to the one from the simulation in the whole $p_{\rm T}$ range. The
statistical significance of the signal, defined as $S/\sqrt{S+B}$ where $S$ is
the raw signal yield obtained by integrating the Gaussian function and $B$ is
the background under the peak (within 3$\sigma$), varies from 4 to 60,
depending on the D-meson species, the $p_{\rm T}$ interval, and the centrality
class. In the ${\rm D^{0}}$ analysis without decay-vertex reconstruction, the
$S/B$ in the $p_{\rm T}$ range 0–1 $\mathrm{GeV}/c$ is approximately $\sim$
8.5$\times$10-6 before the mixed-event background subtraction.

###  3.2 Yield corrections and beauty feed-down subtraction

The D-meson raw yields were corrected to obtain the $p_{\rm T}$-differential
yields of prompt D-mesons according to:

$\frac{\mathrm{d^{2}}N}{\mathrm{d}y\leavevmode\nobreak\
\mathrm{d}p_{\mathrm{T}}}=\frac{1}{2}\frac{1}{\Delta
p_{\mathrm{T}}}\frac{f_{\mathrm{prompt}}(p_{\mathrm{T}})\times
N^{\mathrm{D+\overline{D}}}_{\mathrm{raw}}(p_{\mathrm{T}})\Big{\rvert}_{|y|<y_{\mathrm{fid}}}}{\alpha_{y}(p_{\mathrm{T}})\times(\mathrm{Acc}\times\epsilon)_{\mathrm{prompt}}(p_{\mathrm{T}})\times\mathrm{BR}\times
N_{\mathrm{events}}}.$ (2)

The raw yield values $N^{\rm D+\overline{D},raw}$, which contain the
contribution of feed-down from beauty-hadron decays, were multiplied by the
fraction of promptly produced D mesons $f_{\rm{prompt}}$ and divided by a
factor of two to obtain the charged-averaged yields. Furthermore, they were
divided by the product of prompt D-meson acceptance and efficiency
($\mathrm{Acc}\times\epsilon)_{\mathrm{prompt}}$, the branching ratio BR of
the decay channel, the width of the $p_{\rm T}$ interval ($\Delta p_{\rm T}$),
and by the number of events $N_{\mathrm{events}}$. The factor
$\alpha_{y}(p_{\mathrm{T}})$ normalises the corrected yields measured in
$\Delta y=$ 2$y_{\mathrm{fid}}(p_{\mathrm{T}})$ to one unit of rapidity. It
was computed as the ratio between the generated D-meson yield in
$|y|<y_{\mathrm{fid}}(p_{\mathrm{T}})$ and that in $|y|<0.5$ using a data-
driven $p_{\rm T}$ shape and a rapidity distribution from fixed order plus
next-to-leading logarithms (FONLL) [53, 54] perturbative QCD calculations,
after verifying that the D-meson rapidity distributions in FONLL are
consistent with those from PYTHIA 8 [55].

The ($\mathrm{Acc}\times\epsilon$) correction was obtained separately for
prompt and feed-down D mesons using simulations with the GEANT3 transport
package [56], including a detailed description of the detector geometry and
response as well as of the LHC beam conditions. The Pb–Pb collisions at
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$
were produced with the HIJING v1.36 [57] event generator and D-meson signals
were added by injecting ${\rm c\overline{c}}$ or ${\rm b\overline{b}}$ pairs
generated with the PYTHIA v8.2.43 event generator with Monash tune [58]. The
D-meson $p_{\rm T}$ distributions from the simulations were reweighted in
order to use realistic momentum distributions in the determination of the
acceptance and the efficiency, which depend on $p_{\rm T}$. The weights were
defined with an iterative procedure to match the $p_{\rm T}$ dependence
measured for ${\rm D^{0}}$ mesons in the intervals used in the analysis, for
both centrality classes and in the momentum interval 1$<p_{\rm T}<$
50$\leavevmode\nobreak\ \mathrm{GeV}/c$.

Figure 2: Product of acceptance and efficiency ($\mathrm{Acc}\times\epsilon$)
as a function of $p_{\rm T}$ for prompt (red circles) and feed-down (blue
squares) D mesons in Pb–Pb collisions for the 0–10% centrality class.

Figure 2 shows the ($\mathrm{Acc}\times\epsilon$) for prompt and feed-down
${\rm D^{0}}$ (top panels), ${\rm D^{+}}$ (bottom left panel), and ${\rm
D^{*+}}$ (bottom right panel) mesons with $|y|<y_{\mathrm{fid}}$ as a function
of $p_{\rm T}$ in the 0–10% centrality class. For the analysis that does not
exploit the selections on the ${\rm D^{0}}$-meson decay vertex (top left
panel), the efficiency is higher by a factor of about 100 (5) at low (high)
$p_{\rm T}$ as compared to the analysis with decay-vertex reconstruction (top
right panel). The $({\rm Acc}\times\epsilon)$ is almost independent of $p_{\rm
T}$ (the mild increase with increasing $p_{\rm T}$ is mainly determined by the
geometrical acceptance of the detector) and is the same for prompt and feed-
down ${\rm D^{0}}$, as expected when no selection is made on the displacement
of the ${\rm D^{0}}$-meson decay vertex from the interaction point.

The efficiencies for the analysis with geometrical selections on the displaced
decay-vertex topology range from about 10-3 at low $p_{\rm T}$ to 0.1–0.3 at
high $p_{\rm T}$, depending on the D-meson species. The decreasing efficiency
with decreasing $p_{\rm T}$ is due to the fact that more stringent selection
criteria are needed at low $p_{\rm T}$ because of the larger background and
the smaller average displacement of the decay vertex from the interaction
point. The trend of the efficiency with $p_{\rm T}$ is not completely smooth,
especially for ${\rm D^{0}}$ and ${\rm D^{*+}}$, reflecting the $p_{\rm T}$
intervals in which the optimisation of the selection criteria was performed.
The ($\mathrm{Acc}\times\epsilon$) is higher for semicentral collisions than
for central collisions, since less stringent selections are applied because of
the lower combinatorial background. The difference between the efficiencies
for prompt and feed-down D mesons is due to the geometrical selections applied
since the latter are more displaced from the primary vertex due to large
B-meson lifetime ($\tau\approx$ 500 $\mu$m) and are therefore more efficiently
selected by the majority of the selection criteria applied in the analysis.
However, the requirement on the difference between measured and expected decay
track impact parameter rejects efficiently D mesons coming from beauty-hadron
decays. This is particularly effective for ${\rm D^{+}}$ mesons (Fig. 2 bottom
left panel) whose feed-down efficiencies become lower than the prompt
efficiencies after a selection on the aforementioned variable is applied.

The fraction of selected prompt D mesons $f_{\rm prompt}$ was obtained with
the procedure introduced in Refs. [59, 30] to account for the contribution of
D mesons from beauty-hadron decays in the measured raw yield. The $f_{\rm
prompt}$ factor was estimated in each $p_{\rm T}$ interval using perturbative
QCD calculations for the cross section of feed-down D mesons, efficiencies
from the simulations, and a hypothesis on the $R_{\rm AA}$ of feed-down D
mesons. In detail, the expression for $f_{\rm prompt}$ is:

$\begin{split}f_{\rm prompt}&=1-\frac{N^{\rm D+\overline{D}\,\textnormal{feed-
down}}_{\rm raw}}{N^{\rm D+\overline{D}}_{\rm raw}}\\\ &=1-R_{\rm
AA}^{\textnormal{feed-down}}\times\langle T_{\rm
AA}\rangle\times\left(\frac{{\rm d}\sigma}{{\rm d}p_{\rm
T}}\right)^{\textnormal{FONLL,\,PYTHIA 8}}_{{\textnormal{ feed-
down},\,|y|<0.5}}\times\frac{2\Delta p_{\rm T}\times\alpha_{y}\times({\rm
Acc}\times\epsilon)_{\textnormal{feed-down}}\times{\rm BR}\times N_{\rm
events}}{N^{\rm D+\overline{D}}_{\rm raw}}\,,\end{split}$ (3)

where $N^{\rm D+\overline{D}}_{\rm raw}$ is the measured raw yield and
$N^{{\rm D+\overline{D}}\,\textnormal{feed-down}}_{\rm raw}$ is the estimated
raw yield of $\rm D$ mesons from beauty-hadron decays. The beauty-hadron
production cross section in pp collisions at
$\sqrt{s}=5.02\leavevmode\nobreak\ \mathrm{TeV}$, estimated with FONLL
calculations [60], was folded with the beauty-hadron${\rightarrow\rm D}+X$
decay kinematics from PYTHIA 8 and multiplied by the $({\rm
Acc}\times\varepsilon)$ of feed-down $\rm D$ mesons, by the $\langle T_{\rm
AA}\rangle$ of the corresponding centrality class, and by the other factors
introduced in Eq. 2. Finally, the nuclear modification factor of D mesons from
beauty-hadron decays ($R_{\rm AA}^{\textnormal{feed-down}}$) was accounted
for.

The comparison of the $R_{\rm AA}$ of prompt D mesons ($R_{\rm AA}^{\rm
prompt}$) at $\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\
\mathrm{TeV}$ [30] with that of $\rm J/\psi$ from B-meson decays at the same
energy measured by the CMS [61] and the ATLAS [62] collaborations indicates
that prompt charmed hadrons are more suppressed than non-prompt charmed
hadrons. The difference between the $R_{\rm AA}$ values is about a factor two
in central collisions at a $p_{\rm T}$ value of about $10\leavevmode\nobreak\
\mathrm{GeV}/c$ [30] and it is described by model calculations with parton-
mass-dependent energy loss. Therefore, for the centrality classes 0–10% and
30–50%, the value $R_{\rm AA}^{\textnormal{feed-down}}=2\times R_{\rm AA}^{\rm
prompt}$ was assumed to compute the $f_{\rm prompt}$ factor for D mesons with
$3<p_{\rm T}<24\leavevmode\nobreak\ \mathrm{GeV}/c$. To estimate a systematic
uncertainty, this hypothesis was varied in the range $1<R_{\rm
AA}^{\textnormal{feed-down}}/R_{\rm AA}^{\rm prompt}<3$ considering the data
uncertainties and model variations. For $p_{\rm T}<3\leavevmode\nobreak\
\mathrm{GeV}/c$ and $24<p_{\rm T}<50\leavevmode\nobreak\ \mathrm{GeV}/c$,
since model calculations predict a lower difference between $R_{\rm AA}^{\rm
prompt}$ and $R_{\rm AA}^{\textnormal{feed-down}}$ [63, 64], the hypothesis
$R_{\rm AA}^{\textnormal{feed-down}}=1.5\times R_{\rm AA}^{\rm prompt}$ was
used, with a variation in the range $1<R_{\rm AA}^{\textnormal{feed-
down}}/R_{\rm AA}^{\rm prompt}<2$ for estimating the systematic uncertainty.

The resulting values of the $f_{\rm prompt}$ factor, for the central
hypotheses on $R_{\rm AA}^{\textnormal{feed-down}}/R_{\rm AA}^{\rm prompt}$,
range from about 0.80 to 0.95, depending on the D-meson species, centrality
class, and $p_{\rm T}$ interval. The values of $f_{\rm prompt}$ for the ${\rm
D^{0}}$ analysis without decay-vertex reconstruction are larger compared to
those from the analysis with geometrical selections since the feed-down
component is not enhanced by the topological selection criteria. The
systematic uncertainties obtained from the variation of the hypotheses are
discussed in Section 4.

##  4 Systematic uncertainties

The systematic uncertainties on the corrected yields of ${\rm D^{0}}$, ${\rm
D^{+}}$, and ${\rm D^{*+}}$ mesons were estimated as a function of $p_{\rm T}$
for the 0–10% and 30–50% centrality classes considering the following sources:
(i) extraction of the raw yield from the invariant-mass distributions; (ii)
reconstruction efficiency of the decay-particle tracks; (iii) D-meson
selection efficiency; (iv) PID efficiency; (v) generated D-meson $p_{\rm T}$
and rapidity shape; (vi) subtraction of the contribution originating from
beauty-hadron decays. In addition, a global normalisation uncertainty due to
the branching-ratio uncertainty and the centrality interval determination was
considered. The estimated values of the systematic uncertainties are
summarised in Table 2 for the three D-meson species in representative $p_{\rm
T}$ intervals.

Table 2: Relative systematic uncertainties on prompt D-meson yields in Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$ for representative $p_{\rm T}$ intervals. The first $p_{\rm T}$ interval of the ${\rm D^{0}}$ corresponds to the analysis without decay vertex reconstruction. For the uncertainties on the correction factors, values below 0.5% are considered negligible. Particle | ${\rm D^{0}}$ | ${\rm D^{+}}$ | ${\rm D^{*+}}$
---|---|---|---
$p_{\rm T}$ interval ($\mathrm{GeV}/c$) | 0–1 | 1.5–2 | 7.5–8 | 3–3.5 | 7.5–8 | 3–3.5 | 7.5–8
0–10% centrality | Raw-yield extraction  | 9% | 7% | 2% | 5% | 5% | 7% | 3%
Correction factor |  |  |  |  |  |  |
Tracking efficiency | 8% | 9% | 9% | 14% | 15.5% | 13% | 9%
Selection efficiency | negl. | 8% | 2% | 5% | 5% | 13% | 4%
PID efficiency | negl. | negl. | negl. | negl. | negl. | 1% | 1%
MC $p_{\rm T}$ shape | negl. | 0.5% | negl. | negl. | negl. | 0.5% | negl.
Feed-down from beauty | ${}^{+2.7}_{-2.9}\%$ | ${}^{+7.6}_{-8.3}\%$ | ${}^{+11.6}_{-12}\%$ | ${}^{+3.4}_{-3.2}\%$ | ${}^{+7.1}_{-6.6}\%$ | ${}^{+8.4}_{-8.1}\%$ | ${}^{+10.3}_{-10.3}\%$
Branching ratio  | 0.8% | 1.7% | 1.1%
Centrality limits  | <0.1%
30–50% centrality | Raw-yield extraction  | 11% | 4% | 2% | 4% | 3% | 5% | 4%
Correction factor |  |  |  |  |  |  |
Tracking efficiency | 7% | 7.5% | 7% | 11% | 11% | 10% | 9%
Selection efficiency | negl. | 8% | 2% | 4% | 3% | 9% | 4%
PID efficiency | negl. | negl. | negl. | negl. | negl. | negl. | negl.
MC $p_{\rm T}$ shape | negl. | 1% | negl. | negl. | negl. | negl. | negl.
Feed-down from beauty | ${}^{+2.7}_{-2.9}\%$ | ${}^{+9.6}_{-10.4}\%$ | ${}^{+11.4}_{-11.7}\%$ | ${}^{+3.4}_{-3.2}\%$ | ${}^{+6.6}_{-6.3}\%$ | ${}^{+8.5}_{-8.2}\%$ | ${}^{+9.9}_{-9.9}\%$
Branching ratio  | 0.8% | 1.7% | 1.1%
Centrality limits  | 2%

The systematic uncertainty on the D-meson raw yield was estimated in each
$p_{\rm T}$ interval by varying the invariant-mass interval considered in the
fit and the functional form of the background fit function. The sensitivity to
the line shape of the signal was tested by considering the raw yield values
obtained by counting the candidates in the invariant-mass region of the signal
peak after subtracting the background estimated from the side bands. In the
case of the ${\rm D^{0}}$-meson analysis without decay-vertex reconstruction,
an additional contribution related to the line shape of the signal was
estimated by varying the width of the Gaussian function by $\pm 10\%$ with
respect to the Monte Carlo expectation, based on the deviations between the
Gaussian width values observed in data and simulations for the analysis with
decay-vertex reconstruction. The systematic uncertainty was defined as the RMS
of the distribution of the signal yields obtained from all these variations.
In the case of ${\rm D^{0}}$ mesons, an additional contribution due to the
modelling of the reflection contribution in the fit was estimated by varying
(by $\pm 20\%$) the ratio of the integral of the reflections to the integral
of the signal and the shape of the templates used in the invariant-mass fits.
The assigned uncertainty ranges from 2% to 11% depending on the D-meson
species and $p_{\rm T}$ interval, being on average smaller in the 30–50%
centrality class due to the larger signal-to-background ratio in this class as
compared to central collisions.

The systematic uncertainty on the track reconstruction efficiency has two
contributions. The first one was estimated by varying the track quality
selection criteria. The second contribution originates from possible
differences in the probability to match the TPC tracks to the ITS hits in data
and in simulations. It was estimated by comparing the matching efficiency in
data and simulations after weighting the relative abundances of primary and
secondary particles in the simulation to match those observed in data. The
uncertainty was estimated for the single track and propagated to the
reconstructed D mesons using the decay kinematics. The estimated uncertainty
is the quadratic sum of the two contributions and it depends on the D-meson
species and $p_{\rm T}$, ranging from 3% to 10.5% for the two-body decay of
${\rm D^{0}}$ mesons and from 5.5% to 15.5% for the three-body decays of ${\rm
D^{+}}$ and ${\rm D^{*+}}$ mesons.

The uncertainty on the D-meson selection efficiency originates from
imperfections in the description of the D-meson kinematic properties and of
the detector resolution and alignment in the simulation. For the analyses with
decay-vertex reconstruction, the systematic uncertainty was estimated by
comparing the corrected yields obtained by repeating the analysis with
different sets of selection criteria resulting in a significant modification
of the raw yields, signal-to-background ratios, and efficiencies. The
estimated uncertainty depends on the D-meson species and $p_{\rm T}$ interval,
being larger at low $p_{\rm T}$ and for central collisions, where more
stringent selections are used to obtain a good statistical significance of the
signal. The values obtained in these analyses range from 2% (3%) to 8% (13%)
for the two-(three-) body decay channel. In the case of the ${\rm
D^{0}}$-meson analysis without decay-vertex reconstruction, no geometrical
selections on the displaced decay-vertex topology are applied, and the
efficiencies are higher than those of the analysis with decay-vertex
reconstruction and almost independent of $p_{\rm T}$. The stability of the
corrected yield was tested against variations of the single-track $p_{\rm T}$
selection and no systematic effect was observed.

The uncertainty on the PID selection efficiency was estimated by repeating the
analyses with decay-vertex reconstruction without applying the PID selections.
The resulting corrected yields were found to be compatible with those obtained
with the PID selection and therefore no systematic uncertainty was assigned.
For the ${\rm D^{0}}$-meson analysis without decay-vertex reconstruction, the
analysis without applying PID selections could not be performed due to the
insufficient statistical significance of the signal. More stringent PID
criteria (at 2$\sigma$ level on TPC, or TOF or both) were tested and
compatible values for the corrected yields were obtained. Based on this result
and on the fact that the PID selections are the same as used in the analysis
with decay-vertex reconstruction, no uncertainty due to PID was assigned. In
the case of the ${\rm D^{*+}}$ analysis, more stringent PID selection criteria
were used in the 0–10% centrality class. A systematic uncertainty of 1% was
estimated by comparing the pion and kaon PID selection efficiencies in the
data and in the simulation and combining the observed differences using the
${\rm D^{*+}}$-meson decay kinematics. For this study, pure pion samples were
selected from strange-hadron decays, while pure kaon samples in the TPC (TOF)
were obtained using a tight PID selection in the TOF (TPC).

An additional contribution to the systematic uncertainty on the efficiency
originates from a possible difference between the real and simulated D-meson
$p_{\rm T}$ and rapidity distributions. The effect of the $p_{\rm T}$ shape
was estimated by calculating the efficiency using alternative D-meson
transverse momentum shapes via a reweighting technique. In particular, the
$p_{\rm T}$ distributions from FONLL calculations with and without hot-medium
effects parametrised based on the $R_{\rm AA}$ in central collisions from
different models were considered. For the analyses with decay-vertex
reconstruction, the resulting uncertainty, which also includes the effect of
the $p_{\rm T}$ dependence of the nuclear modification factor, was estimated
to be negligible for $p_{\rm T}>5\leavevmode\nobreak\ \mathrm{GeV}/c$ and to
increase to 1–1.5% in the lowest $p_{\rm T}$ intervals considered in the
analysis, where the efficiency varies steeply with $p_{\rm T}$. Instead, no
sensitivity to the generated ${\rm D^{0}}$ $p_{\rm T}$ shape was observed in
the results of the ${\rm D^{0}}$-meson analysis without decay-vertex
reconstruction. The simulated rapidity shape of D mesons affects the
$\mathrm{Acc}\times\epsilon$ and the $\alpha_{y}$ factors in the calculation
of the cross section. It was verified that the D-meson rapidity distributions
in the PYTHIA 8 simulations and FONLL calculations are similar, resulting in a
negligible effect on the $\mathrm{Acc}\times\epsilon$ and $\alpha_{y}$
correction factors. For the latter factor, an extreme assumption of a flat
rapidity shape was tested and the difference with respect to the FONLL case
was found to be smaller than 1% for $p_{\rm T}<$ 10 $\mathrm{GeV}/c$ and to be
smaller than 2% at higher $p_{\rm T}$. Considering that at high $p_{\rm T}$
the assumption of a flat rapidity shape is an extreme variation, the effect of
the generated ${\rm d}N/{\rm d}y$ shape was considered to be negligible and no
systematic uncertainty was assigned.

The systematic uncertainty on the $f_{\rm prompt}$ correction factor was
estimated by varying (i) the FONLL parameters (b-quark mass and factorisation
and renormalisation scales, according to the prescription in Ref. [60]) in the
calculation of the $p_{\rm T}$-differential production cross section of feed-
down D mesons, and (ii) the ratio between the feed-down and prompt D-meson
$R_{\rm AA}$, as described above. The resulting uncertainty ranges between 2%
and 14%, depending on the D-meson species and $p_{\rm T}$ and centrality
interval.

The normalisation uncertainty due to the centrality interval definition was
estimated from the variation of raw yield observed when varying the limits of
the centrality classes to account for the uncertainty on the fraction of the
hadronic cross section used in the Glauber fit to determine the centrality
percentiles (see Ref. [29] for details).

In the calculation of the nuclear modification factor, the systematic
uncertainties on the D-meson yield in Pb–Pb collisions and on the pp reference
cross section were propagated as uncorrelated, except for the uncertainty on
the BR, which cancels out in the ratio, and the contribution to the feed-down
uncertainty originating from the variation of the parameters of the FONLL
calculation, which was considered to be fully correlated between the Pb–Pb and
pp measurements. In particular, the contributions of the uncertainties on (i)
the normalisation of the pp cross section (due to the luminosity
determination), (ii) the centrality limits of the Pb–Pb samples, and (iii) the
$\left\langle T_{\rm AA}\right\rangle$ estimated with the Glauber model are
common to all the $p_{\rm T}$ intervals and therefore they constitute a
normalisation uncertainty on the $R_{\rm AA}$, which is shown separately from
the other sources when displaying the results.

##  5 Results

###  5.1 Transverse-momentum-differential yields

In this section, the results on the $p_{\rm T}$-differential (${\rm
d}^{2}N/{\rm d}y{\rm d}p_{\rm T}$) and $p_{\rm T}$-integrated (${\rm d}N/{\rm
d}y$) yields of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ mesons
at midrapidity in Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$ in the 0–10% and 30–50% centrality
classes are presented. In the case of the ${\rm D^{0}}$ mesons, two results
for the ${\rm d}^{2}N/{\rm d}y{\rm d}p_{\rm T}$ were obtained from the
analyses with and without selections on the decay-vertex topology. They are
compared in Fig. 3, with the inset showing their ratio in the common $p_{\rm
T}$ range. In the calculation of the ratio, the results in the narrower
$p_{\rm T}$ intervals of the analysis with decay-vertex reconstruction were
merged together to match the binning of the analysis without vertexing, and
the systematic uncertainties were propagated treating the raw yield extraction
as uncorrelated between the two analyses, and all the other sources of
uncertainty as correlated. The vertical error bars represent the statistical
uncertainties and the systematic uncertainties are depicted as boxes except
for the BR uncertainty, which is reported separately. The symbols,
representing the data points, are positioned horizontally at the center of
each interval and the horizontal bars represent the width of the $p_{\rm T}$
interval. This convention is adopted in all the figures shown in this section.

The two results for the $p_{\rm T}$-differential yield of prompt ${\rm D^{0}}$
mesons are found to be consistent within statistical uncertainties, which are
independent between the two analysis techniques because of the largely
different signal-to-background ratios and efficiencies. The usage of these two
techniques allows the measurement of the ${\rm D^{0}}$ $p_{\rm
T}$-differential yields in a wide transverse momentum range extending down to
$p_{\rm T}=$ 0 for the first time in Pb–Pb collisions. The most precise
measurement of the prompt ${\rm D^{0}}$-meson $p_{\rm T}$ spectrum, which will
be shown and used for comparisons throughout the paper, is obtained by using
the results of the analysis without decay-vertex reconstruction in the $p_{\rm
T}$ interval $0<p_{\rm T}<1\leavevmode\nobreak\ \mathrm{GeV}/c$ and those from
the analysis with decay-vertex reconstruction for $p_{\rm
T}>1\leavevmode\nobreak\ \mathrm{GeV}/c$.

Figure 3: Transverse momentum distributions ${\rm d}^{2}N/{\rm d}y{\rm
d}p_{\rm T}$ of prompt ${\rm D^{0}}$ mesons from the analysis with (open
marker) and without (full marker) decay-vertex reconstruction in the 0–10% and
30–50% centrality classes in Pb–Pb collisions at
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$.
Statistical uncertainties (bars) and systematic uncertainties (boxes) are
shown.

The $p_{\rm T}$-differential yields ${\rm d}^{2}N/{\rm d}y{\rm d}p_{\rm T}$ of
prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ mesons are shown in
Fig. 4 for the 0–10% and 30–50% centrality classes. They are compared with the
reference yields from pp collisions, which are computed as $\langle T_{\rm
AA}\rangle\times{\rm d}\sigma_{\rm pp}/{\rm d}p_{\rm T}$, where ${\rm
d}\sigma_{\rm pp}/{\rm d}p_{\rm T}$ is the D-meson $p_{\rm T}$-differential
cross section measured in pp collisions at $\sqrt{s}=5.02\leavevmode\nobreak\
\mathrm{TeV}$ [37, 65], and $\langle T_{\rm AA}\rangle$ is the average nuclear
overlap function [47]. The D-meson production cross sections in pp collisions
are measured up to 36 $\mathrm{GeV}/c$ and they are extrapolated with FONLL
towards higher $p_{\rm T}$, with the method described in Refs. [30, 29]. The
spectra in the 30–50% centrality class are scaled by the factor reported in
the legend for visibility. A clear suppression of the production yield for
$p_{\rm T}>3\leavevmode\nobreak\ \mathrm{GeV}/c$ is visible in Pb–Pb
collisions and it is stronger in central than in semicentral collisions.

Figure 4: Transverse momentum distributions ${\rm d}^{2}N/{\rm d}y{\rm
d}p_{\rm T}$ of prompt ${\rm D^{0}}$ (left), ${\rm D^{+}}$ (middle), and ${\rm
D^{*+}}$ (right) mesons in the 0–10% (cross) and 30–50% (diamond) centrality
classes in Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$. The reference pp distributions
multiplied by $\langle T_{\rm AA}\rangle$ are shown as well. Statistical
uncertainties (bars) and systematic uncertainties (boxes) are shown. The
uncertainties on the BRs are quoted separately and the horizontal bars
represent bin widths. The spectra in the 30–50% centrality class are scaled by
the factor reported in the legend for visibility. Filled and empty markers of
the pp reference indicate the measured [37, 65] and $p_{\rm T}$-extrapolated
values, respectively, of the $p_{\rm T}$-spectrum. Figure 5: ${\rm D^{+}}/{\rm
D^{0}}$ (top left panel) and ${\rm D^{*+}}/{\rm D^{0}}$ (bottom left panel)
ratios as a function of $p_{\rm T}$ in central and semicentral Pb–Pb
collisions compared to the results obtained from pp collisions [37]. The right
panels show the ratios in central Pb–Pb collisions compared to the predictions
from the statistical hadronisation model (SHMc) [66, 67] with a core–corona
approach and the FastReso package [68] for resonance decays. Statistical
(bars) and systematic (boxes) uncertainties are shown.

Figure 5 shows the $p_{\rm T}$-dependent ratios of the production yields of
prompt ${\rm D^{+}}$ and ${\rm D^{*+}}$ mesons relative to the ${\rm D^{0}}$
ones, in the left panels, compared with the values measured in pp collisions
at $\sqrt{s}=5.02\leavevmode\nobreak\ \mathrm{TeV}$ [37]. The systematic
uncertainties were propagated to the ratios, considering the contributions
from the tracking efficiency and the beauty-hadron feed-down subtraction as
fully correlated among the different D-meson species. The values obtained in
Pb–Pb and pp collisions are compatible within uncertainties and this indicates
there is no modification of the relative abundances of ${\rm D^{0}}$, ${\rm
D^{+}}$, and ${\rm D^{*+}}$ mesons as a function of $p_{\rm T}$ and for
different centrality classes. The ${\rm D^{+}}/{\rm D^{0}}$ and ${\rm
D^{*+}}/{\rm D^{0}}$ ratios are described within uncertainties by the
GSI–Heidelberg statistical hadronisation model (SHMc) [66, 67] (right panels
of Fig. 5), which predicts the $p_{\rm T}$ spectra of charm hadrons with a
core–corona approach. The low-$p_{\rm T}$ region is dominated by the core
contribution described with a blast-wave function, while the corona
contribution is more relevant at high $p_{\rm T}$ and is parametrised from pp
measurements. The rise predicted by the model, especially for ${\rm
D^{*+}}/{\rm D^{0}}$, are due to the different D-meson masses and to the
collective radial expansion of the system.

The visible production yields of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and
${\rm D^{*+}}$ mesons in the two centrality classes were evaluated by
integrating the ${\rm d}^{2}N/{\rm d}y{\rm d}p_{\rm T}$ over the $p_{\rm T}$
intervals of the measurements. The results are reported in Table 3. For the
yield integration, the systematic uncertainty was propagated assuming all
sources of uncertainty as fully correlated among $p_{\rm T}$ intervals, except
for the one on the raw-yield extraction, which was treated as uncorrelated
because of the bin-to-bin variations of $S/B$ and background invariant-mass
shape.

The production yield of prompt D mesons per unit of rapidity, ${\rm d}N/{\rm
d}y$, was obtained by extrapolating (where needed) the visible cross section
to the full $p_{\rm T}$ range. Since the $p_{\rm T}$-differential yields of
${\rm D^{0}}$ mesons in the 0–10% and 30–50% centrality classes were measured
down to $p_{\rm T}=$ 0, it was possible to obtain the ${\rm d}N/{\rm d}y$ at
midrapidity for the first time in Pb–Pb collisions at the LHC without using
models or assumptions to extrapolate the $p_{\rm T}$-differential yield to
unmeasured momentum range at low $p_{\rm T}$. Conversely, the ${\rm d}N/{\rm
d}y$ of ${\rm D^{+}}$ and ${\rm D^{*+}}$ mesons were obtained via an
extrapolation procedure exploiting the measured $p_{\rm T}$-differential
ratios relative to the ${\rm D^{0}}$ meson. The measured $p_{\rm
T}$-differential ${\rm D^{+}}/{\rm D^{0}}$ and ${\rm D^{*+}}/{\rm D^{0}}$
ratios, shown in Fig. 5, were fit with two different functions (a constant and
a constant plus logarithmic function) in order to extrapolate the ratios in
the interval $0<p_{\rm T}<2\leavevmode\nobreak\ \mathrm{GeV}/c$, where the
yields could not be measured. The extrapolation in the high-$p_{\rm T}$ range
is negligible because the fraction of D-meson yield in central (semicentral)
collisions with $p_{\rm T}>$ 50 (36)$\leavevmode\nobreak\ \mathrm{GeV}/c$ is
negligible. The yields of ${\rm D^{+}}$ and ${\rm D^{*+}}$ mesons in the
interval $0<p_{\rm T}<2\leavevmode\nobreak\ \mathrm{GeV}/c$ were obtained by
multiplying the measured ${\rm D^{0}}$ yield in this $p_{\rm T}$ interval by
the extrapolated ratio from the fit function. For the ${\rm D^{+}}/{\rm
D^{0}}$ ratios, which show a flat trend as a function of $p_{\rm T}$, the fit
to a constant was used in the determination of the central value of the ${\rm
D^{+}}$ yield, while for the ${\rm D^{*+}}/{\rm D^{0}}$ ratio the fit with a
constant plus a logarithmic function was used. The systematic uncertainty due
to the extrapolation was estimated by evaluating the yields using the other
fit function and the shape of the SHMc predictions with normalisation
parameter left free. The difference between the central value and the yields
obtained with other functions for the extrapolation was assigned as systematic
uncertainty and summed in quadrature to the uncertainty on the ${\rm D^{0}}$
yield in the aforementioned $p_{\rm T}$ interval. The results for the prompt
D-meson production yields per unit of rapidity ${\rm d}N/{\rm d}y$ in $|y|<$
0.5 are reported in Table 4. In the case of the ${\rm D^{0}}$ meson, the
result coincides with the visible yield because the measurement extends down
to $p_{\rm T}=$ 0 and the contribution of ${\rm D^{0}}$ mesons with $p_{\rm
T}>$ 50 (36)$\leavevmode\nobreak\ \mathrm{GeV}/c$ is negligible.

Table 3: Visible production yield of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ mesons in $|y|<$ 0.5 in the 0–10% and 30–50% centrality classes of Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$. | Kinematic range | Visible production yield
---|---|---
0–10% centrality
${\rm D^{0}}$ | 0 $<p_{\rm T}<$ 50 $\mathrm{GeV}/c$ | 6.819 $\pm$ 0.457 (stat.) ${}^{+0.912}_{-0.936}$ (syst.) $\pm$ 0.054 (BR)
${\rm D^{+}}$ | 2 $<p_{\rm T}<$ 50 $\mathrm{GeV}/c$ | 0.992 $\pm$ 0.073 (stat.) ${}^{+0.154}_{-0.155}$ (syst.) $\pm$ 0.017 (BR)
${\rm D^{*+}}$ | 3 $<p_{\rm T}<$ 50 $\mathrm{GeV}/c$ | 0.438 $\pm$ 0.037 (stat.) ${}^{+0.084}_{-0.085}$ (syst.) $\pm$ 0.005 (BR)
30–50% centrality
${\rm D^{0}}$ | 0 $<p_{\rm T}<$ 50 $\mathrm{GeV}/c$ | 1.275 $\pm$ 0.099 (stat.) ${}^{+0.167}_{-0.173}$ (syst.) $\pm$ 0.010 (BR)
${\rm D^{+}}$ | 2 $<p_{\rm T}<$ 50 $\mathrm{GeV}/c$ | 0.179 $\pm$ 0.008 (stat.) ${}^{+0.024}_{-0.024}$ (syst.) $\pm$ 0.003 (BR)
${\rm D^{*+}}$ | 2 $<p_{\rm T}<$ 36 $\mathrm{GeV}/c$ | 0.230 $\pm$ 0.023 (stat.) ${}^{+0.038}_{-0.039}$ (syst.) $\pm$ 0.002 (BR)
Table 4: Measured $p_{\rm T}$-integrated yields of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ in $|y|<$ 0.5 in the 0–10% and 30–50% centrality classes of Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$. The right column reports the ${\rm d}N/{\rm d}y$ predicted by the GSI–Heidelberg statistical hadronisation model [66, 67]. | Measured ${\rm d}N/{\rm d}y$ | SHMc ${\rm d}N/{\rm d}y$
---|---|---
0–10% centrality
${\rm D^{0}}$ | 6.819 $\pm$ 0.457 (stat.) ${}^{+0.912}_{-0.936}$ (syst.) $\pm$ 0.054 (BR) | 6.42 $\pm$ 1.07
${\rm D^{+}}$ | 3.041 $\pm$ 0.073 (stat.) ${}^{+0.154}_{-0.155}$ (syst.) $\pm$ 0.052 (BR) ${}^{+0.352}_{-0.618}$ (extrap.) | 2.84 $\pm$ 0.47
${\rm D^{*+}}$ | 3.803 $\pm$ 0.037 (stat.) ${}^{+0.084}_{-0.085}$ (syst.) $\pm$ 0.041 (BR) ${}^{+0.854}_{-1.175}$ (extrap.) | 2.52 $\pm$ 0.42
30–50% centrality
${\rm D^{0}}$ | 1.275 $\pm$ 0.099 (stat.) ${}^{+0.167}_{-0.173}$ (syst.) $\pm$ 0.010 (BR) | 1.06 $\pm$ 0.15
${\rm D^{+}}$ | 0.552 $\pm$ 0.008 (stat.) ${}^{+0.024}_{-0.024}$ (syst.) $\pm$ 0.009 (BR) ${}^{+0.068}_{-0.114}$ (extrap.) | 0.471 $\pm$ 0.069
${\rm D^{*+}}$ | 0.663 $\pm$ 0.023 (stat.) ${}^{+0.038}_{-0.039}$ (syst.) $\pm$ 0.007 (BR) ${}^{+0.149}_{-0.165}$ (extrap.) | 0.419 ${}^{+0.065}_{-0.061}$

In the right column the ${\rm d}N/{\rm d}y$ values predicted by the SHMc [66,
67] are reported. In the SHMc approach, the yield of hadrons containing charm
quarks can be calculated utilising as input values (i) the temperature
$T_{\mathrm{chem}}$, (ii) the volume of the fireball at the chemical freeze-
out estimated from a fit to light-flavour hadron yields, and (iii) the number
of ${\rm c\overline{c}}$ pairs produced in the Pb–Pb collision. The latter was
calculated from the charm-quark production cross section estimated from
measurements in pp and p–Pb collisions along with guidance from
parameterisations of the nuclear modification of the PDFs [69]. The SHMc
calculations agree with the measurements within uncertainties, even though the
data lie on the upper edge of the uncertainty band of the theoretical
predictions, which is related to the uncertainty on the total charm production
cross section used in the calculation. It is useful to note that the charm
production cross section in pp collisions, which was used to determine the
charm content of the fireball in the SHMc calculations, is lower than the
recent measurement reported in Ref. [70] and therefore the yields from the
SHMc would increase if the measured cross section were used as the input to
the calculations. However, a firmer conclusion on the SHMc predictions for
charm hadrons will be drawn when the measured ${\rm d}N/{\rm d}y$ of ${\rm
D_{s}^{+}}$, $\Lambda^{+}_{\rm{c}}$, and ${\rm J}/\psi$ will be available in
Pb–Pb collisions and included in the comparison, since it is expected that the
modifications of the hadronisation mechanisms in the presence of a QGP affect
the relative abundances of different charm-hadron species [71, 25, 72, 73, 74,
75, 76].

###  5.2 Nuclear modification factor

The $R_{\rm AA}$ of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$
mesons was computed, as defined in Eq. 1, using the ${\rm d}^{2}N/{\rm d}y{\rm
d}p_{\rm T}$ measured in Pb–Pb collisions and the pp reference at the same
centre of mass energy reported in Fig. 4. The obtained results are shown in
Fig. 6 for the two centrality classes. The nuclear modification factors of the
three D-meson species are compatible among each other within statistical
uncertainties.

Figure 6: $R_{\rm AA}$ of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm
D^{*+}}$ mesons as a function of $p_{\rm T}$ for the 0–10% (left panel) and
30–50% (right panel) centrality classes. Statistical (bars), systematic
(boxes), and normalisation (shaded box around unity) uncertainties are shown.
Filled markers are obtained with the measured pp reference [37, 65], empty
markers with the $p_{\rm T}$-extrapolated reference.

The average $R_{\rm AA}$ of ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$
mesons was computed as a weighted average using the inverse of the quadratic
sum of the relative statistical and uncorrelated systematic uncertainties as
weights, in the $p_{\rm T}$ intervals where more than one D-meson $R_{\rm AA}$
value is available. The systematic uncertainties due to the raw-yield
extraction and selection efficiency were considered as uncorrelated among
different D-meson species and therefore they were used in the definition of
the weights and propagated through the weighted average as uncorrelated, while
the other sources of uncertainty (tracking efficiency, generated $p_{\rm T}$
shape, and beauty hadron feed-down) were treated as fully correlated.

The prompt D-meson average nuclear modification factors in the 0–10% and
30–50% centrality classes are shown in Fig. 7 together with the $R_{\rm AA}$
in the 60–80% centrality class taken from Ref. [30], which was measured using
the sample of Pb–Pb collisions collected in 2015. The suppression increases
from peripheral to central collisions. The $R_{\rm AA}$ shows a minimum value
at $p_{\rm T}=$ 6–8 $\mathrm{GeV}/c$, corresponding to a suppression of the
yields by a factor 5 and 2.5 with respect to the binary-scaled pp reference in
the 0–10% and 30–50% classes, respectively. The stronger suppression observed
in central collisions is due to the increasing medium density, size, and
lifetime of the fireball from peripheral to central collisions. Also shown in
Fig. 7 is the nuclear modification factor $R_{\rm pPb}$ measured in p–Pb
collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\
\mathrm{TeV}$ taken from Ref. [33], which is compatible with unity within
uncertainties, confirming that the suppression observed in Pb–Pb collisions is
to due final-state effects induced by the formation of a hot and dense QGP
medium.

Figure 7: Left panel: prompt D-meson $R_{\rm AA}$ (average of ${\rm D^{0}}$,
${\rm D^{+}}$, and ${\rm D^{*+}}$) as a function of $p_{\rm T}$ measured in
Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$ (2018 data sample) in the 0–10%
and 30–50% centrality classes compared with published results in the 60–80%
centrality class (2015 data sample) [30] and in p–Pb collisions at the same
centre-of-mass energy [33]. Statistical (bars), systematic (boxes), and
normalisation (shaded box around unity) uncertainties are shown. Right panel:
prompt D-meson $R_{\rm AA}$ in the 10% most central Pb–Pb collisions at
$\sqrt{s_{\rm\scriptscriptstyle NN}}=\leavevmode\nobreak\ $5.02 TeV and 2.76
TeV [29] compared to the ${\rm D^{0}}$ $R_{\rm AA}$ measured by the STAR
collaboration in Au–Au collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=\leavevmode\nobreak\ $200 GeV [32].

The right panel of Fig. 7 shows the $R_{\rm AA}$ of prompt D mesons measured
by the ALICE collaboration in Pb–Pb collisions at
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$
and 2.76 TeV [29] and of ${\rm D^{0}}$ mesons measured at RHIC by the STAR
collaboration in Au–Au collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=\leavevmode\nobreak\ $200 GeV [32]. The results from ALICE at the two
different energies are compatible within uncertainties. This similarity
originates from the interplay between the higher medium temperature and
density at 5.02 TeV, which causes a larger energy loss and a lower $R_{\rm
AA}$, and the harder $p_{\rm T}$ distribution of charm quarks at 5.02 TeV,
which would increase the $R_{\rm AA}$ if the medium temperature were the same
at the two collision energies [63]. The nuclear modification factor of ${\rm
D^{0}}$ mesons at RHIC energies shows a trend with $p_{\rm T}$ similar to the
one measured at the LHC, with a possible hint for a smaller $R_{\rm AA}$ at
low $p_{\rm T}$ ($<2\leavevmode\nobreak\ \mathrm{GeV}/c$) and a larger $R_{\rm
AA}$ at high $p_{\rm T}$ ($>4\leavevmode\nobreak\ \mathrm{GeV}/c$) in
collisions at $\sqrt{s_{\rm\scriptscriptstyle NN}}=\leavevmode\nobreak\ $200
GeV compared to LHC energies. This difference could be due to dependence with
collision energy of the charm-quark $p_{\rm T}$ distributions, the
initial-/final-state effects, and the medium properties. However, the rather
large uncertainties prevent a firm conclusion on a
$\sqrt{s_{\rm\scriptscriptstyle NN}}$ dependence of the $R_{\rm AA}$ from
being drawn.

In addition, the $p_{\rm T}$-integrated $R_{\rm AA}$ of prompt ${\rm D^{0}}$
mesons at midrapidity was calculated from the $p_{\rm T}$-integrated yield in
Pb–Pb collisions reported above and the ${\rm D^{0}}$ production cross section
measured in pp collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$ [37]. The results in the 0–10% and
30–50% centrality classes are:

$\begin{split}R^{{\mathrm{prompt\leavevmode\nobreak\
D^{0}}}}_{{\mathrm{AA}}}(0\text{--}10\%)=0.689\pm 0.054\leavevmode\nobreak\
({\mathrm{stat.}})^{+0.104}_{-0.106}\leavevmode\nobreak\
({\mathrm{syst.}}),\\\ R^{{\mathrm{prompt\leavevmode\nobreak\
D^{0}}}}_{{\mathrm{AA}}}(30\text{--}50\%)=0.775\pm 0.069\leavevmode\nobreak\
({\mathrm{stat.}})^{+0.117}_{-0.120}\leavevmode\nobreak\
({\mathrm{syst.}}).\end{split}$ (4)

Figure 8 shows the results obtained in Pb–Pb collisions compared with the
nuclear modification factor $R_{\rm pPb}$ measured in p–Pb collisions at the
same centre-of-mass energy [33]. The total charm cross section is expected to
scale with the number of binary collisions $N_{\rm coll}$, as introduced in
Section 1, and thus the $R_{\rm AA}$ should be equal to one. However, the
nuclear shadowing effect reduces the charm production in Pb–Pb (and p–Pb )
collisions with respect to pp interactions. In addition, the possible enhanced
production of ${\rm D_{s}^{+}}$ and $\Lambda^{+}_{\mathrm{c}}$ due to the
hadronisation via recombination is expected to further decrease the fraction
of charm quarks that hadronise into ${\rm D^{0}}$ mesons in Pb–Pb collisions
compared to pp collisions [71, 25, 72, 73, 74, 75, 76]. The measured $p_{\rm
T}$-integrated $R_{\rm AA}$ is significantly below unity and this confirms the
suppression of the ${\rm D^{0}}$-meson yield in Pb–Pb collisions with respect
to the binary-scaled pp reference due to shadowing and the possible
modifications in the hadronisation mechanism. Conversely, the $p_{\rm
T}$-integrated ${\rm D^{0}}$ $R_{\rm pPb}$ is closer to unity, as expected
from the smaller shadowing effects in p–Pb compared to Pb–Pb collisions (where
it affects the nucleons of both the projectile and the target nuclei). The
integrated $R_{\rm AA}$ is also compared with perturbative QCD calculations of
${\rm D^{0}}$-meson production including only initial-state effects modeled
using two different sets of nuclear PDFs, namely nCTEQ15 [77, 78, 79, 80, 81]
and EPPS16 [82, 83]. The calculations with EPPS16 do not include the
dependence of the shadowing on the impact parameter of the Pb–Pb collision and
therefore they are the same in the central and semicentral event classes. The
predictions with nCTEQ15 are obtained applying a Bayesian reweighting of the
nuclear PDFs, which is constrained by measurements of heavy-flavour production
in p–Pb collisions at the LHC [78], and are labelled as nCTEQ15rwHF in Fig. 8.
They include a modelling of the centrality dependence of the nuclear
modification of the PDFs. The uncertainty bands represent a 90% confidence
level uncertainty. In the nCTEQ15 case they are determined by considering
three different factorisation scales in addition to the PDF uncertainties,
with the scale variation constituting the main source of uncertainty, as
described in Ref. [78]. Both model calculations include only the initial-state
effects due to the nuclear PDFs, but not the possible modifications of the
relative abundances of different charm hadron species due to hadronisation via
recombination. The measured $R_{\rm AA}$ values lie on the upper edge of the
theoretical predictions and this could be due to a smaller shadowing effect in
the data with respect to the model expectations.

Figure 8: $p_{\rm T}$-integrated nuclear modification factors of prompt ${\rm
D^{0}}$ mesons measured in p–Pb [33] and Pb–Pb collisions at
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$.
Statistical (bars) and systematic (brackets) uncertainties are shown. The
results are compared with calculations at 90% of confidence level of
theoretical models nCTEQ15 (with Bayesian reweighting, see text for details)
[77, 78, 79, 80, 81] and EPPS16 [82, 83] that include only initial-state
effects.

###  5.3 Discussion: energy loss regime

The comparison of the nuclear modification factor of prompt D mesons with that
of pions and of particles originating from beauty-hadron decays can provide
essential insights into the characteristics of the in-medium parton energy
loss, in particular on its predicted dependence on the colour charge and the
quark mass. To investigate possible differences with respect to light-flavour
particles, the average $R_{\rm AA}$ of prompt ${\rm D^{0}}$, ${\rm D^{+}}$,
and ${\rm D^{*+}}$ mesons in the 0–10% centrality class is compared in Fig. 9
with that of charged pions [84] and charged particles [85]. The charged-
particle $R_{\rm AA}$ is shown for $p_{\rm T}>20\leavevmode\nobreak\
\mathrm{GeV}/c$ in order to extend the comparison to light-flavour particles
up to the highest $p_{\rm T}$ interval of the D-meson measurement. The $R_{\rm
AA}$ of charged particles can be used in this comparison at high $p_{\rm T}$
in place of the pion one because for $p_{\rm T}>$ 8–10$\leavevmode\nobreak\
\mathrm{GeV}/c$. In this range, the nuclear modification factors of different
light-flavour hadron species are found to be consistent with each other and
the particle composition is compatible with that measured in pp collisions,
that is dominated by pions [86, 84].

Figure 9: Average $R_{\rm AA}$ of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and
${\rm D^{*+}}$ mesons in the 0–10$\%$ centrality class compared to $R_{\rm
AA}$ of charged pions [84], charged particles [85], inclusive J/$\psi$
measured by ALICE [87], and of prompt and non-prompt J/$\psi$ from CMS [61] in
Pb–Pb collisions at $\sqrt{s_{\rm\scriptscriptstyle
NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$.

The $R_{\rm AA}$ of D mesons is larger than that of pions for $p_{\rm
T}<8\leavevmode\nobreak\ \mathrm{GeV}/c$, providing a clear evidence for a
different nuclear modification factor of D mesons and light-flavour particles
at low and intermediate $p_{\rm T}$. This difference originates from the
interplay of several effects that concur in determining the magnitude and the
$p_{\rm T}$-dependence of the $R_{\rm AA}$, so the intepretation in terms of
different in-medium parton energy loss of charm quarks, light quarks, and
gluons is not straighforward. As pointed out in Ref. [19], also in the
presence of a colour-charge and quark-mass dependent energy loss, similar
values of D-meson and pion $R_{\rm AA}$ are expected at high $p_{\rm T}$
($\gtrsim$ 8$\leavevmode\nobreak\ \mathrm{GeV}/c$) due to the harder $p_{\rm
T}$ distribution and the harder fragmentation function of charm quarks
compared to those of light quarks and gluons. At low $p_{\rm T}$, where a
large difference is observed between the D-meson and pion $R_{\rm AA}$’s, it
should be considered that the pion yield can have a significant contribution
from soft production up to transverse momenta of about
3–4$\leavevmode\nobreak\ \mathrm{GeV}/c$ due to the strong radial flow at LHC
energies. This soft component, which is not present in the D-meson production,
does not scale with the number of binary nucleon–nucleon collisions. In
addition, effects due to the radial flow and the hadronisation can affect
D-meson and light-hadron yields differently at a given $p_{\rm T}$, playing a
role in the interpretation of their different $R_{\rm AA}$ at low and
intermediate $p_{\rm T}$. For instance, an enhanced production of ${\rm
D_{s}^{+}}$ mesons and charm baryons in Pb–Pb collisions due to recombination
would imply a reduction of the fraction of charm quarks hadronising into non-
strange meson species compared to pp collisions.

A quantitative understanding of the parton in-medium energy-loss from the
$R_{\rm AA}$ measurements needs therefore to be carried out via comparisons
with model calculations. In particular, in the high-$p_{\rm T}$ region, where
effects due to radial flow and modifications in the hadronisation mechanisms
are expected to be negligible, models based on perturbative QCD (pQCD)
calculations of high-$p_{\rm T}$ parton energy loss are expected to describe
the data. The $R_{\rm AA}$ and the elliptic flow $v_{2}$ of prompt D mesons
(taken from Ref. [36]) and pions (taken from Ref. [84, 88]) in the 0–10% and
30-50% centrality classes are compared in Fig. 10 to three of these models,
namely CUJET3.1 [89, 90, 91], DREENA-A [92, 93, 94], and ${\rm SCET_{M,G}}$
[95]. The DREENA-A is a numerical framework that calculates the medium-
modified distribution of high-$p_{\rm T}$ hadrons based on a dynamical energy-
loss formalism coupled with a modelling of initial parton momentum
distributions, a full 3+1D hydrodynamic evolution of the medium, and
fragmentation functions. The dynamical energy-loss formalism [96, 97, 98] is a
model of jet-medium interactions incorporating collisional and radiative
energy loss mechanisms in a QCD medium of finite size and temperature composed
of dynamical scattering centers. It includes a colour-charge and quark-mass
dependence as well as finite magnetic mass effects and running strong coupling
constant. The CUJET3.1 framework provides a calculation of jet energy loss in
a QCD medium described with relativistic viscous hydrodynamics. The jet-medium
interactions are based on the DGLV opacity expansion model [99, 100, 22]
including both inelastic and elastic scatterings with their colour-charge and
quark-mass dependence, and taking into account interactions with both chromo-
electric and magnetic charges of the medium. The ${\rm SCET_{\rm M,G}}$ model
implements medium-induced gluon radiation via modified splitting functions. It
is based on a soft-collinear effective theory (SCET [101, 102]) describing the
parton shower formed in the vacuum via soft and collinear splittings. This
effective theory was extended to the case of massless [103, 104] and massive
[95] quarks propagating in strongly-interacting matter by including additional
splitting processes induced by the interactions of the incident parton with
the QCD medium mediated by Glauber gluon exchange. These medium-induced
splittings are calculated to first order in the opacity series expansion, thus
limiting the applicability of these calculations to the high $p_{\rm T}$
region.

The considered models provide a fair description of both the $R_{\rm AA}$ and
the $v_{2}$ of D mesons and pions for $p_{\rm T}>10\leavevmode\nobreak\
\mathrm{GeV}/c$, where radiative energy loss is expected to be the dominant
interaction mechanism. This suggests that the dependences of radiative energy
loss on the colour charge and the quark mass of the hard-scattered parton, as
well as on the path length in the hot and dense medium are reasonably well
described in these calculations.

Figure 10: Average $R_{\rm AA}$ (top panels) and $v_{2}$ (bottom panels) of
prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ mesons [36] and pions
(charged hadrons) [84, 85, 88] compared with predictions from models based on
pQCD calculations, namely CUJET3.1 [89, 90, 91], DREENA-A [92], and ${\rm
SCET_{\rm M,G}}$ [95].

The comparison of the nuclear modification factors of particles originating
from charm and beauty quarks is also reported in Fig. 9, to provide insight
into the predicted quark-mass dependence of parton energy loss. In particular,
the $R_{\rm AA}$ of inclusive ${\rm J}/\psi$ mesons measured by ALICE [87] at
midrapidity in the 0–20% centrality class and those of prompt and non-prompt
${\rm J}/\psi$ from the CMS collaboration [61] are shown.

The $R_{\rm AA}$ of prompt D mesons in the 0–10% centrality class is lower
than that of non-prompt ${\rm J}/\psi$ mesons from beauty hadron decays
measured in the same centrality interval. This difference provides an
indication for the predicted quark-mass dependence of in-medium energy loss,
even though a proper interpretation of the different nuclear modification
factors in terms of energy loss requires a full modelling of the initial
momentum distributions of charm and beauty quarks, of the heavy-quark
fragmentation, and of the beauty-hadron decay kinematics. The model from Ref.
[105], which includes all these effects, provides a good description of the
measurements of prompt D and non-prompt ${\rm J}/\psi$ mesons. In this model,
the large difference in the $R_{\rm AA}$’s in the $p_{\rm T}$ interval around
10 $\mathrm{GeV}/c$ is mostly due to the different in-medium energy loss of c
and b quarks, as the effects of initial distributions, quark fragmentation,
and beauty-hadron decays are found to be small.

The comparison of the nuclear modification factors of prompt D mesons and
prompt (inclusive) ${\rm J}/\psi$ mesons, also shown in Fig. 9, provides
insight into the interplay of different QGP medium effects in the charm
sector, where they are expected to affect differently the production of open
charm and charmonium states. At high momentum ($p_{\rm T}\gtrsim
10\leavevmode\nobreak\ \mathrm{GeV}/c$), the $R_{\rm AA}$ of prompt D mesons
is compatible within uncertainties with that of prompt ${\rm J}/\psi$ mesons,
as well as with that of light-flavour hadrons. This may suggest that in this
$p_{\rm T}$ region the yield of charmonia has a significant contribution from
production within the parton shower originating from the splitting of a hard-
scattered gluon [106, 107], which experiences in-medium energy loss before
fragmenting [108, 109]. At lower $p_{\rm T}$ ($2\lesssim p_{\rm T}\lesssim
4\leavevmode\nobreak\ \mathrm{GeV}/c$), the results for inclusive ${\rm
J}/\psi$, which are dominated by the prompt contribution, show a magnitude and
a trend of the $R_{\rm AA}$ similar to that of D mesons. In this region, ${\rm
J}/\psi$ production is likely dominated by recombination of c and
$\overline{\rm c}$ quarks in the medium either at the phase boundary or during
the QGP expansion [110, 87], similarly to the recombination of charm and light
quarks forming D mesons. The similar $R_{\rm AA}$ of ${\rm J}/\psi$ and D
mesons in this momentum region may signal the dominant contribution of
hadronisation via recombination after the interactions of the charm quarks
with the QGP medium constituents [111, 112, 113, 114], with possible small
differences arising from the different kinematics involved in the charm-quark
coalescence with $\overline{\rm c}$ and light antiquarks.

###  5.4 Discussion: transport models and an estimate of the spatial diffusion
coefficient

As discussed in Section 1, the charm-hadron yields and angular distributions
at low and intermediate $p_{\rm T}$ are sensitive to the diffusion and the
possible thermalisation of charm quarks in the medium. The comparison of the
measured D-meson $R_{\rm AA}$ and $v_{2}$ to models implementing charm-quark
transport in a hydrodynamically expanding QGP could therefore provide insight
into the interactions of heavy quarks with the medium, constraining in
particular the spatial diffusion coefficient, which is the relevant transport
coefficient in the diffusion regime. In the left panels of Fig. 11 the
measured nuclear modification factor of prompt D mesons for the 0–10$\%$ (top)
and 30–50$\%$ (bottom) centrality class, respectively, is compared with
various predictions from models implementing charm-quark transport in a
hydrodynamically expanding medium [115, 116, 117, 118, 119, 120, 76, 121, 122,
123, 124, 125]. In particular, the TAMU [76], POWLANG-HTL [117, 118], PHSD
[125], and Catania [122, 123] models describe the interactions of the charm
quarks with the medium constituents solely via collisional processes, while
the MC@sHQ+EPOS2 [115], DAB-MOD [116], LBT [119, 120], LGR [121], and LIDO
[124] calculations include also radiative processes. All the models, except
for DAB-MOD, include initial-state effects by using nuclear PDFs (nPDFs) in
the calculation of the initial $p_{\rm T}$ distributions of charm quarks. A
contribution of hadronisation via quark recombination, in addition to charm-
quark fragmentation, is included in all theoretical predictions. The
theoretical uncertainties, where available, are displayed with a coloured
band.

Figure 11: Average $R_{\rm AA}$ (left) and elliptic flow $v_{2}$ [36] (right)
of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ mesons in the
0–10$\%$ (top) and 30–50$\%$ (bottom) centrality classes compared with
predictions of models implementing the charm-quark transport in a
hydrodynamically expanding medium [115, 116, 117, 118, 119, 120, 76, 121, 122,
123, 124, 125].

Most of the models capture the measured $p_{\rm T}$ trend and magnitude of the
nuclear modification factor and provide a good description of the data in the
$p_{\rm T}$ interval 6 $<p_{\rm T}<$ 10$\leavevmode\nobreak\ \mathrm{GeV}/c$
where the minimum of the $R_{\rm AA}$ is observed, but many of them show
significant deviations from the measurements at low $p_{\rm T}$ ($p_{\rm
T}\lesssim$ 4–5$\leavevmode\nobreak\ \mathrm{GeV}/c$). In particular, the LBT
[119, 120], LIDO [124], MC@sHQ+EPOS2 [115], and Catania [122, 123]
calculations tend to underestimate the results for 2 $<p_{\rm T}<$
4$\leavevmode\nobreak\ \mathrm{GeV}/c$ in both centrality classes, while the
PHSD [125] prediction is above the measured $R_{\rm AA}$ in semicentral
collisions for the same $p_{\rm T}$ range. The POWLANG-HTL [117, 118]
calculation, instead, predicts a narrower and more pronounced radial flow peak
at low $p_{\rm T}$ as compared to the measured one for the 10% most central
Pb–Pb collisions. Most of the models underestimate the measured $R_{\rm AA}$
in both centrality classes for $p_{\rm T}<$ 1.5$\leavevmode\nobreak\
\mathrm{GeV}/c$, with the exception of Catania [122, 123] calculations, which
tend to overpredict the measured $R_{\rm AA}$ at low $p_{\rm T}$ in central
and semicentral collisions and PHSD predictions, which overestimate the
$R_{\rm AA}$ in the 30–50% centrality class. Nevertheless, it is important to
remark that the $R_{\rm AA}$ at low $p_{\rm T}$ is sensitive not only to the
modelling of the charm-quark interactions with the medium but also to the
parametrisation of the nPDFs used in the calculations and to the
hydrodynamical description of the underlying medium.

More stringent constraints to the implementation of the heavy-quark
interactions with the medium constituents can be provided by the simultaneous
comparison of $R_{\rm AA}$ and $v_{2}$ measurements in Pb–Pb collisions at
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$
[36] with models, as reported in Fig. 11. All the models describe reasonably
well the $v_{2}$ in the most central collisions, while they tend to
underestimate the measured points in the 2 $<p_{\rm T}<$
6$\leavevmode\nobreak\ \mathrm{GeV}/c$ interval for the 30–50% centrality
class, except for LBT which reproduces well the measured $v_{2}$ but misses
completely the $R_{\rm AA}$ in the same $p_{\rm T}$ range. In this $p_{\rm T}$
region, the measured $v_{2}$ originates predominantly from the charm-quark
interactions with the QGP constituents, which impart the flow of the medium to
the heavy quarks, and from the hadronisation via recombination, which enhances
the charm-hadron $v_{2}$ with respect to the one of the charm quark because
the D meson picks up the $v_{2}$ of the light quark. Similarly, the peak
observed in the $R_{\rm AA}$ for 2 $<p_{\rm T}<$ 6$\leavevmode\nobreak\
\mathrm{GeV}/c$ is also due to the interplay between the diffusion and
recombination of the charm quarks with the medium constituents. Thus, the
measurements of the $R_{\rm AA}$ and $v_{2}$ in this $p_{\rm T}$ region are
particularly sensitive to quark diffusion and thermalisation with the medium,
and to the hadronisation mechanisms.

The simultaneous description of $R_{\rm AA}$ and $v_{2}$ is challenging for
the models and therefore the data have the potential to constrain the model
ingredients and parameters. The global agreement between the measured $R_{\rm
AA}$ and the theoretical models was evaluated by computing a $\chi^{2}/$ndf,
as done in Refs. [126, 36]. The statistical and systematic uncertainties
(treating separately the contributions correlated and uncorrelated among
$p_{\rm T}$ intervals) were considered in the calculation together with the
theoretical ones, when available. Since the upper $p_{\rm T}$ limit of the
predictions is different for each model, the $\chi^{2}/$ndf was computed in
the 0 $<p_{\rm T}<$ 8$\leavevmode\nobreak\ \mathrm{GeV}/c$ interval, which is
common among all the models except for the LIDO predictions which start from
1.5$\leavevmode\nobreak\ \mathrm{GeV}/c$. Therefore, the $\chi^{2}/$ndf
computation for LIDO was performed excluding the first two $p_{\rm T}$
intervals of the $R_{\rm AA}$. This low-$p_{\rm T}$ range provides high
sensitivity to charm-quark diffusion and hadronisation in the QGP. The
$\chi^{2}/$ndf values are reported in Table 5. The large spread in the
computed $\chi^{2}/$ndf is not only due to the improved precision of the
measurement, but also to the differences among the theoretical models. They do
not only differ in terms of the interaction of charm quarks with the medium,
as previously highlighted, but also in terms of the considered nuclear PDFs,
the bulk evolution of the medium (i.e. ideal or viscous hydrodynamics), the
charm hadronisation mechanism (i.e. fragmentation and/or recombination), and
whether a hadronic phase is included or not. In particular, the charm-quark
recombination mechanism is implemented with different approaches in the
various models: most of the models use an instantaneous recombination based on
the Wigner function formalism [127], while TAMU implements a resonance-
recombination model [128], PHSD a dynamical coalescence via a Monte Carlo
approach [125], and POWLANG an in-medium string formation approach [117]. The
$R_{\rm AA}$ predictions are deeply influenced by these additional model
ingredients, other than the transport properties of the medium. Therefore, a
rather mild requirement on the data-to-model consistency, namely
$\chi^{2}/$ndf $<$ 5, was applied to select the models considered for the
estimation of the heavy-quark spatial diffusion coefficient. With this
criterion, the selected models are TAMU [76], MC@sHQ+EPOS2 [115], LIDO [124],
LGR [121], and Catania [122, 123]. A similar analysis was performed for the
elliptic and triangular flow [36]. In order to further constrain the
description of heavy-quark transport in the medium and the spatial diffusion
coefficient $D_{s}$, the results of the $\chi^{2}/$ndf analysis for the
$R_{\rm AA}$ model calculations were combined with those obtained for the
elliptic and triangular flow in Ref. [36]. Note that some theoretical
predictions for the $v_{2}$ were updated after the publication of Ref. [36],
namely LBT [119, 120], LIDO [124], LGR [121], PHSD [125], and TAMU [76].
Therefore the $\chi^{2}/$ndf of these predictions with respect to the measured
$v_{2}$ and $v_{3}$ was re-computed. The outcome did not change significantly
with respect to Ref. [36], with the only exception of LIDO [124] for which the
updated predictions provide a better description of the data with respect to
the old ones, resulting in a $\chi^{2}/$ndf $<$ 2\. Thus, the LIDO model is
included among those providing a good description of $v_{2}$ and $v_{3}$,
while it was not considered in Ref. [36]. The models that describe reasonably
well both $R_{\rm AA}$ (with $\chi^{2}/$ndf $<$ 5) and $v_{2}$ and $v_{3}$
(with $\chi^{2}/$ndf $<$ 2) are TAMU [76], MC@sHQ+EPOS2 [115] , LIDO [124],
LGR [121], and Catania [122, 123]. These models use a value of heavy-quark
spatial diffusion coefficient in the range $1.5<2\pi
D_{\mathrm{s}}T_{\mathrm{c}}<4.5$ at the pseudocritical temperature
$T_{\mathrm{pc}}=155\leavevmode\nobreak\ \mathrm{MeV}$ [3]. According to lQCD
calculations the $2\pi D_{\mathrm{s}}T_{\mathrm{c}}$ lies between 2 and 6
[129, 130, 131], and these results are also in agreement with the $D_{s}$
range estimated from the $v_{2}$ measurements by ALICE (1.5–7 from Ref. [36]),
STAR (2–12 from Ref. [132]), and PHENIX (value of $\sim$3 from Ref. [133]).
The inclusion of the $R_{\rm AA}$ in the $\chi^{2}/$ndf improved the
constraint on the spatial diffusion coefficient with respect to the range
reported in Ref. [36]. It is however important to remark that this coefficient
is not the only key parameter of the models describing the heavy-quark
transport in an expanding medium, but there are other ingredients (such as the
parameters of the underlying hydrodynamics, the modelling of the
hadronisation, the description of the interactions in the hadronic phase,
amongst others, see e.g. Refs. [17, 134]) playing an important role.

Table 5: Summary of the $\chi^{2}/\mathrm{ndf}$ values obtained in $0<p_{\rm T}<8\leavevmode\nobreak\ \mathrm{GeV}/c$ for the different model predictions compared with the measured D-meson $R_{\rm AA}$. Model | $\chi^{2}/\mathrm{ndf}$
---|---
$R_{\rm AA}$ | global
0–10% | 30–50%
Catania [122, 123] | 94.9$/$15 | 48.9$/$15 | 143.8$/$30
DAB-MOD [116] | 110.2$/$15 | 123.9$/$15 | 234.1$/$30
LBT [119, 120] | 342.6$/$15 | 69.2$/$15 | 411.8$/$30
LIDO [124] | 31.8$/$13 | 14.5$/$13 | 46.4$/$26
LGR [121] | 4.7$/$15 | 4.5$/$15 | 9.2$/$30
MC@sHQ+EPOS2 [115] | 31.8$/$15 | 24.8$/$15 | 56.6$/$30
PHSD [125] | 103.2$/$15 | 191.4$/$15 | 294.7$/$30
POWLANG-HTL [117, 118] | 331.0$/$15 | 137.6$/$15 | 468.6$/$30
TAMU [76] | 16.7$/$15 | 13.5$/$15 | 30.2$/$30

A deeper insight on the impact of the different implementations of the charm-
quark interaction and hadronisation in the QGP can be obtained from Figs. 12
and 13.

Figure 12: Prompt D-meson $R_{\rm AA}$ in the 0–10% centrality class (left
panel) and $v_{2}$ in the 30–50% centrality class (right panel) compared with
the LIDO [124] and LGR [121] predictions obtained with and without including
radiative processes in the charm-quark interactions with the medium.

In Fig. 12 the $R_{\rm AA}$ and the $v_{2}$, measured in the 0–10% and 30–50%
centrality classes, respectively, are compared with two different
implementations of the LIDO [124] and LGR [121] models, in order to assess the
role of elastic and radiative processes in the charm-quark interactions with
medium constituents. In particular, the first implementation is the standard
one including both elastic and radiative processes, while the second
prediction was obtained by switching off the radiative processes. At low
$p_{\rm T}$ (i.e. $\lesssim$ 3–4$\leavevmode\nobreak\ \mathrm{GeV}/c$), the
collisional processes are expected to be the dominant interaction and this is
confirmed by the similarity of the predictions for prompt D-meson $R_{\rm AA}$
and $v_{2}$ with (full band) and without (striped band) radiative processes.
This is further supported by the agreement, in the same $p_{\rm T}$ range,
between the experimental data and other theoretical models which implement
only elastic processes, like TAMU [76], PHSD [125], and POWLANG-HTL [117,
118]. Radiative processes, instead, become dominant at intermediate and high
$p_{\rm T}$. This can be observed in Fig. 12 where the predictions without
these processes overestimate (underestimate) the measured $R_{\rm AA}$
($v_{2}$) for $p_{\rm T}\gtrsim$ 5–6$\leavevmode\nobreak\ \mathrm{GeV}/c$.

Figure 13: Prompt D-meson $R_{\rm AA}$ in the 0–10% centrality class (left
panel) and $v_{2}$ in the 30–50% centrality class (right panel) compared with
the PHSD [125], POWLANG [117, 118], and DAB-MOD [116] predictions obtained
with and without including hadronisation via recombination.

Similarly, Fig. 13 shows the comparison of the experimental data for $R_{\rm
AA}$ and $v_{2}$ with two different versions of the PHSD [125], POWLANG-HTL
[117, 118], and DAB-MOD [116] models, in order to investigate the effects of
the hadronisation mechanism, and in particular the role of recombination. This
plays a key role in the predictions for $R_{\rm AA}$ and $v_{2}$, since the
relation between the momentum of the charm hadron and that of the parent charm
quark is different for the fragmentation process, where the hadron inherits a
fraction of the initial quark momentum, and the recombination, where the
D-meson $p_{\rm T}$ is larger than the one of the charm quark and the charm
hadron inherits also the collective flow of the light quark. In Fig. 13, the
predictions are provided with (solid line) and without (dashed line) the
implementation of the recombination process in the hadronisation mechanism.
The calculations performed including only the fragmentation process
underestimate both the $R_{\rm AA}$ and $v_{2}$, while the inclusion of the
recombination of charm quarks with light quarks pushes the predictions closer
to the experimental data. This indicates that recombination with light quarks
from the medium plays a relevant role in the hadronisation of charm quarks at
the QGP phase boundary.

##  6 Conclusions

We have reported the measurement of the $p_{\rm T}$-differential production
yields of prompt ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm D^{*+}}$ mesons and
charge conjugates at midrapidity ($|y|<$ 0.5) in central and semicentral Pb–Pb
collisions at a centre-of-mass energy per nucleon pair
$\sqrt{s_{\rm\scriptscriptstyle NN}}=5.02\leavevmode\nobreak\ \mathrm{TeV}$.
The results were obtained with the data sample collected at the end of 2018
with the ALICE detector. The $p_{\rm T}$-spectra were measured in finer
$p_{\rm T}$ intervals with respect to the previous measurements at the same
centre-of-mass energy [30], providing a more precise description of the
$p_{\rm T}$ distribution. The large data sample allowed for the first
measurement of the ${\rm D^{0}}$-meson yield down to $p_{\rm T}=$ 0 in Pb–Pb
collisions at the LHC. This enabled the determination of the $p_{\rm
T}$-integrated production yields of prompt ${\rm D^{0}}$, which is obtained in
a model independent way, of prompt ${\rm D^{+}}$ and ${\rm D^{*+}}$, and the
comparison with predictions from the statistical hadronisation model [66, 67].
The average $R_{\rm AA}$ of the ${\rm D^{0}}$, ${\rm D^{+}}$, and ${\rm
D^{*+}}$ mesons reaches a maximum suppression value of 5.5 (i.e. $R_{\rm
AA}\sim$ 0.18) and 2.7 (i.e. $R_{\rm AA}\sim$ 0.4) in the 0–10% and 30–50%
centrality classes, respectively, at $p_{\rm T}=$ 6–8$\leavevmode\nobreak\
\mathrm{GeV}/c$. This suppression becomes less pronounced in peripheral
collisions with a minimum value of 0.7–0.8, as observed in the 60–80%
centrality class in Ref. [30], and it is due to final-state effects, since it
is not observed in minimum bias p–Pb collisions [33]. However, it was pointed
out in Ref. [135, 136] that event selection and geometry biases can cause an
apparent suppression of $R_{\rm AA}$, especially in peripheral collisions,
even in the absence of nuclear effect. The $R_{\rm AA}$ of prompt D mesons at
$p_{\rm T}>$ 8$\leavevmode\nobreak\ \mathrm{GeV}/c$ is well described by
models which include both collisional and radiative energy loss processes. In
addition, the $p_{\rm T}$-integrated $R_{\rm AA}$ of prompt ${\rm D^{0}}$
mesons was obtained and compared with the $R_{\rm pPb}$ measured in p–Pb
collisions at the same centre-of-mass energy [33]. The results lie on the
upper edge of the calculations which consider nuclear modification of the
PDFs.

In central collisions, the $p_{\rm T}$-differential $R_{\rm AA}$ of prompt D
mesons is compatible within uncertainties with that of charged particles for
$p_{\rm T}\gtrsim$ 8 $\mathrm{GeV}/c$ and prompt ${\rm J}/\psi$ mesons for
$p_{\rm T}\gtrsim$ 10 $\mathrm{GeV}/c$, while an ordering in $R_{\rm AA}$ is
observed at lower $p_{\rm T}$. The latter is due to the interplay of different
effects, such as the diffusion of charm quarks in the medium possibly leading
to their thermalisation in the QGP and the hadronisation via recombination
with light quarks from the medium, which along with energy loss contribute to
the modification of the $p_{\rm T}$-distribution of the hadrons. The
comparison of the $R_{\rm AA}$ of prompt D mesons and non-prompt ${\rm
J}/\psi$ for $p_{\rm T}>$ 6$\leavevmode\nobreak\ \mathrm{GeV}/c$ together with
model calculations supports the predictions of quark-mass dependent energy
loss. At low and intermediate $p_{\rm T}$ ($\lesssim$ 6–8$\leavevmode\nobreak\
\mathrm{GeV}/c$), the $R_{\rm AA}$ is well described by different transport
model calculations for the two centrality classes. However, most of them fail
in describing simultaneously both $R_{\rm AA}$ and $v_{2}$ in central and
semicentral Pb–Pb collisions. By considering the few models that are in fair
agreement with both observables, the heavy-quark spatial diffusion coefficient
was estimated to be in the range $1.5<2\pi D_{\mathrm{s}}T_{\mathrm{c}}<4.5$
at the pseudocritical temperature $T_{\mathrm{pc}}=155\leavevmode\nobreak\
\mathrm{MeV}$, which is a narrower interval as compared to estimations based
on previous D-meson measurements at LHC energies [126, 36]. Therefore, the
simultaneous comparison of the data with the theoretical predictions for these
observables allowed for more stringent constraints to the heavy-quark spatial
diffusion coefficient and for significant progress in the understanding of the
interaction processes and the hadronisation mechanisms of charm quarks in the
high-density QCD medium.

## Acknowledgements

The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector: A. I. Alikhanyan National Science Laboratory (Yerevan Physics Institute) Foundation (ANSL), State Committee of Science and World Federation of Scientists (WFS), Armenia; Austrian Academy of Sciences, Austrian Science Fund (FWF): [M 2467-N36] and Nationalstiftung für Forschung, Technologie und Entwicklung, Austria; Ministry of Communications and High Technologies, National Nuclear Research Center, Azerbaijan; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Financiadora de Estudos e Projetos (Finep), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Universidade Federal do Rio Grande do Sul (UFRGS), Brazil; Ministry of Education of China (MOEC) , Ministry of Science & Technology of China (MSTC) and National Natural Science Foundation of China (NSFC), China; Ministry of Science and Education and Croatian Science Foundation, Croatia; Centro de Aplicaciones Tecnológicas y Desarrollo Nuclear (CEADEN), Cubaenergía, Cuba; Ministry of Education, Youth and Sports of the Czech Republic, Czech Republic; The Danish Council for Independent Research | Natural Sciences, the VILLUM FONDEN and Danish National Research Foundation (DNRF), Denmark; Helsinki Institute of Physics (HIP), Finland; Commissariat à l’Energie Atomique (CEA) and Institut National de Physique Nucléaire et de Physique des Particules (IN2P3) and Centre National de la Recherche Scientifique (CNRS), France; Bundesministerium für Bildung und Forschung (BMBF) and GSI Helmholtzzentrum für Schwerionenforschung GmbH, Germany; General Secretariat for Research and Technology, Ministry of Education, Research and Religions, Greece; National Research, Development and Innovation Office, Hungary; Department of Atomic Energy Government of India (DAE), Department of Science and Technology, Government of India (DST), University Grants Commission, Government of India (UGC) and Council of Scientific and Industrial Research (CSIR), India; Indonesian Institute of Science, Indonesia; Istituto Nazionale di Fisica Nucleare (INFN), Italy; Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) KAKENHI and Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT)of Applied Science (IIST), Japan; Consejo Nacional de Ciencia (CONACYT) y Tecnología, through Fondo de Cooperación Internacional en Ciencia y Tecnología (FONCICYT) and Dirección General de Asuntos del Personal Academico (DGAPA), Mexico; Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), Netherlands; The Research Council of Norway, Norway; Commission on Science and Technology for Sustainable Development in the South (COMSATS), Pakistan; Pontificia Universidad Católica del Perú, Peru; Ministry of Education and Science, National Science Centre and WUT ID-UB, Poland; Korea Institute of Science and Technology Information and National Research Foundation of Korea (NRF), Republic of Korea; Ministry of Education and Scientific Research, Institute of Atomic Physics and Ministry of Research and Innovation and Institute of Atomic Physics, Romania; Joint Institute for Nuclear Research (JINR), Ministry of Education and Science of the Russian Federation, National Research Centre Kurchatov Institute, Russian Science Foundation and Russian Foundation for Basic Research, Russia; Ministry of Education, Science, Research and Sport of the Slovak Republic, Slovakia; National Research Foundation of South Africa, South Africa; Swedish Research Council (VR) and Knut & Alice Wallenberg Foundation (KAW), Sweden; European Organization for Nuclear Research, Switzerland; Suranaree University of Technology (SUT), National Science and Technology Development Agency (NSDTA) and Office of the Higher Education Commission under NRU project of Thailand, Thailand; Turkish Energy, Nuclear and Mineral Research Agency (TENMAK), Turkey; National Academy of Sciences of Ukraine, Ukraine; Science and Technology Facilities Council (STFC), United Kingdom; National Science Foundation of the United States of America (NSF) and United States Department of Energy, Office of Nuclear Physics (DOE NP), United States of America.

## References

  * [1] J. C. Collins and M. J. Perry, “Superdense matter: Neutrons or asymptotically free quarks?”, Phys. Rev. Lett. 34 (1975) 1353–1356.
  * [2] F. Karsch, “Lattice simulations of the thermodynamics of strongly interacting elementary particles and the exploration of new phases of matter in relativistic heavy ion collisions”, J. Phys. Conf. Ser. 46 (2006) 122–131, arXiv:hep-lat/0608003 [hep-lat].
  * [3] Wuppertal-Budapest Collaboration, S. Borsanyi et al., “Is there still any $T_{c}$ mystery in lattice QCD? Results with physical masses in the continuum limit III”, JHEP 1009 (2010) 073, arXiv:1005.3508 [hep-lat].
  * [4] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, “Full result for the QCD equation of state with 2+1 flavors”, Phys. Lett. B730 (2014) 99–104, arXiv:1309.5258 [hep-lat].
  * [5] A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. Ding, et al., “The chiral and deconfinement aspects of the QCD transition”, Phys. Rev. D85 (2012) 054503, arXiv:1111.1710 [hep-lat].
  * [6] W. Busza, K. Rajagopal, and W. van der Schee, “Heavy Ion Collisions: The Big Picture, and the Big Questions”, Ann. Rev. Nucl. Part. Sci. 68 (2018) 339–376, arXiv:1802.04801 [hep-ph].
  * [7] ALICE Collaboration, K. Aamodt et al., “Two-pion Bose-Einstein correlations in central Pb-Pb collisions at $\sqrt{{s}_{NN}}=$ 2.76 TeV”, Phys. Lett. B696 (2011) 328–337, arXiv:1012.4035 [nucl-ex].
  * [8] A. Andronic et al., “Heavy-flavour and quarkonium production in the LHC era: from proton-proton to heavy-ion collisions”, Eur. Phys. J. C76 no. 3, (2016) 107, arXiv:1506.03981 [nucl-ex].
  * [9] F.-M. Liu and S.-X. Liu, “Quark-gluon plasma formation time and direct photons from heavy ion collisions”, Phys. Rev. C89 no. 3, (2014) 034906, arXiv:1212.6587 [nucl-th].
  * [10] M. H. Thoma and M. Gyulassy, “Quark damping and energy loss in the high temperature QCD”, Nucl. Phys B351 (1991) 491–506.
  * [11] E. Braaten and M. H. Thoma, “Energy loss of a heavy fermion in a hot plasma”, Phys. Rev. D44 (1991) 1298–1310.
  * [12] E. Braaten and M. H. Thoma, “Energy loss of a heavy quark in the quark - gluon plasma”, Phys. Rev. D 44 no. 9, (1991) R2625.
  * [13] R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff, “Radiative energy loss and p(T) broadening of high-energy partons in nuclei”, Nucl. Phys. B484 (1997) 265–282, arXiv:hep-ph/9608322 [hep-ph].
  * [14] M. Gyulassy and M. Plumer, “Jet quenching in dense matter”, Phys. Lett. B243 (1990) 432–438.
  * [15] S. Batsouli, S. Kelly, M. Gyulassy, and J. L. Nagle, “Does the charm flow at RHIC?”, Phys. Lett. B557 (2003) 26–32, arXiv:nucl-th/0212068 [nucl-th].
  * [16] V. Greco, C. M. Ko, and R. Rapp, “Quark coalescence for charmed mesons in ultrarelativistic heavy ion collisions”, Phys. Lett. B595 (2004) 202–208, arXiv:nucl-th/0312100 [nucl-th].
  * [17] A. Beraudo et al., “Extraction of heavy-flavor transport coefficients in QCD matter”, Nucl. Phys. A979 (2018) 21–86, arXiv:1803.03824 [nucl-th].
  * [18] F. Prino and R. Rapp, “Open heavy flavor in QCD matter and in nuclear collisions”, J. Phys. G43 no. 9, (2016) 093002, arXiv:1603.00529 [nucl-ex].
  * [19] M. Djordjevic, “Heavy flavor puzzle at LHC: a serendipitous interplay of jet suppression and fragmentation”, Phys. Rev. Lett. 112 no. 4, (2014) 042302, arXiv:1307.4702 [nucl-th].
  * [20] Y. L. Dokshitzer and D. Kharzeev, “Heavy quark colorimetry of QCD matter”, Phys. Lett. B519 (2001) 199–206, arXiv:hep-ph/0106202 [hep-ph].
  * [21] N. Armesto, C. A. Salgado, and U. A. Wiedemann, “Medium induced gluon radiation off massive quarks fills the dead cone”, Phys. Rev. D69 (2004) 114003, arXiv:hep-ph/0312106 [hep-ph].
  * [22] M. Djordjevic and M. Gyulassy, “Heavy quark radiative energy loss in QCD matter”, Nucl. Phys A733 (2004) 265–298, arXiv:nucl-th/0310076 [nucl-th].
  * [23] B.-W. Zhang, E. Wang, and X.-N. Wang, “Heavy quark energy loss in nuclear medium”, Phys. Rev. Lett. 93 (2004) 072301, arXiv:nucl-th/0309040 [nucl-th].
  * [24] H. van Hees, V. Greco, and R. Rapp, “Heavy-quark probes of the quark-gluon plasma at RHIC”, Phys. Rev. C73 (2006) 034913, arXiv:nucl-th/0508055 [nucl-th].
  * [25] A. Andronic, P. Braun-Munzinger, K. Redlich, and J. Stachel, “Statistical hadronization of charm in heavy-ion collisions at SPS, RHIC and LHC”, Phys. Lett. B571 (2003) 36–44, arXiv:nucl-th/0303036 [nucl-th].
  * [26] R. Glauber and G. Matthiae, “High-energy scattering of protons by nuclei”, Nucl. Phys B21 (1970) 135–157.
  * [27] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, “Glauber modeling in high energy nuclear collisions”, Ann.Rev.Nucl.Part.Sci. 57 (2007) 205–243, arXiv:nucl-ex/0701025 [nucl-ex].
  * [28] C. Loizides, J. Nagle, and P. Steinberg, “Improved version of the PHOBOS Glauber Monte Carlo”, SoftwareX 1-2 (2015) 13–18, arXiv:1408.2549 [nucl-ex].
  * [29] ALICE Collaboration, J. Adam et al., “Transverse momentum dependence of D-meson production in Pb-Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}=$ 2.76 TeV”, JHEP 03 (2016) 081, arXiv:1509.06888 [nucl-ex].
  * [30] ALICE Collaboration, S. Acharya et al., “Measurement of D0, D+, D∗+ and D${}_{s}^{+}$ production in Pb-Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}=5.02$ TeV”, JHEP 10 (2018) 174, arXiv:1804.09083 [nucl-ex].
  * [31] CMS Collaboration, A. M. Sirunyan et al., “Nuclear modification factor of D0 mesons in PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV”, Phys. Lett. B782 (2018) 474–496, arXiv:1708.04962 [nucl-ex].
  * [32] STAR Collaboration, J. Adam et al., “Centrality and transverse momentum dependence of $D^{0}$-meson production at mid-rapidity in Au+Au collisions at ${\sqrt{s_{\rm NN}}=\rm{200}}$ GeV”, Phys. Rev. C99 no. 3, (2019) 034908, arXiv:1812.10224 [nucl-ex].
  * [33] ALICE Collaboration, S. Acharya et al., “Measurement of prompt D0, D+, D∗+, and ${\mathrm{D}}_{\mathrm{S}}^{+}$ production in p–Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}$ = 5.02 TeV”, JHEP 12 (2019) 092, arXiv:1906.03425 [nucl-ex].
  * [34] S. Voloshin and Y. Zhang, “Flow study in relativistic nuclear collisions by Fourier expansion of Azimuthal particle distributions”, Z. Phys. C 70 (1996) 665–672, arXiv:hep-ph/9407282.
  * [35] A. M. Poskanzer and S. A. Voloshin, “Methods for analyzing anisotropic flow in relativistic nuclear collisions”, Phys. Rev. C58 (1998) 1671–1678, arXiv:nucl-ex/9805001.
  * [36] ALICE Collaboration, S. Acharya et al., “Transverse-momentum and event-shape dependence of D-meson flow harmonics in Pb-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV”, Phys. Lett. B813 (2021) 136054, arXiv:2005.11131 [nucl-ex].
  * [37] ALICE Collaboration, S. Acharya et al., “Measurement of ${{\mathrm{D}}^{0}}$ , ${{\mathrm{D}}^{+}}$ , ${{\mathrm{D}}^{*+}}$ and ${{\mathrm{D}}^{+}_{\mathrm{s}}}$ production in pp collisions at ${\sqrt{{\textit{s}}}\leavevmode\nobreak\ =\leavevmode\nobreak\ 5.02\leavevmode\nobreak\ {\text{TeV}}}$ with ALICE”, Eur. Phys. J. C 79 no. 5, (2019) 388, arXiv:1901.07979 [nucl-ex].
  * [38] G. D. Moore and D. Teaney, “How much do heavy quarks thermalize in a heavy ion collision?”, Phys. Rev. C71 (2005) 064904, arXiv:hep-ph/0412346.
  * [39] ALICE Collaboration, K. Aamodt et al., “The ALICE experiment at the CERN LHC”, JINST 3 (2008) S08002.
  * [40] ALICE Collaboration, B. Abelev et al., “Performance of the ALICE Experiment at the CERN LHC”, Int.J.Mod.Phys. A29 (2014) 1430044, arXiv:1402.4476 [nucl-ex].
  * [41] ALICE Collaboration, K. Aamodt et al., “Alignment of the ALICE Inner Tracking System with cosmic-ray tracks”, JINST 5 (2010) P03003, arXiv:1001.0502 [physics.ins-det].
  * [42] J. Alme et al., “The ALICE TPC, a large 3-dimensional tracking device with fast readout for ultra-high multiplicity events”, Nucl. Instrum. Methods A 622 no. 1, (2010) 316 – 367.
  * [43] A. Akindinov et al., “Performance of the ALICE Time-Of-Flight detector at the LHC”, Eur. Phys. J. Plus 128 (2013) 44.
  * [44] ALICE Collaboration, E. Abbas et al., “Performance of the ALICE VZERO system”, JINST 8 (2013) P10016, arXiv:1306.3130 [nucl-ex].
  * [45] R. Arnaldi et al., “The zero degree calorimeters for the ALICE experiment”, Nucl. Instrum. Meth. A 581 (2007) 397–401. [Erratum: Nucl.Instrum.Meth.A 604, 765 (2009)].
  * [46] ALICE Collaboration, J. Adam et al., “Centrality dependence of the charged-particle multiplicity density at midrapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV”, Phys. Rev. Lett. 116 no. 22, (2016) 222302, arXiv:1512.06104 [nucl-ex].
  * [47] ALICE Collaboration, “Centrality determination in heavy ion collisions”, ALICE-PUBLIC-2018-011 (2018) .
  * [48] C. Loizides, J. Kamin, and D. d’Enterria, “Improved Monte Carlo Glauber predictions at present and future nuclear colliders”, Phys. Rev. C97 no. 5, (2018) 054910, arXiv:1710.07098 [nucl-ex]. [Erratum: Phys. Rev. C 99, 019901 (2019)].
  * [49] Particle Data Group Collaboration, P. Zyla et al., “Review of Particle Physics”, PTEP 2020 no. 8, (2020) 083C01.
  * [50] ALICE Collaboration, J. Adam et al., “Measurement of D${}_{s}^{+}$ production and nuclear modification factor in Pb-Pb collisions at $\sqrt{{\mathrm{s}}_{\mathrm{NN}}}=$ 2.76 TeV”, JHEP 03 (2016) 082, arXiv:1509.07287 [nucl-ex].
  * [51] ALICE Collaboration, S. Acharya et al., “Measurement of D-meson production at mid-rapidity in pp collisions at ${\sqrt{s}=7}$ TeV”, Eur. Phys. J. C 77 no. 8, (2017) 550, arXiv:1702.00766 [hep-ex].
  * [52] ALICE Collaboration, J. Adam et al., “D-meson production in p-Pb collisions at $\sqrt{s_{\rm NN}}=$5.02 TeV and in pp collisions at $\sqrt{s}=$7 TeV”, Phys. Rev. C94 no. 5, (2016) 054908, arXiv:1605.07569 [nucl-ex].
  * [53] M. Cacciari, M. Greco, and P. Nason, “The $p_{\rm T}$ spectrum in heavy flavor hadroproduction”, JHEP 9805 (1998) 007, arXiv:hep-ph/9803400 [hep-ph].
  * [54] M. Cacciari, S. Frixione, and P. Nason, “The $p_{\rm T}$ spectrum in heavy flavor photoproduction”, JHEP 0103 (2001) 006, arXiv:hep-ph/0102134 [hep-ph].
  * [55] T. Sjöstrand et al., “An introduction to PYTHIA 8.2”, Comput. Phys. Commun. 191 (2015) 159–177, arXiv:1410.3012 [hep-ph].
  * [56] R. Brun, F. Bruyant, F. Carminati, S. Giani, M. Maire, A. McPherson, G. Patrick, and L. Urban, GEANT: Detector Description and Simulation Tool; Oct 1994.  CERN Program Library. CERN, Geneva, 1993.  Long Writeup W5013.
  * [57] X. N. Wang and M. Gyulassy, “HIJING: A Monte Carlo model for multiple jet production in pp, pA and AA collisions”, Phys. Rev. D44 (1991) 3501–3516.
  * [58] P. Skands, S. Carrazza, and J. Rojo, “Tuning PYTHIA 8.1: the Monash 2013 Tune”, Eur. Phys. J. C 74 no. 8, (2014) 3024, arXiv:1404.5630 [hep-ph].
  * [59] ALICE Collaboration, B. Abelev et al., “Suppression of high transverse momentum D mesons in central Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV”, JHEP 1209 (2012) 112, arXiv:1203.2160 [nucl-ex].
  * [60] M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason, et al., “Theoretical predictions for charm and bottom production at the LHC”, JHEP 1210 (2012) 137, arXiv:1205.6344 [hep-ph].
  * [61] CMS Collaboration, A. M. Sirunyan et al., “Measurement of prompt and nonprompt charmonium suppression in PbPb collisions at 5.02 TeV”, Eur. Phys. J. C 78 no. 6, (2018) 509, arXiv:1712.08959 [nucl-ex].
  * [62] ATLAS Collaboration, M. Aaboud et al., “Prompt and non-prompt $J/\psi$ and $\psi(2\mathrm{S})$ suppression at high transverse momentum in $5.02\leavevmode\nobreak\ \mathrm{TeV}$ Pb+Pb collisions with the ATLAS experiment”, Eur. Phys. J. C 78 no. 9, (2018) 762, arXiv:1805.04077 [nucl-ex].
  * [63] M. Djordjevic and M. Djordjevic, “Predictions of heavy-flavor suppression at 5.1 TeV Pb+Pb collisions at the CERN Large Hadron Collider”, Phys. Rev. C92 no. 2, (2015) 024918, arXiv:1505.04316 [nucl-th].
  * [64] M. He, R. J. Fries, and R. Rapp, “Heavy flavor at the Large Hadron Collider in a strong coupling approach”, Phys. Lett. B735 (2014) 445–450, arXiv:1401.3817 [nucl-th].
  * [65] ALICE Collaboration, S. Acharya et al., “Measurement of beauty and charm production in pp collisions at $\sqrt{s}=5.02$ TeV via non-prompt and prompt D mesons”, JHEP 05 (2, 2021) 220, arXiv:2102.13601 [nucl-ex].
  * [66] A. Andronic, P. Braun-Munzinger, M. K. Köhler, K. Redlich, and J. Stachel, “Transverse momentum distributions of charmonium states with the statistical hadronization model”, Phys. Lett. B797 (2019) 134836, arXiv:1901.09200 [nucl-th].
  * [67] A. Andronic, P. Braun-Munzinger, M. K. Köhler, A. Mazeliauskas, K. Redlich, J. Stachel, and V. Vislavicius, “The multiple-charm hierarchy in the statistical hadronization model”, JHEP 07 (2021) 035, arXiv:2104.12754 [hep-ph].
  * [68] A. Mazeliauskas, S. Floerchinger, E. Grossi, and D. Teaney, “Fast resonance decays in nuclear collisions”, Eur. Phys. J. C 79 no. 3, (2019) 284, arXiv:1809.11049 [nucl-th].
  * [69] A. Andronic, P. Braun-Munzinger, K. Redlich, and J. Stachel, “Decoding the phase structure of QCD via particle production at high energy”, Nature 561 no. 7723, (2018) 321–330, arXiv:1710.09425 [nucl-th].
  * [70] ALICE Collaboration, S. Acharya et al., “Charm-quark fragmentation fractions and production cross section at midrapidity in pp collisions at the LHC”, arXiv:2105.06335 [nucl-ex].
  * [71] I. Kuznetsova and J. Rafelski, “Heavy flavor hadrons in statistical hadronization of strangeness-rich QGP”, Eur. Phys. J. C 51 (2007) 113–133, arXiv:hep-ph/0607203.
  * [72] P. R. Sorensen and X. Dong, “Suppression of non-photonic electrons from enhancement of charm baryons in heavy ion collisions”, Phys. Rev. C74 (2006) 024902, arXiv:nucl-th/0512042.
  * [73] S. H. Lee, K. Ohnishi, S. Yasui, I.-K. Yoo, and C.-M. Ko, “$\Lambda_{\rm{c}}$ enhancement from strongly coupled quark-gluon plasma”, Phys. Rev. Lett. 100 (2008) 222301, arXiv:0709.3637 [nucl-th].
  * [74] Y. Oh, C. M. Ko, S. H. Lee, and S. Yasui, “Heavy baryon/meson ratios in relativistic heavy ion collisions”, Phys. Rev. C79 (2009) 044905, arXiv:0901.1382 [nucl-th].
  * [75] M. He, R. J. Fries, and R. Rapp, “$D_{s}$-meson as quantitative probe of diffusion and hadronization in nuclear collisions”, Phys. Rev. Lett. 110 no. 11, (2013) 112301, arXiv:1204.4442 [nucl-th].
  * [76] M. He and R. Rapp, “Hadronization and charm-hadron ratios in heavy-ion collisions”, Phys. Rev. Lett. 124 no. 4, (2020) 042301, arXiv:1905.09216 [nucl-th].
  * [77] J. P. Lansberg and H. S. Shao, “Towards an automated tool to evaluate the impact of the nuclear modification of the gluon density on quarkonium, D and B meson production in proton–nucleus collisions”, Eur. Phys. J. C 77 no. 1, (2017) 1, arXiv:1610.05382 [hep-ph].
  * [78] A. Kusina, J. P. Lansberg, I. Schienbein, and H. S. Shao, “Gluon shadowing in heavy-flavor production at the LHC”, Phys. Rev. Lett. 121 no. 5, (2018) 052004, arXiv:1712.07024 [hep-ph].
  * [79] A. Kusina, J. P. Lansberg, I. Schienbein, and H. S. Shao, “Reweighted nuclear PDFs using heavy-flavor production data at the LHC: nCTEQ15rwHF & EPPS16rwHF”, Phys. Rev. D104 no. 1, (2021) 014010, arXiv:2012.11462 [hep-ph].
  * [80] H. S. Shao, “HELAC-Onia: An automatic matrix element generator for heavy quarkonium physics”, Comput. Phys. Commun. 184 (2013) 2562–2570, arXiv:1212.5293 [hep-ph].
  * [81] H. S. Shao, “HELAC-Onia 2.0: an upgraded matrix-element and event generator for heavy quarkonium physics”, Comput. Phys. Commun. 198 (2016) 238–259, arXiv:1507.03435 [hep-ph].
  * [82] K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, “EPPS16: Nuclear parton distributions with LHC data”, Eur. Phys. J. C 77 no. 3, (2017) 163, arXiv:1612.05741 [hep-ph].
  * [83] I. Helenius and H. Paukkunen, “Revisiting the D-meson hadroproduction in general-mass variable flavour number scheme”, JHEP 05 (2018) 196, arXiv:1804.03557 [hep-ph].
  * [84] ALICE Collaboration, S. Acharya et al., “Production of charged pions, kaons, and (anti-)protons in Pb-Pb and inelastic pp collisions at $\sqrt{s_{NN}}$ = 5.02 TeV”, Phys. Rev. C101 no. 4, (2020) 044907, arXiv:1910.07678 [nucl-ex].
  * [85] ALICE Collaboration, S. Acharya et al., “Transverse momentum spectra and nuclear modification factors of charged particles in pp, p–Pb and Pb–Pb collisions at the LHC”, JHEP 11 (2018) 013, arXiv:1802.09145 [nucl-ex].
  * [86] ALICE Collaboration, B. Abelev et al., “Production of charged pions, kaons and protons at large transverse momenta in pp and Pb–Pb collisions at $\sqrt{s_{\rm NN}}$ =2.76 TeV”, Phys. Lett. B736 (2014) 196–207, arXiv:1401.1250 [nucl-ex].
  * [87] ALICE Collaboration, S. Acharya et al., “Centrality and transverse momentum dependence of inclusive J/$\psi$ production at midrapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV”, Phys. Lett. B805 (2020) 135434, arXiv:1910.14404 [nucl-ex].
  * [88] ALICE Collaboration, S. Acharya et al., “Anisotropic flow of identified particles in Pb-Pb collisions at ${\sqrt{s}}_{\mathrm{NN}}=5.02$ TeV”, JHEP 09 (2018) 006, arXiv:1805.04390 [nucl-ex].
  * [89] J. Xu, J. Liao, and M. Gyulassy, “Consistency of perfect fluidity and jet quenching in semi-quark-gluon monopole plasmas”, Chin. Phys. Lett. 32 no. 9, (2015) 092501, arXiv:1411.3673 [hep-ph].
  * [90] J. Xu, J. Liao, and M. Gyulassy, “Bridging soft-hard transport properties of quark-gluon plasmas with CUJET3.0”, JHEP 02 (2016) 169, arXiv:1508.00552 [hep-ph].
  * [91] S. Shi, J. Liao, and M. Gyulassy, “Global constraints from RHIC and LHC on transport properties of QCD fluids in CUJET/CIBJET framework”, Chin. Phys. C 43 no. 4, (2019) 044101, arXiv:1808.05461 [hep-ph].
  * [92] S. Stojku, B. Ilic, M. Djordjevic, and M. Djordjevic, “Extracting the temperature dependence in high-$p_{\perp}$ particle energy loss”, Phys. Rev. C103 no. 2, (2021) 024908, arXiv:2007.07851 [nucl-th].
  * [93] D. Zigic, I. Salom, J. Auvinen, M. Djordjevic, and M. Djordjevic, “DREENA-B framework: first predictions of $R_{AA}$ and $v_{2}$ within dynamical energy loss formalism in evolving QCD medium”, Phys. Lett. B791 (2019) 236–241, arXiv:1805.04786 [nucl-th].
  * [94] D. Zigic, I. Salom, J. Auvinen, M. Djordjevic, and M. Djordjevic, “DREENA-C framework: joint $R_{AA}$ and $v_{2}$ predictions and implications to QGP tomography”, J. Phys. G 46 no. 8, (2019) 085101, arXiv:1805.03494 [nucl-th].
  * [95] Z.-B. Kang, F. Ringer, and I. Vitev, “Effective field theory approach to open heavy flavor production in heavy-ion collisions”, JHEP 03 (2017) 146, arXiv:1610.02043 [hep-ph].
  * [96] M. Djordjevic, “Theoretical formalism of radiative jet energy loss in a finite size dynamical QCD medium”, Phys. Rev. C80 (2009) 064909, arXiv:0903.4591 [nucl-th].
  * [97] M. Djordjevic and U. W. Heinz, “Radiative energy loss in a finite dynamical QCD medium”, Phys. Rev. Lett. 101 (2008) 022302, arXiv:0802.1230 [nucl-th].
  * [98] M. Djordjevic, “Collisional energy loss in a finite size QCD matter”, Phys. Rev. C74 (2006) 064907, arXiv:nucl-th/0603066.
  * [99] M. Gyulassy, P. Levai, and I. Vitev, “Non-abelian energy loss at finite opacity”, Phys. Rev. Lett. 85 (2000) 5535–5538, arXiv:nucl-th/0005032.
  * [100] M. Gyulassy, P. Levai, and I. Vitev, “Reaction operator approach to non-abelian energy loss”, Nucl. Phys. B594 (2001) 371–419, arXiv:nucl-th/0006010.
  * [101] C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart, “An effective field theory for collinear and soft gluons: Heavy to light decays”, Phys. Rev. D63 (2001) 114020, arXiv:hep-ph/0011336.
  * [102] C. W. Bauer, D. Pirjol, and I. W. Stewart, “Soft collinear factorization in effective field theory”, Phys. Rev. D65 (2002) 054022, arXiv:hep-ph/0109045.
  * [103] A. Idilbi and A. Majumder, “Extending soft-collinear-effective-theory to describe hard jets in dense QCD media”, Phys. Rev. D80 (2009) 054022, arXiv:0808.1087 [hep-ph].
  * [104] G. Ovanesyan and I. Vitev, “An effective theory for jet propagation in dense QCD matter: jet broadening and medium-induced bremsstrahlung”, JHEP 06 (2011) 080, arXiv:1103.1074 [hep-ph].
  * [105] M. Djordjevic, B. Blagojevic, and L. Zivkovic, “Mass tomography at different momentum ranges in quark-gluon plasma”, Phys. Rev. C94 no. 4, (2016) 044908, arXiv:1601.07852 [nucl-th].
  * [106] LHCb Collaboration, R. Aaij et al., “Study of J/$\psi$ Production in Jets”, Phys. Rev. Lett. 118 no. 19, (2017) 192001, arXiv:1701.05116 [hep-ex].
  * [107] R. Bain, L. Dai, A. Leibovich, Y. Makris, and T. Mehen, “NRQCD confronts LHCb data on quarkonium production within jets”, Phys. Rev. Lett. 119 no. 3, (2017) 032002, arXiv:1702.05525 [hep-ph].
  * [108] M. Spousta, “On similarity of jet quenching and charmonia suppression”, Phys. Lett. B767 (2017) 10–15, arXiv:1606.00903 [hep-ph].
  * [109] F. Arleo, “Quenching of hadron spectra in heavy ion collisions at the LHC”, Phys. Rev. Lett. 119 no. 6, (2017) 062302, arXiv:1703.10852 [hep-ph].
  * [110] ALICE Collaboration, J. Adam et al., “Differential studies of inclusive J/$\psi$ and $\psi$(2S) production at forward rapidity in Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV”, JHEP 05 (2016) 179, arXiv:1506.08804 [nucl-ex].
  * [111] K. Zhou, N. Xu, Z. Xu, and P. Zhuang, “Medium effects on charmonium production at ultrarelativistic energies available at the CERN Large Hadron Collider”, Phys. Rev. C89 no. 5, (2014) 054911, arXiv:1401.5845 [nucl-th].
  * [112] X. Du and R. Rapp, “Sequential regeneration of charmonia in heavy-ion collisions”, Nucl. Phys. A943 (2015) 147–158, arXiv:1504.00670 [hep-ph].
  * [113] M. He, R. J. Fries, and R. Rapp, “Heavy-quark diffusion and hadronization in quark-gluon plasma”, Phys. Rev. C86 (2012) 014903, arXiv:1106.6006 [nucl-th].
  * [114] S. Cao, K. J. Sun, S. Q. Li, S. Y. F. Liu, W. J. Xing, G. Y. Qin, and C. M. Ko, “Charmed hadron chemistry in relativistic heavy-ion collisions”, Phys. Lett. B807 (2020) 135561, arXiv:1911.00456 [nucl-th].
  * [115] M. Nahrgang, J. Aichelin, P. B. Gossiaux, and K. Werner, “Influence of hadronic bound states above $T_{c}$ on heavy-quark observables in Pb + Pb collisions at at the CERN Large Hadron Collider”, Phys. Rev. C89 no. 1, (2014) 014905, arXiv:1305.6544 [hep-ph].
  * [116] R. Katz, C. A. G. Prado, J. Noronha-Hostler, J. Noronha, and A. A. P. Suaide, “Sensitivity study with a $D$ and $B$ mesons modular simulation code of heavy flavor $R_{\mathrm{AA}}$ and azimuthal anisotropies based on beam energy, initial conditions, hadronization, and suppression mechanisms”, Phys. Rev. C102 no. 2, (2020) 024906, arXiv:1906.10768 [nucl-th].
  * [117] A. Beraudo, A. De Pace, M. Monteno, M. Nardi, and F. Prino, “Heavy flavors in heavy-ion collisions: quenching, flow and correlations”, Eur. Phys. J. C 75 no. 3, (2015) 121, arXiv:1410.6082 [hep-ph].