Non-linear flow modes of identified particles in Pb-Pb collisions at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV

The pT-differential non-linear flow modes, v4,22, v5,32, v6,33 and v6,222 for π±, K±, KS0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathrm{K}}_{\mathrm{S}}^0 $$\end{document} , p + p¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{p}} $$\end{document}, Λ + Λ¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\Lambda} $$\end{document} and ϕ-meson have been measured for the first time at sNN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s_{\mathrm{NN}}} $$\end{document} = 5.02 TeV in Pb-Pb collisions with the ALICE detector at the Large Hadron Collider. The results were obtained with a multi-particle technique, correlating the identified hadrons with reference charged particles from a different pseudorapidity region. These non-linear observables probe the contribution from the second and third order initial spatial anisotropy coefficients to higher flow harmonics. All the characteristic features observed in previous pT-differential anisotropic flow measurements for various particle species are also present in the non-linear flow modes, i.e. increase of magnitude with increasing centrality percentile, mass ordering at low pT and particle type grouping in the intermediate pT range. Hydrodynamical calculations (iEBE-VISHNU) that use different initial conditions and values of shear and bulk viscosity to entropy density ratios are confronted with the data at low transverse momenta. These calculations exhibit a better agreement with the anisotropic flow coefficients than the non-linear flow modes. These observations indicate that non-linear flow modes can provide additional discriminatory power in the study of initial conditions as well as new stringent constraints to hydrodynamical calculations.


Introduction
Lattice quantum chromodynamics (QCD) calculations [1,2] suggest that at extremely high temperature and energy density a state of matter is produced in which quarks and gluons are no longer confined into hadrons. This state of matter is called the quark-gluon plasma (QGP) [3][4][5]. The main goal of heavy-ion collision experiments is to study the properties of the QGP, such as the speed of sound, the equation of state and its shear and bulk viscosities.
One of the observables sensitive to these properties is the azimuthal angular distribution of particles emitted in the plane perpendicular to the beam axis. In a heavy-ion collision, the overlap region of the colliding nuclei exhibits an irregular shape [6][7][8][9][10][11][12]. This spatial irregularity is a superposition of the geometry, i.e. centrality [13] of the collision reflected in the value of the impact parameter, and the initial energy density in the transverse plane which fluctuates from event to event. Through interactions between partons -1 -

JHEP06(2020)147
and at later stages between the produced particles, this spatial irregularity is transferred into an anisotropy in momentum space. The latter is usually decomposed into a Fourier expansion of the azimuthal particle distribution [14] according to (1.1) where N , p T , η and ϕ are the particle yield, transverse momentum, pseudorapidity and azimuthal angle of particles, respectively, and Ψ n is the azimuthal angle of the n th -order symmetry plane [7][8][9][10]12]. The coefficient v n is the magnitude of the n th -order flow vector coefficient V n , defined as V n = v n e inΨn , and can be calculated according to where the angle brackets denote an average over all particles in all events. Since the symmetry planes are not accessible experimentally, the flow coefficients are estimated solely from the azimuthal angles of the particles emitted in the transverse plane. Measurements of different anisotropic flow coefficients at both the Relativistic Heavy Ion Collider (RHIC) [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] and the Large Hadron Collider (LHC) [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46] not only confirmed the production of a strongly coupled quark-gluon plasma (sQGP) but also contributed in constraining the value of the ratio between shear viscosity and entropy density (η/s) which is very close to the lower limit of 1/4π conjectured by AdS/CFT [47]. In addition, the comparison between experimental data [41] and viscous hydrodynamical calculations [48] showed that higher order flow coefficients and more importantly their transverse momentum dependence are more sensitive probes than lower order coefficients, i.e. v 2 and v 3 , to the initial spatial irregularity and its fluctuations [10]. This initial state spatial irregularity is usually quantified with the standard (momentdefined) anisotropy coefficients, n . In the Monte Carlo Glauber model, n and its corresponding initial symmetry plane, Φ n can be calculated from the transverse positions of the nucleons participating in a collision according to [9,49] n e inΦn = r n e inϕ r n (for n > 1), (1.3) where the brackets denote an average over the transverse position of all participating nucleons that have an azimuthal angle ϕ and a polar distance from the centre r. Model calculations show that v 2 and to a large extent, v 3 are for a wide range of impact parameters linearly proportional to their corresponding initial spatial anisotropy coefficients, 2 and 3 , respectively [9], while for larger values of n, v n scales with n , a cumulant-based definition of initial anisotropic coefficients. As an example, the fourth order spatial anisotropy is given by [50,51] 4 e i4Φ 4 = 4 e i4Φ 4 + 3 r 2 2 r 4 2 2 e i4Φ 2 , (1.4) where the second term in the right hand side of eq. (1.4) reveals a non-linear dependence of 4 on the lower order 2 . This further supports the earlier ideas that the higher order -2 -JHEP06(2020)147 flow vector coefficients, V n (n > 3) obtain contributions not only from the linear response of the system to n , but also a non-linear response proportional to the product of lower order initial spatial anisotropies [52,53]. In particular, for a single event, V n with n = 4, 5, 6 can be decomposed to the linear (V L n ) and non-linear (V NL n ) modes according to where χ n,mk , known as non-linear flow mode coefficients, quantify the contributions of the non-linear modes to the total V n [53,54]. For simplicity, the magnitude of the total V n will be referred to as anisotropic flow coefficient (v n ) in the rest of this article. The magnitude of the p T -differential non-linear modes for higher order flow coefficients, v NL n , can be written as: where brackets denote an average over all events. The approximation is valid assuming a weak correlation between the lower (n = 2, 3) and higher (n > 3) order flow coefficients [52,55]. Various measurements of the p T -differential anisotropic flow, v n (p T ), of charged particles [33,38,43,45,46,56] provided a testing ground for model calculations that attempt to describe the dynamical evolution of the system created in heavy-ion collisions. Early predictions showed that the p T -differential anisotropic flow for different particle species can reveal more information about the equation of state, the role of the highly dissipative hadronic rescattering phase as well as probing particle production mechanisms [57,58]. In order to test these predictions, v n (p T ) coefficients were measured for different particle species at RHIC [15][16][17][18] and at the LHC [39,40,42,44]. These measurements reveal a characteristic mass dependence of v n (p T ) in the low transverse momentum region (p T < 3 GeV/c), a result of an interplay between radial and anisotropic flow, and mass dependent thermal velocities [57,58]. In the intermediate p T region (3 p T 8 GeV/c) the measurements indicate a particle type grouping where baryons have a larger v n than the one of mesons. This feature was explained in a dynamical model where flow develops at the partonic level followed by quark coalescence into hadrons [59,60]. In this picture the invariant spectrum of produced particles is proportional to the product of the spectra of their constituents and, in turn, the flow coefficient of produced particles is the sum of the v n values of their -3 -

