Measurement of prompt $\rm{D_{s}^{+}}$-meson production and azimuthal anisotropy in Pb-Pb collisions at $\sqrt{s_{\rm NN}} = 5.02$ TeV

The production yield and angular anisotropy of prompt ${D_s^+}$ mesons were measured as a function of transverse momentum ($p_{ T}$) in Pb-Pb collisions at a centre-of-mass energy per nucleon pair $\sqrt{s_{ NN}} = 5.02$ TeV collected with the ALICE detector at the LHC. ${D_s^+}$ mesons and their charge conjugates were reconstructed at midrapidity ($|y|<0.5$) from their hadronic decay channel ${D_s^+ \to \phi \pi^+}$, with ${\phi \to K^-K^+}$, in the $p_{ T}$ intervals $2<p_{ T}<50$ GeV/$c$ and $2<p_{ T}<36$ GeV/$c$ for the 0-10% and 30-50% centrality intervals. For $p_{ T}>10$ GeV/$c$, the measured ${D_s^+}$-meson nuclear modification factor $R_{ AA}$ is consistent with the one of non-strange D mesons within uncertainties, while at lower $p_{ T}$ a hint for a ${D_s^+}$-meson $R_{ AA}$ larger than that of non-strange D mesons is seen. The enhanced production of ${D_s^+}$ relative to non-strange D mesons is also studied by comparing the $p_{ T}$-dependent ${D_s^+/D^0}$ production yield ratios in Pb-Pb and in pp collisions. The ratio measured in Pb-Pb collisions is found to be on average higher than that in pp collisions in the interval $2<p_{ T}<8$ GeV/$c$ with a significance of 2.3$\sigma$ and 2.4$\sigma$ for the 0-10% and 30-50% centrality intervals. The azimuthal anisotropy coefficient $v_2$ of prompt ${D_s^+}$ mesons was measured in Pb-Pb collisions in the 30-50% centrality interval and is found to be compatible with that of non-strange D mesons. The main features of the measured $R_{ AA}$, ${D_s^+/D^0}$ ratio, and $v_2$ as a function of $p_{ T}$ are described by theoretical calculations of charm-quark transport in a hydrodynamically expanding quark-gluon plasma including hadronisation via charm-quark recombination with light quarks from the medium. The $p_{ T}$-integrated production yield of ${D_s^+}$ mesons is compatible with the prediction of the statistical hadronisation model.


Abstract
The production yield and angular anisotropy of prompt D + s mesons were measured as a function of transverse momentum (p T ) in Pb-Pb collisions at a centre-of-mass energy per nucleon pair √ s NN = 5.02 TeV collected with the ALICE detector at the LHC. D + s mesons and their charge conjugates were reconstructed at midrapidity (|y| < 0.5) from their hadronic decay channel D + s → φπ + , with φ → K − K + , in the p T intervals 2 < p T < 50 GeV/c and 2 < p T < 36 GeV/c for the 0-10% and 30-50% centrality intervals. For p T > 10 GeV/c, the measured D + s -meson nuclear modification factor R AA is consistent with the one of non-strange D mesons within uncertainties, while at lower p T a hint for a D + s -meson R AA larger than that of non-strange D mesons is seen. The enhanced production of D + s relative to non-strange D mesons is also studied by comparing the p T -dependent D + s /D 0 production yield ratios in Pb-Pb and in pp collisions. The ratio measured in Pb-Pb collisions is found to be on average higher than that in pp collisions in the interval 2 < p T < 8 GeV/c with a significance of 2.3σ and 2.4σ for the 0-10% and 30-50% centrality intervals. The azimuthal anisotropy coefficient v 2 of prompt D + s mesons was measured in Pb-Pb collisions in the 30-50% centrality interval and is found to be compatible with that of non-strange D mesons. The main features of the measured R AA , D + s /D 0 ratio, and v 2 as a function of p T are described by theoretical calculations of charm-quark transport in a hydrodynamically expanding quark-gluon plasma including hadronisation via charmquark recombination with light quarks from the medium. The p T -integrated production yield of D + s mesons is compatible with the prediction of the statistical hadronisation model.

