Anisotropic flow of identified particles in Pb-Pb collisions at $\mathbf{\sqrt{{\textit s}_{\rm NN}}}=5.02$ TeV

The elliptic ($v_2$), triangular ($v_3$), and quadrangular ($v_4$) flow coefficients of $\pi^{\pm}$, ${\rm K}^{\pm}$, ${\rm p+\overline{p}}$, ${\Lambda+\overline{\Lambda}}$, ${\rm K}^{\rm 0}_{\rm S}$, and the $\phi$-meson are measured in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. Results obtained with the scalar product method are reported for the rapidity range $\vert y \vert<$ 0.5 as a function of transverse momentum, $p_\text{T}$, at different collision centrality intervals between 0-70%, including ultra-central (0-1%) collisions for $\pi^{\pm}$, ${\rm K}^{\pm}$, and ${\rm p+\overline{p}}$. For $p_\text{T}<3$ GeV$\kern-0.05em/\kern-0.02em c$, the flow coefficients exhibit a particle mass dependence. At intermediate transverse momenta ($3<p_\text{T}<$~8-10 GeV$\kern-0.05em/\kern-0.02em c$), particles show an approximate grouping according to their type (i.e., mesons and baryons). The $\phi$-meson $v_2$, which tests both particle mass dependence and type scaling, follows ${\rm p+\overline{p}}$ $v_2$ at low $p_\text{T}$ and $\pi^{\pm}$ $v_2$ at intermediate $p_\text{T}$. The evolution of the shape of $v_{\rm n}(p_{\mathrm{T}})$ as a function of centrality and harmonic number $n$ is studied for the various particle species. Flow coefficients of $\pi^{\pm}$, ${\rm K}^{\pm}$, and ${\rm p+\overline{p}}$ for $p_\text{T}<3$ GeV$\kern-0.05em/\kern-0.02em c$ are compared to iEBE-VISHNU and MUSIC hydrodynamical calculations coupled to a hadronic cascade model (UrQMD). The iEBE-VISHNU calculations describe the results fairly well for $p_\text{T}<2.5$ GeV$\kern-0.05em/\kern-0.02em c$, while MUSIC calculations reproduce the measurements for $p_\text{T}<1$ GeV$\kern-0.05em/\kern-0.02em c$. A comparison to $v_{\rm n}$ coefficients measured in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV is also provided.


Introduction
Ultra-relativistic heavy-ion collisions are used to study the properties of the quark-gluon plasma (QGP), a state of deconfined quarks and gluons expected at high temperatures or baryon densities [1]. Measurements of anisotropies in particle azimuthal distributions relative to the collision symmetry planes at the Relativistic Heavy Ion Collider (RHIC) [2][3][4][5] and the Large Hadron Collider (LHC) [6][7][8] have shown that the produced hot and dense matter behaves as a strongly-interacting QGP. Comparisons to predictions from hydrodynamic models indicate that the QGP has a shear viscosity to entropy density ratio (η/s) close to the theoretical lower limit from the anti-de Sitter/conformal field theory (AdS/CFT) correspondence of 1/4π forh = k B = 1 [9].
Azimuthal anisotropies in particle production relative to the collision symmetry planes, often referred to as anisotropic flow, arise from the asymmetry in the initial geometry of the collision combined with the initial inhomogeneities of the system's energy density [10]. Anisotropic flow depends on the equation of state and transport coefficients of the system, such as η/s and bulk viscosity to entropy density ratio (ζ /s). Its magnitude is quantified via the coefficients v n in a Fourier decomposition of the particle azimuthal distribution [11] E d 3 N dp 3 = 1 2π d 2 N p T dp T dy (1 + 2 where E is the energy, p the momentum, p T the transverse momentum, ϕ the azimuthal angle, η the pseudorapidity of the particle, and Ψ n the n-th harmonic symmetry plane angle. The second order flow coefficient v 2 , called elliptic flow, is the largest contribution to the asymmetry of non-central collisions because of the almond-like geometry of the overlap region between the colliding nuclei in the plane perpendicular to the beam direction. The third-order flow coefficient v 3 , named triangular flow, is generated by fluctuations in the initial distribution of nucleons and gluons in the overlap region [12][13][14][15]. The fourth-order flow coefficient v 4 , called quandrangular flow, is generated both by initial geometry, fluctuations, and is in addition sensitive to the non-linear hydrodynamic response of the medium [16,17]. It has been shown that higher-order flow coefficients are more sensitive to η/s than v 2 [18,19]. In addition to probing η/s and ζ /s, anisotropic flow constrains the initial spatial density (e.g. energy and entropy density), freeze-out conditions of the system, and particle production mechanisms in different p T regions. Stronger constraints are achieved by studying anisotropic flow of identified particles. To guide interpretation of the results in the context of these processes, three kinematic 'regions of interest' are defined in the p T -differential v n measurements, v n (p T ). For p T 3 GeV/c, anisotropic flow is a remnant of the collective dynamics during the hydrodynamic expansion of the system. The interplay between the isotropic expansion (radial flow) and anisotropic flow leads to a characteristic mass ordering of v n (p T ) [20][21][22][23][24][25][26][27][28], meaning that heavier particles have smaller v n (p T ). At intermediate p T (3 p T 8 GeV/c), the values of v n for different particles tend to separate mesons and baryons [27][28][29][30][31][32][33]. The flow of baryons is larger than that of mesons in this p T range, supporting the hypothesis of hadronization through quark coalescence [34], where it is assumed that the invariant spectrum of produced particles is proportional to the product of the spectra of their constituents [35,36]. However, the scaling only holds approximately at RHIC [32] and at the level of ±20% in Pb-Pb collisions at √ s NN = 2.76 TeV [27,28]. This behaviour is also qualitatively consistent with a scenario in which particle production includes interactions of jet fragments with bulk matter [37,38]. For p T 8 GeV/c, anisotropic flow is generated when hard partons that propagate through the system lose energy via (multiple) scattering and gluon radiation [39,40], resulting in v n that remain non-zero up to very high p T [41][42][43][44].
Anisotropic flow of identified particles is an important observable when studying the characteristics of the QGP. However, since particles can scatter and be regenerated in between the chemical and kinetic freeze-out of a collision (the hadronic phase), information about the QGP phase imprinted in v n (p T ) can be altered by late-stage interactions and resonance decays, which can affect both v n and p T [45], leading to a deviation in mass ordering in v n (p T ) at low p T [46]. The φ -meson has been suggested as a particularly sensitive probe of the early collision phase as its production rate via regeneration in the hadronic phase is negligible [47] and it is theorized to have a low hadronic cross section [48][49][50], making it insensitive to the dissipative effects of the hadronic phase of the collision (although it should be noted that there is no consensus on the exact value of the cross section between the φ -meson and nucleons in heavy-ion collisions [51][52][53][54]). Recent experimental studies [27, 55,56] suggest that the φ -meson may be more sensitive to the hadronic phase than anticipated.
In this article, we present measurements of p T -differential elliptic, triangular, and quadrangular flow coefficients of π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson in Pb-Pb collisions at √ s NN = 5.02 TeV, extending greatly, and improving in precision upon, the previous measurements of identified particle v n in Pb-Pb collisions at √ s NN = 2.76 TeV as carried out by ALICE [27,28,33]. The results are reported for a wide range of particle transverse momenta within the rapidity range |y| < 0.5 at different collision centralities between 0-70% range. To isolate the fraction of anisotropic flow that is generated by initial-state fluctuations rather than geometry, the flow coefficients are also studied in ultra-central collisions (0-1% collision centrality). Centrality estimates the degree of overlap between the two colliding nuclei and is expressed as percentiles of the inelastic hadronic cross section, with low percentage values corresponding to head-on collisions. The measurements are performed using the scalar product method [57][58][59] with a (pseudo-)rapidity gap of |∆η| > 2.0 between the identified particles under study and the charged reference particles. The flow coefficients are measured separately for particles and anti-particles and are found to be compatible within the statistical uncertainties for most p T and centrality intervals. Any residual differences are included in the systematic uncertainties, and v n denotes the average between results for particles and anti-particles. This paper is organized as follows. Analysis details, particle identification, reconstruction methods, and flow measurement techniques are outlined in Sec. 2. The evaluation of systematic uncertainties is discussed in Sec. 3. The flow coefficients of π ± , K ± , p+p (v 2 , v 3 , and v 4 ), Λ+Λ, K 0 S (v 2 and v 3 ), and the φ -meson (v 2 ) are reported and compared to model calculations in Sec. 4. Finally, the results are summarized in Sec. 5.
2 Experimental setup and data analysis ALICE [60][61][62] is a dedicated heavy-ion experiment at the LHC optimized to study the properties of strongly interacting matter produced in heavy-ion collisions. A full overview of the detector layout and its performance can be found in [62,63]. The main subsystems used in this analysis are the Inner Tracking System (ITS) [64], Time Projection Chamber (TPC) [65], Time Of Flight detector (TOF) [66], and V0 [67]. The ITS, TPC, and TOF detectors cover full azimuth within pseudorapidity range |η| < 0.9 and lie within a homogeneous magnetic field of up to 0.5 T. The ITS consists of six layers of silicon detectors used for tracking and vertex reconstruction. The TPC is the main tracking detector and is also used to identify particles via specific ionization energy loss, dE/dx. The TOF in conjunction with the timing information from the T0 detector [68] provide particle identification based on flight time. The T0 is made up of two arrays of Cherenkov counters T0C and T0A, located at -3.3 < η < -3.0 and 4.5 < η < 4.9, respectively. Two scintillator arrays (V0), which cover the pseudorapidity ranges −3.7 < η < −1.7 (V0C) and 2.8 < η < 5.1 (V0A), are used for triggering, event selection, and the determination of centrality [69] and Q n -vectors (see Sec. 2.5). Both V0 detectors are segmented in four rings in the radial direction with each ring divided into eight sectors in the azimuthal direction. In addition, two tungsten-quartz neutron Zero Degree Calorimeters (ZDCs), installed 112.5 meters from the interaction point on each side, are used for event selection.

