Azimuthal anisotropy of D meson production in Pb-Pb collisions at $\sqrt{s_{\rm NN}} = 2.76$ TeV

The production of the prompt charmed mesons $D^0$, $D^+$ and $D^{*+}$ relative to the reaction plane was measured in Pb-Pb collisions at a centre-of-mass energy per nucleon-nucleon collision of $\sqrt{s_{\rm NN}} = 2.76$ TeV with the ALICE detector at the LHC. D mesons were reconstructed via their hadronic decays at central rapidity in the transverse momentum ($p_{\rm T}$) interval of 2-16 GeV/$c$. The azimuthal anisotropy is quantified in terms of the second coefficient $v_2$ in a Fourier expansion of the D meson azimuthal distribution, and in terms of the nuclear modification factor $R_{\rm AA}$, measured in the direction of the reaction plane and orthogonal to it. The $v_2$ coefficient was measured with three different methods and in three centrality classes in the interval 0-50%. A positive $v_2$ is observed in mid-central collisions (30-50% centrality class), with an mean value of $0.204_{-0.036}^{+0.099}$ (tot.unc.) in the interval $2<p_{\rm T}<6$ GeV/$c$, which decreases towards more central collisions (10-30% and 0-10% classes). The positive $v_2$ is also reflected in the nuclear modification factor, which shows a stronger suppression in the direction orthogonal to the reaction plane for mid-central collisions. The measurements are compared to theoretical calculations of charm quark transport and energy loss in high-density strongly-interacting matter at high temperature. The models that include substantial elastic interactions with an expanding medium provide a good description of the observed anisotropy. However, they are challenged to simultaneously describe the strong suppression of high-$p_{\rm T}$ yield of D mesons in central collisions and their azimuthal anisotropy in non-central collisions.