JHEP06(2020)147
and a length of 500 cm and it is positioned around the ITS. It provides full azimuthal coverage in the pseudorapidity range |η| < 0.9.
Charged particles were identified using the information from the TPC and the TOF detectors [63]. The TPC allows for a simultaneous measurement of the momentum of a particle and its specific energy loss ( dE/dx ) in the gas. The detector provides a separation more than two standard deviations (2σ) for different hadron species at p T < 0.7 GeV/c and the possibility to identify particles on a statistical basis in the relativistic rise region of dE/dx (i.e. 2 < p T < 20 GeV/c) [64]. The dE/dx resolution for the 5% most central Pb-Pb collisions is 6.5% and improves for more peripheral collisions [64]. The TOF detector is situated at a radial distance of 3.7 m from the beam axis, around the TPC and provides a 3σ separation between π-K and K-p up to p T = 2.5 GeV/c and p T = 4 GeV/c, respectively [64]. This is done by measuring the flight time of particles from the collision point with a resolution of about 80 ps. The start time for the TOF measurement is provided by the T0 detectors, two arrays of Cherenkov counters positioned at opposite sides of the interaction points covering 4.6 < η < 4.9 (T0A) and −3.3 < η < −3.0 (T0C). The start time is also determined using a combinatorial algorithm that compares the timestamps of particle hits measured by the TOF to the expected times of the tracks, assuming a common event time t ev . Both methods of estimating the start time are fully efficient for the 80% most central Pb-Pb collisions [64].
A set of forward detectors, the V0 scintillator arrays [67], were used in the trigger logic and for the determination of the collision centrality. The V0 consists of two detectors, the V0A and the V0C, positioned on each side of the interaction point, covering the pseudorapidity intervals of 2.8 < η < 5.1 and −3.7 < η < −1.7, respectively.
For more details on the ALICE apparatus and the performance of the detectors, see refs. [63,64].
3 Event sample, track selection and particle identification

Trigger selection and data sample
The analysis is performed on minimum bias Pb-Pb collision data at √ s NN = 5.02 TeV collected by the ALICE detector in 2015. These events were triggered by the coincidence between signals from both V0A and V0C detectors. An offline event selection, exploiting the signal arrival time in V0A and V0C, measured with a 1 ns resolution, was used to discriminate beam induced-background (e.g. beam-gas events) from collision events. This led to a reduction of background events in the analysed samples to a negligible fraction (< 0.1%) [64]. Events with multiple reconstructed vertices were rejected by comparing multiplicity estimates from the V0 detector to those from the tracking detectors at midrapidity, exploiting the difference in readout times between the systems. The fraction of pileup events left after applying these dedicated pileup removal criteria is negligible. All events selected for the analysis had a reconstructed primary vertex position along the beam axis (z vtx ) within 10 cm from the nominal interaction point. After all the selection criteria, a filtered data sample of approximately 40 million Pb-Pb events in the 0-60% centrality interval was analysed to produce the results presented in this article.