Event and track selection
The data sample recorded by ALICE during the 2015 LHC Pb-Pb run at √ s NN = 5.02 TeV is used for this analysis. The minimum-bias trigger requires signals in both V0A and V0C detectors. An offline event selection is applied to remove beam-induced background (i.e. beam-gas events) and pileup events. The former is rejected utilizing the V0 and ZDC timing information. The remaining contribution of such interactions is found to be smaller than 0.02% [63]. Pileup events, which constitute about 0.25% of recorded sample, are removed by comparing multiplicity estimates from the V0 detector to those of tracking detectors at mid-rapidity, exploiting the difference in readout times between the systems. The fraction of pileup events left after applying the dedicated pileup removal criteria is found to be negligible. The primary vertex position is determined from tracks reconstructed in the ITS and TPC as described in Ref. [63]. Only events with a primary vertex position within ±10 cm from the nominal interaction point along the beam direction are used in the analysis. Approximately 67 × 10 6 Pb-Pb events in the 0-70% centrality interval pass these selection criteria. Centrality is estimated from the energy deposition measured in the V0 detector [69].
Charged-particle tracks, used to measure the v n of π ± , K ± , p+p and the φ -meson, are reconstructed using the ITS and TPC within |η| < 0.8 and 0.5 < p T < 16.0 GeV/c with a track-momentum resolution better than 4% for the considered range [63]. Additional quality criteria are used to reduce the contamination from secondary charged particles (i.e., particles originating from weak decays, γ-conversions, and secondary hadronic interactions in the detector material) and fake tracks (random associations of space points). Only tracks with at least 70 space points, out of a maximum of 159, with a χ 2 per degree-offreedom for the track fit lower than 2, are accepted. Moreover, each track is required to cross at least 70 TPC pad rows and to be reconstructed from at least 80% of the number of expected TPC space points, in addition to having at least one hit in the two innermost layers of the ITS. Furthermore, tracks with a distance of closest approach (DCA) to the reconstructed event vertex smaller than 2 cm in the longitudinal direction (z) and (0.0105 + 0.0350 (p T c/GeV) −1.1 ) cm in the transverse plane (xy) are selected. Relevant selection criteria for tracks used for the reconstruction of K 0 S and Λ+Λ are given in Sec. 2.3.