Introduction
Strongly-interacting matter at temperatures exceeding the pseudo-critical value of T pc ≈ 154-158 MeV and at vanishing baryon density is predicted to behave as a plasma of deconfined quarks and gluons (QGP) [1,2]. A QGP is formed and studied in ultrarelativistic heavy-ion collisions at the CERN Large Hadron Collider (LHC) and existing measurements indicate that it behaves as a strongly-coupled liquidlike system [3]. The lifetime of the QGP produced at the energy densities reached at the LHC is of the order of 10 fm/c [4]. Heavy quarks (charm and beauty) are sensitive probes to investigate the properties of the medium formed in these collisions. Due to their large masses, heavy quarks are produced predominantly in hard partonic scattering processes occurring during the early stages of the collision (i.e. on timescales shorter than the QGP formation time) and therefore experience the entire evolution of the medium. Heavy quarks propagate through the expanding hot and dense medium, interacting and exchanging energy and momentum with its constituents via both inelastic and elastic quantum chromodynamic (QCD) processes. At high momentum, the main effect of these interactions is the energy loss of the heavy quarks in the QGP due to medium-induced gluon radiation and collisional processes. On the other hand, low-momentum heavy quarks, including those shifted to low momentum by the energy loss, probe the diffusion regime dominated by elastic interactions. Since the charm and beauty quark masses are large compared to the medium temperature, the propagation of low-momentum heavy quarks through the fireball can be treated as a "Brownian motion", characterised by many elastic collisions with relatively small momentum transfers [5,6]. As a consequence of the large number of soft collisions with the medium constituents, heavy quarks can acquire significant collective flow when diffusing through the expanding fireball. Due to their large masses, charm quarks have a thermalisation time which is comparable to the fireball lifetime [5,7], and therefore they carry sensitive information on their coupling strength to the expanding medium, preserving memory of the thermalisation process. The process of hadronisation is also predicted to be modified in the presence of the QGP. Once the fireball approaches the pseudo-critical temperature for the transition to a hadron gas, a significant fraction of low-and intermediate-momentum heavy quarks could hadronise via recombination with other quarks from the medium [8][9][10][11], in competition with the fragmentation mechanism, which describes quark-tohadron transitions in pp, e ± p, and e + e − collisions [12,13].
The effects of the interaction of heavy quarks with the medium are commonly quantified by two main observables: the nuclear modification factor R AA and the elliptic flow v 2 . The R AA is defined as the ratio of the transverse-momentum (p T ) differential yields in nucleus-nucleus (AA) collisions and the cross section in proton-proton collisions, scaled by the average nuclear overlap function T AA R AA (p T ) = 1 T AA × dN AA /dp T dσ pp /dp T , where the yield in nucleus-nucleus collisions dN AA /dp T is measured in a given centrality interval and the T AA value is proportional to the average number of nucleon-nucleon collisions [14]. The T AA can be estimated via Glauber-model calculations tuned to match the measured multiplicity distribution of charged particles [15]. The elliptic flow v 2 is the second coefficient of the Fourier expansion of the particle-yield distribution in the azimuthal direction ϕ relative to the initial-state symmetry plane angle Ψ 2 : v 2 = cos[2(ϕ − Ψ 2 )] , where indicates the average over all particles and all events [16,17].
Measurements of non-strange D-meson production in heavy-ion collisions at RHIC [18] and LHC [19][20][21] energies show a substantial suppression of the D-meson yields compared to pp collisions at intermediate and high p T . In central nucleus-nucleus collisions, the R AA exhibits a pronounced drop for p T > 4-5 GeV/c, reaches a minimum around p T ≈ 8 GeV/c, and slightly increases at higher p T . This trend is described by different state-of-the-art model calculations of charm-quark energy loss in the QGP [22][23][24]. A positive D-meson v 2 is measured at p T > 8-10 GeV/c for semicentral Pb-Pb collisions at the LHC [25,26], and it is understood as originating from the path-length dependence of the Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration charm-quark energy loss in the geometrically anisotropic medium created in collisions with finite impact parameter. At lower p T , larger values of D-meson v 2 are observed accompanied by a bump-like structure in the R AA reflecting the radial flow of the fireball [19][20][21]. In particular, the D-meson v 2 for semicentral collisions shows a maximum value at p T ≈ 3 GeV/c, a clear mass ordering v 2 (D) < v 2 (p) < v 2 (π) at low p T (p T < 3 GeV/c), and a similar magnitude as the v 2 of charged pions at intermediate p T (3 < p T < 6 GeV/c) [25,26]. These features are consistent with a scenario in which low-momentum charm quarks acquire a significant collective flow when diffusing through the expanding QGP and hadronise via recombination with light quarks from the medium. The measured R AA and v 2 in this p T region are described qualitatively, and to some extent also quantitatively, by transport models including charmquark interactions in a hydrodynamically expanding QGP and hadronisation via both fragmentation and recombination [27][28][29][30][31][32][33][34][35][36][37][38]. However, a simultaneous description of the nuclear modification factor and the anisotropic flow of D mesons is still a challenge for theoretical models.
Studies of the production of different charm-hadron species, dubbed heavy-flavour hadrochemistry, can provide information on the hadronisation mechanism of charm quarks. In particular, an enhancement of the ground-state charm-strange meson yield relative to that of non-strange D mesons is expected in nucleus-nucleus collisions at low and intermediate momenta as compared to pp interactions, if the dominant process for D-meson formation is the recombination of charm quarks with light quarks from the medium, due to the large abundance of strange quarks in the QGP [39][40][41][42][43]. It was also pointed out in Ref. [43] that the comparison of the v 2 of D + s mesons to that of D mesons without strange-quark content (D 0 , D + , and D * + ) could provide sensitivity to the transport properties of the hadronic phase, since D + s mesons are expected to decouple early from the hadron gas and therefore do not pick up significant additional v 2 in the hadronic phase.
The production of D + s mesons was measured at RHIC [44] and the LHC [20,45] in Au-Au and Pb-Pb collisions at different centralities. So far, the results have shown that at low and intermediate p T the D + s /D 0 ratio in central, semicentral, and peripheral collisions is larger than the value measured in pp collisions, though the relatively large uncertainties do not allow firm conclusions. The magnitude and the p T dependence of the D + s /D 0 ratio are captured, at least qualitatively, by models including hadronisation via quark coalescence along with strangeness enhancement in the QGP [33,43,46,47], suggesting a relevant role of recombination processes in the hadronisation of low-momentum charm quarks in the QGP.
In this Letter, we report the measurements of the p T -differential yield and the nuclear modification factor of prompt D + s mesons in central (0-10%) and semicentral (30-50%) Pb-Pb collisions at √ s NN = 5.02 TeV, together with the measurement of the prompt D + s -meson elliptic flow in semicentral collisions. D + s mesons and their charge conjugates were reconstructed at midrapidity, |y| < 0.5, through their hadronic decay channel D + s → φπ + with a subsequent decay φ → K − K + . Prompt D + s mesons are defined as those produced directly in the hadronisation of charm quarks or originating from the decays of directly-produced excited open-charm and charmonium states, hence excluding weak decays of beauty hadrons. The data sample used for the analysis reported in this paper was collected with the ALICE detector at the end of 2018 and is larger by a factor of about 8 (4) for central (semicentral) collisions with respect to the sample collected in 2015, used for the previous publications of D + s -meson R AA and v 2 [20, 48].

