Measurement of D$^0$, D$^+$, D$^{*+}$ and D$^+_{\rm s}$ production in Pb-Pb collisions at $\mathbf{\sqrt{s_{\rm NN}}}= 5.02$ TeV

We report measurements of the production of prompt D$^0$, D$^+$, D$^{*+}$ and D$^+_{\rm s}$ mesons in Pb-Pb collisions at the centre-of-mass energy per nucleon-nucleon pair $\sqrt{s_{\rm NN}}=5.02$ TeV, in the centrality classes 0-10%, 30-50% and 60-80%. The D-meson production yields are measured at mid-rapidity ($|y|<0.5$) as a function of transverse momentum ($p_{\rm T}$). The $p_{\rm T}$ intervals covered in central collisions are: $1<p_{\rm T}<50$ Gev/$c$ for D$^0$, $2<p_{\rm T}<50$ GeV/$c$ for D$^+$, $3<p_{\rm T}<50$ GeV/$c$ for D$^{*+}$, and $4<p_{\rm T}<16$ GeV/$c$ for D$^+_{\rm s}$ mesons. The nuclear modification factors ($R_{\rm AA}$) for non-strange D mesons (D$^0$, D$^+$, D$^{*+}$) show minimum values of about 0.2 for $p_{\rm T}$ = 6-10 GeV/$c$ in the most central collisions and are compatible within uncertainties with those measured at $\sqrt{s_{\rm NN}}=2.76$ TeV. For D$^+_{\rm s}$ mesons, the values of $R_{\rm AA}$ are larger than those of non-strange D mesons, but compatible within uncertainties. In central collisions the average $R_{\rm AA}$ of non-strange D mesons is compatible with that of charged particles for $p_{\rm T}>8$ GeV/$c$, while it is larger at lower $p_{\rm T}$. The nuclear modification factors for strange and non-strange D mesons are also compared to theoretical models with different implementations of in-medium energy loss.


Introduction
Ultra-relativistic collisions of heavy nuclei produce a state of strongly-interacting matter characterised by high energy density and temperature. According to Quantum Chromodynamics (QCD) on the lattice, in these extreme conditions matter undergoes a phase transition to a Quark-Gluon Plasma (QGP) state in which quarks and gluons are deconfined and chiral symmetry is partially restored [1][2][3][4].
Heavy quarks (such as charm and beauty) are predominantly produced in the early stage of the collision in hard scattering processes between partons of the incoming nuclei. Because of their large masses, their production time (∼ 0.1 and 0.02 fm/c for charm and beauty, respectively [5]) is shorter than the formation time of the QGP, which is about 0.3-1.5 fm/c at Large Hadron Collider (LHC) energies [6]. In contrast, the thermal production and annihilation rates are negligible [7]. Heavy quarks therefore experience the full evolution of the hot and dense QCD medium.
During their propagation through the medium, heavy quarks are exposed to interactions with the medium constituents and lose part of their energy via inelastic (gluon radiation) [8,9] or elastic scatterings (collisional processes) [10][11][12]. The colour-charge dependence of the strong interaction and parton-massdependent effects are predicted to influence the amount of energy loss (see [5,13] for recent reviews). Low-momentum heavy quarks can participate in the collective expansion of the system as a consequence of multiple interactions with the medium [14,15]. It was also suggested that low-momentum heavy quarks could hadronise not only via fragmentation in the vacuum, but also via the mechanism of recombination with other quarks in the medium [15,16]. In this scenario, the large abundance of strange quarks in nucleus-nucleus collisions with respect to proton-proton collisions is expected to lead to an increased production of D + s mesons relative to non-strange D mesons [17]. The effects of energy loss and the dynamics of heavy-quark hadronisation can be studied using the nuclear modification factor R AA , which compares the transverse-momentum (p T ) differential production yields in nucleus-nucleus collisions (dN AA /dp T ) with the cross section in proton-proton collisions (dσ pp /dp T ) scaled by the average nuclear overlap function T AA : R AA (p T ) = 1 T AA · dN AA /dp T dσ pp /dp T .
The average nuclear overlap function T AA is defined as the average number of nucleon-nucleon collisions N coll , which can be estimated via Glauber model calculations [18][19][20][21], divided by the inelastic nucleon-nucleon cross section.
Measurements of prompt D-meson production by the ALICE Collaboration in Pb-Pb collisions at √ s NN = 2.76 TeV [22][23][24][25] showed a strong suppression of the D-meson yields by a factor of 5-6 for 8 < p T < 12 GeV/c in the 10% most central collisions. Recent results from the CMS Collaboration on D 0 production in the p T range 2-100 GeV/c show a similar suppression for 6 < p T < 10 GeV/c in the 10% most central Pb-Pb collisions at √ s NN = 5.02 TeV, decreasing with increasing p T [26]. In contrast, the D-meson nuclear modification factor in p-Pb collisions at √ s NN = 5.02 TeV, where an extended QGP phase is not expected to be formed, was found to be consistent with unity within uncertainties for 0 < p T < 24 GeV/c [27]. These results indicate that the strong suppression is due to substantial final-state interactions of charm quarks with the QGP formed in Pb-Pb collisions.
In this article, we present the measurement of p T -differential yields and the nuclear modification factor for prompt D 0 , D + , D * + and D + s mesons (including their antiparticles), in Pb-Pb collisions at √ s NN = 5.02 TeV collected with the ALICE detector during the LHC Run 2 in 2015. Prompt D mesons are defined as those produced by the hadronisation of charm quarks or from the decay of excited open charm and charmonium states, hence excluding the decays of beauty hadrons. The experimental apparatus is briefly presented in Section 2, together with the data sample used for the analysis. The reconstruction of D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration Centrality class T AA (mb −1 ) N events 0-10% 23.42 ± 0.75 10.4 × 10 6 30-50% 3.82 ± 0.14 20.8 × 10 6 60-80% 0.404 ± 0.017 20.8 × 10 6 Table 1: Average nuclear overlap function and number of events for the three centrality classes used in the analysis.
D-meson hadronic decays and all corrections applied to the raw yields are presented in Section 3. The procedure used to obtain the proton-proton reference cross section at √ s = 5.02 TeV and the estimation of the systematic uncertainties are described in Section 4 and Section 5, respectively. The results for the central (0-10%), semi-central (30-50%) and peripheral (60-80%) collisions are presented in Section 6. A comparison with charged-pion and charged-particle R AA is reported in the same Section, along with detailed comparisons with model calculations, including a simultaneous comparison of the R AA and elliptic flow v 2 . Conclusions are drawn in Section 7.