Identification of π ± , K ± and p+p
Particle identification is performed using the specific ionization energy loss, dE/dx, measured in the TPC and the time of flight obtained from the TOF. A truncated-mean procedure is used to estimate the dE/dx (where the 40% highest-charge clusters are discarded), which yields a dE/dx resolution around 6.5% in the 0-5% centrality class [63]. At least 70 clusters are used for the dE/dx estimation. The TOF measures the time that a particle needs to travel from the primary vertex to the detector itself with a time resolution of ≈ 80 ps [63]. The start time for the TOF measurement is provided by the T0 detector or from a combinatorial algorithm which uses the particle arrival times at the TOF detector itself [63,66].
Expressing the difference between the expected dE/dx and the time of flight for π ± , K ± and p+p, and the measured signals in both TPC and TOF, in units of the standard deviations from the most probable value for both detectors (nσ TPC , nσ TOF ), and applying a selection on the number of accepted nσ , allows for particle identification on a track-by-track basis. The TPC dE/dx of different particle species are separated by at least 4σ for p T < 0.7 GeV/c, while in the relativistic rise region of the dE/dx (p T > 2 GeV/c) particle identification is still possible but only on a statistical basis [63]. The TOF detector provides 3σ separation between π ± and K ± for p T < 2.5 GeV/c, and between K ± and p+p for p T < 4 GeV/c [63].
The information from the TPC and TOF is combined using a quadratic sum nσ PID = nσ 2 TPC + nσ 2 TOF for 0.5 < p T ≤ 4 GeV/c. Particles are selected by requiring nσ PID < 3 for each species. The smallest nσ PID is used to assign the identity when the selection criterion is fulfilled by more than one species. When measuring p+p v n (p T ), only p are considered for p T < 2 GeV/c to exclude secondary protons from detector material and to reduce the contributions from weak decays. At high transverse momenta (p T > 4 GeV/c), K ± cannot reliably be identified. Identification of π ± and p+p for p T > 4 GeV/c is done utilizing the TPC dE/dx signal only. Pions (protons) are selected from the upper (lower) part of the expected pion (proton) dE/dx distribution. For example, proton selection typically varies in the range from 0 to −3σ TPC or from −1.5σ TPC to −4.5σ TPC depending on the momentum.
Secondary contamination from weak decays decreases from about 30% to 5% for p+p in the p T range 0.7-4.0 GeV/c and from about 5% to 0.5% for π ± in the p T range 0.5-4.0 GeV/c, while it is negligible for K ± .
The v n coefficients are not corrected for these contaminations; their effect on v n is at maximum ≈ 8%, for p + p v 2 at p T < 1 GeV/c for central collisions, and negligible for K ± , π ± v n . The contamination from other particle species is below 3% and 20% at p T > 4.0 GeV/c for π ± and p+p, respectively, and contamination from fake tracks is negligible. The v n results are reported for 0.5 < p T < 16.0(12.0, 6.0) GeV/c for π ± v 2 (v 3 , v 4 ), 0.7 < p T < 16.0(12.0, 6.0) GeV/c for p+p v 2 (v 3 , v 4 ), and 0.5 < p T < 4.0 GeV/c for K ± v n , all within |y| < 0.5.
2.3 Reconstruction of K 0 S and Λ + Λ The K 0 S and Λ+Λ are reconstructed in the K 0 S → π + + π − and Λ → p + π − (Λ → p + π + ) channels with branching ratios of 69.2% [70] and 63.9% [70] respectively. Reconstruction of K 0 S and Λ+Λ is based on identifying secondary vertices from which two oppositely-charged particles originate, called V 0 s. Topological selection criteria pertaining to the shape of the V 0 decay can be imposed, as well as requirements on the species identity of the decay products (called daughter particles).
The V 0 candidates are selected to have an invariant mass between 0.4 and 0.6 GeV/c 2 and 1.07 and 1.17 GeV/c 2 for K 0 S and Λ+Λ, respectively. The invariant mass of the V 0 is calculated based on the assumption that the daughter particles are either a π + π − pair, or a pπ − (pπ + ) pair. The daughter particles have been identified over the entire p T range using the TPC following the nσ approach detailed in Sec. 2.2 (|nσ TPC | < 3). The daughter tracks were reconstructed within |η| < 0.8, while the criteria on the number of TPC space points, the χ 2 per TPC space point per degree-of-freedom, the number of crossed TPC pad rows, and the percentage of the expected TPC space points used to reconstruct a track are identical to those applied for primary particles. In addition, the minimum DCA of daughter tracks to the primary vertex is 0.1 cm. Furthermore, the maximum DCA of daughter tracks to the secondary vertex is 0.5 cm to ensure that they are products of the same decay.
To reject secondary vertices arising from decays into more than two particles, the cosine of the pointing angle θ p is required to be larger than 0.998. This angle is defined as the angle between the momentumvector of the V 0 assessed at its decay position and the line connecting the V 0 decay vertex to the primary vertex and has to be close to 0 as a result of momentum conservation. In addition, the V 0 candidates are only accepted when they are produced at a distance between 5 and 100 cm from the nominal primary vertex in the radial direction. The lower value is chosen to avoid any bias from the efficiency loss when secondary tracks are being wrongly matched to clusters in the first layer of the ITS. To assess the systematic uncertainty related to contaminations from Λ+Λ and electron-positron pairs coming from γ-conversions to the K 0 S sample, a selection in the Armenteros-Podolanski variables [71] is applied for the K 0 S candidates, rejecting ones with q ≤ |α|/5. Here q is the momentum projection of the positively charged daughter track in the plane perpendicular to the V 0 momentum and α = (p + L − p − L )/(p + L + p − L ), with p ± L the projection of the positive or negative daughter tracks' momentum onto the momentum of the V 0 .
To obtain the p T -differential yield of K 0 S and Λ+Λ (which, together with background yields, are used for the v n extraction cf. Eq. 4), invariant mass distributions at various p T intervals are parametrized as a sum of a Gaussian distribution and a second-order polynomial function. The latter is introduced to account for residual contaminations (background yield) that are present in the K 0 S and Λ+Λ signals after the topological and daughter track selections. The K 0 S and Λ+Λ yields are extracted by integration of the Gaussian distribution. Obtained yields have not been corrected for feed-down from higher mass baryons (Ξ ± , Ω ± ) as earlier studies have shown that these have a negligible effect on the measured v n [27]. The v n (p T ) results are reported within |y| < 0.5 and 0.5 < p T < 10 GeV/c for K 0 S and 0.8 < p T < 10 GeV/c for Λ+Λ.

Reconstruction of φ -mesons
The φ -meson is reconstructed in the φ → K + +K − channel with a branching ratio of 48.9% [70]. Its reconstruction proceeds by first identifying all primary K ± tracks in an event, following the procedure for primary charged K ± outlined in Sec. 2.2. The K ± identification criterion nσ PID < 3 is chosen as it improves the significance of the φ -meson yield, while retaining a sufficient reconstruction efficiency. The vector sums of all possible K ± pairs are called φ -meson candidates, the yield of which is obtained as function of invariant mass M K + K − in various p T intervals. The p T -differential φ -meson yield is obtained by first subtracting a background yield from the candidate yield. This background yield is estimated using an event-mixing technique [72], in which K ± from different collisions are paired into background tracks, and is normalized to the candidate yield for 1.04 < M K + K − < 1.09 GeV/c. Collisions with similar characteristics (vertex position, centrality) are used for this mixing. To obtain the p T -differential yield of φ -mesons, the invariant mass distributions of the candidate yield is, after the aforementioned subtraction, parametrized as a sum of a Breit-Wigner distribution and a second-order polynomial function, the latter introduced to account for residual contaminations. The φ -meson yields are extracted by integration of the Breit-Wigner distribution and, together with background yields, used for the v n extraction (see Eq. 4).
The v 2 (p T ) results are reported for 0.9 < p T < 6.5 GeV/c within |y| < 0.5.

