Anisotropic flow and flow fluctuations of identified hadrons in Pb$-$Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.02 TeV

The first measurements of elliptic flow of $\pi^\pm$, ${\rm K}^\pm$, p+$\overline{\rm p}$, ${\rm K_{S}^0}$, $\Lambda$+$\overline{\Lambda}$, $\phi$, $\Xi^-$+$\Xi^+$, and $\Omega^-$+$\Omega^+$ using multiparticle cumulants in Pb$-$Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV are presented. Results obtained with two- ($v_2\{2\}$) and four-particle cumulants ($v_2\{4\}$) are shown as a function of transverse momentum, $p_{\rm T}$, for various collision centrality intervals. Combining the data for both $v_2\{2\}$ and $v_2\{4\}$ also allows us to report the first measurements of the mean elliptic flow, elliptic flow fluctuations, and relative elliptic flow fluctuations for various hadron species. These observables probe the event-by-event eccentricity fluctuations in the initial state and the contributions from the dynamic evolution of the expanding quark-gluon plasma. The characteristic features observed in previous $p_{\rm T}$-differential anisotropic flow measurements for identified hadrons with two-particle correlations, namely the mass ordering at low $p_{\rm T}$ and the approximate scaling with the number of constituent quarks at intermediate $p_{\rm T}$, are similarly present in the four-particle correlations and the combinations of $v_2\{2\}$ and $v_2\{4\}$. In addition, a particle species dependence of flow fluctuations is observed that could indicate a significant contribution from final state hadronic interactions. The comparison between experimental measurements and CoLBT model calculations, which combine the various physics processes of hydrodynamics, quark coalescence, and jet fragmentation, illustrates their importance over a wide $p_{\rm T}$ range.


Introduction
The primary goal of the ultra-relativistic heavy-ion collision programme at the Large Hadron Collider (LHC) is to study the properties of the quark-gluon plasma (QGP), a novel state of strongly interacting matter at high temperatures and energy densities [1,2]. Studies of the azimuthal anisotropy of particle production have contributed significantly to the characterization of the system created in heavy-ion collisions [3,4]. Anisotropic flow reflects the conversion of the initial state spatial anisotropy into final state anisotropies in momentum space. This translation is facilitated by interactions between the constituents of the quark-gluon plasma (QGP) [5][6][7] and at later stages, after hadronisation, between the produced particles. Anisotropic flow is quantified by studying the azimuthal distribution of particles emitted in the plane transverse to the beam direction [3]. This is usually expressed in terms of a Fourier series in the azimuthal angle ϕ [8,9] according to E d 3 N dp 3 = 1 2π where E, N, p, p T , ϕ, and η are the energy, yield, momentum, transverse momentum, azimuthal angle, and pseudorapidity of particles, respectively, and Ψ n is the azimuthal angle of the symmetry plane of order n [10,11]. The v n coefficients are given by where ⟨⟩ denote an average over all particles in a single event. The second Fourier coefficient, v 2 , is usually referred to as elliptic flow. It is the dominant harmonic in heavy-ion collisions with large values of impact parameter (i.e. non-central collisions). Its value is sensitive to some of the basic transport coefficients of the QGP, e.g. the shear viscosity over entropy density ratio (η/s). The study of elliptic flow has been instrumental in establishing the strongly-coupled QGP paradigm, first in collisions at the Relativistic Heavy Ion Collider (RHIC) [12][13][14][15] and, since 2010, in collisions at the Large Hadron Collider (LHC) [16][17][18].
Around the start of the LHC heavy-ion program, it was realised that elliptic flow is also a sensitive probe of the initial state of heavy-ion collisions [19]. Its magnitude fluctuates from one event to the other, reflecting the event-by-event fluctuating energy-density profiles of the nuclear overlap region prior to the formation of the QGP. Initial event-by-event geometry fluctuations lead to the fluctuations of evenharmonic anisotropic flow and generate non-zero odd harmonics. In fact, the initial geometry fluctuations lead to ⟨v k n ⟩ ̸ = ⟨v n ⟩ k and the development of different order symmetry planes Ψ n in different kinematic regions in p T or η [20][21][22]. Thus, a comprehensive investigation of the final state flow fluctuations is crucial for understanding the event-by-event initial geometry fluctuations and their impact on the system dynamic evolution. Studies of charged particles in Pb-Pb collisions at LHC energies indicated non-Gaussian initial state fluctuations and, consequently, made it possible to constrain their probability distribution function (p.d.f.) [23,24]. Studies of flow fluctuations have so far been performed both experimentally and in theoretical model calculations for the measurements integrated over a large kinematic range [25][26][27][28]. On the other hand, a p T -differential study, and in particular with identified hadrons, has not been done before. These studies can provide insights on the interplay between the expansion of the system and its late-stage, highly-dissipative hadronic phase, as well as particle production mechanisms.
Similar studies in the past for various flow coefficients have been pivotal in establishing the need to include viscous corrections in hydrodynamic models and, consequently, in constraining the value of η/s to be very close to the conjectured lower limit of 1/4π calculated for infinitely strongly coupled gauge theories via the AdS/CFT correspondence [29]. Detailed studies of how anisotropic flow develops for different particle species as a function of p T for various centrality intervals (i.e. an estimate of the degree of overlap between the two colliding nuclei) of Pb-Pb collisions at LHC energies [26][27][28] confirmed a number of qualitative features already observed at RHIC [12][13][14][15]: the mass ordering of v n sis for p T > 2 GeV/c [42]. The Time of Flight detector (TOF) [44] is located around the TPC and is used for particle identification by measuring the flight time of particles from the collision point with a resolution of about 80 ps [42]. The start time for the TOF measurement is provided by the T0 detector with a resolution of about 25 ps [42,45], two arrays of Cherenkov counters covering the pseudorapidity ranges −3.3 < η < −3.0 (T0C) and 4.6 < η < 4.9 (T0A), or from a combinatorial algorithm that uses the particle arrival times at the TOF detector itself [42,44]. Both methods of estimating the start time are fully efficient for the 60% most central Pb-Pb collisions. The TOF provides a 3σ separation between π ± -K ± and K ± -p+p up to p T = 2.5 GeV/c and p T = 4 GeV/c, respectively [42]. The ITS, TPC, and TOF detectors cover the full azimuth within |η| < 0.9.
In the forward region, two scintillator arrays (V0) [46] are used for triggering, event selection, and the determination of the collision centrality [47]. The V0 consists of two systems, the V0C and V0A, positioned at −3.7 < η < −1.7 and 2.8 < η < 5.1, respectively. 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.
More details on the ALICE setup and the performance of the detectors can be found in Refs. [41,42].
3 Analysis procedure