Experimental apparatus and data sample
A description of the ALICE experimental apparatus and its performance in pp, p-Pb and Pb-Pb collisions can be found in [28,29]. The main detectors used in the present analysis are the V0 detector, the Inner Tracking System (ITS) [30], the Time Projection Chamber (TPC) [31] and the Time Of Flight (TOF) detector [32], located inside a large solenoidal magnet providing a uniform magnetic field of 0.5 T parallel to the LHC beam direction (z axis in the ALICE reference system), and the Zero Degree Calorimeters (ZDC) [33], located at z = ±112.5 m from the nominal interaction point. The analysed sample consists of Pb-Pb collision data recorded with a minimum-bias interaction trigger that required coincident signals in both scintillator arrays of the V0 detector [34]. The V0 detector consists of two scintillator arrays, which cover the full azimuth in the pseudorapidity intervals −3.7 < η < −1.7 and 2.8 < η < 5.1. Events produced by 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 reconstructed interaction point (primary vertex) within ±10 cm from the centre of the ITS detector along the beam line were used in the analysis. For the data sample considered in this paper, the probability of in-bunch collision pileup (i.e. collisions with two or more simultaneous interactions per bunch crossing) was negligible, while the request of at least a hit in one of the two innermost layers of the ITS rejected tracks produced in out-of-bunch pileup collisions.
Collisions were divided into centrality classes, determined from the sum of the V0 signal amplitudes and defined in terms of percentiles of the hadronic Pb-Pb cross section. In order to relate the centrality classes to the collision geometry, the distribution of the V0 summed amplitudes was fitted with a function based on the Glauber model [18][19][20][21] combined with a two-component model for particle production [35], which decomposes particle production in nucleus-nucleus collisions into the contributions due to soft and hard interactions. The centrality classes used in the present analysis, together with the corresponding average nuclear overlap function T AA [36] and the number of events (N events ) in each class, are summarised in Table 1. The corresponding integrated luminosity is about L int ≈ 13 µb −1 [37].

