Event-shape engineering for the D-meson elliptic flow in mid-central Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

The production yield of prompt D mesons and their elliptic flow coefficient $v_2$ were measured with the Event-Shape Engineering (ESE) technique applied to mid-central (10-30% and 30-50% centrality classes) Pb-Pb collisions at the centre-of-mass energy per nucleon pair $\sqrt{s_{\rm NN}} =5.02$ TeV, with the ALICE detector at the LHC. The ESE technique allows the classification of events, belonging to the same centrality, according to the azimuthal anisotropy of soft particle production in the collision. The reported measurements give the opportunity to investigate the dynamics of charm quarks in the Quark-Gluon Plasma and provide information on their participation in the collective expansion of the medium. D mesons were reconstructed via their hadronic decays at mid-rapidity, $|\eta|<0.8$, in the transverse momentum interval $1<p_{\rm T}<24$ GeV/$c$. The $v_2$ coefficient is found to be sensitive to the event-shape selection confirming a correlation between the D-meson azimuthal anisotropy and the collective expansion of the bulk matter, while the per-event D-meson yields do not show any significant modification within the current uncertainties.


Introduction
Quantum Chromo-Dynamics (QCD) calculations on the lattice predict the existence of a plasma of deconfined quarks and gluons, known as the Quark-Gluon Plasma (QGP) [1][2][3][4]. The transition from the hadronic phase to the QGP state occurs at high temperatures and energy densities, which can be reached in collisions of heavy nuclei at ultra-relativistic energies. The QGP created in ultra-relativistic heavy-ion collisions was found to behave as a nearly ideal fluid (i.e. with a small shear viscosity over entropy density ratio, η/s), undergoing an expansion that can be described by relativistic hydrodynamics [5][6][7][8][9][10].
Heavy flavours (charm and beauty quarks), due to their large masses, m c ≈ 1.3 GeV/c 2 and m b ≈ 4.5 GeV/c 2 , are predominantly produced in hard-scattering processes characterised by timescales shorter than the QGP formation time [11][12][13][14]. Thus, they experience the entire evolution of the medium interacting with its constituents via inelastic (gluon radiation) [15][16][17] and elastic (collisional) [18] QCD processes. Such interactions with the medium constituents can also lead to a modification of the hadronisation mechanism with respect to the fragmentation in vacuum: a significant fraction of low-and intermediate-momentum charm and beauty quarks can hadronise via recombination with other quarks from the medium [19][20][21].
Heavy-flavour hadrons are effective probes of the properties of the medium produced in heavy-ion collisions. A strong modification of their transverse-momentum (p T ) distributions in heavy-ion collisions with respect to pp collisions was observed at RHIC [22][23][24][25] and LHC energies [26][27][28][29][30][31][32]. In particular, the observed 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 provides compelling evidence of the heavy-quark energy loss in deconfined strongly-interacting matter [13,17].
Further insight into the interactions of heavy quarks with the medium can be obtained through measurements of the azimuthal distributions of heavy-flavour hadrons in heavy-ion collisions. The initial spatial anisotropy present in the early stages of nucleus-nucleus collisions is converted via multiple interactions into an azimuthally anisotropic distribution in momentum space of the produced particles [33,34]. This anisotropy can be characterised in terms of the Fourier coefficients v n of the azimuthal distribution of particle momenta relative to the symmetry-plane angles Ψ n (for the n th harmonic) [34,35]. The values of the Fourier coefficients depend on the geometry of the collision, the fluctuations in the distributions of nucleons within the nuclei [36], and the dynamics of the expansion. The second-order coefficient v 2 = cos[2(ϕ − Ψ 2 )] , where ϕ is the particle momentum azimuthal angle and the brackets indicate the average over all the measured particles in the considered events, is usually denoted as elliptic flow. In non-central heavy-ion collisions, it represents the dominant term in the Fourier expansion [33,35]. The measurement of the azimuthal anisotropy of heavy-flavour hadrons at low p T is sensitive to whether charm quarks take part in the collective expansion of the medium [37], as well as to the fraction of heavy-flavour hadrons hadronising via recombination with flowing light quarks [38,39]. At high p T , it can constrain the path-length dependence of heavy-quark in-medium energy loss [40,41]. A positive v 2 in the heavy-flavour sector was observed at RHIC in Au-Au collisions at a centre-of-mass energy per nucleon pair √ s NN = 200 GeV [22,42,43 . The D-meson results are described by theoretical calculations including mechanisms that impart a positive v 2 to charm quarks through the interactions with the hydrodynamically-expanding medium, namely collisional processes, and recombination of charm and light quarks [50][51][52][53][54][55][56][57][58][59]. According to these model calculations, the same mechanisms affect the beauty-quark propagation in the medium, although the beauty-hadron v 2 is expected to be smaller than that of charm hadrons and to have a different transverse momentum dependence due to the large mass of the b quarks. Precise measurements of v 2 of heavy-flavour hadrons help to constrain model parameters, e.g., the heavy-quark spatial diffusion coefficient D s in the QGP, ESE for the D-meson v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration which is related to the relaxation time (or the time scale for equilibration) of the heavy quarks inside the QGP [48,60].
The Event Shape Engineering (ESE) technique [61] can be used to further investigate the dynamics of heavy quarks in the medium. This technique has already been exploited in the light-flavour sector to study the interplay between the initial geometry of the nucleus-nucleus collisions and the subsequent evolution of the system [62,63], and to investigate the Chiral Magnetic Effect (CME) [64,65]. The ESE technique is based on the observation of a large event-by-event v n variation at fixed collision centrality [66]. Hydrodynamic calculations show a linear correlation between the final state v 2 (and v 3 ) and the corresponding eccentricities in the initial state ε 2 (and ε 3 ) for small values of η/s [67][68][69]. These observations suggest the possibility to select heavy-ion collisions with different initial geometrical shape on the basis of the magnitude of the average bulk flow.
The ESE technique provides a tool to investigate the correlation between the flow coefficients of D mesons and soft hadrons: measuring the D-meson v 2 in classes of events in a given centrality interval, but with different magnitude of the average event flow can be useful to study the interplay between the anisotropic flow of heavy quarks and that of the bulk matter. In addition, it could provide insights on how the fluctuations in the initial geometry of the system affect the path-length-dependent energy loss experienced by the heavy quarks in the QGP. For these reasons, the application of the ESE technique to the D-meson v 2 measurements could be exploited to infer more information on the dynamics of the charm quark in the QGP and has the potential to set additional constraints on parameters of model calculations implementing heavy-quark transport in an hydrodynamically expanding medium [70]. Model calculations for the correlation between the v 2 values of soft hadrons and heavy-flavour mesons on an event-by-event basis have recently become available. A linear correlation between the high-p T D-meson v 2 , which originates from the path-length dependence of in-medium energy loss, and the elliptic flow of charged hadrons is predicted in [71], based on a model for charm-quark energy loss in a medium described event-by-event with viscous hydrodynamics. Within the heavy-quark transport model of [72], an almost linear correlation is obtained between the v 2 of pions and that of D 0 mesons with p T > 2 GeV/c, which is dominated by low-p T mesons and is therefore sensitive to the degree of thermalisation of charm quarks with the collectively expanding medium. According to these calculations, 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.
Finally, the measurement of the D-meson yields at low and intermediate p T in ESE-selected events allows the investigation of a possible interplay between elliptic and radial flow, already observed for charged and identified particles [62]. This observation is possibly related to the correlation between the density of participant nucleons and the initial eccentricity of the collision. For high-p T D mesons, the measurement of the yields in collisions with different initial eccentricity via the ESE technique could further constrain in-medium energy loss models.
In this paper, the D 0 , D + and D * + meson v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV in the 10-30% centrality class are presented and compared to the results in the 30-50% centrality class published in [48]. The D 0 and D + v 2 obtained with ESE and the measurement of D-meson yield ratios in ESE-selected events in the 10-30% and 30-50% centrality classes are reported as well.

Data analysis
The D 0 , D + and D * + mesons were reconstructed at mid-rapidity, exploiting the tracking and particle identification capabilities of the ALICE detector at the LHC. A detailed description of the ALICE experimental apparatus and its performance can be found in [73,74]. The main detectors used for the analysis presented in this paper are the Inner Tracking System (ITS), a six-layer silicon detector ESE for the D-meson v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration used to track charged particles and for the reconstruction of primary and secondary vertices; the Time Projection Chamber (TPC), which provides track reconstruction as well as particle identification via the measurement of the specific ionisation energy loss dE/dx; the Time-Of-Flight (TOF) detector, an array of Multi-Gap Resistive Plate Chambers that provides particle identification via the measurement of the flight time of the particles. These detectors cover the pseudorapidity interval |η| < 0.9 and are located in a large solenoidal magnet providing a uniform magnetic field of 0.5 T parallel to the LHC beam direction. In addition, two detectors were used for the event selection and classification: the V0 detector, which consists of two arrays of 32 scintillators each, covering the full azimuth in the pseudorapidity intervals −3.7 < η < −1.7 (V0C) and 2.8 < η < 5.1 (V0A); and the Zero Degree Calorimeters (ZDC), located at 112.5 m from the interaction point on either side, to detect spectator neutrons and protons of the colliding nuclei.
The analysed data sample consists of Pb-Pb collisions at √ s NN = 5.02 TeV collected using a minimumbias interaction trigger that required coincident signals in both scintillator arrays of the V0 detector. Events were selected offline by using the V0 and the neutron ZDC timing information, to remove contaminations produced by the interaction of the beams with residual gas in the vacuum pipe. Only events with a reconstructed primary vertex within ±10 cm from the centre of the detector along the beam line were analysed. Events satisfying the aforementioned selections were divided in centrality classes, defined in terms of percentiles of the hadronic Pb-Pb cross section. This classification was based on a fit to the sum of the signal amplitudes measured in the V0 detectors. The fit function assumes the Glauber model [75,76] combined with a two-component model for particle production [77]. The number of events in each centrality class considered for this analysis (10-30% and 30-50%) is about 20.7 × 10 6 , corresponding to an integrated luminosity of about 13 µb −1 . The events in each centrality class were further divided in 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 [69,78], defined as where M is the multiplicity (number of tracks used in the q 2 calculation) and is the second-order flow vector, which is built starting from the azimuthal angles (ϕ i ) of the considered particles. The denominator in Eq. 1 is introduced to remove the dependence of |Q Q Q 2 | on √ M in the absence of flow [69].
The Q Q Q 2 vector was measured using charged tracks reconstructed in the TPC (q TPC 2 ), with |η| < 0.8 and 0.2 < p T < 5 GeV/c, to exploit the good ϕ resolution of the TPC and the large multiplicity at midrapidity, which are crucial to maximise the selectivity of q 2 . In order to remove autocorrelations between D mesons and q 2 , the tracks used to form the D-meson candidates were excluded from the computation of q 2 . However, with this definition of q 2 , some residual non-flow correlations (i.e. correlations among particle emission angles not induced by the collective expansion but rather by particle decays and jet production) could still be included. As shown in [62], the introduction of a pseudorapidity separation of more than one unit between the region used to calculate q 2 and the region used to measure the observables would suppress unwanted non-flow contributions. Therefore, to investigate a possible effect induced by non-flow contaminations, q 2 was also measured using the V0A detector (q V0A 2 ), allowing for a pseudorapidity separation of at least 2 units between the D-meson decay tracks and the particles used for the q 2 determination. In this case, the Q Q Q 2 vector was calculated from the azimuthal distribution of the energy deposition measured in the V0A detector, and its components are given by where the sum runs over the 32 sectors (N sectors ) of the V0A detector, ϕ i is the angle of the centre of the sector i and w i is the amplitude measured in sector i, once the gain equalisation method [79] is applied to correct effects of non-uniform acceptance. The comparison between the two ESE selections is discussed in Sec. 2.1.
The left panel of Fig. 1 shows the centrality dependence of the q TPC 2 distribution. As expected in case of large initial-state fluctuations, the q 2 distribution is broad and reaches values larger than twice the mean value [61]. Moreover, because of the different average elliptic flow and multiplicity, the q 2 distribution changes as a function of centrality. Hence, a selection on a fixed value of q 2 would induce a non-flat centrality distribution, that would spoil the event-shape selection. For this reason, the selection of the events according to their q 2 was performed by defining q 2 percentiles in 1%-wide centrality intervals. The results presented in the following sections are obtained in two ESE-selected classes, corresponding to the 60% and the 20% of events with smallest and largest q 2 , respectively. The q TPC 2 distributions for these classes in the 30-50% centrality interval are displayed in the right panel of Fig. 1. In the following, we will refer to these two classes as "small-q 2 " and "large-q 2 ". In case of no event-shape selection, we will use the "unbiased" term.
The D mesons, together with their charge conjugates, were reconstructed via their hadronic decay channels D 0 → K − π + , D + → K − π + π + and D * + → D 0 π + → K − π + π + . The D-meson candidates were built combining pairs and triplets of tracks with proper charge sign, |η| < 0.8, p T > 0.4 GeV/c, a minimum number of 70 (out of 159) associated space points in the TPC and no less than two hits (out of six) in the ITS, with at least one in the two innermost layers. For the soft pion produced in the D * + decay, also tracks reconstructed only in the ITS, with at least three associated hits and with p T > 0.1 GeV/c, were considered. These selections limit the D-meson rapidity acceptance, which drops steeply to zero for |y| > 0.6 for p T = 1 GeV/c and |y| > 0.8 for p T > 5 GeV/c. Therefore, a p T -dependent fiducial acceptance selection, |y D | < y fid (p T ), was applied. The selection value, y fid (p T ), was defined according to a second-order polynomial function, increasing from 0.6 to 0.8 in the range 1 < p T < 5 GeV/c, and fixed to a constant value of 0.8 for p T > 5 GeV/c.  0.7597 ± 0.0001 Table 1: Event-plane resolution R 2 in the 10-30% and 30-50% centrality classes for the unbiased, small-q 2 and large-q 2 samples. The quoted uncertainty is statistical only. on the reconstruction of secondary vertices with a separation of a few hundred microns from the primary vertex (cτ 123 and 312 µm for D 0 and D + , respectively [80]). The main variables used to enhance the statistical significance and the signal-to-background ratio are the displacement of the decay tracks from the interaction vertex, the separation between the secondary and primary vertices and the pointing angle of the reconstructed D-meson momentum to the primary vertex. In the case of the strong decay of the D * + meson, the secondary vertex cannot be resolved from the primary vertex, and therefore geometrical selections were applied on the displaced decay-vertex topology of the produced D 0 mesons. In addition, for D 0 and D + mesons, the normalised difference between the measured and expected transverse-plane impact parameters of each of the decay particles and the transverse-plane impact parameter to the primary vertex (d xy 0 ) of the D + -meson candidates were applied to suppress the fraction of D mesons coming from beauty-hadron decays (denoted by "feed-down" in the following) and hence reduce the associated systematic uncertainty. These selections were found to be especially effective for D + mesons, for which a rejection of the feed-down contribution up to 50% at high p T was achieved. The selection criteria for each D-meson species were optimised as a function of p T independently for the two centrality classes, because of the different combinatorial background. Within a given centrality class, the same selection criteria were applied in the different ESE-selected samples. In order to further reduce the combinatorial background, a particle identification for charged pions and kaons with the TPC and TOF detectors was applied, using a selection in units of resolution (at ±3 σ ) around the expected mean values of dE/dx and time of flight, respectively.
Monte Carlo simulations with a detailed description of the detector and its response, based on the GEANT3 transport package [81], were used to study the signal invariant-mass distributions and the reconstruction efficiencies, as described in the following. In the Monte Carlo simulation, the underlying Pb-Pb events at √ s NN = 5.02 TeV were simulated using the HIJING v1.383 generator [82] and cc or bb pairs were added with the PYTHIA v6.421 generator [83] with Perugia-2011 tune [84]. The generated D mesons were forced to decay into the hadronic channels of interest for the analysis.
The D-meson elliptic flow, v 2 , is measured using the Event-Plane (EP) method [35]. This analysis technique relies on the event-by-event estimate of the second-order harmonic symmetry plane Ψ 2 using the so-called event-plane angle For the measurements of v 2 in the unbiased and q TPC 2 -selected samples, Q Q Q 2 was estimated with the full ESE for the D-meson v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration   V0 detector using Eq. 3 (with N sectors corresponding to the 64 sectors of the full V0 detector). In case of the ESE selection based on q V0A 2 , only the 32 sectors of the V0C were used for the ψ 2 determination, to avoid autocorrelations with the q 2 measurement.
After the topological and kinematical selections, the D-meson candidates were divided in two samples, according to their azimuthal angle relative to the event-plane angle ∆ϕ = ϕ D − ψ 2 , namely in-plane (] − π 4 , π 4 ] and ] 3π 4 , 5π 4 ]) and out-of-plane (] π 4 , 3π 4 ] and ] 5π 4 , 7π 4 ]). The separation of at least 0.9 units of pseudorapidity (|∆η| > 0.9) between the D mesons and the particles used to measure ψ 2 , naturally ensured by the selection of D-meson decay tracks and the V0 (V0C) acceptance, suppresses non-flow contributions. The v 2 can therefore be expressed by the following equation where N in-plane and N out-of-plane are the D-meson raw yields in the two ∆ϕ intervals. The raw yields can be directly used in Eq. 5, without an efficiency correction, since simulations showed that the D-meson reconstruction and selection efficiencies do not depend on ∆ϕ [45]. The factor 1/R 2 in Eq. 5 is the correction due to the finite resolution of the estimated ψ 2 angle. In case of v 2 measurements in the unbiased and q TPC 2 -selected samples, the event-plane resolution R 2 was determined by correlating three sub-events of charged particles reconstructed in the V0 itself, in the positive (0 < η < 0.8) and negative (−0.8 < η < 0) semivolumes of the TPC [35]. In case of ESE selection based on q V0A 2 , the three subevents considered were the charged particles reconstructed in the V0C, in the V0A, and in the full volume of the TPC (|η| < 0.8). The values of R 2 estimated in the 10-30% and 30-50% centrality classes, for the unbiased, small-q 2 and large-q 2 samples are reported in Tab. 1. The R 2 factor is higher (lower) in the large-q 2 (small-q 2 ) class with respect to that evaluated for the unbiased sample and similarly in the V0 case than in the V0C one, since the event-plane resolution The in-plane and out-of-plane raw yields were obtained by fitting the invariant-mass distributions M(Kπ) for D 0 candidates, M(Kππ) for D + candidates and the mass-difference ∆M = M(Kππ) − M(Kπ) distributions for D * + candidates in each centrality class. The fit function was composed by a Gaussian distribution to describe the signal and an exponential term for the background of D 0 and D + candidates or by a threshold function multiplied by an exponential function, a √ ∆m − m π · e b(∆m−m π ) , for the D * + background. Since the invariant-mass resolution does not exhibit any dependence on ∆ϕ or q 2 , the width of the Gaussian, for each D-meson species and p T interval, was fixed to that obtained from a fit to the invariant-mass distribution integrated over ∆ϕ and q 2 , where the signal has higher statistical significance.   In addition, for the determination of the D 0 -meson yield, the contribution of signal candidates present in the invariant-mass distribution with the wrong K-π mass assignment was taken into account by including an additional term in the fit function, parametrised with a double-Gaussian shape [45] determined with Monte Carlo simulations. The contribution of the reflected signal, 2-5% under the D 0 -peak region depending on p T , was considered as background and therefore not included in the raw yield. Examples of invariant-mass fits for the three D-meson species in the unbiased sample in the 10-30% centrality class and for D 0 and D + mesons in the ESE-selected samples in the 30-50% centrality class are shown in Fig. 2 and Fig. 3, respectively.
The measured raw D-meson yields contain a feed-down contribution which, depending on the D-meson species, p T and the topological selections, can vary between 5% and 20%. The strategy adopted to correct the observed v 2 for the fraction of prompt D mesons in the measured raw yields is the same as the one used in [48], and it is described in the following. The observed v 2 can be expressed as a linear combination of the prompt (D mesons coming directly from the hadronisation of a c-quark or from the decay of an excited open charm or charmonium state) and the feed-down contributions v obs 2 = f prompt v prompt 2 where f prompt is the fraction of promptly produced D mesons estimated as a function of p T with the same method used in [32]. In particular, it is computed using (i) FONLL calculations [85,86] for the ESE for the D-meson v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration production cross-section of beauty hadrons, (ii) the beauty-hadron decay kinematics from the EvtGen package [87], (iii) the product of efficiency and acceptance (Acc × ε) from Monte Carlo simulations and (iv) an hypothesis on the nuclear modification factor of feed-down D mesons. The nuclear modification factor is defined as R AA = (dN AA /dp T )/( T AA dσ pp /dp T ), where dN AA /dp T and dσ pp /dp T are the p Tdifferential yield and production cross section of D mesons in nucleus-nucleus (AA) and pp collisions, respectively, and T AA is the average nuclear overlap function in the considered centrality class [77]. The hypothesis R feed-down AA = 2R prompt AA was used to estimate the central value of f prompt . This choice is motivated by the comparison of the R AA of prompt D mesons at √ s NN = 2.76 TeV [88] with that of J/ψ from beauty-hadron decays at the same energy measured by the CMS Collaboration [30], which indicates that the charm-hadron production yield is more suppressed than that of the beauty hadrons by about a factor of two. This difference is described by model calculations with parton-mass-dependent energy loss [53]. The selection efficiency and therefore f prompt are different in the 10-30% and 30-50% centrality classes, because of the different geometrical selections applied on the displaced decayvertex topology. In the case of the ESE selection, the (Acc × ε) is the same for the large-q 2 and small-q 2 samples, because the same selection criteria were used in the two ESE-selected classes and the efficiency was found not to depend on local particle density. Therefore, considering also the same can be attributed to different contributions of non-flow correlations, but also to the different eccentricity discriminating power of q 2 measured with the two detectors. This discriminating power depends on the magnitude of the elliptic flow, on the multiplicity used in the q 2 calculation and on the performance of the detector (i.e. the angular resolution or the linearity of the response as a function of charged-particle multiplicity). To disentangle the two effects, the selectivity of q TPC 2 was artificially reduced by rejecting randomly 85% of tracks used for the calculation of q TPC 2 . A similar strategy was used in [62]. Figure 4 illustrates the comparison of the effect of the ESE selection on the D 0 -meson v 2 obtained using q TPC 2 (left-hand panels), q TPC 2 with random rejection of 85% of the tracks (middle panels) and q V0A 2 (right-hand panels), for the 10-30% (top panels) and the 30-50% (bottom panels) centrality classes. The separation between the measurements in the ESE-selected sample with respect to the unbiased one is reduced in the case of the random rejection of the tracks with respect to the default q TPC 2 , confirming the reduced q TPC 2 selectivity. The results obtained with q V0A 2 are similar to those obtained reducing artificially the selectivity of q TPC 2 , although they are compatible within uncertainties with both q TPC 2 measurements. This indicates that the statistical precision of the measurement is not sufficient to draw a firm conclusion about non-flow contaminations in the measurement performed by selecting the events according to q TPC 2 . The q TPC 2 -based selection was thus chosen for the evaluation of the results presented in the following sections, except for the comparison of the effect of the ESE selection on the D-mesons and the chargedparticle v 2 , for which the q  Figure 4: Comparison between the D 0 v 2 values measured in the unbiased sample and in the two event-shape classes obtained using TPC and V0A to compute q 2 , for the 10-30% (top row) and 30-50% (bottom row) centrality classes. Only statistical uncertainties are shown.

