Measurement of flow harmonics with multi-particle cumulants in Pb+Pb collisions at \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}}}=2.76$$\end{document}sNN=2.76 TeV with the ATLAS detector

ATLAS measurements of the azimuthal anisotropy in lead–lead collisions at \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}}}=2.76$$\end{document}sNN=2.76 TeV are shown using a dataset of approximately 7 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upmu $$\end{document}μb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{-1}$$\end{document}-1 collected at the LHC in 2010. The measurements are performed for charged particles with transverse momenta \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.5<p_{\mathrm {T}}<20$$\end{document}0.5<pT<20 GeV and in the pseudorapidity range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\eta |<2.5$$\end{document}|η|<2.5. The anisotropy is characterized by the Fourier coefficients, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {v}_n$$\end{document}vn, of the charged-particle azimuthal angle distribution for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n = 2$$\end{document}n=2–4. The Fourier coefficients are evaluated using multi-particle cumulants calculated with the generating function method. Results on the transverse momentum, pseudorapidity and centrality dependence of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {v}_n$$\end{document}vn coefficients are presented. The elliptic flow, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {v}_2$$\end{document}v2, is obtained from the two-, four-, six- and eight-particle cumulants while higher-order coefficients, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {v}_3$$\end{document}v3 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {v}_4$$\end{document}v4, are determined with two- and four-particle cumulants. Flow harmonics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {v}_n$$\end{document}vn measured with four-particle cumulants are significantly reduced compared to the measurement involving two-particle cumulants. A comparison to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {v}_n$$\end{document}vn measurements obtained using different analysis methods and previously reported by the LHC experiments is also shown. Results of measurements of flow fluctuations evaluated with multi-particle cumulants are shown as a function of transverse momentum and the collision centrality. Models of the initial spatial geometry and its fluctuations fail to describe the flow fluctuations measurements.


Introduction
The anisotropy of charged-particle azimuthal angle distributions in heavy-ion collisions has been a subject of extensive experimental studies at RHIC [1][2][3][4][5][6] and more recently at the LHC [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]. The results provide conclusive evidence that the hot and dense matter produced in these collisions behaves collectively and has properties resembling those of a nearly perfect fluid [25]. e-mail: atlas.publications@cern.ch The final-state anisotropy is a global property of particle production that arises from the initial spatial asymmetry of the collision region in a plane transverse to the beam axis for heavy-ion collisions with a non-zero impact parameter. It is characterized by the coefficients, v n , of the Fourier expansion of the measured azimuthal angle distributions [1,26]: where n is the order of the Fourier harmonic, referred to as flow harmonic, φ is the azimuthal angle of the outgoing particle, n defines the azimuthal angle of the nth-order symmetry plane of the initial geometry [15], and the angled brackets denote an average over charged particles in an event.
Due to the symmetry in the azimuth of the plane defined by n , all sine terms of the Fourier expansion vanish. For evaluation of the coefficients v n in the "event-plane" method, the initial plane of symmetry is estimated from the measured correlations between particles, using the so-called subevent method [26]. As a consequence, only the two-particle correlations are exploited in the determination of v n (see Eq. 1). This leads to a problem of disentangling all-particle flow and contributions from particle correlations unrelated to the initial geometry, known as non-flow correlations. These non-flow effects include correlations due to energy and momentum conservation, resonance decays, quantum interference phenomena and jet production. They generally involve only a small number of produced particles. In order to suppress non-flow correlations, methods that use genuine multi-particle correlations, estimated using cumulants, were proposed [27][28][29][30].
Calculating multi-particle correlations in large-multiplicity heavy-ion collisions at high energies is limited by the computing requirements needed to perform nested loops over thousand of particles per event to analyse all particle multiplets. To avoid this problem, the generating function formalism [27][28][29] is exploited to calculate multi-particle cumulants, and the results obtained are presented in this paper. An alternative approach was proposed in Ref. [30] to express multi-particle correlations in terms of the moments of the flow vector, Q n , and is used in this paper as a cross-check of multi-particle cumulants obtained with the generating function method. The cumulant approach to measure flow harmonics also provides the possibility to study event-to-event fluctuations in the amplitudes of different harmonics, which can be related to the fluctuations in the initial transverse shape of the interaction region [31][32][33].
The cumulant method has been used to measure the anisotropic flow in NA49 [34], STAR [35] and recently also at the LHC experiments [7,9,20,23]. The results show that the Fourier coefficients determined with four-particle cumulants are smaller than those derived with two-particle cumulants due to the suppression in the former of non-flow two-particle correlations. In this paper, the method is used to measure flow harmonics in lead-lead collisions at √ s NN = 2.76 TeV with the ATLAS detector. The elliptic flow v 2 is measured using two-, four-, six-and eight-particle cumulants. For v 3 and v 4 measurements the two-and four-particle cumulants are exploited.
This paper is organized as follows. Section 2 describes the ATLAS detector, trigger, and offline event selections. Section 3 contains a description of additional selection criteria for events and charged-particle tracks. Section 4 gives details of the Monte Carlo simulation samples used to derive the tracking efficiency and fake-track rates. The analysis method and procedure is outlined in Sect. 5. Section 6 contains a discussion of the systematic errors. Results are presented in Sect. 7. Section 8 is devoted to summary and conclusions.