Flow analysis techniques
The flow coefficients v n are measured using the scalar product method [57][58][59], written as where u n,k = exp(inϕ k ) is the unit flow vector of the particle of interest k with azimuthal angle ϕ k , Q n is the event flow vector, and n is the harmonic number. Brackets · · · denote an average over all events, the double brackets · · · an average over all particles in all events, and * the complex conjugate.
The vector Q n is calculated from the azimuthal distribution of the energy deposition measured in the V0A. Its x and y components are given by where the sum runs over the 32 channels j of the V0A detector, ϕ j is the azimuthal angle of channel j defined by the geometric center, and w j is the amplitude measured in channel j. The vectors Q A n and Q B n are determined from the azimuthal distribution of the energy deposition measured in the V0C and the azimuthal distribution of the tracks reconstructed in the ITS and TPC, respectively. The amplitude measured in each channel of the V0C (32 channels as for the V0A) is used as weight in the case of Q A n , while unity weights are applied for Q B n . Tracks used for Q B n are selected following the procedure for primary charged tracks outlined in Sec. 2.1 for 0.2 < p T < 5.0 GeV/c. In order to account for a non-uniform detector response, the components of the Q n and Q A n vectors are recalibrated using a recentering procedure (i.e. subtraction of the Q n -vector averaged over many events from the Q n -vector of each event) [73]. The large gap in pseudorapidity between u n,k and Q n (|∆η| > 2.0) greatly suppresses short-range correlations unrelated to the azimuthal asymmetry in the initial geometry, commonly referred to as 'non-flow'. These correlations largely come from the inter-jet correlations and resonance decays. The v n of the K 0 S , Λ+Λ, and φ -meson cannot directly be measured using Eq. 2 as K 0 S , Λ+Λ and the φ -meson cannot be identified on a particle-by-particle basis. Therefore, the v tot n of V 0 s and φ -meson candidates is measured as function of both invariant mass, M d + d − , and candidate p T . This v tot n can be written [74] as the weighted sum of v n (p T ) of the particle of interest, v sig n , and that of background tracks, v bg where signal and background yields N sig and N bg are obtained for each p T interval from the K 0 S , Λ+Λ and φ -meson reconstruction procedures outlined in Secs 2.3 and 2.4. The formalism of Eq. 2 is used to

Systematic uncertainties
The systematic uncertainties on v n fall into the following categories: those arising from event selection, those arising from charged particle tracking, uncertainties in particle identification, uncertainties in V 0 finding, and those coming from the extraction of v n (p T ).
For p T ≤ 4 GeV/c, a p T -dependent systematic uncertainty is assigned to v 2 , v 3 , and v 4 of π ± , K ± , p+p, Λ+Λ, K 0 S and the φ -meson. Per measured point, the difference between the nominal measurement and a variation on the nominal measurement is calculated. If this difference between the nominal data point and the systematic variation is significant (where significance is evaluated based on the recommendations in [75]), it is considered to be a systematic uncertainty. When various checks are performed to quantify the effect of one systematic uncertainty (e.g. using three different centrality estimators to estimate the uncertainty in centrality determination), the maximum significant deviation that is found between the Table 1: Summary of systematic uncertainties for the v 2 of π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson. The uncertainties depend on p T and centrality range; minimum and maximum values are listed here. Empty fields indicate that a given check does not apply to the particle of interest. If an uncertainty has been tested but cannot be resolved within statistical precision, the field is marked negl for negligible. Horizontal lines are used to separate the different categories of systematic uncertainties as explained in Sec. 3. nominal measurement and the systematic variations is assigned as a systematic uncertainty. For each particle species, a p T -independent average uncertainty is reported for p T > 4 GeV/c in order to suppress sensitivity to statistical fluctuations. The uncertainty is obtained by fitting a zeroth-order polynomial to the significant p T -dependent relative uncertainties.
The systematic uncertainties are evaluated (if applicable) for each particle species, v n (p T ) and centrality intervals. A quadratic sum of the systematic uncertainties from the independent sources is reported as final systematic uncertainty on the measurements. An overview of the magnitude of the relative systematic uncertainties per particle species is given in Tabs. 1, 2, and 3 for v 2 , v 3 , and v 4 , respectively.

Event selection
The nominal event selection criteria and centrality determination are discussed in Sec. 2.1. Event selection criteria are varied by (i) changing the default centrality estimator from energy deposition in the V0 scintillator to either an estimate based on the number of hits in the first or second layer of the ITS; n fit ranges 0-2% 0-2% Table 2: Summary of systematic uncertainties for the v 3 of π ± , K ± , p+p, Λ+Λ, and K 0 S . The uncertainties depend on p T and centrality range; minimum and maximum values are listed here. Empty fields indicate that a given check does not apply to the particle of interest. If an uncertainty has been tested but cannot be resolved within statistical precision, the field is marked negl for negligible. Horizontal lines are used to separate the different categories of systematic uncertainties as explained in Sec. 3.
(ii) performing the v n analysis of π ± , K ± , and p+p in 1% wide centrality intervals to test the effect of multiplicity fluctuations (a test not possible for K 0 S , Λ+Λ v 3 ); (iii) not rejecting events with tracks caused by pileup or imposing a stricter than default pileup rejection by requiring a tighter correlation between the V0 and central barrel multiplicities; (iv) requiring the reconstructed primary vertex of a collision to lie alternatively within ±12 cm and ±5 cm from the nominal interaction point along the beam axis; (v) analyzing events recorded under different magnetic field polarities independently; (v) analyzing events recorded at different collision rates independently.

Charged particle tracking
The nominal charged particle track selection criteria are outlined in Sec. 2.1. Charged particle track selection criteria are varied by (i) requiring the third layer of the ITS to be part of the track reconstruction rather than the first two layers only; (ii) using only tracks that have at least three hits per track in the ITS, complemented by tracks without hits in the first two layers of the ITS (in which case the primary interaction vertex is used as an additional constraint for the momentum determination); (iii) changing the requirement on the minimum number of TPC space points that are used in the reconstruction from 70 to 60, 80, and 90; (iv) an additional systematic uncertainty is evaluated combining the following checks of the track quality: rejecting tracks that are reconstructed close to the sector boundaries of the TPC to which the sensitive pad rows do not extend, varying the minimum number of crossed TPC pad rows from 70 to 120, and requesting at least 90% instead of 80% of the expected TPC space points to reconstruct a track. Variations (i) and (ii) are referred to as tracking mode in Tabs. 1, 2, and 3.