Introduction
Collisions of heavy nuclei at ultra-relativistic energies are expected to lead to the formation of a highdensity colour-deconfined state of strongly-interacting matter.According to calculations of Quantum Chromo-Dynamics (QCD) on the lattice (see e.g.[1][2][3][4]), a phase transition to the Quark-Gluon Plasma (QGP) state can occur in these collisions, when conditions of high energy density and temperature are reached.Heavy quarks (charm and beauty), with large masses m c ≈ 1.3 and m b ≈ 4.5 GeV/c 2 , are produced in pairs predominantly at the initial stage of the collision [5] in hard scattering processes characterized by timescales shorter than the medium formation time.They traverse the medium and interact with its constituents via both inelastic (medium-induced gluon radiation, i.e. radiative energy loss) [6,7] and elastic (collisional) [8] QCD processes.Heavy-flavour hadrons are thus effective probes of the properties of the medium formed in the collisions.
Compelling evidence for heavy-quark energy loss in strongly-interacting matter is provided by the observation of a modification of the transverse momentum (p T ) distributions of heavy-flavour hadrons.This modification is quantified by the nuclear modification factor R AA (p T ) = dN AA /dp T T AA dσ pp /dp T , where dN AA /dp T is the differential yield in nucleus-nucleus collisions in a given centrality class, dσ pp /dp T is the cross section in pp collisions, and T AA is the average nuclear overlap function [9].In central nucleus-nucleus collisions at RHIC and LHC energies, R AA values significantly below unity were observed for heavy-flavour hadrons with p T values larger than a few GeV/c [10][11][12][13][14][15].A suppression by a factor up to 3-5 (R AA ≈ 0.25) at p T ≃ 5 GeV/c was measured in central collisions for inclusive electrons and muons from heavy-flavour hadron decays, both at RHIC ( √ s NN = 200 GeV), by the PHENIX and STAR Collaborations [10,11], and at the LHC ( √ s NN = 2.76 TeV), by the ALICE Collaboration [14].
At the LHC, the effect was also measured separately for charm, via D mesons by the ALICE Collaboration [13], and for beauty, via non-prompt J/ψ particles from B hadron decays by the CMS Collaboration [15].
The D meson suppression at RHIC and at the LHC is described (see [12,13]) by model calculations that implement a combination of mechanisms of heavy-quark interactions with the medium, via radiative and collisional processes, as well as in-medium formation and dissociation of charmed hadrons [16][17][18][19][20][21][22].
Model comparisons with more differential measurements can provide important insights into the relevance of the various interaction mechanisms and the properties of the medium.In particular, the dependence of the partonic energy loss on the in-medium path length is expected to be different for each mechanism (linear for collisional processes [8] and close to quadratic for radiative processes [7]).In addition, it is an open question whether low-momentum heavy quarks participate, through interactions with the medium, in the collective expansion of the system and whether they can reach thermal equilibrium with the medium constituents [23,24].It was also suggested that low-momentum heavy quarks could hadronize not only via fragmentation in the vacuum, but also via the mechanism of recombination with other quarks from the medium [24,25].
These questions can be addressed with azimuthal anisotropy measurements of heavy-flavour hadron production with respect to the reaction plane, defined by the beam axis and the impact parameter of the collision.For non-central collisions, the two nuclei overlap in an approximately lenticular region, the short axis of which lies in the reaction plane.Hard partons are produced at an early stage, when the geometrical anisotropy is not yet reduced by the system expansion.Therefore, partons emitted in the direction of the reaction plane (in-plane) have, on average, a shorter in-medium path length than partons emitted orthogonally (out-of-plane), leading a priori to a stronger high-p T suppression in the latter case.
In the low-momentum region, the in-medium interactions can also modify the parton emission directions, thus translating the initial spatial anisotropy into a momentum anisotropy of the final-state particles.Both effects cause a momentum anisotropy that can be characterized with the coefficients v n and the symmetry planes Ψ n of the Fourier expansion of the p T -dependent particle distribution d 2 N/dp T dϕ in azimuthal angle ϕ.The elliptic flow is the second Fourier coefficient v 2 , which can also be expressed as the average over all particles in all events of the angular correlation cos[2(ϕ − Ψ 2 )].The symmetry planes Ψ n for all harmonics would coincide with the reaction plane if nuclei were spherically symmetric with a matter density depending only on the distance from the centre of the nucleus.Due to fluctuations in the positions of the participant nucleons, the plane of symmetry fluctuates event-by-event around the reaction plane, independently for each harmonic, so that the Ψ n directions no longer coincide.
A path-length dependent energy loss, which gives a positive v 2 , is considered to be the dominant contribution to the azimuthal anisotropy of charged hadrons in the high p T region, above 8-10 GeV/c [29,30].At low p T , a large v 2 is considered as an evidence for the collective hydrodynamical expansion of the medium [31,32].Measurements of light-flavour hadron v 2 over a large p T range at RHIC and LHC are generally consistent with these expectations [18,[33][34][35][36][37][38][39].In contrast to light quarks and gluons, which can be produced or annihilated during the entire evolution of the medium, heavy quarks are produced predominantly in initial hard scattering processes and their annihilation rate is small [5].Thus, the final state heavy-flavour hadrons at all transverse momenta originate from heavy quarks that experienced each stage of the system evolution.High-momentum heavy quarks quenched by in-medium energy loss are shifted towards low momenta and, while participating in the collective expansion, they may ultimately thermalize in the system.In this context, the measurement of D meson v 2 is also important for the interpretation of recent results on J/ψ anisotropy [26], because J/ψ mesons formed from cc recombination would inherit the azimuthal anisotropy of their constituent quarks [27,28].
An azimuthal anisotropy in heavy-flavour production was observed in Au-Au collisions at RHIC with v 2 values of up to about 0.13 for electrons from heavy-flavour decays [40].The measured asymmetry is reproduced by several models [19][20][21][41][42][43][44][45][46] implementing heavy-quark transport within a medium that undergoes a hydrodynamical expansion.The transport properties, i.e. the diffusion coefficients, of heavy quarks in the medium can be related to its shear viscosity [41].For LHC energies these models predict a large v 2 (in the range 0.10-0.20 in semi-central collisions) for D mesons at p T ≈ 2-3 GeV/c and a decrease to a constant value v 2 ≈ 0.05 at high p T .The models described in Refs.[20,[43][44][45][46] include, at the hadronization stage, a contribution from the recombination of charm quarks with light quarks from the medium, which enhances v 2 at low p T .
The measurement of the D meson v 2 in the centrality class 30-50% in Pb-Pb collisions at √ s NN = 2.76 TeV, carried out using the ALICE detector, was presented in [47].The v 2 coefficient was found to be significantly larger than zero in the interval 2 < p T < 6 GeV/c and comparable in magnitude with that of charged particles.
Here the measurement is extended to other centrality classes and accompanied with a study of the azimuthal dependence of the nuclear modification factor with respect to the reaction plane.The decays D 0 → K − π + , D + → K − π + π + and D * + → D 0 π + and charge conjugates were reconstructed.The v 2 coefficient was measured with various methods in the centrality class 30-50% as a function of p T .For the D 0 meson, which has the largest statistical significance, the centrality dependence of v 2 in the range 0-50% is presented and the anisotropy is also quantified in terms of the nuclear modification factor R AA in two 90 • -wide azimuthal intervals centred around the in-plane and out-of-plane directions.
The experimental apparatus is presented in Section 2. The data analysis is described in Section 3, including the data sample, the D meson reconstruction and the anisotropy measurement methods.Systematic uncertainties are discussed in Section 4. The results on v 2 and R AA are presented in Section 5 and compared with model calculations in Section 6.

Experimental apparatus
The ALICE apparatus is described in [48].In this section, the characteristics of the detectors used for the D meson analyses are summarized.The z-axis of the ALICE coordinate system is defined by the beam direction, the x-axis lies in the horizontal plane and is pointing towards the centre of the LHC accelerator ring and the y-axis is pointing upward.
Charged-particle tracks are reconstructed in the central pseudo-rapidity1 region (|η| < 0.9) with the Time Projection Chamber (TPC) and the Inner Tracking System (ITS).For this analysis, charged hadron identification was performed using information from the TPC and the Time Of Flight (TOF) detectors.These detectors are located inside a large solenoidal magnet that provides a field with a strength of 0.5 T, parallel to the beam direction.Two VZERO scintillator detectors, located in the forward and backward pseudo-rapidity regions, are used for online event triggering, collision centrality determination and, along with the Zero Degree Calorimeter (ZDC), for offline event selection.
The ITS [49] includes six cylindrical layers of silicon detectors surrounding the beam vacuum tube, at radial distances from the nominal beam line ranging from 3.9 cm for the innermost layer to 43 cm for the outermost one.The two innermost layers consist of Silicon Pixel Detectors (SPD) with a pixel size of 50 × 425 µm 2 (rϕ × z, in cylindrical coordinates), providing an intrinsic spatial resolution of 12 µm in rϕ and 100 µm in z.The third and fourth layers use Silicon Drift Detectors (SDD) with an intrinsic spatial resolution of 35 µm and 25 µm in rϕ and z, respectively.The two outermost layers of the ITS contain double-sided Silicon Strip Detectors (SSD) with an intrinsic spatial resolution of 20 µm in rϕ and 830 µm in the z-direction.The alignment of the ITS sensor modules is crucial for the precise space point recontruction needed for the heavy-flavour analyses.It was performed using survey information, cosmic-ray tracks and pp data.A detailed description of the employed methods can be found in [49].
The TPC [51] covers the pseudo-rapidity interval |η| < 0.9 and extends in radius from 85 cm to 247 cm.Charged-particle tracks are reconstructed and identified with up to 159 space points.The transverse momentum resolution for tracks reconstructed with the TPC and the ITS ranges from about 1% at p T = 1 GeV/c to about 2% at 10 GeV/c, both in pp and Pb-Pb collisions.The TPC also provides a measurement of the specific energy deposition dE/dx, with up to 159 samples.The truncated mean method, using only the lowest 60% of the measured dE/dx samples, gives a Gaussian distribution with a resolution (ratio of sigma over centroid) of about 6%, which is slightly dependent on the track quality and on the detector occupancy.
The TOF detector [52] is positioned at a radius of 370-399 cm and it has the same pseudo-rapidity coverage as the TPC (|η| < 0.9).The TOF provides an arrival time measurement for charged tracks with an overall resolution, including the measurement of the event start time, of about 80 ps for pions and kaons at p T = 1 GeV/c in the Pb-Pb collision centrality range used in this analysis [52].
The VZERO detector [53] consists of two arrays of scintillator counters covering the pseudo-rapidity regions −3.7 < η < −1.7 (VZERO-C) and 2.8 < η < 5.1 (VZERO-A).Each array is composed of 8 × 4 segments in the azimuthal and radial directions, respectively.This detector provides a low-bias interaction trigger (see Section 3.1).For Pb-Pb collisions, the signal amplitude from its segments is used to classify events according to centrality, while the azimuthal segmentation allows for an estimation of the reaction plane.
The ZDCs are located on either side of the interaction point at z ≈ ±114 m.The timing information from the neutron ZDCs was used to reject parasitic collisions between one of the two beams and residual nuclei present in the vacuum tube.
Table 1: Number of events and integrated luminosity for the considered centrality classes, expressed as percentiles of the hadronic cross section.The uncertainty on the integrated luminosity derives from the uncertainty of the hadronic Pb-Pb cross section from the Glauber model [9,54].

Data sample and event selection
The analysis was performed on a data sample of Pb-Pb collisions recorded in November and December 2011 at a centre-of-mass energy per nucleon-nucleon collision of √ s NN = 2.76 TeV.The events were collected with an interaction trigger based on information from the VZERO detector, which required coincident signals recorded in the detectors at forward and backward pseudo-rapidities.An online selection based on the VZERO signal amplitude was used to enhance the sample of central and midcentral collisions through two separate trigger classes.Events were further selected offline to remove background coming from parasitic beam interactions by using the time information provided by the VZERO and the neutron ZDC detectors.Only events with a reconstructed interaction point (primary vertex), determined by extrapolating charged-particle tracks, within ±10 cm from the centre of the detector along the beam line were used in the analysis.
Collisions were classified in centrality classes, determined from the sum of the amplitudes of the signals in the VZERO detector and defined in terms of percentiles of the total hadronic Pb-Pb cross section.In order to relate the centrality classes to the collision geometry, the distribution of the VZERO summed amplitudes was fitted by a model based on the Glauber approach for the geometrical description of the nuclear collision [9] complemented by a two-component model for particle production [54].The centrality classes used in the analysis are reported in Table 1, together with the number of events in each class and the corresponding integrated luminosity.
D 0 and D + candidates were defined from pairs and triplets of tracks within the fiducial acceptance |η| < 0.8, selected by requiring at least 70 associated space points in the TPC, χ 2 /ndf < 2 for the momentum fit, and at least two associated hits in the ITS, with at least one of them in the SPD.A transverse momentum threshold p T > 0.4 GeV/c was applied in order to reduce the combinatorial background.D * + candidates were obtained by combining the D 0 candidates with tracks selected with the same requirements as described above, but with a lower transverse momentum threshold p T > 0.1 GeV/c and at least three associated hits in the ITS, with at least one of them in the SPD.The lower p T threshold was used because the momentum of the pions from D * + decays is typically low, as a consequence of the small mass difference between D * + and D 0 .
The selection of tracks with |η| < 0.8 introduces a steep drop in the acceptance of D mesons for rapidities larger than 0.7-0.8,depending on p T .A fiducial acceptance region was, therefore, defined as: |y| < y fid (p T ), with y fid (p T ) increasing from 0.7 to 0.8 in 2 < p T < 5 GeV/c and taking a constant value of 0.8 for p T > 5 GeV/c [13].The D meson v 2 results are not expected to be affected by this small variation in rapidity acceptance.
The D meson yields were measured with an invariant mass analysis of reconstructed decays, using kinematic and geometrical selection criteria, and particle identification (PID).The selection of D 0 and D + decays was based on the reconstruction of secondary vertices with a separation of a few hundred microns from primary vertex.In the case of the D * + decay, the secondary vertex of the produced D 0 was reconstructed.The coordinates of the primary vertex and of the secondary vertices, as well as the corresponding covariance matrices, were computed using a χ 2 minimization method [56].
The selection strategy is the same as in previous pp [56,57] and Pb-Pb [13] analyses.It exploits the displacement of the decay tracks from the primary vertex (transverse impact parameter, d 0 ), the separation between the secondary and primary vertices (decay length, L) and the pointing of the reconstructed meson momentum to the primary vertex.
The transverse impact parameter d 0 of a given track is defined as the signed distance of closest approach of the extrapolated track to the primary vertex in the (x, y) plane.The sign of d 0 is attributed based on the position of the primary vertex with respect to the curve of the (x, y) projection of the track.In Pb-Pb collisions, the impact parameter resolution in the transverse direction is better than 65 µm for tracks with a transverse momentum larger than 1 GeV/c and reaches 20 µm for p T > 20 GeV/c [13].This includes the contribution from the primary vertex precision, which is better than 10 µm in the central and semi-central Pb-Pb events used in this analysis.The impact parameter measurement is significantly less precise along the longitudinal direction, e.g.170 µm at p T = 1 GeV/c.
A pointing condition was applied via a selection on the angle ϑ pointing between the direction of the reconstructed momentum of the candidate and the straight line connecting the primary and secondary vertices.For Pb-Pb collisions, two additional selection variables were introduced with respect to pp analyses, namely the projection of the pointing angle and of the decay length onto the transverse plane (ϑ xy pointing and L xy ).The selection requirements were tuned so as to provide a large statistical significance for the signal and to keep the selection efficiency as high as possible.The chosen selection values depend on the p T of the D meson and become more stringent from peripheral to central collisions.
The selection criteria for the centrality class 30-50% are described in the following.The D 0 candidates were selected by requiring the decay tracks to have an impact parameter significance |d 0 |/σ d 0 > 0.5 (σ d 0 is the uncertainty on the track impact parameter), and to form a secondary vertex with a track-to-track distance of closest approach smaller than 250-300 µm, depending on p T , and a decay length larger than 100 µm.The product of the decay track impact parameters, which are of opposite sign for welldisplaced signal topologies, was required to be below −(200 µm) 2 at low p T (2-3 GeV/c) and below −(120 µm) 2 for high p T candidates (12)(13)(14)(15)(16) GeV/c), with a smooth variation between these values in 2-12 GeV/c.A significance of the projection of the decay length in the transverse plane L xy /σ L xy (where σ L xy is the uncertainty on L xy ) larger than 5 was also required.A selection on the angle ϑ * between the kaon momentum in the D 0 rest frame and the boost direction was used to reduce the contamination from background candidates that do not represent real two-body decays and typically have large values of | cos ϑ * |.The selection | cos ϑ * | < 0.8 was applied.The pointing of the D 0 momentum to the primary vertex was implemented by requiring cos ϑ pointing > 0.95 and cos ϑ xy pointing > 0.998 at low p T (2-3 GeV/c).Since the background is lower at high p T , the cuts were progressively made less stringent for increasing p T .In the 0-10% and 10-30% centrality classes the combinatorial background is larger than in 30-50%.Therefore, the selections were made more stringent and they are similar to those used for the 0-20% centrality class in [13].
The D + candidates were selected by requiring a decay length larger than 1200-1600 µm, depending on p T , and cos ϑ pointing larger than 0.998 (0.990) in the p T interval 3-4 (8-12) GeV/c, with a smooth variation in-between.Further requirements to reduce the combinatorial background were cos ϑ xy pointing > 0.993-0.998and L xy /σ L xy > 9-11, depending on the candidate p T .In general, the D + selection criteria are more stringent than those of the D 0 because of the larger combinatorial background.
In the D * + analysis, the selection of the decay D 0 candidates was similar to that used for the D 0 analysis.Only D 0 candidates with invariant mass within 2.5 σ of the PDG mass value [55] were used, where σ is the p T -dependent Gaussian sigma of the invariant mass distribution observed in data.The decay pion was selected with the same track quality criteria as for the D 0 and D + decay tracks.
Pions and kaons were identified with the TPC and TOF detectors, on the basis of the difference, expressed in units of the resolution (σ ), between the measured signal and that expected for the considered particle species.Compatibility regions at ±3 σ around the expected mean energy deposition dE/dx and timeof-flight were used.Tracks without a TOF signal were identified using only the TPC information.This particle identification strategy provided a reduction by a factor of about three of the combinatorial background in the low-p T range, while preserving most of the signal (see Section 3.4).
The D 0 and D + raw yields were obtained with a fit to the invariant mass M distribution of the D meson candidates.For the D * + signal the mass difference ∆M = M(K − π + π + ) − M(K − π + ) was considered.The fit function is the sum of a Gaussian to describe the signal and a term describing the background, which is an exponential for D 0 and D + and has the form f (∆M) = a (∆M − m π ) b for the D * + , where m π is the charged pion mass and a and b are free parameters.The centroids and the widths of the Gaussian functions were found to be in agreement, respectively, with the D meson PDG mass values [55] and with the simulation results, confirming that the background fluctuations were not causing a distortion in the signal line shape.An example of invariant mass distributions will be shown in Section 3.3.

Azimuthal anisotropy analysis methods
The p T -differential azimuthal distribution of produced particles can be described by a Fourier series: where Ψ n is the initial state spatial plane of symmetry of the n-th harmonic, defined by the geometrical distribution of the nucleons participating in the collision.In order to determine the second harmonic coefficient v 2 , the Q vector is defined from the azimuthal distribution of charged particles, where ϕ i are the azimuthal angles and N is the multiplicity of charged particles.The weights w i are discussed later in the text.The charged particles used for the Q vector determination are indicated in the following as reference particles (RFP).
The azimuthal angle of the Q vector is called event plane angle and it is an estimate of the second harmonic symmetry plane Ψ 2 .
The charged particle tracks used for the Q vector determination were selected with the following criteria: at least 50 associated space points in the TPC; χ 2 /ndf < 2 for the momentum fit in the TPC; a distance of closest approach to the primary vertex smaller than 3.2 cm in z and 2.4 cm in the (x, y) plane.In order to minimize the non-uniformities in the azimuthal acceptance, no requirement was applied on the number of ITS points associated to the track.To avoid auto-correlations between the D meson candidates and the event plane angles, the Q vector was calculated for each candidate excluding from the set of reference particles the tracks used to form that particular candidate.Tracks with p T > 150 MeV/c were considered and the pseudo-rapidity interval was limited to the positive region 0 < η < 0.8, where the TPC acceptance and efficiency were more uniform as function of the azimuthal angle for this data set.The remaining azimuthal non-uniformity was corrected for using weights w i in Eq. ( 2), defined as the inverse of the ϕ distribution of charged particles used for the Q vector determination, 1/(dN/dϕ i ), multiplied by a function f . This function mimics the p T -dependence of the charged particle v 2 and it improves the estimate of Ψ 2 by enhancing the contribution of particles with a stronger flow signal (see e.g.Ref. [36]).The distribution of the event plane angle ψ 2 obtained for this set of reference particles is shown in Fig. 1 (a), for the centrality range 30-50%.The distribution, divided by its integral, exhibits a residual non-uniformity below 1%.
An additional study was performed with the Q vector determined from the azimuthal distribution of signals in the segments of the VZERO detectors, which are sensitive to particles produced at forward and backward rapidities.The Q vector was calculated with Eq. ( 2), with the sum running over the eight azimuthal sectors of each VZERO detector, where ϕ i was defined by the central azimuth of the i-th sector, and w i equal to the signal amplitude in the i-th sector for the selected event, which is proportional to the number of charged particles crossing the sector.Non-uniformities in the VZERO acceptance and efficiency were corrected for using the procedure described in [61].The residual non-uniformity is about 1%, as shown in Fig. 1 (a).
For the event plane method, the measured anisotropy v obs 2 was divided by the event plane resolution correction factor R 2 according to the equation v 2 = v obs 2 /R 2 , with R 2 being smaller than one.This resolution depends on the multiplicity and v 2 of the RFP [58].For the event plane computed using TPC tracks, R 2 was determined from the correlation of the event plane angles reconstructed from RFP in the two sides of the TPC, −0.8 < η < 0 and 0 < η < 0.8, i.e. two samples of tracks (called sub-events) with similar multiplicity and v 2 .R 2 is shown in Fig. 1 (b) as a function of collision centrality.The average R 2 values in the three centrality classes used in this analysis are 0.6953 (0-10%), 0.8503 (10-30%) and 0.8059 (30-50%).The statistical uncertainty on R 2 is negligible (∼ 10 −4 ).The systematic uncertainty on R 2 was estimated by using the three-sub-event method described in [62].In this case, the event planes reconstructed in the TPC (0 < η < 0.8), VZERO-A (2.8 < η < 5.1) and VZERO-C (−3.7 < η < −1.7) were used.This method yielded R 2 values smaller than those obtained from the two-sub-events method by 6.9%, 2.0% and 2.3% for the centrality classes 0-10%, 10-30% and 30-50%.A part of this difference can be attributed to the presence of short-range non-flow correlations that are suppressed when the three sub-events with a pseudo-rapidity gap are used.Non-flow correlations can originate from resonance or cascade-like decays and from jets.The resolution of the event plane determined from the VZERO detector (summing the signals in VZERO-A and VZERO-C) is also shown in Fig. 1 (b).In this case, R 2 was measured with three sub-events, namely the signals in the VZERO detector (both A and C sides) and the tracks in the positive and negative η regions of the TPC.The systematic uncertainty was estimated from the difference with the results obtained with two TPC sub-events separated by 0.4 units in pseudorapidity (η gap).The event plane determination has a poorer resolution with the VZERO detector than with the TPC tracks.As a consequence, the v 2 measurement is expected to be more precise with the TPC event plane.
In the event plane method, the D meson yield was measured in two 90 • -wide intervals of ∆ϕ = ϕ D − ψ 2 : in-plane (− π 4 < ∆ϕ ≤ π 4 and 3π 4 < ∆ϕ ≤ 5π 4 ) and out-of-plane ( π 4 < ∆ϕ ≤ 3π 4 and 5π 4 < ∆ϕ ≤ 7π 4 ).ϕ D is defined as the azimuthal angle of the D meson momentum vector at the primary vertex.The invariant mass distributions for the three meson species are shown in Fig. 2 in three p T intervals for the 30-50% centrality class, along with the fits used for the yield estimation (Section 3.2).When fitting the invariant mass distribution in the two ∆ϕ intervals, the centroid and the width of the Gaussian functions were fixed, for each meson species and for each p T interval, to those obtained from a fit to the invariant mass distribution integrated over ϕ, where the statistical significance of the signal is larger.Integrating Eq. ( 1) and including the correction for the event plane resolution 1/R 2 yields: The contribution of higher harmonics to the v 2 value calculated with this equation can be evaluated by integrating the corresponding terms of the Fourier series.All odd harmonics, as well as v 4 and v 8 , induce the same average contribution to N in-plane and N out-of-plane due to symmetry, and therefore they do not affect v 2 calculated with Eq. ( 4).The contribution of v 6 , v 10 and higher harmonics is assumed to be negligible based on the values measured for light-flavour hadrons [63,64].
The measurement of the elliptic flow with the scalar product method is given by [58]: where indicates an average over D meson candidates in all events.The vector u is defined as u = (cos 2ϕ D , sin 2ϕ D ), where ϕ D the D meson candidate momentum azimuthal direction.The Q a,b and u a,b vectors were computed from charged particles and D meson candidates, respectively, in two separate pseudo-rapidity regions: a) 0 < η < 0.8 and b) −0.8 < η < 0. The elliptic flow was computed by correlating D mesons from the positive η region with the charged particles in the negative η region, and vice versa.This separation in pseudo-rapidity suppresses two-particle correlations at short distance that are due to decays (D * → D + X and B → D ( * ) + X ).The denominator in Eq. ( 5) plays a similar role as the resolution correction in the event plane method.Since the resolution is proportional to the number of used RFP, the vectors Q a and Q b were normalized by N a and N b , respectively, before averaging over all events.The azimuthal non-uniformity of the TPC response, which results in non-zero average values of Q a and Q b , was corrected for using a re-centering procedure [58]:   The two-particle cumulant is defined by the equation [60,65,66]: For this method, the azimuthal non-uniformity of the detector acceptance and efficiency was corrected for with the aforementioned re-centering procedure.In contrast to the scalar product method, there is no pseudo-rapidity gap between the D mesons and the RFP for the two-particle cumulant method.
For both the scalar product and two-particle cumulant methods, the v 2 of D meson candidates was computed in narrow intervals of invariant mass M for D 0 and D + and mass difference ∆M for the D * + .
In each invariant mass interval, the measured v 2 is the weighted average of the D meson v 2 (v S 2 ) and the background v 2 (v B 2 ) with the weights given by the relative fractions of signal (S) and background (B) in  that interval.In order to extract the values of v S 2 and v B 2 , a simultaneous fit of the distributions of counts and v 2 as a function of invariant mass M was performed.The invariant mass distribution was fitted with a sum of two terms for signal and background, as explained in Section 3.2.The v 2 (M) distribution was fitted with a function: The background contribution v B 2 was parametrized by a linear function of M.An example of the corresponding distributions and fits is shown in Fig. 3 for D 0 mesons in the interval 4 < p T < 6 GeV/c with the two-particle cumulants method (a) and D * + mesons in the interval 2 < p T < 4 GeV/c with the scalar product method (b).The values of v S 2 , hereafter indicated as v 2 {2} and v 2 {SP}, are also reported in the figure .Since the measured D meson yield has a feed-down contribution from B meson decays, the measured v 2 is a combination of v 2 of promptly produced and feed-down D mesons.In fact, the contribution of D mesons from B meson decays is enhanced by the applied selection criteria, because the decay vertices of the feed-down D mesons are, on average, more displaced from the primary vertex.The elliptic flow of promptly produced D mesons, v prompt 2 , can be obtained from the measured v all 2 (v 2 {EP}, v 2 {2} or v 2 {SP}) as: where f prompt is the fraction of promptly produced D mesons in the measured raw yield and v feed-down 2 is the elliptic flow of D mesons from B decays, which depends on the dynamics of beauty quarks in the medium.These two quantities have not been measured.According to Eq. ( 8), the value of v all 2 is independent of f prompt and equal to v . The central value of the prompt D meson elliptic flow was defined under this assumption, which removes the need to apply the feed-down correction.Because of the larger mass of the b quark, the v 2 of B mesons is expected to be lower than that of D mesons.Therefore, the choice of v feed-down 2 = v prompt 2 as central value is the most conservative for the observation of D meson v 2 > 0. The details of the systematic uncertainty related to this assumption are discussed in Section 4.
Azimuthal anisotropy of charm production in Pb-Pb collisions at √ s NN = 2.76 TeV 13