JHEP06(2020)147
Events were classified according to fractions of the inelastic hadronic cross section. The 0-5% interval represents the most central interactions (i.e. smallest impact parameter) and is referred to as most central collisions. On the other hand, the 50-60% centrality interval corresponds to the most peripheral (i.e. largest impact parameter) collisions in the analysed sample. The centrality of the collision was estimated using the signal amplitude measured in the V0 detectors which is related to the number of particles crossing their sensitive areas. Details about the centrality determination can be found in [68].
3.2 Selection of primary π ± , K ± and p + p In this analysis, tracks are reconstructed using the information from the TPC and the ITS detectors. The tracking algorithm, based on the Kalman filter [69,70], starts from a collection of space points (referred to as clusters) inside the TPC and provides the quality of the fit by calculating its χ 2 value. Each space point is reconstructed at one of the TPC pad rows [63], where the deposited ionisation energy is also measured. The specific ionisation energy loss dE/dx is estimated using a truncated mean, excluding the 40% highest-charge clusters associated to the track. The obtained dE/dx has a resolution, which we later refer to as σ TPC . The tracks are propagated to the outer layer of the ITS, and the tracking algorithm attempts to identify space points in each of the consecutive layers, reaching the innermost ones (i.e. SPD). The track parameters are then updated using the combined information from both the TPC and the ITS detectors.
Primary charged pions, kaons and (anti-)protons were required to have at least 70 reconstructed space points out of the maximum of 159 in the TPC. The average distance between space point and the track fit per TPC space point per degree of freedom (see [64] for details) was required to be below 2. These selections reduce the contribution from short tracks, which are unlikely to originate from the primary vertex. To reduce the contamination by secondary tracks from weak decays or from the interaction with the material, only particles within a maximum distance of closest approach (DCA) between the tracks and the primary vertex in both the transverse plane (DCA xy < 0.0105 + 0.0350(p T c/GeV) −1.1 cm) and the longitudinal direction (DCA z < 2 cm) were analysed. Moreover, the tracks were required to have at least two associated ITS clusters in addition to having a hit in either of the two SPD layers. This selection leads to an efficiency of about 80% for primary tracks at p T ∼ 0.6 GeV/c and a contamination from secondaries of about 5% at p T = 1 GeV/c [71]. These values depend on particle species and transverse momentum [71].
The particle identification (PID) for pions (π ± ), kaons (K ± ) and protons (p + p) used in this analysis relies on the two-dimensional correlation between the number of standard deviations in units of the resolution from the expected signals of the TPC and the TOF detectors similar to what was reported in [39,40,42]. In this approach particles were selected by requiring their standard deviations from the dE/dx and t TOF values to be less than a p T -dependent value, maintaining a minimum purity of 90% for π ± and 75% for K ± and 80% for p + p. In order to further reduce the contamination from other species, the standard deviation of a given track was required to be the minimum among other candidate species.

JHEP06(2020)147
In addition, for the evaluation of systematic effects (see section 5) the minimum purity was varied to more strict values, a condition that becomes essential with increasing transverse momentum where the relevant detector response for different particle species starts to overlap. The results for all three particle species were extrapolated to 100% purity and the uncertainty from the extrapolation was also considered in the estimation of the total systematic uncertainty.
3.3 Reconstruction of K 0 S , Λ + Λ and φ meson In this analysis, the K 0 S and Λ + Λ are reconstructed via the following fully hadronic decay channels: K 0 S → π + + π − and Λ(Λ) → p(p) + π − (π + ) with branching ratios of 69.2% and 63.9% [72], respectively. The reconstruction is performed by identifying the candidates of secondary vertices, denoted as V 0 s, from which two oppositely-charged decay products originate. Such candidates are obtained during data processing by looking for a characteristic V-shaped decay topology among pairs of reconstructed tracks.
The daughter tracks were reconstructed within |η| < 0.8, while the criteria on the number of TPC space points, 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 the daughter tracks to the primary vertex is 0.1 cm. Furthermore, the maximum DCA of the daughter tracks is 0.5 cm to ensure that they are products of the same decay. To suppress the combinatorial background, the PID is applied for the daughter particles in the whole p T region by requiring the particle to be within 3σ TPC for a given species hypothesis.
To reject combinatorial background, the cosine of the pointing angle, θ p , was required to be larger than 0.998. This angle is defined as the angle between the momentum vector of the V 0 candidate assessed at its decay vertex and the line connecting the V 0 decay vertex to the primary vertex and has to be close to 1 as a result of momentum conservation. In addition, only the candidates reconstructed between 5 and 100 cm from the nominal primary vertex in radial direction were accepted. The lower value was 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, where the occupancy is the largest. To assess the systematic uncertainty related to the contamination from Λ + Λ and electron-positron pairs coming from γ-conversions to the K 0 S sample, a selection in the Armenteros-Podolanski variables [73] was applied for the K 0 S candidates, rejecting the ones with q ≤ 0.2|α|. 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 track momentum onto the momentum of the V 0 .
The reconstruction of φ meson candidates is done via the hadronic decay channel: φ → K + + K − with a branching ratio of 48.9% [72]. The φ meson candidates were reconstructed from the charged tracks passing all criteria for charged kaons. These kaon daughters were identified utilising the Bayesian PID approach [74] with a minimum probability threshold of 85% using the TPC and TOF detectors. Additionally, to reduce combinatorial background, a track was identified as a kaon if it had the highest probability among all considered species (e ± , µ ± , π ± , K ± and p + p). The vector sum of all possible pairs of charged kaons are -7 -JHEP06(2020)147 called φ candidates. The invariant mass distribution (M K + K − inv ) of φ candidates was then obtained in various p T intervals by subtracting a combinatorial background yield from the candidate yield. This combinatorial background yield was estimated from like-sign kaon pairs (unphysical φ state with total charge of ±2) normalised to the candidate yield.

