Higher-order moments of the elliptic flow distribution in PbPb collisions at $\sqrt{s_\mathrm{NN}}$ = 5.02 TeV

The hydrodynamic flow-like behavior of charged hadrons in high-energy lead-lead collisions is studied through multiparticle correlations. The elliptic anisotropy values based on different orders of multiparticle cumulants, $v_{2}\{2k\}$, are measured up to the tenth order ($k$ = 5) as functions of the collision centrality at a nucleon-nucleon center-of-mass energy of $\sqrt{s_\mathrm{NN}}$ = 5.02 TeV. The data were recorded by the CMS experiment at the LHC and correspond to an integrated luminosity of 0.607 nb$^{-1}$. A hierarchy is observed between the coefficients, with $v_{2}\{2\}>v_{2}\{4\} \gtrsim v_{2}\{6\} \gtrsim v_{2}\{8\} \gtrsim v_{2}\{10\}$. Based on these results, centrality-dependent moments for the fluctuation-driven event-by-event $v_{2}$ distribution are determined, including the skewness, kurtosis and, for the first time, superskewness. Assuming a hydrodynamic expansion of the produced medium, these moments directly probe the initial-state geometry in high-energy nucleus-nucleus collisions.


Introduction
In high-energy nucleus-nucleus collisions, a hot and dense state of strongly interacting quarks and gluons is created, the so-called quark-gluon plasma (QGP).Clear evidence for this state was obtained at the BNL RHIC [1][2][3][4].Experiments at the CERN LHC, performed with nucleusnucleus collisions at nucleon-nucleon center-of-mass energies up to more than twenty times greater than available for the RHIC studies, confirmed the main RHIC results with much larger event samples and kinematic ranges.Unexpectedly, a similar behavior has been observed in smaller systems formed in high-multiplicity proton-proton and proton-nucleus collisions [5,6].One important feature of the QGP is its apparent collective hydrodynamic expansion.The initial spatial geometry of the overlap region of the colliding nuclei results in anisotropic pressure gradients that, in turn, are reflected in the azimuthal angle distribution of outgoing particles.This anisotropy has been used to perform detailed studies of the QGP formed at the LHC [7][8][9][10][11][12][13][14][15][16][17][18][19].Fluctuations in nuclear densities that appear due to the fluctuations in the position of nucleons within the incident nuclei also have a significant influence on the QGP expansion [8,12,[20][21][22].The effect of these geometry fluctuations on the observed particle anisotropy can give insight into the early-stage dynamics of the collisions [23][24][25].
Figure 1 schematically depicts a collision of two nuclei in the transverse X ′ − Y ′ plane that is perpendicular to the beam axis Z, in the laboratory frame.The impact parameter vector b, which is not experimentally accessible, connects the centers of the colliding nuclei in the transverse plane.The symmetry plane of the collision, denoted X − Y, defined in terms of the beam direction and b, is randomly oriented with respect to the X ′ − Y ′ plane.In this schematic figure, ϕ denotes the azimuthal angle of one of the outgoing particles.The anisotropic flow is quantified in terms of coefficients in the Fourier expansion of the ϕ dependence of the invariant yield of particles relative to the event-plane angle Ψ EP , where E is the energy of the particle, p T its transverse momentum, and y its rapidity.The eventplane angle Ψ EP is defined to be in the direction of maximum outgoing particle density [26] and is, on average, the same as the rotation angle of the symmetry plane with respect to the laboratory frame.The direction of Ψ EP fluctuates about the symmetry plane rotation angle because of fluctuations in the initial state geometry and resolution effects arising from finite particle multiplicities.The elliptic flow harmonic v 2 is the leading term of the Fourier series expansion of the azimuthal angle distribution in the event plane frame.
Many methods have been developed to study the anisotropic expansion of the system formed in nuclear collisions [26][27][28].The cumulant method was introduced in Refs.[29,30] based on multiparticle correlations which do not require an explicit measure of the event-plane angle.In this method, v n values can be extracted from c n {2k} cumulant coefficients, where n is the order of the Fourier expansion term and 2k is the number of particles used to determine the correlation.For example, v 2 {2}, v 2 {4}, and v 2 {6} are referred to as the 2-, 4-, and 6-particle cumulant based values for the elliptic flow harmonic, respectively.Nonflow, short-range correlations arising from jets and resonance decays, can be suppressed by correlating four or more particles.
It is then possible to determine the central moments of the event-by-event, fluctuation-driven, v 2 distribution based on the v 2 {2k} values.The 3rd, 4th, and 5th central moments are called the skewness (s), kurtosis (κ), and superskewness (p) respectively.The skewness describes the overall asymmetry of the v 2 distribution, while the superskewness gives a measure of the asymmetry of the tails.The kurtosis describes the peakedness of the center, or equivalently, the relative contribution (heaviness) of the distribution tails.In a hydrodynamic expansion the v 2 is proportional to the initial-state eccentricity ϵ 2 , and therefore ϵ 2 fluctuations produce non-Gaussian v 2 distributions [31].Thus, the central moments can be related to the initial state geometry and reveal its detailed structure.
In this paper, the cumulant based values are determined using the Q-cumulant method in which one calculates multiparticle cumulants in terms of moments of Q-vectors [32].Determining the cumulant based v 2 {2k} (k = 1, 2, 3 . . . ) values of the v 2 distribution is a way to study the fluctuation behavior related to the initial-state geometry [30].For Gaussian fluctuations, the higher central moments vanish and the v 2 {2k} values with k > 1 are all expected to have the same value [33].The measured higher-order v 2 {2k} (k = 2, 3, 4) values, however, show a fine splitting [34,35], , that can be attributed to a non-Gaussian behavior [31].In Ref. [36], it is noted that the main signature of non-Gaussian fluctuations is a nonzero skewness of the v 2 distribution.This is used to suggest, as described in Ref. [36], a hydrodynamic probe based on the ratio ).The basic premise of a hydrodynamic probe is that the observed azimuthal angle correlations that appear in the bulk medium can be related to the initial-state geometry.
We present for the first time a measurement of v 2 {10}, derived from ten-particle cumulants, for lead-lead (PbPb) collisions at a center-of-mass energy per nucleon pair of √ s NN = 5.02 TeV.The data used in this analysis correspond to an integrated luminosity of 0.607 nb −1 and were recorded with the CMS detector at the LHC.This enables the development of a new hydrodynamic probe using the v 2 {10} value and the corresponding lower-order v 2 {2k} values.Both the original and the new hydrodynamic probes should be independent of collision centrality if moments of the v 2 distribution with order higher than the skewness are negligible [36].Centrality refers to the percentage of the total inelastic hadronic nucleus-nucleus cross section [37], with 0% corresponding to the maximum overlap of the colliding nuclei.This connection between the centrality dependence of the hydrodynamic probes and higher order moments of the v 2 distribution is investigated by measuring both the centrality dependence itself and also by using the cumulants to extract the higher order moments directly.This paper is organized as follows.In Section 2 and in Appendix A, we derive formulas for the v 2 {10} harmonic and its statistical uncertainty.Section 3 and Appendices B, and C present derivations of the formulas used in this analysis to calculate the hydrodynamic probe ratios.In Section 4, we introduce and derive formulas for the standardized moments (denoted as γ exp i , i = 1, 2, 3) and the corrected moments (denoted as γ exp i, corr , i = 1, 2, 3) of the fluctuation-driven, event-by-event v 2 distribution.Section 5, 6, and 7 present the experimental details, results, and summary, respectively.The numerical values of the results for this analysis are tabulated in the HEPDATA record [38].