Azimuthal dependence of the nuclear modification factor
The in-plane and out-of-plane nuclear modification factors of prompt D 0 mesons are defined as: where dN in (out) AA /dp T are the D 0 meson per-event yields, integrated over the two 90 • -wide intervals used to determine v 2 with the event plane method.The factor 2 in Eq. ( 9) accounts for the fact that the D meson yields for Pb-Pb collisions are integrated over half of the full azimuth.R in (out) AA was measured in the 30-50% centrality class for D 0 mesons, which have the highest signal significance, using the yields relative to the event plane defined with TPC tracks in 0 < η < 0.8.The average value of the nuclear overlap function in this centrality class, T AA = 3.87 ± 0.18 mb −1 , was determined with the procedure described in [54].
The yields of prompt D 0 mesons in the two azimuthal intervals were obtained as: The raw yields N D 0 +D 0 raw were divided by a factor of two to obtain the charge (particle and antiparticle) averaged yields.The factor c refl (p T ) was introduced to correct the raw yields for the contribution of signal candidates that are present in the invariant mass distribution both as D 0 → K − π + and as D 0 → π − K + (the combination with wrong mass hypothesis assignment is called 'reflection').To correct for the contribution of B meson decay feed-down, the raw yields were multiplied by the prompt fraction f prompt , whose determination is described later in this section.Furthermore, they were divided by the product of prompt D meson acceptance and efficiency (Acc × ε) prompt , normalized by the decay channel branching ratio (BR), the transverse momentum (∆p T ) and rapidity (∆y = 2 y fid ) interval widths and the number of events (N events ).The normalization by ∆y gives the corrected yields in one unit of rapidity |y| < 0.5.
The (Acc × ε) correction was determined, as a function of p T , using Monte Carlo simulations with a detailed description of the ALICE detector geometry and the GEANT3 particle transport package [69].The simulation was tuned to reproduce the (time-dependent) position and width of the interaction vertex distribution, as well as the number of active electronic channels and the accuracy of the detector calibration.The HIJING v1.383 [68] generator was used to simulate Pb-Pb collisions at √ s NN = 2.76 TeV and all the produced particles were transported through the detector simulation.Prompt and feed-down D meson signals were added using pp events from the PYTHIA v6.4.21 [67] event generator with the Perugia-0 tune [70].Each simulated pp event contained a cc or bb pair with D mesons decaying into the hadronic channels of interest for the analysis.Out of all the particles produced in these PYTHIA pp events, only the heavy-flavour decay products were kept and transported through the detector simulation together with the particles produced by HIJING.In order to minimize the bias on the detector occupancy, the number of D mesons injected into each HIJING event was adjusted according to the Pb-Pb collision centrality.
The efficiencies were evaluated from simulated events that had the same average charged-particle multiplicity, corresponding to the same detector occupancy, as observed for real events in the centrality class 30-50%.Figure 4 shows (Acc × ε) for prompt and feed-down D 0 mesons within the rapidity interval |y| < y fid .The magnitude of (Acc × ε) increases with p T , starting from about 1% and reaching about 10-15% at high p T .Also shown in Fig. 4 are the (Acc × ε) values for the case where no PID was applied.
The relative difference with respect to the (Acc × ε) obtained using also the PID selection is only about 5%, thus illustrating the high efficiency of the applied PID criteria.The (Acc × ε) for D mesons from B decays is larger than for prompt D mesons by a factor of about 1.5, because the decay vertices of the The values for prompt (solid lines) and feed-down (dotted lines) D 0 mesons are shown.Also displayed, for comparison, are the values for prompt D 0 mesons without PID selection (dashed lines).The lower panel shows the ratio of the efficiencies for prompt D 0 mesons in the in-plane and out-of-plane regions used for the analysis.This ratio was estimated using simulation samples with a difference in particle multiplicity similar to that observed in data for the two azimuthal regions.
feed-down D mesons are more displaced from the primary vertex and are, therefore, more efficiently selected by the analysis cuts.
The possible difference in the reconstruction and selection efficiency between in-plane and out-of-plane D 0 mesons was studied using simulations.This difference could arise from the variation of the particle density, and consequently of the detector occupancy, induced by the azimuthal anisotropy of bulk particle production.The difference in occupancy was estimated in data using the multiplicity of SPD tracklets in the two considered azimuthal intervals.Tracklets are defined as combinations of two hits in the two SPD layers that are required to point to the primary vertex.They can be used to measure the multiplicity of charged particles with p T > 50 MeV/c and |η| < 1.6.The SPD tracklet multiplicity in the 30-50% centrality class was found to be larger in-plane than out-of-plane by about 12%.In order to study the efficiency variation, two sets of simulated events with 12% difference in average multiplicity were used.The ratio of the two efficiencies was found to be consistent with unity (see lower panel of Fig. 4) and therefore no correction was applied.
The correction factor c refl for the contribution of reflections to the raw yield was determined by including in the invariant mass fit procedure a template of the distribution of reflected signal candidates, which was obtained from the simulation for each p T interval.This distribution has a centroid close to the D 0 mass and has typical r.m.s.values of about 100 MeV/c 2 , i.e. about one order of magnitude larger than the signal invariant mass resolution.The distribution from the simulation was parametrized with the sum of two Gaussians, in order to remove the statistical fluctuations.In the fit with the template, the ratio of the integrals of the total distribution of reflections and of the Gaussian used for the signal were fixed to the value obtained from the simulation.This ratio is mostly determined by the PID selection, which limits the probability that a true K − π + pair can be also compatible with the π − K + mass hypothesis.For the v 2 analysis described in the previous section, the PID selection was used only for tracks with p < 4 GeV/c.Since the contribution of the reflections does not depend on the angle relative to the event plane, it is not necessary to apply the c refl correction for v 2 .For the R AA analysis, in order to minimize the correction, the PID selection was extended to tracks with p > 4 GeV/c, requiring the compatibility of the TOF and TPC signals with the expectations for kaons and pions within 3 σ .It was verified that this change results in a variation of v 2 well within the uncertainties.The correction factor c refl was determined as the ratio of the signal yield from the fit including the reflections template and the signal yield from the fit without the template.It was computed using the sum of the in-plane and out-of-plane invariant mass distributions, in order to have a more precise value, and it was applied as in Eq. ( 10) for both the in-plane and out-of-plane yields.The procedure was validated using the simulation, where the signal yield obtained from the fit with the template can be compared with the true signal yield.The numerical value of c refl ranges from 0.98 in the interval 2 < p T < 3 GeV/c to 0.90 in the interval 4 < p T < 16 GeV/c.Figure 5 shows an example of the fits without (a) and with (b) template for the interval 4-6 GeV/c.The fraction f prompt of promptly produced D mesons in the measured raw yields was obtained, following the procedure introduced in [13], as: In this expression, where the symbol of the p T -dependence has been omitted for brevity, N D 0 raw is the measured raw yield (corrected by the c refl factor) and N D 0 feed-down raw is the contribution of D 0 mesons from B decays to the raw yield, estimated on the basis of the FONLL calculation of beauty production [71].In detail, the B meson production cross section in pp collisions at √ s = 2.76 TeV was folded with the B → D 0 + X decay kinematics using EvtGen [72] and multiplied by: the average nuclear overlap function T AA in the 30-50% centrality class, the acceptance-times-efficiency for feed-down D 0 mesons, and the other factors introduced in Eq. (10).In addition, the nuclear modification factor R feed-down AA of D mesons from B decays was accounted for.The comparison of the R AA of prompt D mesons [73] with that of J/ψ from B decays [74] measured in the CMS experiment indicates that charmed hadrons are more suppressed than beauty hadrons.Therefore, it was assumed that the ratio of the nuclear modification factors for feed-down and prompt D mesons lies in the range 1 was used to compute the correction, and the variation over the full range, which also accounts for possible centrality and p T dependences, was used to assign a systematic uncertainty.The hypothesis on the nuclear modification of feed-down D mesons was changed with respect to the assumption used in [13], based on the most recent results on the R AA of prompt D meson and nonprompt J/ψ mentioned above.As it was done for the v 2 measurement, the feed-down contribution was computed assuming v feed-down . Therefore, the ratio R feed-down AA /R prompt AA is the same in-plane and out-of-plane.The systematic uncertainty related to this assumption is discussed in Section 4. For D 0 mesons, assuming R feed-down AA = 2 • R prompt AA , the resulting f prompt ranges from about 0.80 in the lowest transverse momentum interval (2 < p T < 3 GeV/c) to about 0.75 at high p T .
The D 0 yields in the two azimuthal regions with respect to the event plane, obtained from Eq. ( 10), were corrected for the event plane resolution using the correction factor R 2 (Section 3.3) and the relation given in Eq. (4).For example, the correction factor for the in-plane 2 )/2, where N in (out) is the D 0 raw yield.The value R 2 = 0.8059 ± 0.0001 for the 30-50% centrality class and the typical N in /N out magnitude result in a correction of approximately +4 (−6)% for the inplane (out-of-plane) yields.
The prompt D 0 meson production cross section in pp collisions used in the calculation of the nuclear modification factor was obtained by scaling the p T -differential cross section in |y| < 0.5 at √ s = 7 TeV, measured using a data sample of L int = 5 nb −1 [56].The p T -dependent scaling factor was defined as the ratio of the cross sections obtained from FONLL calculations [71] at √ s = 2.76 and 7 TeV [75].The scaled D 0 meson p T -differential cross section is consistent with that measured at √ s = 2.76 TeV using a smaller statistics data sample with L int = 1.1 nb −1 [57], which only covered a reduced p T interval with a statistical uncertainty of 20-25% and was therefore not used as pp reference.The correction for reflections was not applied for the D 0 cross section in pp collisions.It was verified that the resulting signal bias is smaller than 5% (c refl > 0.95), which is less than the systematic uncertainty assigned for the yield extraction (10-20% [56]).