Analysis method
In this article the p T -differential non-linear flow modes are calculated based on eqs. (1.6)-(1.9). Each event is divided into two subevents "A" and "B", covering the ranges −0.8 < η < 0.0 and 0.0 < η < 0.8, respectively. Thus v n,mk (p T ) is a weighted average of v A n,mk (p T ) and v B n,mk (p T ). The measured v n,mk (p T ) coefficients are calculated using d n,mk (p T ) and c mk,mk multi-particle correlators given by These correlators were obtained using the Generic Framework with sub-event method originally used in [54,75,76], which allows precise non-uniform acceptance and efficiency corrections. In this analysis, d n,mk (p T ) is measured by correlating the azimuthal angle of the particle of interest (ϕ 1 (p T )) from subevent "A"("B") with that of reference particles 1 from subevent "B"("A") and c mk,mk by selecting half of the reference particles from subevent "A" and the other half from "B". Thus, eqs. , , where denotes an average over all particles and events. This multi-particle correlation technique by construction removes a significant part of non-flow correlations. In order to further reduce residual non-flow contributions, a pseudorapidity gap was applied between the two pseudorapidity regions (|∆η| > 0.4). In addition, particles with like-sign charges were correlated. These two variations do not significantly affect the results but any variation was included in the final systematics in table 1.
For charged hadrons, i.e. π ± , K ± and p + p, the d n,mk correlators are calculated on a track-by-track basis as a function of p T for each centrality percentile. For particle species JHEP06(2020)147 reconstructed on a statistical basis from the decay products, i.e. K 0 S , Λ+Λ and φ meson, the selected sample contains both signal and the background. Therefore, the d n,mk correlators are measured as a function of invariant mass (M inv ) and p T for each centrality percentile. The d n,mk vs. M inv method is based on the additivity of correlations and is a weighted sum of the d sig n,mk and d bkg n,mk according to where N sig and N bkg are signal and background yields obtained for each p T interval and centrality percentile from fits to the K 0 S , Λ + Λ and φ meson invariant mass distributions. To obtain the p T -differential yield of K 0 S and Λ + Λ, the invariant mass distributions at various p T intervals were parametrised as a sum of two Gaussian distributions and a thirdorder polynomial function. The latter was introduced to account for residual contamination (background yield) that is present in the K 0 S and Λ + Λ signals after the topological and daughter track selections. The K 0 S and Λ + Λ yields were extracted by integration of the Gaussian distribution. The obtained yields were not corrected for feed-down from higher mass baryons (Ξ ± ,Ω ± ) as earlier studies have shown that these have a negligible effect on the measured v n [39]. Similarly, to obtain the p T -differential yield of φ-mesons, the invariant mass distributions of the candidate yield was parametrized as a sum of a Breit-Wigner distribution and a third-order polynomial function, the latter introduced to account for residual contamination.
To extract d sig n,mk in a given p T range, d total n,mk (M inv ) was fitted together with the fit values from the invariant mass distribution and parametrising d bkg n,mk (M inv ) with a first order polynomial function.

