
\section{Analysis techniques}
\label{sec:analysis}

\subsection{D-meson raw yields}
\label{sec:first_method}
 
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\pm0.031\%$), ${\rm D}^+\to {\rm K}^-\pi^+\pi^+$ ($\mathrm{BR} = 9.38\pm0.16\%$), and ${\rm D}^{*+} \to {\rm D}^0 \pi^+$ ($\mathrm{BR}=67.7\pm0.5\%$)~\cite{Zyla:2020zbs}. 
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~\cite{Zyla:2020zbs}). 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~\mathrm{GeV}/c$ in both centrality classes and of ${\rm D^+}$ (${\rm D^{*+}}$) candidates for $p_{\rm T}>$ 2.5 (3)$~\mathrm{GeV}/c$ and for $p_{\rm T}>$ 2 (2)$~\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 \mbox{$\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~\mathrm{GeV}/c$, and fixed to 0.8 for $p_{\rm T}>5~\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.~\cite{Adam:2015jda, Acharya:2018hre}. 
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.~\cite{Acharya:2017jgo}. 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 ($\pm2 \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 $\mbox{p--Pb}~$ collisions, in order to extend the measurement down to $p_{\rm T}= 0$~\cite{Adam:2016ich, Acharya:2019mgn}. This technique, which is used for the first time in the analysis of \mbox{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 \textit{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.~\ref{fig:Dmass} 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~\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~\ref{fig:Dmass} 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.
 
\begin{figure*}[!t]
\begin{center}
\includegraphics[width=1.\textwidth]{figures/performance_plot/Dmesons_InvMassFits_PbPb5TeV_010.eps}
\caption{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 \textit{S} are also reported. Top row: ${\rm D^0}$-meson candidates with $0<p_{\rm T}<1~\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~\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~\mathrm{GeV}/c$ (left) and ${\rm D^{*+}}$-meson candidates with $16<p_{\rm T}<24~\mathrm{GeV}/c$ (right).} 
\label{fig:Dmass} 
\end{center}
\end{figure*}
 
 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~\mathrm{GeV}/c$ for ${\rm D^0}$ and $p_{\rm T}<5~\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$~\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}$~\cite{Acharya:2018hre}.
 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 ~\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.

\subsection{Yield corrections and beauty feed-down subtraction}
\label{sec:corrected_yields}
 The D-meson raw yields were corrected to obtain the $p_{\rm T}$-differential yields of prompt D-mesons according to:
 \begin{equation}
  \frac{\mathrm{d^2}N}{\mathrm{d}y~\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}} } .
  \label{eq:dNdpt}
 \end{equation}

 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)~\cite{Cacciari:1998it,Cacciari:2001td} perturbative QCD calculations, after verifying that the D-meson rapidity distributions in FONLL are consistent with those from PYTHIA 8~\cite{Sjostrand:2014zea}.

The ($\mathrm{Acc}\times\epsilon$) correction was obtained separately for prompt and feed-down D mesons using simulations with the GEANT3 transport package~\cite{Brun:1994aa}, including a detailed description of the detector geometry and response as well as of the LHC beam conditions. The \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$ were produced with the HIJING v1.36~\cite{Wang:1991hta} 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~\cite{Skands:2014pea}. 
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$~\mathrm{GeV}/c$.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=0.4\textwidth]{./figures/performance_plot/DzeroLowPtAccEffVsPt_PbPb_5TeV_010.eps}
\includegraphics[width=0.4\textwidth]{./figures/performance_plot/DzeroAccEffVsPt_PbPb_5TeV_010.eps}
\includegraphics[width=0.4\textwidth]{./figures/performance_plot/DplusAccEffVsPt_PbPb_5TeV_010.eps}
\includegraphics[width=0.4\textwidth]{./figures/performance_plot/DstarAccEffVsPt_PbPb_5TeV_010.eps}
\caption{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 \mbox{Pb--Pb}~ collisions for the 0--10\% centrality class.  
} 
\label{fig:Deff}
\end{center}
\end{figure}
Figure~\ref{fig:Deff} 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.~\ref{fig:Deff} 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.~\cite{ALICE:2012ab, Acharya:2018hre} 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{equation}
  \label{eq:fcNbMethod}
