Energy dependence and fluctuations of anisotropic flow in Pb-Pb collisions at $\mathbf{\sqrt{{\textit s}_\text{NN}}} = \mathbf{5.02}$ and $\mathbf{2.76}$ TeV

Measurements of anisotropic flow coefficients with two- and multi-particle cumulants for inclusive charged particles in Pb-Pb collisions at $\sqrt{{\textit s}_\text{NN}} = 5.02$ and 2.76 TeV are reported in the pseudorapidity range $|\eta|<0.8$ and transverse momentum $0.2<p_\text{T}<50$ GeV/$c$. The full data sample collected by the ALICE detector in 2015 (2010), corresponding to an integrated luminosity of 12.7 (2.0) $\mu$b$^{-1}$ in the centrality range 0-80%, is analysed. Flow coefficients up to the sixth flow harmonic ($v_6$) are reported and a detailed comparison among results at the two energies is carried out. The $p_\text{T}$ dependence of anisotropic flow coefficients and its evolution with respect to centrality and harmonic number $n$ are investigated. An approximate power-law scaling of the form $v_n(p_\text{T}) \sim p_\text{T}^{n/3}$ is observed for all flow harmonics at low $p_\text{T}$ ($0.2<p_\text{T}<3$ GeV/$c$). At the same time, the ratios $v_n/v_m^{n/m}$ are observed to be essentially independent of $p_\text{T}$ for most centralities up to about $p_\text{T} = 10$ GeV/$c$. Analysing the differences among higher-order cumulants of elliptic flow ($v_2$), which have different sensitivities to flow fluctuations, a measurement of the standardised skewness of the event-by-event $v_2$ distribution $P(v_2)$ is reported and constraints on its higher moments are provided. The Elliptic Power distribution is used to parametrise $P(v_2)$, extracting its parameters from fits to cumulants. The measurements are compared to different model predictions in order to discriminate among initial-state models and to constrain the temperature dependence of the shear viscosity to entropy-density ratio.


