Elliptic and triangular flow of (anti)deuterons in Pb-Pb collisions at $\sqrt {s_{\rm NN}} = 5.02$ TeV

The measurements of the (anti)deuteron elliptic flow (v2) and the first measurements of triangular flow (v3) in Pb-Pb collisions at a center-of-mass energy per nucleon-nucleon collision √ sNN = 5.02 TeV are presented. A mass ordering at low transverse momentum (pT) is observed when comparing these measurements with those of other identified hadrons, as expected from relativistic hydrodynamics. The measured (anti)deuteron v2 lies between the predictions from the simple coalescence and blast-wave models, which provide a good description of the data only for more peripheral and for more central collisions, respectively. The mass number scaling, which is violated for v2, is approximately valid for the (anti)deuterons v3. The measured v2 and v3 are also compared with the predictions from a coalescence approach with phase-space distributions of nucleons generated by IEBEVISHNU with AMPT initial conditions coupled with URQMD, and from a dynamical model based on relativistic hydrodynamics coupled to the hadronic afterburner SMASH. The model predictions are consistent with the data within the uncertainties in midcentral collisions, while a deviation is observed in the most central collisions.


I. INTRODUCTION
The production mechanism of light (anti)nuclei in highenergy hadronic collisions is still not fully clear and is under intense debate in the scientific community [1][2][3][4][5]. The understanding of the production of loosely bound multibaryon states in heavy-ion collisions has additional complications due to the fact that the phase transition is followed by a hadron gas phase with intense rescattering of hadrons. At the Large Hadron Collider (LHC) energies, the lifetime of the hadronic phase between chemical and kinetic freeze-out is in the range 4-7 fm/c [6] and the kinetic freeze-out temperature, when elastic interactions cease, is of the order of 100 MeV [7,8]. The binding energy of multibaryon systems such as light (anti)nuclei typically does not exceed a few MeV, which is almost two orders of magnitude smaller than the temperature of the system. Considering the high density of hadrons in the posthadronization stage and the large dissociation cross sections of light (anti)nuclei, it is not clear how such loosely bound systems can survive under these extreme conditions.
Existing phenomenological models provide very different interpretations for this observation. In the statistical hadronization model [1][2][3]9,10], light (anti)nuclei as well as all other hadron species are assumed to be emitted by a source in local thermal and hadrochemical equilibrium. Their abundances are fixed at the chemical freeze-out, occurring at a temperature of T chem = 156 ± 4 MeV for Pb-Pb collisions at the LHC [11]. This model provides a good description of the measured hadron yields in central nucleus-nucleus collisions [1]. However, the mechanism of hadron production and the propagation of loosely bound states through the hadron gas phase are not addressed by this model. In the context of the statistical hadronization model, it has been conjectured that such objects could be produced at the phase transition as compact colorless quark clusters with the same quantum numbers of the final-state hadrons. The survival of these states at high temperatures is interpreted as due to the low interaction cross section with the surrounding medium [1].
In the coalescence approach, multibaryon states are assumed to be formed by the coalescence of baryons at the kinetic freeze-out. In the simplest versions of this model [12,13], baryons are treated as pointlike particles and the coalescence happens instantaneously if the momentum difference between nucleons is smaller than a given threshold, which is typically of the order of 100 MeV/c, while spatial coordinates are ignored. In the state-of-the-art implementations of the coalescence approach [4,14], the quantum-mechanical properties of baryons and their bound states are taken into account and the coalescence probability is calculated from the overlap between the wave functions of baryons and the Wigner density of the final-state cluster. All light (anti)nuclei produced at the phase transition are assumed to be destroyed by the interactions in the hadron gas phase and regenerated with the same amount only at the latest stage of the system evolution.
To address the open question of the survival of loosely bound multibaryon states in the hadron gas phase with intense rescattering, models based on relativistic hydrodynamics coupled to a hadronic afterburner have been recently developed [4,5]. In these models, nucleons and light nuclei are produced at the phase transition using the Cooper-Frye formula [15], which describes the hadron production based on the local energy density of the fireball, and their yields are fixed to the value predicted by the thermal model at the chemical freeze-out temperature. Their propagation through the hadronic medium is simulated based on known interaction cross sections and resonant states using different transport codes. Existing calculations are based on URQMD [16,17], with light nuclei being produced by nucleon coalescence, and SMASH [5], where (anti)deuterons are assumed to be destroyed and regenerated with equal rates in the hadronic stage. The model based on URQMD with nucleon coalescence [4] provides a good description of the elliptic flow of (anti)deuterons measured in Pb-Pb collisions at √ s NN = 2.76 TeV [18] and of that of (anti) 3 He measured in Pb-Pb collisions at √ s NN = 5.02 TeV [19]. The model is able to describe the lowp T spectra of deuterons, but overpredicts the deuteron data above 2.5 GeV/c and the (anti) 3 He spectra in the full momentum interval. The hybrid model based on SMASH successfully describes the measured (anti)deuteron p T spectra and coalescence parameter B 2 , defined as the ratio of the invariant yield of deuterons and that of protons squared, measured in Pb-Pb collisions at √ s NN = 2.76 TeV [18].
A conceptually similar approach, based on the analogy between the evolution of the early universe after the Big Bang and the space-time evolution of the system created in heavy-ion collisions, has recently been developed [20]. The production of light (anti)(hyper)nuclei in heavy-ion collisions at the LHC is considered in the framework of the Saha equation assuming that disintegration and regeneration reactions involving light nuclei proceed in relative chemical equilibrium after the chemical freeze-out of hadrons.
The existing models depict radically different pictures of the posthadronization stage for loosely bound states. Considering this scenario, the measurements of radial and anisotropic flow of light (anti)nuclei, i. e., the harmonics (v n ) of the Fourier decomposition of their azimuthal production distribution with respect to a symmetry plane of the collision, are relevant to study their propagation through the hadron gas phase and the dynamics of their interactions with other particles. Compared to the elliptic flow, the triangular flow of light (anti)nuclei has a better sensitivity to the fluctuating initial conditions as well as the properties of the created systems. Therefore, tighter constraints on the theoretical model that describe the production mechanism of light (anti)nuclei can be set.
The elliptic flow of (anti)deuterons was measured as a function of the transverse momentum (p T ) for different centrality classes in Pb-Pb collisions at √ s NN = 2.76 TeV [18]. A clear mass ordering is observed at low p T (p T < 3 GeV/c) when this measurement is compared to that of other hadron species [21], as expected from relativistic hydrodynamics. The simple coalescence model, based on the assumption that the (anti)deuteron invariant yield is proportional to the invariant yield of (anti)protons squared, is found to overestimate the measured v 2 in all centrality intervals. The data are better described by the blast-wave model, a simplified version of the relativistic hydrodynamic approach in which the collective expansion is described using a parametrized hydrodynamic flow field. The elliptic flow of (anti) 3 He was measured in Pb-Pb collisions at √ s NN = 5.02 TeV [19]. Also in the case of (anti) 3 He, the mass ordering is observed for p T < 3 GeV/c and the measured elliptic flow lies between the predictions of the blast-wave [22] and the simple coalescence model. A better description of the measurement is provided by a more sophisticated coalescence model where the phase-space distributions of protons and neutrons are generated by the IEBE-VISHNU hybrid model with AMPT initial conditions [4].
The picture that has emerged so far, regarding the elliptic flow of (anti)nuclei measured at LHC energies, is that the simple coalescence and blast-wave models represent the upper and lower edges of a region where the data are mostly located. Recent developments in the coalescence approach, which take into account momentum-space correlations of nucleons and their quantum-mechanical properties, provide a better description of the data [4,5].
In this paper, a precision measurement of the (anti)deuteron elliptic flow and the first ever measurement of (anti)deuteron triangular flow for different p T and centrality intervals in Pb-Pb collisions at √ s NN = 5.02 TeV are presented. Thanks to the large data sample collected at higher energy, the elliptic flow measurement is performed in wider p T and up to a higher centrality interval compared to that in Pb-Pb collisions at √ s NN = 2.76 TeV, allowing for a more differential comparison with the theoretical models.