Systematic uncertainties
The values of v 2 are affected by systematic uncertainties related to (i) the signal extraction from the invariant-mass distributions, (ii) the correction for the beauty feed-down contribution, (iii) the presence of non-flow effects, and (iv) the centrality dependence of the event-plane resolution correction R 2 .
The uncertainty on the D-meson raw yield extraction from the invariant-mass distributions of candidates in the in-plane and out-of-plane azimuthal angle intervals was estimated with a multi-trial approach by repeating the fits several times with different configurations. In particular, the lower and upper limits of the fit range and the background fit function were varied, while the Gaussian width was kept fixed to the one extracted from the fits to the invariant mass distributions integrated over q 2 and ∆ϕ. For each fit configuration, the D-meson v 2 was calculated from the in-plane and out-of-plane yields. The absolute systematic uncertainties were assigned as the r.m.s. of the v 2 distribution resulting from the different fits. They range from 0.005 to 0.040 in the 30-50% centrality class and from 0.008 to 0.040 in the 10-30% centrality class, depending on the p T interval and the D-meson species. Further checks on the stability of the results were performed by repeating the procedure leaving the Gaussian width as a free parameter in the fits and by using a bin-counting method for the definition of the raw yield. With the latter method, the signal yield was obtained by counting the histogram entries in the invariant-mass region of the signal (|M − M peak | < 3.5σ ), after subtracting the background contribution estimated from a fit to the side bands (|M − M peak | > 4σ ). The v 2 values obtained from these checks were found to be within the uncertainty estimated by varying the fit conditions and therefore no additional systematic uncertainty was assigned. For the analysis with ESE selection, further studies were carried out by comparing the output of the multiple-trial fit procedure described above in the small-q 2 , large-q 2 and q 2 -integrated samples for each of the tested fit configurations. These checks indicated that this contribution to the systematic uncertainty is uncorrelated between the event samples selected based on the q 2 value. ESE for the D-meson v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration The contribution of the beauty feed-down correction to the systematic uncertainty was estimated varying (i) the quark mass and the renormalisation and factorisation scales in the FONLL calculations; (ii) the R feed-down AA hypothesis; and (iii) the v feed-down 2 hypothesis as described in Sec. 2. The value of the corresponding absolute systematic uncertainty ranges from 0.001 to 0.030 depending on the D-meson species and p T as well as on the ESE-selected class.
The systematic uncertainty on the event-plane resolution correction factor R 2 has two contributions, which are common to the unbiased, small-q 2 and large-q 2 samples. The first one originates from possible non-flow effects affecting the estimation of R 2 , when the particles reconstructed in the two semivolumes of the TPC are used as sub-events. It was estimated by comparing the value of R 2 obtained by introducing two different pseudorapidity gaps (∆η = 0.2 and ∆η = 0.4) between the sub-events of the TPC tracks with positive/negative η. The second contribution is due to the centrality dependence of R 2 within the classes used in the analysis. The central value of R 2 was computed from the three subevent correlations averaged over the events in the 10-30% and 30-50% intervals. The uncertainty was estimated by comparing this value with those obtained as weighted averages of the R 2 values in narrow centrality intervals, using as weights either the D-meson yields or the average number of nucleonnucleon collisions. A systematic uncertainty of 2% on R 2 was assigned based on these studies for all centrality and ESE-selected classes.
For the ESE-selected samples, an additional bias on the resolution correction factor can originate from autocorrelations because of the usage of TPC tracks (V0A signals) for both q TPC 2 (q V0A 2 ) and R 2 determination. In particular, the selection on q TPC 2 can bias the correlation between the sub-events of charged particles reconstructed in the TPC with 0 < η < 0.8 and with −0.8 < η < 0 used in the three sub-event calculation of R 2 . To estimate this systematic uncertainty, an alternative approach to compute R 2 was utilised, which is based on (i) the correlations between the sub-events reconstructed with the V0 and half of the TPC tracks (with η < 0) and (ii) the assumption that the ratio of the variables χ V0 and χ TPC,η<0 governing the event plane resolution (see Ref. [35] for its definition) is the same in the unbiased and ESE-selected samples. The difference between the R 2 values obtained with this approach and the three sub-event method, which amounts to 3% and 5% in the 10-30% and 30-50% centrality classes, respectively, was assigned as systematic uncertainty on the ESE-selected samples. The same procedure was adopted for the samples selected using the q V0A 2 . In this case, the systematic uncertainty was estimated to be of the order of 1% for the large-q V0A 2 sample, while negligible for the small-q V0A 2 sample, for both the centrality classes.
As discussed in Sec. 2.1, a further bias in the analyses with q TPC 2 -based selection could be induced by non-flow correlations between the D meson and the sample of tracks used for the q 2 measurement, which can include charged particles originating from the fragmentation of the charm quarks. To further study this effect, the analysis with q TPC 2 -based selection was repeated introducing a "jet-veto" pseudorapidity gap of |∆η| = 0.1 units between each D-meson candidate and the tracks used to measure q TPC 2 . Since no significant difference was observed, no systematic uncertainty was assigned.