The ATLAS detector and trigger
The results presented in this paper were obtained from a sample of minimum-bias lead-lead collisions at √ s NN = 2.76 TeV recorded by ATLAS in 2010 and corresponding to an integrated luminosity of approximately 7 µb −1 . The measurements were performed using the ATLAS inner detector and forward calorimeters [36]. The inner detector covers the complete azimuthal range and extends over the pseudorapidity region |η| < 2.5. 1 The inner detector silicon tracker, used in this analysis for track reconstruction, consists of layers of pixel and microstrip detectors (SCT) immersed in a 2 T axial magnetic field. The forward calorimeters (FCal) use liquid argon with copper-tungsten absorbers to perform both the electromagnetic and hadronic energy measurements 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upward. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the beam pipe. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). with copper-tungsten/liquid argon technology, and also provide complete coverage in azimuth for 3.2 < |η| < 4.9. The trigger system was used to select minimum-bias lead-lead collisions. It required a coincidence of signals recorded in both zero-degree calorimeters (ZDC), located symmetrically at z = ±140 m, and in the minimum-bias trigger scintillator (MBTS) counters at z = ±3.6 m.

Event and track selections
Additional offline event selections were also applied, requiring a time difference between the two MBTS counters of less than 3 ns and at least one primary vertex reconstructed using charged-particle tracks. Events satisfying the abovedescribed selections were also required to have a reconstructed primary vertex within 100 mm of the nominal centre of the ATLAS detector.
The precision silicon tracking detectors were used to reconstruct charged-particle trajectories with a minimum p T of 0.5 GeV. Special track-quality criteria are imposed to deal with high particle densities in Pb+Pb collisions. Tracks are required to have at least eight hits in the SCT, at least two pixel hits and a hit in the pixel layer closest to the interaction point if expected. A track must have no missing pixel hits and no missing SCT hits, when a hit is expected. The transverse and longitudinal impact parameters with respect to the vertex, |d 0 | and |z 0 sin θ | respectively, were each required to be less than 1 mm. Specifically for this analysis it was also required that |d 0 /σ d 0 | < 3 and |z 0 sin θ/σ z | < 3, where σ d 0 and σ z are the uncertainties on d 0 and z 0 sin θ , respectively, as obtained from the covariance matrix of the track fit. The latter requirements improve both the tracking performance at high p T and the purity of the track sample. The number of reconstructed tracks per event is denoted N rec ch . For this analysis, the additional requirement of N rec ch ≥ 10 for tracks with 0.5 < p T < 5 GeV was imposed to allow the measurement of correlations involving as many as eight particles.
The correlation between the summed transverse energy ( E FCal T ) measured in the FCal and N rec ch was investigated in order to identify background events. Events having an N rec ch vs. E FCal T correlation distinctly different from that for the majority of Pb+Pb collisions were removed. The removed events, less than 0.01 % of the sample, were found to contain multiple Pb+Pb collisions. After applying all selection requirements, the data sample consists of about 35 × 10 6 Pb+Pb collision events.
The summed transverse energy is used to define the centrality of the collision. A detailed analysis of the E FCal T distribution [15] showed that the fraction of the total inelastic cross-section sampled by the trigger and event selection requirements is (98 ± 2) %. The E FCal centile fraction of all events after accounting for the 2 % inefficiency in recording the most peripheral collisions. The analysis is performed in narrow centrality intervals: 1 % centrality bins for the 20 % of events with the largest E FCal T , and 5 % centrality bins for the remaining events. These narrow centrality intervals are then combined into wider bins to ensure sufficiently small statistical uncertainties on the measured flow harmonics. The 20 % of events with the smallest E FCal T (most peripheral collisions) are not considered in this analysis, due to the inefficiency in the event triggering and the correspondingly large uncertainties of measurements performed for these low-multiplicity collisions. For each centrality interval, a standard Glauber Monte Carlo model [37,38] is used to estimate the average number of participating nucleons, N part , which provides an alternative measure of the collision centrality.