Data analysis
The D mesons and their charge conjugates were reconstructed in the decay channels D 0 → K − π + (with branching ratio, BR, of (3.93 ± 0.04)%), D + → K − π + π + (BR of (9.46 ± 0.24)%), D * + → D 0 π + (BR of (67.7 ± 0.5)%) and D + s → φ π + → K + K − π + (BR of (2.27 ± 0.08)%) [38]. D 0 , D + and D + s candidates were defined using pairs and triplets of tracks with proper charge-sign combination having |η| < 0.8, p T > 0.4 GeV/c, a minimum number of 70 (out of 159) associated space points in the TPC and at least two hits (out of six) in the ITS, with at least one in the two innermost layers. D * + candidates were formed D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration by combining D 0 candidates with tracks having |η| < 0.8, p T > 0.1 GeV/c and at least three associated hits in the ITS. For D + s candidate selection, one of the two pairs of opposite-sign tracks was required to have an invariant mass compatible with the φ mass (m φ = 1019.461 ± 0.019 MeV/c 2 [38]). In particular, the difference between the reconstructed K + K − invariant mass and φ mass was required to be less than 5-10 MeV/c 2 depending on the D + s p T interval. This selection preserves 70-85% of the D + s signal. The selection of tracks with |η| < 0.8 limits the D-meson acceptance in rapidity, which, depending on p T , varies from |y| < 0.6 for p T = 1 GeV/c to |y| < 0.8 for p T > 5 GeV/c. A p T -dependent fiducial acceptance cut, |y D | < y fid (p T ), was therefore applied to the D-meson rapidity. The value of y fid (p T ) increases from 0.6 to 0.8 in the range 1 < p T < 5 GeV/c, and the variation can be described according to a second-order polynomial function. For p T > 5 GeV/c one has y fid = 0.8.
The selection strategy is similar to the one used in previous analyses [25,39] and is mainly based on the separation between primary and secondary vertex, the displacement of the tracks from the primary vertex and the pointing of the reconstructed D-meson momentum to the primary vertex. In comparison to previous analyses, additional selection criteria were exploited. In particular, the normalised difference between the measured and expected transverse-plane impact parameters of each of the decay particles (already introduced in [40]) and the transverse-plane impact parameter to the primary vertex (d xy 0 ) of the D-meson candidates were used. Besides the rejection of the combinatorial background, a selection based on the latter two variables has the advantage to suppress significantly the fraction of D mesons coming from beauty-hadron decays (feed-down) and hence reduce the associated systematic uncertainty. The cut values on the selection variables were optimised in each centrality class independently, in order to obtain a large statistical significance of the D-meson signals, while keeping the selection efficiency of promptly produced D mesons as large as possible. Further background reduction was obtained by applying particle identification for charged pions and kaons with the TPC and TOF detectors. A ±3 σ window around the expected mean values of specific ionisation energy loss dE/dx in the TPC gas and time-of-flight from the interaction point to the TOF detector was used for the identification, where σ is the resolution on these two quantities. In central collisions, a 2 σ selection was used for D * + and D + (for p T < 3 GeV/c) candidates. For D + s candidates, tracks without a TOF signal (mostly at low momentum) were identified using only the TPC information and requiring a 2 σ compatibility with the expected dE/dx. The stricter PID selection strategy was needed due to the large background of track triplets and, in case of D + s , because of its short lifetime, which limits the effectiveness of the geometrical selections on the displaced decay-vertex topology.
The D 0 , D + and D + s raw yields were obtained from binned maximum-likelihood fits to the candidate invariant-mass (M) distributions, while for the D * + the mass difference ∆M = M(Kππ) − M(Kπ) distributions were used. Examples for these distributions are shown in Fig. 1 for the centrality class 0-10%. The D 0 , D + and D + s candidate invariant-mass distributions were fitted with a function composed of a Gaussian term for the signal and an exponential function to describe the background shape, with the exception of the D 0 p T intervals 1-2 GeV/c and 2-3 GeV/c, where the background was found to be better described by a second-order polynomial function (a fourth-order polynomial was used in 1-2 GeV/c for the 0-10% centrality class). The ∆M distribution of D * + candidates was fitted with a Gaussian function for the signal and a threshold function multiplied by an exponential for the background (a √ ∆M − m π · e b(∆M−m π ) , where m π is the pion mass and a and b are free parameters). The contribution of signal candidates that are present in the invariant-mass distribution of the D 0 meson with the wrong decay-particle mass assignment (reflection), was parametrised by fitting the simulated reflection invariant-mass distributions with a double Gaussian function, and it was included in the total D 0 fit function. The ratio between the reflected signal and the yields of the D 0 was taken from simulations (typically 2-5% of the raw yield, depending on p T ) [39]. The Monte Carlo simulation used for this study is the same one used to determine the reconstruction efficiency, as described in the following dedicated paragraph. In addition, given the critical signal extraction induced by the small signal-to-background      Figure 1: Invariant-mass distributions for the four D-meson species in selected p T intervals for the centrality class 0-10%. Fitted values for the meson mass µ, width σ and raw yield S are also given. Top row: D 0 mesons with 1 < p T < 2 GeV/c, before (left) and after (right) subtraction of the background fit function. Middle row: D + mesons with 8 < p T < 10 GeV/c and D * + mesons (difference of M(Kππ) and M(Kπ)) with 24 < p T < 36 GeV/c. Bottom row: D + s mesons with 4 < p T < 6 GeV/c and 12 < p T < 16 GeV/c; the D + → K + K − π + signal is visible on the left of the D + s signal.
ratio of the D 0 meson in 1 < p T < 2 GeV/c, the width of the Gaussian used to describe the signal was fixed to the value obtained in the simulations. The Gaussian widths obtained from the simulations were found to be consistent with those extracted from the data in the full p T range, for all measured centrality classes, with deviations of at most 10-15%. In the fit to the D + s -candidate invariant-mass distribution, an additional Gaussian was used to describe the D + → K + K − π + signal on the left of the D + s signal. The statistical significance S/ √ S + B of the observed signals, estimated within 3 standard deviations, varies from 5 to 33 depending on the D-meson species, the p T interval, and the centrality class.
The D-meson raw yields were corrected in order to obtain the p T -differential yields of prompt D mesons dN D dp T |y|<0. 5 = (2) √ s NN = 5.02 TeV ALICE Collaboration The raw yields N D+D raw were divided by a factor of two to obtain the charge-averaged (particle and antiparticle) yields. To correct for the contribution of feed-down from beauty-hadron decays, the raw yields were multiplied by the fraction of promptly produced D mesons, f prompt (see Eq. 3). Furthermore, they were divided by the product of prompt D-meson acceptance and efficiency (Acc × ε) prompt , by the branching ratio BR of the decay channel, by the transverse momentum interval width ∆p T and by the number of events N events . The factor α y (p T ) = y fid (p T )/0.5 normalises the corrected yields measured in |y| < y fid (p T ) to one unit of rapidity |y| < 0.5, assuming a flat rapidity distribution for D mesons in |y| < y fid (p T ). This assumption was validated to the 1% level with simulations for pp collisions [41,42] and it is justified also for Pb-Pb collisions. For example, measurements of the prompt and non-prompt J/ψ R AA in Pb-Pb collisions at √ s NN = 2.76 TeV do not exhibit a significant rapidity dependence [43].
The correction for acceptance and efficiency (Acc × ε) prompt was determined using Monte Carlo simulations with a detailed description of the detector and its response, based on the GEANT3 transport package [44]. The underlying Pb-Pb events at √ s NN = 5.02 TeV were simulated using the HIJING v1.383 generator [45] and D-meson signals were added using the PYTHIA v6.421 generator [46] with Perugia-2011 tune. Each simulated PYTHIA pp event contained a cc or bb pair, and D mesons were forced to decay into the hadronic channels of interest for the analysis. In the most central event class, the p T distribution of D mesons in the MC simulation for p T > 2 GeV/c was weighted in order to match the shape measured in data for D 0 mesons in finer p T intervals with respect to those used in the analysis. In the centrality classes and p T ranges where an analysis in finer p T intervals was not possible, the simulated D-meson p T distribution was weighted to match the shape given by model calculations. In particular, fixed-order plus next-to-leading-log perturbative QCD calculations (FONLL) [47,48] multiplied by the R AA (p T ) of D mesons computed using the BAMPS model (which implements both elastic and radiative processes) for the 30-50% centrality class [49][50][51] were used for the corresponding centrality class. For the p T intervals 1-2 GeV/c and 16-50 GeV/c in the 0-10% centrality class and for the 60-80% centrality class, where the R AA is nearly flat in the measured p T interval, only the FONLL calculations were used. Figure 2 shows the acceptance-times-efficiency (Acc × ε) for prompt and feed-down D mesons with rapidity |y| < y fid (p T ) in the centrality class 0-10%, after the aforementioned p T -distribution weighting procedure. The difference between the (Acc × ε) factor for prompt and feed-down D-mesons arises from the geometrical selections applied, given the different decay topology of D mesons coming from B decays. In particular, the feed-down D mesons are on average more displaced from the primary vertex due to the large B-meson lifetime (cτ ≈ 500 µm [38]) and therefore are more efficiently selected by the majority of the analysis cuts (e.g. for D 0 and D * + in most of the p T intervals). On the contrary, the selections on the difference between measured and expected decay-track impact parameters and on the D-meson impact parameter reject more feed-down D mesons, thus reducing the feed-down efficiencies as compared to the previous analyses (e.g. for D + and D + s ). The (Acc × ε) is higher for more peripheral collisions, by up to a factor larger than two at low p T , since less stringent selections can be applied because of the lower combinatorial background.
The f prompt factor was obtained, following the procedure introduced in [22], by subtracting the contribution of D mesons from beauty-hadron decays from the measured raw yield in each p T interval. It was estimated using perturbative QCD calculations, efficiencies from MC simulations, and an hypothesis on the R AA of feed-down D mesons. The expression for f prompt reads: In this expression, N D+D raw is the measured raw yield and N D+Dfeed-down raw is the estimated raw yield of D D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration mesons from beauty-hadron decays. In detail, the beauty-hadron production cross section in pp collisions at √ s = 5.02 TeV, estimated with FONLL calculations [52], was folded with the beauty-hadron→ D + X decay kinematics using the EvtGen package [53] and multiplied by T AA of the corresponding centrality class, by the (Acc × ε) for feed-down D mesons, and by the other factors introduced in Eq. (2). In addition, the nuclear modification factor of D mesons from beauty-hadron decays was accounted for. The comparison of the R AA of prompt D mesons (R prompt AA ) at √ s NN = 2.76 TeV [24] with that of J/ψ from B-meson decays [43] at the same energy measured by the CMS Collaboration indicates that prompt charmed hadrons are more suppressed than non-prompt charmed hadrons. The R AA values differ by a factor of about two in central collisions at a transverse momentum of about 10 GeV/c [24] and this difference is described by model calculations with parton-mass-dependent energy loss. Therefore, for the centrality classes 0-10% and 30-50%, the value R feed-down AA = 2 · R prompt AA was used to compute the correction for non-strange D mesons with 3 < p T < 24 GeV/c. This hypothesis was varied in the range 1 < R feed-down AA /R prompt AA < 3 considering the data uncertainties and model variations to estimate a systematic uncertainty. For 1 < p T < 3 GeV/c and 24 < p T < 50 GeV/c, where model calculations predict a reduced difference between the R AA values of prompt and non-prompt charm hadrons [54], the hypothesis R feed-down 4 Proton-proton reference for R AA The p T -differential cross sections of prompt D mesons with |y| < 0.5 in pp collisions at √ s = 5.02 TeV, used as reference for the nuclear modification factor, were obtained by scaling the measurements at √ s = 7 TeV [40] to √ s = 5.02 TeV with FONLL calculations [52]. These measurements reach up to p T = 36 GeV/c for D 0 , 24 GeV/c for D + and D * + , and 12 GeV/c for D + s mesons. The uncertainties on the p T -dependent scaling factor from √ s = 7 TeV to √ s = 5.02 TeV were determined by varying the FONLL parameters (charm-quark mass, factorisation and renormalisation scales) as described in [55]. The uncertainties range from +17 − 4 % for 1 < p T < 2 GeV/c to about ±3% for p T > 10 GeV/c. At high D-meson p T (36 < p T < 50 GeV/c for D 0 , 24 < p T < 50 GeV/c for D + and D * + , and 12 < p T < 16 GeV/c for D + s ), the FONLL calculation at √ s = 5.02 TeV [52] was used as a reference by scaling the values for each meson species to match the central value of the scaled data at lower p T . This procedure is described in Ref. [23]. As an example, the total systematic uncertainties on the pp reference for D 0 mesons with 36 < p T < 50 GeV/c is +38 −28 %.