Systematic uncertainties
Several sources of systematic uncertainty were considered for both v 2 and R AA analyses.The uncertainties on v 2 are described first.Afterwards, the systematic uncertainties affecting R AA in-plane and out-of-plane are discussed.The uncertainties for the 30-50% centrality class are summarized in Tables 2 and 3.In the following, the quoted uncertainties are symmetric around the central value of the measurement, unless the upper and lower parts are reported separately.

Uncertainties on v 2
One of the main sources of uncertainty originates from the D meson yield extraction using a fit to the invariant mass distributions.This uncertainty was estimated by repeating the fits under different conditions and by utilizing alternative methods for the yield determination.For the v 2 analysis with the event plane method, the fit ranges and the functional forms for the combinatorial background were varied.Polynomial and exponential functions were tried for D 0 and D + background shapes, while a threshold function multiplied by an exponential was considered for the D * + : a , with a and b as free parameters.The D meson yield was also extracted by counting the entries in the invariant mass distributions after background subtraction.For this procedure the background was estimated with a fit to the left and right sides of the D meson invariant mass peak (side-band regions), using the fit functions described in Section 3.2.The v 2 analysis employing the event plane method was performed by fixing the Gaussian centroids and widths of the in-plane and out-of-plane invariant mass distributions to the values obtained from a fit of the ϕ-integrated distribution.The analysis was repeated with free Gaussian parameters in the fit.The systematic uncertainty due to the yield measurement was estimated as the maximum variation of the v 2 values obtained from the described tests.It amounts to 10-20% for the D 0 meson, depending on the p T and centrality intervals, and 20-50% for the D + and D * + mesons, depending on the p T interval.The same procedure was applied for the two-particle correlation methods (scalar product and two-particle cumulants), except for the bin counting method and the fixed Gaussian centroids and widths.Instead, the parametrization of the background v B 2 (M) was varied from a first order to a second order polynomial.The resulting uncertainty is in the range 15-30%.
For the event plane method, two alternative procedures were considered to extract v 2 , which are not directly based on the measurement of the signal yields from the invariant mass distribution.These procedures use the distribution of cos(2∆ϕ) versus invariant mass (where ∆ϕ = ϕ D − ψ 2 ) and the relation v 2 = cos(2∆ϕ) .In the first procedure, the distribution of cos(2∆ϕ) is considered for the signal region (|M − m D | < 3 σ ) and the two side-band regions (4 < |M − m D | < 7 σ ).The distribution of cos(2∆ϕ) for the background is obtained by averaging, bin-by-bin, the distributions of cos(2∆ϕ) in the two side bands.This background distribution is then rescaled to the integral of the background fit function in the invariant mass signal region and it is subtracted from the total cos(2∆ϕ) distribution in the signal region.In this way, the distribution of cos(2∆ϕ) of the signal is obtained.Its mean value gives the D meson v 2 .In the second procedure, a distribution of cos(2∆ϕ) as a function of invariant mass is used for a simultaneous fit of the v 2 and the yield, as in the case of the two-particle correlation methods.These two alternative procedures result in D meson v 2 values that are consistent with those obtained from the event plane method with two ∆ϕ bins.Therefore, no systematic uncertainty is taken for the v 2 extraction procedure.
The v 2 analysis was repeated with different sets of cuts for the selection of D meson candidates.A set of tighter and a set of looser cuts with respect to those described in Section 3.2 were considered for each D meson species, thus varying the signal yield by about 30-50% and, consequently, the significance and the signal-to-background ratio.The resulting v 2 values were found to be consistent within statistical uncertainties.Consequently, this contribution to the systematic uncertainty was neglected.
The uncertainty due to the event plane resolution was estimated with the two and three sub-event methods with an η gap.The three sub-events were defined using the TPC tracks and the signals in the two VZERO detectors.The resolutions estimated with these two methods differ by 6.9%, 2.0% and 2.3% in the 0-10%, 10-30% and 30-50% centrality classes, respectively (see Fig. 1 (b)).A symmetric systematic uncertainty equal to the relative difference between R 2 values obtained with the two and three sub-event methods was assigned to the D meson v 2 .
The uncertainty due to the centrality dependence of the event plane resolution was estimated from the difference between two ways to define the average resolution in the centrality classes used in the analysis, starting from the resolutions in fine centrality intervals (see Fig. 1 (b)).Namely, a plain arithmetic average and an average weighted with the D meson yield measured in smaller centrality classes (2.5% wide).The latter was estimated using D 0 meson raw yields in wide p T intervals and the sum of the two ∆ϕ intervals, in order to reduce the statistical fluctuations.The difference between these averages was found to be about 2%, 0.5% and 2% for the 0-10%, 10-30% and 30-50% centrality classes, respectively.The resulting total uncertainties on R 2 amount to 7%, 2% and 3% for the three centrality classes.
The distribution of collision impact parameters selected in a given centrality class slightly depends on the pseudo-rapidity coverage of the detector used for the centrality determination.The analysis was repeated using the number of tracks in the TPC as a centrality estimator, instead of the total signal measured in the VZERO detector.A relative systematic uncertainty of 10% was assigned to the v 2 values measured with the scalar product and two-particle cumulant methods, on the basis of the difference of the resulting v 2 values.This difference could originate from the dependence of the RFP multiplicity fluctuations on the centrality estimator.No significant difference was observed for the event plane method when using the TPC, instead of the VZERO, for the centrality determination.
As explained at the end of section 3.3, the central value of the prompt D meson v 2 was obtained without applying a correction for the feed-down from B meson decays, on the basis of the assumption (see Eq. ( 8)).The systematic uncertainty associated with this assumption was estimated by varying it in the interval 0 . This range covers all model predictions for v 2 of charm and beauty hadrons [20,21,43].The lower limit of the variation range, v feed-down 2 = 0, gives v prompt 2 = v all 2 / f prompt .The f prompt values for each of the D meson species and each p T interval were obtained using FONLL calculations [71] (see section 3.4).Under the assumption R feed-down AA = 2 • R prompt AA , the f prompt values change from 0.8 to 0.75 (0.85 to 0.8) from low to high p T for D 0 (D + and D * + ) mesons (the feed-down contribution is larger for D 0 mesons because of the stronger constraint on the separation between the secondary and the primary vertex).A set of f prompt values was computed by varying the heavy-quark masses and the perturbative scales in the FONLL calculation as prescribed in [71], and the ratio R feed-down AA /R prompt AA in the range 1 < R feed-down AA /R prompt AA < 3. The smallest value of f prompt was used to assign the uncertainty related to the B feed-down contribution to the elliptic flow of prompt D mesons.The maximum relative uncertainty is about +45 −0 %.