Analysis procedure
A version of the cumulant method, the Q-cumulant method, was introduced in Ref. [32] and is based on the flow vector Q n , with Q n = ∑ M j=1 e inϕ j .Here, n is the flow harmonic order, M denotes the event multiplicity, i.e., the number of analyzed tracks in the given event, and ϕ j is the azimuthal angle of the jth track measured in the laboratory frame.
The 2m-particle azimuthal angle correlators ⟨2m⟩ [32] are expressed through the flow vector ) .( 2) Here the summation is over all unique track indices.The weighted mean of the 2m-particle azimuthal angle correlator ⟨⟨2m⟩⟩ is defined as where the double brackets ⟨⟨. ..⟩⟩ denote a weighted average over all events within a given class [39].The weight is taken as the number of distinct 2m-particle combinations that can be formed for an event [39]: The Q-cumulant c n {2k} values are multivariate polynomial functions of the ⟨⟨2m⟩⟩ values.
Based on Eqs. ( 5) and ( 6), one finds where the a 2k coefficients can be determined using the recursion relationship In the ideal case, all 2k-particle cumulants c n {2k} with k > 1 result in the same flow harmonic value v n .But, non-Gaussian flow fluctuations and non-flow correlations result in a splitting of the c n {2k} values with k > 1.In this case, the v n {2k} values will depend on the cumulant orders as Details concerning the Q-cumulant method up to the k = 4 order can be found in Refs.[32,39].In the case where k = 5, the above equations give the relationship between the v n {10} magnitude and the corresponding c n {10} value, with and Formulas for the analytic calculations of the statistical uncertainties in the v n {2k} (k = 1, . . ., 4) magnitudes are given in Ref. [40].Appendix A discusses the analytical calculation of the statistical uncertainty for the v 2 {10} value.