Event and track selection
The data sample used in this analysis consists of Pb-Pb collisions at √ s NN = 5.02 TeV recorded by the ALICE detector in the years 2015 and 2018 LHC data-taking campaigns. A minimum bias trigger was provided by requiring signals in both V0A and V0C scintillator arrays. In addition, the sample of semicentral collisions was enhanced by an online selection based on the V0 signal amplitudes. Beam-induced background events (i.e. beam-gas interactions) were removed offline utilizing the V0 and ZDC timing information. Pileup of collisions from different bunch crossings in the TPC was rejected by comparing multiplicity estimates from the V0 detector to those of tracking detectors at midrapidity, exploiting the difference in readout times between the systems. The primary vertex position, determined from tracks reconstructed in the ITS and TPC, was required to be within ±10 cm from the nominal interaction point along the beam direction. These selection criteria were met by approximately 245 million events in the 10-60% centrality interval. The collision centrality was estimated from the amplitudes of the signals measured in the V0 detector [47].
Charged-particle tracks, used to measure the v 2 of π ± , K ± , p+p, φ -mesons, and inclusive charged particles, were reconstructed using the ITS and TPC within |η| < 0.8 and 0.2 < p T < 10 GeV/c. Each track was required to have a minimum number of 70 TPC space points (out of a maximum of 159) with a χ 2 per TPC space point lower than 4 and at least 2 hits in the ITS with a χ 2 per ITS hit smaller than 36. Moreover, tracks with a distance of closest approach (DCA) larger than 2 cm in the longitudinal direction were rejected. In the transverse plane, a p T -dependent DCA selection of the form 0.0105 + 0.0350 p −1.1 T cm was applied. These selection criteria lead to an efficiency of about 80% for primary tracks at p T > 0.5 GeV/c and contamination from fake tracks (random associations of space points) and secondary charged particles (i.e. particles originating from weak decays, conversions, and secondary hadronic interactions in the detector material) of about 5% at p T ≈ 1 GeV/c.
3.2 Selection of π ± , K ± , and p+p Particle identification of π ± , K ± , and p+p is performed using the dE/dx from the TPC and the time of flight from the TOF system, if available. The identification is based on the normalised difference between the measured and the expected signal for a given species (σ TPC and σ TOF , respectively). It uses the correlation between nσ TPC and nσ TOF in a Bayesian approach [48], where the signals converted into probabilities are folded with the expected abundances (priors) of each particle species. The minimal probability threshold has been set to 0.95 for π ± and 0.85 for K ± and p+p. In addition, particles are selected by requiring |nσ TPC | < 3 and |nσ TOF | < 3 for each species in the whole p T range. This procedure ensures a high purity of the studied sample, thus reducing the uncertainties due to particle misidentification. The resulting purity, estimated using Monte Carlo (MC) simulations, is higher than 95% for π ± for 0.2 < p T < 10 GeV/c, above 80% for K ± for 0.3 < p T < 6 GeV/c, and reaches values larger than 90% for p+p for 0.5 < p T < 6 GeV/c.