\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\Deltap_{\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}
\end{equation}
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~\mathrm{TeV}$, estimated with FONLL calculations~\cite{Cacciari:2012ny}, 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.~\ref{eq:dNdpt}. 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~\mathrm{TeV}$~\cite{Acharya:2018hre} with that of $\rm J/\psi$ from B-meson decays at the same energy measured by the CMS~\cite{Sirunyan:2017isk} and the ATLAS~\cite{Aaboud:2018quy} 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~\mathrm{GeV}/c$~\cite{Acharya:2018hre} 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\timesR_{\rm AA}^{\rm prompt}$ was assumed to compute the $f_{\rm prompt}$ factor for D mesons with $3<p_{\rm T}<24~\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~\mathrm{GeV}/c$ and $24<p_{\rm T}<50~\mathrm{GeV}/c$, since model calculations predict a lower difference between $R_{\rm AA}^{\rm prompt}$ and $R_{\rm AA}^{\textnormal{feed-down}}$~\cite{Djordjevic:2015hra,He:2014cla}, the hypothesis $R_{\rm AA}^{\textnormal{feed-down}}=1.5\timesR_{\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~\ref{sec:systematics}.

\section{Conclusions}
\label{sec: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 \mbox{Pb--Pb}~ collisions at a centre-of-mass energy per nucleon pair $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\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~\cite{Acharya:2018hre}, 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 \mbox{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~\cite{Andronic:2019wva,Andronic:2021wva}.
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 \mbox{$p_{\rm T} =$ 6--8$~\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.~\cite{Acharya:2018hre}, and it is due to final-state effects, since it is not observed in minimum bias \mbox{p--Pb}~ collisions~\cite{Acharya:2019mno}. However, it was pointed out in Ref.~\cite{Loizides:2017sqq, ALICE:2018ekf} 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$~\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 \mbox{p--Pb}~ collisions at the same centre-of-mass energy~\cite{Acharya:2019mno}. 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$~\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$~\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 \mbox{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~\mathrm{MeV}$, which is a narrower interval as compared to estimations based on previous D-meson measurements at LHC energies~\cite{Acharya:2017qps,Acharya:2020pnh}.
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. 

\section{Detector and data sample}
\label{sec:detector}

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.~\cite{Aamodt:2008zz,Abelev:2014ffa}. 

The main detectors used for the track reconstruction and particle identification in the analysis are the Inner Tracking System (ITS)~\cite{Aamodt:2010aa}, the Time Projection Chamber (TPC)~\cite{ALME2010316}, and the Time-Of-Flight (TOF) detector~\cite{Akindinov:2013tea}. 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~\cite{Abbas:2013taa} and the Zero Degree Calorimeter (ZDC)~\cite{Arnaldi:1999zz}. 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 \mbox{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 $\pm10$~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.~\cite{Acharya:2018hre}. 
 The centrality classes used in the present analysis, together with the corresponding average nuclear overlap function $\av{T_{\rm AA}}$~\cite{Adam:2015ptt, ALICE-PUBLIC-2018-011}, the number of events ($N_{\rm events}$), and the recorded integrated luminosity $L_{\rm int}$~\cite{Loizides:2017ack} in each class, are summarised in Table~\ref{tab:Nevents}.
 

	\begin{table}[!t]
	\centering
	\caption{Average nuclear overlap function, number of events, and recorded integrated luminosities for the two centrality classes used in the analysis.}
	\begin{tabular}{cccc}
	\hline
        \noalign{\vskip 0.1cm}
	Centrality class & $\langleT_{\rm AA}\rangle$ (mb$^{-1}$) & $N_{\rm events}$ & $L_{\rm int} (\mu \mathrm{b}^{-1})$\\[0.1cm]
	\hline
        \noalign{\vskip 0.1cm}
	\phantom{0}0--10\% & $23.26\pm0.17$ & $100 \times 10^6$ & 130.5 $\pm$ 0.5 \\
	30--50\% & $3.92\pm0.06$ & $85 \times 10^6$ & 55.5 $\pm$ 0.2\\
	\hline
	\end{tabular}		
	\label{tab:Nevents}
	\end{table}

\section*{Acknowledgements}
\input{fa_2021-09-21.tex}   
\end{acknowledgement}

\bibliographystyle{utphys}  


\section{Introduction}
\label{sec: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~\cite{Collins:1975us, Karsch:2006xs, Borsanyi:2010bp, Borsanyi:2013bia, Bazavov:2011nk}. 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~\cite{Busza:2018rrf} and its lifetime at the energy densities reached at the LHC is of the order of 10 $\mathrm{fm}/c$~\cite{ALICE:2011dyt}.

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~\cite{Andronic:2015wma}) are shorter than the QGP formation time, which is about 0.3--1.5 $\mathrm{fm}/c$ at LHC energies~\cite{Liu:2012ax}. 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~\cite{Thoma:1990fm, Braaten:1991jj, Braaten:1991we} and inelastic~\cite{Baier:1996sk, Gyulassy:1990ye} 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~\cite{Batsouli:2002qf,Greco:2003vf}. 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~\cite{Rapp:2018qla}. 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~\cite{Prino:2016cni}, 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~\cite{Gyulassy:1990ye,Baier:1996sk}. 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~\cite{Djordjevic:2013pba}. 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~\cite{Dokshitzer:2001zm,Armesto:2003jh,Djordjevic:2003zk,Zhang:2003wk}. Similarly, collisional energy loss is predicted to depend on the quark mass and to be smaller for heavy quarks~\cite{vanHees:2005wb}.

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~\cite{Prino:2016cni}. 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~\cite{Greco:2003vf, Andronic:2003zv}.

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}\pt$) and the production cross section in proton--proton collisions (${\rm d}\sigma_{\rm pp}/{\rm d}\pt$) scaled by the average nuclear overlap function $\langleT_{\rm AA}\rangle$
\begin{equation}
\label{eq:Raa}
R_{\rm AA}(p_{\rm T})=
{1\over \av{T_{\rm AA}}} {{\rm d} N_{\rm AA}/{\rm d}p_{\rm T} \over
{\rm d}\sigma_{\rm pp}/{\rm d}p_{\rm T}},
\end{equation}

where $\langleT_{\rm AA}\rangle$ is defined as the average number of binary nucleon--nucleon collisions ($\langleN_{\rm coll}\rangle$), which can be estimated via Glauber model calculations~\cite{Glauber:1970jm,Miller:2007ri,Loizides:2014vua}, divided by the inelastic nucleon--nucleon cross section.

Measurements of prompt D-meson production were performed by the ALICE Collaboration in \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} =$ 2.76 TeV~\cite{Adam:2015sza} and 5.02~$\mathrm{TeV}$~\cite{Acharya:2018hre} 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/\textit{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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$~\cite{Sirunyan:2017xss} 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 \mbox{Au--Au}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} =$ 200 GeV~\cite{Adam:2018inb} 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 \mbox{p--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$~\cite{Acharya:2019mno}, 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~\cite{Voloshin:1994mz,Poskanzer:1998yz}. 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~\cite{Acharya:2020pnh} in \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\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$\%$) \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$ are reported.
The statistical and systematic uncertainties have been reduced due to the \mbox{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~\cite{Acharya:2018hre}, and to the production cross section measured in pp collisions at the same centre-of-mass energy~\cite{Acharya:2019mgn}.
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~\cite{Acharya:2018hre} 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~\cite{Acharya:2019mgn}.

The paper is structured as follows. The data sample and the experimental apparatus are briefly introduced in Section~\ref{sec:detector}. The D-meson reconstruction procedure and the corrections applied to the raw yields are presented in Section~\ref{sec:analysis}, while the estimation of the systematic uncertainties is described in Section~\ref{sec:systematics}. The results are presented in Section~\ref{sec:results}, together with a comparison with the charged-pion, proton, and prompt and non-prompt ${\rm J}/\psi~R_{\rm AA}$, along with model calculations.
In addition, a simultaneous comparison of the measured $R_{\rm AA}$ and elliptic flow $v_2$~\cite{Acharya:2020pnh} 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~\cite{Moore:2004tg}, as well as on the role of the different hadronisation mechanisms. Conclusions and perspectives are drawn in Section~\ref{sec:conclusions}.




\section{Results}
\label{sec:results}

\subsection{Transverse-momentum-differential yields}
In this section, the results on the $p_{\rm T}$-differential (${\rm d}^2N/{\rm d}y {\rm d}\pt$) 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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\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}^2N/{\rm d}y {\rm d}\pt$ were obtained from the analyses with and without selections on the decay-vertex topology. They are compared in Fig.~\ref{fig:Dzero_spectra}, 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 \mbox{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~\mathrm{GeV}/c$ and those from the analysis with decay-vertex reconstruction for $p_{\rm T}>1~\mathrm{GeV}/c$. 

\begin{figure*}[!t]
\begin{center}
\includegraphics[width=0.495\textwidth]{figures/pTspectra_raa_plot/PromptD0_pTspectra_woVertexing_010_PbPb5TeV.eps}
\includegraphics[width=0.495\textwidth]{figures/pTspectra_raa_plot/PromptD0_pTspectra_woVertexing_3050_PbPb5TeV.eps}
\caption{Transverse momentum distributions ${\rm d}^2N/{\rm d}y {\rm d}\pt$ 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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$. Statistical uncertainties (bars) and systematic uncertainties (boxes) are shown.} 
\label{fig:Dzero_spectra} 
\end{center}
\end{figure*}

The $p_{\rm T}$-differential yields ${\rm d}^2N/{\rm d}y {\rm d}\pt$ of prompt ${\rm D^0}$, ${\rm D^+}$, and ${\rm D^{*+}}$ mesons are shown in Fig.~\ref{fig:Dspectra} for the 0--10\% and 30--50\% centrality classes. 
They are compared with the reference yields from pp collisions, which are computed as $\langleT_{\rm AA}\rangle\times{\rm d}\sigma_{\rm pp}/{\rm d}\pt$, where ${\rm d}\sigma_{\rm pp}/{\rm d}\pt$ is the D-meson $p_{\rm T}$-differential cross section measured in pp collisions at $\sqrt{s} = 5.02~\mathrm{TeV}$~\cite{Acharya:2019mgn,Acharya:2021cqv}, and $\langleT_{\rm AA}\rangle$ is the average nuclear overlap function~\cite{ALICE-PUBLIC-2018-011}. 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.~\cite{Acharya:2018hre,Adam:2015sza}. 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~\mathrm{GeV}/c$ is visible in \mbox{Pb--Pb}~ collisions and it is stronger in central than in semicentral collisions. 

\begin{figure*}[!t]
\begin{center}
\includegraphics[width=1.\textwidth]{figures/pTspectra_raa_plot/DzeroDplusDstar_dNdpt_010_3050.eps}
\caption{Transverse momentum distributions ${\rm d}^2N/{\rm d}y {\rm d}\pt$ 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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$. The reference pp~distributions multiplied by $\langleT_{\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~\cite{Acharya:2019mgn,Acharya:2021cqv} and $p_{\rm T}$-extrapolated values, respectively, of the $p_{\rm T}$-spectrum.} 
\label{fig:Dspectra} 
\end{center}
\end{figure*}

\begin{figure}[!t]
\begin{center}
\includegraphics[width=1.05\textwidth]{./figures/ratios_plot/Dplus_over_D0_Dstar_over_D0_PbPb2018_SHM_010_3050_4pads.eps}
\caption{${\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 \mbox{Pb--Pb}~ collisions compared to the results obtained from pp collisions~\cite{Acharya:2019mgn}. The right panels show the ratios in central \mbox{Pb--Pb}~ collisions compared to the predictions from the statistical hadronisation model (SHMc)~\cite{Andronic:2019wva,Andronic:2021wva} with a core--corona approach and the FastReso package~\cite{Mazeliauskas:2018irt} for resonance decays. Statistical (bars) and systematic (boxes) uncertainties are shown.} 
\label{fig:Dmeson_ratios}
\end{center}
\end{figure}

Figure~\ref{fig:Dmeson_ratios} 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~\mathrm{TeV}$~\cite{Acharya:2019mgn}. 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 \mbox{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)~\cite{Andronic:2019wva,Andronic:2021wva} (right panels of Fig.~\ref{fig:Dmeson_ratios}), 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}^2N/{\rm d}y {\rm d}\pt$ over the $p_{\rm T}$ intervals of the measurements. The results are reported in Table~\ref{tab:dNdy_visible}.
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 \mbox{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.~\ref{fig:Dmeson_ratios}, 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~\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)$~\mathrm{GeV}/c$ is negligible.
The yields of ${\rm D^+}$ and ${\rm D^{*+}}$ mesons in the interval $0<p_{\rm T}<2~\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~\ref{tab:dNdy_010_integrated}. 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)$~\mathrm{GeV}/c$ is negligible.

\begin{table}[!t]
\renewcommand*{\arraystretch}{1.2}
\begin{center}
\caption{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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$.}
\begin{tabular}{c|c|c}
\toprule
 & Kinematic range & Visible production yield \\
\hline
\multicolumn{3}{c}{\rule{0pt}{2.6ex} 0--10\% centrality} \\
\hline
 \rule{0pt}{2.6ex} 
 ${\rm D^0}$            & 0 $< p_{\rm T} <$ 50 $\mathrm{GeV}/c$  & 6.819 $\pm$ 0.457 {\rm(stat.)} $^{+0.912}_{-0.936}$ {\rm(syst.)} $\pm$ 0.054 {\rm(BR)}\\
 \rule{0pt}{2.6ex} 
 ${\rm D^+}$            & 2 $< p_{\rm T} <$ 50 $\mathrm{GeV}/c$  & 0.992 $\pm$ 0.073 {\rm(stat.)} $^{+0.154}_{-0.155}$ {\rm(syst.)} $\pm$ 0.017 {\rm(BR)}\\
 \rule{0pt}{2.6ex}\rule[-1.2ex]{0pt}{0pt} 
 ${\rm D^{*+}}$            & 3 $< p_{\rm T} <$ 50 $\mathrm{GeV}/c$  & 0.438 $\pm$ 0.037 {\rm(stat.)} $^{+0.084}_{-0.085}$ {\rm(syst.)} $\pm$ 0.005 {\rm(BR)} \\
\hline
\multicolumn{3}{c}{\rule{0pt}{2.6ex} 30--50\% centrality} \\
\hline
 \rule{0pt}{2.6ex} 
 ${\rm D^0}$            & 0 $< p_{\rm T} <$ 50 $\mathrm{GeV}/c$  & 1.275 $\pm$ 0.099 {\rm(stat.)} $^{+0.167}_{-0.173}$ {\rm(syst.)} $\pm$ 0.010 {\rm(BR)}\\
 \rule{0pt}{2.6ex} 
 ${\rm D^+}$            & 2 $< p_{\rm T} <$ 50 $\mathrm{GeV}/c$  & 0.179 $\pm$ 0.008 {\rm(stat.)} $^{+0.024}_{-0.024}$ {\rm(syst.)} $\pm$ 0.003 {\rm(BR)}\\
 \rule{0pt}{2.6ex}\rule[-1.2ex]{0pt}{0pt} 
 ${\rm D^{*+}}$            & 2 $< p_{\rm T} <$ 36 $\mathrm{GeV}/c$  & 0.230 $\pm$ 0.023 {\rm(stat.)} $^{+0.038}_{-0.039}$ {\rm(syst.)} $\pm$ 0.002 {\rm(BR)}\\
\bottomrule
\end{tabular}
\label{tab:dNdy_visible}
\end{center}
\end{table}

\begin{table}[!tb]
\renewcommand*{\arraystretch}{1.2}
\begin{center}
\caption{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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$. The right column reports the ${\rm d}N/{\rm d}y$ predicted by the GSI--Heidelberg statistical hadronisation model~\cite{Andronic:2019wva,Andronic:2021wva}.}
\begin{tabular}{c|l|c}
\toprule
 & \multicolumn{1}{|c|}{Measured ${\rm d}N/{\rm d}y$} & SHMc ${\rm d}N/{\rm d}y$ \\
\hline
\multicolumn{3}{c}{\rule{0pt}{2.6ex} 0--10\% centrality} \\
\hline
 \rule{0pt}{2.6ex} 
 ${\rm D^0}$            & 6.819 $\pm$ 0.457 {\rm(stat.)} $^{+0.912}_{-0.936}$ {\rm(syst.)} $\pm$ 0.054 {\rm(BR)} & 6.42 $\pm$ 1.07 \\
 \rule{0pt}{2.6ex} 
 ${\rm D^+}$            & 3.041 $\pm$ 0.073 {\rm(stat.)} $^{+0.154}_{-0.155}$ {\rm(syst.)} $\pm$ 0.052 {\rm(BR)} $^{+0.352}_{-0.618}$ {\rm(extrap.)} & 2.84 $\pm$ 0.47\\
 \rule{0pt}{2.6ex}\rule[-1.2ex]{0pt}{0pt} 
 ${\rm D^{*+}}$            & 3.803 $\pm$ 0.037 {\rm(stat.)} $^{+0.084}_{-0.085}$ {\rm(syst.)} $\pm$ 0.041 {\rm(BR)} $^{+0.854}_{-1.175}$ {\rm(extrap.)} & 2.52 $\pm$ 0.42\\
 \hline
\multicolumn{3}{c}{\rule{0pt}{2.6ex} 30--50\% centrality} \\
\hline
 \rule{0pt}{2.6ex} 
 ${\rm D^0}$            & 1.275 $\pm$ 0.099 {\rm(stat.)} $^{+0.167}_{-0.173}$ {\rm(syst.)} $\pm$ 0.010 {\rm(BR)} & 1.06 $\pm$ 0.15 \\
 \rule{0pt}{2.6ex} 
 ${\rm D^+}$            & 0.552 $\pm$ 0.008 {\rm(stat.)} $^{+0.024}_{-0.024}$ {\rm(syst.)} $\pm$ 0.009 {\rm(BR)} $^{+0.068}_{-0.114}$ {\rm(extrap.)} & 0.471 $\pm$ 0.069\\
 \rule{0pt}{2.6ex}\rule[-1.2ex]{0pt}{0pt} 
 ${\rm D^{*+}}$            & 0.663 $\pm$ 0.023 {\rm(stat.)} $^{+0.038}_{-0.039}$ {\rm(syst.)} $\pm$ 0.007 {\rm(BR)} $^{+0.149}_{-0.165}$ {\rm(extrap.)} & 0.419 $^{+0.065}_{-0.061}$\\
\bottomrule
\end{tabular}
\label{tab:dNdy_010_integrated}
\end{center}
\end{table}

In the right column the ${\rm d}N/{\rm d}y$ values predicted by the SHMc~\cite{Andronic:2019wva,Andronic:2021wva} 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 \mbox{p--Pb}~\ collisions along with guidance from parameterisations of the nuclear modification of the PDFs~\cite{Andronic:2017pug}.
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.~\cite{ALICE:2021dhb} 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~\cite{Kuznetsova:2006bh,Andronic:2003zv,Sorensen:2005sm,Lee:2007wr,Oh:2009zj,He:2012df,He:2019vgs}. 


\subsection{Nuclear modification factor}
The $R_{\rm AA}$ of prompt ${\rm D^0}$, ${\rm D^+}$, and ${\rm D^{*+}}$ mesons was computed, as defined in Eq.~\ref{eq:Raa}, using the ${\rm d}^2N/{\rm d}y {\rm d}\pt$ measured in \mbox{Pb--Pb}~ collisions and the pp reference at the same centre of mass energy reported in Fig.~\ref{fig:Dspectra}. The obtained results are shown in Fig.~\ref{fig:Dmeson_raa} for the two centrality classes. The nuclear modification factors of the three D-meson species are compatible among each other within statistical uncertainties.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=0.49\textwidth]{./figures/raa_plots/DmesonComparison_010.eps}
\includegraphics[width=0.49\textwidth]{./figures/raa_plots/DmesonComparison_3050.eps}
\caption{$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~\cite{Acharya:2019mgn,Acharya:2021cqv}, empty markers with the $p_{\rm T}$-extrapolated reference.} 
\label{fig:Dmeson_raa}
\end{center}
\end{figure}
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.~\ref{fig:Dmeson_avg_raa} together with the $R_{\rm AA}$ in the 60--80\% centrality class taken from Ref.~\cite{Acharya:2018hre}, which was measured using the sample of \mbox{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 \mbox{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.~\ref{fig:Dmeson_avg_raa} is the nuclear modification factor $R_{\rm pPb}$ measured in $\mbox{p--Pb}~$ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$ taken from Ref.~\cite{Acharya:2019mno}, which is compatible with unity within uncertainties, confirming that the suppression observed in \mbox{Pb--Pb}~ collisions is to due final-state effects induced by the formation of a hot and dense QGP medium.
\begin{figure}[!t]
\begin{center}
\includegraphics[width=0.49\textwidth]{./figures/raa_plots/DmesonAverage_010_3050_6080_pPb_comparison_.eps}
\includegraphics[width=0.49\textwidth]{./figures/raa_plots/DmesonAverage_5dot02_2dot76_200_comparison_.eps}
\caption{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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\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)~\cite{Acharya:2018hre} and in \mbox{p--Pb}~ collisions at the same centre-of-mass energy~\cite{Acharya:2019mno}.  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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}}=~$5.02 TeV and 2.76 TeV~\cite{Adam:2015sza} compared to the ${\rm D^0}$ $R_{\rm AA}$ measured by the STAR collaboration in $\mbox{Au--Au}~$ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}}=~$200 GeV~\cite{Adam:2018inb}.} 
\label{fig:Dmeson_avg_raa}
\end{center}
\end{figure}