Introduction
The primary goal of ultra-relativistic heavy-ion collisions is to study the properties of QCD matter at extremely high temperatures and/or densities and to understand the microscopic dynamics from which these properties arise, especially in the non-perturbative regime. The study of anisotropies in the azimuthal distribution of produced particles, commonly called anisotropic flow, has contributed significantly to the characterization of the system created in heavy-ion collisions [1][2][3][4][5]. According to the current paradigm of bulk particle production, anisotropic flow is determined by the response of the system to its initial spatial anisotropies. Initial-state spatial anisotropies come in turn from both the geometry of the collision and fluctuations in the wave function of the incident nuclei [3][4][5][6][7][8]. The significant magnitude of anisotropic flow is interpreted as evidence of the formation of a strongly-coupled system, which can effectively be described as a fluid with very low shear viscosity to entropy-density ratio (η/s) [9].
Anisotropic flow is quantified by the coefficients v n of a Fourier series decomposition of the distribution in azimuthal angle ϕ of final-state particles [10] where Ψ n corresponds to the symmetry plane angle of order n. The dominant flow coefficient in noncentral heavy-ion collisions is the second flow harmonic (v 2 ), called elliptic flow, which is mostly a result of the average ellipsoidal shape of the overlapping area between the colliding nuclei, whereas higher harmonics originate from initial-state fluctuations. For transverse momenta p T 3 GeV/c, anisotropic flow is thought to be quantitatively determined by the whole evolution of the system, including the phase of hadronic rescatterings that takes place after chemical freeze-out [11]. Flow coefficients have been shown to be sensitive not only to initial-state anisotropies, but also to the transport parameters (such as shear and bulk viscosity [12,13]) and the equation of state of the system, and they have therefore been used to constrain these properties [14,15]. However, given the different heterogeneous phases that the system is believed to undergo, it has not been possible so far to simultaneously constrain the large number of model parameters, although attempts have been made [16,17].
In this regard, the energy dependence of anisotropic flow has been shown to provide additional discriminating power over initial-state models and temperature dependence of transport parameters [18,19]. In fact, some theoretical uncertainties in the determination of anisotropic flow coefficients are expected to partially cancel in the ratios of v n coefficients measured at different collision energies, such as those on the choice of initial-state model or on the absolute value of η/s. These ratios would then effectively constrain the variations with collision energy and, therefore, system temperature of the parameters to which anisotropic flow is most sensitive.
It is known that the magnitude of anisotropic flow, being approximately proportional to the initial-state spatial anisotropy [20], fluctuates from collision to collision even for fixed centrality [6,[21][22][23][24], and that its probability distribution function (p.d.f.) P(v n ) is to a first approximation Bessel-Gaussian [1,25], i.e. the product of a modified Bessel function and a Gaussian function. It has been pointed out that small deviations from a Bessel-Gaussian shape are to be expected independently from the details of initial-state fluctuations [26][27][28]. Evidence of such small deviations has been previously reported [29]. These deviations are due to first order to the flow p.d.f. having a finite skewness. Its quantitative determination would therefore improve the characterization of these deviations. For dimensional reasons, it is convenient to use a standardised skewness (γ 1 ), defined as [30] where v n {RP} refers to the anisotropic flow with respect to the reaction plane Ψ RP , i.e. the plane spanned by the impact parameter and the beam axis, and the brackets · · · indicate an average over all events. It is worthwhile to note that the symmetry planes Ψ n do not generally coincide with Ψ RP because of initial-state fluctuations.
A robust experimental method to quantify flow fluctuations is to measure v n with multi-particle cumulants, which have different sensitivities to the moments of the underlying flow p.d.f.
The number in curly brackets indicates the order of the cumulant.
For elliptic flow, a large difference between v 2 {2} and v 2 {4} and approximately equal values of the higher order cumulants ) have been previously observed [29,31], which is indeed consistent with an approximately Bessel-Gaussian flow p.d.f.. However, a fine-splitting of a few percent among the higher order cumulants has also been reported [29], which is thought to be determined by the residual deviations from Bessel-Gaussian shape, in particular a non-zero skewness. A negative value of γ 1 , which corresponds to P(v 2 ) being left skewed, is expected [27] from the necessary condition on the initial-state eccentricity ε 2 < 1, which acts as a right cutoff on P(v 2 ). The Elliptic Power distribution, proposed in [26,27], was motivated mainly by this observation and it was shown to provide a good description of P(v 2 ) in a wide centrality range [32]. Moreover, γ 1 has been predicted to increase in absolute value from central to peripheral collisions [30], being roughly proportional to v 2 {RP} and being inversely proportional to the square root of the system size [28]. γ 1 can be estimated from the fine-splitting among two-and multi-particle cumulants [30] It is denoted as γ exp 1 to emphasize that it does not exactly match the definition of γ 1 given in Eq. 2, although the two have been estimated to coincide within a few percents [30]. The derivation of Eq. 7 relies on a Taylor expansion of the generating function in powers of the moments, truncated at the order of the skewness. It is experimentally possible to test the validity of this approximation through the universal equality that it implies [30] The precision up to which this equality holds depends on the residual contribution of higher central moments of the flow p.d.f., e.g. the kurtosis, to the multi-particle cumulants.
At high p T (p T 10 GeV/c) the dominant mechanism that determines azimuthal anisotropies of the produced final-state particles is thought to be path-length dependent energy-loss of highly energetic partons [33][34][35]. Although several experimental observations, such as jet azimuthal anisotropies [36,37], are consistent with this hypothesis, the details of the process are largely unconstrained and measurements of anisotropic flow of high-p T particles can help in this regard. Although the mechanism that determines it is fundamentally different, the origin of anisotropic flow at high p T is common to the one at low p T : initial-state geometry and its event-by-event fluctuations. Measurements reported in [38] seem to confirm this interpretation.
Recent CMS results on non-Gaussian elliptic flow fluctuations [39] appeared during the writing of this article. Numerical data are not yet available, but the observations seem to be essentially compatible with our measurements and their conclusions agree with those of this article.