Experimental apparatus and data sample
The ALICE apparatus comprises a central barrel, which is composed of a set of detectors for charged particle reconstruction and identification at midrapidity, a forward muon spectrometer, and various forward and backward detectors for triggering and event characterisation. A detailed description of the detectors and an overview of their typical performances can be found in Refs. [49,50].
Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration The D + s -meson decay candidates and charged conjugates were reconstructed and identified with the central barrel detectors, which cover the full azimuth in the pseudorapidity interval |η| < 0.9 and are embedded in a large solenoidal magnet providing a homogeneous magnetic field B = 0.5 T parallel to the beam direction. Charged-particle trajectories are reconstructed from their hits in the Inner Tracking System (ITS) and the Time Projection Chamber (TPC). The ITS is the innermost detector of the ALICE central barrel, it consists of six cylindrical layers of silicon detectors, allowing a precise determination of the track parameters in the vicinity of the interaction point. The TPC provides track reconstruction with up to 159 three-dimensional space points along the trajectory of a charged particle and provides particle identification via the measurement of the specific ionisation energy loss dE/dx. The Time-Of-Flight (TOF) detector, positioned at a radial distance of about 4 m from the beam axis, extends the particle-identification capabilities of the TPC by measuring the flight time of the charged particles from the interaction point to the TOF. The V0 detector is used for triggering and event selection, as well as for the estimation of the collision centrality and the reference plane for the elliptic flow measurement. It consists of two scintillator arrays, located on both sides of the nominal interaction point and covering the full azimuth in the pseudorapidity intervals −3.7 < η < −1.7 (V0C) and 2.8 < η < 5.1 (V0A). The neutron Zero Degree Calorimeters (ZDC), located along the beam axis on both sides of the central barrel at about 110 m distance from the interaction point, are used for event selection, along with the V0 detector.
The events used in the analysis were recorded with a minimum bias (MB) trigger which required coincident signals in the V0A and V0C detectors. Two additional trigger classes were used to enrich the sample of central and semicentral collisions via an online event selection based on the V0-signal amplitude. Background events due to the interaction of one of the beams with residual gas in the vacuum tube and other machine-induced backgrounds were rejected offline using the V0 and the ZDC timing information [50]. In order to have a uniform acceptance in pseudorapidity, only events with a primary vertex reconstructed within ±10 cm from the centre of the detector along the beam-line direction were considered in the analysis. Collisions were classified into centrality intervals, defined in terms of percentiles of the hadronic Pb-Pb cross section, based on the V0 signal amplitude, as described in detail in Ref. [51]. Central and semicentral collisions were considered in the analysis of the D + s -meson production. The sample of central collisions consists of about 100 × 10 6 events in the 0-10% centrality interval, corresponding to an integrated luminosity L int ≃ 130 µb −1 . For semicentral collisions, a sample of about 85 × 10 6 events in the 30-50% interval was utilised, corresponding to L int ≃ 56 µb −1 .
The average values of the nuclear overlap function, T AA , for the considered central and semicentral event intervals were estimated via Glauber-model simulations anchored to the measured charged-particle multiplicity distribution, and are 23.26 ± 0.17 mb −1 and 3.92 ± 0.06 mb −1 [52], respectively.
The Monte Carlo samples utilised in the analysis were obtained simulating Pb-Pb collisions with the HIJING 1.36 [53] event generator. In each simulated event, additional cc-and bb-quark pairs were injected using the PYTHIA 8.243 event generator [54,55] (Monash-13 tune [56]) and D + s mesons were forced to decay into the hadronic channel of interest for the analysis. The generated particles were propagated through the detector using the GEANT3 transport package [57]. The conditions of all the ALICE detectors in terms of active channels, gain, noise level, and alignment, and their evolution with time during the data taking period, were taken into account in the simulations.
3 Analysis technique D + s mesons and their charge conjugates were reconstructed via the decay channel D + s → φπ + → K − K + π + with branching ratio BR = (2.24 ± 0.08)% [58]. The analysis was based on the reconstruction of decayvertex topologies displaced from the interaction vertex. The separation induced by the weak decays of prompt D + s mesons is typically a few hundred of µm, cτ ≃ 151 µm [58].
Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration D + s -meson candidates were built combining triplets of tracks with the proper charge signs, each with |η| < 0.8, at least 70 out of 159 crossed TPC pad rows, a fit quality χ 2 /ndf < 1.25 in the TPC (where ndf is the number of degrees of freedom involved in the track fit procedure), and a minimum of two (out of six) hits in the ITS, with at least one in either of the two innermost layers, which provide the best pointing resolution. Moreover, at least 50 clusters available for particle identification (PID) in the TPC were required and only tracks with p T above 0.6 (0.4) GeV/c were considered for central (semicentral) collisions. These track-selection criteria limit the D + s -meson acceptance in rapidity, which drops steeply to zero for |y| > 0.5 at low p T and for |y| > 0.8 at p T > 5 GeV/c. Thus, only D + s -meson candidates within a p T -dependent fiducial acceptance region, |y| < y fid (p T ), were selected. The y fid (p T ) value was defined as a second-order polynomial function, increasing from 0.5 to 0.8 in the transverse-momentum range 0 < p T < 5 GeV/c, and as a constant term, y fid = 0.8, for p T > 5 GeV/c.
Unlike previous D-meson analyses based on linear selections [20, 21,26], a machine-learning (ML) approach based on Boosted Decision Trees (BDT) was adopted for the candidate selection to reduce the large combinatorial background [59]. In particular, the implementation of the BDT algorithm provided by the XGBoost [60] library was employed. Signal samples of prompt D + s mesons for the BDT training were obtained from Monte Carlo simulations as described in Section 2. The background samples were obtained from the sidebands of the candidate invariant-mass distributions in the data. Before the training, loose kinematic and topological selections were applied to the D + s -meson candidates together with the PID of decay-product tracks. Pions and kaons were selected by requiring compatibility with the respective particle hypothesis within three times the detector resolution between the measured and the expected signals for either the TPC dE/dx or the time of flight. Tracks without TOF hits were identified using only the TPC information. In addition, the absolute difference between the reconstructed K + K − invariant mass and the PDG average mass for the φ meson [58] (∆M KK ) was required to be below 15 MeV/c 2 . The candidate information provided to the BDTs, as an input for the models to distinguish among prompt D + s mesons and background candidates, was mainly based on the displacement of the tracks from the primary vertex, the distance between the D + s -meson decay vertex and the primary vertex, the D + s -meson impact parameter, and the cosine of the pointing angle between the D + s -meson candidate line of flight (the vector connecting the primary and secondary vertex) and its reconstructed momentum vector. The value of ∆M KK and additional variables related to the PID of decay tracks were also included. Independent BDTs were trained in the different p T intervals of the analysis and for the different centrality intervals. Subsequently, they were applied to the real data sample in which the belonging class, i.e., prompt D + s meson or combinatorial background, of particle candidates is unknown. Selections on the BDT output, which is related to the candidate's probability to be a prompt D + s meson, were optimised to reject a large fraction of the combinatorial background while maintaining high signal-selection efficiency.