Monte Carlo simulations
A Monte Carlo (MC) sample was used in the analysis to determine tracking efficiencies and rates of falsely reconstructed tracks (fake-track rates). The HIJING event generator [39] was used to produce minimum-bias Pb+Pb collisions. Events were generated with the default parameters, except for the jet quenching, which was turned off. Flow harmonics were introduced into HIJING at the generator level by changing the azimuthal angle of each particle [1] in order to produce an anisotropic azimuthal angle distribution consistent with previous ATLAS v n (n = 2-6) measurements [15,16]. The detector response simulation [40] uses the GEANT4 package [41] with data-taking conditions corresponding to those of the 2010 Pb+Pb run and simulated events are reconstructed in the same way as the data.
The tracking efficiency, ( p T , η), and the fake-track rate f ( p T , η) are determined [42] using the Monte Carlo sample described above. The MC reproduces the measured centrality dependence of the track-quality parameters. The efficiency is found to depend weakly on the collision centrality. For the lowest transverse momenta (0.5-0.6 GeV), the efficiency at |η| < 1 is of the order of 50 % and falls to about 30 % at high |η|. For higher transverse momenta it reaches about 70 % at |η| < 1 and drops to about 50 % at high |η|. The rate of falsely reconstructed tracks (the fake-track rate) is typically below 1 %. It increases to 3-7 % for the lowest transverse momenta in the most central collisions.

Analysis procedure
Fourier coefficients, v n , are measured using 2k-particle correlations [27][28][29] defined as: where the notation v n {2k} is used for the v n flow harmonic derived from the 2k-particle correlations, and k is an integer. Azimuthal angles of particles forming a 2k-particle cluster are denoted by φ l , where l = 1, . . . , 2k. The double angled brackets denote an average, first over charged particles in an event, and then over events, while the single angled brackets denote averaging over events. The multi-particle correlation, corr n {2k} , includes contributions from the collective anisotropic flow and from non-flow effects (see Sect. 1). It was proposed in Refs. [27][28][29] to exploit the cumulant expansion of multi-particle correlations in order to reduce the nonflow contribution. The anisotropic flow related to the initial geometry is a global, collective effect involving correlations between all outgoing particles. Thus, in the absence of nonflow effects, v n {2k} is expected to be independent of k. On the other hand, most of the non-flow correlations, such as resonance decays or interference effects, contribute only to correlations between small numbers of particles. The idea of using 2k-particle cumulants is to suppress the non-flow contribution by eliminating the correlations which act between fewer than 2k particles. More specifically, the cumulant of e.g. the four-particle correlations, defined as: measures the genuine four-particle correlations. So, if the non-flow contribution is only due to the two-particle correlations, then c n {4} directly measures flow harmonics. Similarly, using the cumulant of the six-particle correlations allows one to remove contributions from two-and four-particle correlations. The different cumulants provide independent estimates of the same flow harmonic v n , with the estimate based on correlations among many particles being more precise due to the suppressed non-flow correlations. In the absence of non-flow correlations, cumulants of different order should give the same estimate of v n . The generating function formalism for calculating 2kparticle cumulants (GFC method) was proposed in Ref. [29]. With this method, the number of required computing operations is proportional to the number of particles per event. The cumulant generating function of multi-particle azimuthal correlations, C n (z), is defined in the plane of a complex variable z as: where the angled brackets represent the average over events in a given centrality interval, and the product runs over the N particles within a given Pb+Pb event [27][28][29]. The weighting factors, w j , are used in this analysis to correct for any non-uniformity in the azimuthal angle distribution of reconstructed tracks. The weights are obtained from the data using the two-dimensional distribution in the η-φ plane of all reconstructed tracks. For each bin j in (δη, δφ) = (0.1, 2π/64) a weight is calculated as w j = N (δη) /N (δη, δφ), where N (δη) is the average number of tracks in the δη slice to which this bin belongs, while N (δη, δφ) is the number of tracks in the (δη, δφ) bin.
The expansion of the cumulant generating function in powers of |z| provides the cumulant c n {2k}, which is equal to the coefficient of the term |z| 2k /k! 2 of this expansion. In practice, to construct the c n {2k} cumulant the power series is truncated to order |z| 2k and C n (z) is computed at a discrete set of interpolating points z p,q = x p,q + iy p,q [29], where: For this analysis, the parameters p = 1, . . . , 5 and q = 0, . . . , q max − 1 with q max = 11 were chosen as recommended in Ref. [29]. The r 0 parameter (r 0 ≡ |z|/ √ p) should be as small as possible, chosen such that the results remain stable under its variation. The r 0 values used were chosen to be 4.0, 2.2, 1.6, 1.1 and 1.0 for centrality intervals 0-5 %, 5-10 %, 10-20 %, 20-30 % and 30-80 %, respectively. For these values, the cumulants are found to be stable when varying r 0 between r 0 /2 and 2r 0 . The only differences, up to about 2 %, were seen when using the eight-particle cumulants to calculate the elliptic flow harmonic and are accounted for in the systematic uncertainty on v 2 {8}. An alternative method to calculate multi-particle correlations and cumulants in a single pass over all particles in each event, referred to as the QC method, was proposed in Ref. [30]. In this method, the expressions for the multiparticle correlations are derived in terms of the moments of the flow vector Q n , defined as Q n = N j=1 w j e inφ j , where the index n denotes the order of the flow harmonic, the sum runs over all N particles in an event and w j are weights as defined above. The QC method is used to calculate the cumulants, c n {2k}, which are compared with the cumulants obtained from the GFC method.
A practical application of the cumulant method involves two main steps [27][28][29]. First, the reference 2k-particle cumulants, c n {2k}, are derived from the cumulant generating function calculated from particles measured over a broad range of transverse momentum and pseudorapidity. This step is equivalent to the event-plane estimate in the standard method (see Eq. 1) and the reference cumulants play a similar role to the event-plane resolution correction [26]. In Fig. 1 Multi-particle cumulants for the second-order flow harmonic, c 2 {2k} for k = 1, 2, 3, 4, obtained with the GFC [29] and QC [30] methods shown as a function of centrality. The horizontal axis ranges from central collisions (0-20 %) to more peripheral collisions (60-80 %) the next step, the differential flow is calculated in p T and η bins using cumulants, denoted d n {2k}, computed from a differential generating function. To determine the d n {2k} cumulants, each charged particle from a p T and η bin is correlated with 2k − 1 reference particles. The differential flow harmonics v n {2k}( p T , η), are then calculated with respect to the reference cumulants as prescribed in Refs. [28,29]: In order to calculate the reference cumulants, c n {2k}, all charged particles with pseudorapidities |η| < 2.5 and transverse momenta 0.5 < p T < 5 GeV are used in this analysis. The results for c 2 are shown as a function of centrality in Fig. 1 for two-, four-, six-and eight-particle cumulants obtained from the GFC and QC methods. The figure shows that the two methods yield consistent results over a wide range of collision centralities. Differences, up to ∼ 20 %, are observed only for the most peripheral collisions. For the most central (0-2 %) Pb+Pb collisions, the cumulants c n {2k} for k > 1 are, within sizeable statistical errors, consistent with zero. However, they have incorrect signs, which prevents the calculation of flow harmonics due to the square-root function in the denominator of Eqs. (8), (9) and (10).
For higher flow harmonics, the cumulants c 3 {2k} and c 4 {2k} obtained from both the GFC and QC methods are consistent with zero for k > 2. Therefore, only two-and fourparticle cumulants can be used to derive third-and fourthorder flow coefficients. Figure 2 shows the centrality dependence of the two-and four-particle cumulants, obtained from the GFC and QC methods, for n = 3 and n = 4. The figure demonstrates an overall good agreement between the cumulants calculated using the two different methods. In the case of four-particle cumulants, the centrality range of the method's applicability is limited to 0-60 % for n = 3 and 0-25 % for n = 4.
The differential flow harmonics, v n {2k}( p T , η), are determined using the differential cumulants d n {2k} and Eqs. (7)-(10) in bins of transverse momentum and pseudorapidity for events from a given centrality interval. The pseudorapidity range |η| < 2.5 is divided into 50 bins of width 0.1 each. In transverse momentum, 28 bins of variable width, covering the p T range from 0.5 GeV to 20 GeV, are used. These differential flow harmonics can then be integrated over wider phase-space bins or the full range in either pseudorapidity or transverse momentum, or both. In this integration procedure, the harmonics v n {2k}( p T , η) measured in each small bin are weighted by the charged-particle multiplicity in that bin, corrected for tracking efficiency and fake-track rate, using the MC-determined corrections ( p T , η) and f ( p T , η) as described in Sect. 4.