Uncertainties on R AA
For the analysis of the D 0 meson R AA in-plane and out-of-plane, the same sources of systematic uncertainty as for the v 2 measurement with the event plane method were considered.Additional systematic uncertainties, which are specific to the R AA measurement, stem from the tracking, selection and particle identification efficiencies, and from the uncertainty of the proton-proton reference yield.The evaluation of these uncertainties is similar as in [13] and it is described in the following.
In order to reduce the statistical fluctuations, the uncertainty of the D 0 yield extraction was estimated using the ϕ-integrated invariant mass distributions.The fit procedure was varied, as described for the v 2 analysis.The resulting uncertainty is 7% for 2 < p T < 8 GeV/c and 10% for 8 < p T < 16 GeV/c.The systematic uncertainty on the correction factor for signal reflections, c refl , was estimated by changing by ±50% the ratio of the integral of the reflections over the integral of the signal obtained from the simulation and used in the invariant mass fit with the reflections template.In addition, the shape of reflections invariant mass distribution template was varied using a polynomial parametrization of the distribution from the simulation, instead of a double-Gaussian parametrization.These variations resulted in an uncertainty of 1-2% for 2 < p T < 4 GeV/c and of 5% for 4 < p T < 16 GeV/c on the c refl factor.
The systematic uncertainty of the tracking efficiency was estimated by comparing the probability to match the TPC tracks extrapolated to the ITS hits in data and simulation, and by varying the track quality selection criteria (for example, the minimum number of associated hits in the TPC and in the ITS and maximum χ 2 /ndf of the momentum fit).The efficiency of the track matching and the association of hits in the silicon pixel layers was found to be described by the simulation with maximal deviations on the level of 5% in the p T range relevant for this analysis (0.5-15 GeV/c).The effect of misassociating ITS hits to tracks was studied using simulations.It was found that the fraction of D mesons with at least one decay track with a wrong hit associated increases with centrality, due to the higher detector occupancy, and vanishes at high p T , where the track extrapolation between ITS layers is more precise.In the centrality class 30-50%, this fraction is about 2% in the transverse momentum interval 2 < p T < 16 GeV/c.It was verified that the signal selection efficiencies are the same for D mesons with and without wrong hit associations.The total systematic uncertainty of the track reconstruction procedure amounts to 5% for single tracks, which results in a 10% uncertainty for D 0 mesons (two-track final state).
The uncertainty of the correction for the selection on the decay topology was evaluated by repeating the analysis with different sets of cuts and was defined as the variation of the resulting corrected yields with respect to the value corresponding to the baseline cuts.This resulted in a variation up to 10% in the p T intervals used in the analysis.The analysis was repeated without applying the PID selection and the resulting corrected yields were found to be consistent within 5% with those obtained with the PID selection.Therefore, a systematic uncertainty of 5% was assigned for the PID efficiency correction in the simulation.
The uncertainty of the efficiencies arising from the difference between the real and simulated D meson momentum distributions depends on the width of the p T intervals and on the variation of the efficiencies within them.This uncertainty includes also the effect of the p T dependence of the nuclear modification factor.The mean efficiency in a given p T interval was computed by re-weighting the simulated D 0 meson yield according to the p T distribution measured for D 0 mesons in central Pb-Pb collisions [13].The systematic uncertainty was defined as the difference with respect to the efficiency computed using the p T distribution from a FONLL calculation [71] multiplied by the R AA value from one of the models [21] that closely describe the central value of the measurement (see Section 6).This uncertainty is of 2% in the interval 2 < p T < 3 GeV/c, where the efficiency increases steeply with p T , and below 1% for The uncertainty of 3% on the event plane resolution correction factor R 2 in the 30-50% centrality class was propagated to the R AA observables, resulting in an uncertainty in the range 0.5-2%, depending on the p T interval.
The systematic uncertainty due to the subtraction of feed-down D mesons from B meson decays was estimated following the procedure described in [13].The contribution of the uncertainties inherent in the FONLL perturbative calculation was included by varying the heavy-quark masses and the factorization and renormalization scales in the ranges proposed in [71].This contribution partly cancels in the R AA ratio, because these variations are done simultaneously for the Pb-Pb yield and for the pp reference cross section.The uncertainty introduced by the hypothesis on the value of the feed-down D meson R AA was estimated from the variation 1 < R feed-down AA /R prompt AA < 3. The total uncertainty due to the feed-down correction, which is common to the in-plane and out-of-plane R AA , ranges between +9 −13 % at low p T and +14 −12 % at high p T .The hypothesis on the value of v 2 for D mesons from B decays, that was varied in the range 0 , introduces an additional contribution to the systematic uncertainty, which is anti-correlated between R in-plane AA and R out-of-plane AA .This uncertainty is typically of +5 −0 % for in-plane and +0 −5 % for out-of-plane.The uncertainty of the pp reference used for the calculation of R AA has two contributions.The first is due to the systematic uncertainty of the measured D 0 meson p T -differential yield at √ s = 7 TeV and it is about 17%, approximately constant with p T [56].The second contribution is due to the scaling to √ s = 2.76 TeV.It ranges from +31 −10 % at low p T to about 5% at high p T [13].The uncertainties on the pp cross section normalization (3.5%) [56] and the average nuclear overlap function T AA (4.7% for the class 30-50%) were also included.The contribution due to the 1.1% relative uncertainty on the fraction of the hadronic cross section used in the Glauber fit to determine the centrality classes [54] was obtained by estimating the variation of the D meson dN/dp T when the limits of the centrality classes are shifted by ±1.1% (e.g.instead of 30-50%, 30.3-50.6% and 29.7-49.5%)[13].
The resulting uncertainty, common to all p T intervals, is 2% for the 30-50% centrality class.The total normalization uncertainty, computed taking the quadratic sum of these three contributions, is 6.2%.
The systematic uncertainties of R AA were grouped in three categories, depending on their correlation between the in-plane and out-of-plane measurements.The uncorrelated systematic uncertainties affect the two R AA independently; this category includes only the yield extraction uncertainty.The correlated systematic uncertainties affect the two R AA in the same way and do not affect their relative difference.The uncertainties on the correction efficiencies (for track reconstruction, selection cuts, particle identification and D 0 p T distribution in the simulation) and on the correction for reflections, as well as those on the pp reference, the variation of pQCD scales and the R feed-down AA hypothesis used for the feed-down subtraction are included in this category.Another correlated uncertainty is due to the normalization ( T AA and centrality class definition), which is quoted separately.The anti-correlated systematic uncertainties could shift the two R AA in opposite directions, affecting their difference.This category includes the contribution from the unknown azimuthal anisotropy of feed-down D mesons (variation of v feed-down 2 ) and the contribution from the event plane resolution correction factor.Within each category, the uncertainties from different sources were added in quadrature.with the event plane (from Ref. [47]), scalar product and two-particle cumulant methods (columns).For the first method, the event plane was estimated with TPC tracks in 0 < η < 0.8; for the other methods, TPC tracks in −0.8 < η < 0.8 were used as RFP.The symbols are positioned at the average p T measured within each interval.