Data sample and analysis methods
The sample of Pb-Pb collisions used for this measurement was recorded with the ALICE detector [40,41] in November and December 2015 (2010), during the Run 2 (Run 1) of the LHC, at a centre of mass energy per nucleon of √ s NN = 5.02 (2.76) TeV. The detectors used in the present analysis are the Inner Tracking System (ITS) and Time Projection Chamber (TPC), for primary vertex determination and charged particle tracking, and the V0 detector, for symmetry plane determination, centrality estimation [42] and trigger. The trigger conditions are described in [40]. About 78.4 × 10 6 (12.6 × 10 6 ) minimumbias events in the centrality range 0-80%, corresponding to an integrated luminosity of 12.7 µb −1 (2.0 µb −1 ), with a reconstructed primary vertex position along the beam direction (z vtx ) within ±10 cm from the nominal interaction point, passed offline selection criteria [40] for the data sample at √ s NN = 5.02 (2.76) TeV. Centrality is determined from the measured amplitude in the V0, which is proportional to the number of charged tracks in the corresponding acceptance (2.8 < η < 5.1 for V0A and −3.6 < η < −1.7 for V0C).
Charged tracks with transverse momentum 0.2 < p T < 50 GeV/c and pseudorapidity |η| < 0.8 are used in the present analysis. These tracks are reconstructed using combined information from the ITS and TPC. A minimum number of TPC space points of 70 (out of 159) is required for all tracks, together with a χ 2 per TPC space point (χ 2 TPC ) in the range 0.1 < χ 2 TPC < 4. A minimum number of 2 ITS hits, of which at least one in the two innermost layers, is required, together with a χ 2 per ITS hit per degree of freedom (χ 2 ITS ) smaller than 36. Only tracks with a distance of closest approach (DCA) to the primary vertex position less than 3.2 cm in the beam direction and 2.4 cm transverse to it are used. These track selection criteria ensure an optimum rejection of secondary particles and a p T resolution better than 5% in the p T range used in the present analysis [40].
Anisotropic flow coefficients are measured with the Q-cumulant method [43], using the implementation proposed in [44]. Track weights (w) are used in the construction of the Q-vectors, in order to correct for non-uniform reconstruction efficiency and acceptance where M is the charged track multiplicity, n the harmonic and m an integer exponent of the weights. After applying track weights, the effects due to non-uniformities in azimuthal acceptance, which would introduce a bias in the measured flow coefficients, are observed to be negligible. This is evaluated by measuring the event-averaged values of the real and imaginary part of Q n , which are consistent with zero. Multi-particle cumulants are measured on an event-by-event basis and then, in order to minimise statistical fluctuations, averaged over all events using the corrected charged track multiplicity as a weight, following the procedure proposed in [43]. All observables are computed in small centrality bins (1%) and then integrated, when limited size of the data sample makes it necessary, in wider centrality intervals using the charged particle yield in each 1% centrality bin as weight. This avoids that the event weighting procedure, based on multiplicity, would introduce a bias in the average centrality within a large centrality bin, since multiplicity varies with centrality.
For p T -integrated results, the m-particle cumulants are calculated using all tracks within given p T range, while for p T -differential results one particle at a given p T is correlated with m − 1 particles in the full p T range (0.2 < p T < 50 GeV/c). In terms of reference (c n {m}) and differential (d n {m}) cumulants, as defined in [43], the flow coefficients are measured as For two-particle correlations, a separation in pseudorapidity between the correlated particles (∆η) is applied in order to suppress short-range azimuthal correlations which are not associated to the symmetry planes, usually called 'non-flow'. These correlations arise from jets, mini-jets and resonance decays. For flow coefficients of higher order (v n {m > 2}), non-flow contribution has been previously found to be negligible in Pb-Pb collisions [24,31]. Results corresponding to |∆η| > 1 (denoted with v n {2, |∆η| > 1}) are obtained with the two-particle cumulant correlating tracks from opposite sides of the TPC acceptance, −0.8 < η < −0.5 and 0.5 < η < 0.8. Results corresponding to |∆η| > 2 (and reported as v n {2, |∆η| > 2}) are obtained with the scalar product method [45], correlating all tracks at mid-rapidity (|η| < 0.8) with the n-th harmonic Q-vector Q V0A n calculated from the azimuthal distribution of the energy deposition measured in the V0A detector [2] v n {2, |∆η| > 2} = u n,0 Q V0A* where u n,0 = e inϕ is the unit flow vector from charged particle tracks at mid-rapidity and Q n,1 is computed from the same type of tracks according to Eq. 9. Both methods have their own limitations and thus are complementary to each other: v n {2, |∆η| > 2} can be reliably employed only up to the fourth harmonic, because of the finite azimuthal segmentation of the V0 detectors (8 sectors in 2π), while v n {2, |∆η| > 1} suffers from bigger statistical uncertainties, due to the limited acceptance from which tracks are selected, and bigger non-flow contribution for p T > 10 GeV/c.
The systematic uncertainties are evaluated by varying the track and event selection criteria and comparing the variation in the flow coefficients relative to the default results. The absolute value of the variation itself is assigned as a systematic uncertainty if it is considered significant according to the Barlow criterion [46]. Different track quality variables are varied: number of TPC space points, χ 2 TPC and χ 2 ITS , fraction of shared TPC space points and number of ITS hits. For each of these, the default values are varied in order to increase the fraction of excluded tracks as much as 5 times. No significant differences are observed in the reported measurements between positively and negatively charged particles. Concerning the event selection criteria, the following are investigated: polarity of the magnetic field, reconstructed primary vertex position along the beam direction (selecting only events with z vtx within ±8 cm from the nominal interaction point), pile-up rejection (imposing stronger or weaker constraints on the consistency of different event multiplicity estimators) and variations in the instantaneous luminosity delivered to the ALICE detector by the LHC. The uncertainty on centrality determination is evaluated using an alternative  estimator based on the number of hits in the second ITS layer (|η| < 1.4), which is directly proportional to the number of charged particles in the corresponding acceptance. Among the aformentioned sources, for all observables in this article, track quality and centrality determination are the dominant sources. The total systematic uncertainties are evaluated summing in quadrature the systematic uncertainties coming from each of the sources, i.e. considering the different sources to be uncorrelated. Figure 1 presents the centrality dependence of the flow coefficients v n (n = 2, . . . 6) averaged in the p T range 0.2 < p T < 3.0 GeV/c, where collective effects are expected to dominate azimuthal anisotropies. Measurements are performed with the two-and four-particle cumulant method, denoted with v n {2, |∆η| > 1} and v n {4}, respectively. Results at both √ s NN = 5.02 and 2.76 TeV are shown. A clear hierarchy is observed among flow coefficients, with the second harmonic (elliptic flow) being the dominant one and the higher harmonics progressively smaller. The centrality-averaged (0-50%) values of harmonics v 3 − v 6 are decreasing as v n+1 /v n ∼ 0.5. In contrast to a strong increase in v 2 from central to mid-central collisions and a decrease after about 45% centrality towards peripheral collisions, a weak centrality dependence is observed for the higher harmonics. This holds true at both energies and is consistent with previous observations [24,47]. The characteristic centrality dependence of the elliptic flow was observed already at RHIC energies [48]. Compared to previous measurements in the p T range 0.2 < p T < 5 GeV/c [47], the differences in v n coefficients arising from the different choice of p T range are of about 1% and 2% for v 2 and v 3 − v 6 , respectively.  Figure 2 shows the ratio of v n {2, |∆η| > 1} (n = 2, 3, 4) and v 2 {4} between √ s NN = 5.02 and 2.76 TeV, i.e. the relative variation of these flow coefficients between those two energies. Since the systematic uncertainties of the measurements at different energies are partially correlated, the resulting systematic uncertainty on the ratio is reduced. All harmonics are observed to increase with energy, between about 2 and 10%. A hint of a centrality dependence is observed only for v 2 , with the increase growing slightly from mid-central towards more peripheral collisions. No significant difference is observed in the increase of v 2 measured with two-or four-particle correlations. Since the difference between v 2 {2, |∆η| > 1} and v 2 {4} is directly related to flow fluctuations, this observation suggests that the fluctuations of elliptic flow do not vary significantly between the two energies, within experimental uncertainties. The ratios are compared to hydrodynamical calculations with EKRT initial conditions [49] and different parametrisations of the temperature dependence of η/s [18]. The p-values for the comparison between data and models are also shown in Tab. 1. Among the two parametrisations that provide the best description of RHIC and LHC data [50], both are consistent with the measurements, except for v 3 {2, |∆η| > 1}, albeit the one with constant η/s = 0.2 agrees slightly better. These comparisons take into account the correlation between systematic uncertainties of data points in different centrality intervals. This observation might indicate little or no temperature dependence of η/s within the temperature range at which anisotropic flow develops at the two center of mass energies. As a reference, the p-values for the comparison between data and unity in the same centrality range (5-50%) are also reported in Tab. 1. Figure 3 shows the p T -differential measurements of v n (n = 2, . . . 6), with two-and four-particle cumulants, in the p T range 0.2 < p T < 50 GeV/c and in wide centrality bins, between 0% and 70%. The

