Transverse-momentum and event-shape dependence of D-meson flow harmonics in Pb-Pb collisions at $\sqrt{s_{NN}} = 5.02$ TeV

The elliptic and triangular flow coefficients $v_2$ and $v_3$ of prompt D$^{0}$, D$^{+}$, and D$^{*+}$ mesons were measured at midrapidity ($|y|<0.8$) in Pb-Pb collisions at the centre-of-mass energy per nucleon pair of $\sqrt{s_{NN}} = 5.02$ TeV with the ALICE detector at the LHC. The D mesons were reconstructed via their hadronic decays in the transverse momentum interval $1<p_{\rm T}<36$ GeV/$c$ in central (0-10%) and semi-central (30-50%) collisions. Compared to pions, protons, and J/$\psi$ mesons, the average D-meson $v_{n}$ harmonics are compatible within uncertainties with a mass hierarchy for $p_{\rm T} \lesssim 3$ GeV/$c$, and are similar to those of charged pions for higher $p_{\rm T}$. The coupling of the charm quark to the light quarks in the underlying medium is further investigated with the application of the event-shape engineering (ESE) technique to the D-meson $v_2$ and $p_{\rm T}$-differential yields. The D-meson $v_2$ is correlated with average bulk elliptic flow in both central and semi-central collisions. Within the current precision, the ratios of per-event D-meson yields in the ESE-selected and unbiased samples are found to be compatible with unity. All the measurements are found to be reasonably well described by theoretical calculations including the effects of charm-quark transport and the recombination of charm quarks with light quarks in a hydrodynamically expanding medium.


Introduction
The formation of a strongly coupled colour-deconfined medium in ultra-relativistic heavy-ion collisions, called quark-gluon plasma (QGP), has been established both at RHIC and LHC energies [1,2]. The QGP behaves as a near-perfect fluid with small shear viscosity over entropy density ratio, η/s, undergoing an expansion that can be described by relativistic hydrodynamics [3].
In heavy-ion collisions, heavy quarks (charm and beauty) are predominantly produced via hard-scattering processes on a time scale shorter than the QGP formation time [4,5], and therefore they experience all the stages of the system evolution, interacting with the medium constituents via both elastic (collisional) [6] and inelastic (gluon radiation) [7][8][9] processes. The measurement of the suppression of the yield of heavy-flavour hadrons in central nucleus-nucleus collisions relative to pp collisions scaled by the number of nucleon-nucleon collisions at both RHIC [10][11][12][13][14] and LHC energies [15][16][17][18][19][20][21] provides compelling evidence of heavy-quark energy loss in deconfined strongly interacting matter.
Additional insights into the QGP properties can be obtained by measuring the azimuthal anisotropy of heavy-flavour hadrons. In non-central nucleus-nucleus collisions the initial spatial anisotropy of the overlap region is converted via multiple interactions into an azimuthally anisotropic distribution in the momentum space of the produced particles [22,23]. This anisotropy is characterised in terms of the Fourier coefficients v n = cos[n(ϕ − Ψ n )] , where ϕ is the azimuthal angle of the particle and Ψ n is the azimuthal angle of the symmetry plane for the nth-order harmonic [23,24]. The values of the Fourier coefficients depend on the geometry of the collision, the fluctuations in the distributions of nucleons and gluons within the nuclei [25], and the dynamics of the expansion. The second order flow coefficient v 2 , called elliptic flow, is related to the almond-shaped geometry of the overlap region between the colliding nuclei and, consequently, is the largest contribution to the anisotropy in non-central collisions. The third harmonic coefficient v 3 , named triangular flow, originates from event-by-event fluctuations in the initial distribution of nucleons and gluons in the overlap region [26]. In particular, the measurement of the azimuthal anisotropy of heavy-flavour hadrons at low p T can help quantify the extent to which charm and beauty quarks participate in the collective expansion of the medium [27], as well as the fraction of heavy-flavour hadrons hadronising via recombination with flowing light quarks [28,29]. At high p T , instead, the charm hadron azimuthal anisotropy can constrain the path-length dependence of heavy-quark in-medium energy loss [30,31]. Precise measurements of heavy-flavour v n coefficients are useful to constrain the parameters of models that implement the heavy-quark transport in the QGP. In this context, the heavy-quark spatial diffusion coefficient D s in the QGP is particularly interesting, since it is related to the relaxation (equilibration) time of heavy quarks τ Q = (m Q /T )D s , where m Q is the quark mass and T is the medium temperature [32].
Further investigation into the dynamics of heavy quarks in the medium can be performed with the event-shape engineering (ESE) technique [33], which allows for selection of events with the same centrality but different initial geometry on the basis of the magnitude of the average bulk flow. In fact, hydrodynamic calculations show that the average flow of the bulk of soft hadrons is proportional to the initial-state eccentricities [34] for small values of η/s [3,35,36]. By classifying the events with the ESE technique it is possible to investigate the correlation between the flow coefficients of D mesons and soft hadrons. According to the available calculations [34,37,38], the initial system ellipticity is converted into parton flow with a similar efficiency for bulk and charm quarks, despite the different production mechanisms, dynamics, and hadronisation of heavy quarks and light partons forming the bulk of the medium. Moreover, the measurement of the D-meson spectra in events with different average eccentricity provides information about the possible correlation between the radial and elliptic flows at low-intermediate p T , and the charm-quark energy loss and the elliptic flow at high p T . The correlation with the radial flow is expected to be present from the observation of the scaling of the flow harmonics with the particle mass [39], while the correlation with the in-medium energy loss would be motivated by the different path traversed by the charm quark in the medium in the case of an isotropic or an eccentric system.
A positive D-meson v 2 is observed at both RHIC [10,40,41] and the LHC [42][43][44][45][46][47][48]. The comparison of the D-meson v 2 with the charged-pion v 2 and with theoretical models [49][50][51][52][53][54][55][56][57] indicates that charm quarks participate in the collective expansion of the medium and that both collisional processes and the recombination of charm and light quarks contribute to the observed elliptic flow. Furthermore, a positive D 0 -meson v 3 was measured by the CMS Collaboration [47]. The p T -differential yields and v 2 of D mesons were measured by the ALICE Collaboration in samples of events selected on the basis of the average bulk elliptic flow with the ESE technique [48]. A correlation between the D-meson v 2 and the v 2 of the bulk of light hadrons was observed, while the ratio of the p T -differential yields in ESE-selected samples to the yields measured without any ESE selection was found to be compatible with unity within the large uncertainties.
In this Letter, the measurement of the non-strange D-meson flow harmonics performed on a large sample of Pb-Pb collisions at √ s NN = 5.02 TeV collected by ALICE in 2018 is reported. With this data sample, the D-meson v 2 is measured with the Scalar Product (SP) method in an extended p T interval and with smaller uncertainties with respect to the previous results obtained with the Event Plane (EP) method described in [46,48] in the 30-50% centrality class. The results obtained with the SP and the EP method were found to be compatible between each other, as reported in previous publications [43]. The measurement of the D-meson v 2 coefficient in the 0-10% centrality class and v 3 coefficient in the 0-10% and 30-50% centrality classes are also presented. In addition, the measurement of the v 2 and the modification of the p T distributions in the ESE-selected samples is reported in narrower classes of the average event flow with respect to [48]. The measurements are compared to theoretical calculations in order to assess information about the participation of the charm quark in the collective motion of the system and its interactions with the QGP constituents.