The right panel of Fig.~\ref{fig:Dmeson_avg_raa} shows the $R_{\rm AA}$ of prompt D mesons measured by the ALICE collaboration in \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$ and 2.76 TeV~\cite{Adam:2015sza} and of ${\rm D^0}$ mesons measured at RHIC by the STAR collaboration in \mbox{Au--Au}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}}=~$200 GeV~\cite{Adam:2018inb}. 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~\cite{Djordjevic:2015hra}. 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~\mathrm{GeV}/c$) and a larger $R_{\rm AA}$ at high $p_{\rm T}$ ($>4~\mathrm{GeV}/c$) in collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}}=~$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 \mbox{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~\mathrm{TeV}$~\cite{Acharya:2019mgn}. The results in the 0--10\% and 30--50\% centrality classes are:
\begin{equation}
\begin{split}
	R^{{\mathrm{prompt~D^0}}}_{{\mathrm{AA}}} (0\text{--}10\%) = 0.689 \pm 0.054~ ({\mathrm{stat.}}) ^{+0.104}_{-0.106}~({\mathrm{syst.}}), \\
	R^{{\mathrm{prompt~D^0}}}_{{\mathrm{AA}}} (30\text{--}50\%) = 0.775 \pm 0.069~({\mathrm{stat.}}) ^{+0.117}_{-0.120}~({\mathrm{syst.}}).	
\end{split}
\end{equation}
Figure~\ref{fig:D0_integrated_raa} shows the results obtained in \mbox{Pb--Pb}~ collisions compared with the nuclear modification factor $R_{\rm pPb}$ measured in \mbox{p--Pb}~~collisions at the same centre-of-mass energy~\cite{Acharya:2019mno}. The total charm cross section is expected to scale with the number of binary collisions $N_{\rm coll}$, as introduced in Section~\ref{sec:introduction}, and thus the $R_{\rm AA}$ should be equal to one. However, the nuclear shadowing effect reduces the charm production in \mbox{Pb--Pb}~ (and \mbox{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 \mbox{Pb--Pb}~ collisions compared to pp collisions~\cite{Kuznetsova:2006bh,Andronic:2003zv,Sorensen:2005sm,Lee:2007wr,Oh:2009zj,He:2012df,He:2019vgs}.
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 \mbox{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 \mbox{p--Pb}~~compared to \mbox{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~\cite{Lansberg:2016deg,Kusina:2017gkz,Kusina:2020dki,Shao:2012iz,Shao:2015vga} and EPPS16~\cite{Eskola:2016oht,Helenius:2018uul}. 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 \mbox{p--Pb}~ collisions at the LHC~\cite{Kusina:2017gkz}, and are labelled as nCTEQ15$_{\mathrm{rwHF}}$ in Fig.~\ref{fig:D0_integrated_raa}. 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.~\cite{Kusina:2017gkz}. 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.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=0.65\textwidth]{./figures/Dzero-RAA-PtIntegr-withEPPS16mbnCTEQ15rwHF90_2.eps}
\caption{$p_{\rm T}$-integrated nuclear modification factors of prompt ${\rm D^0}$ mesons measured in \mbox{p--Pb}~~\cite{Acharya:2019mno} and \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\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)~\cite{Lansberg:2016deg,Kusina:2017gkz,Kusina:2020dki,Shao:2012iz,Shao:2015vga} and EPPS16~\cite{Eskola:2016oht,Helenius:2018uul} that include only initial-state effects.} 
\label{fig:D0_integrated_raa}
\end{center}
\end{figure}









\subsection{Discussion: energy loss regime}
\label{sec:results_pQCD}


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.~\ref{fig:D_pi_jpsi} with that of charged pions~\cite{Acharya:2019yoi} and charged particles~\cite{Acharya:2018qsh}.
The charged-particle $R_{\rm AA}$ is shown for $p_{\rm T}>20~\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$~\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~\cite{Abelev:2014laa,Acharya:2019yoi}.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=0.7\textwidth]{./figures/raa_pQCD_models/D-pi-JpsiPub-Raa.eps}
\caption{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~\cite{Acharya:2019yoi}, charged particles~\cite{Acharya:2018qsh}, inclusive J/$\psi$ measured by ALICE~\cite{Acharya:2019lkh}, and of prompt and non-prompt J/$\psi$ from CMS~\cite{Sirunyan:2017isk} in \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$.}
\label{fig:D_pi_jpsi}
\end{center}
\end{figure}

The $R_{\rm AA}$ of D mesons is larger than that of pions for $p_{\rm T}<8~\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.~\cite{Djordjevic:2013pba}, 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$~\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$~\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 \mbox{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.~\cite{Acharya:2020pnh}) and pions (taken from Ref.~\cite{Acharya:2019yoi, Acharya:2018zuq}) in the 0--10\% and 30-50\% centrality classes are compared in Fig.~\ref{fig:D_pi_models} to three of these models, namely CUJET3.1~\cite{Xu:2014tda,Xu:2015bbz,Shi:2019nyp}, DREENA-A~\cite{Stojku:2020tuk,Zigic:2018ovr,Zigic:2018smz}, and ${\rm SCET_{M,G}}$~\cite{Kang:2016ofv}.
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~\cite{Djordjevic:2009cr,Djordjevic:2008iz,Djordjevic:2006tw} 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~\cite{Gyulassy:2000fs,Gyulassy:2000er, Djordjevic:2003zk} 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~\cite{Bauer:2000yr,Bauer:2001yt}) describing the parton shower formed in the vacuum via soft and collinear splittings.
This effective theory was extended to the case of massless~\cite{Idilbi:2008vm,Ovanesyan:2011xy} and massive~\cite{Kang:2016ofv} 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~\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.



\begin{figure}[!t]
\begin{center}
\includegraphics[width=\textwidth]{./figures/raa_pQCD_models/D-pi-Raa-v2-models_uncBand.eps}
\caption{Average $R_{\rm AA}$ (top panels) and $v_2$ (bottom panels) of prompt ${\rm D^0}$, ${\rm D^+}$, and ${\rm D^{*+}}$ mesons~\cite{Acharya:2020pnh} and pions (charged hadrons)~\cite{Acharya:2019yoi,Acharya:2018qsh,Acharya:2018zuq}  compared with predictions from models based on pQCD calculations, namely CUJET3.1~\cite{Xu:2014tda,Xu:2015bbz,Shi:2019nyp}, DREENA-A~\cite{Stojku:2020tuk}, and ${\rm SCET_{\rm M,G}}$~\cite{Kang:2016ofv}.}
\label{fig:D_pi_models}
\end{center}
\end{figure}

The comparison of the nuclear modification factors of particles originating from charm and beauty quarks is also reported in Fig.~\ref{fig:D_pi_jpsi}, 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~\cite{Acharya:2019lkh} at midrapidity in the 0--20\% centrality class and those of prompt and non-prompt ${\rm J}/\psi$ from the CMS collaboration~\cite{Sirunyan:2017isk} 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.~\cite{Djordjevic:2016vfo}, 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.~\ref{fig:D_pi_jpsi}, 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~\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~\cite{Aaij:2017fak,Bain:2017wvk}, which experiences in-medium energy loss before fragmenting~\cite{Spousta:2016agr,Arleo:2017ntr}.
At lower $p_{\rm T}$ ($2 \lesssim p_{\rm T} \lesssim 4~\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~\cite{Adam:2015isa,Acharya:2019lkh}, 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~\cite{Zhou:2014kka,Du:2015wha,He:2011qa,Cao:2019iqs}, with possible small differences arising from the different kinematics involved in the charm-quark coalescence with $\overline{\rm c}$ and light antiquarks.



\subsection{Discussion: transport models and an estimate of the spatial diffusion coefficient}
\label{sec:results_transport}
As discussed in Section~\ref{sec:introduction}, 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.~\ref{fig:Dmeson_transport} 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~\cite{Nahrgang:2013xaa,Katz:2019fkc,Beraudo:2014boa,Beraudo:2017gxw,Cao:2016gvr,Cao:2017hhk,He:2019vgs,Li:2019lex,Scardina:2017ipo,Plumari:2019hzp,Ke:2018jem,Song:2015sfa}. In particular, the TAMU~\cite{He:2019vgs}, POWLANG-HTL~\cite{Beraudo:2014boa,Beraudo:2017gxw}, PHSD~\cite{Song:2015sfa}, and Catania~\cite{Scardina:2017ipo,Plumari:2019hzp} models describe the interactions of the charm quarks with the medium constituents solely via collisional processes, while the MC@sHQ+EPOS2~\cite{Nahrgang:2013xaa}, DAB-MOD~\cite{Katz:2019fkc}, LBT~\cite{Cao:2016gvr,Cao:2017hhk}, LGR~\cite{Li:2019lex}, and LIDO~\cite{Ke:2018jem} 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.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=1.05\textwidth]{./figures/raa_transport_models/D_Raa_V2_vs_transportModel_010.eps}
\includegraphics[width=1.05\textwidth]{./figures/raa_transport_models/D_Raa_V2_vs_transportModel_3050.eps}

\caption{Average $R_{\rm AA}$ (left) and elliptic flow $v_{2}$~\cite{Acharya:2020pnh} (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~\cite{Nahrgang:2013xaa,Katz:2019fkc,Beraudo:2014boa,Beraudo:2017gxw,Cao:2016gvr,Cao:2017hhk,He:2019vgs,Li:2019lex,Scardina:2017ipo,Plumari:2019hzp,Ke:2018jem,Song:2015sfa}.} 
\label{fig:Dmeson_transport}
\end{center}
\end{figure}

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$~\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$~\mathrm{GeV}/c$). In particular, the LBT~\cite{Cao:2016gvr,Cao:2017hhk}, LIDO~\cite{Ke:2018jem}, MC@sHQ+EPOS2~\cite{Nahrgang:2013xaa}, and Catania~\cite{Scardina:2017ipo,Plumari:2019hzp} calculations tend to underestimate the results for 2 $< p_{\rm T} <$ 4$~\mathrm{GeV}/c$ in both centrality classes, while the PHSD~\cite{Song:2015sfa} prediction is above the measured $R_{\rm AA}$ in semicentral collisions for the same $p_{\rm T}$ range. The POWLANG-HTL~\cite{Beraudo:2014boa,Beraudo:2017gxw} 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 \mbox{Pb--Pb}~ collisions. Most of the models underestimate the measured $R_{\rm AA}$ in both centrality classes for $p_{\rm T} <$ 1.5$~\mathrm{GeV}/c$, with the exception of Catania~\cite{Scardina:2017ipo,Plumari:2019hzp} 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 \mbox{Pb--Pb}~ collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\mathrm{TeV}$~\cite{Acharya:2020pnh} with models, as reported in Fig.~\ref{fig:Dmeson_transport}. 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$~\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$~\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.~\cite{Acharya:2017qps, Acharya:2020pnh}. 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$~\mathrm{GeV}/c$ interval, which is common among all the models except for the LIDO predictions which start from 1.5$~\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~\ref{tab:chi2}. 
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~\cite{Dover:1991zn}, while TAMU implements a resonance-recombination model~\cite{Ravagli:2007xx}, PHSD a dynamical coalescence via a Monte Carlo approach~\cite{Song:2015sfa}, and POWLANG an in-medium string formation approach~\cite{Beraudo:2014boa}.
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~\cite{He:2019vgs}, MC@sHQ+EPOS2~\cite{Nahrgang:2013xaa}, LIDO~\cite{Ke:2018jem}, LGR~\cite{Li:2019lex}, and  Catania~\cite{Scardina:2017ipo,Plumari:2019hzp}. 
A similar analysis was performed for the elliptic and triangular flow~\cite{Acharya:2020pnh}. 
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.~\cite{Acharya:2020pnh}. Note that some theoretical predictions for the $v_2$ were updated after the publication of Ref.~\cite{Acharya:2020pnh}, namely LBT~\cite{Cao:2016gvr,Cao:2017hhk}, LIDO~\cite{Ke:2018jem}, LGR~\cite{Li:2019lex}, PHSD~\cite{Song:2015sfa}, and TAMU~\cite{He:2019vgs}. 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.~\cite{Acharya:2020pnh}, with the only exception of LIDO~\cite{Ke:2018jem} 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.~\cite{Acharya:2020pnh}.  
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~\cite{He:2019vgs}, MC@sHQ+EPOS2~\cite{Nahrgang:2013xaa}	, LIDO~\cite{Ke:2018jem}, LGR~\cite{Li:2019lex}, and Catania~\cite{Scardina:2017ipo,Plumari:2019hzp}. 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~\mathrm{MeV}$~\cite{Borsanyi:2010bp}. 
According to lQCD calculations the $2\pi D_\mathrm{s} T_\mathrm{c}$ lies between 2 and 6~\cite{Kaczmarek:2014jga,Ding:2012sp,Banerjee:2011ra}, 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.~\cite{Acharya:2020pnh}), STAR (2--12 from Ref.~\cite{Adamczyk:2017xur}), and PHENIX (value of $\sim$3 from Ref.~\cite{Adare:2010de}).
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.~\cite{Acharya:2020pnh}.
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.~\cite{Rapp:2018qla,Cao:2018ews}) playing an important role.