Collision energy, transverse momentum and centrality dependence
0.484 0.468 0.022 Table 1: p-values for the comparison among ratios of v n {2, |∆η| > 1} (n = 2, 3, 4) and v 2 {4} between √ s NN = 5.02 and 2.76 TeV and model calculations using different parametrisations of η/s(T ) [18], shown in Fig. 2, and unity, in the centrality range 5-50%.  p T dependence is qualitatively similar for all harmonics: v n increases with increasing p T up to about 3-4 GeV/c, after which it starts decreasing. Comparing the different harmonics, they seem to follow the hierarchy observed in the p T -integrated results in the whole p T range: In the centrality range 10-40% a significant non-zero value of v 2 {2, |∆η| > 1} and v 2 {4} is observed up to p T ≈ 30 GeV/c; for the higher harmonics, significant values are only measured for p T ≤ 10 GeV/c.
Looking at the p T dependence in more detail, the flow harmonics are found to follow an approximate power-law scaling up to around the maximum, with exponents being proportional to the harmonic number n, v n (p T ) ∼ p n/3 T , as shown by the dashed fitted lines in Fig. 3. In ideal hydrodynamics, the p T dependence of anisotropic flow for massive particles should follow a power-law function v n (p T ) ∼ p n T in the region of p T /M up to order one, where M is the particle's mass, and at higher momenta it has been  predicted to be linear in p T for all n, v n (p T ) ∼ p T [51,52]. This p T dependence is notably different from the one observed in the data. At very low p T this is presumably because the relevant momentum region for inclusive particles, mostly pions, is below the range of our measurements, and at higher p T ideal hydro is not expected to hold because of momentum dependent viscous corrections at freeze out [53] and/or non-linear mode mixing for n 4 [20,54]. The emergence of this simple power-law dependence however remains unexpected and surprising. Figure 4 shows the p T -differential measurements of v n (n = 2, 3, 4) calculated with the scalar product method with respect to V0A. The same p T and centrality range as in Fig. 3 is shown. A significant v 2 {2, |∆η| > 2} is observed up to p T ≈ 40 GeV/c in the centrality range 10-40%. v n {2, |∆η| > 2} and v n {2, |∆η| > 1} (n = 2, 3, 4) are found to be compatible within 2% in the p T range 0.2 < p T < 10 GeV/c, while a systematic difference (with v 2 {2, |∆η| > 1} > v 2 {2, |∆η| > 2}) is observed for 10 < p T < 50 GeV/c, ranging from about 3% in centrality 0-5% to about 10% in centrality 40-50%. This difference is attributed to small residual non-flow contributions which are suppressed by the larger pseudorapidity gap. Such an interpretation is consistent with the observed centrality dependence: the relative contribution of non-flow increases with decreasing multiplicities, therefore from central to peripheral collisions. Possible differences among v n {2, |∆η| > 1} and v n {2, |∆η| > 2} (n = 2, 3, 4) arising from the decorrelation of event planes at different pseudorapidities have been estimated to be less than 1% and 3% for v 2 and v 3−4 , respectively, based on η-dependent factorization ratios [55] measured at 2.76 TeV [56]. This estimation assumes that such decorrelation only depends on |∆η| and not η in the pseudorapidity range under consideration (|η| < 5.1). Figure 5 shows the ratios of p T -differential v n {2, |∆η| > 1} (n = 2, 3, 4) and v 2 {4} between √ s NN = 5.02 and 2.76 TeV. Overall, the ratios are consistent with unity, indicating that p T -differential anisotropic flow does not change significantly across collision energies and that the increase observed in the p T -integrated values can be mostly attributed to an increase of p T , as previously noted [47]. This observation is also consistent with little or no variation of η/s between the two collision energies, as already shown in Fig.  2. The possible variations in p T -integrated values arising from the differences in the p T -differential ones have been estimated to be less than 1%, by integrating v n (p T ) at √ s NN = 5.02 (2.76) TeV with charged particle spectra at 2.76 (5.02) TeV. Figure 6 shows the comparison of p T -differential flow measurements with different models, in three centrality intervals: 5-10% (top panel), 20-30% (middle panel) and 40-50% (bottom panel). At low p T (p T < 2 GeV/c), flow coefficients are expected to be mostly determined by the collective expansion of the system, which is commonly described by hydrodynamic models. The measurements are compared to three calculations, one employing IP-Glasma initial conditions [57] matched to the MUSIC viscous hydrodynamic code [58] and two calculations using iEBE-VISHNU viscous hydrodynamic code [59] with AMPT [60] or TRENTo [61] initial conditions. The parameters of TRENTo were tuned to reproduce previous measurements in Pb-Pb collisions at √ s NN = 2.76 TeV [16]; with such tuning TRENTo has been shown [61] to effectively mimic IP-Glasma initial conditions and therefore the two calculations TRENTo+iEBE-VISHNU and IP-Glasma+MUSIC are expected to be based on similar initial conditions. All models employ a transport cascade model (UrQMD [62]) to describe the hadronic phase after freezeout. Compared to data, all models are found to underestimate the data for p T < 0.5 GeV/c. For 1 < p T < 2 GeV/c the predictions from IP-Glasma+MUSIC and TRENTo+iEBE-VISHNU overestimate the data, while those from AMPT-IC+iEBE-VISHNU are found to be still in agreement. Overall, all models qualitatively describe the p T dependence of flow coefficients in this low p T range.
At high p T (p T > 10 GeV/c), azimuthal anisotropies are on the contrary expected to be determined by path-length dependent parton energy-loss. The measurements are compared to predictions from [63], which combine an event-by-event hydrodynamic description of the medium (v-USPhydro [64]) with a jet quenching model (BBMG [65]). Two sets of predictions for v 2 {2}, v 2 {4} and v 3 {2}, assuming a linear (dE/dx ∼ L) and a quadratic (dE/dx ∼ L 2 ) dependence of the energy loss on the path length L, are compared to data. Other parameters of the model, such as η/s, are expected to have a minor contribution within the presented centrality ranges [63]. For v 2 {2, |∆η| > 2}, the linear case is compatible with the data, while the quadratic one can be excluded within 95% confidence level. For v 3 {2, |∆η| > 2}, neither of the two sets of predictions can be excluded within uncertainties. Our results are found to be in good agreement with CMS data [38].
The evolution of the shape of p T -differential v n coefficients with respect to centrality is investigated by calculating the ratios of v n (p T ) in a given centrality range and v n (p T ) in centrality 20-30%, normalised by the corresponding ratio of p T -integrated v n v n (p T ) ratio to 20 In order to reduce statistical fluctuations, a parametrisation of v n (p T )[20-30%] fitted to data is employed. If the shape of v n (p T ) does not change with centrality, v n (p T ) ratio to 20-30% is identical to 1 in the full p T range. The results are shown in Fig. 7: deviations from unity up to about 10% are observed at low p T (p T < 3 GeV/c) and up to about 30% at intermediate p T (3 < p T < 6 GeV/c), where v n (p T ) reaches its maximum. These variations are observed to be larger for higher harmonics (v 3−4 ), in particular for central collisions. The effects due to a change in particle composition of the inclusive charged particle sample with centrality are estimated to be negligible. These deviations are attributed mostly to the combined effect of radial flow and parton density which, in the coalescence model picture [66], decrease from central to peripheral collision shifting the maximum of v n (p T ) from higher to lower p T . At high p T (p T > 10 GeV/c), results on v 2 {2, |∆η| > 2} are consistent with those at low p T , suggesting a common origin of the centrality evolution of elliptic flow in the two regimes, presumably initial-state geometry and its fluctuations. This interpretation is consistent with the findings of [38]. The attribution of the scaling of v n (p T ) up to p T = 8 GeV/c to initial-state geometries agrees with studies [67,68] using the Event Shape Engineering technique [69] and p T -dependent elliptic flow fluctuations [70]. Finally, the models using hydrodynamic calculations [59] and jet energy loss ones [63] are observed to be in good agreement with the v 2 data at low and high p T , respectively.
At RHIC [71,72] and LHC [7] it had been observed that the ratios of harmonics follow a power-law scaling, i.e. v  Figure 8 shows these ratios for n = 3, 4 and m = 2, 3, as a function of p T . These ratios are indeed observed to be independent of p T , in most of the p T range and for most centrality ranges, except for centrality 0-5%. Up to about the maximum of v n (p T ), the scaling is numerically related to, but actually significantly more precise than, the observed approximate power-law dependences v n (p T ) ∼ p n/3 T pointed out in Fig. 3. Surprisingly however, the scaling extends much further, in particular v 3 (p T )/v 2 (p T ) 3/2 is constant to better than about 10%, out to the highest measured p T in excess of 10 GeV/c. The ratio v 4 (p T )/v 2 (p T ) 4/2 shows stronger deviations at high p T , starting at around the maximum of v 2 (p T ). A separation of v 4 into linear and non-linear components would be required to see if the v 4 /v 2 scaling at low p T , and/or its violations at high p T , is related to the mode mixing, which is particularly strong for the 4th harmonic and at high p T , or possibly also to quark coalescence [51,73,74].
As noted in the context of Fig. 3, the observed ratio scaling is not expected in ideal hydrodynamics. While not all viscous hydrodynamical models shown in Fig. 6 describe the data up to the highest p T very well, they all do exhibit the same power-law scaling in the ratio of harmonics over the p T range 0.5 < p T < 3 GeV/c, with a precision comparable to the one seen in the data, while they strongly deviate for p T < 0.5 GeV/c. The scaling may be related to viscosity, as also postulated in [75,76], in particular to the large and p T -dependent viscous corrections appearing at hadronisation [53]. However, a harmonic number dependence of these viscous corrections which could reproduce the scaling observed in the data, has so far, to the best of our knowledge, never been quantitatively investigated.    Fig. 7: Ratios v n (p T ) ratio to 20-30% of inclusive charged particles for Pb-Pb collisions at √ s NN = 5.02 TeV, in different centrality classes, measured with the scalar product method with respect to the V0A Q-vector. Hydrodynamic calculations [59,63] are shown for comparison. Figure 9 shows the integrated v 2 in the p T range 0.2 < p T < 3 GeV/c as a function of centrality, measured with two-, four-, six-and eight-particle cumulants at √ s NN = 5.02 and 2.76 TeV. The corresponding cumulants (c 2 {2, 4, 6, 8}) are reported in Fig. 10. The centrality dependence is similar for all multiparticle cumulants and similar to what is shown in Fig. 1. The differences between v 2 {2} (shown in Fig.  9) and v 2 {2, |∆η| > 1} (shown in Fig. 1), which range from about 4% in mid-central collisions to about 20% in peripheral ones, are mostly attributed to non-flow contributions, which are suppressed in the case of results with a pseudorapidity gap. The possible differences arising from the decorrelation of event planes at different pseudorapidities are expected to be less than 1%, as previously argued.