Hydrodynamic probes
A compact formulation for the cumulant expansion is given through the formalism of generating functions.The Fourier-Laplace transform of the symmetry plane elliptic harmonic vector v 2 = v x e x + v y e y is ⟨e l•v 2 ⟩, with a vector variable l = l x e x + l y e y .Here, e x is the unit vector in the direction of the impact parameter, while e y is the unit vector perpendicular to b.The angle brackets ⟨⟩ indicate an event average that, for this analysis, involves an average over events within a centrality class.As the distribution of (v x , v y ) is symmetric under v y → −v y , the mean value of the x− and y−component of v 2 are ⟨v x ⟩ ≡ v 2 , and ⟨v y ⟩ = 0.The variances are σ 2 x = ⟨(v x − ⟨v x ⟩) 2 ⟩ and σ 2 y = ⟨v 2 y ⟩.Higher-order central moments up to the kurtosis are defined in Ref. [42] in terms of the components of the elliptic harmonic vector.The expressions found in Ref. [42] are reproduced here, expanded to include terms related to the superskewness: x s 12 , and In Ref. [36], the generating function of the cumulant is defined as ln⟨e l•v 2 ⟩.If one expands this generating function up to the 4th power in (l x , l y ) one gets the expression in terms of central moments of the (v x , v y ) distribution, that includes the σ, s, and κ terms [42].By extending this expansion up to the 5th power in (l x , l y ), we obtain One can follow the procedure for expanding the generating function of the v 2 {2k} cumulants [36] to derive expressions for the hydrodynamic probes that include higher-order moments of the v 2 distribution.Instead of using Eq. ( 7) from Ref. [36], in this analysis Eq. ( 13) is used.The difference is that in the earlier work the expansion is performed up to the skewness, while the current analysis additionally covers the 4th and 5th moments.Then, for the higher-order v 2 {2k} (k = 2, . . ., 5) values, we find , and ( 16) If one neglects the 5th-order moments p, as well as the terms with the multiplier (σ 2 y − σ 2 x ), then the hydrodynamic probe proposed with Eq. ( 14) in Ref. [36] is given as Here, one can see that neglecting the kurtosis κ moment, with respect to the skewness, would reduce the right-hand side of Eq. ( 18) to a constant value of 1/11.Appendix C shows the relation that corresponds to Eq. ( 18) if none of the terms in Eqs. ( 14)-( 16) are neglected.
The approximation given by Eq. ( 18) can be tested by expressing the right-hand side in terms of the measurable quantities v 2 {2k}, with k = 2, 3, and 4. Based on Taylor expansion of the v 2 {2k} 2 values, as shown in Appendix B, one finds if terms with (σ 2 y − σ 2 x ) are neglected.