Reconstruction of φ mesons
The φ meson is reconstructed in the decay channel φ → K + + K − with a branching ratio of 49.2% [49]. Its decay products are selected using the same criteria for primary K ± (see Sec. 3.2). The φ meson yield is obtained from the invariant mass (M K + K − ) reconstructed from all possible K ± pairs after subtracting the combinatorial background evaluated using the like-sign kaon pairs in each p T and centrality interval. The resulting M K + K − distribution is parametrised as a sum of a Breit-Wigner (BW) distribution and a third-order polynomial function that accounts for residual contamination within the invariant mass range of 0.99 < M K + K − <1.07 GeV/c 2 . The p T -differential yield of φ mesons is extracted by integration of the BW distribution and used for the v 2 extraction together with the background yield (see Eq. 16). The procedure of the reconstruction of φ meson is identical to the previous measurements [28], while the extraction of v 2 {4}(p T ) is slightly different and explained in Section 3.7.

Reconstruction of K 0
S and Λ+Λ The reconstruction of K 0 S and Λ+Λ is based on identifying their secondary vertices called V 0 s in the decay channels K 0 S → π + + π − and Λ → p + π − (Λ → p + π + ) with branching ratios of 69.2% and 63.9% [49], respectively. Selection criteria related to the distinctive V-shaped decay topology and requirements on the characteristics of the daughter particles are applied to suppress the large combinatorial background. The invariant mass is calculated assuming that the daughter particles, identified using the TPC (|nσ TPC | < 3) over the entire p T range, are either a π + π − pair or a pπ − (pπ + ) pair. The V 0 candidates are selected with an invariant mass between 0.4 and 0.6 GeV/c 2 for K 0 S and 1.08 and 1.16 GeV/c 2 for Λ+Λ. The daughter tracks are reconstructed within |η| < 0.8 using the same TPC track quality requirements described in Section 3.1 for charged tracks. In addition, the ratio between the number of space points and the number of crossed rows in the TPC is required to be larger than 0.8, the minimum DCA of daughter tracks to the primary vertex is 0.1 cm, and the maximum DCA of daughter tracks to the secondary vertex is 0.5 cm. Only V 0 candidates produced at a radial distance between 5 and 100 cm from the beam line and with a cosine of the pointing angle (the angle between the line connecting the primary and V 0 vertices and the V 0 momentum vector) larger than 0.998 are accepted. To reduce the contamination from Λ+Λ and electron-positron pairs coming from γ conversions, an additional selection is applied in the Armenteros-Podolanski variables [50] of the K 0 S candidates, similar to what is done in Ref. [28]. To obtain the p T -differential yield of K 0 S and Λ+Λ, the invariant mass distributions in various p T intervals are parametrised as a sum of two Gaussian distributions with the same mean and a third-order polynomial function which accounts for residual background. The K 0 S and Λ+Λ yields are extracted by integration of the Gaussian distributions and are not corrected for feed-down from higher mass baryons (e.g. Ξ ± , Ω ± ), but these have a negligible effect on v 2 [26]. The Ξ − +Ξ + and Ω − +Ω + are reconstructed through the cascade topology of the following weak decays: and Ω − → Λ + K − (Ω + → Λ + K + ) with branching ratios of 99.9% and 67.8% [49], respectively, with a subsequent Λ (Λ) decay. Candidates are found by applying topological and kinematic criteria first to select the V 0 with an invariant mass between 1.08 and 1.16 GeV/c 2 and then to match it with one of the remaining secondary tracks. They are selected by requiring the DCA between the V 0 and the track to be less than 0.3 cm, the cosine of the pointing angle to be at least 0.999 and 0.998 for the cascade and V 0 , respectively, the DCA between the V 0 and primary vertex to be larger than 0.05 cm, the minimum DCA of V 0 daughter tracks to the primary vertex to be 0.1 cm, the maximum DCA of V 0 daughter tracks to be 1.0 cm, and the minimum DCA of the daughter track to the primary vertex to be 0.03 cm. Only Ξ − +Ξ + and Ω − +Ω + candidates produced at a radial distance between 0.9 and 100 cm from the beam line with the same radial distance reported in Section 3.4 for V 0 are accepted. Each of the three daughter tracks is also required to have p T > 0.15 GeV/c within |η| < 0.8 and to pass the TPC track quality criteria detailed above for charged tracks. In addition, the daughter tracks are checked for compatibility with the pion, kaon, or proton hypotheses by selecting particles with |nσ TPC | < 3 for each species. The p T -differential yield of Ξ − +Ξ + and Ω − +Ω + is obtained by fitting the invariant mass distributions with a sum of two Gaussian distributions with the same mean and a third-order polynomial function that describe the signal and the background, respectively.