\begin{table}[!tb]
\caption{Summary of the $\chi^2/\mathrm{ndf}$ values obtained in $0<p_{\rm T}<8~\mathrm{GeV}/c$ for the different model predictions compared with the measured D-meson $R_{\rm AA}$.}
\centering
\renewcommand*{\arraystretch}{1.2}
\begin{tabular}[t]{l|c c c}
\toprule
\multirow{3}{*}{Model} & \multicolumn{3}{c}{$\chi^2/\mathrm{ndf}$}\\

\cline{2-4}
& \multicolumn{2}{c}{$R_{\rm AA}$} & \multirow{2}{*}{global} \\
\cline{2-3}
&  0--10\% & 30--50\% & \\
\hline
\hline
Catania~\cite{Scardina:2017ipo,Plumari:2019hzp}		& 94.9$/$15 		& 48.9$/$15 	& 143.8$/$30 \\
DAB-MOD~\cite{Katz:2019fkc}							& 110.2$/$15 	& 123.9$/$15  	& 234.1$/$30 \\
LBT~\cite{Cao:2016gvr,Cao:2017hhk}					& 342.6$/$15 	& 69.2$/$15  	& 411.8$/$30 \\
LIDO~\cite{Ke:2018jem}								& 31.8$/$13		& 14.5$/$13  	& 46.4$/$26 \\
LGR~\cite{Li:2019lex}								& 4.7$/$15		& 4.5$/$15   	& 9.2$/$30  \\
MC@sHQ+EPOS2~\cite{Nahrgang:2013xaa}					& 31.8$/$15		& 24.8$/$15  	& 56.6$/$30  \\
PHSD~\cite{Song:2015sfa}								& 103.2$/$15		& 191.4$/$15 	& 294.7$/$30 \\
POWLANG-HTL~\cite{Beraudo:2014boa,Beraudo:2017gxw}	& 331.0$/$15		& 137.6$/$15  	& 468.6$/$30 \\
TAMU~\cite{He:2019vgs}								& 16.7$/$15		& 13.5$/$15  	& 30.2$/$30  \\
\bottomrule
\end{tabular}
\label{tab:chi2}	
\end{table}

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.~\ref{fig:Dmeson_collrad} and~\ref{fig:Dmeson_fragcoal}.
\begin{figure}[!t]
\begin{center}
\includegraphics[width=1\textwidth]{./figures/raa_transport_models/D_Raa010_V23050_RadColl.eps}
\caption{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~\cite{Ke:2018jem} and LGR~\cite{Li:2019lex} predictions obtained with and without including radiative processes in the  charm-quark interactions with the medium.} 
\label{fig:Dmeson_collrad}
\end{center}
\end{figure}
In Fig.~\ref{fig:Dmeson_collrad} 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~\cite{Ke:2018jem} and LGR~\cite{Li:2019lex} 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$~\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~\cite{He:2019vgs}, PHSD~\cite{Song:2015sfa}, and POWLANG-HTL~\cite{Beraudo:2014boa,Beraudo:2017gxw}. Radiative processes, instead, become dominant at intermediate and high $p_{\rm T}$. This can be observed in Fig.~\ref{fig:Dmeson_collrad} where the predictions without these processes overestimate (underestimate) the measured $R_{\rm AA}$ ($v_2$) for $p_{\rm T} \gtrsim$ 5--6$~\mathrm{GeV}/c$.