In the above relation, h 1 denotes the hydrodynamic probe and h
Taylor 1 denotes the corresponding Taylor expansion of this probe as expressed in terms of the measured v 2 {2k} values.
By employing the v 2 {10} harmonic, we define a new hydrodynamic probe h 2 as In the derivation of Eq. ( 20), none of the terms in Eqs. ( 15)-( 17) have been neglected.If one neglects the 5th p moment, then the right-hand side of the Eq. ( 20) will reduce to the constant value of 3/19.A deviation from this constant value indicates that the p moment cannot be neglected.
Again, as in the case with the derivation of the Eq. ( 19), one can express the right-hand side of Eq. ( 20) through the second power of the higher-order v 2 {2k} terms, with k = 3, 4, 5 (as shown in Appendix B): If one neglects the term with the multiplier (σ 2 y − σ 2 x ), Eq. ( 21) reduces to where h 2 denotes the new hydrodynamic probe and h Taylor 2 denotes the Taylor expansion of this probe, as expressed in terms of the measured v 2 {2k} values.
In Section 6, it will be shown that the term proportional to the asymmetry of the fluctuations (σ 2 y − σ 2 x ) is negligible.The same conclusion is found in Ref. [36].

The standardized and corrected moments
The standardized skewness, γ 1 , and standardized kurtosis, γ 2 , are defined as The corresponding expressions given in terms of the measured v 2 {2k} values, are derived in Ref. [42] as and The standardized superskewness, γ 3 , is defined as In order to express the standardized superskewness through the measured v 2 {2k} values with k = 1, . . ., 5, one can use the expansion of the fifth power of the v 2 {2k} value (as discussed in Appendix B).The resulting expression for the standardized superskewness, expressed in terms of measurable quantities, is The standardized moments, γ exp i (i = 1, 2, 3), have contributions from higher-order moments of the v 2 distribution that are not negligible (as shown by Eq. ( 13) in Ref. [42]).For example, if higher-order moments are considered, then for the standardized skewness we obtain additional terms O N and O D , with where Corrected moments, which are free of contributions from higher-order moments, can be defined if These conditions are satisfied with the elliptic power distribution [43] if the eccentricity parameter ϵ 0 in this model satisfies the condition ϵ 0 ≤ 0.15 and a linear flow response is assumed.Mathematical details that support the validity of the conditions given by Eq. ( 29) are presented in Appendix D. With these conditions, the corrected moments can be expressed in terms of the measured v n {2k} (k = 1, . . ., 5) values.As an example, the corrected skewness can be expressed as If one replaces the powers that appear in Eq. ( 30) with the expansions given in Appendix B, one finds where O N and O D are given by We note that the contribution from the term proportional to (σ 2 y − σ 2 x ) as well as the contributions from moments of order higher than 5 are negligible.
The same correction procedure can also be applied to the standardized kurtosis and the standardized superskewness.The corrected form of these moments are given by