Flow observables
A common approach to study the event-by-event flow fluctuations for a given flow coefficient is by using the two-and multiparticle cumulants [51,52], which have different sensitivities to effects stemming from non-flow and flow fluctuations, where δ n denotes the two-particle non-flow effects.
Assuming that flow fluctuation σ v n is relatively small compared to v n , which was found to be true for non-central heavy-ion collisions at the LHC [53][54][55], and also assuming that non-flow effects can be experimentally removed (or largely suppressed, e.g. by using appropriate η gaps) in two-particle correlation measurements, the v n can be written as [56] v where ⟨v n ⟩ and σ v n are the mean and fluctuations of the anisotropic flow coefficient, respectively. These two quantities correspond to the first and second moments of the event-by-event v n distribution. The observable ⟨v n ⟩ is expected to be free from flow fluctuations and to only reflect the true elliptic flow from the flow symmetry plane.
Both ⟨v n ⟩ and σ v n can be calculated using the measured v n {2} and v n {4} as Furthermore, the relative flow fluctuations F(v n ) are defined as

Flow extraction methods
The measurement of the p T -differential v n coefficients of identified hadrons is performed using two-and four-particle cumulant method [51], according to Here c n {m} and d n {m} are the reference and differential m-particle cumulants, respectively, which can be obtained from m-particle correlations. For the specific case of two and four particles, they are given by where ⟨⟨m⟩⟩ and ⟨⟨m ′ ⟩⟩ are the event-averaged reference and differential m-particle correlations, respectively. In order to suppress two-particle non-flow correlations, a pseudorapidity gap of |∆η| > 0.8 between the two particles (subevents) is applied when computing the correlations of Eqs. 12 and 13. The relevant flow coefficients will be denoted as v 2 {2, |∆η| > 0.8} later in the text.
The multiparticle correlation technique with the generic framework implementation [38] allows for correcting for detector inefficiencies and non-uniformities in the azimuthal particle distribution using weights. Using this method, one can measure the p T -differential flow with two-and four-particle cumulants of inclusive charged hadrons, π ± , K ± , and p+p for each centrality percentile. For K 0 S , Λ+Λ, φ , Ξ − +Ξ + , and Ω − +Ω + , the identification on a particle-by-particle basis is not possible. Besides a signal component (true particles), the candidates contain a non-negligible combinatorial background. Considering the fact that the corresponding p T -differential multiparticle cumulants of both signal and background might not be easily decomposed into individual contributions, the invariant mass method [57] was developed exploiting the additivity of correlation. This method has been used in previous anisotropic flow measurements for reconstructed particles with the two-particle correlation method [26, 28, 33] and is also used in this analysis. The two-and four-particle correlations are measured as a function of both invariant mass and candidate p T . The relation between the signal and background components is given for each p T interval by where the total m-particle correlation ⟨⟨m ′ ⟩⟩ total n can be regarded as the sum of the m-particle correlations of signal particles ⟨⟨m ′ ⟩⟩ signal n and the correlations of the background ⟨⟨m ′ ⟩⟩ bg n . Here the signal function is weighted by its corresponding fraction f signal defined as where N signal (m inv ) and N bg (m inv ) are signal and background yields, respectively. Correspondingly, the weight of the background function is determined with f bg = 1− f sig . Both N signal (m inv ) and N bg (m inv ) are obtained from the invariant mass distribution of K 0 S , Λ+Λ, φ , Ξ − +Ξ + , and Ω − +Ω + for each p T interval and centrality class following the procedures outlined in Sections 3.3, 3.4, and 3.5. After measuring ⟨⟨m ′ ⟩⟩ total n (m inv ) via multiparticle correlations, the ⟨⟨m ′ ⟩⟩ signal n can be determined for a given centrality class and p T interval using a simultaneous fit to the ⟨⟨m ′ ⟩⟩ total n (m inv ) and N total (m inv ) distributions, where ⟨⟨m ′ ⟩⟩ bg n (m inv ) is parametrised as a first-order polynomial function. An example of such a procedure with distributions of m inv , ⟨⟨2 ′ ⟩⟩, and ⟨⟨4 ′ ⟩⟩ together with fit functions is shown in Fig. 1