Particle identification
The nominal particle identification approach for π ± , K ± , and p+p is outlined in Sec. 2.2. Particle identification criteria are varied by (i) changing the minimum number of clusters in the TPC that are used to estimate the dE/dx from 70 to 60, 80, and 90; (ii) rejecting tracks that satisfy the particle identification criterion for more than one particle species simultaneously for p T < 4 GeV/c; (iii) varying the particle identification criterion from nσ PID < 3 to nσ PID < 1, nσ PID < 2, and nσ PID < 4; (iv) varying the nσ TPC ranges that are used for particle identification for p T > 4 GeV/c. (viii) varying the particle identification criterion of the V 0 daughter tracks from |nσ TPC | < 3 to |nσ TPC | < 1 and |nσ TPC | < 4; (ix) changing the minimum value of the cos θ p from 0.998 to 0.98; (x) varying the minimum radial distance to the primary vertex at which the V 0 can be produced from 5 cm to 1 cm and 15 cm; (xi) varying the maximum radial distance to the beam pipe at which the V 0 can be produced from 100 cm to 50 cm and 150 cm; (xii) the contamination from Λ+Λ decays and γ-conversions to the K 0 S sample is checked by only selection V 0 daughter tracks with a dE/dx value 2σ away from the expected electron dE/dx, effectively excluding electrons, and limiting the value of the Armenteros-Podolanski variables α and q.
The yield extraction, as explained in Sec 2.3 for the K 0 S and Λ+Λ, and Sec 2.4 for the φ -meson, is varied by: (i) using a third-order polynomial as parametrization of residual background in the invariant mass spectra; (ii) using for the φ -meson a Voigtian distribution (a convolution of a Gaussian distribution and Breit-Wigner distribution, where the width of the Breit-Wigner distribution is set to the natural width of the φ -meson, allowing for the Gaussian distribution to describe the smearing of the φ -meson width due to the detector resolution) for the parametrization of the φ -meson invariant mass yield; using for the K 0 S and Λ+Λ a sum of two Gaussian distributions with the same mean for the parametrization of the K 0 S , Λ+Λ invariant mass yield; (iii, for the φ -meson only) using the yield of like-sign kaon pairs, in which two kaons with equal charge from the same event are used as candidate, for background yield description instead of event mixing.

Extraction of the v n (p T )
The nominal approach of measuring v n (p T ) is outlined in Sec. 2.5, and is varied by: (i) performing flow analysis for π ± , K ± , and p+p for positive and negative charges independently; (ii) performing flow analysis for positive and negative rapidities independently; (iii) performing flow analysis for π ± , K ± , and p+p in 1% centrality intervals and merging the result rather than measuring in wider centrality intervals directly; (iv) suppressing the signal from a specific V0A channel in the evaluation of the Q n -vector (see Eq. 3), which, on average, measures a lower energy deposition with respect to the ones reported by the other channels from the same ring; (v) performing flow analysis with the Q n -vector calculated from the V0A or V0C separately; (vi) testing various testing the assumption made on v bg n by changing the parametrization from a second-order polynomial to a first-order polynomial function.

Results and discussion
The flow coefficients v 2 , v 3 , and v 4 of identified particles are presented for various centrality classes in Sec. 4.1; scaling properties are discussed in Sec. 4.2. Comparisons to various model calculations, studies on the shape evolution of v n (p T ) with centrality, and comparisons to v n measured at √ s NN = 2.76 TeV are shown in Secs. 4.3, 4.4, and 4.5, respectively. Figure 2 shows the v 2 (p T ) of π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson for various centrality classes in the range 0-70%. For the π ± , K ± and p+p, measurements performed in ultra-central collisions (0-1%) are also presented. For the φ -meson, the results are reported in the 5-60% centrality range, where v 2 can be measured accurately. The magnitude of v 2 increases strongly with decreasing centrality up to the 40-50% centrality interval for all particle species. This evolution is expected, since the eccentricity of the overlap zone of the colliding nuclei increases for peripheral collisions and v 2 scales approximately linearly with eccentricity [76]. For more peripheral collisions (i.e. 50-60% and 60-70%), the value of v 2 is smaller than in the previous centrality intervals for all particle species except the φ -meson. This suggests that the system has a shorter lifetime in more peripheral collisions, which does not allow for the generation of large v 2 [77]. Furthermore, the reduced contribution of eccentricity fluctuations and hadronic interactions might play an important role in these centrality ranges [78]. A non-zero, positive v 2 is found in the 0-1% centrality interval for p T < 6 GeV/c for π ± , K ± , and p+p, which mostly reflects  the contribution from event-by-event fluctuations in the initial nucleon and gluon density as the system shape is almost spherical at vanishing impact parameter.

Centrality and p T dependence of flow coefficients
The third-order flow coefficent v 3 is generated by inhomogeneities in the initial nucleon and gluon density and not by the collision geometry [12][13][14][15], while v 4 arises from initial collision geometry, fluctuations, and the non-linear hydrodynamic response of the medium [16,17]. Higher-order flow harmonics are more sensitive to transport coefficients than v 2 [15], as the dampening effect of η/s leads to a stronger decrease of these coefficients [18,19]. Figures 3 and 4 present the v 3 (p T ) of π ± , K ± , p+p, Λ+Λ, and K 0 S and v 4 (p T ) of π ± , K ± , and p+p for various centrality classes in the 0-50% range. Non-zero, positive v 3 and v 4 are observed for particle species throughout the entire p T ranges up to ≈ 8 GeV/c. Unlike v 2 , the coefficients v 3 and v 4 increase weakly from ultra-central to peripheral collisions. This observation illustrates that higher-order flow coefficients are mainly generated by event-by-event fluctuations in the initial nucleon and gluon density. All flow coefficients increase monotonically with increasing p T up to 3-4 GeV/c where a maximum is reached. The position of this maximum depends on centrality and particle species as it takes place at higher p T for heavier particles for various centrality classes. This behaviour can be explained by the centrality dependence of radial flow combined with the parton density, which will be detailed in Sec. 4.4. Figure 5 presents the evolution of v n (p T ) of different particle species for various centrality classes. In the most central collisions, initial nucleon-density fluctuations are expected to be the main contributor to the generation of v n . For the 0-1% centrality interval, v 3 is the dominant flow coefficient for 1.5 < p T < 6.0 GeV/c, 2.0 < p T < 4 GeV/c, and 2.5 < p T < 6 GeV/c for π ± , K ± , and p+p, respectively. Furthermore, v 4 becomes equal to v 2 at p T ≈ 2.0 GeV/c (2.2, 2.5) for π ± (K ± , p+p), after which it increases gradually and reaches a magnitude similar to v 3 at around 3.5 GeV/c. A similar trend is observed in the 0-5% centrality class for all particle species. However, the crossing between flow coefficients (the p T value at which they reach a similar magnitude), which also depends on the particle mass, takes place at different p T values than for the 0-1% centrality interval. This dependence of the crossing between different flow coefficients can be attributed to the interplay of elliptic, triangular, and quadrangular flow with radial flow. Upwards of 5% collision centrality, v 2 is larger than v 3 and v 4 , confirming the hypothesis that collision geometry dominates the generation of flow coefficients. Figure 6 shows the v 2 (p T ) of π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson in a given centrality interval arranged into panels of various centrality classes, which allows for further illustration of the interplay between elliptic and radial flow. For p T < 2-3 GeV/c, v 2 of the different particle species is mass-ordered, meaning that lighter particles have a larger v 2 than heavier particles at the same p T . This behaviour is indicative of strong radial flow which imposes an equal, isotropic velocity boost to all particles in addition to the anisotropic expansion of the medium [20][21][22]. For 3 < p T < 8-10 GeV/c, particles are grouped according to their number of constituent quarks, which supports the hypothesis of particle production via quark coalescence [34]. Particle type scaling and mass ordering are most directly tested by the φ -meson v 2 , as its mass is close to the proton mass. Figure 6 demonstrates that the φ -meson v 2 follows proton v 2 at low p T , but pion v 2 at intermediate p T in all centrality classes. The crossing between meson and baryon v 2 , which depends on the particle species, happens at higher p T values for central than peripheral collisions as a result of the larger radial flow in the former. Lastly, it is seen that the v 2 of baryons is higher than that of mesons up to p T ≈ 10 GeV/c, indicating that particle type dependence of v 2 persists up to high p T . For p T > 10 GeV/c, where v 2 depends only weakly on transverse momentum, the magnitude of p+p v 2 is compatible with that for π ± within statistical and systematic uncertainties. Furthermore, the nuclear modification factor in this high p T region is found to be the same for the two particle species within uncertainties [79].    The p T -differential v 3 of π ± , K ± , p+p, Λ+Λ, and K 0 S for various centrality classes. Statistical and systematic uncertainties are shown as bars and boxes, respectively.
Both v 3 and v 4 show a clear mass ordering at p T < 2-3 GeV/c, confirming the interplay between triangular and quadrangular flow and radial flow. For 3 < p T < 8 GeV/c, particles are grouped into mesons and baryons and, analogous to the trend of v 2 in this p T region, the flow of baryons is larger than that of mesons. The crossing between meson and baryon v 3 and v 4 also exhibits a centrality and particle mass dependence. Figures 6 and 7 also show a comparison between K ± and K 0 S v 2 and v 3 as a function of p T for various centrality classes. A difference in v n (p T ) is found between the K ± and K 0 S measurements: the magnitude of K 0 S v n is systematically smaller than the magnitude of K ± v n . This difference in v n exhibits no p T dependence, but changes with centrality for v 2 . For 0.8 < p T < 4.0 GeV/c, the difference in v 2 ranges from 7% ± 3.5%(syst) ± 0.7%(stat) in the most central collisions to 1.5% ± 1.5%(syst) ± 0.4%(stat) in peripheral collisions. In the same kinematic range, a deviation in v 3 of 6.5% ± 5%(syst) ± 1.7%(stat) is found in the most central collisions and of 6% ± 4.5%(syst) ± 1%(stat) in peripheral collisions. This difference is similar in magnitude and centrality dependence as the one reported by ALICE in Pb-Pb collisions at √ s NN = 2.76 TeV in [27].