CMS detector, event selection, and systematic uncertainties
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections.Forward calorimeters (HF), made of steel and quartz-fibres, extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors.Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid.A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [44].
Events of interest are selected using a two-tiered trigger system.The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz [45].The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to about 5 to 8 kHz before data storage [46].
The data analyzed in this paper, before applying the selection described below, consist of 4.27×10 9 minimum bias lead-lead (PbPb) collision events at √ s NN = 5.02 TeV, collected in 2018 with an integrated luminosity of 0.607 nb −1 [47,48].The minimum bias events are triggered by requiring signals above readout thresholds of 3 GeV in each of the two HF calorimeters [46].Further selections are applied to reject events from beam-gas interactions and nonhadronic collisions [49].The events are also required to have at least one reconstructed primary vertex (based on two or more tracks) within a distance (|z vtx |) of 15 cm from the nominal interaction point along the beam axis.If multiple vertices are found in an event then the primary vertex is selected as the one with the highest track multiplicity in the event.The total energy deposited in both HF calorimeters is used for determination of the event centrality [50].
Only tracks that satisfy the high-purity selection criteria described in Ref. [49] are used in the analysis.In addition, a reconstructed track is only considered as a candidate track from the primary vertex if the separation along the beam axis between the track and the vertex (d z ), and the track-vertex impact parameter measured transverse to the beam (d xy ), are both less than three times their respective uncertainties σ z and σ xy .The relative uncertainty in the transverse momentum (σ p T /p T ) measurement is required to be less than 10%.In addition, tracks must have 11 or more hits (n hits ) along their trajectory in the pixel and strip tracking detectors combined.The relative χ 2 of the track fit (χ 2 /N dof /N layers ) is required to be less than 0.18, where N dof is the number of degrees of freedom during the track fitting and N layers is the number of tracking layers used.
Each reconstructed track is weighted by the inverse of the efficiency factor w(η, p T , centrality) calculated based on Monte Carlo HYDJET 1.9 [51].The efficiency factor accounts for the recon-struction efficiency ϵ and the fraction of misreconstructed tracks f misrec (η, p T , centrality), with w = ϵ/(1 − f misrec ) [52].The response of the CMS detector to the generated events is simulated with GEANT4 [53].
The systematic uncertainties are determined by varying the vertex selection, the track selection, the centrality determination, and the efficiency correction.For each of the indicated sources of systematic uncertainty, an absolute difference between results from the varied and nominal values are used to obtain the uncertainty.The vertex selection uncertainty is found by performing the analysis with |z vtx | < 3 cm and with 3 < |z vtx | < 15 cm.The uncertainty related to the track selection is determined by comparing tight (d z /σ z < 2, d xy /σ xy < 2, σ p T /p T < 0.05 and χ 2 /N dof /N layers < 0.15) criterion and loose (d z /σ z < 5, d xy /σ xy < 5, σ p T /p T < 0.15 and χ 2 /N dof /N layers < 0.18) criterion for tracks with the results obtained with the nominal values.The finite HF energy resolution results in small migration of events across the centrality bin boundaries.The centrality calibration is varied to calculate the effect of the minimum bias event selection efficiency of the HF calorimeters [54].Differences in the p T spectra between the data and the Monte Carlo simulation can affect the efficiency factor.To estimate the systematic uncertainty arising from these differences, the efficiency factor w is varied by ± 5%.The final systematic uncertainty is taken as a quadratic sum of the four indicated sources, with the track selection found to have the dominant contribution.