\begin{figure}[!t]
\begin{center}
\includegraphics[width=1\textwidth]{./figures/raa_transport_models/D_Raa010_V23050_FragCoal_3models.eps}
\caption{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~\cite{Song:2015sfa}, POWLANG~\cite{Beraudo:2014boa,Beraudo:2017gxw}, and DAB-MOD~\cite{Katz:2019fkc} predictions obtained with and without including hadronisation via recombination.} 
\label{fig:Dmeson_fragcoal}
\end{center}
\end{figure}
Similarly, Fig.~\ref{fig:Dmeson_fragcoal} shows the comparison of the experimental data for $R_{\rm AA}$ and $v_2$ with two different versions of the 
PHSD~\cite{Song:2015sfa}, POWLANG-HTL~\cite{Beraudo:2014boa,Beraudo:2017gxw}, and DAB-MOD~\cite{Katz:2019fkc} 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.~\ref{fig:Dmeson_fragcoal}, 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.

\section{Systematic uncertainties}
\label{sec:systematics}

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~\ref{tab:syst} for the three D-meson species in representative $p_{\rm T}$ intervals.


\begin{table}[!tb]
\begin{center}
\caption{Relative systematic uncertainties on prompt D-meson yields in Pb--Pb collisions at $\sqrt{s_{\rm \scriptscriptstyle NN}} = 5.02~\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.}
\begin{tabular}{|c|l|ccc|cc|cc|}
\hline
\multicolumn{2}{|r|}{\rule{0pt}{2.6ex} Particle} & \multicolumn{3}{c|}{${\rm D^0}$} & \multicolumn{2}{c|}{${\rm D^+}$}  & \multicolumn{2}{c|}{${\rm D^{*+}}$} \\[0.5ex]
\hline
\multicolumn{2}{|r|}{\rule{0pt}{2.6ex} $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.5ex]
\hline
\hline
\multirow{9}{*}{\parbox{2 cm}{\centering 0--10\%\\ centrality}} 
& Raw-yield extraction \rule{0pt}{2.6ex}            & 9\%  & 7\%  & 2\%             & 5\%  & 5\%                     & 7\%  & 3\%\\
& Correction factor & & & & & & &\\
& $\qquad$ Tracking efficiency   & 8\%  & 9\%  & 9\%         & 14\%  & 15.5\%               & 13\%  & 9\%\\
& $\qquad$ Selection efficiency  & negl.    & 8\%  & 2\%             & 5\%  & 5\%                     & 13\%  & 4\%\\
& $\qquad$ PID efficiency        & negl.    & negl.    & negl.             & negl.  & negl.                       & 1\%  & 1\%\\
& $\qquad$ MC $p_{\rm T}$ shape        & negl.    & 0.5\%  & negl.               & negl.  & negl.                       & 0.5\%  & negl. \\[1ex]
& 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}\%$\\[1ex]	
\cline{3-9} 
& Branching ratio \rule{0pt}{2.6ex}      & \multicolumn{3}{c|}{0.8\%}      &\multicolumn{2}{c|}{1.7\%}        &\multicolumn{2}{c|}{1.1\%}\\
\cline{3-9} 
& Centrality limits \rule{0pt}{2.6ex}    & \multicolumn{7}{c|}{<0.1\%}\\
\hline
\hline
\multirow{9}{*}{\parbox{2 cm}{\centering 30--50\%\\ centrality}} 
& Raw-yield extraction \rule{0pt}{2.6ex}           & 11\%  & 4\%  & 2\%             & 4\%  & 3\%                     & 5\%  & 4\%\\
& Correction factor & & & & & & &\\
& $\qquad$ Tracking efficiency   & 7\%  & 7.5\%  & 7\%           & 11\%  & 11\%               & 10\%  & 9\%\\
& $\qquad$ Selection efficiency  & negl.    & 8\%  & 2\%             & 4\%  & 3\%                     & 9\%  & 4\%\\
& $\qquad$ PID efficiency        & negl.    & negl.    & negl.           & negl.    & negl.                   & negl.    & negl. \\
& $\qquad$ MC $p_{\rm T}$ shape        & negl.    & 1\%  & negl.               & negl.  & negl.                     & negl.    & negl. \\[1ex]
& 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}\%$\\[1ex]	
\cline{3-9} 
& Branching ratio \rule{0pt}{2.6ex}      & \multicolumn{3}{c|}{0.8\%}      &\multicolumn{2}{c|}{1.7\%}        &\multicolumn{2}{c|}{1.1\%}\\
\cline{3-9} 
& Centrality limits \rule{0pt}{2.6ex}    & \multicolumn{7}{c|}{2\%}\\
\hline
\end{tabular}
\label{tab:syst}
\end{center}
\end{table}

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~\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.~\cite{Cacciari:2012ny}) 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.~\cite{Adam:2015sza} 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 $\av{T_{\rm AA}}$ 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. 