Systematic uncertainties
The systematic uncertainties were estimated by varying the selection criteria for all particle species as well as the topological reconstruction requirements for K 0 S , Λ + Λ and φ. The contributions from different sources were extracted from the relative ratio of the p Tdifferential v n,mk between the default selection criteria described in section 3 and their variations summarised in this section. Sources with a statistically significant contribution (where significance is evaluated as recommended in [77]) were added in quadrature to form the final value of the systematic uncertainties on the non-linear flow modes. An overview of the magnitude of the relative systematic uncertainties per particle species is given in table 1 for π ± , K ± and p + p and table 2 for K 0 S , Λ + Λ and the φ-meson. The systematic uncertainties are grouped into five categories, i.e. event selection, tracking, particle identification, topological cuts and non-flow contribution and are described below.
The effects of event selection criteria on the measurements were studied by: (i) varying the primary vertex position along the beam axis (z vtx ) from a nominal ±10 cm to ±8 cm and ±6 cm; (ii) changing the centrality estimator from the signal amplitudes in the V0 scintillator detectors to the number of clusters in the first or second layer of the SPD, (iii) analysing events recorded for different magnetic field polarities independently; (iv) not rejecting all events with tracks caused by pileup. Systematic uncertainties induced by the selection criteria imposed at the track level were investigated by: (i) changing the tracking from global mode, where combined track information from both TPC and ITS detectors are used, to what is referred to as hybrid mode. In the latter mode, track parameters from the TPC are used if the algorithm is unable to match the track reconstructed in the TPC with associated ITS clusters; (ii) increasing the number of TPC space points from 60 up to 90 and (iii) decreasing the value of the χ 2 per TPC space point per degree of freedom from 4 to 3; (iv) varying the selection criteria on both the transverse and longitudinal components of the DCA to estimate the impact of secondary particles from a strict p T -dependent cut to 0.15 cm and 2 cm to 0.2 cm, respectively.
Systematic uncertainties associated with the particle identification procedure were studied by varying the PID method from a p T -dependent one described in section 3.2 to an even stricter version where the purity increases to higher than 95% (π ± ), 80% (K ± ) and 80% (p + p) across the entire p T range of study. The second approach relied on the Bayesian method with a probability of at least 80% which gives an increase in purity to at least 97% (π ± ), 87% (K ± ) and 90% (p + p) across the entire p T range of study. To further check the effect of contamination the purity of these species was extrapolated to 100%. Table 1. List of the maximum relative systematic uncertainties of each individual source for v n,mk of π ± , K ± and p + p. The uncertainties depend on the transverse momenta. Percentage ranges are given to account for all centrality intervals.
The topological cuts were also varied to account for the V 0 and φ-meson reconstruction. These selection criteria were varied by (i) changing the reconstruction method for V 0 particles to an alternate technique that uses raw tracking information during the Kalman filtering stage (referred to as online V0 finder); (ii) varying the minimum radial distance from the primary vertex at which the V 0 can be produced from 5 cm to 10 cm; (iii) changing the minimum value of the cosine of pointing angle from 0.998 to 0.99; (iv) varying the minimum number of crossed TPC pad rows by the V 0 daughter tracks from 70 to 90; (v) changing the requirement on the minimum number of TPC space points that are used in the reconstruction of the V 0 daughter tracks form 70 to 90; (vi) requesting a minimum ratio of crossed to findable TPC clusters from 0.8 to 1.0; (vii) changing the minimum DCA of the V 0 daughter tracks to the primary vertex from 0.1 cm to 0.3 cm; (viii) changing the maximum DCA of the V 0 daughter tracks from 0.5 cm to 0.3 cm; (ix) requiring a minimum p T of the V 0 daughter tracks of 0.2 GeV/c.
In addition, the non-flow contribution was studied by (i) selecting like sign pairs of particles of interest and reference particles to decrease the effect from the decay of resonance particles; (ii) applying pseudorapidity gaps between the two subevents from |∆η| > 0.0 to |∆η| > 0.4.
Tables 1 and 2 summarise the maximum relative systematic uncertainties for each individual systematic source described above for all transverse momenta. The systematic uncertainties are expressed for each non-linear mode and particle species in a range to account for all centrality intervals in this article.