II. THE ALICE DETECTOR
A detailed description of the ALICE detector can be found in Ref. [23] and references therein. The main subdetectors used for the present analysis are the V0 detector, the inner tracking system (ITS), the time projection chamber (TPC), and the time-of-flight (TOF) detector, which are located inside a solenoidal magnet that provides a uniform field of 0.5 T directed along the beam direction. The V0 detector [24] consists of two arrays of scintillation counters placed around the beam vacuum tube on either side of the interaction point: one covering the pseudorapidity interval 2.8 < η < 5.1 (V0A) and the other one covering −3.7 < η < −1.7 (V0C). Each V0 array consists of four rings in the radial direction, with each ring composed of eight cells with the same azimuthal size. The scintillator arrays have an intrinsic time resolution better than 0.5 ns, and their timing information is used in coincidence for offline rejection of events produced by the interaction of the beams with residual gas in the vacuum pipe. The V0 scintillators are used to determine the collision centrality from the measured charged-particle multiplicity [25][26][27] and to measure the orientation of the symmetry plane of the collision.
The ITS [28], designed to provide high-resolution track points in the vicinity of the nominal vertex position, is composed of three subsystems of silicon detectors placed around the interaction region with a cylindrical symmetry. The silicon pixel detector (SPD) is the subsystem closest to the beam vacuum tube and it is made of two layers of pixel detectors. The third and the fourth layers are formed by silicon drift detectors, while the outermost two layers are equipped with double-sided silicon strip detectors. The ITS covers the pseudorapidity interval |η| < 0.9.
The same pseudorapidity interval is covered by the TPC, which is the main tracking detector, consisting of a hollow cylinder the axis of which coincides with the nominal beam axis.
The active volume of 90 m 3 is filled with a gas mixture containing 88% Ar and 12% CO 2 .
The trajectory of a charged particle is estimated using up to 159 space points. The charged-particle tracks are then built by combining the hits in the ITS and the reconstructed space points in the TPC. The TPC is also used for particle identification (PID) by measuring the specific energy loss (dE/dx) in the TPC gas.
The TOF detector [29] covers the full azimuth in the pseudorapidity interval |η| < 0.9. The detector is based on multigap resistive plate chamber technology and it is located, with cylindrical symmetry, at an average radial distance of 380 cm from the beam axis. The TOF allows for PID, based on the difference between the measured time of flight and its expected value, computed for each mass hypothesis from the track momentum and length. The resolution on the measurement of the time of flight is about 60 ps in heavy-ion collisions.