Systematic uncertainties
The systematic uncertainties on the measurements presented in this paper are evaluated by varying different aspects of the analysis and comparing the results obtained to the baseline results for the transverse momentum, pseudorapidity and centrality dependence of v 2 , v 3 and v 4 . The following sources are considered as potential contributors to the systematic uncertainty on the measured flow harmonics.
An overall scale uncertainty on flow harmonics comes from the uncertainty in the fraction of the total inelastic crosssection accepted by the trigger and the event selection criteria. It is evaluated by varying the centrality bin definitions, using the modified selections on E FCal T , which account for the 2 % uncertainty in the sampled fraction of the crosssection.
All formulas applied in the analysis are valid under the assumption that the sine terms in the Fourier expansion vanish due to the azimuthal symmetry of the initial geometry. However, due to some distortions in the detector acceptance and non-uniformities in the measured azimuthal angle distributions, a residual sine term may be present. The magnitude of the sine term is calculated as the imaginary part of the differential generating function. The deviation from zero of the average sine term with respect to v n is treated as the systematic uncertainty. Some detector distortions may lead to an asymmetry between positive and negative η hemispheres, and the difference between the flow harmonics measured at positive and negative pseudorapidities is also considered as a systematic uncertainty.
A small contribution to the systematic uncertainty, only for v 2 {8}, comes from the stability of the results with respect to the assumed value of the r 0 parameter (see discussion in the previous section). The correction applied to ensure the uniformity of the azimuthal angle distribution of reconstructed tracks (via weights w j ) is also checked by comparing the baseline results to those obtained with w j ≡ 1. The contribution to the systematic uncertainty related to the trackquality definition is evaluated by comparing results obtained with more restrictive or less restrictive requirements. Both the transverse and longitudinal impact parameter cuts, |d 0 | and |z 0 sin θ |, are changed by ±0.5 mm with respect to the nominal value of 1 mm and the significance cuts, |d 0 /σ d 0 | < 3 and |z 0 sin θ/σ z | < 3, are changed by ±1.
The analysis procedure is also checked through MC studies by comparing the observables at the generator/particle level with those obtained in the MC simulated sample for which the same analysis chain and correction procedure is used as for the data. The measured flow harmonics in data agree qualitatively with the reconstructed MC harmonics and show similar trends as a function of η. In the phase-space region where tracking performance suffers from low efficiency and high fake-track rates ( p T < 1.5 GeV and |η| > 1), systematic differences are observed between the flow harmonics calculated at the generator level and at the reconstruction level after the corrections. In general, in this phase-space region, the reconstructed flow harmonics are smaller than the generator-level ones and show an η dependence, not present at the generator level. To account for this η-dependent bias,    3 and v 4 measured with four-particle cumulants averaged over the accessible centrality ranges 0-60 % and 0-25 %, respectively. A single entry is given where the uncertainty only varies by a small amount over the p T range from 0.5 to 20 GeV or η range from −2.5 to 2.5. Otherwise the range in uncertainties is provided corresponding to the range in p T or η For higher-order flow harmonics, the sine term and MC closure each contribute about 3 %, while all other sources contribute less than 1 %. The total systematic uncertainty is the sum in quadrature of all the individual contributions. For illustration, Table 1 shows the total systematic uncertainties on v 2 {4} for three representative centrality intervals: 2-5 %, 15-20 % and 60-80 %. For higher-order flow harmonics, the systematic uncertainties are listed in Table 2. For comparison, the statistical uncertainties on v n are also listed. It can be seen that the uncertainties on the measured v 2 at high p T and on v 3 and v 4 over the whole kinematic range, are dominated by large statistical errors.
In addition to the evaluation of the systematic uncertainty, further cross-checks are performed. The comparison between cumulants calculated with the GFC and QC methods is discussed in detail in Sect. 5. The analysis is performed separately for negatively and positively charged particles and the resulting v n coefficients are found to be consistent within their statistical and systematic uncertainties. Furthermore, the consistency of v 2 {2k} for k > 1 measured for same-sign particles and all combinations of charged particles confirms the global collective feature of the measured effect. Consistency is also observed between measurements obtained from sub-samples collected early and late during the data-taking period. The analysis also evaluates a potential bias which may be due to the large spread in charged-particle multiplicities in centrality intervals defined by E FCal T . For this purpose, in a given centrality bin selected by E FCal T , the analysis is restricted to events with multiplicities limited to a very narrow range (corresponding to ±RMS/2 around the mean multiplicity) and compared to the analysis performed for the full range of multiplicities. Both give v n harmonics consistent with each other within their statistical and systematic uncertainties.