Results and discussion
In this section, the results of the p T -dependent non-linear flow modes v 4,22 , v 5,32 , v 6,33 and v 6,222 of identified particles are presented for various centrality intervals in Pb-Pb collisions at √ s NN = 5.02 TeV. We first present the centrality and p T dependence of v n,mk in section 6.1. The scaling properties of the non-linear flow modes are also discussed Table 2. List of the maximum relative systematic uncertainties of each individual source for v n,mk of K 0 S , Λ + Λ and φ-meson. The uncertainties depend on the transverse momenta and centrality interval. Percentage ranges are given to account for all centrality intervals. "N/A" indicates that a certain check was not applicable to the given particle of interest. If a source was checked and proved to have a negligible effect, the field is marked as "-". in this section. These results are compared with v n measurements for the same particle species in section 6.2. Finally, the comparison with two model calculations is shown in section 6.3. Note that in some of the following sections the same data are used in different representations to highlight the various physics implications of the measurements in each section.
6.1 Centrality and p T dependence of non-linear flow modes Figure 2 presents the magnitude of the non-linear mode for the fourth order flow coefficient, v 4,22 (p T ), for π ± , K ± , K 0 S , p + p, Λ + Λ and the φ-meson in a wide range of centrality intervals, i.e. 0-5% up to 50-60%. For the φ-meson, the results are reported from the 10-20% up to the 40-50% centrality interval, where v 4,22 can be measured accurately. The magnitude of this non-linear flow mode rises steeply with increasing centrality interval from 0-5% to 40-50% for all particle species. This increase is expected as v 4,22 reflects the contribution of the second order eccentricity, ε 2 , which increases from central to peripheral collisions, in v 4 [9,54]. For more peripheral collisions (i.e. 50-60%), the magnitude of v 4,22 does not increase further with respect to the neighbouring centrality interval (40-50%). This effect that was observed also in v n measurements [39,42] is probably due to the shorter lifetime of the produced system in more peripheral collisions, which prevents v 4,22 from developing further. Figure 3 presents the non-linear mode for the fifth order flow coefficient, i.e. v 5,32 (p T ), of π ± , K ± , K 0 S , p + p, and Λ + Λ for the same range of centrality intervals, i.e. 0-5% up to 50-60%. Statistical precision limits extending the measurements of non-linear flow modes of the φ-meson for n > 4. The measurements show a significant increase in the magnitude of this non-linear flow mode with increasing centrality percentile. This is due to the fact that v 5,32 (p T ) has a contribution from both ε 2 and ε 3 . It is shown in MC studies that ε 2 and to a smaller extent, ε 3 increase for peripheral collisions [9].       Figures 4 and 5 present the non-linear terms for the sixth order flow coefficient, i.e. v 6,33 (p T ) for π ± , K ± , K 0 S , p + p and Λ + Λ for the 0-5% up to 40-50% centrality intervals and v 6,222 (p T ) for π ± , K ± , p + p for the 0-5% up to 50-60% centrality intervals. As expected, measurements of v 6,222 (p T ) which probe the contribution of ε 2 , show an increase in the magnitude of this non-linear flow mode with increasing centrality percentile. On the other hand, the v 6,33 (p T ) measurements, which probe the contribution of ε 3 , present little to no dependence on centrality as previously observed for charged particles in [54]. In figure 6 the same data points are grouped by centrality interval to highlight how v 4,22 develops for a given centrality for various particle species as a function of p T . A clear mass ordering can be seen in the low p T region (i.e. p T < 2.5 GeV/c) for all collision centralities. This mass ordering arises from the interplay between radial flow and the initial spatial anisotropy, generated from both the geometry and the fluctuating initial energy density profile. This creates a depletion in the particle spectra at lower p T values which becomes larger in-plane than out-of plane due to the velocity profile. This naturally leads to lower v 4,22 (p T ) values for heavier particles [57,58,78]. Similarly, figures 7, 8 and 9 show the p T -differential v 5,32 , v 6,33 and v 6,222 , respectively, of different particle species for each centrality interval. A clear mass ordering is seen in the low p T region, (i.e. p T < 2.5 GeV/c), for v 5,32 (p T ) and to a smaller extent for v 6,33 (p T ) as well as for some centrality intervals of v 6,222 (p T ).
In addition, in the intermediate p T region (for p T > 2.5 GeV/c) the data points of figures 6-9 exhibit a particle type grouping. In particular, the data points form two groups, one for mesons and one for baryons with the values of v n,mk of the latter being larger. This particle type grouping was previously observed in v n measurements of various particle species [15-18, 39, 40, 42]. This grouping was explained in ref. [60] in the picture of particle production via quark coalescence indicating that flow develops at the partonic stage. In this picture, known as NCQ scaling, the flow of mesons (baryons) is roughly twice (thrice) the flow of their constituent quarks in the intermediate transverse momentum region [59,60]. The ALICE measurements show that this scaling at the LHC energies holds at an approximate level of 20% for v n [39,40,42]. Figures 10, 11, 12 and 13 present v 4,22 , v 5,32 , v 6,33 and v 6,222 , respectively, scaled by the number of constituent quarks (n q ) as a function of p T /n q for π ± , K ± , K 0 S , p + p, Λ + Λ and the φ-meson grouped in different centrality intervals. The scaling is consistent with the observations reported for higher order anisotropic flow coefficients [42]. It is seen that for the non-linear flow modes this scaling holds at an approximate level (±20%) for p T > 1 GeV/c, where quark coalescence is expected to be the dominant process.

Comparison with v n of identified particles
The comparison of the features discussed before i.e. mass ordering and particle type grouping between the non-linear and the anisotropic flow coefficient is of particular interest. Based on a naive expectation the mass ordering should develop quantitatively in a different way between the non-linear (i.e. due to the dependence on ε 2 2 ) and the anisotropic flow coefficient. In parallel, if coalescence is the dominant particle production mechanism in the intermediate p T region, one expects a similar grouping between v NL n and v n . Such a comparison could only be performed for v 4,22 (p T ) (this study) and the v 4 (p T ) measurements [42] and was done by taking the difference between pions and protons at a given p T in both modes and normalising it by the integrated flow of the corresponding mode for charged particles [41] . This comparison is shown in figure 14 for the 0-5% up to the 40-50% centrality interval. It can be seen that in the low p T region (p T < 2.5 − 3 GeV/c) where the mass ordering is prominent, the data points exhibit a general agreement for all centrality intervals. However, there is a hint that the relative ratio for v 4,22 is smaller than the one of the v 4 for p T < 0.8 GeV/c and for the centrality intervals 0-30%. If this difference and its centrality dependence persist for low values of p T , it could indicate that the hydrodynamic evolution is reflected differently in v 4 and v 4,22 and could be explained by the contribution of ε 2 2 . As stated earlier, the mass splitting is a result of an interplay of radial and anisotropic flow, leading to a stronger in-plane expansion compared to out-of-plane, and the particle thermal motion. Particles with larger mass have smaller thermal velocities, and are thus affected stronger by the difference between in-and out-ofplane expansion velocities, thus leading to the mass splitting of v n (p T ). The comparison of the p T dependence of v NL n and v n can therefore provide a unique opportunity to test this picture, as it would allow results for the cases of exactly the same average radial flow and temperature, but differing in anisotropic flow to be compared. On the other hand, in the intermediate p T region (p T > 2.5 GeV/c), the same comparison shows that the results are compatible in all centrality intervals within one standard deviation. This implies a similar particle type grouping between v 4 and v 4,22 which is in line with the expectation that quark coalescence affects both flow modes similarly.