A. Event and track selections
The data sample used for the measurements presented in this paper was recorded by ALICE in 2015 during the LHC Pb-Pb run at √ s NN = 5.02 TeV. A minimum bias trigger was used during the data taking, which required coincident signals in both V0 detectors. An offline event selection is applied to remove beam-gas collisions using the timing information provided by the V0 detectors and the zero-degree calorimeters [23]. Events with multiple primary vertices identified with the SPD are tagged as pileup and removed from the analysis. In addition, events with significantly different charged-particle multiplicities measured by the V0 detector and by the tracking detectors at midrapidity, which have different readout times, are rejected. After the offline event selection, the remaining contribution of beam-gas events is smaller than 0.02% [23] and the fraction of pileup events is found to be negligible. The primary vertex position is determined from tracks reconstructed in the ITS and TPC as described in Ref. [23] and only events with a reconstructed primary vertex position along the beam axis within 10 cm from the nominal interaction point are selected. The total number of events selected for the analysis for centrality 0-70% is about 73 × 10 6 .
Deuteron (d) and antideuteron (d) candidates are selected from charged-particle tracks reconstructed in the ITS and TPC in the kinematic range |η| < 0.8 and 0.8 < p T < 6 GeV/c. Only tracks with at least 70 clusters out of a maximum of 159 and with a χ 2 per degree of freedom for the track fit lower than 2 are accepted. In addition, in order to guarantee a trackmomentum resolution of 2% in the measured p T range and a dE/dx resolution of about 6%, each track is required to be reconstructed from at least 80% of the number of expected TPC clusters and to have at least one hit in either of the two innermost layers of the ITS. The distances of closest approach (DCA) to the primary vertex in the plane perpendicular and parallel to the beam axis for the selected tracks are determined with a resolution better than 300 μm [23]. To suppress the contribution of secondary particles, the reconstructed tracks are required to have a longitudinal DCA smaller than 2 cm and a transverse DCA smaller than 0.0105 + 0.0350/p 1.1 T cm, with p T in units of GeV/c. The latter corresponds to approximately 7σ DCA (p T ), where σ DCA (p T ) is the transverse DCA resolution in the corresponding p T interval.