Elliptic flow fluctuations
A fine-splitting of less than 1% is observed among v 2 {4}, v 2 {6} and v 2 {8}, as it can be seen from their ratios, shown in Fig. 11 for both collision energies.
TeV show a significant centrality dependence: the deviations of the ratios from unity is about 0.2% in central and increases up to about 1% for mid-central collisions. A further increase seems to be observed for more peripheral collisions, up to about 2% for centralities above 50%. The systematic uncertainties on these ratios are greately reduced with respect to those on v 2 {m} (m = 2, 4, 6, 8), since the dominant sources of systematic uncertainty (track quality variables and centrality determination) among the two-and multi-particle cumulants are highly correlated. This fine-splitting is consistent with non-Bessel-Gaussian behaviour of event-by-event flow fluctuations, as previously explained. These ratios are found to be independent of the choice of p T range within 0.2 < p T < 3 GeV/c, indicating that the characterization of flow fluctuations at low p T does not depend on p T , even for such fine-splitting. Results at √ s NN = 2.76 TeV are found to be compatible, indicating that these ratios do not change significantly across collision energies. Compared to calculations [30] employing MC-Glauber initial conditions [77] and viscous hydrodynamics (v-USPhydro) for Pb-Pb collisions at √ s NN = 2.76 TeV, the  [29], as shown in Fig. 12. Figure 13 shows the ratio between v 2 {8} and v 2 {6} at √ s NN = 5.02 TeV. A hint of a further fine-splitting between these two, of the order of 0.05%, is observed. The results suggest little or no centrality dependence within centrality 10-50%. This difference is also consistent with non-Bessel-Gaussian elliptic flow fluctuations, and can be attributed to different contributions of the skewness to these higher-order cumulants [30]. Corresponding results at √ s NN = 2.76 TeV, here and in the following, are not shown because of the large statistical uncertainties. Figure 14 shows s NN = 5.02 TeV: these two are observed to be in agreement, which demonstrates the validity of Eq. 8. This observation sets an upper limit of 4 × 10 −4 at 95% confidence level for possible contributions to multi-particle cumulants from higher moments of the flow p.d.f. (kurtosis and beyond) in the centrality range 10-50%. This estimate is obtained assuming gaussian systematic uncertainties and summing them in quadrature with the statistical ones.  [30] for Pb-Pb collisions at √ s NN = 2.76 TeV, the results are found to be compatible for the entire centrality range. This observation is consistent with the elliptic flow p.d.f. being progressively more left-skewed going from central to peripheral collisions. We attribute this feature to the combination of an increase in ε 2 and the geometrical constrain ε 2 < 1, as previously argued.
In order to report the full p.d.f. of elliptic flow P(v 2 ), which can be compared to previous experimental  results and theoretical predictions, it is parametrised with the Elliptic Power distribution [26,27] and its three free parameters (α, ε 0 and k 2 ) are extracted from fits to the elliptic flow cumulants c 2 {2, |∆η| > 1} and c 2 {m} (m = 4, 6, 8) at √ s NN = 5.02 TeV. The parameter α quantifies the magnitude of elliptic flow fluctuations, ε 0 the mean eccentricity in the reaction plane and k 2 is the proportionality coefficient between initial-state eccentricity and v 2 coefficient: v 2 = k 2 ε 2 . The relation between cumulants and Elliptic Power parameters is given by [27] where and 2 F 1 is the hypergeometric function. The results are shown in Fig. 16. The systematic uncertainties are assigned varying the fit ranges and initial values of the parameters and shifting the data points according to the corresponding systematic uncertainties. An additional source of uncertainty, which is investigated, is a possible cubic response coefficient k 2 , defined as v 2 = k 2 ε 2 + k 2 ε 3 2 . This coefficient is introduced to quantify the possible increase of flow fluctuations that the hydrodynamic expansion of the medium Pb-Pb 2.76 TeV Fig. 12: Ratios of elliptic flow coefficients v 2 of inclusive charged particles between measurements with different multi-particle cumulant methods, as a function of centrality, at √ s NN = 2.76 TeV. Hydrodynamic calculations [30] and ATLAS measurements [29] are shown for comparison. introduces with respect to geometrical fluctuations in the initial state and was argued to be non-zero in mid-central and peripheral collisions due to general properties of the hydrodynamic phase [78]. In particular, k 2 is expected to be ≤ 0.15 in the centrality range 0-60% [78]. The residual differences in α, ε 0 and k 2 when including k 2 as an additional free parameter are considered in the systematic uncertainties. The statistical uncertainties are evaluated using the subsampling method: the analysed dataset is divided into 10 sub-samples and v 2 {m} is measured in each of them. The Elliptic Power parameters are then extracted in each subsample and their dispersion is used to estimate the statistical uncertainties.