Transverse momentum dependence of flow harmonics
The transverse momentum dependence of v 2 {2}, v 2 {4}, v 2 {6} and v 2 {8} is shown in Fig. 3 in 14 centrality intervals as indicated in the plots. The v 2 coefficients are integrated over the full pseudorapidity range |η| < 2.5. The elliptic flow mea-surements, v 2 {EP}, obtained with the event-plane method [26] are also shown. For this comparison, the measurements from Ref. [15] were reanalysed with the same track-quality requirements and centrality intervals as for the present analysis. The event-plane v 2 is systematically smaller than v 2 {2} since it is less affected by short-range two-particle correlations, which are partially removed from v 2 {EP} due to the separation between the phase-space region where the event plane angle is determined (3.2 < |η| < 4.9) and the phase space where charged-particle momenta are reconstructed (|η| < 2.5). This effect is particularly pronounced at high transverse momenta, where v 2 {2} is strongly influenced by jet-related correlations. At lower transverse momenta, the difference between v 2 {2} and v 2 {EP} can also be attributed to flow fluctuations. The difference between v 2 {EP} and v 2 {4} is mainly due to flow fluctuations. The v 2 {2k} for k > 1 is systematically smaller than v 2 {2}, consistent with the expected suppression of non-flow effects in v 2 obtained with cumulants of more than two particles. The results for v 2 {2k} with k > 1 agree with each other, within the uncertainties, for all centrality intervals, indicating that already the four-particle  and v 2 {8} at high transverse momenta may reflect both the anisotropy of the initial geometry and the path-length dependence of the parton energy loss in the dense, strongly interacting medium [45]. Figure 4 shows the comparison of our results for v 2 {4} integrated over |η| < 0.8 as a function of p T , to these coefficients measured by the CMS [20] and ALICE [9] experiments in several centrality intervals. The results on the elliptic flow harmonic measured with four-particle cumulants are consistent within uncertainties for the three experiments.
The transverse momentum dependence of the higher-order harmonics, v 3 and v 4 , is shown in Fig. 5 and compared to the results obtained with the event-plane method. Due to the large uncertainties on the harmonics measured with fourparticle cumulants, especially for events with low multiplicities, the results are shown in wide centrality ranges: for v 3 in the two broad centrality intervals, 0-25 % and 25-60 %, and for v 4 in the full accessible centrality range, 0-25 %. In addition, the results for v n {4} are shown in fine p T bins