Results
In Fig. 5 the elliptic flow coefficient v 2 of prompt D 0 , D + , and D * + mesons is reported as a function of p T in the centrality class 10-30%. The symbols are positioned at the average p T of the reconstructed D mesons, which is determined as the average of the p T distribution of candidates in the signal invariantmass region, after subtracting the contribution of the background candidates estimated from the side bands. The systematic uncertainty of the feed-down correction is displayed separately in the figure. The v 2 of D 0 , D + and D * + mesons is consistent among the various species and larger than zero in the interval 2 < p T < 8 GeV/c. inverse of the squared absolute statistical uncertainties as weights and is reported in the left panel of Fig. 6. The systematic uncertainties were propagated by considering the contribution from the eventplane resolution R 2 and the feed-down correction as correlated among the D-meson species. In the right panel of Fig. 6, the average v 2 of D 0 , D + , and D * + as a function of p T in the centrality class 30-50% taken from [48] is reported. The measurements in both centrality classes are compatible within uncertainties with the D 0 -meson v 2 measured with the Scalar Product (SP) method [69,78] by the CMS Collaboration in |y| < 1 [49]. The charged-pion v 2 measured in |y| < 0.5 by the ALICE Collaboration using the SP method [91] is also superimposed for comparison. The D-meson v 2 is similar in magnitude to that of π ± for 4 < p T < 10 GeV/c. In the region p T < 4 GeV/c, where a mass ordering for light hadrons is observed and described by hydrodynamical calculations [92], the values of the D-meson v 2 are slightly lower than those of π ± , but compatible within uncertainties. Figure 7 shows the prompt D 0 and D + v 2 as a function of p T in the small-q TPC 2 and large-q TPC 2 samples, in the centrality classes 10-30% (top row) and 30-50% (bottom row). The measurement of the D * + v 2 in the ESE-selected samples was not possible due to the small statistical significance, while the measurements of D 0 and D + mesons were performed in wider p T intervals compared to the unbiased v 2 measurement and in the reduced range 2 < p T < 12 GeV/c, due to the limited size of the data sample. The measurements of the v 2 of the two different D-meson species in the ESE-selected classes are compatible with each other within uncertainties. Also in this case, the symbols are positioned at the average D-meson p T determined as described above.
The average v 2 of D 0 and D + mesons has been calculated in the small-q TPC 2 and large-q TPC 2 samples with the same weighted average procedure described above. It is shown for the two considered centrality classes in the top panels of Fig. 8 together with the v 2 measured in the unbiased sample, recalculated in the same p T intervals of the ESE analysis. In the bottom panels of the same figure, the ratio of the average √ s NN = 5.02 TeV ALICE Collaboration (see text for details) and unbiased samples, in the 10-30% (left) and 30-50% (right) centrality classes. Bottom panels: ratios of the measured v 2 in the ESE-selected classes to the one obtained from the unbiased sample. D-meson v 2 from the ESE-selected samples with respect to that of the unbiased samples is illustrated. The statistical uncertainties on the ratio were propagated taking into account the degree of correlation between the measured yields in the small-q TPC 2 (large-q TPC 2 ) and the unbiased sample. The systematic uncertainties were propagated considering the contribution from the centrality dependence and the nonflow contaminations among sub-events of R 2 as well as the feed-down correction as correlated between the measurements in the ESE-selected and the unbiased samples.
The observation of a flat ratio as a function of p T for light hadron v 2 with ESE-selection at √ s NN = 2.76 TeV indicated that the q 2 value is connected to a global property of the event [62]. For D mesons, the modification of the v 2 in the q TPC 2 -selected samples is compatible within uncertainties with a flat behaviour as a function of p T for both the 10-30% and the 30-50% centrality classes. However, the current precision of the measurement does not allow to exclude a p T dependence which would indicate the presence of non-flow contaminations.
Selecting the 20% (60%) highest (lowest) q TPC 2 sample leads to a change of about 40% (25%) in the measured v 2 . The corresponding variation of the average q TPC 2 in the ESE-selected classes was found to be about 65% and 30% in the large-q TPC  (see text for details) and unbiased samples, in the 10-30% (left) and 30-50% (right) centrality classes. The charged-particle v 2 obtained at the same energy, centrality classes and ESE samples are superimposed for comparison. Bottom panels: ratios of the measured v 2 in the ESE-selected classes to the one obtained from the unbiased sample.
(decrease) of the D-meson v 2 and the average q TPC 2 observed in the large-q TPC 2 (small-q TPC 2 ) sample with respect to the unbiased one is similar within uncertainties in the two centrality intervals considered. Considering as null hypothesis v 2 (large-q TPC 2 ) = v 2 (small-q TPC 2 ), the probability to observe the measured positive ∆v 2 = v 2 (large-q TPC 2 ) − v 2 (small-q TPC 2 ) in the full p T range of the measurement, corresponds to a significance of about 4σ , taking into account both statistical and systematic uncertainties in each centrality class. It is however important to keep in mind that part of the observed effect could be slightly enlarged by non-flow contaminations, as previously mentioned.
The effect of the ESE selection on the D-meson v 2 was compared to that observed for charged particles. For this comparison, the ESE selection was performed using q V0A 2 , in order to avoid autocorrelations and non-flow contaminations. In the top panels of Fig. 9, the average D 0 and D + v 2 in the ESE-selected and unbiased samples in the 10-30% (left panel) and 30-50% (right panel) centrality classes are depicted together with the charged-particle v 2 measured at the same energy, centrality classes and ESE-selected samples. The charged-particle v 2 was measured with the SP method considering reconstructed tracks with |η| < 0.8 and 0.2 < p T < 12 GeV/c, selected as in Ref. [91]. The   figure show the ratios of the v 2 measured in the ESE-selected samples with respect to the unbiased one. The ratios between the charged-particle v 2 show almost no p T dependence, confirming that the usage of the q V0A 2 provides a selection of a global property of the collision. The relative variation of the charged-particle v 2 in the large-q V0A 2 and small-q V0A 2 samples was found to be of about 14-15% and 7-8%, respectively. These values reflect the reduced sensitivity of the ESE selection obtained using the V0A with respect to that based on TPC tracks. The ratios of the average D-meson v 2 in the ESE-selected samples with respect to the unbiased one were found to be compatible within uncertainties with those of charged particles in the corresponding samples, suggesting that the response to the ESE selection is similar for D mesons and the bulk of light hadrons. However, given the reduced selectivity of q V0A 2 and the current experimental uncertainties, the ratios of the average D-meson v 2 are also compatible with unity, and therefore a firm conclusion cannot be drawn. Nevertheless, the comparison between D mesons and charged particles will be crucial for future larger data samples, to better asses the magnitude of the correlation between the D-meson and the soft-hadron v 2 .
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 yields of prompt D 0 and D +  The D-meson raw yields integrated over ∆ϕ were extracted from the fits to the invariant-mass distributions in the ESE-selected and unbiased classes and normalised to the corresponding number of events in the considered sample. As described in Sec. 2, the selection and reconstruction efficiencies of prompt D mesons do not show any dependence on q 2 within the ESE selections considered in this analysis, therefore no correction to the raw yields was applied. The fraction of prompt D mesons, f prompt , was estimated using the same strategy adopted for the v 2 measurement and it is the same in the ESE-selected and the unbiased samples.
The ratio of the D-meson yields in the small-q TPC 2 (large-q TPC 2 ) sample to those in the unbiased sample are shown in Fig. 10 as a function of p T in the 10-30% (top row) and 30-50% (bottom row) centrality classes. The systematic uncertainty on the raw D-meson yield extraction was evaluated directly on the ratio of the yields, applying the same strategy used for the v 2 (see Sec. 3). The systematic uncertainty on the reconstruction and selection efficiency, arising from a possible imperfect description of the data in the Monte Carlo simulations, cancels out in the ratio, since the efficiency is the same in the two ESE-selected classes.
The average of the ratio of the D 0 and D + yields in the small-q TPC 2 (large-q TPC 2 ) sample to those in the unbiased sample is depicted in Fig. 11. It was computed by using the inverse of the squared relative statistical uncertainties as weights.
In the 10-30% centrality class, the ratio between the D-meson yields in ESE-selected samples to those in the unbiased sample was found to be compatible with unity in the measured p T range. In the 30-50% centrality class, the central values of the D-meson per-event yields in the large-q TPC 2 (small-q TPC 2 ) samples were found to be higher (lower) than those in the unbiased sample in all the measured p T intervals in the range 3 < p T < 12 GeV/c. However, the ratios between the yields in the ESE-selected samples to the unbiased yields are compatible with unity within about one standard deviation.
In the light-hadron sector, the effect induced by the correlation between radial and elliptic flow, attributed to a larger initial density in more anisotropic events, was observed to be of the order of 5% for charged ESE for the D-meson v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV ALICE Collaboration pions with p T ≈ 4 GeV/c in mid-central Pb-Pb collisions at √ s NN = 2.76 TeV [62]. Since the ratio between the D-meson yields in ESE-selected samples to those in the unbiased sample was found to be compatible with unity, a possible similar effect is expected to be smaller than the present experimental uncertainties, which do not allow for any conclusion.