Nuclear modification factor measurement
The raw yields of D + s mesons, including both particles and antiparticles, were extracted from binned maximum-likelihood fits to the invariant-mass distributions. The raw yields could be extracted in transverse-momentum intervals in the ranges 2 < p T < 50 GeV/c and 2 < p T < 36 GeV/c for the 0-10% and the 30-50% centrality intervals, respectively. The fit function was composed of a Gaussian for the description of the signal and of an exponential term for the background. An additional Gaussian was used to describe the peak due to the decay D + → K − K + π + , with a branching ratio of (9.68 ± 0.18) × 10 −3 [58], present at a lower invariant-mass value than the D + s -meson signal peak. The statistical significance of the observed signals S/ √ 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 standard deviations, varies from 4 to 24 depending on the p T and centrality intervals.
The p T -differential corrected yield of prompt D + s mesons was computed for each p T interval according The raw-yield values N D+D,raw , which contain the contribution of non-prompt D + s mesons from beautyhadron decays, were divided by a factor of two and multiplied by the prompt fraction f prompt to obtain the charge-averaged yields of prompt D + s mesons. Furthermore, they were divided by the acceptance times efficiency of prompt D + s mesons (Acc × ε) prompt , the BR of the decay channel, the width of the p T interval ∆p T , the correction factor for the rapidity coverage c ∆y , and the number of analysed events N evt . The correction factor for the rapidity acceptance c ∆y was computed with FONLL perturbative QCD calculations [61,62]. It was defined as the ratio between the generated D-meson yield in ∆y = 2 y fid and that in |y| < 0.5. The resulting values were in agreement within 1% with PYTHIA 8 simulations for pp collisions. To account for possible differences in Pb-Pb collisions and as an extreme variation, a flat rapidity distribution was also considered. The discrepancies with respect to FONLL calculations were negligible in comparison to other sources of systematic uncertainty described in Section 4.
The (Acc × ε) correction was obtained from the simulations described in Section 2 using samples not employed in the BDT training. The D + s -meson p T distributions from simulations were reweighted in order to use realistic momentum distributions in the determination of the (Acc × ε) factor, which depends on p T . In particular, the weights were defined to match the shape given by FONLL calculations multiplied by the R AA of D + s mesons predicted by the TAMU [33] model. The (Acc × ε) factors as a function of p T for prompt and non-prompt D + s mesons in the 0-10% and 30-50% centrality intervals are shown in Fig. 1. The difference between the (Acc × ε) factor for prompt and non-prompt D + s mesons arises from the BDT selections applied, given the different decay topology of D + s mesons coming from beauty-hadron decays. In particular, the non-prompt D + s mesons are on average more displaced from the primary vertex due to the large beauty-hadron lifetime, cτ ≃ 500 µm [58], and therefore are more efficiently selected in the low-p T region. At high p T , where the candidate decay length is less important to separate signal from background, the BDT selections are able to suppress the non-prompt efficiency with respect to the prompt one. The (Acc × ε) is higher for semicentral collisions, by up to a factor two at low p T , since less stringent selections can be applied thanks to the lower combinatorial background.
The f prompt fraction in each p T interval was obtained following the procedure employed in Refs. [20,21,63]. The calculation was based on the beauty-hadron production cross sections in pp collisions at √ s = 5.02 TeV from FONLL calculations, the beauty hadron to D + X decay kinematics from the PYTHIA 8 decayer, the (Acc × ε) correction factor for non-prompt D + s mesons, and the T AA for the corresponding centrality interval. In addition, the nuclear modification factor of D + s mesons from beautyhadron decays was accounted for and R prompt AA = R non-prompt AA was assumed as in Ref. [20]. The values of f prompt range between 0.86 and 0.91 depending on the p T interval and the centrality interval.
The prompt D + s -meson nuclear modification R AA factor was computed following Eq. 1. The measurement of the p T -differential cross section of prompt D + s mesons with |y| < 0.5 in pp collisions at √ s = 5.02 TeV from Ref. [64], which reaches up to p T = 24 GeV/c, was used as a reference for the R AA computation. At higher D + s -meson p T , 24 < p T < 50 GeV/c, FONLL calculations were used as a reference by scaling the predictions to match the measured values at lower p T . The p T -extrapolation procedure is the same as in Ref. [63]. As an example, the total systematic uncertainty of the pp reference in the 36 < p T < 50 GeV/c interval is +42 −33 %.