Pseudorapidity dependence of flow harmonics
The pseudorapidity dependence of v n {2k} is studied as a function of centrality for flow coefficients integrated over the p T range from 0.5 GeV to 20 GeV. Figure 6 shows v 2 {2},  The centrality dependence of the elliptic flow harmonic, integrated over the full range in η and p T and obtained with cumulants of various orders, is shown as a function of N part in Fig. 8. The coefficients v 2 {2k}, and in general v n {2k}, can also be calculated from the moments of the distribution, p(v n ), of the event-by-event (EbyE) measured flow harmonics, v k n = v k n p(v n ) as: ATLAS has measured p(v n ) for n = 2, 3, 4 [17]. The comparison of v 2 {2k} obtained with the cumulant method to v calc 2 {2k, EbyE} is shown in Fig. 8. Good agreement between the two independent measurements is seen. The cumulant method gives v 2 values larger than those calculated from the p(v 2 ) distribution only for v 2 {2} measured in the most peripheral collisions, due to contributions from short-range two-particle correlations in the former. The ratios of v 2 {6} and v 2 {8} to v 2 {4} are shown in Fig. 9. The left panel shows results from the cumulant method. The ratios are systematically below unity, most significantly at low N part . This effect, which is of the order of 1-2 %, is significant for the ratio v 2 {6}/v 2 {4} while it is within the present uncertainty of the cumulant measurements for v 2 {8}/v 2 {4}. Better precision is achieved for v calc 2 {2k, EbyE} (right panel of Fig. 9). The difference between v 2 {4} and v 2 {6} or v 2 {8} is attributed to the non-Bessel-Gaussian character of the p(v 2 ) distribution measured in peripheral collisions [17].
It is interesting to compare flow harmonic measurements obtained with different methods, which have different sensi- The v n {2} values are the largest (predominantly) due to large contributions from short-range two-particle correlations, which are suppressed in the event-plane v n measurements. The v n coefficients measured with the event-plane method are systematically larger than the mean values of the event-by-event measurement of flow harmonics. This difference is naturally attributed to the flow fluctuations, which contribute to v n {EP} but are suppressed in v n {EbyE}. The flow coefficients measured with the four-particle cumulant method are the smallest, mainly due to the contribution from flow fluctuations, which is negative for v n {4} and positive for v n measured with the event-plane method. In addition, some residual two-particle correlations unrelated to the azimuthal asymmetry in the initial geometry contribute to v n {EP}, but are negligibly small in the case of v n {4}.
The centrality dependence of v n {4} is shown in Fig. 11 for n = 2, 3 and 4. The elliptic flow v 2 {4} shows a strong centrality dependence, rising with N part until reaching a maximum at N part ≈ 100, and then decreasing for more central collisions. This strong centrality dependence is not seen for the higher flow harmonics v 3 {4} and v 4 {4}. In addition, the magnitude of the third-and fourth-order flow coefficients is much smaller than the magnitude of the elliptic flow; e.g. for N part of about 300, v n ≈ 0.05, 0.02, and 0.01 for v 2 , v 3 and v 4 , respectively. For smaller N part this difference is even larger, with v 2 reaching more than 0.1 and v 3 and v 4 staying at the same level as at higher N part . Figure 11 also shows v calc 2,3 {4, EbyE} obtained from the measured p(v 2 ) and p(v 3 ) distributions. The measured p(v 4 ) in Ref. [17] is truncated at large values of v 4 and therefore is not used here for the comparison. Good agreement between the two independent measurements is also seen for the third-order flow harmonics.