Systematic uncertainties
The systematic uncertainties were evaluated by varying the selection criteria for each particle species, in every p T range and centrality interval, such as track quality criteria for π ± , K ± , and p+p or topological reconstruction requirements for φ , K 0 S , Λ+Λ, Ξ − +Ξ + , and Ω − +Ω + . Only statistically significant differences between the nominal data points and the systematic variations, where significance is evaluated based on the recommendations in Ref. [58], were assigned as systematic uncertainties. The uncertainties from the independent sources were added in quadrature to obtain the final systematic uncertainties on the measurements. For each particle species, a p T -dependent systematic uncertainty is reported for p T < 3 GeV/c, while a p T -independent average uncertainty is assigned at higher transverse momenta to suppress statistical fluctuations. Tables 1, 2, and 3 summarise the minimum and maximum values of the relative systematic uncertainties per particle species for all p T and centrality ranges.
For the event selection criteria, the primary vertex position along the beam line was varied from 10 cm to 8 cm, and the centrality determination was changed from energy deposition in the V0 scintillator to the number of hits in the first layer of the ITS. Additionally, for π ± , K ± , and p+p, the event sample was separated based on the polarity of the ALICE solenoid, and the two sub-samples were studied independently.
Systematic uncertainties arising from the selection criteria imposed at the track level were investigated by requiring 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). In addition, systematics uncertainties in the measurements were investigated by increasing the minimum number of TPC space points from 70 to 90, and by varying the DCA from the strict p T -dependent selection to 0.15 cm in the transverse plane and from 2 cm to 0.2 cm in the longitudinal direction, in order to estimate the impact of secondary particles. These variations are referred as "tracking mode" in Tables 1 and 3. Furthermore, the minimal probability threshold in the Bayesian particle identification was increased from 0.95 to 0.98 for π ± and from 0.85 to 0.9 for K ± and p+p.
0-2% 0-5% 0-5% 0-1% 0-1% 0-2% Bayesian particle identification 0-5% 0-5% 0-4% 0-5% 0-4% 0-4% For K 0 S and Λ+Λ, the topological requirements on the V 0 s themselves were varied by changing the maximum DCA of the V 0 daughter tracks to the secondary vertex from 0.5 cm to 0.3 cm and the minimum DCA of the V 0 daughter tracks to the primary vertex from 0.1 cm to 0.3 cm. In addition, the minimum radial distance to the primary vertex at which the V 0 can be produced was changed from 5 cm to 1 cm and 10 cm. The selection criteria imposed on the daughter tracks were varied by increasing the minimum number of TPC space points from 70 to 90, requiring the ratio between the number of space points and the number of crossed rows in the TPC to be larger than 0.9 or 1.0 instead of 0.8 (denoted as "track quality" in Table 2), and requiring a minimum p T of 0.2 GeV/c. Finally, the strategy for reconstructing V 0 was changed from online, where V 0 s were reconstructed during the track fitting procedure, to offline, which took place after all the tracks have been reconstructed.
The Ξ − +Ξ + and Ω − +Ω + finding criteria were varied by changing the maximum DCA between the V 0 and bachelor track from 0.3 cm to 0.25 cm, increasing the minimum DCA between the V 0 and primary vertex from 0.05 cm to 0.06 cm, decreasing the minimum DCA of V 0 bachelor tracks to the primary vertex from 0.1 cm to 0.08 cm. In addition, the criteria were changed by requiring the maximum DCA of V 0 bachelor tracks to be 0.8 cm instead of 1.0 cm, increasing the minimum DCA of the bachelor track to the primary vertex from 0.03 cm to 0.035 cm, and decreasing the minimum value of the cosine of the pointing angle for the cascade from 0.999 to 0.995. For the V 0 , the invariant mass range was changed from (1.08-1.16) GeV/c 2 to (1.10-1.14 ) GeV/c 2 , while the minimum radial distance was varied from 5 cm to 1 cm and 10 cm. A minimum value of 0.995 instead of 0.998 was used for the cosine of the pointing angle. For each of the three daughter tracks, the PID criterion was varied from |nσ TPC | < 3 to |nσ TPC | < 2 and the minimum p T was increased to 0.2 GeV/c.
An additional contribution from fitting parameter variations was studied for all the reconstructed particles, following the same approach used in previous works [26,33]. The resulting systematic uncertain- The p T -differential v 2 measured with two-particle correlations with a pseudorapidity gap of |∆η| > 0.8 for different particle species and centralities in Pb-Pb collisions at √ s NN = 5.02 TeV. The vertical error bars and the filled boxes represent statistical and systematic uncertainties, respectively. ties, negligible for K 0 S , Λ+Λ, and φ mesons, but significant for Ξ − +Ξ + and Ω − +Ω + , are summarized in Table 3.