Systematic uncertainties
Systematic uncertainties on the D-meson yield in Pb-Pb collisions were estimated considering the following sources: (i) extraction of the raw yield from the invariant-mass distributions; (ii) track reconstruction efficiency; (iii) D-meson selection efficiency; (iv) PID efficiency; (v) generated D-meson p T shape in the simulation; (vi) subtraction of the feed-down from beauty-hadron decays. In addition, the uncertainties on the branching ratios [38] were considered. A procedure similar to that described in [22][23][24][25] and outlined in what follows was used to estimate the uncertainties as a function of p T and centrality. The systematic uncertainties on the raw yield extraction were evaluated for each D-meson species and in each p T interval by varying the lower and upper limits of the fit range, and the background fit function. In addition, the same approach was used with a bin-counting method, in which the signal yield was obtained by integrating the invariant-mass distribution after subtracting the background estimated from a fit to the side-bands. It ranges between 2% and 15% depending on the D-meson species and p T interval. In the case of D 0 , an additional contribution due to signal reflections in the invariant-mass distribution was estimated by varying the ratio of the integral of the reflections over the integral of the signal and the shape of the templates used in the invariant-mass fits. For the D 0 meson in the interval 1 < p T < 2 GeV/c, the signal line shape was varied by using Gaussian functions with the widths fixed to ±15% with respect to the value expected from Monte Carlo simulations, based on the deviations between the Gaussian width values observed in data and simulations. For the four D mesons, further checks on the stability of the results were performed by repeating the fits varying the invariant-mass bin width.
The systematic uncertainty on the track reconstruction efficiency was estimated by varying the trackquality selection criteria and by comparing the probability to match the TPC tracks to the ITS hits in data and simulation. The comparison of the matching efficiency in data and simulations was made after weighting the relative abundances of primary and secondary particles in the simulation to match those observed in data, which were estimated via fits to the inclusive track impact parameter distributions. The estimated uncertainty depends on the D-meson p T and ranges from 3% to 8% for the two-body decay of D 0 mesons and from 6% to 12% for the three-body decays of D + , D * + and D + s mesons. To estimate the uncertainty on the PID selection efficiency, for the three non-strange D-meson species the analysis was repeated without PID selection. The resulting cross sections were found to be compatible D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration with those obtained with the PID selection and therefore no systematic uncertainty was assigned. For the D + s meson, the lower signal yield and the larger combinatorial background prevented a signal estimation without particle identification. In this case, a 3% uncertainty was estimated by repeating the analysis with a 3 σ PID selection, similar to that used for non-strange D-mesons for which no systematic effects were observed. This value was also verified by comparing the pion and kaon PID selection efficiencies in the data and in the simulation and combining the observed differences using the D + s decay kinematics (for this test, pure pion samples were selected using strange hadron decays, while kaon samples in the TPC were obtained using a tight PID selection in the TOF).
The uncertainty on the D-meson selection efficiency (see Cut efficiency in Table 2) originates from imperfections in the description of the D-meson kinematic properties and of the detector resolutions and alignments in the simulation. It 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 efficiencies, raw yield and background values. The assigned uncertainty for non-strange D mesons is 5% in most of the p T intervals and it increases to 10-15% in the lowest p T intervals, where the efficiencies are low and vary steeply with p T , because of the tighter selections. A larger uncertainty of 10% in all p T intervals was estimated for D + s mesons, for which more stringent selection criteria were utilized in the analysis as compared to non-strange D mesons.
The systematic effect on the efficiency due to a possible difference between the real and simulated D-meson transverse momentum distributions was estimated by using alternative D-meson p T distributions. In particular, the p T distributions from FONLL calculations with and without hot-medium effects parametrised based on the R AA in central collisions from the BAMPS [56], LBT [57] and TAMU [58] models were used in this study. The uncertainty, which also includes the effect of the p T dependence of the nuclear modification factor, was estimated to be, for non-strange D mesons in central collisions, about 10% in the lowest p T intervals and decreasing to zero for p T > 5 GeV/c. For D + s mesons the uncertainty was estimated as 7% in 4-6 GeV/c, 2% in 6-8 GeV/c and 1% at higher p T .
The systematic uncertainty on the subtraction of feed-down from beauty-hadron decays (i.e. the calculation of the f prompt fraction) was estimated by varying i) the p T -differential feed-down D-meson cross section from the FONLL calculation within the theoretical uncertainties, ii) the ratio of the feed-down and prompt D-meson R AA in the ranges described at the end of Section 3. The resulting uncertainty ranges between 2% and 15%, depending on D-meson species, centrality classes and p T intervals.
The systematic uncertainties on the p T -differential spectra and R AA in the two extreme centrality classes are listed for all D-meson species in Table 2 for the lowest p T interval accessible as well as for the intermediate range 7 < p T < 8 GeV/c (6 < p T < 8 GeV/c for the D + s meson). The systematic uncertainties on the R AA measurement include those on the D-meson corrected yields described above, those on the proton-proton reference cross section, and the uncertainties on the average nuclear overlap function.
The systematic uncertainty on the pp reference used for the calculation of R AA has two contributions. The first one is the systematic uncertainty on the measured p T -differential D-meson cross section at √ s = 7 TeV. The second contribution is the scaling to √ s = 5.02 TeV, which has been discussed in Section 4.
In the calculation of the nuclear modification factor, the systematic uncertainty on the feed-down subtraction deriving from the variation of the parameters of the FONLL calculation was considered to be correlated in the Pb-Pb and pp measurements, while all the other sources of systematic uncertainties were treated as uncorrelated.
The uncertainties on the R AA normalisation are the quadratic sum of (i) the pp normalisation uncertainty (3.5%), (ii) the uncertainty on T AA , which ranges from 3.2% to 4.2% depending on the centrality, and  Table 2: Relative systematic uncertainties on the dN/dp T in Pb-Pb collisions, on the extrapolated dN/dp T in pp collisions and on the R AA of D 0 , D * + , D + , and D + s in two centrality classes considered in the analysis for the lowest accessible p T intervals and for the intermediate range 7 < p T < 8 GeV/c (6 < p T < 8 GeV/c for the D + s meson). D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration (iii) the variation of raw yield (< 0.1%, 2% and 3% for the 0-10%, 30-50% and 60-80% centrality classes, respectively) obtained when the centrality intervals are varied to account for the uncertainty on the fraction of the hadronic cross section used in the Glauber fit to determine the centrality [23], and the branching ratio uncertainty cancels out in the ratio.