Elliptic flow
The elliptic flow v 2 measured with the event plane method is shown as a function of p T in the left column of Fig. 6 for D 0 , D + and D * + mesons in the 30-50% centrality class.The event plane was estimated from TPC tracks in the range 0 < η < 0.8.The symbols are positioned horizontally at the average p T of reconstructed D mesons.This value was determined as the average of the p T distribution of candidates in the signal invariant mass region, after subtracting the contribution of the background candidates, which was estimated from the side bands.This average p T of the reconstructed D mesons is larger than that of the produced D mesons, because the efficiency increases with increasing p T (see Fig. 4).The vertical error bars represent the statistical uncertainty, the open boxes are the systematic uncertainties from the anisotropy determination and the event plane resolution, and the filled boxes are the uncertainties due to the B feed-down contribution.The elliptic flow of the three D meson species is consistent within statistical uncertainties and ranges between 0.1 and 0.3 in the interval 2 < p T < 8 GeV/c.For p T > 12 GeV/c, v 2 is consistent with zero within the large statistical uncertainties.The central and right-most panels of the same figure show the v 2 results obtained with the scalar product and two-particle cumulant methods, respectively.The results from the three methods are consistent within statistical uncertainties for the three meson species.
Figure 7 shows the v 2 of the D 0 mesons measured with the event plane (a) and scalar product (b) methods using reference particles from the TPC detector (i.e. in a η range that overlaps with the D meson acceptance) or from the VZERO detectors at −3.7 < η < −1.7 and 2.8 < η < 5.1 (i.e. with a large η gap with respect to the D mesons).The agreement between the results with and without η gap indicates that the bias due to non-flow correlations is within the statistical precision of the measurement.
For the 30-50% centrality class an average v 2 of D 0 , D + and D * + was already computed in [47] from the event plane method results, using the statistical uncertainties as weights.The resulting D meson v 2 has a value 0.204 ± 0.030 (stat) ± 0.020 (syst) +0.092 −0 (B feed-down), averaged over the p T intervals 2-3, 3-4, 4-6 GeV/c.This value is larger than zero with a significance, calculated from the combined statistical and systematic uncertainties, of 5.7 σ .
Figure 8 shows the D 0 meson v 2 in the three centrality classes 0-10%, 10-30% and 30-50% as a function of p T .The D 0 meson v 2 is compared with that of charged particles [36], for the same centrality classes.D meson and charged particle results are obtained with the event plane method using TPC and the VZERO detectors, respectively.The magnitude of v 2 is similar for charmed hadrons and light-flavour hadrons, which dominate the charged-particle sample.
The centrality dependence of the D 0 elliptic flow is shown in Fig. 9 for three transverse momentum intervals in the range 2 < p T < 6 GeV/c.A decreasing trend of v 2 towards more central collisions is observed, as expected because of the decreasing initial geometrical anisotropy.