Results
In this section, the results for the p T -differential v 2 and its relative fluctuations measured in various collision centrality intervals of Pb-Pb collisions at √ s NN = 5.02 TeV are presented. The v 2 measured with two-and four-particle cumulants and the corresponding results for flow fluctuations for different particle species are reported in Sec. 5.1 and Sec. 5.2, respectively. In Section 5.3, the experimental data are compared with a state-of-the-art hydrodynamic model calculations, namely, the coupled linear  Boltzmann transport (CoLBT) [59], that applies T R ENTo initial state model and incorporates the bulk expansion of the medium with a specific shear viscosity η/s = 0.10 and interactions of energetic partons with it, as well as a coalescence mechanism for particle production. Note that the same data will be shown in different representations to highlight the various physics implications of the measurements. The data points will be drawn together with their statistical and systematic uncertainties represented by the error bars and shaded boxes around each marker, respectively. Figure 2 presents the p T -differential v 2 {2, |∆η| > 0.8} measurements for unidentified charged hadrons (h ± ) as well as for π ± , K ± , p+p, K 0 S , Λ+Λ, φ , Ξ − +Ξ + , and Ω − +Ω + from Pb-Pb collisions at √ s NN = 5.02 TeV. The five panels show different centrality intervals, with the most central and the most peripheral, i.e. 10-20% and 50-60%, drawn in the top left and bottom right, respectively. This analysis profits from the data samples collected by ALICE in 2015 and 2018 which allow extending the previous results from two-particle correlations [25-27, 33] to higher-order cumulants. Figure 3 presents the first p T -differential v 2 measurements using four-particle cumulants (i.e. v 2 {4}) for the same particle species as reported in Fig. 2 from Pb-Pb collisions at √ s NN = 5.02 TeV. In both cases, similar features of the p T -differential measurements as reported and discussed in detail in Refs. [26-28, 33] are confirmed. The progressive increase of v 2 with the centrality of the collision for a given p T interval illustrates the final-state anisotropy that originates from the initial-state ellipsoidal geometry in non-central collisions, quantified by the spatial eccentricity ε 2 . Furthermore, both the effects known in the literature as mass ordering and the meson-baryon particle type grouping are present in these new measurements. The former originates from the radial flow of the system, while the latter is explained in a dynamical picture where flow develops at the partonic level followed by quark coalescence into hadrons [32].

Mass ordering and scaling properties
The meson-baryon grouping is generally attributed to hadron production via coalescence in the intermediate p T region [32], where the direct contribution from hydrodynamic expansion may no longer be dominant and the path-length dependence of energy loss might not play a significant role yet [60]. This      Figure 6: The p T dependence of the standard deviation of v 2 (σ v 2 ) for different particle species and centralities in Pb-Pb collisions at √ s NN = 5.02 TeV. The vertical error bars and the filled boxes represent statistical and systematic uncertainties, respectively.