Scaling properties
To test the hypothesis of particle production via quark coalescence [34], which would lead to a grouping of v n of mesons and baryons at intermediate p T , both v n and p T are divided by the number of constituent quarks (n q ) independently for each particle species. The v 2 /n q , v 3 /n q , and v 4 /n q of π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson, plotted as a function of p T /n q , are reported in Figs. 9, 10, and 11 for various centrality classes.
For p T /n q > 1 GeV/c, the scaling is only approximate. To quantify the degree to which the measurements deviate from the n q scaling, the p T /n q dependence of v n /n q has been divided by a cubic spline fit to the p+p v n /n q . In the region where quark coalescence is hypothesized to be the dominant process (≈ 1 < p T /n q < 3 GeV/c) [34,80], a deviation from the exact scaling of ± 20% is found for v 2 for central collisions, which decreases to ±15% for the most peripheral collisions. For higher harmonics, a ±20% deviation is found for all centrality classes. This deviation is in agreement with earlier observations [27, 28, 32].

Comparison with model calculations
To test the validity of the hydrodynamic description of the QGP evolution, the v n measurements in the 0-5%, 10-20% and 40-50% centrality intervals are compared to hydrodynamical calculations in Figs. 12, 13, and 14 for π ± , K ± , and p+p, respectively. Predictions from MUSIC [81] and iEBE-VISHNU [82] simulations are depicted by the different coloured curves. The first calculation is based on MUSIC [83], an event-by-event 3+1 dimensional viscous hydrodynamic model, coupled to a hadronic cascade model (UrQMD) [84,85], which allows the influence of the hadronic phase on the anisotropic flow to be studied for different particle species. The IP-Glasma model [86,87] is used to simulate the initial conditions of the collision. MUSIC uses a starting time for the hydrodynamic evolution of  The p T /n q dependence of v 3 /n q of π ± , K ± , p+p, Λ+Λ, and K 0 S for various centrality classes. Statistical and systematic uncertainties are shown as bars and boxes, respectively. τ 0 = 0.4 fm/c, a switching temperature between the macroscopic hydrodynamic description and the microscopic transport evolution of T sw = 145 MeV, a value of η/s = 0.095, and a temperature dependent ζ /s. The second calculation employs the iEBE-VISHNU hybrid model [88], which is an event-byevent version of the VISHNU hybrid model [89], and couples 2+1 dimensional viscous hydrodynamics VISH2+1 [77] to UrQMD. The T R ENTo [90] and AMPT [91] models are used to describe the initial conditions. For both configurations, τ 0 = 0.6 fm/c and T sw = 148 MeV are set from [92], where these values have been obtained utilizing Bayesian statistics from a simultaneous fit of final charged-particle density, mean transverse momentum, and integrated flow coefficients v n in Pb-Pb collisions at √ s NN = 2.76 TeV.   The p T -differential v 2 (top), v 3 (middle), and v 4 (bottom) of π ± for the 0-5%, 10-20%, and 40-50% centrality classes compared to hydrodynamical calculations from MUSIC model using IP-Glasma initial conditions (magenta) [81] and the iEBE-VISHNU hybrid model using AMPT (orange) or T R ENTo (cyan) initial conditions [82]. Statistical and systematic uncertainties of the data points are shown as bars and boxes, respectively. The uncertainties of the hydrodynamical calculations are depicted by the thickness of the curves. The ratios of the measured v n to a fit to the hydrodynamical calculations are also presented for clarity.  Fig. 13: (Colour online) The p T -differential v 2 (top), v 3 (middle), and v 4 (bottom) of K ± for the 0-5%, 10-20%, and 40-50% centrality classes compared to hydrodynamical calculations from MUSIC model using IP-Glasma initial conditions (magenta) [81] and the iEBE-VISHNU hybrid model using AMPT (orange) or T R ENTo (cyan) initial conditions [82]. Statistical and systematic uncertainties of the data points are shown as bars and boxes, respectively. The uncertainties of the hydrodynamical calculations are depicted by the thickness of the curves. The ratios of the measured v n to a fit to the hydrodynamical calculations are also presented for clarity.
The temperature-dependent η/s and ζ /s extracted in [92] are used for T R ENTo initial conditions, while η/s = 0.08 and ζ /s = 0 are taken for AMPT. Figures 12,13,and 14 show that the hydrodynamical calculations qualitatively reproduce the v n measurements. The differences between the data points and models are visualized in Figs. 12, 13, and 14 as the ratios of the measured v n to a fit to the theoretical calculations. The iEBE-VISHNU calculations using AMPT initial conditions describe the p T -differential v n of π ± , K ± , and p+p more accurately than T R ENTo based and MUSIC calculations for p T > 1 GeV/c. Using AMPT initial conditions, there is good agreement between π ± and K ± v n and iEBE-VISHNU calculations for p T < 2 GeV/c, while p+p v n is described fairly well up to p T = 3 GeV/c. The T R ENTo based predictions follow π ± and K ± v n up to slightly lower transverse momenta (p T <1-2 GeV/c) and to p T < 3 GeV/c) for p+p, depending on the considered centrality interval. The MUSIC calculations are in agreement with the measured v n for p T < 1 GeV/c in central collisions, however they overestimate v 2 at lower p T in more peripheral collisions. The p T -differential v 2 (top), v 3 (middle), and v 4 (bottom) of p+p for the 0-5%, 10-20%, and 40-50% centrality classes compared to hydrodynamical calculations from MUSIC model using IP-Glasma initial conditions (magenta) [81] and the iEBE-VISHNU hybrid model using AMPT (orange) or T R ENTo (cyan) initial conditions [82]. Statistical and systematic uncertainties of the data points are shown as bars and boxes, respectively. The uncertainties of the hydrodynamical calculations are depicted by the thickness of the curves. The ratios of the measured v n to a fit to the hydrodynamical calculations are also presented for clarity.