Comparison with models
The comparison of various anisotropic flow measurements and hydrodynamic calculations are presented and discussed in great detail in [79][80][81]. A recent comparison between v n measurements reported by the ALICE collaboration [42] and two hydrodynamic calcula- tions from [81] shed new light on the initial conditions and the transport properties of the created system in Pb-Pb collisions. Both hydrodynamic calculations are based on iEBE-VISHNU [82], an event-by-event version of the VISHNU hybrid model [83] coupling 2 + 1 dimensional viscous hydrodynamics (VISH2+1) [84] to a hadronic cascade model (UrQMD). The initial conditions used for these calculations are described by AMPT [85] and T R ENTo [86], both with τ 0 =0.6 fm/c and T sw =148 MeV [87]. For AMPT initial conditions, constant values of specific shear viscosity over entropy density (η/s = 0.08, the lower limit conjectured by AdS/CFT) and bulk viscosity over entropy density (ζ/s = 0) are utilised. The version of the model that uses T R ENTo [86] initial conditions incorporates temperature dependent specific shear and bulk viscosities extracted from the global bayesian analysis [87]. The comparison between v n measurements and these two hydrodynamic calculations illustrates a qualitative agreement. This agreement between the data and the models depends on the particle species, transverse momentum range and centrality percentile. Overall, the AMPT model reproduces the measurements more accurately than the T R ENTo model [42]. In order to further investigate the performance of these two models in reproducing the v n measurements, and provide a quantitative comparison, the relative ratios between each model and the measurements of π ± , K ± and p + p are obtained. Table 3 summarises these relative ratios. The values represent the ranges across all centralities that each model is able to describe the measurements of v n for each particle species. Comparisons between the performance of the two models show that the AMPT calculations reproduce v 2 slightly better that T R ENTo. Both models reproduce the v 3 measurements relatively better than the v 2 , however AMPT performs better than T R ENTo. Finally, the comparison between the models and the v 4 measurements show that AMPT has an absolute better performance compared to T R ENTo. These values should be taken with caution as v 4 has larger uncertainties with respect to v 3 and v 2 .  π ± K ± p + p π ± K ± p + p π ± K ± p + p AMPT calculations 3-13% 0-16% 0-20% 0-8% 5-12% 0-4% 6-12% 5-12% 0-4% T R ENTo calculations 6-17% 0-19% 3-19% 2-15% 7-22% 0-11% 7-25% 16-28% 0-21% Table 3. List of minimum and maximum values of the fit with a constant function to relative ratios between data and each model for v n (n = 2, 3, 4) of π ± , K ± and p + p. The percentages show deviations of the fit from unity obtained for the 0-5% up to 40-50% centrality intervals.
To achieve additional constraints on the initial conditions and transport properties of the system and test the validity of these hydrodynamic models, a comparison is performed between the measured p T -dependent non-linear flow modes for π ± , K ± , p + p, K 0 S and Λ + Λ with the same two hydrodynamical calculations reported in [81]. Figures 15-18 present the comparison between the measurements and the two model predictions for the p T -differential v 4,22 , v 5,32 , v 6,33 and v 6,222 , respectively, for π ± , K ± and p + p and figures 19-21 present these comparisons for the p T -differential v 4,22 , v 5,32 and v 6,33 for K 0 S and Λ + Λ for the 0-10% up to 50-60% centrality interval (40-50% centrality interval for of v 4,22 reveals that T R ENTo reproduces the data very well from the 0-10% up to 30-40% centrality interval and fails to reproduce the measurements for the remaining more peripheral centrality intervals. On the other hand, AMPT overestimates the measurements from the 0-10% up to 30-40% centrality interval. For the 40-50% centrality interval, it reproduces the measurements for all particle species except π ± , where it slightly underestimates the results. For more peripheral collisions, it reproduces the K ± , p + p and Λ + Λ measurements and underestimates the results for π ± and K 0 S . In a similar attempt to the comparison between the v n measurement and the model calculation in table 3, the performance of these models were further studied for v n,mk by taking the relative ratios between each model and the measurements of π ± , K ± and p + p.  π ± K ± p + p π ± K ± p + p π ± K ± p + p π ± K ± p + p AMPT calculations 5-32% 2-30% 3-30% 3-28% 5-29% 1-65% 0-46% 0-46% 0-97% 6-52% 0-80% 0-118% T R ENTo calculations 0-30% 4-33% 0-21% 24-49% 33-97% 12-58% 0-43% 0-46% 0-95% 0-20% 0-34% 0-78% Table 4. List of minimum and maximum values of the fit with a constant function to relative ratios between the data and each model for v n,mk of π ± , K ± and p + p. The percentages show deviations of the fit from unity obtained for the 0-10% up to 50-60% (40-50% for v 6,33 ) centrality intervals.
to stress, however, that the non-linear flow modes have smaller magnitudes with respect to v n and any discrepancy between the models and the data becomes magnified in the ratios reported in table 4.
For v 5,32 , the comparison is different, with the T R ENTo predictions overestimating the measurements for all centrality intervals, and AMPT reproducing the data better than T R ENTo. The AMPT model overestimates the measurements from the 0-10% to 20-30% centrality interval. It underestimates the measurements of π ± , K ± and p + p for more peripheral collisions while it reproduces the measurements of K 0 S and Λ + Λ relatively well up to the 40-50% centrality interval. These comparisons are reflected in table 4 where AMPT performs on average 20-27% better than T R ENTo for π ± , K ± and p + p.
For v 6,33 , both models reproduce the data for the 0-10% centrality interval. For the 10-20% up to 30-40% centrality interval, AMPT reproduces the data while T R ENTo over- estimates the measurements. Finally, the comparison with v 6,222 shows an agreement between both models and the measurements of π ± , K ± and p + p at 0-10% up to 30-40% centrality intervals. 3 All in all, this study shows larger differences between the model calculations and the v n,mk measurements with respect to that of v n , indicating a larger sensitivity to the initial conditions and transport properties for the non-linear flow modes. As a result, it is useful to tune the input parameters of hydrodynamic models considering also the non-linear flow measurements.

Summary
In this article, the measurements of the non-linear flow modes, v 4,22 , v 5,32 , v 6,222 and v 6,33 are for the first time reported as a function of transverse momentum for different particle species, i.e. π ± , K ± , K 0 S , p + p, Λ+Λ and φ-meson. The results are presented in a wide range of centrality intervals from 0-5% up to 50-60% in Pb-Pb collisions at √ s NN = 5.02 TeV.
The magnitude of the non-linear flow modes, v n,mk , were obtained with a multi-particle correlation technique, namely the generic framework, selecting the identified hadron under study and the reference flow particles from different, non-overlapping pseudorapidity regions. The measured v 4,22 , v 5,32 and v 6,222 exhibit a distinct centrality dependence. This centrality dependence originates from the contribution of initial state eccentricity, ε 2 , as shown in eq. (1.5). As expected, v 6,33 does not exhibit a considerable centrality dependence since ε 3 quantifies primarily the event-by-event fluctuations of the initial energy density profile. This is supported by the relatively large magnitude of v 6,33 in the most-central collisions (0-5%). A clear mass ordering is observed in the low p T region (p T < 2.5 GeV/c). A closer comparison between v 4 and v 4,22 shows that this mass ordering seems slightly larger for v 4,22 than v 4 at very low p T (p T < 0.8 GeV/c). In the intermediate p T region (p T > 2.5 GeV/c), a particle type grouping is observed where the magnitude of the non-linear modes for baryons is larger than for mesons similar to observations in v n measurements. The NCQ scaling holds at an approximate level of ±20% within the current level of statistical and systematic uncertainties, similar to that of the anisotropic flow coefficients [42].  Figure 18. The p T -differential v 6,222 of π ± , K ± and p + p in the 0-10% up to 50-60% centrality intervals of Pb-Pb collisions at √ s NN = 5.02 TeVcompared with iEBE-VISHNU hybrid models with two different sets of initial parameters: AMPT initial conditions (η/s= 0.08 and ζ/s = 0) shown as solid bands and T R ENTo initial conditions (η/s(T) and ζ/s(T)) as hatched bands. The bottom panels show the difference between the measurements and each model. Statistical and systematic uncertainties are shown as bars and boxes, respectively.
The comparison of two models based on the iEBE-VISHNU hybrid model, with two different initial conditions (AMPT and T R ENTo) and transport properties shows that neither of the models is able to fully describe the measurements. The quality of the model description depends on the centrality percentile and particle species similar to the modeldata comparisons of the anisotropic flow coefficients [42]. The measurements are better predicted by the models in more central collisions. All in all, the model using AMPT initial conditions (η/s = 0.08 and ζ/s = 0) exhibits a magnitude and shape closer to the measurements. As a result, in order to further constrain the values of the transport properties and the initial conditions of the system, it is necessary to tune the input parameters of future hydrodynamic calculations attempting to describe these measurements.     Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.