Results on flow fluctuations
The measurements of v 2 with two-and four-particle cumulants provide the first opportunity to investigate the first moments of the v 2 distribution for different particle species. Figure 5 presents the mean value of v 2 , denoted as ⟨v 2 ⟩, as a function of p T for the same combination of hadrons and centrality intervals as in the previous figures. Assuming that the non-flow contribution in v 2 2 {2, |∆η| > 0.8} is negligible [21,36], this mean value is calculated according to Eq. 7 by replacing v 2 2 {2} with v 2 2 {2, |∆η| > 0.8} measurement. Figure 6 presents the transverse momentum dependence of the second moment of the v 2 distribution, i.e. the standard deviation σ v 2 , measured for the first time for different particle species. As in the previous case, assuming negligible non-flow contribution in v 2 2 {2, |∆η| > 0.8}, σ v 2 is approximated according to Eq. 8. The data points of both ⟨v 2 ⟩ and σ v 2 show, as expected, the same qualitative features as in the previous cases of Figs. 2 and 3: namely the mass ordering developing at low values of p T and the particle type grouping that is evident at higher p T .
Combining ⟨v 2 ⟩ and σ v 2 , one can quantify the relative v 2 fluctuations (F(v 2 )) according to Eq. 9. This quantity is displayed in Fig. 7 as a function of p T and centrality intervals for the various particle species presented in this article. It can be seen that for central events, there is no significant p T or particle species dependence. However, for more peripheral collisions, and in particular starting from the interval 30-40% and above, the data points exhibit a non-monotonic transverse momentum dependence, with a minimum in F(v 2 ) that lies at higher values of p T for baryons than for mesons. In addition, the F(v 2 ) for baryons in 1 < p T < 3 GeV/c is systematically lower than for mesons. Interestingly, the momentum region where this apparent particle type grouping develops for F(v 2 ) does not coincide with the region where a similar grouping is reported for measurements of v 2 . This could point to a different origin for these two observations. For p T > 3 GeV/c, all data points converge into a universal band within the uncertainties. The origin of this characteristic behaviour of F(v 2 ) is studied using hydrodynamical models in the following section.
Finally, to further study the nature of flow fluctuations, Fig. 8 presents the p T -dependence of the ratio v 2 {4}/v 2 {2, |∆η| > 0.8}. This ratio is expected to be sensitive to the fluctuations within the picture of initial state models. Within these models, the spatial eccentricity ε 2 fluctuates from event to event. These fluctuations are transferred through the low viscosity QGP to the final state and are imprinted in how v 2 fluctuates. Since v 2 ∝ ε 2 , the ratio v 2 {4}/v 2 {2, |∆η| > 0.8} is expected to reflect the ratio between ε 2 {4} and ε 2 {2}, which have positive and negative contributions from the initial eccentricity fluctuations, and thus can provide strong constraints on initial state models. It can be seen that for central collisions, this ratio does not exhibit any significant p T or particle species dependence. Starting from the 20-30% centrality interval, however, a decrease in the ratio can be seen between 1 and 5 GeV/c. It becomes progressively more pronounced for more peripheral events. In addition, starting from the 30-40% centrality interval, and similar to the picture that develops for F(v 2 ), the data points indicate a particle type grouping, with the values for baryons being systematically larger than the ones for mesons. This apparent dependence on particle species highlights that final state effects play a significant role in these observables.