Fluctuations of flow harmonics
Measurements of elliptic flow dynamic fluctuations have attracted much interest, since flow fluctuations can be traced back to fluctuations of the initial collision zone. Experimentally, flow fluctuations are difficult to measure due to unavoidable contamination by non-flow effects. The reported elliptic flow fluctuation measurements from RHIC [31][32][33] are affected by non-flow correlations, despite the attempts made to estimate their contribution. Interestingly, RHIC results indicate that flow fluctuations are mostly determined by initial-state geometry fluctuations, which thus seem to be preserved throughout the system evolution.
The relative flow harmonic fluctuations, defined as σ v n / v n , can be calculated using the width and mean value of the p(v n ) distributions and compared to the predictions for fluctuations of the initial geometry. The latter can be characterized by the eccentricities, ε n , which can be estimated from the transverse positions (r, φ) of nucleons participating in the collision: ε n = r n cos nφ 2 + r n sin nφ 2 r n .
Such a comparison of σ v n / v n , derived from the event-byevent measurement of v n , to the Glauber model [37] and MC-KLN model [46], which combines the Glauber approach with saturated low-x gluon distribution functions, is discussed in Ref. [17]. In general, none of the considered models of the relative fluctuations of ε n gives a consistent description of the relative flow fluctuations over the entire range of collision centralities.
In this analysis, the measure of relative flow fluctuations, F(v n ), is defined as: The above formula provides a valid estimate of σ v n / v n under the assumptions that non-flow correlations are absent in v n {2} and v n {4}, and that flow fluctuations are small compared to v n (σ v n << v n ). The first assumption is obviously not fulfilled by v n {2}, which is strongly contaminated by non-flow correlations. Therefore, v n {EP} is used instead of v n {2}, following the approach proposed in Ref. [9]. The second assumption is not valid for fluctuations of third-and fourth-order flow harmonics, and also for elliptic flow harmonics measured in the most central Pb+Pb collisions. Nevertheless, it is interesting to study this alternative measure of flow fluctuations and to compare it to the same quantity predicted by initial-state models. For head-on nucleus-nucleus collisions, N part ≈ 400, the prediction for eccentricity fluctuations σ ε n / ε n reaches the limit of √ 4/π − 1 ≈ 0.52 for the fluctuations-only scenario when the ε n distribution is described by a two-dimensional Gaussian function in the transverse plane [47]. In Ref. [17] it was shown that this limit is indeed reached by σ v n / v n for v 2 measured in the most central Pb+Pb collisions and for v 3 and v 4 over the entire centrality range. For the fluctuations-only scenario, the estimate F(v n ) should be close to one since then v n {4} ≈ 0. Thus, it is interesting to compare this alternative measure of flow fluctuations to the same quantity derived from the initial eccentricity distributions, F(ε n ). It can be seen by comparison with Eq. (16) that the quantity F depends not only on the second moment of the ε n distribution (as does σ ε n / ε n ), but also on the fourth moment and, therefore, can provide a more sensitive test of model assumptions. Figure 12 shows the p T dependence of the relative elliptic flow fluctuations calculated for different centrality intervals with Eq. (16), where v 2 {2} is replaced by v 2 {EP}. For all centrality intervals, except 2-5 %, the relative elliptic flow fluctuations depend only weakly on p T over the whole p T range, indicating that they are predominantly associated with fluctuations of the initial geometry. A similar p T dependence of the relative elliptic flow fluctuations was recently reported by the ALICE collaboration [9], although the ALICE results for the 0-5 % most central collisions show a much stronger p T dependence than the present measurement for the centrality interval 2-5 %. This discrepancy may be due to different contributions of non-flow effects to v 2 {EP} measured in the two experiments.
The quantity F(v n ) is further investigated as a function of the collision centrality using flow harmonics averaged  [17] is shown as open circles (marked in the legend as "EbyE") with the error bars denoting statistical and systematic uncertainties added in quadrature. The same quantity calculated from the initial eccentricity distributions obtained from the Glauber [37] and MC-KLN [46] models is shown by curves. Open squares show the CMS measurement of F(v 2 ) [23] over p T and η. The dependence of F(v 2 ) on N part is shown in Fig. 13. Two sets of measurements are shown:  [17]. The two measurements show similar centrality dependence, but the estimate based on the cumulant method is systematically smaller (by up to about 15 %) than that calculated from p(v 2 ).
F(v 2 ) can also be compared to σ v 2 / v 2 determined from the p(v 2 ) distribution. It was shown in Ref. [17] that the two measures of elliptic flow fluctuations agree for the most peripheral collisions. For semi-central collisions, σ v 2 / v 2 is systematically larger than F(v 2 ). A significant discrepancy F(v 2 ) shows a strong centrality dependence. Starting with the most peripheral collisions, it decreases with N part , reaching a minimum at N part ≈ 200 and then rises steeply with centrality up to about 1.0 for the most central collisions. A comparison to F( 2 ) calculated from the eccentricity distributions predicted by the Glauber and MC-KLN models is shown in Fig. 13. One can see that the MC-KLN model describes the measurements for Pb+Pb collisions with N part above 150 reasonably well, while the Glauber model significantly over-predicts the measured F(v 2 ) across the entire centrality range. Figure 13 also shows that our results are consistent with the CMS estimate of F(v 2 ) [23].
The study of F(v n ) is also performed for higher-order flow harmonics, n = 3, 4. Figure 14 shows F(v 3 ) (top plot) and F(v 4 ) (bottom plot) obtained using the cumulant method with v 3 {2} and v 4 {2} replaced by v 3 {EP} and v 4 {EP}, respectively. Large relative fluctuations of the third-and fourthorder harmonics, of the order of 0.7-0.8, are measured over the whole accessible centrality range, with a relatively weak centrality dependence. The results are consistent with the CMS measurement of F(v 3 ) [23] as well as with F calculated using v calc 3 {4, EbyE} and v calc 3 {2, EbyE} [17]. F( 3 ) and F( 4 ) obtained from the eccentricity distributions predicted by the Glauber and MC-KLN models are also shown. One can see that none of the models gives a consistent description of F(v 3 ) and F(v 4 ).