Results
The transverse-momentum distributions dN/dp T of prompt D 0 , D + , D * + and D + s mesons are shown in Fig. 3 for the 0-10%, 30-50% and 60-80% centrality classes. The vertical bars represent the statistical uncertainties and the empty boxes the systematic uncertainties. The uncertainty on the branching ratios is quoted separately. Their average was computed using the inverse of the quadratic sum of the relative statistical and uncorrelated systematic uncertainties as weights, in the p T intervals where more than one D-meson species is available, (Fig. 5, right-hand panels). The systematic uncertainties were propagated through the averaging procedure, considering the contributions from the tracking efficiency, the beauty-hadron feed-down subtraction and the FONLL-based √ s-scaling of the pp cross section from √ s = 7 TeV to √ s = 5.02 TeV as fully correlated uncertainties among the three D-meson species. The average nuclear modification factors in the 0-10% and 30-50% centrality classes (top and middle righthand panels of Fig. 5) show a suppression that is maximal at p T = 6-10 GeV/c, where a reduction of the yields by a factor of about 5 and 2.5 with respect to the binary-scaled pp reference is observed in the two centrality classes, respectively. The suppression gets smaller with decreasing p T for p T < 6 GeV/c, and R AA is compatible with unity in the interval 1 < p T < 3 GeV/c. The average R AA in the 60-80% centrality class shows a suppression by about 20-30%, without a pronounced dependence on p T .
The R AA of prompt D + s mesons is shown in the right-hand panels of Fig. 5, where it is compared with the average R AA of non-strange D mesons: the values are larger for D + s mesons, but the two measurements are compatible within one standard deviation of the combined uncertainties, as is the case for the ratios shown in Fig. 4. In Fig. 6, the nuclear modification factor of D mesons at √ s NN = 5.02 TeV in the 0-10% centrality class is compared with the same measurement at √ s NN = 2.76 TeV [23]. The measurement at √ s NN = 5.02 TeV have total uncertainties reduced by a factor of about two and extended p T coverage from 36 to 50 GeV/c. The suppression is compatible within uncertainties at the two energies, as also observed for charged particles [59].
The close similarity of the R AA measurements at the two energies was predicted by the Djordjevic model [54] (Fig. 6, right panel), and it results from the combination of a higher medium temperature at 5.02 TeV (estimated to be about 7% higher than at 2.76 TeV), which would decrease the R AA by about D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration 10%, with a harder p T distribution of charm quarks at 5.02 TeV, which would increase the R AA by about 5% if the medium temperature were the same as at 2.76 TeV.
As explained in Section 1, the measurement of the R AA of open-charm mesons is essential to understand in-medium parton energy loss, in particular its colour-charge and quark-mass dependence. In Fig. 7, the R AA of prompt D mesons is compared with that of charged particles in the same p T intervals, at the same energy and in the same centrality classes [59]. The ratio of their nuclear modification factors is displayed in the bottom panels, for the three centrality classes. The R AA of D mesons and charged particles differ D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration by more than 2 σ of the combined statistical and systematic uncertainties in all the p T intervals within 3 < p T < 8 GeV/c in central collisions. The difference is less than 2 σ in this range for semi-central collisions, while the two R AA are the same within 1 σ for p T > 10 GeV/c in both central and semicentral collisions. In the 60-80% class the measurements are compatible in the common p T interval. The interpretation of the difference observed for p T < 8 GeV/c in central and semi-central collisions is not straightforward, because several factors can play a role in defining the shape of the R AA .
In presence of a colour-charge and quark-mass dependent energy loss, the harder p T distribution and the harder fragmentation function of charm quarks compared to those of light quarks and gluons should lead to similar values of D-meson and pion R AA , as discussed in [60]. Since the pions are the dominant contribution in the inclusive charged-particle yields, this statement is expected to be still valid for the comparison of the D-meson and the charged particle R AA . In addition, it should be considered that the yield of light-flavour hadrons could have a substantial contribution up to transverse momenta of about 2-3 GeV/c from soft production processes, such as the break-down of participant nucleons into quarks and D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration gluons that subsequently hadronise. This component scales with the number of participants rather than the number of binary collisions. Finally, the effects of radial flow and hadronisation via recombination, as well as initial-state effects, could affect D-meson and light-hadron yields differently at a given p T .
The average R AA of the three non-strange D-meson species in the three centrality classes are compared with theoretical models in Fig. 8. Models based on heavy-quark transport and models based on perturbative QCD calculations of high-p T parton energy loss are shown in the left and in the right panels, respectively. Transport models in the left panels include: BAMPS el. [56], POWLANG [61] and TAMU [58], in which the interactions are only described by collisional (i.e. elastic) processes; BAMPS el.+rad. [56], LBT [57], MC@sHQ+EPOS2 [62] and PHSD [63], in which also energy loss from medium-induced gluon radiation is considered, in addition to collisional process. In the right panels, the CUJET3.0 [64] and Djordjevic [54] models include both radiative and collisional energy loss processes, while the SCET [65] model implements medium-induced gluon radiation via modified splitting functions with finite quark masses 1 . All models, with the exception of BAMPS and CUJET3.0, include a nuclear modification of the parton distribution functions. The LBT, MC@sHQ, PHSD, POWLANG and TAMU models include a contribution of hadronisation via quark recombination, in addition to independent fragmentation. Most of the models provide a fair description of the data in the region p T < 10 GeV/c in central collisions (except for BAMPS el., where the radiative term is missing), but many of them (LBT, PHSD, POWLANG and SCET) provide a worse description of non-central collisions. In the high-p T region above 10 GeV/c only the BAMPS el.+rad., CUJET3.0, Djordjevic, MC@sHQ+EPOS2 and SCET models can describe the data in central collisions. The CUJET3.0 and Djordjevic models provide a fair description of the R AA in all three centrality classes for p T > 10 GeV/c, where radiative energy loss is expected to be the dominant interaction mechanism, suggesting that the dependence of radiative energy loss on the path length in the hot and dense medium is well understood.
In Fig. 9, the non-strange and strange D-meson R AA are compared with the models that provide both observables. An increase of the D + s R AA is expected in the two models, PHSD and TAMU, in particular 1 The SCET curves reported here differ from those of Ref. [65] because the latter used an extrapolation of the charged-particle multiplicity at for p T < 5 GeV/c, with respect to non-strange D mesons. This increase is induced by hadronisation via quark recombination in the QGP, as well as by different interaction cross sections for non-strange D and for D + s in the hadronic phase of the system evolution. In the transverse momentum interval covered by the D + s measurement (p T > 4 GeV/c), the PHSD model predicts the effect to be very small, while the TAMU model predicts a sizeable difference of about 30% up to about 8 GeV/c, similar to the trend shown by the data.
The simultaneous comparison of R AA and elliptic flow v 2 measurements at √ s NN = 5.02 TeV [66] with models can provide more stringent constraints to the implementation of the interaction and hadronisation processes for heavy quarks in the QGP. The comparison with models that compute both observables is shown in Fig. 10 for the R AA and v 2 , in the 0-10% and 30-50% centrality classes, respectively. The TAMU model overestimates R AA and underestimates v 2 at high p T , probably because it does not include radiative energy loss. The BAMPS el. model overestimates the maximum flow while underestimating the R AA value at high p T . The radiative energy loss contribution in BAMPS el.+rad. improves the description of R AA but gives v 2 values lower than the data. The LBT, PHSD, POWLANG and MC@sHQ models provide instead a fair description of v 2 . Nevertheless, energy loss is overestimated at high p T in the 0-10% centrality classes (but also in semi-central events) by PHSD, POWLANG and LBT, while at low p T the measured R AA is slightly higher than what predicted within LBT, PHSD and MC@sHQ. D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration  Figure 8: Average R AA of D 0 , D + and D * + mesons compared with model calculations. The three rows refer to the 0-10%, 30-50% and 60-80% centrality classes. The left panels show models based on heavy-quark transport, while the right panels show models based on pQCD energy loss. Model nomenclature and references: BAMPS [56], CUJET3.0 [64], Djordjevic [54], LBT [57], MC@sHQ+EPOS2 [62], PHSD [63] POWLANG [61], SCET [65], TAMU [58]. Some of the models are presented with two lines with the same style and colour, which encompass the model uncertainty band.