Detector and data sample
A detailed description of the ALICE apparatus and data acquisition framework can be found in [58,59]. The main detectors used for this analysis are the Inner Tracking System (ITS) [60], the Time Projection Chamber (TPC) [61], and the Time-Of-Flight (TOF) detector [62]. The ITS is a six-layer silicon detector which provides the event selection, the reconstruction of primary and secondary vertices, and the tracking of charged particles. The TPC detector is used for the track reconstruction and the particle identification (PID) via the measurement of the specific energy loss dE/dx, while the TOF detector provides PID via the measurement of the flight time of the particles. These detectors are located inside a solenoid providing a uniform magnetic field of 0.5 T parallel to the LHC beam direction and cover the pseudorapidity interval |η| < 0.9. A minimum-bias interaction trigger was provided by the coincidence of signals in the two scintillator arrays of the V0 detector [63], covering the full azimuth in the pseudorapidity regions −3.7 < η < −1.7 (V0C) and 2.8 < η < 5.1 (V0A). An online selection based on the V0 signal amplitudes was applied in order to enhance the sample of central and mid-central collisions through two separate trigger classes. Background events from beam-gas interactions were removed offline using the time information provided by the V0 and the neutron Zero-Degree Calorimeters (ZDC) [64]. Only events with a primary vertex reconstructed within ±10 cm from the centre of the detector along the beam line were considered in the analysis.
Events were divided into centrality classes, defined in terms of percentiles of the hadronic Pb-Pb cross section, using the amplitudes of the signals in the V0 arrays. The number of events in each centrality class considered for this analysis (0-10% and 30-50%) is about 100 × 10 6 and 85 × 10 6 , corresponding to an integrated luminosity of ≃ 130 µb −1 and ≃ 56 µb −1 , respectively [65]. In order to apply the ESE technique, the events in each centrality class were further divided into samples with different average elliptic anisotropy of final-state particles, selected according to the magnitude of the second-order harmonic reduced flow vector q 2 [36], defined as where M is the number of tracks used in the |Q Q Q 2 | calculation selected as described below, and is a vector built from the azimuthal angles (ϕ k ) of the considered particles. The Q Q Q 2 vector was measured using charged tracks reconstructed in the TPC with |η| < 0.8 and 0.2 < p T < 5 GeV/c to exploit the ϕ resolution of the TPC and the large multiplicity at midrapidity, which are crucial to maximise the selectivity of q 2 with respect to the final state flow eccentricity [48,66]. The denominator in Eq. 1 is introduced to remove the dependence of |Q Q Q 2 | on √ M in the absence of flow [36]. The tracks used to form the D-meson candidates were excluded from the computation of q 2 to partially remove autocorrelations between D mesons and q 2 . The effect of residual autocorrelations between the D mesons and q 2 was studied in [48] by computing q 2 from the azimuthal distribution of the energy deposition measured in the V0A detector, and hence introducing a pseudorapidity gap of two units between the D mesons and q 2 . The v 2 values obtained with the q 2 calculated with TPC tracks and using the V0 detector were found to be compatible with a reduction of the eccentricity discriminating power of the two detectors without allowing for a firm conclusion on the magnitude of non-flow contamination. The same study was repeated for the data sample used for this analysis, leading to the same conclusions.
The selection of the events according to the average elliptic flow of the event was performed by defining q 2 percentiles in 1%-wide centrality intervals as described in [48] and [67] to avoid the introduction of biases in the centrality (multiplicity) distribution of the selected events. The ESE-selected classes defined for the analyses presented in this paper correspond to the 20% of events with smallest and largest q 2 , respectively, and will be indicated as "small-q 2 " and "large-q 2 ". In case of no ESE selection, the term "unbiased" will be used.