Shape evolution of v n (p T ) as function of centrality
The evolution of the shape of v n (p T ) as function of centrality is quantified by taking the ratio of v n (p T ) in a given centrality interval to the v n (p T ) measured in the 20-30% centrality interval where the second fraction on the right-hand side of the equation serves as a normalization factor which is constructed from the p T -integrated v n values obtained in the 20-30% centrality interval (v n | 20−30% ) and the centrality interval of interest (v n ). Centrality-dependent variations in the shape of v n (p T ) will present themselves as deviations from unity of the observed v n (p T ) ratio to 20−30% .
The shape evolution of elliptic and triangular flow is shown in Figs. 15 and 16 for π ± , K ± , p+p, and inclusive charged particles (the latter taken from [44]). For inclusive charged particles, variations in shape of about 10% are observed for p T < 3 GeV/c, which increase to about 30% for p T < 6 GeV/c. The shape evolution of v 2 (p T ) shows different trends for π ± , K ± , and p+p. While π ± v 2 (p T ) ratio to 20−30% follows inclusive charged particle over the considered p T range, the elliptic flow of p+p (K ± ) varies between 20% (10%) to 250% (55%) at low p T from most central to peripheral collisions. The variations are more pronounced for v 3 (p T ), in particular for central collisions. The mass dependence found in the shape evolution of both v 2 and v 3 for p T < 4 GeV/c can be attributed to variations of the magnitude 60-70% Fig. 15: (Colour online) Centrality dependence of v 2 (p T ) ratio to 20−30% for π ± , K ± , p+p, and inclusive charged particles [44]. Statistical and systematic uncertainties are shown as bars and boxes, respectively. 40-50% Fig. 16: (Colour online) Centrality dependence of v 3 (p T ) ratio to 20−30% for π ± , K ± , p+p, and inclusive charged particles [44]. Statistical and systematic uncertainties are shown as bars and boxes, respectively.
of radial flow and quark density, both being larger for central than peripheral collisions. Radial flow has a stronger effect on the v n of heavier particles than that of lighter particles at low p T , while the quark density influences the peak value of v n (p T ) in the coalescence model picture [35,36,93]. For p T > 4 GeV/c, the shape evolution shows little (if any) particle type dependence.
The shape evolution of v 2 (p T ) for π ± , K ± , and p+p is compared to calculations from the MUSIC and 40-50% Fig. 17: (Colour online) Centrality dependence of v 2 (p T ) ratio to 20−30% for π ± (upper panels), K ± (middle panels), and p+p (lower panels) compared to hydrodynamical calculations from the MUSIC model using IP-Glasma initial conditions (magenta) [81], the iEBE-VISHNU hybrid model using AMPT (orange) or T R ENTo (cyan) initial conditions [82]. Statistical and systematic uncertainties of the data points are shown as bars and boxes, respectively.
iEBE-VISHNU hybrid models in Fig. 17. Both models describe the shape evolution for p+p over the p T range 0.7 < p T < 3 GeV/c. The iEBE-VISHNU model reproduces the shape evolution for π ± and K ± for p T < 1.5 GeV/c. Calculations from the MUSIC model deviate strongly from the observed shape evolution for π ± and K ± in peripheral collisions.
As quark density depends on centrality, the maximum v n is expected to be found at higher p T in more central collisions. To further quantify this aspect of the shape evolution of v n (p T ), the p T of π ± , p+p, Λ+Λ, and K 0 S where v 2 (p T ) and v 3 (p T ) reach a maximum, divided by number of constituent quarks n q , is reported in Fig. 18 as a function of centrality. The φ -meson and K ± are not included since the kinematic range and granularity of the measurements do not allow for a reliable extraction of a maximum. The left panel of Fig. 18 shows that the p T /n q at which v 2 (p T ) reaches a maximum, denoted as p T | v max 2 /n q , decreases with increasing centrality percentile for collision centralities larger than 5-10%, following the expectations from the hypothesis of hadronization through coalescence. The systematic uncertainties as presented in Fig. 18 have been evaluated directly on p T | v max n /n q to accurately take into account that some systematic uncertainties can be point-by-point correlated in p T . In the 0-5% centrality interval, there is a hint of a lower p T | v max 2 /n q than in the 5-10% centrality class for all particle species. The observed p T | v max 2 /n q is compatible among all particle species with the exception of the p+p p T | v max 2 /n q , which is slightly lower in the 0-20% centrality range. The right panel of Fig. 18 presents p T | v max 3 /n q , which shows, within the large uncertainties, a weak (if any) centrality dependence for π ± and K 0 S and no centrality dependence for p+p and Λ+Λ. The p T | v max 3 /n q is the same for the different particle species within uncertainties.
In the scenario of ideal hydrodynamics, v n is a power law function of the radial expansion velocity of the medium [94,95] so that v n ∝ p T n up to p T ∼ M for particles with mass M. Figure 19 shows |v n | 1/n /p T   : (Colour online) Centrality dependence of |v n | 1/n /p T of inclusive charged particles [44], π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson for n = 2 (upper panels) and n = 3 (lower panels). Statistical and systematic uncertainties are shown as bars and boxes, respectively. as function of p T for n = 2 and n = 3 in various centrality intervals for inclusive charged particles [44], π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson (n = 2 only). When v n ∝ p T n , the observable |v n | 1/n /p T should be a constant. For π ± and the inclusive charged particles, the v n ∝ p T n scaling is broken both for v 2 and v 3 for all centrality intervals, as is also hypothesized in [96]. It should be noted however that the kinematic constraints imposed on the measurement preclude testing the scaling hypothesis in the full relevant momentum region. The scaling holds up to p T ≈ 1 GeV/c for K ± and K 0 S , and up to p T ≈ 2 GeV/c for p+p, Λ+Λ, and the φ -meson for the 0-5% and 10-20% centrality intervals. Similar qualitative observations are found in the three hydrodynamical calculations [81,82].
If v n indeed exhibits a power law dependence on p T n , ratios of the form of v 1/n n /v 1/m m are p T -independent. Previous measurements at RHIC [97,98] and the LHC [99,100] have shown that the ratios v  be due to fluctuations in the initial geometry [98]. The ratios v 3 /|v 2 | 3/2 , v 4 /|v 2 | 4/2 , and v 4 /|v 3 | 4/3 , which probe the same scaling but are in practice more sensitive, are shown in Figs. 20, 21, and 22, respectively. For each figure, v n /|v m | n/m is shown for inclusive charged particles [44], π ± , K ± and p+p in various centrality intervals. For v 3 /|v 2 | 3/2 and v 4 /|v 2 | 4/2 , no obvious p T dependence is found for inclusive charged particles between 5-50% collision centrality. For the 0-5% centrality class, the ratios are flat for p T < 3 GeV/c and rise monotonically for higher momenta. No particle type dependence of the ratios is found for p T > 1.5 GeV/c, below which the ratios for p+p v n rise. This rise of the p+p v n ratios can be attributed to an increase of radial flow which affects the independent harmonics differently. For the ratio v 4 /|v 3 | 4/3 , no p T dependence is observed over the full centrality range. Large statistical uncertainties do not allow conclusions to be drawn on the behaviour of p+p v n in the v 4 /|v 3 | 4/3 ratio. The transport properties and initial condition models can be further constrained by studying the energy dependence of anisotropic flow. Figure 23 presents the v 2 (p T ), v 3 (p T ), and v 4 (p T ) of π ± , K ± , and p+p compared to ALICE measurements performed at √ s NN = 2.76 TeV [28].
The v n coefficients at √ s NN = 2.76 TeV have been measured using the scalar product method, taking the particle of interest under study and the charged reference particles from different, non-overlapping pseudorapidity regions between |η| < 0.8. Assuming no anisotropic flow in minimum bias pp collisions at the same collision energy, the non-flow contributions are estimated from minimum bias pp collisions and subtracted from the measured v n coefficients. Ratios of the measurements presented in this paper to a cubic spline fit to the ones performed at √ s NN = 2.76 TeV are given in the figure for each presented centrality interval and flow coefficient. The uncertainties in these ratios are obtained by summing the statistical and systematic uncertainties on the independent measurements in quadrature, and propagating the obtained uncertainties as uncorrelated.
An increase of radial flow with increasing collision energy is expected to lead to a suppression of v n at low p T , an effect which would be most pronounced for heavier particles. Although a possible suppression of p+p v n at √ s NN = 5.02 TeV can be seen between 1 p T 3 GeV/c in central collisions and additionally for v 2 (p T ) of π ± and K ± at the same centrality interval, the precision of the results does not allow for conclusions to be drawn as the measurements at different collision energies are compatible within uncertainties. Figure 24 shows the v 2 (p T ) of Λ+Λ, K 0 S , and the φ -meson compared to ALICE measurements performed at √ s NN = 2.76 TeV [27], where the v 2 coefficients at √ s NN = 2.76 TeV have been measured using the scalar product method with an |∆η| > 0.9 gap to suppress non-flow. No differences are observed between the K 0 S and Λ+Λ v 2 (p T ) measured at two different collision energies. The strongly improved precision of the φ -meson measurement at √ s NN = 5.02 TeV, both in terms of statistical uncertainty and granularity in p T , shows that the v 2 (p T ) follows a mass ordering at low p T and groups with mesons after p T ≈ 3 GeV/c for all centrality intervals.