Nuclear modification factor in and out of the event plane
The nuclear modification factors of D 0 mesons in the 30-50% centrality class are shown in Fig. 10 for the in-plane and out-of-plane directions with respect to the event plane.The event plane was estimated with TPC tracks in 0 < η < 0.8.The error bars represent the statistical uncertainties, which are to a large extent independent for the two azimuthal intervals, since they are dominated by the statistical uncertainties of the Pb-Pb data.The uncorrelated (empty boxes), correlated (brackets) and anti-correlated (shaded boxes) systematic uncertainties are shown separately.The normalization uncertainty, shown as a box at R AA = 1, is common to both measurements.A large suppression is observed in both directions with respect to the event plane for p T > 4 GeV/c.At lower transverse momentum, the suppression appears to be reduced, especially in the in-plane direction, where R AA reaches unity at a p T of 2-3 GeV/c.Overall, a stronger suppression in the out-of-plane  is equivalent to the observation of v 2 > 0 (as shown in the top-left panel of Fig. 6), since Eq. ( 4) can be expressed also as

Comparison with model calculations
A number of theoretical model calculations are available for the elliptic flow coefficient v 2 and the nuclear modification factor R AA of heavy-flavour hadrons.Figure 11 shows a comprehensive comparison of these models to measurements of the R AA of D 0 mesons in-plane and out-of-plane in the 30-50% centrality class, of the average R AA of D 0 , D + and D * + in the 0-20% centrality class [13], and of the v 2 averaged over the D meson species in the centrality class 30-50% [47].
The following models are considered and compared to data: I WHDG [18].This is a perturbative QCD calculation of parton energy loss, including both radiative (DGLV [76]) and collisional processes.A realistic collision geometry based on the Glauber model [9] is used, without hydrodynamical expansion, so that the anisotropy results only from path-length dependent energy loss.Hadronization is performed using vacuum fragmentation functions.The medium density is constrained on the basis of the π 0 R AA in central collisions at √ s NN = 200 GeV and scaled to LHC energy according to the increase of the charged-particle multiplicity.The model describes well the D meson R AA in the centrality interval 0-20% (slightly overestimating the suppression, as it does also for charged particles [13]), and gives an almost p T -independent v 2 ≈ 0.06, which is smaller than the measured values in the range 2 < p T < 6 GeV/c.Consequently, the difference between the in-plane and out-of-plane R AA suppression is underestimated: the model describes well the out-of-plane R AA and lies below the in-plane R AA .
II MC@sHQ+EPOS, Coll+Rad(LPM) [77].This pQCD model includes collisional and radiative (with Landau-Pomeranchuk-Migdal correction [78]) energy loss mechanisms for heavy quarks with running strong coupling constant.The medium fluid dynamical expansion is based on the EPOS model [79].A component of recombination of heavy quarks with light-flavour quarks from the QGP is also incorporated in the model.This model yields a substantial anisotropy (v 2 ≈ 0.12-0.08 from low to high p T ), which is close to that observed in data.The nuclear modification factor is substantially overestimated below p T ≈ 5 GeV/c and correctly described at higher p T .
only.The heavy-quark transport coefficient is calculated within a non-perturbative T -matrix approach, where the interactions proceed via resonance formation that transfers momentum from the heavy quarks to the medium constituents.The model includes hydrodynamic medium evolution, constrained by light-flavour hadron spectra and elliptic flow data, and a component of recombination of heavy quarks with light-flavour quarks from the QGP.Diffusion of heavy-flavour hadrons in the hadronic phase is also included.The model provides a good description of the observed suppression of D mesons over the entire p T range.The maximum anisotropy, v 2 of about 0.13 at 2 < p T < 4 GeV/c, is close to that observed in the data.Towards larger p T , the model tends to underestimate v 2 , as well as the difference of the in-plane and out-of-plane R AA .
IV POWLANG [19].This transport model is based on collisional processes treated within the framework of Langevin dynamics, within an expanding deconfined medium described by relativistic viscous hydrodynamics.The transport coefficients entering into the relativistic Langevin equation are evaluated by matching the hard-thermal-loop calculation of soft collisions with a perturbative QCD calculation for hard scatterings.Hadronization is implemented via vacuum fragmentation functions.This model overestimates the high-p T suppression, it yields a value for v 2 significantly smaller than observed in data and also underestimates the difference between the in-plane and out-of-plane suppression.
V BAMPS [21].This partonic transport model is based on the Boltzmann approach to multiparton scattering.Heavy quarks interact with the medium via collisional processes computed with running strong coupling constant.Hadronization is performed using vacuum fragmentation functions.The lack of radiative processes is accounted for by scaling the binary cross section with a correction factor, which is tuned to describe the heavy-flavour decay electron elliptic flow and nuclear modification factor at RHIC.When applied to calculations for LHC energy, this correction factor results in an underestimation of the D meson R AA for p T > 5 GeV/c and a large azimuthal anisotropy, with v 2 values up to 0.20, similar to those observed in the data.The nuclear modification factors in-plane and out-of-plane are well described up to 5 GeV/c, while for higher p T the in-plane R AA is underestimated.
VI UrQMD [45].The Langevin approach for the transport of heavy quarks is in this case implemented within the UrQMD model [80].This model includes a realistic description of the medium evolution by combining hadronic transport and ideal hydrodynamics.The transport of heavy quarks is calculated on the basis of a resonance model with a decoupling temperature of 130 MeV.Hadronization via quark coalescence is included.The calculation parameters are tuned to reproduce the heavy-flavour measurements at RHIC ( √ s NN = 200 GeV) and kept unchanged for calculations at the LHC energy.The model describes the measured D meson v 2 , as well as R AA in the interval 4 < p T < 8 GeV/c, but it fails to reproduce the significant suppression measured for R AA at p T of 2-3 GeV/c.VII Cao, Qin, Bass [46].This model is also based on the Langevin approach.In addition to quasielastic scatterings, radiative energy loss is incorporated by treating gluon radiation as an additional force term.The space-time evolution of the medium is modelled using a viscous hydrodynamic simulation.The hadronization of heavy quarks has a contribution based on the recombination mechanism.With respect to [46], the curves shown in the figure were obtained with a more recent parametrization for the nuclear shadowing of the parton distribution functions.This model provides a good description of the R AA data in central collisions, but it yields a value of v 2 significantly smaller than the measured one (similarly to the WHDG and POWLANG models) and also underestimates the difference between the in-plane and out-of-plane suppression.
Overall, the anisotropy is qualitatively described by the models that include both charm quark energy loss in a geometrically anisotropic medium and mechanisms that transfer to charm quarks the elliptic flow induced during the system expansion.These mechanisms include collisional processes (MC@sHQ+EPOS, Coll+Rad(LPM) [77], BAMPS [21]) and resonance scattering with hadronization via recombination (TAMU elastic [44], UrQMD [45]) in a hydrodynamically expanding QGP.Models that do not include a collective expansion of the medium or lack a contribution to the hadronization of charm quarks from recombination with light quarks from the medium predict in general a smaller anisotropy than observed in the data.The comparison for R AA and v 2 shows that it is challenging to simultaneously describe the large suppression of D mesons in central collisions and their anisotropy in non-central collisions.In general, the models that are best in describing R AA tend to underestimate v 2 and the models that describe v 2 tend to underestimate the measured R AA at high p T .It is also worth noting that most of the calculations do reproduce the RHIC measurements of heavy-flavour decay electron R AA and v 2 .