B. (Anti)deuterons identification
The (anti)deuteron identification technique used in this analysis is similar to that used in the previous measurement in Pb-Pb collisions at √ s NN = 2.76 TeV [18]. For transverse momenta up to 1.4 GeV/c (anti)deuterons are identified using only the TPC information by requiring that the average dE/dx is within 3σ from the expected average value for the (anti)deuteron mass hypothesis. For p T > 1.4 GeV/c the 3σ TPC identification is complemented by the signal provided by the TOF detector. The number of (anti)deuterons in each where m TOF is the particle mass calculated using the time of flight measured by the TOF and m d pdg is the nominal mass of deuterons taken from [30]. In the left panel of Fig. 1 the M distribution for (anti)deuterons with 2.2 p T <2.4 GeV/c in the centrality interval 20-30% is shown. The d + d signal is fitted using a Gaussian with an exponential tail, while the background, originating from TOF hits incorrectly associated to tracks extrapolated from the TPC, is modeled with an exponential function.
Deuterons and antideuterons are summed together (d + d) in all the centrality intervals and for p T larger than 1.4 GeV/c. This is possible since the v 2 and v 3 measured for v 2 and v 3 for d and d are consistent within the statistical uncertainties. At lower p T , deuterons produced by spallation in interactions between particles and the detector material or in the beam vacuum tube constitute a significant background. For this reason, for p T < 1.4 GeV/c only antideuterons, which are not affected by this background, are used in the analysis. Since no difference is expected for the v 2 and v 3 of d and d, hereafter deuterons will denote results for antideuterons for p T < 1.4 GeV/c and the sum of d and d elsewhere. The contribution of secondary deuterons produced in weak decays of hypertritons is negligible considering that the production rate of (hyper)nuclei with mass number A = 3 is suppressed compared to that of A = 2 by a factor of approximately 300 in Pb-Pb collisions at √ s NN = 2.76 TeV [31]. A similar suppression is expected in Pb-Pb collisions at √ s NN = 5.02 TeV.