Elliptic flow measurement
The elliptic flow of prompt D + s mesons was measured for semicentral events in transverse-momentum intervals in the range 2 < p T < 24 GeV/c. The same ML models trained for the R AA measurement in the 30-50% centrality interval were used and the same selections on the BDT output were applied. The  Figure 1: Acceptance-times-efficiency factor for D + s mesons as a function of p T . The (Acc × ε) factors for prompt (red) and non-prompt (blue) D + s mesons in Pb-Pb collisions for the 0-10% centrality interval are shown, together with those for prompt (orange) and non-prompt (green) D + s mesons for the 30-50% centrality interval.
analysis procedure for the v 2 determination followed closely with what was done in Ref. [26] for the measurement of the non-strange D-meson elliptic flow. The D + s -meson v 2 coefficients were measured using the Scalar Product (SP) method [65,66] and can be expressed as where u 2 = e i2ϕ D is the unit flow vector of the D-meson candidate with azimuthal angle ϕ D , Q Q Q k 2 is the subevent 2 nd -harmonic flow vector for the subevent k, and M k represents the subevent multiplicity. The SP denominator was calculated with the formula introduced in Ref. [66], where the three subevents, indicated as A, B, and C, are defined by the particles measured in the V0C, V0A, and TPC detectors, respectively. For the TPC detector, the Q Q Q 2 vector was computed from the azimuthal angles of charged tracks reconstructed with |η| < 0.8 and M was the number of measured tracks. For the V0A and V0C detectors, the Q Q Q 2 vectors were calculated from the azimuthal distribution of the energy deposition in the detector sectors and M was the sum of the amplitudes measured in each channel [26]. The Q Q Q 2 vectors were corrected for detector effects arising from the non-uniform acceptance [67]. The single bracket in Eq. 3 refers to an average over all the events, while the double brackets denote the average over all particles in the considered p T interval and all events. The SP denominator was obtained as a function of the collision centrality.
The elliptic flow of D + s mesons cannot be directly measured using Eq. 3 as signal candidates cannot be identified on a particle-by-particle basis. The measured anisotropic flow coefficient v tot 2 can be written as a weighted sum of the v 2 of candidates reconstructed from true D + s -meson decays (v sig 2 ) and that of the where N sig and N bkg are the raw signal and background yields, respectively. An additional v D + 2 free parameter and the corresponding raw yield N D + were included to account for the D + → K − K + π + Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration  Figure 2: Simultaneous fit to the invariant-mass spectrum and v 2 (M D ) of D + s -meson candidates in the 4 < p T < 6 GeV/c interval for the 30-50% centrality interval. The solid blue and the dotted red curves represent the total and combinatorial-background fit functions, respectively. contribution to the measured v tot 2 distribution. A simultaneous fit to the invariant-mass spectrum and the v tot 2 distribution as a function of the invariant mass was performed in each p T interval to extract the elliptic flow coefficients. The fit function for the invariant-mass distributions was composed of two Gaussian terms to describe the signal and the peak due to the decay D + → K − K + π + , and an exponential distribution for the background, as for the R AA measurement of Section 3.1. The v sig 2 was measured from a fit to the v tot 2 distribution with the function of Eq. 4, where v bkg 2 (M D ) was described by a linear function. Figure 2 shows the simultaneous fit to the invariant-mass spectrum and v tot 2 (M D ) of D + s mesons in the 4 < p T < 6 GeV/c interval for the 30-50% centrality interval.
The reconstructed D + s -meson signal is a mixture of prompt D + s mesons and non-prompt D + s mesons from beauty-hadron decays. Therefore, the v = v prompt 2 /2. This hypothesis is based on the v 2 measurements of the non-prompt J/ψ performed by ATLAS and CMS [69,70], and on the available model calculations [71][72][73] 4 Systematic uncertainties