Summary
A measurement of flow harmonics of charged particles in Pb+Pb collisions at √ s NN = 2.76 TeV from the ATLAS experiment at the LHC is presented using a dataset of approximately 7 µb −1 collected in 2010. The analysis is based on the cumulant expansion of multi-particle azimuthal correlations, which suppresses correlations not related to the initialstate geometry. Another advantage of the cumulant method is that it provides several different measurements of the same harmonic, v n , allowing for estimation of the non-flow contributions and for consistency checks. The need for huge computing power in calculating multi-particle correlations is overcome by using the generating function formalism. Flow coefficients v n for n = 2, 3, 4 were obtained using two-and four-particle cumulants. In addition, for the elliptic flow (n = 2), the analysis is for the first time extended to the six-and eight-particle cumulants. The transverse momentum, pseudorapidity and centrality dependence of flow harmonics is presented. An attempt is also made to estimate the flow harmonic fluctuations using the measured v n {4} and v n obtained with the event-plane method.
The transverse momentum dependence of v 2 {2} shows significant non-flow contributions. This contribution is reduced in v 2 {EP}. The elliptic flow obtained with the fourparticle cumulants provides a measure of v 2 with nonflow correlations strongly suppressed. Using six-and eightparticle cumulants gives results consistent, within the errors, with those obtained with four-particle cumulants, indicating that four-particle cumulants efficiently suppress nonflow correlations. Similar conclusions can be drawn from the study of the η-dependence of the p T -integrated v 2 as well as from the centrality dependence of v 2 averaged over p T and η. As for v 2 , the higher-order flow harmonics determined using four-particle cumulants are significantly reduced compared to the measurement involving two-particle cumulants.
The flow harmonics, v n {4}, determined from the fourparticle cumulants increase sharply with p T reaching a maximum at 2-3 GeV. At higher transverse momenta, v 2 {4} decreases and beyond p T of about 7 GeV, it plateaus at the level of about 0.04 up to the highest accessible transverse momenta. The higher-order harmonics also decrease above 3 GeV and reach a value of about 0.03 when integrated over p T from 4 GeV to 20 GeV. The four-particle harmonics, v n {4}, are found to depend weakly on the pseudorapidity over the range |η| < 2.5.
The centrality dependence of the p T -and η-integrated v n {4} reveals a clear distinction between v 2 and the higherorder harmonics: v 2 strongly depends on the collision centrality, reflecting its sensitivity to the varying shape of the initial collision geometry, while v 3 and v 4 show only a weak centrality dependence, predominantly attributed to geometry fluctuations. Over the studied centrality range, except for the most central collisions, the measured v 3 and v 4 are much smaller than v 2 . For the most central collisions, a similar magnitude is measured for v 2 , v 3 and v 4 because all are dominated by geometry fluctuations. The measured v 2 {2k} for k = 2, 3, 4 and v 3 {4} are found to agree with the same coefficients calculated from the moments of the measured p(v 2 ) and p(v 3 ) distributions.
The relative flow harmonic fluctuations, F(v n ), defined in Eq. (16), are estimated using v n {EP} and v n {4}. For the elliptic flow harmonic, a strong centrality dependence is observed, following a trend similar to that exhibited by F as estimated from the p(v 2 ) distribution. In contrast, F(v 3 ) and F(v 4 ) show a weak centrality dependence. The large magnitudes of F obtained for third-and fourth-order harmonics, and also for the elliptic flow harmonic measured in the most central collisions, indicate the dominant role of initial-state fluctuations. The comparison to the same quantity derived from the initial-state eccentricity distributions, modelled by the Glauber and MC-KLN models, shows that none of these models can describe the flow harmonic fluctuations well, particularly for higher-order flow harmonics. Therefore, the measurements presented in this paper provide valuable constraints on models of initial spatial anisotropy and subsequent hydrodynamic evolution of systems produced in ion-ion collisions with nucleon-nucleon centre-of-mass energies at the TeV energy scale.