C. Flow analysis techniques
The particle azimuthal distribution of charged particles with respect to the nth-order flow symmetry plane n [32][33][34][35] can be expressed as a Fourier series: where E is the energy of the particle, p is the momentum, ϕ is the azimuthal angle, y is the rapidity, and (2) The second coefficient of the Fourier series (v 2 ) is called elliptic flow and is related to the initial geometrical anisotropy of the overlap region of the colliding nuclei. The third-order flow coefficient (v 3 ), called triangular flow, is generated by fluctuations in the initial distribution of nucleons and gluons in the overlap region [34,36,37]. The same fluctuations are responsible for the v 2 measured in most central collisions (centrality < 5%) [38]. The v n coefficients are measured using the scalar product (SP) method [32,39]. This is a two-particle correlation technique based on the scalar product of the unit flow vector of the particle of interest, k, and the Q vector. The unit flow vector is denoted by u n,k = exp(inϕ k ), where ϕ k is the azimuthal angle of the particle k.
The Q-vector is computed from a set of reference flow particles and is defined as where, in general, ϕ i is the azimuthal angle for the ith reference flow particle, n is the order of the harmonic, and w i is a weight applied to correct for reference flow.
The v n flow coefficients are calculated as Single brackets ... denote an average over all events, while double brackets ... indicate an average over all particles in all events, and the asterisk denotes the complex conjugate. The denominator is a correction factor that is introduced to take into account the resolution of the Q n vector. In this analysis, the Q n vector is calculated from the azimuthal distribution of the energy deposition measured in the V0A, while the Q A n and Q B n vectors are determined from the azimuthal distribution of the energy deposited in the V0C and the azimuthal distribution of tracks reconstructed in the TPC, respectively. Using these detectors, a pseudorapidity gap | η| > 2 between the particle of interest and the reference flow particles is introduced. Such a pseudorapidity gap reduces nonflow effects, which are correlations not arising from the collective expansion of the system (e. g., resonance decays and jets).
The purity of the sample of deuterons identified using the TPC in the 0.8 < p T < 1.4 GeV/c interval is around 100%. In this transverse momentum interval the v 2 and v 3 coefficients were evaluated on a track-by-track basis and then averaged in each p T interval. For higher p T , the v n coefficients are calculated in different ranges of M. v n ( M) contains contributions from the signal (v sig n ) and from the background (v bkg (5) where N sig is the number of deuterons, N bkg is the number of background particles, and N tot is their sum. The signal v n is extracted from a fit to the observed v n as a function of M, in which v bkg n is described using a first-order polynomial function, and v sig n is a free fit parameter. N sig and N bkg are obtained from the fit to the M distribution using a Gaussian with an exponential tail for the signal and an exponential for the background. The signal extraction procedure is illustrated in Fig. 1 The elliptic and triangular flows of deuterons are measured in centrality intervals of 5% width and then the results in wider centrality intervals are obtained as weighted averages of these measurements using the number of deuteron candidates, in the same centrality interval of 5% width as a weight, similarly to what was performed in Ref. [19].

D. Systematic uncertainties
The sources of systematic uncertainties for the elliptic and triangular flow of deuterons are related to event selection, tracking, (anti-)deuteron identification, and the technique used for the signal extraction. The contribution related to the event selection is estimated by taking into account the differences 055203-4 in the v 2 and v 3 measurements obtained using different eventselection criteria. In particular, the fiducial region for the vertex position along the beam axis is varied from the range [−10, 10] to [−7, 7] cm to probe the magnitude of potential edge effects. To investigate possible effects due to charge asymmetries during tracking and geometrical asymmetries in the detector, the differences between the results obtained by using opposite magnetic-field polarities are included. Analogously, the default centrality estimator is changed to that based on the number of hits in the first or second layer of the ITS. Finally, the effect related to pileup rejection is tested by requiring a stronger correlation between the V0 and central barrel multiplicities. These contributions are assumed to be independent and added in quadrature. The total systematic uncertainty due to event selection is found to be around 1.5% for both v 2 and v 3 .
To estimate the systematic uncertainties due to reconstruction and identification of deuterons, the track selection and the TPC PID criteria are varied with respect to the default choice and the v n measurements are repeated for each of these different settings. The rms of the distribution of v n measurements in each p T interval is considered as systematic uncertainty. To minimize the effect of statistical fluctuations, all variations smaller than 2 |σ 2 0 − σ 2 i | are not included in the estimate of the systematic uncertainties [40], where σ 0 is the statistical uncertainty of the default value while σ i is that corresponding to the ith selection criterion. The probability distribution for the variations of data points due to systematic effects related to tracking and PID is assumed to be uniform in each p T interval and the difference between the maximum and minimum value divided by √ 12 is assigned as systematic uncertainty. This contribution ranges from 1 and 3% depending on p T and centrality.
To estimate the contribution to the systematic uncertainties due to the signal extraction, the function used to describe the v bkg n is changed. In addition to a first-order polynomial, a constant function and a second-order polynomial are also used, and the maximum difference with respect to the default measurement is considered as systematic uncertainty. A contribution up to 5% is observed for central collisions and for p T < 2 GeV/c. Moreover, different functions and fitting ranges are used to describe the signal and the background of Eq. (5). More specifically, besides a Gaussian function with an exponential tail, a Gaussian is also used for the signal, while single and double exponentials and linear functions are also used for the background. This contribution is relevant only for p T > 1.4 GeV/c, where the TOF is used to extract the signal, and is found to vary from 1 to 6% depending on p T and centrality. Table I shows the summary of the different contributions to the systematic uncertainties for the v 2 and v 3 of deuterons. The total uncertainties are given by their sum in quadrature, assuming that all contributions are independent.

IV. RESULTS AND DISCUSSION
The v 2 and v 3 of deuterons measured in Pb-Pb collisions at √ s NN = 5.02 TeV are shown in Fig. 2 as a function of p T for different centrality intervals. In the measured p T interval, an increasing trend is observed with increasing p T and going from central to more peripheral Pb-Pb collisions, as expected based on the relativistic hydrodynamic description of the collective expansion of a hot and dense medium [41]. Initial-state fluctuations of the energy density distribution of partons in the colliding nuclei imply a nonzero v 3 [37].
The measurement presented in this paper shows that these initial-state effects, already observed for other hadron species at LHC energies [42,43], are also visible for deuterons.
The measurement of the deuterons v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV is compared to that in Pb-Pb collisions at √ s NN = 2.76 TeV [18] in Fig. 3 for two centrality intervals. The observed v 2 and their trend are similar at the two center-of-mass energies, but a decrease of the observed elliptic flow for a given p T is observed with increasing centerof-mass energy. This effect is more pronounced in peripheral rather than in central collisions. A similar effect was observed for the proton v 2 measurements [43] and is interpreted as partially due to the increasing radial flow with increasing collision energy, which produces a shift of the v 2 towards higher p T . The effect due to radial flow is assessed quantitatively by comparing the ratio of the deuteron and proton v 2 as a function of p T at the two energies. The ratio between the deuteron v 2 in Pb-Pb collisions at √ s NN = 5.02 TeV to that measured at √ s NN = 2.76 TeV, with v 2 and p T scaled by the mass number A = 2, is shown in Fig. 4 for two centrality intervals in comparison with the same ratio for protons. As indicated by these ratios, the radial flow effects are quantitatively very similar for protons and deuterons. It has to be noted that a mass scaling would lead to the same conclusion since the binding energy of deuterons is 2.2 MeV, i. e., the deuteron   3 He v 2 is measured using the event plane method [19]. Vertical bars and boxes represent the statistical and systematic uncertainties, respectively. mass is approximately equal to 2m p , where m p is the proton mass.
The elliptic flow of deuterons is compared to that of pions, kaons, protons, and (anti) 3 He measured at the same centerof-mass energy [19,43] in Fig. 5. Since the (anti) 3 He elliptic flow is measured in centrality intervals of 20% width due to its rarer production compared to that of lighter hadrons, the v 2 of pions, kaons, protons, and deuterons are recalculated to match the same centrality intervals. This is achieved by averaging the v 2 measurements of these particles in narrower centrality intervals weighted by the corresponding p T spectra [8,44]. A clear mass ordering of v 2 is observed at low p T , as expected for a system expansion driven by the pressure gradient as described by relativistic hydrodynamics [41,45,46].
In Fig. 6, the deuteron v 3 is compared to that of pions, kaons, and protons at the same center-of-mass energy [43] for the centrality intervals 0-20% (left) and 20-40% (right). Also for v 3 , a clear mass ordering is observed for p T 2.5 and 3 GeV/c for the centrality intervals 0-20 and 20-40%, respectively.

A. Comparison with the blast-wave model predictions
The elliptic flow of deuterons is compared with the expectations of the blast-wave model [22,47,48], which is based on the assumption that the system produced in heavy-ion collisions is locally thermalized and expands collectively with a common velocity field. The system is assumed to undergo an instantaneous kinetic freeze-out at the temperature T kin and to be characterized by a common transverse radial flow velocity at the freeze-out surface. A simultaneous fit of the v 2 and the p T spectra of pions, kaons, and protons [8,43] with the blast-wave model is performed in the transverse-momentum ranges 0.5 p π T < 1 GeV/c, 0.7 p K T < 2 GeV/c, and 0.7 p p T < 2.5 GeV/c. The four free parameters of the blast-wave function are the kinetic freeze-out temperature (T kin ), the variation in the azimuthal density of the source (s 2 ), the mean transverse expansion rapidity (ρ 0 ), and the amplitude of its azimuthal variation (ρ a ), as described in Ref. [47]. The values of these parameters extracted from the fits are reported in Table II for   assumption that the same kinetic freeze-out conditions apply for all particles produced in the collision. The deuteron mass is taken from [30]. The blast-wave fits to the v 2 of pions, kaons, and protons and the predictions for the deuterons v 2 are reported in Fig. 7 for the centrality intervals 10-20 and 40-50%. In the lower panels, the data-to-fit ratios for pions, kaons, and protons and the ratios of the deuterons v 2 to the model are shown. Because of the finite size of the p T intervals, the average of the blast-wave function within the interval, weighted by the p T spectrum of the corresponding particle species, is considered in the calculation of these ratios.
The predictions of the blast-wave model underestimate the deuteron elliptic flow experimental values in semiperipheral collisions for p T > 1.4 GeV/c, while they are close to the measurements for central events in the measured p T interval. This is better observed in Fig. 8, which shows the centrality evolution of the data-to-model ratios.

B. Test of the coalescence hypothesis
The deuterons v 2 and v 3 are compared to the expectations of a coalescence approach based on mass number scaling and isospin symmetry, for which the proton and neutron v 2 (v 3 ) are identical. In particular, the v 2 (v 3 ) measured for protons [43] was used to predict the v 2 (v 3 ) of deuterons using the following relation [49]: The results of this calculation for different centrality intervals for v 2 are shown in the left panel of Fig. 9. The measured elliptic flow in 10-20 and 40-50% centrality intervals of deuterons is compared with coalescence model predictions from Eq. (6) using the measured v n of protons. Similarly, the right panel of Fig. 9 shows a comparison between the FIG. 7. Blast-wave fits to the v 2 (p T ) of pions, kaons, and protons [43] and predictions of the deuterons v 2 (p T ) for the centrality intervals 10-20% (left) and 40-50% (right). In the lower panels, the data-to-fit ratios are shown for pions, kaons, and protons as well as the ratio of the deuterons v 2 to the blast-wave predictions. Vertical bars and boxes represent the statistical and systematic uncertainties, respectively. The dashed line is to guide the eye. The coalescence model overestimates the deuteron v 2 by about 20 to 30% in central collisions and is close to the data for semiperipheral collisions, as illustrated in Fig. 10, which shows the centrality evolution of the data-to-model ratio. The coalescence approach seems to have a slightly better agreement with deuterons v 3 ; however, the large statistical uncertainties on the v 3 measurements do not allow for conclusive statements.

C. Comparison with IEBE-VISHNU and coalescence calculations
In Fig. 11, the deuterons v 2 and v 3 are compared to a model [4] implementing light nuclei formation via coalescence of nucleons originating from a hydrodynamical evolution of the fireball coupled to a URQMD simulation of the hadronic cascade [16,17]. In this model, the coalescence probability is calculated as the superposition of the wave functions of protons and neutrons and the Wigner function of the deuterons. The coalescence happens in a flowing  [5] based on the JETSCAPE generator [52]. The model predictions, based on the SMASH afterburner and which used TRENTo [51] initial conditions, are shown as bands. The width of the band represents the statistical uncertainty associated with the model. In the lower panel the data-to-model differences are shown. Vertical bars and boxes represent the statistical and systematic uncertainties, respectively.

D. Comparison with hybrid (hydrodynamics plus transport) approach expectations
The deuterons v 2 measured in the centrality intervals 10-20, 20-30, and 30-40% are compared in Fig. 12 with the predictions from a hybrid model based on relativistic viscous hydrodynamics, with fluctuating initial conditions generated by TRENTO [51], coupled to the hadronic afterburner SMASH [5]. The simulations are obtained by using the JETSCAPE 1.0 event generator [52]. The parameters of this model, including the shear and bulk viscosities, are tuned to the measurements of p T spectra and azimuthal flow of pions, kaons, and protons obtained by ALICE in Pb-Pb collisions at √ s NN = 2.76 TeV [7,21] and by PHENIX and STAR in Au-Au collisions at √ s NN = 200 GeV [53][54][55]. The interactions of deuterons with other hadrons in the hadron gas phase are simulated using SMASH in which all known resonances and the experimentally known cross sections, most importantly π d → π np and its inverse reaction, are included.
In this model, the number of deuterons at the kinetic freeze-out is independent from their primordial abundance at the Cooper-Frye hypersurface. It was found that even when their initial number is set to zero the number of deuterons regenerated in the hadronic phase converges towards the equilibrium value, which is the same as that predicted by the statistical hadronization model. Considering that in this model only ≈1% of the primordial deuterons survive the hadronic stage, the elliptic flow of deuterons observed after the kinetic freeze-out is almost identical to that of the regenerated ones. For this reason, deuterons are not sampled at the Cooper-Frye hypersurface for these predictions.
The model predictions are consistent with the measured v 2 within the uncertainties in the centrality intervals 20-30 and 30-40% for 0.8< p T < 4 GeV/c, while the data are overestimated by up to 30% in the centrality interval 10-20% for p T > 2 GeV/c.

V. SUMMARY
The measurements of the deuterons v 2 and the first measurement of v 3 in Pb-Pb collisions at √ s NN = 5.02 TeV are presented. The observed centrality and p T dependence are consistent with the expectations from relativistic hydrodynamics. A mass ordering is observed for p T < 3 GeV/c when comparing these results with the measured v 2 and v 3 of pions, kaons, and protons. The shift of the deuterons v 2 towards higher p T with respect to the measurement in Pb-Pb collisions at √ s NN = 2.76 TeV, mainly due to a stronger radial flow at higher center-of-mass energy, is consistent with that observed for the proton v 2 measurement. The results of this measurement are compared with the expectations from the simple coalescence approach, in which the deuteron v 2 is obtained from that of protons assuming that the deuteron invariant yield is proportional to that of protons squared, and with the predictions of the blast-wave model. The deuteron v 2 is overestimated by a simple coalescence approach, which describes the data only in peripheral (centrality >50%) collisions. On the other hand, the blast-wave model underestimates the peripheral measurements and it is close to the data in central collisions. These results are consistent with the scenario previously seen for deuterons and 3 He elliptic flow: these simplified models bracket a region where the light nucleus v 2 is located and describe reasonably the data in different multiplicity regimes, indicating that neither of these two models is able to describe the deuteron production measurement from low to high multiplicity environments.
Similar considerations are valid for the deuterons v 3 with some limitations due to the rather large statistical uncertainties. This specific aspect will be addressed with the larger data sample that will be collected in run 3 following the ALICE upgrade, where a significant improvement of the statistical precision is expected. This measurement will be crucial to better constrain models that describe the production of light nuclei in heavy-ion collisions.
A more advanced coalescence model coupled to hydrodynamics and the hadronic afterburner URQMD, which takes into account the quantum-mechanical properties of nucleons and nuclei and space-momentum correlations of nucleons, provides a good description of the deuterons v 2 and v 3 for p T > 2.5 GeV/c. The model predictions deviate from the data at lower p T , in particular for the centrality interval 10-20%. The same model provides a good description of the deuteron v 2 measured in Pb-Pb collisions at √ s NN = 2.76 TeV and that of 3 He at √ s NN = 5.02 TeV. The deuteron v 2 is also compared to the predictions from a hybrid model based on relativistic hydrodynamics coupled to the hadronic afterburner SMASH.

055203-11
The model predictions are consistent with the data within the uncertainties in the centrality intervals 20-30 and 30-40%, while a deviation of up to 30% is observed in the centrality interval 10-20% for 2 < p T < 3 GeV/c.
In general, the state-of-the-art implementations of coalescence and the hybrid approach based on hydrodynamics coupled to hadronic afterburners provide better descriptions of the data compared to the simple coalescence and blast-wave models. Further efforts, on both the experimental and the theoretical side, are needed to have a more comprehensive understanding of dynamics and production of light nuclei.

ACKNOWLEDGMENTS
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centers and the Worldwide LHC Computing Grid collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the AL-