Analysis technique
The charmed mesons were reconstructed at midrapidity via the decay channels D 0 → K − π + (with branching ratio, BR = 3.89 ± 0.04%), D + → K − π + π + (BR = 8.98 ± 0.28%), and D * + → D 0 π + (BR = 67.7 ± 0.5%) and their charge conjugates [68]. D 0 and D + candidates were built combining pairs and triplets of tracks with the proper charge, p T > 0.4 GeV/c, |η| < 0.8, a fit quality χ 2 /ndf < 2 in the TPC (where ndf is the number of degrees of freedom involved in the track fit procedure), at least 70 (out of 159) associated space points in the TPC, and a minimum number of two hits in the ITS, with at least one in the two innermost layers. D * + candidates were formed by combining D 0 candidates with low-p T tracks, referred to here as "soft pions", which were required to have p T > 0.1 GeV/c, |η| < 0.8, and at least three associated hits in the ITS. These selections limit the D-meson acceptance in rapidity, which drops to zero for |y| > 0.6 for p T = 1 GeV/c and |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, defined as a second-order polynomial function increasing from 0.6 to 0.8 in the range 1 < p T < 5 GeV/c, and fixed to 0.8 for p T > 5 GeV/c. The D-meson candidate selection approach adopted to reduce the combinatorial background is similar to that used in previous analyses [43,46]. The analysis procedure searches for decay vertices displaced from the primary vertex, exploiting the mean proper decay lengths of about 123 and 312 µm for D 0 and D + mesons, respectively [68]. The variables mainly used to distinguish between signal and background candidates are based on the separation between the primary and decay vertices, the displacement of the tracks from the primary vertex, and the pointing angle of the reconstructed D-meson momentum to the primary vertex, and are the same as described in [21,69]. In the strong decay of the D * + meson the primary vertex cannot be differentiated from the secondary vertex. Therefore the geometrical selections were applied on the secondary vertex topology of the produced D 0 mesons. The optimisation of the selection criteria for each D-meson species was performed as a function of p T and independently for the two centrality classes. Further reduction of the combinatorial background was obtained by applying PID for the daughter tracks with the TPC and TOF detectors. A selection in units of resolution (±3 σ ) was applied on the difference between the measured and expected signals of pions and kaons for both dE/dx and time-of-flight. The same selections are applied both for the unbiased and the ESE-selected measurements.
The D-meson elliptic and triangular flow coefficients, v 2 and v 3 , were measured using the Scalar Product (SP) method [36,70,71]. For each D-meson candidate, the v n coefficients can be expressed in terms of the Q n vectors, introduced in Sec. 2, as where u n = e inϕ D is the unit flow vector of the D-meson candidate with azimuthal angle ϕ D and the symbol "*" denotes the complex conjugation. The denominator is called the resolution (R n ) and is calculated with the formula introduced in [36], where the three subevents, indicated as A, B, and C, are defined by the particles measured in the V0C, V0A, and TPC detectors, respectively. Q Q Q k n is the subevent flow vector corresponding to the nth-order harmonic for the subevent k, and M k represents the subevent multiplicity. This is defined as the sum of the amplitudes measured in each channel for the V0A and the V0C. For the V0A and V0C detectors, the Q n vectors were calculated from the azimuthal distribution of the energy deposition, and their components are given by where the sum runs over the 32 sectors (N sectors ) of the V0A or V0C detector, ϕ k is the azimuthal angle of the centre of the sector k, and w k is the amplitude measured in sector k, once the gain of the single channels is equalised and the recentering is applied to correct effects of non-uniform acceptance [72]. For the TPC detector, the Q n vectors were computed using Eq. 2. 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 R n is obtained as a function of the collision centrality.
The v n of the D mesons cannot be directly measured using Eq. 3 as D 0 , D + , and D * + cannot be identified on a particle-by-particle basis. Therefore, a simultaneous fit to the invariant-mass spectrum and the v tot n distribution as a function of the invariant mass (M D ) was performed in each p T interval for D 0 and D + candidates separately in order to measure the raw yields and the flow coefficients. For the D * + case the distributions are studied as a function of the mass difference ∆M = M(Kππ) − (Kπ). The measured anisotropic flow coefficient, v tot n , can be written as a weighted sum of the v n of the D-meson candidate, v sig n , and that of background, v bkg n [73] as where N sig and N bkg are the raw signal and background yields, respectively. The fit function for the invariant-mass distributions was composed of a Gaussian term to describe the signal and an exponential distribution for the background for D 0 and D + candidates, while for the D * + candidates the background was described by the function a √ ∆M − m π e b(∆M−m π ) , where a and b are free parameters. In the case of the D 0 invariant mass the contribution of signal candidates with the reflected K-π mass assignment was taken into account with an additional term. Its invariant-mass distribution was parameterised with 0.14 0.145 0.15 ) and charge conj.  Figure 1: Simultaneous fits to the invariant-mass spectrum and v 2 (M D ) of D 0 (left panel), D + (middle panel), and D * + (right panel) meson candidates in the 3 < p T < 4 GeV/c, 5 < p T < 6 GeV/c, and 8 < p T < 10 GeV/c intervals, respectively, for the 30-50% centrality class. The solid blue and the dotted red curves represent the total and the combinatorial-background fit functions, respectively. For the D 0 candidates, the green dashed curve represents the contribution of the reflected signal. a double-Gaussian distribution based on Monte Carlo (MC) simulations [43,46,48,69,74]. In the MC simulation, the underlying Pb-Pb events at √ s NN = 5.02 TeV were simulated using the HIJING v1.383 generator [75] and cc or bb pairs were added with the PYTHIA 6.4.25 generator [76] with Perugia-2011 tune [77]. In the simultaneous fit, the v n parameter for the candidates with wrong K-π mass assignment was set to be equal to v sig n , provided that the origin of these candidates are real D 0 mesons. The v sig n is measured from the fit to the v tot n distribution with the function of Eq. 5, where v bkg n is a linear function for D + and D * + mesons, and D 0 mesons with p T > 4 GeV/c. For the D 0 candidates with p T < 4 GeV/c, a second-order polynomial function was used instead. Figure 1 shows the simultaneous fit to the invariantmass spectrum and v tot 2 (M D ) in the p T intervals 3-4 GeV/c for D 0 , 5-6 GeV/c for D + , and 8-10 GeV/c for D * + in the 30-50% centrality class.
The reconstructed D-meson signal is a mixture of prompt D mesons from c-quark hadronisation or strong decays of excited charmonium or open-charm states, and D mesons from beauty-hadron decays, called "feed-down" in the following. The v sig n is therefore a linear combination of prompt (v prompt n ) and feeddown (v feed-down n ) contributions, and can be expressed as where f prompt is the fraction of promptly produced D mesons estimated as a function of p T with the theory-driven method described in [21]. This method uses (i) FONLL calculations [78,79] for the production cross section of beauty hadrons, (ii) the beauty-hadron decay kinematics from the EvtGen package [80], (iii) the product of efficiency and acceptance (Acc × ε) from Monte Carlo simulations, and (iv) a hypothesis on the nuclear modification factor of feed-down D mesons.
The anisotropic flow coefficients of promptly produced D mesons were obtained assuming v feed-down The hypothesis is based on the measurement of the non-prompt J/ψ performed by CMS [19] and on the available model calculations [49,81,82], For the measurement of the modification of the p T -differential distributions of D mesons in the ESEselected samples compared to the unbiased sample, the raw yields were extracted via fits to the invariantmass distributions of D 0 , D + , and D * + candidates and normalised to the corresponding number of events in the corresponding ESE-selected sample. The same functions adopted in the simultaneous fits for the invariant-mass distributions were used. The extracted raw yields were not corrected for the efficiency in the ratio calculation, under the assumption that the reconstruction and selection efficiencies do not depend on q 2 . This assumption was previously verified in [48].