Comparison with models
The comparison of results from anisotropic flow studies with hydrodynamic calculations has been instrumental in constraining some of the basic transport coefficients of the QGP. However, such comparisons were limited until now to the low p T region, i.e. in ranges where the mass ordering discussed in the previous section is prominent. One of the first attempts to provide a unified physics picture throughout the entire transverse momentum range for different particle species was presented recently in Ref. [60]. In this article, the authors used the CoLBT hydrodynamic model [59] which allows for the simultaneous description of the evolution of parton showers and the bulk medium. The latter is prescribed by a (3+1)-D viscous hydrodynamic model that is initialized at τ 0 = 0.6 fm/c and uses a value of specific shear viscosity η/s = 0.10. The freeze-out temperature is set to T fo = 150 MeV, beyond which a hadronic after-burner describes the interactions between hadrons. The remaining parameters of the model were adjusted to reproduce the measured yields, p T spectra, and integrated v n of unidentified charged hadrons in Pb-Pb collisions. One of the important ingredients which is introduced in this model is the way ) p p (   Figure 9: The p T -differential v 2 {4} for π ± , K ± , and p+p measured in Pb-Pb collisions at √ s NN = 5.02 TeV compared with expectations of the same quantity from the CoLBT hydrodynamic model with quark coalescence [60]. The left and right panels present the comparison for the 10-20% and 40-50% centrality intervals, respectively. The vertical error bars and the filled boxes represent statistical and systematic uncertainties of the data, respectively. The thickness of the model curves reflect the uncertainties of the hydrodynamic calculations. hadrons emerge, with the typical hydrodynamic freeze-out at low p T being complemented by a quark coalescence prescription at intermediate p T and fragmentation at high p T [60]. Figure 9 presents the evolution of v 2 {4} as a function of p T for π ± , K ± , and p+p for two characteristic centrality intervals, 10-20% (central) and 40-50% (peripheral), in the left and right panels, respectively. The measurements are compared with the expectations for the same particle species from the CoLBT hydrodynamic model, represented by the shaded bands. It can be seen that the model describes the p T dependence of v 2 {4} over the entire p T range. In particular, at low values of p T (< 2-3 GeV/c) where the hydrodynamic expansion of the medium plays a dominant role, the model describes both the increase as a function of p T and the mass ordering. The v 2 {4} reaches a peak value at around p T ≈ 3 GeV/c for pions and kaons and at p T ≈ 4 GeV/c for protons, before decreasing at high p T . This can be naturally explained by the interplay between the hydrodynamical expansion, hadron production through quark coalescence, and jet fragmentation [60].
Within the CoLBT model, the hydrodynamic contribution to v 2 is dominant for all particle species up to p T = 4 GeV/c, whereas the jet fragmentation plays an increasingly important role for p T > 6 GeV/c. In the intermediate p T region (4-6 GeV/c), quark coalescence contributes in CoLBT to the development of the value of v 2 , even though in the model this mechanism accounts for less than 25% of the total particle yield. This is because the value of v 2 from coalescence is significantly larger than the v 2 from fragmentation up to 6 GeV/c. The additional mechanism of the coalescence prescription in the model is important to reproduce the experimental results quantitatively and to provide the proper connection between the low and high p T regions. In the former, the mass ordering develops, while in the latter the fragmentation is the dominant particle production mechanism and no significant particle species dependence is observed.
It is known that neither the hydrodynamic expansion nor the fragmentation alone leads to the precise NCQ scaling development. Such contributions in the final v 2 could consequently give a natural explanation for the significant deviation from a universal NCQ scaling observed in Fig. 4. On the other hand, the model calculation with only contributions from hydrodynamic expansion and fragmentation but without the contribution from quark coalescence significantly underestimates v 2 {4} for p T above 3 GeV/c for π ± , K ± , shown in Fig. 10. Nevertheless, the crossing of v 2 of pions, kaons, and protons can develop according to CoLBT with a combination of the hydrodynamic expansion coupled only to jet fragmentation (for details refer to Fig. 4 in Ref. [60]). This might challenge the prevailing idea discussed in the literature (see, e.g., Ref. [61]) that the crossing can be attributed to quark coalescence. It is important to note that the development of the crossing point in the absence of coalescence in CoLBT arises from a particle species-dependent p T value where fragmentation becomes dominant over hydrodynamics. hydrodynamics and jet fragmentation, respectively [60]. The 40-50% centrality interval was chosen as representative for these comparisons. The model without quark coalescence contribution describes qualitatively the features and the p T dependence of the measurements, but significantly underestimates ⟨v 2 ⟩ of pion and kaon for p T above 3 GeV/c. This is very different from what has been observed in Figs. 9 and 11. Despite the sizable uncertainties of CoLBT calculations, the contribution from quark coalescence seems non-negligible for σ v 2 , v 2 {4}/v 2 {2}, and F(v 2 ), when comparing the calculations of hydro+coal+frag (shown in Fig. 11) and hydro+frag (shown in Fig. 12).

Summary
In summary, the first measurement of p T -differential elliptic flow using two-and four-particle cumulants for π ± , K ± , p+p, K 0 S , Λ+Λ, φ , Ξ − +Ξ + , and Ω − +Ω + in Pb-Pb collisions at √ s NN = 5.02 TeV is presented. The mean elliptic flow, elliptic flow fluctuations, and relative elliptic flow fluctuations are obtained for various particle species. Differences in the value of relative flow fluctuations for different particle species are observed, suggesting that final state hadronic interactions further modify the flow fluctuations. A distinct mass ordering is found in the 10-60% centrality interval for p T < 3 GeV/c, which arises from the interplay between the elliptic and radial flow. In the intermediate p T range, the magnitude of v 2 {4}, ⟨v 2 ⟩, and σ v 2 for baryons is larger than that for mesons by about 50%. In addition, particles show an approximate constituent quark scaling. This scaling is tested for v 2 {4}, which is expected to measure flow with little (or no) non-flow contamination. NCQ scaling describes the data no better than ±20%, an accuracy similar to what was reported for the v 2 using two-particle correlations. Furthermore, the relative flow fluctuation F(v 2 ) for the identified hadrons shows an apparent splitting between baryons and mesons for centrality above 30%, which suggests a significant role for final-state interactions in developing this observable. Last but not least, CoLBT hydrodynamic calculations with the implementation of quark coalescence describe the measurements over a large p T range, which confirms the relevance of the quark coalescence hadronization mechanism in the particle production in Pb-Pb collisions at the LHC.