Nuclear modification factor measurement
The measurement of the D + s -meson corrected yield is affected by the following sources of systematic uncertainties: (i) the raw-yield extraction from the invariant-mass distributions, (ii) the track-reconstruction efficiency, (iii) the PID and selection efficiency, (iv) the generated D + s -meson p T shape in the simulation, and (v) the prompt fraction estimation. In addition, the uncertainty due to the branching ratio of Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration  [58], and that due to the centrality-interval definition were considered. This last contribution arises from the uncertainty of the fraction of the hadronic cross section used in the Glauber fit to determine the centrality, and was estimated to be < 0.1% and 2% for the 0-10% and 30-50% centrality intervals, respectively [63]. A procedure similar to that described in Refs. [20,21] was used to estimate the uncertainties as a function of the p T interval and the centrality interval. The estimated values of the systematic uncertainties are summarised in Table 1 for representative p T intervals, together with the total systematic uncertainty obtained from the sum in quadrature of the different contributions.
The systematic uncertainty of the raw-yield extraction was evaluated by repeating the fit of the invariantmass distribution varying the lower and upper limits of the fit range, the bin width, and the functional form of the background fit function. The systematic uncertainty was defined as the RMS of the distribution of the signal yields obtained from all these variations and ranges from 2% to 8% depending on the centrality interval and the p T interval.
The systematic uncertainty of the track-reconstruction efficiency was estimated by varying the trackquality selection criteria and by comparing the prolongation probability of the TPC tracks to the ITS hits in data and simulation. The comparison was performed after weighting the relative abundances of primary and secondary particles in the simulation to match those observed in data [74]. The estimated uncertainty ranges from 5% to 14%.
The systematic uncertainty of the selection efficiency originates from imperfections in the description of the detector resolutions and alignments in the simulation. It was estimated by comparing the corrected yields obtained by repeating the analysis with different selections on the BDT output, which resulted in up to 50% higher and lower efficiencies with respect to the central values. The assigned systematic uncertainty ranges from 3% to 9%. Possible systematic effects due to the loose PID selection, applied prior to the machine-learning one, were investigated comparing pion and kaon PID selection efficiencies in data and in simulations. A pure sample of pions was selected from K 0 S and Λ decays, while samples of kaons in the TPC (TOF) were obtained applying a strict PID selection using the TOF (TPC) information. Since no significant differences were observed, no systematic uncertainty was assigned.
An additional contribution to the systematic uncertainty of the efficiency originates from possible differences between the real and simulated D + s -meson p T distributions. It was estimated by calculating Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration Table 2: Systematic uncertainties of the prompt D + s -meson v 2 in Pb-Pb collisions for the 30-50% centrality interval in representative p T intervals. The uncertainties of the fitting procedure and non-prompt contribution subtraction are quoted as absolute uncertainties, while that of the SP denominator as relative uncertainty. the efficiency using alternative D + s -meson p T shapes obtained by re-weighting the p T distribution from MC simulations to match those predicted by theoretical models. The p T distributions from FONLL calculations including or not hot-medium effects, parametrised using the p T -differential R AA from the LGR [34], PHSD [75], TAMU [33], and Catania [35] models, were considered. The resulting uncertainty was estimated to be about 5% and 3% for the 0-10% and 30-50% centrality intervals, respectively, in the lowest p T intervals where the efficiency varies steeply with p T , and to decrease to zero above 12 GeV/c.
The systematic uncertainty of the prompt fraction was estimated by varying the FONLL parameters (bquark mass, factorisation, and renormalisation scales, according to the prescription reported in Ref. [76]) in the calculation of the p T -differential production cross section of non-prompt D + s mesons. In addition, the ratio of the non-prompt and prompt D + s -meson R AA was varied in the range 1 3 < R non-prompt AA /R prompt AA < 3 as done in Ref. [20]. The resulting uncertainty ranges between +8 −16 % and +12 −23 %. In the R AA calculation, the BR uncertainty of the D + s -meson yield in Pb-Pb collisions and of the pp reference cross section cancels out in the ratio. The contribution due to the prompt fraction uncertainty, estimated by the variation of the parameters of the FONLL calculation, was considered to be fully correlated and the remaining systematic uncertainties were propagated as uncorrelated. The uncertainties of the R AA normalisation are the quadratic sum of the pp normalisation uncertainty, 2.1% [77], the T AA uncertainty, 0.7% (1.5%) for the 0-10% (30-50%) centrality interval [52], and the one related to the centrality-interval definition described above.

Elliptic flow measurement
The systematic uncertainties of the measurement of the D + s -meson v 2 coefficients were estimated with procedures similar to those detailed in Ref. [26]. They include the following sources: (i) the signal extraction from the invariant-mass and v tot 2 distributions, (ii) the non-prompt D + s contribution, and (iii) the centrality dependence of the SP denominator. The selection efficiency was observed to be independent of the D + s -meson azimuthal direction, therefore no contribution to the systematic uncertainty was assigned. The non-flow effects are naturally suppressed due to the pseudorapidity gap of at least 0.9 units between the pseudorapidity interval used for the D + s -meson reconstruction, and the V0C used for the Q Q Q 2 -vector determination. The estimated values of the systematic uncertainties are summarised in Table 2 for representative p T intervals.
The uncertainty due to the simultaneous fit was estimated by repeating the fit several times with different configurations, as done for the R AA measurement. The RMS of the v 2 distribution obtained from the different trials, separately for each p T interval, was assigned as systematic uncertainty. The absolute systematic uncertainty values due to the signal extraction range between 0.01 and 0.03 depending on p T .
The systematic uncertainty related to the correction for the contribution of non-prompt D + s to the measured v 2 has two main sources. The first one is due to the f prompt calculation and it was treated as described in Section 4.1 for the R AA measurement. The second source is due to the assumption  −0.005 and +0.039 −0.009 for the different p T intervals. The contribution of the SP denominator to the systematic uncertainty is due to the centrality dependence. The uncertainty was evaluated as the difference of the centrality integrated value, computed from the events in the 30-50% interval, with that obtained as weighted average of SP-denominator values in narrow centrality intervals using the D + s -meson yields as weights. A systematic uncertainty of 0.5% was assigned.