Systematic uncertainties
The D-meson v n coefficients are affected by the systematic uncertainties due to (i) the signal extraction from the invariant-mass and v tot n distributions, (ii) the beauty feed-down contribution, (iii) on the selection of the centrality interval in which R n is calculated, and (iv), for the ESE-selected samples, the uncertainties due to possible bias caused by the presence of auto-correlation effects between the subevents used for R n and q 2 calculations.
The uncertainty due to the simultaneous fit was estimated by repeating the fit several times with different configurations. In particular, the lower and upper limits of the fit range, the bin width, and the background fit functions used for the invariant-mass and v tot n distributions were varied. For each configuration the Dmeson v n was calculated and the absolute systematic uncertainty for each p T interval was assigned as the r.m.s. of the v n distribution obtained from the different trials. The absolute systematic uncertainty values on the v n are reported in Table 1 and they depend on the D-meson species, the p T interval and the ESE-selected class. This source of uncertainty was considered as uncorrelated among the p T intervals and the centrality classes for the two harmonics. The correlation between the small-q 2 /large-q 2 and the unbiased case was investigated and the outcome indicated that this uncertainty source is uncorrelated between the different q 2 -selected samples.
For the p T -differential yield ratios in ESE-selected samples, the uncertainty for the signal extraction was estimated using the same approach described above, directly on the ratio of the yields in the ESE-selected and unbiased samples, leading to a systematic uncertainty value from 0.7% to 5%, depending on the p T and the D-meson species.
The systematic uncertainty source related to the beauty feed-down correction has two main contributions. The first is due to the f prompt calculation and it was studied by varying the quark mass and the renormalisation and factorisation scales in the FONLL calculations, the R feed-down  Table 1 and they depend on the D-meson species, the p T interval and the ESE-selected class. The uncertainty due to the beauty feed-down correction was assumed to be fully correlated among the p T bins for the measured v n coefficients in the same centrality class.
The non-flow effects are naturally suppressed because of the pseudorapidity gap of at least 0.9 units between the pseudorapidity interval used for the D-meson reconstruction, and the V0C used for the Q n -vector determination. Furthermore, the auto-correlation effect due to the usage of the TPC tracks for the q 2 estimate has been discussed in Sec. 2 and the related systematic uncertainty was found to be negligible, as described in [48].
The contribution of the R n to the systematic uncertainty is due to the centrality dependence. The central value of R n was estimated using the three subevent formula, as described in Sec. 3, averaged over the events in the 0-10% and 30-50% intervals. The uncertainty was evaluated as the difference of the centrality integrated R n values with those obtained as weighted averages of R n values in narrow centrality intervals using the D-meson yields as weights. A systematic uncertainty of 3.5% and 0.5% was assigned on R 2 in the 0-10% and 30-50% centrality classes and for all ESE-selected samples. For the R 3 , an uncertainty of 0.5% was assigned in the 30-50% interval while it was found to be negligible for the 0-10% class. The uncertainty associated with the resolution factor is smaller for the third harmonic than for the second harmonic, due to the milder centrality dependence of R 3 compared with that of R 2 .
For the ESE-selected samples an additional source of systematic uncertainty on the resolution originates from auto-correlations due to the usage of the TPC tracks both for q 2 and R 2 determination. This potential bias is assessed by replacing the ratio Q Q Q V0C n /M V0C ·Q Q Q TPC * n /M TPC / Q Q Q V0A n /M V0A ·Q Q Q TPC * n /M TPC in Eq. 3 with the one from the q 2 -integrated analysis, following the same approach used for the J/ψ azimuthal anisotropy measurement [83]. In this case, the systematic uncertainty was estimated to be 3.5% for the small-q 2 and 1% for the large-q 2 samples in the 0-10% centrality class, and 0.5% for both q 2selected classes in the 30-50% centrality class, as reported in Table 1. The last two sources of systematic uncertainty, related to the resolution, are considered to be fully correlated among the different p T intervals.
For the analysis of the p T -differential yield ratios in ESE-selected and unbiased samples the reconstruction efficiency was verified to be independent of q 2 . Consequently, it cancels out in the ratio of the two ESE-selected classes. by using the inverse squared absolute statistical uncertainties as weights, after having compared their compatibility [67]. The systematic uncertainties were propagated to the average by considering the contributions from the centrality dependence of the R n resolution and the correction for the beauty feeddown component in the D-meson yields as correlated among the D-meson species. The D-meson v n harmonics are compared to the corresponding coefficients measured for charged pions and protons at midrapidity (|y| < 0.5) [39] as well as to inclusive J/ψ mesons at forward rapidity (2.5 < y < 4) [84].