Summary
We have presented a comprehensive set of results on the azimuthal anisotropy of charm production at central rapidity in Pb-Pb collisions at √ s NN = 2.76 TeV, obtained by reconstructing the decays D 0 → K − π + , D + → K − π + π + and D * + → D 0 π + .
Azimuthal anisotropy of charm production in Pb-Pb collisions at √ s NN = 2.76 TeV 27 The azimuthal anisotropy parameter v 2 was measured with the event plane, scalar product and twoparticle cumulant methods, as a function of transverse momentum for semi-central collisions in the 30-50% quantile of the hadronic cross section.The measured anisotropy was found to be consistent among D meson species, as well as for the three methods.The average v 2 of the three mesons in the interval 2 < p T < 6 GeV/c is larger than zero with a significance of 5.7 σ , combining statistical and systematic uncertainties.With a smaller significance, a positive v 2 is also observed for p T > 6 GeV/c, likely to originate from a path-length dependence of the partonic energy loss.The azimuthal anisotropy of D 0 mesons, which have larger statistical significance than D + and D * + , was also measured in the centrality classes 0-10% and 10-30%.For all three centrality classes, the D 0 meson v 2 is comparable in magnitude to that of inclusive charged particles.An indication for a decrease of v 2 towards more central collisions was observed for 3 < p T < 6 GeV/c.
The anisotropy was also quantified in terms of the D 0 meson nuclear modification factor R AA , measured in the direction of the event plane and orthogonal to it.For p T > 3 GeV/c, a stronger suppression relative to proton-proton collisions is observed in the out-of-plane direction, where the average path length of heavy quarks through the medium is larger.
The results indicate that, during the collective expansion of the medium, the interactions between its constituents and charm quarks transfer to the latter information on the azimuthal anisotropy of the system.
The new results for v 2 and R AA measured in and out of the event plane, as well as previously published R AA in the most central collisions [13], were compared with model calculations.The anisotropy is best described by the models that include mechanisms, like collisional energy loss, that transfer to charm quarks the elliptic flow induced during the system expansion.In some of these models the charmed meson v 2 is further enhanced by charm quark recombination with light quarks from the medium.However, it is challenging for models to describe simultaneously the large suppression of D mesons in central collisions and their anisotropy in non-central collisions.The results reported in this article provide important constraints on the mechanisms of heavy-quark energy loss and on the transport properties of the expanding medium produced in high-energy heavy-ion collisions.

Figure 1 :
Figure 1: (Color online) (a) Distribution of event plane angle ψ 2 , estimated from TPC tracks with 0 < η < 0.8 (solid line) or with the VZERO detector signals (dashed line) in the centrality range 30-50%.The distributions are normalized by their integral.(b) Event plane resolution correction factor R 2 as a function of centrality for the TPC and VZERO detectors.The boxes represent the systematic uncertainties estimated from the variation of R 2 when changing the sub-events used for its determination.

Figure 2 :
Figure 2: (Color online) Distributions of the invariant mass for D 0 (upper panels) and D + (central panels) candidates and of the mass difference for D * + candidates (lower panels) in the two ∆ϕ intervals used in the event plane method, for Pb-Pb collisions in the 30-50% centrality class.The rapidity interval is |y| < y fid (see Section 3.2 for details).For each meson species three p T intervals are shown, along with the fits used to extract the signal yield.The definition of the two ∆ϕ intervals is sketched in the top-left panel.

Figure 3 :
Figure 3: (Color online) Examples of v 2 extraction with two-particle correlation methods in a selected p T interval for Pb-Pb collisions in the 30-50% centrality range: the two-particle cumulants method for D 0 (a) and the scalar product method for D * + (b).The lower panels report the D meson v 2 values obtained with the simultaneous fit procedure, as described in the text.The rapidity interval is |y| < y fid (see Section 3.2 for details).

Figure 4 :
Figure 4: (Color online) Product of acceptance and efficiency for D 0 mesons in Pb-Pb collisions for 30-50% centrality class (upper panel).The rapidity interval is |y| < y fid (see Section 3.2 for details).The values for prompt (solid lines) and feed-down (dotted lines) D 0 mesons are shown.Also displayed, for comparison, are the values for prompt D 0 mesons without PID selection (dashed lines).The lower panel shows the ratio of the efficiencies for prompt D 0 mesons in the in-plane and out-of-plane regions used for the analysis.This ratio was estimated using simulation samples with a difference in particle multiplicity similar to that observed in data for the two azimuthal regions.

Figure 5 :
Figure 5: (Color online) Invariant mass distribution of D 0 candidates with 4 < p T < 6 GeV/c in the centrality class 30-50%: (a) fit without template for reflections and (b) with template for reflections (dotted line).The raw yield obtained as integral of the signal Gaussian function is reported.

Figure 6 :
Figure 6: (Color online) v 2 as a function of p T in the 30-50% centrality class, for D 0 , D + and D * + mesons (rows)with the event plane (from Ref.[47]), scalar product and two-particle cumulant methods (columns).For the first method, the event plane was estimated with TPC tracks in 0 < η < 0.8; for the other methods, TPC tracks in −0.8 < η < 0.8 were used as RFP.The symbols are positioned at the average p T measured within each interval.

Figure 7 :
Figure 7: (Color online) D 0 meson v 2 as a function of p T in the 30-50% centrality class, with the reference particles from the TPC or from the VZERO detectors (−3.7 < η < −1.7 and 2.8 < η < 5.1).(a) Event plane method.(b) Scalar product method.For visibility, the symbols for the VZERO case are slightly displaced horizontally.

Figure 8 :Figure 9 :
Figure 8: (Color online) Comparison of prompt D 0 meson and charged-particle v 2[36] in three centrality classes as a function of p T .Both measurements are done with the event plane method.For charged particles a gap of two η units is used.

Figure 10 :
Figure 10: (Color online) Nuclear modification factor R AA of D 0 mesons in the 30-50% centrality class in two 90 •wide azimuthal intervals centred on the in-plane and on the out-of-plane directions.The correlated, uncorrelated, and anti-correlated contributions to the systematic uncertainty are shown separately.

Table 2 :
Systematic uncertainties on the measurement of v 2 in the 30-50% centrality class for the interval 4 < p T < 6 GeV/c.The uncertainties are comparable in the other p T intervals.

Table 3 :
Systematic uncertainties on the measurement of the D 0 meson R AA in-plane and out-of-plane in the 30-50% centrality class for two p T intervals.The uncertainties are grouped according to the type of correlation between the in-plane and out-of-plane cases. )