Centrality (%)
The resulting p.d.f., constructed using the Elliptic Power distribution (Eq. 17) with the parameters shown in Fig. 16 and scaled by its mean ( v 2 ), is reported in Fig. 17, for centralities 5-10%, 25-30% and 45-50%. The systematic uncertainties take into account the correlation of the uncertainties of the Elliptic Power parameters. Other centrality ranges that are not shown here are reported in the appendix A. Scaling by v 2 allows a comparison of our data with results by the ATLAS collaboration [67] obtained in different p T ranges. The observed agreement is also consistent, as previously noted, with elliptic flow fluctuations at low p T not depending on p T and not changing significantly between collision energies, except for the trivial increase in p T -integrated v 2 due to the change in p T . Comparison with iEBE-VISHNU model calculations with AMPT and TRENTo initial conditions [59] indicates that TRENTo initial conditions are better at describing the experimental data. The data are found to be in agreement also with predictions from the IP-Glasma+MUSIC model [58] (with initial conditions very similar to the TRENTo ones), although the uncertainties on the theoretical predictions do not allow to draw firm conclusions.

Conclusions
Anisotropic flow coefficients are measured up to the sixth harmonic for inclusive charged particles at mid-rapidity (|η| < 0.8), in a wide centrality (0-80%) and p T (0.2 < p T < 50 GeV/c) ranges, for Pb-Pb collisions at √ s NN = 5.02 and 2.76 TeV. Comparing the results at √ s NN = 5.02 and 2.76 TeV the energy dependence of anisotropic flow at the LHC is investigated. Comparison with different model calculations demonstrates that these measurements have the potential to constrain initial-state fluctuations, transport parameters of the medium and path-length dependence of energy loss of high-p T partons. The evolution of v n (p T ) with respect to centrality and harmonic number n is also investigated. Flow coefficient of all harmonics are observed to follow an approximate power-law scaling of the form v n (p T ) ∼ p n/3 T in the p T range 0.2 < p T < 3 GeV/c. The ratios v n /v n/m m n = 3, 4 and m = 2, 3 are also observed to be independent of p T within the same p T range and show deviations of about 10% for 3 < p T < 10 GeV/c.   [58,59] and previous measurements from ATLAS [67] at lower energies are shown for comparison.
The fluctuations of elliptic flow are investigated through the fine-splitting of the higher-order multiparticle cumulants (v 2 {4}, v 2 {6}, v 2 {8}), from which the standardised skewness (γ exp 1 ) of the flow p.d.f. is extracted. Results are found to be compatible both with predictions from hydrodynamical models and with previous ATLAS results at lower energies. It is concluded that the characterization of elliptic flow fluctuations at low p T does not depend on the p T range and on the collision energy, except for the increase in p T -integrated v 2 due to the change in p T . Direct constraints on the contribution of higher moments to the multi-particle cumulants are also reported. Finally, the full elliptic flow p.d.f., parametrised with the Elliptic Power distribution, is reported in the centrality ranges 0-60%. These results are also found to be in agreement with previous experimental results. Overall, calculations including initial conditions matching the IP-Glasma description are observed to better reproduce the elliptic flow p.d.f. while failing to describe the p T dependence of anisotropic flow coefficients, whereas the opposite situation is observed for calculations that employ AMPT initial conditions. [35] E. V. Shuryak, "The Azimuthal asymmetry at large p(t) seem to be too large for a 'jet quenching'," Phys. Rev. C66 (2002) 027902, arXiv:nucl-th/0112042 [nucl-th].

A Appendix
The elliptic flow p.d.f. P(v 2 ), constructed as explained in Sec. 4, in the centrality ranges not shown in Fig. 17 are reported in Fig. A.1-A.3.