Summary
In summary, the elliptic, triangular, and quadrangular flow coefficients of π ± , K ± , p+p, Λ+Λ, K 0 S , and the φ -meson have been measured in Pb-Pb collisions at √ s NN = 5.02 TeV over a broad range of transverse momentum and in various centrality ranges. The precision of these measurements provide constraints for initial-state fluctuations and transport coefficients of the medium. The magnitude of v n increases with decreasing centrality up to the 40-50% centrality interval for all particle species. This increase is stronger for v 2 than for v 3 and v 4 , which indicates that collision geometry dominates the generation of elliptic flow while higher flow coefficients are mainly generated by event-by-event fluctuations in the initial nucleon and gluon densities. This interpretation is also supported by the non-zero, positive v n found in the 0-1% centrality interval. In most central collisions (i.e. 0-1% and 0-5%), v 3 and v 4 reach a similar magnitude as v 2 at different p T values depending on particle mass, after which they increase gradually. For p T < 3 GeV/c, the v n coefficients show a mass ordering consistent with an interplay between anisotropic flow and the isotropic expansion (radial flow) of the collision system. In this transverse momentum range, the iEBE-VISHNU hydrodynamical calculations describe the measured v n of π ± , K ± , and p+p fairly well for p T < 2.5 GeV/c, while MUSIC reproduces the measurements for p T < 1 GeV/c. It should be noted that neither of the presented hydrodynamical models is able to fully describe the measurements. At intermediate transverse momenta (3 < p T < 8-10 GeV/c), particles show an approximate grouping by the number of constituent quarks at the level of ±20% for all flow coefficients in the 0-50% centrality range. The φ -meson v 2 , which tests both particle mass dependence and type scaling, follows p+p v 2 at low p T and π ± v 2 at intermediate p T . The baryon v n has a magnitude larger than that of mesons for p T < 8-10 GeV/c, indicating that the particle type dependence persists up to high p T . For p T > 10 GeV/c, the v 2 of p+p is compatible with that of π ± within uncertainties. The shape evolution of v 2 (p T ) as function of centrality shows different trends for π ± , K ± , and p+p and varies between 20% (10%) to 250% (55%) for p+p (K ± ) at low p T from most central to peripheral collisions; variations are more pronounced for v 3 (p T ), in particular for central collisions. Ratios v 3 /|v 2 | 3/2 and v 4 /|v 2 | 4/2 are flat for p T < 3 GeV/c and rise monotonically for higher momenta for the 0-5% centrality class. No particle type dependence of the ratios is found for p T > 1.5 GeV/c, below which the ratios for p+p v n rise, which can be attributed to an increase of radial flow which affects the independent harmonics differently. For the ratio v 4 /|v 3 | 4/3 , no p T dependence is observed over the full centrality range. The measurements are compatible with those performed in Pb-Pb collisions at √ s NN = 2.76 TeV within uncertainties.