Unbiased flow harmonics
The D-meson elliptic flow increases significantly from central to semi-central collisions, as expected from the increasing eccentricity of the interaction region. Conversely, the triangular flow is compatible in the two centrality classes within the large uncertainties, following the milder centrality dependence of the third flow harmonic observed for light-flavour particles [39]. For p T < 3-4 GeV/c (p T < 4-5 GeV/c) the measured D-meson v 2 (v 3 ) is lower than that of pions and protons. This observation is consistent within uncertainties with the hypothesis of a mass hierarchy, v n (D) < v n (p) < v n (π), in the low p T region (p T 3 GeV/c). In semi-central events the v n coefficients of J/ψ mesons seem to follow the mass hierarchy (v n (J/ψ) < v n (D)). In central events the data suggests a similar behaviour, however within the current uncertainties no firm conclusions can be drawn. This observation can be explained  by the interplay between the anisotropic flow and the isotropic expansion of the system (radial flow), which imposes an equal velocity boost to all particles. For 4 p T 6-8 GeV/c, the D-meson v n coefficients are similar to those of charged pions and lower than those of protons. This observation is consistent with a scaling of the v n coefficients with the number of constituent quarks, which supports the hypothesis of particle production via quark coalescence [85]. In the same p T interval, for the 30-50% centrality class the larger values of v n for D mesons compared to J/ψ mesons can be explained by (i) the hadronisation via coalescence together with the larger flow coefficients of up and down quarks compared to that of charm quarks [28] and (ii) the fraction of J/ψ mesons coming from beauty-hadron decays [86,87], which are expected to have lower v 2 and v 3 than charmed mesons [26,88]. In the 0-10% centrality class, the current experimental uncertainties do not allow for firm conclusions on the expected difference for J/ψ and D mesons. The measured v n coefficients for all the hadron species are compatible within uncertainties for p T 8 GeV/c. Similar values of v n coefficients are expected, because in this kinematic range the charm-quark mass is small compared to the momentum, and because the path-length dependence of the in-medium parton energy loss is similar for high-p T charm quarks and gluons.
In Fig. 3, the average D-meson v n coefficients are compared to theoretical calculations that include the charm-quark transport in a hydrodynamically expanding medium. The theoretical uncertainties,  [34,56], PHSD [53], Catania [94,95], and BAMPS el [49] calculations the interactions between the charm quarks and the medium constituents are modelled with collisional processes, while the MC@sHQ+EPOS2 [26], LBT [57,90], LIDO [91,92], BAMPS el+rad [55], DAB-MOD(M&T) [37,93], and LGR [96] models include also radiative processes. The difference in the variants of the BAMPS model indicates that in this model elastic collisions are the dominant process that imparts a positive D-meson v 2 in the low and intermediate p T region. All the models except for BAMPS include the hadronisation of the charm quark via coalescence, in addition to the fragmentation mechanism. Initial-state event-by-event fluctuations are included in the POWLANG HTL, LIDO, PHSD, MC@sHQ+EPOS2, LBT, and DAB-MOD(M&T) models, which are therefore the only ones that provide predictions for the triangular flow. Although the models differ in several aspects related to the interactions both in the QGP and in the hadronic phase as well as to the medium expansion, most of them provide a fair description of the measured v n harmonics. The largest difference is observed in the 2 < p T < 6 GeV/c interval for the v 2 in the 30-50% centrality class, where most of the models provide a prediction lower than the measured points. This is more evident for the LIDO model, which shows a deviation of 5.4 σ , and BAMPS el+rad , which underestimates the measured v 2 by about a factor two with more than 10 σ significance. In contrast to this, BAMPS el overestimates the measurement by about 3 σ . The underestimation of the data by the BAMPS el+rad model can be eventually due to the missing implementation of the charm-quark coalescence with light quarks from the medium, which seems to be necessary in the description of the measured v 2 . In the same p T range, the DAB-MOD model overestimates the measured v 2 in the 0-10% centrality class by 3.7 σ . These discrepancies expressed in number of standard deviations were computed combining the probability to observe a deviation from the null hypothesis (i.e. the model prediction) for all the measured points in the 2 < p T < 6 GeV/c interval, considering both the experimental (statistical and systematic) and the theoretical uncertainties, when available. The global agreement between the data and the theoretical models was evaluated by computing the χ 2 /ndf, as done in [46]. The values are reported in Table 2. All the centrality classes and v n harmonics were considered when the model predictions were available. Compared to the results in [46], for almost all the models the χ 2 /ndf is found to be higher than unity, most likely because of the improved precision of the measurement. The models that describe p T and ESE dependence of D-meson v n harmonics ALICE Collaboration the data with χ 2 /ndf < 2 are MC@sHQ+EPOS2, LBT, LGR, PHSD, POWLANG, Catania, and TAMU which is more in agreement with the data compared to [46], thanks to the improved description of the charm-quark coalescence in its latest version [89]. These models use a value of heavy-quark spatial diffusion coefficient in the range 1.5 < 2πD s T c < 7 at the critical temperature T c = 155 MeV [97], which is consistent with the interval obtained in [46]. It is however important to consider that not all the theoretical models provide predictions for all the v n harmonics in all the centrality classes reported in this article, hence the global interpretation of these comparisons could not be conclusive.
5.2 Event-shape engineered flow harmonics and p T p T p T -differential yields The average v 2 of prompt D 0 , D + , and D * + mesons measured in the ESE-selected samples is shown in Fig. 4 for the 0-10% (top row) and the 30-50% (bottom row) centrality classes. The measurements in the small-q 2 sample are reported in the left column, those in the large-q 2 sample in the right column, while the measurements in the unbiased samples recomputed in the same p T intervals of the ESE analysis are in the middle column. A reduced p T range (2 < p T < 16 GeV/c) and wider p T intervals compared to the unbiased v 2 measurement were adopted due to the limited size of the ESE-selected samples. The average v 2 among the three D-meson species was computed as described in Sec. 5.1. In Fig. 5 the ratio between the average D-meson v 2 measured in the ESE-selected samples with respect to that in the unbiased sample is depicted. The statistical uncertainties of the ratio were calculated taking into consideration the degree of correlation between the measurements in the ESE-selected and unbiased samples. The systematic uncertainties arising from the centrality dependence of R n , the non-flow contaminations among subevents, and the correction for the beauty feed-down contribution were considered as fully correlated.
The D-meson v 2 was found to be on average about 50% higher (lower) in the 20% of the events with largest (smallest) q 2 in both the 0-10% and 30-50% centrality classes. No significant centrality dependence was found within the current uncertainties. The corresponding variation of the average q 2 in the small-q 2 (large-q 2 ) sample with respect to the unbiased one was found to be about 65% (75%) and 60% (65%) for the 0-10% and 30-50% centrality class, respectively. This confirms the correlation between the D-meson azimuthal anisotropy and the collective expansion of the bulk matter already observed in [48]. This modification of the v 2 coefficient was found to be independent of p T within uncertainties, which might suggest that the ESE selection is related to a global property of the events (i.e. a property that is independent of the measured particle and is related to the entire event). A similar trend was also observed for light-flavour particles [66].
Figures 4 and 5 also compare the measured v 2 and v 2 ratios between ESE-selected and unbiased samples to the POWLANG, LIDO, DAB-MOD, and Catania theoretical predictions. For the POWLANG model, both the predictions obtained with the transport coefficients from weak coupling (Hard Thermal Loop, HTL [98]) and from lattice QCD calculations (lQCD [99]) are reported. For the DAB-MOD model, a version based on the heavy-quark transport (M&T [32]) and a parametric model for the heavy-quark energy loss (E loss [100]) were considered. In the LIDO, DAB-MOD, and Catania models the ESE selection is performed with a q 2 estimator computed starting from generated quantities [91,93,95], while in the POWLANG model the elliptic eccentricity ε 2 is directly used [34]. The v 2 measured in the small-q 2 sample is described by all the available models within the uncertainties. On the contrary, in the 30-50% centrality class the LIDO, DAB-MOD, and Catania models underestimate the measurement in the large-q 2 sample, which is instead well described by the POWLANG HTL prediction.
In the case of POWLANG lQCD, the theoretical prediction is compatible with the measured v 2 for p T < 4 GeV/c and lower for higher p T . The DAB-MOD calculations give a better description of the experimental data with the M&T approach for p T < 5 GeV/c and in the E loss case for p T > 5 GeV/c. When the ratios between the v 2 in the ESE-selected and the unbiased samples are considered, the models seem to better describe the measured values, owing to similar discrepancies between theoretical predictions and experimental data in the ESE-selected and unbiased samples which lead to similar ratio values in the different models. In the small-q 2 samples the model predictions are more similar TeV in the small-q 2 , large-q 2 (see text for details), and unbiased samples, for the 0-10% (top panels) and 30-50% (bottom panels) centrality classes, compared to model calculations [34,37,91,[93][94][95]. In the LIDO, DAB-MOD, and Catania predictions, the ESE selection is performed with a q 2 estimator, while in the POWLANG model the elliptic eccentricity ε 2 is used.
to each other and the discrepancies are less significant, also due to the larger experimental uncertainties. Interestingly, different implementations of the same model with the studied transport parameterisations (i.e. POWLANG HTL vs. POWLANG lQCD, and DAB-MOD(M&T) vs. DAB-MOD(E loss )) give similar predictions, suggesting that the effect of the ESE selection is more related to the initial geometry and the underlying hydrodynamic expansion rather than the dynamic evolution of the heavy quarks in the medium.
To study a possible interplay between the azimuthal anisotropy of the event and the charm-quark radial flow (at low/intermediate p T ) and in-medium energy loss (at high p T ), the ratio of the measured per-event yields of prompt D 0 , D + , and D * + mesons in the ESE-selected and unbiased samples has been calculated as a function of p T in the range 2 < p T < 24 GeV/c. The average D-meson ratios, computed by using the inverse of the squared relative statistical uncertainties as weights, are compared to the POWLANG and LIDO models in Fig. 6. The POWLANG model predicts a hardening (softening) of the p T distributions in the large (small)-q 2 class of events due to an interplay between the radial and elliptic flows, while no significant modification is predicted by the LIDO model. Within the current precision, the measured per-event yield ratios and are found compatible with unity, and hence to the LIDO model predictions, and with the POWLANG model in the case of lQCD, while the measured effect seems to be lower than the effect predicted with HTL transport coefficients.  Figure 5: Ratio of the average prompt D 0 , D + , and D * + meson v 2 coefficients measured in the small-q 2 (left panels) and large-q 2 (right panels) selected samples with respect to that of the unbiased sample as a function of p T in Pb-Pb collisions at √ s NN = 5.02 TeV for the 0-10% (top panels) and 30-50% (bottom panels) centrality classes, compared to model calculations [34,37,91,93,95].

Conclusions
The elliptic and triangular flow of D 0 , D + , and D * + mesons was measured with the SP method at midrapidity (|y| < 0.8) as a function of p T in central (0-10%) and semi-central (30-50%) Pb-Pb collisions at √ s NN = 5.02 TeV.
Compared to other particle species, the average D-meson v n harmonics were found to be compatible with the hypothesis of a mass hierarchy for p T 3 GeV/c as observed for light-flavour hadrons [39]. At intermediate p T , the D-meson v n is similar to those of charged pions, lower than those of protons, and higher than those of J/ψ mesons, supporting the hypothesis of charm-quark hadronisation via coalescence. Moreover, the contribution to the hadronisation of charm quarks from coalescence with light quarks from the medium seems to be necessary in the theoretical models to quantitatively reproduce the measured D-meson v n . For p T 8 GeV/c, the D-meson v 2 and v 3 are compatible within uncertainties with the values measured for the other particle species, indicating a similar path-length dependence of   Figure 6: Average of the ratio of p T -differential D 0 , D + , and D * + yields measured in the ESE-selected samples to those in the unbiased sample in Pb-Pb collisions at √ s NN = 5.02 TeV for the 0-10% (top panels) and 30-50%
the energy loss of high-p T charm quarks and gluons. The comparison of the measured D-meson v n with theoretical calculations suggests that the interactions with the hydrodynamically expanding medium impart a positive v 2 and v 3 to the charm quarks.
The elliptic flow and the modification of the p T distributions of D 0 , D + , and D * + mesons were also investigated with the event-shape engineering technique. The D-meson v 2 was found to be larger (smaller) in events with larger (smaller) q 2 , confirming the correlation with average bulk elliptic flow. The ratios of the p T -differential yields measured in the ESE-selected samples and the unbiased sample were found to be compatible with unity. The measurements in the ESE-selected samples are qualitatively described by theoretical calculations and provide new constraints to models based on charm-quark transport in a hydrodynamically expanding medium and charm-quark energy loss in the QGP.

Acknowledgements
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The