Results
It was previously shown in Ref. [36] that the presence of non-Gaussian fluctuations in the initial-state energy density leads to a fine splitting between the higher-order v n {2k} cumulants.The v 2 {2k} (k = 1, . . ., 5) values are presented in Fig. 2 as functions of centrality in PbPb collisions at √ s NN = 5.02 TeV.These results are obtained from charged particles with |η| < 2.4 and 0.5 < p T < 3.0 GeV/c.The statistical uncertainties are about two orders of magnitude smaller than the systematic ones.A clear splitting between the v 2 {2} and higher-orders cumulant based v 2 {2k} (k = 2, . . ., 5) values is visible.The difference is attributed to flow fluctuations, with v 2 {2} 2 ≈ v 2 {2k} 2 + 2σ 2 v for k > 1 [36,55], where σ 2 v is the v 2 variance.From Fig. 2, it is clear that the flow fluctuations become larger going to more peripheral collisions (higher centrality percentages).
In Fig. 2, the ordering and a fine splitting between the v 2 {2k} (k = 2, . . ., 5) values are not clearly visible.The v 2 {2k} (k = 1, . . .5) values need to satisfy the ordering relationship v 2 {2k} > v 2 {2(k + 1)} in order to define the hydrodynamic probes given by Eqs.(19) and (22).The splitting between the higher-order cumulant based becomes much finer than the splitting with respect to the lowest-order cumulant based v 2 {2} value.Figure 3 shows the corresponding relative differences (v 2 {2k} − v 2 {10})/v 2 {10} (k = 1, . . ., 4) as functions of centrality.The magnitudes of these differences are plotted using a logarithmic scale.The fine splitting between cumulants of different orders as well as the expected hydrodynamic ordering of the splittings is clearly seen.The relative difference between adjacent v 2 {2k} − v 2 {10} values decreases by about an order of magnitude for each increment in k.
The measured cumulants of different orders are used to calculate the hydrodynamic probes given by the left-hand sides of Eqs. ( 19) and (22).Figure 4 displays these distributions with closed symbols.Open symbols in the same figure show the right-hand sides of Eqs. ( 19) and (22).These are also constructed using the measured v 2 {2k} (k = 1, . . ., 5) cumulants.v 2 {2k} values, statistical and systematic uncertainties quickly increase as peripheral collisions are approached.Also, in the region of very central collisions, because of the small magnitudes of the cumulant based v 2 {2k} values, the statistical and systematic uncertainties are larger and, therefore, statistical and systematic uncertainties of the hydrodynamic probes are increased.The magnitudes of the statistical uncertainties of the h 2 distribution are larger with respect to those corresponding to the h 1 distribution because higher-order cumulants are involved.After the 0-5% centrality bin, where the measurement uncertainties are large, the h 1 and h 2 values show an increasing trend going to more peripheral events.Based on an event-by-event measurement of the v 2 distribution, it was reported in Ref. [34] that the h 1 has a value of 0.143 ± 0.008 (stat) ± 0.014 (syst) for 20-25% central events, and increases to 0.185 ± 0.005 (stat) ± 0.012 (syst) as the centrality increases to 55-60%.Using the Q-cumulant method, Ref. [35]   Taking into consideration differences in acceptance, the current analysis is consistent with the results of previous analyses [34,35], although with smaller uncertainties.There is a very good agreement between the distributions obtained by the expressions given on left-hand and righthand sides of Eq. ( 19) and Eq. ( 22) in Fig. 4.This indicates the importance of including higherorder terms in the Taylor expansion of the v 2 generating function.Figure 5 shows ratios between the left-and right-hand sides of the hydrodynamic probes given by Eq. ( 19) (left plot) and Eq. ( 22) (right plot) as functions of centrality.The two sides differ by less than 1.3% and 0.3% in the cases of Eqs. ( 19) and (22), respectively.Thus, based on Eqs.(C.2) and ( 21), which differ from the Eqs.( 19) and (22) only by the term that is proportional to (σ 2 y − σ 2 x ), one can conclude that the contribution of this term is quite small and can be neglected.Although with different acceptance and differences in the applied methodology, the obtained results for the experimentally measured standardized skewness are in a fair agreement with the results obtained in Ref. [34].The γ exp 1 is negative over the entire analyzed centrality range.This is a consequence of the v 2,x distribution having a long tail on the low v 2,x side (as shown in Fig. 1 of Ref. [36]).The v 2,y distribution is symmetric and thus its skewness is equal to zero.The γ exp 1 had been predicted to become more negative as the centrality percentile increases [36].This measurement confirms the prediction.