Summary
We have presented measurements of the p T -differential production yields of prompt D 0 , D + , D * + and D + s mesons at central rapidity in Pb-Pb collisions in the three centrality classes 0-10%, 30-50% and 60-80% at a centre-of-mass energy per nucleon pair √ s NN = 5.02 TeV.
The average R AA of the three non-strange D-meson species shows minimum values of 0.2 and 0.4 in the centrality classes 0-10% and 30-50%, respectively, at p T of 6-10 GeV/c. R AA increases for p T < 6 GeV/c, and it is compatible with unity at 1 < p T < 3 GeV/c. The average R AA values are compatible with those measured at √ s NN = 2.76 TeV and they have smaller uncertainties by a factor of about two, as well as extended p T coverage up to 50 GeV/c in central collisions. The similarity of the R AA values at the two energies was predicted by the Djordjevic model, and it results from the combination of a higher medium temperature at 5.02 TeV (estimated to be about 7% higher than at 2.76 TeV) with a harder p T distribution of charm quarks at 5.02 TeV. In central and semi-central collisions the average R AA of non-strange D mesons is compatible with that of charged particles for p T > 6 GeV/c, while it is larger at lower p T . The R AA of D + s mesons have generally larger central values than those of the average of non-strange D mesons, but the two measurements are compatible within about one standard deviation of the combined uncertainties.
The R AA of non-strange D mesons at high p T (above 10 GeV/c) is fairly described in the three centrality classes by model calculations that include both radiative and collisional energy loss. This indicates that the centrality dependence of radiative energy loss, which is the dominant contribution at high p T , is under good theoretical control. The R AA in the transverse momentum region below 10 GeV/c is described by several transport model calculations in central collisions, but most models fail in describing the centrality dependence of R AA and in describing simultaneously R AA and the elliptic flow coefficient v 2 . Therefore, the measurements provide significant constraints for the understanding of the interaction of charm quarks with the high-density QCD medium, especially at low and intermediate p T , where the R AA is the result of a more complex interplay among several effects.