Results
The p T -differential production yields dN/dp T of prompt D + s mesons measured in the 0-10% and 30-50% centrality intervals are shown in Fig. 3. For the semicentral class of events, the measurements are scaled by 10 −1 for better visibility. The results are compared with the pp reference cross section multiplied by the corresponding average nuclear overlap function T AA . The larger data sample and the improved analysis technique enable an extended p T coverage and finer p T intervals in the measured dN/dp T of prompt D + s mesons compared to the previous measurement by the ALICE Collaboration in Pb-Pb collisions at √ s NN = 5.02 TeV [20]. A strong suppression of the D + s yields compared to the binary-scaled pp reference is observed for both centrality intervals for p T > 3-4 GeV/c, similarly as for the non-strange D mesons [21]. This suppression is understood in terms of modification of the charmquark momentum spectra due to the interactions within the QGP. D 0 , D + , and D * + mesons in Fig. 4 for the 0-10% and 30-50% centrality intervals, in the left and right panels, respectively. The systematic uncertainties related to the tracking efficiency and the promptfraction estimation are considered as fully correlated between the different D-meson species, and are reported separately from the other sources of systematic uncertainty which are uncorrelated. The R AA of D + s and non-strange D mesons show a minimum value of about 0.2 (0.4) around p T ≈ 10 GeV/c in the 0-10% (30-50%) centrality interval. For lower p T , the R AA increases with decreasing p T reaching about unity around p T ≈ 2-3 GeV/c. In both centrality intervals, the R AA of prompt D + s and non-strange D mesons are compatible within uncertainties for p T 10 GeV/c. In this p T region, the hadronisation is expected to occur mainly via fragmentation and the dominant effect leading to the observed suppression is the charm-quark energy loss in the QGP. For lower p T , the measured R AA of prompt D + s mesons is systematically higher than that of non-strange D mesons but compatible within about one standard deviation of the combined statistical and systematic uncertainties.
In the left and right panels of Fig. 5, the R AA of prompt D + s and non-strange D mesons in the 0-10% centrality interval are compared with theoretical calculations implementing charm-quark transport in the QGP [78]. All the models include an enhancement of the strangeness content of the QGP and the hadronisation of charm quarks is implemented either via fragmentation, which is dominant at high p T , or via coalescence with light quarks in the QGP. In the Catania [35,47] and LGR [34] models the coalescence occurs instantaneously at the phase boundary and is implemented through the Wigner formalism [79]. In the PHSD model [38,75], the hadronisation in heavy-ion collisions is described via a Monte Carlo simulation of the coalescence process in competition to fragmentation. In the TAMU [33] model, the hadronisation via coalescence proceeds via formation of resonant states when approaching the (pseudo)critical temperature within the formalism of a Resonance Recombination Model [11]. For the description of the D-meson p T spectra in pp collisions, all the models use as starting point FONLL calculations [61,62,76]. Charm quarks are hadronised in pp collisions with fragmentation in the PHSD and LGR models, while in the Catania model the charm-quark hadronisation via coalescence is also implemented in addition to that via fragmentation [80]. In pp collisions, the hadronisation in the TAMU model is instead determined with a statistical hadronisation approach, in which the strangeness production is suppressed in pp with respect to heavy-ion collisions. This is described with a suppression factor for strange particles of γ s = 0.6 [81], which is instead unity in heavy-ion collisions. All the models TeV compared with theoretical calculations based on charm-quark transport in a hydrodynamically expanding QGP implementing strangeness enhancement and hadronisation of charm quarks via coalescence in addition to fragmentation in the vacuum [33-35, 38, 75]. The boxes represent the total systematic uncertainties. The colour bands represent the theoretical uncertainty when available.
reproduce qualitatively the measured R AA of prompt D + s and non-strange D mesons. The Catania model underestimates both measurements for 2 < p T < 5 GeV/c by about 2 σ of the combined statistical and systematic uncertainties of the measured points, while it overestimates the non-strange D-meson R AA for p T < 1.5 GeV/c, where no measurement is available for strange mesons. In contrast, the PHSD model describes well the measured nuclear modification factors for p T < 5 GeV/c and underestimates them by about 2 σ for higher p T . The TAMU model describes the measurements within uncertainties, with a tension of about 2 σ of the combined statistical and systematic uncertainties of the D + s -meson measurement in 2 < p T < 3 GeV/c. These three models do not include charm-quark interactions with medium constituents via radiative processes, hence are not expected to describe the R AA of strange and non-strange D mesons for p T > 6-8 GeV/c. The LGR model, which instead includes gluon-radiation processes, provides a good description of the R AA up to high p T . All the models predict a smaller suppression of the D + s -meson R AA compared to non-strange D mesons at low and intermediate p T . The possible enhancement of the yield of D mesons with strange-quark content with respect to that of non-strange D mesons was further investigated by computing the ratio between the p T -differential production yields of prompt D + s mesons and those of prompt D 0 mesons [21]. The systematic uncertainty related to the determination of the tracking efficiency and the contribution due to the subtraction of the component from beauty-hadron decays were propagated as fully correlated in the ratios, while all the other sources of systematic uncertainties were considered as uncorrelated between the measurements of D + s and D 0 mesons. The top row of Fig. 6 shows the D + s /D 0 yield ratios in the 0-10% (left panel) and 30-50% (middle panel) centrality intervals compared to the same quantity measured in minimumbias pp collisions [64]  Pb-Pb collisions are described within uncertainties by the Catania (PHSD) model. The TAMU model significantly overestimates the measured D + s /D 0 by a similar amount in the two colliding systems, leading to a good description of the ratio of the D + s /D 0 measured in Pb-Pb and pp collisions, as shown in the bottom panels of Fig. 6. While the Catania and PHSD models predict a D + s /D 0 ratio almost p T independent for p T < 3 GeV/c and then mildly decreasing towards the pp value at higher p T , the TAMU and LGR models predict a peak around p T ≈ 3-4 GeV/c. The origin of such a peak would be motivated by the different masses of D + s and D 0 mesons and by the collective radial expansion of the system with a common flow-velocity profile, which imposes an equal velocity boost to all particles in case of complete thermalisation. In addition, also the hadronisation via coalescence is expected to modify the p T shape of the D + s /D 0 ratio due to the different masses of u and s quarks. A similar p T shape is predicted by the GSI-Heidelberg statistical hadronisation model (SHMc) [82], which is reported in the top panels of Fig. 6 for central and semicentral Pb-Pb collisions, where the p T spectra of charm hadrons are modelled with a core-corona approach. The low-p T region is dominated by the core contribution described with a Blast Wave function. The corona contribution is instead parametrised from measurements in pp collisions and is relevant at high p T . The p T -spectra modification due to resonance decays is computed using the FastReso package [83]. Within the current uncertainties of the measurement, no firm conclusions can be drawn on the p T shape of the D + s /D 0 ratio in Pb-Pb collisions at low and intermediate p T . These results however provide important indications about the role of the charm-quark hadronisation via coalescence in the QGP, complementary to those obtained via the simultaneous comparison of the measured D-meson R AA and v n coefficients [21,26].
The visible production yield of prompt D + s mesons was evaluated by integrating the p T -differential yield over the narrower p T intervals of the measurement. The systematic uncertainties were propagated as Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration fully correlated among the measured p T intervals, except for the raw-yield extraction uncertainty, which was treated as uncorrelated considering the variations of the signal-to-background ratio and the shape of the combinatorial-background distribution as a function of p T . In order to obtain the p T -integrated production yield, the dN/dp T was extrapolated in 0 < p T < 2 GeV/c. For this purpose, the measured p Tdifferential D + s /D 0 ratio was interpolated using the shape predicted by the PHSD model and leaving the normalisation as a free parameter. The extrapolated D + s /D 0 ratio for p T < 2 GeV/c was then multiplied by the dN/dp T of D 0 mesons measured in the same p T interval [21] to obtain the extrapolated D + s yield, which amounts to about 70% of the total production yield. An additional uncertainty was assigned to the extrapolation procedure, by repeating the computation using the TAMU and Catania transport models, and the SHMc to extrapolate the D + s /D 0 ratio in the unmeasured p T interval. Finally, the p T -integrated production yield was obtained as the sum of the extrapolated one for p T < 2 GeV/c and the measured one. The results for the 0-10% and 30-50% centrality intervals are reported in Table 3. As for the D 0 , D + , and D * + mesons [21], the production yield of prompt D + s mesons at midrapidity is compatible within uncertainties with the one predicted by the SHMc. This suggests that low-p T charm quarks, which determine the total yield, are thermalised in the QGP.
The degree of thermalisation of charm quarks and their hadronisation in the QGP were also studied via the measurement of the azimuthal anisotropy in the prompt D + s -meson production. Figure 7 shows the elliptic flow coefficient v 2 of prompt D + s mesons for the 30-50% centrality interval measured in the transverse-momentum interval 2 < p T < 24 GeV/c, compared with that of prompt non-strange D mesons (left panel) and with theoretical calculations (right panel). The rapidity interval of the measurement, Prompt D + s mesons in Pb-Pb at √ s NN = 5.02 TeV ALICE Collaboration |y| < 0.8, is wider than that quoted for the R AA since no correction for the rapidity acceptance was applied. The measurement was carried out in finer p T intervals and has uncertainties reduced by a factor up to four with respect to the previous measurement [48], thanks to the more advanced D + s -meson selection technique and the larger data sample. Considering as null hypothesis v 2 = 0, the probability to observe the measured positive v 2 in 2 < p T < 8 GeV/c corresponds to a significance of 6.4σ , confirming the participation of the charm quark in the collective motion of the system, as already observed for non-strange D mesons [26,48]. However, within the current uncertainties it is not possible to draw a conclusion about a potential difference between the elliptic flow of strange and non-strange D mesons, which would be motivated by the different mass, the charm-quark hadronisation via recombination with strange quarks in the medium instead of light quarks [84], and possible differences in the hadronic phase [43]. The measured D + s -meson v 2 is compatible within uncertainties with the predictions of the TAMU and PHSD models, which include charm-quark coalescence with flowing strange quarks in the medium.