Centrality (%)
Similarly to the case of the skewness, although with different acceptance and a difference in the applied methodology to obtain the v 2 {2k} values, the γ exp 2 values are in reasonable agreement with the results presented in Ref. [42].The γ exp 2 values are positive over nearly the entire analyzed centrality range.The exception is the centrality range of 10-20%.The sign of γ exp 2 is driven by the mean eccentricity ϵ 0 , and is predicted to be negative for ϵ 0 < 0.28 and positive for ϵ 0 > 0.29 [42].Our results are qualitatively in an agreement with this prediction.If the v 2 values were to reflect a purely linear response to the initial-state eccentricity ϵ 2 , the corresponding kurtosis would equal the kurtosis of the initial ϵ 2 fluctuations.Although a fully linear response is not expected, calculations predict that deviations from a Gaussian behavior for the skewness and kurtosis, although significantly reduced by the hydrodynamic evolution, are still influenced by the early stage hydrodynamics [36,56].
Except for collisions with centrality less than 25%, where it is either positive or vanishes, the γ exp 3 moment is negative, with its absolute magnitude increasing towards more peripheral collisions.This is the first measurement of this moment of the v 2 distribution.Without the superskewness, it would be impossible to describe the centrality dependence observed for the new hydrodynamic probe h 2 .
In addition to the experimental results for the standardized skewness, kurtosis, and superskewness, the corresponding corrected moments are also presented.Except for the superskewness, the corrected skewness and kurtosis have larger slopes with respect to the standardized ones.The corrected moments give additional constraints on models of the initial-state geometry.

Summary
The elliptic anisotropy values based on different orders of multiparticle cumulants, v 2 {2k} (1 ≤ k ≤ 5), are determined as functions of centrality in lead-lead collisions at a center-ofmass energy per nucleon pair of √ s NN = 5.02 TeV, with an integrated luminosity of 0.607 nb −1 .The v 2 {10} value is determined for the first time.A fine splitting is observed between the harmonic values based on different cumulant orders, with The ordering of the v 2 {2k} values is consistent with a hydrodynamic evolution of the quarkgluon plasma (QGP).This splitting is attributed to a non-Gaussian behavior in the event-byevent fluctuations of the v 2 distribution, leading to nonzero values of the skewness, kurtosis, and superskewness.The splitting becomes finer as the k value increases, with the difference between the adjacent v 2 {2k} values decreasing by about an order of magnitude for each increment.The standardized magnitude of the v 2 moments are presented, together with their corrected values, where contributions from higher-order moments (up to the 5th moment) are removed.The large data set of lead-lead collisions collected by the CMS experiment enables a precise measurement of the hydrodynamic probe h 1 as a function of centrality, where ).A strong centrality dependence is observed, with values slowly increasing going to more peripheral collisions.This contrasts with an earlier hydrodynamic expectation that had taken the skewness of the initial-state geometry as the main source of non-Gaussian fluctuations.In that case, the ratio was not expected to depend on centrality.Based on the first v 2 {10} measurements, a new hydrodynamic probe is introduced that gives an even more precise measure of the initial-state geometry assuming a hydrodynamic evolution of the QGP.The new probe, h 2 , defined as ), is also found to have a centrality dependence, with a shape similar to the h 1 results.The centrality dependence of both ratios can be understood in terms of the evolving shape of the interaction region with centrality.Based on these results, centrality-dependent moments for the fluctuation-driven event-by-event v 2 distribution are determined, including the skewness, kurtosis and, for the first time, superskewness.The results provide basic input for a precision test of models that assume a hydrodynamic expansion of the QGP.

C Expanding the h 1 hydrodynamic probe to include terms up to the fifth central moment
The hydrodynamic probe, h 1 , was introduced in Ref. [36], but in this earlier reference its Taylor expansion was only shown to the 3rd moment of the (v x , v y ) distribution.If none of the terms from Eqs. ( 14)-( 16) are neglected, the right-hand side of Eq. ( 18) becomes where , and and 3 F2 is the regularized hypergeometric function.For the case when the ellipticity parameter tends to zero, ϵ 0 → 0, the upper formulas give the following ratios between moments of the initial anisotropy distribution: By assuming linear response between the flow and the initial anisotropy, the relations given by Eq. (D.3) are also well satisfied for the corresponding moments of the v n distribution.Furthermore, these relations are considerably transferred to relations between the corresponding central moments of the v n distribution given by Eq. ( 29).We tested these conditions numerically and found that they are satisfied to a reasonable degree in the range of the eccentricity parameter ϵ 0 ≤ 0.15.