Summary
The first application of the event-shape engineering technique to the measurement of D-meson production in Pb-Pb collisions at √ s NN = 5.02 TeV has been presented.
The elliptic flow of D 0 , D + , and D * + mesons at mid-rapidity in the 10-30% (30-50%) centrality class was measured with the event-plane technique and found to be larger than zero in the transverse momentum interval 2 < p T < 8(10) GeV/c and similar in magnitude to that of charged pions for p T > 4 GeV/c, while slightly lower for p T < 4 GeV/c, in the same centrality class.
The v 2 coefficient of D 0 and D + mesons was measured in events with different magnitude of the average bulk elliptic flow, quantified by the value of q 2 measured using TPC tracks to maximise the selectivity. The observation of a larger (smaller) D-meson v 2 in events with large-q TPC 2 (small-q TPC 2 ) values confirms a correlation between D-meson azimuthal anisotropy and the collective expansion of the bulk of light hadrons. When using the V0A to measure q 2 in order to reduce non-flow contaminations and autocorrelations, the variation of the D-meson v 2 in the small-q V0A 2 and large-q V0A 2 samples was found to be compatible within uncertainties with that of charged particles, suggesting a similar response to the ESE selection.
The ratio of the p T -differential yields measured in the ESE-selected samples with respect to those in the unbiased sample was found to be compatible with unity in both the 10-30% and 30-50% centrality classes, with a possible indication of larger D-meson yield for 3 < p T < 12 GeV/c in events with higherthan-average bulk elliptic flow in the 30-50% centrality class. With the current uncertainties no firm conclusion can be drawn on the possible interplay between the initial spatial anisotropy and the charmquark energy loss and radial flow.
The measurements presented in this paper open the way to the study of heavy-quark production with the Event-Shape Engineering technique, which offers a new possibility to understand the correlation of heavy-quark and bulk properties. An improved precision is expected to be achieved with future data samples that will be collected in 2018 and during Run 3 and 4 of the LHC [93,94].