Summary
In this Letter, a comprehensive and high-precision set of measurements regarding the prompt D + s -meson production at midrapidity in Pb-Pb collisions at √ s NN = 5.02 TeV was reported.
The p T -differential production yields were measured in a wide transverse-momentum interval between 2-50 (2-36) GeV/c in the 0-10% (30-50%) centrality interval. They were used to compute the p Tdifferential R AA and the ratio of D + s -meson production relative to D 0 mesons. The measured R AA shows a strong suppression of the D + s -meson production yield compared to the binary-scaled pp reference, reaching a minimum of about 0.2 (0.4) around p T ≈ 10 GeV/c in the 0-10% (30-50%) centrality interval. For lower p T , the R AA increases reaching about unity for p T ≈ 2-3 GeV/c. The D + s /D 0 yield ratios in Pb-Pb collisions are higher than those measured in pp collisions for p T 8 GeV/c with a significance of 2.3σ and 2.4σ in the 0-10% and 30-50% centrality intervals, respectively. This finding is consistent with the predictions of theoretical calculations implementing the charm-quark transport in a hydrodynamically expanding QGP, which include an enhanced strange-quark production in the medium and the charmquark hadronisation via coalescence. The production yield of prompt D + s mesons, extrapolated down to p T = 0, in the 0-10% centrality interval is compatible with the prediction of the SHMc, suggesting that the bulk of charm quarks are thermalised in the QGP.
The elliptic flow coefficient v 2 of prompt D + s mesons was measured as a function of p T in the 30-50% centrality interval. The D + s -meson v 2 in 2 < p T < 8 GeV/c is positive with a significance of 6.4σ and is compatible within uncertainties with that of non-strange D mesons. The measured v 2 is also described by several transport-model calculations implementing the charm-quark hadronisation via coalescence.
The data reported in this Letter represent the most precise measurements of prompt D + s -meson production in heavy-ion collisions at LHC energies to date, and provide stringent constraints to all models on the production of charm quarks and their hadronisation in the QGP. Future data samples that will be collected with the upgraded ALICE detector in Run 3 will have the potential to further improve and extend to lower p T the measurement of D + s mesons in heavy-ion collisions [85].      [13] ALICE Collaboration, S. Acharya et al., "Charm-quark fragmentation fractions and production cross section at midrapidity in pp collisions at the LHC", arXiv:2105.06335 [nucl-ex].