Figure 1 :
Figure 1: A schematic view of a non-central nucleus-nucleus collision in the transverse plane.

Figure 2 :
Figure 2: The v 2 {2k} (k = 1, . . ., 5) values as functions of centrality in PbPb collisions at √ s NN = 5.02 TeV.The vertical sizes of the open boxes denote the systematic uncertainties.Statistical uncertainties are negligible compared to the marker size.Points are plotted at the center of the respective centrality ranges.The markers are displaced horizontally for better visibility.

Figure 3 :
Figure 3: The relative differences (v 2 {2k} − v 2 {10})/v 2 {10} (k = 1, . . ., 4) as functions of centrality in PbPb collisions at √ s NN = 5.02 TeV.The vertical sizes of the open boxes denote the systematic uncertainties.Statistical uncertainties are negligible compared to the marker size.Points are plotted at the center of the respective centrality ranges.

Figure 5 :
Figure 5: The ratios between the hydrodynamic probes and their Taylor expansions.The ratios are plotted as functions of centrality in PbPb collisions at √ s NN = 5.02 TeV.The vertical sizes of the open boxes denote systematic uncertainties.Statistical uncertainties are negligible compared to the marker size.Points are plotted at the center of the respective centrality ranges.

Figure 6 :
Figure 6: The magnitudes of the measured (closed circles) standardized skewness γ exp 1 (upper), standardized kurtosis γ exp 2 (middle), and standardized superskewness γ exp 3 (lower) as functions of centrality in PbPb collisions at √ s NN = 5.02 TeV.The magnitudes of the corrected skewness γ exp 1 ,corr (upper), corrected kurtosis γ exp 2, corr (middle), and corrected superskewness γ exp 3, corr (lower) are presented with the open circles.The bars (the vertical sizes of the open boxes) denote statistical (systematic) uncertainties.Points are plotted at the center of the respective centrality ranges.

Figure 6
Figure 6 displays the distributions of the measured standardized skewness γ exp 1 , standardized kurtosis γ exp 2 , and standardized superskewness γ exp 3 as functions of centrality.The distributions of the corrected skewness γ exp 1, corr (upper), corrected kurtosis γ exp 2, corr (middle), and corrected superskewness γ exp 3, corr (lower) are presented with the open blue circles.The systematic uncertainties are generally larger than the statistical ones.The values of the corrected moments the NKFIH research grants K 124845, K 124850, K 128713, K 128786, K 129058, K 131991, K 133046, K 138136, K 143460, K 143477, 2020-2.2.1-ED-2021-00181, and TKP2021-NKTA-64 (Hungary); the Council of Science and Industrial Research, India; ICSC -National Research Center for High Performance Computing, Big Data and Quantum Computing, funded by the EU NexGeneration program (Italy); the Latvian Council of Science; the Ministry of Education and Science, project no.2022/WK/14, and the National Science Center, contracts Opus 2021/41/B/ST2/01369 and 2021/43/B/ST2/01552 (Poland); the Fundac ¸ão para a Ciência e a Tecnologia, grant CEECIND/01334/2018 (Portugal); the National Priorities Research Program by Qatar National Research Fund; MCIN/AEI/10.13039/501100011033,ERDF "a way of making Europe", and the Programa Estatal de Fomento de la Investigaci ón Científica y Técnica de Excelencia María de Maeztu, grant MDM-2017-0765 and Programa Severo Ochoa del Principado de Asturias (Spain); the Chulalongkorn Academic into Its 2nd Century Project Advancement Project, and the National Science, Research and Innovation Fund via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, grant B37G660013 (Thailand); the Kavli Foundation; the Nvidia Corporation; the SuperMicro Corporation; the Welch Foundation, contract C-1845; and the Weston Havens Foundation (USA).