Measurement of the distributions of event-by-event flow harmonics in lead--lead collisions at sqrt(s_NN)=2.76 TeV with the ATLAS detector at the LHC

The distributions of event-by-event harmonic flow coefficients v_n for n=2-4 are measured in sqrt(s_NN)=2.76 TeV Pb+Pb collisions using the ATLAS detector at the LHC. The measurements are performed using charged particles with transverse momentum pT>0.5 GeV and in the pseudorapidity range |eta|<2.5 in a dataset of approximately 7 ub^-1 recorded in 2010. The shapes of the v_n distributions are described by a two-dimensional Gaussian function for the underlying flow vector in central collisions for v_2 and over most of the measured centrality range for v_3 and v_4. Significant deviations from this function are observed for v_2 in mid-central and peripheral collisions, and a small deviation is observed for v_3 in mid-central collisions. It is shown that the commonly used multi-particle cumulants are insensitive to the deviations for v_2. The v_n distributions are also measured independently for charged particles with 0.51 GeV. When these distributions are rescaled to the same mean values, the adjusted shapes are found to be nearly the same for these two pT ranges. The v_n distributions are compared with the eccentricity distributions from two models for the initial collision geometry: a Glauber model and a model that includes corrections to the initial geometry due to gluon saturation effects. Both models fail to describe the experimental data consistently over most of the measured centrality range.


Introduction
Heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) create hot, dense matter that is thought to be composed of strongly interacting quarks and gluons. A useful tool to study the properties of this matter is the azimuthal anisotropy of particle emission in the transverse plane [1,2]. This anisotropy is believed to result from pressure-driven anisotropic expansion (referred to as "flow") of the created matter, and is described by a Fourier expansion of the particle distribution in azimuthal angle φ, around the beam direction: where v n and Φ n represent the magnitude and phase of the n th -order anisotropy. These quantities can also be conveniently represented by the per-particle "flow vector" [2]: ⇀ v n = (v n cos nΦ n , v n sin nΦ n ). The angles Φ n are commonly referred to as the event plane (EP) angles.
In typical non-central [2] heavy ion collisions, the large and dominating v 2 coefficient is associated mainly with the "elliptic" shape of the nuclear overlap. However, v 2 in -1 -central (head-on) collisions and the other v n coefficients in general are related to various shape components of the initial state arising from fluctuations of the nucleon positions in the overlap region [3]. The amplitudes of these shape components, characterized by eccentricities ǫ n , can be estimated via a simple Glauber model from the transverse positions (r, φ) of the participating nucleons relative to their centre of mass [4,5]: ǫ n = r n cos nφ 2 + r n sin nφ 2 r n . (1. 2) The large pressure gradients and ensuing hydrodynamic evolution can convert these shape components into v n coefficients in momentum space. Calculations based on viscous hydrodynamics suggest that v n scales nearly linearly with ǫ n , for n < 4 [6]. The proportionality constant is found to be sensitive to properties of the matter such as the equation of state and the ratio of shear viscosity to entropy density [7,8]. In particular, the proportionality constant is predicted to decrease quickly with increasing shear viscosity [9]. Hence detailed measurements of v n coefficients and comparisons with ǫ n may shed light on the collision geometry of the initial state and transport properties of the created matter [10,11].
Significant v n coefficients have been observed for n ≤ 6 at RHIC and the LHC [12][13][14][15][16][17][18]. These observations are consistent with small values for the ratio of shear viscosity to entropy density, and the existence of sizable fluctuations in the initial state. Most of these measurements estimate v n from the distribution of particles relative to the event plane, accumulated over many events. This event-averaged v n mainly reflects the hydrodynamic response of the created matter to the average collision geometry in the initial state. More information, however, can be obtained by measuring ⇀ v n or v n on an event-by-event (EbyE) basis.
Some properties of the v n distributions, such as the mean v n , the standard deviation (hereafter referred to as "width") σ vn , the relative fluctuation σ vn / v n , and the root-meansquare v 2 n ≡ v n 2 + σ 2 vn , were previously estimated from a Monte Carlo template fit method [19], or two-and four-particle cumulant methods [20][21][22]. The value of σ v 2 / v 2 was measured to be 0.3-0.7 in different centrality bins. However, these methods are reliable only for σ vn ≪ v n , and are subject to significant systematic uncertainties. In contrast, v n , σ vn and higher-order moments can be calculated directly from the full v n distributions.
The EbyE distributions of ⇀ v n or v n also provide direct insight into the fluctuations in the initial geometry [23]. If fluctuations of ⇀ v n relative to the underlying flow vector associated with the average geometry, ⇀ v RP n , in the reaction plane 1 (RP) [23,24] are described by a two-dimensional (2D) Gaussian function in the transverse plane, then the probability density of ⇀ v n can be expressed as: Model calculations show that this approximation works well for central and mid-central collisions [23,25]. Integration of this 2D Gaussian over the azimuthal angle gives the 1 The reaction plane is defined by the impact parameter vector and the beam axis.
-2 -one-dimensional (1D) probability density of v n = | ⇀ v n | in the form of the Bessel-Gaussian function [7,26]: (1.4) where I 0 is the modified Bessel function of the first kind. Additional smearing to eq. (1.3) also arises from effects of the finite number of particles produced in the collision. If it is Gaussian, this smearing is expected to increase the observed δ vn value, but the value of v RP n should be stable. The parameters v RP n and δ vn in eq. (1.4) are related to v n and σ vn , and can be estimated directly from a fit of the measured p(v n ) distribution with eq. (1.4). For small fluctuations δ vn ≪ v RP n [23]: For large fluctuations δ vn ≫ v RP n (e.g. in central collisions), eqs. (1.3) and (1.4) can be approximated by: 6) which is equivalent to the "fluctuation-only" scenario, i.e. v RP n = 0. In this case, both the mean and the width are controlled by δ vn [27]: and hence: σ vn v n = 4 π − 1 = 0.523, v 2 n = 2 √ π v n = 1.13 v n . (1.8) In the intermediate case, δ vn ≈ v RP n , a more general approximation to eq. (1.4) can be used via a Taylor expansion of the Bessel function, I 0 (x) = e x 2 /4 1 − x 4 /64 + O(x 6 ) : Defining α ≡ δ vn /v RP n , eqs. (1.9) and (1.10) imply that for v n ≪ 2 √ 2δ vn α, the shape of the distribution is very close to that of eq. 1.6, except for a redefinition of the width. For example, the deviation from the fluctuation-only scenario is expected to be less than 10% over the range v n < 1.6δ vn α. Hence the reliable extraction of v RP selection requires a time difference |∆t| < 3 ns between the MBTS trigger counters on either side of the interaction point to suppress non-collision backgrounds. A coincidence between the ZDCs at forward and backward pseudorapidity is required to reject a variety of background processes, while maintaining high efficiency for non-Coulomb processes. Events satisfying these conditions are required to have a reconstructed primary vertex with z vtx within 150 mm of the nominal centre of the ATLAS detector. About 48 million events pass the requirements for the analysis.
The Pb+Pb event centrality is characterized using the total transverse energy (ΣE T ) deposited in the FCal over the pseudorapidity range 3.2 < |η| < 4.9 measured at the electromagnetic energy scale [34]. A larger ΣE T value corresponds to a more central collision. From an analysis of the ΣE T distribution after all trigger and event selections, the sampled fraction of the total inelastic cross section is estimated to be (98±2)% [35]. The uncertainty associated with the centrality definition is evaluated by varying the effect of the trigger, event selection and background rejection requirements in the most peripheral FCal ΣE T interval [35]. The FCal ΣE T distribution is divided into a set of 5%-wide percentile bins, together with five 1%-wide bins for the most central 5% of the events. A centrality interval refers to a percentile range, starting at 0% for the most central collisions. Thus the 0-1% centrality interval corresponds to the most central 1% of the events; the 95-100% centrality interval corresponds to the least central (i.e. most peripheral) 5% of the events. A standard Glauber model Monte Carlo analysis is used to estimate the average number of participating nucleons, N part , for each centrality interval [35,36]. These numbers are summarized in table 1.
The v n coefficients are measured using tracks reconstructed in the ID that have p T > 0.5 GeV and |η| < 2.5. To improve the robustness of track reconstruction in the highmultiplicity environment of heavy ion collisions, more stringent requirements on track quality, compared to those defined for proton-proton collisions [37], are used. At least 9 hits in the silicon detectors (compared to a typical value of 11) are required for each -5 - track, with no missing Pixel hits and not more than 1 missing SCT hit, after taking into account the known non-operational modules. In addition, at its point of closest approach the track is required to be within 1 mm of the primary vertex in both the transverse and longitudinal directions [33]. The efficiency, ǫ(p T , η), of the track reconstruction and track selection cuts is evaluated using Pb+Pb Monte Carlo events produced with the HIJING event generator [38]. The response of the detector is simulated using GEANT4 [39] and the resulting events are reconstructed with the same algorithms as applied to the data. The absolute efficiency increases with p T by 6% between 0.5 GeV and 0.7 GeV, and varies only weakly for p T > 0.7 GeV. However, the efficiency varies more strongly with η and event multiplicity [40]. For p T > 0.7 GeV, it ranges from 70% at η = 0 to 50% for |η| > 2 in peripheral collisions, while it ranges from 70% at η = 0 to about 40% for |η| > 2 in central collisions. Contributions of fake tracks from random combinations of hits are generally negligible, reaching only 0.1% for |η| < 1 for the highest multiplicity events. This rate increases slightly at large |η|.

Method and data analysis
To illustrate the level of EbyE fluctuations in the data, the top panels of figure 1 show the azimuthal distribution of charged particle tracks with p T > 0.5 GeV for three typical events in the 0-5% centrality interval. The corresponding track-pair ∆φ distributions from the same events are shown in the bottom panels. For each pair of particles two ∆φ entries, |φ 1 − φ 2 | and −|φ 1 − φ 2 |, are made each with a weight of 1/2, and then folded into the [−0.5π, 1.5π] interval. Rich EbyE patterns, beyond the structures in the event-averaged distributions shown by the solid points (arbitrary normalization), are observed. These EbyE distributions are the inputs to the EbyE v n analyses.
The azimuthal distribution of charged particles in an event is written as a Fourier -6 -series, as in eq. (1.1): n,y 2 , v obs n,x = v obs n cos nΨ obs n = cos nφ , v obs n,y = v obs n sin nΨ obs n = sin nφ , (4.2) where the averages are over all particles in the event for the required η range. The v obs n is the magnitude of the observed EbyE per-particle flow vector: ⇀ v obs n = (v obs n,x , v obs n,y ). In the limit of very large multiplicity and in the absence of non-flow effects, it approaches the true flow signal: v obs n → v n . The key issue in measuring the EbyE v n is to determine the response function p(v obs n |v n ), which can be used to unfold the smearing effect due to the finite number of detected particles. Possible non-flow effects from short-range correlations, such as resonance decays, Bose-Einstein correlations and jets, also need to be suppressed.
The rest of this section sets out the steps to obtain the unfolded v n distribution. Since the data-driven unfolding technique has rarely been used in the study of flow phenom-  Figure 1. Single-track φ (top) and track-pair ∆φ (bottom) distributions for three typical events (from left to right) in the 0-5% centrality interval. The pair distributions are folded into the [−0.5π, 1.5π] interval. The bars indicate the statistical uncertainties of the foreground distributions, the solid curves indicate a Fourier parameterization including the first six harmonics: dN/dφ = A(1 + 2 6 i=1 c n cos n(φ − Ψ n )) for single-track distributions and dN/d∆φ = A(1 + 2 6 i=1 c n cos n(∆φ)) for track-pair distributions, and the solid points indicate the eventaveraged distributions (arbitrary normalization). ena, details are provided to facilitate the understanding of the methods and systematic uncertainties. Section 4.1 explains how v obs n and the associated response function can be obtained from the EbyE single-particle distributions, such as those shown in the top panels of figure 1. Section 4.2 describes how v obs n and the response function can be obtained from EbyE two-particle correlations (2PC), similar to those shown in the lower panels of figure 1. In this paper the 2PC approach is used primarily as a consistency check. The Bayesian unfolding procedure, applicable to either the single-particle or 2PC data, is described in section 4.3. The performance of the unfolding is described in section 4.4, while the systematic uncertainties are discussed in section 4.5.

Single-particle method
The azimuthal distribution of particles in figure 1 needs to be corrected for non-uniform detector acceptance. This is achieved by dividing the foreground distribution (S) of a given event by the acceptance function (B) obtained as the φ distribution of all tracks in all events (see top panels of figure 1): , (4.4) where the index i runs over all tracks in an event, ǫ(η, p T ) is the tracking efficiency for a given centrality interval, and v det n,x and v det n,y are Fourier coefficients of the acceptance function in azimuth, also weighted by the inverse of the tracking efficiency. The influence of the structures in the acceptance function can be accounted for by taking the leading-order term of the Taylor expansion of 1/B(φ) in terms of cos nφ and sin nφ: 5) where the values of v det n,x or y are less than 0.007 for n = 2-4. The higher-order corrections to eq. (4.5) involve products of v raw n,x or y and v det n,x or y . They have been estimated and found to have negligible impact on the final v n distributions for n = 2-4. Figure 2 shows the distribution of the EbyE per-particle flow vector ⇀ v obs 2 = (v obs 2,x , v obs 2,y ) and v obs 2 obtained for charged particles with p T > 0.5 GeV in the 20-25% centrality interval. The azimuthal symmetry in the left panel reflects the random orientation of ⇀ v obs 2 because of the random orientation of the impact parameter. Due to the finite track multiplicity, the measured flow vector is expected to be smeared randomly around the true flow vector by a 2D response function p( ⇀ v obs n | ⇀ v n ). In order to determine p( ⇀ v obs n | ⇀ v n ), the tracks in the entire inner detector (referred to as full-ID) for a given event are divided into two subevents with symmetric η range, η > 0 and η < 0 (labelled by a and b and referred to as half-ID). The two half-IDs have the same average track multiplicity to within 1%. The distribution of flow vector differences between the two subevents, p sub ( ⇀ v obs n ) a − ( ⇀ v obs n ) b , is then obtained and is shown in the left panel of figure 3. The physical flow signal cancels in this distribution such that it  Figure 3. Left: The distribution of the difference between the EbyE per-particle flow vectors of the two half-IDs for events in the 20-25% centrality interval for n = 2. Middle: The x-projection overlaid with a fit to a Gaussian. Right: The y-projection overlaid with a fit to a Gaussian. The width from the fit, δ 2SE , and the quality of the fit, χ 2 /DOF, are also shown.
contains mainly the effects of statistical smearing and non-flow. The middle and right panels of figure 3 show the x-and y-projections of the distribution, together with fits to a Gaussian function. The fits describe the data very well (χ 2 /DOF≈ 1) across five orders of magnitude with the same widths in both directions, implying that the smearing effects and any effects due to non-flow short-range correlations are purely statistical. This would be the case if either the non-flow effects are small and the smearing is mostly driven by finite particle multiplicity, or the number of sources responsible for non-flow is proportional to the multiplicity and they are not correlated between the subevents [31]. The latter could be true for resonance decays, Bose-Einstein correlations, and jets. Similar behaviour is observed for all harmonics up to centrality interval 50-55%. Beyond that the distributions are found to be described better by the Student's t-distribution, which is a general probability density -9 -function for the difference between two estimates of the mean from independent samples. The t-distribution approaches a Gaussian distribution when the number of tracks is large.
Denoting the width of these 1D distributions by δ 2SE , the widths of the response functions for the half-ID and the full-ID are δ 2SE / √ 2 and δ 2SE /2, respectively. The response functions themselves can be obtained by rescaling the left panel of figure 3 in both dimensions by constant factors of 2 and √ 2 for the full-ID and half-ID, respectively [1,19,23]: which reduces to a Bessel-Gaussian function in 1D similar to eq. (1.4): The difference between the observed and the true signal, denoted by s = v obs n −v n , accounts for the statistical smearing. The similarity between eq. (4.7) and eq. (1.4) is a direct consequence of the 2D Gaussian smearing. However, the smearing leading to eq. (4.7) is due to the finite charge-particle multiplicity, while the smearing leading to eq. (1.4) is due to the intrinsic flow fluctuations associated with the initial geometry. Hence the smearing in eq. (4.7) is expected to increase the observed δ vn value but the value of v RP n should be relatively stable.
The analytical expression eq. (4.7) can be used to unfold the v obs n distribution, such as that shown in the right panel of figure 2. Alternatively, the measured distribution p(v obs n |v n ), taking into account sample statistics, can be used directly in the unfolding. This measured distribution is obtained by integrating out the azimuthal angle in the 2D response function, and the response function is obtained by rescaling the left panel of figure 3 as described earlier. This approach is more suitable for peripheral collisions where the analytical description using eq. (4.7) is not precise enough.

Two-particle correlation method
The EbyE two-particle correlation (2PC) method starts from the ∆φ information in each event, where ∆φ is calculated for each pair of charged tracks as described at the start of section 4. In order to reduce the effect of short-range correlations in η, the tracks in each pair are taken from different half-IDs. This procedure corresponds to convolving the azimuthal distributions of single particles in the two half-IDs: -10 -where A n = cos n∆φ and B n = sin n∆φ . The parameters A n and B n are calculated by averaging over the pairs in each event, with each track weighted by the tracking efficiency, as in eq. (4.3). An EbyE track-pair variable v obs n,n is subsequently calculated for each event: v obs n,n ≡ A 2 The observed flow signal from the two-particle correlation analysis is then calculated as: where s a = v obsa n − v n and s b = v obs b n − v n are independent variables described by the probability distribution in eq. (4.7) with δ = δ 2SE / √ 2. The response function for v obs,2PC n is very different from that for the single-particle method, but nevertheless can be either calculated analytically via eq. (4.7) or generated from the measured distribution such as that shown in figure 3. For small v n values, the s a s b term dominates eq. (4.10) and the distribution of v obs,2PC n is broader than v obs n . For large v n values, the distributions of s a and s b are approximately described by Gaussian functions and hence: (4.11) where the fact that the average of two Gaussian random variables is equivalent to a Gaussian with a narrower width has been used, and the smearing of the flow vector for the half-IDs (s a and s b ) is taken to be a factor of √ 2 broader than that for the full-ID (s). Hence the distribution of v obs,2PC n is expected to approach the v obs n distribution of the full-ID when v n ≫ δ 2SE / √ 2.

Unfolding procedure
In this analysis, the standard Bayesian unfolding procedure [41], as implemented in the RooUnfold framework [42], is used to obtain the v n distribution. In this procedure, the true v n distribution ("cause"ĉ) is obtained from the measured v obs n or v obs,2PC n distribution ("effect"ê) and the response function A ji ≡ p(e j |c i ) (i,j are bins) as: where the unfolding matrixM 0 is determined from the response function and some initial estimate of the true distributionĉ 0 (referred to as the prior). The matrixM 0 is used to obtain the unfolded distributionĉ 1 andM 1 , and the process is then iterated. More iterations reduce the dependence on the prior and give results closer to the true distributions but with increased statistical fluctuations. Therefore the number of iterations N iter is adjusted according to the sample size and binning. The prior can be chosen to be the v obs n distribution from the full-ID for the single-particle unfolding, or the v obs,2PC n distribution obtained by convolving the two half-IDs (eq. (4.8)) for the 2PC unfolding. However, a -11 -more realistic prior can be obtained by rescaling the v obs n in each event by a constant factor R to obtain a distribution with a mean that is closer to that of the true distribution: where v obs n and σ obs vn are the mean and the standard deviation of the v obs n distribution, respectively, and v EP n is measured using the FCal event plane method from ref.
[16] with the same dataset and the same track selection criteria. The EP method is known to measure a value between the mean and the root-mean-square of the true v n [25, 28] (see figure 13): (4.14) The lower limit is reached when the resolution factor [16] used in the EP method approaches one, and the upper limit is reached when the resolution factor is close to zero. Thus f = 0 (default choice) gives a prior that is slightly broader than the true distribution, f = 1 gives a distribution that has a mean close to the true distribution, and f > 1 typically gives a distribution that is narrower than the true distribution. The unfolding process in this analysis has several beneficial features: 1. The response function is obtained entirely from the data using the subevent method described above (eq. (4.6)).
2. Each event provides one entry for the v obs n distribution and the response function (no efficiency loss), and these distributions can be determined with high precision (from about 2.4 million events for each 5% centrality interval).
3. The unfolded v n distribution is always narrower than the input v obs n distribution, and is fully contained within the available range.

Unfolding performance
This section describes the unfolding based on the single-particle method and summarises a series of checks used to verify the robustness of the results: a) the number of iterations used, b) comparison with results obtained from a smaller η range, c) variation of the priors, d) comparison with the results obtained using the 2PC method, and e) estimation of possible biases due to short-range correlations by varying the η gap between the two half-IDs. Only a small subset of the available figures is presented here; a complementary selection can be found in appendix A.
The left and middle panels of figure 4 show the convergence behaviour of the unfolding based on the single-particle method for v 2 in the 20-25% centrality interval measured with the full-ID. Around the peak of the distribution, the results converge to within a few percent of the final answer by N iter = 4, but the convergence is slower in the tails and there are small, systematic improvements at the level of a few percent for N iter ≥ 8. The -12 -refolded distributions (right panel), obtained by convolving the unfolded distributions with the response function, indicate that convergence is reached for N iter ≥ 8. Figures 5 and 6 show similar distributions for v 3 and v 4 . The performance of the unfolding is similar to that shown in figure 4, except that the tails of the unfolded distributions require more iterations to converge. For example, figure 6 suggests that the bulk region of the v 4 distributions has converged by N iter = 32, but the tails still exhibit some small changes up to N iter = 64. The slower convergence for higher-order harmonics reflects the fact that the physical v n signal is smaller for larger n, while the values of δ 2SE change only weakly with n. These studies are repeated for all centrality intervals. In general, more iterations are needed for peripheral collisions due to the increase in δ 2SE .
The statistical uncertainties in the unfolding procedure are verified via a resampling technique [43]. For small N iter , the statistical uncertainties as given by the diagonal entries of the covariance matrix are much smaller than √ N , where N is the number of entries in each bin, indicating the presence of statistical bias in the prior. However, these uncertainties increase with N iter , and generally approach √ N for 64 ≤ N iter ≤ 128. In this analysis, the centrality range for each harmonic is chosen such that the difference between N iter = 32 and N iter = 128 is less than 10%. The centrality ranges are 0-70% for v 2 , 0-60% for v 3 and 0-45% for v 4 .
The robustness of the unfolding procedure is checked by comparing the results measured independently for the half-ID and the full-ID. The results are shown in figure 7. Despite the large differences between their initial distributions, the final unfolded results agree to within a few percent in the bulk region of the unfolded distribution, and they are nearly indistinguishable on a linear scale. This agreement also implies that the influence due to the slight difference (up to 1%) in the average track multiplicity between the two subevents is small. Systematic differences are observed in the tails of the distributions for v 4 , especially in peripheral collisions, where the half-ID results are slightly broader. This behaviour reflects mainly the residual non-convergence of the half-ID unfolding, since the centrality: 20-25% f=0 as prior  centrality: 20-25% f=0 as prior  centrality: 20-25% f=0 as prior  response function is a factor of √ 2 broader than that for the full-ID. A wide range of priors has been tried in this analysis, consisting of the measured v obs n distribution and the six rescaled distributions defined by eq. (4.13). Figure 8 compares the convergence behaviour of these priors for v 3 in the 20-25% centrality interval. Despite the significantly different initial distributions, the unfolded distributions converge to the same answer, to within a few percent, for N iter ≥ 16. A prior that is narrower than the unfolded distribution leads to convergence in one direction, and a broader prior leads to convergence from the other direction. Taken together, these checks allow a quantitative evaluation of the residual non-convergence in the unfolded distributions. Figure 9 compares the convergence behaviour between unfolding of single-particle v obs   in the tails of the v 4 distribution (bottom-right panel) are due mainly to the remaining non-convergence in the 2PC method, which has a broader response function than the single-particle method. One important issue in the EbyE v n study is the extent to which the unfolded results are biased by non-flow short-range correlations, which may influence both the v obs n distributions and the response functions. This influence should contribute to both the v obs  Figure 9. Comparison of the input distributions (solid symbols) and unfolded distributions for N iter = 128 (open symbols) between the single-particle unfolding and 2PC unfolding in the 20-25% centrality interval for v 2 (left panels), v 3 (middle panels) and v 4 (right panels). The ratios of 2PC to single-particle unfolded results are shown in the bottom panels.
by the consistency between the single-particle and 2PC methods (figure 9), which have different sensitivities to the non-flow effects. Furthermore, both unfolding methods have been repeated requiring a minimum η gap between the two subevents used to obtain the input distributions and the response functions. Six additional cases, requiring η gap = 0.2, 0.4, 0.6, 0.8, 1.0, 2.0, are studied and the results are compared with each other. The unfolded v n distributions are observed to narrow slightly for larger η gap , reflecting the fact that the true v n decreases slowly at large |η| [16] and a larger η gap on average selects particles at large |η|. However, the results are always consistent between the two methods independent of the η gap value used. This consistency supports further the conclusion that the influence of the short-range non-flow correlations on the final unfolded results is not significant.
The dependence of the EbyE v n on the p T of the charged particles has also been checked: the particles are divided into those with 0.5 < p T < 1 GeV and those with p T > 1 GeV, and the EbyE v n measurements are repeated independently for each group of particles. About 60% of detected particles have 0.5 < p T < 1 GeV, and this fraction varies weakly with centrality. The unfolding performance is found to be slightly worse for charged particles with 0.5 < p T < 1 GeV than for those with p T > 1 GeV, due to their much smaller v n signal. Hence the v n range of the unfolded distribution for the final results is chosen separately for each p T range.
The final v n distributions are obtained using the single-particle unfolding with the full-ID and N iter = 128, separately for charged particles in the two aforementioned p T ranges and the combined p T range. The prior is obtained by rescaling the v obs n distribution according to eq. (4.13) with f = 0, and the response function is measured from correlations of the two half-IDs with no η gap in between.

Systematic uncertainties
The systematic uncertainties associated with the unfolding procedure include contributions from the residual non-convergence, dependence on the prior, uncertainty in the response function, the difference between the single-particle method and the 2PC method, and the tracking efficiency. The residual non-convergence is estimated from the difference between the results for N iter = 32 and N iter = 128 or between the results for the half-ID and full-ID. These two estimates are strongly correlated, so in each bin of the unfolded distribution the larger of the two is used. The dependence on the prior is taken as the difference between the results for f = 0 and f = 2.5. Note that the prior for f = 0 (f = 2.5) is broader (narrower) than the final unfolded distributions. The uncertainty of the response function is estimated from the difference between results obtained using the analytical formula eq. (4.7) and results obtained using the measured distribution, as well as the change in the results when the small dependence of δ 2SE on the observed v obs n is taken into account. Results for the p T -dependence of the v n distributions (see section 5 and figure 11) show that the mean values vary with p T , but, after they are rescaled to a common mean, the resulting shapes are almost identical. Motivated by this finding, every source of systematic uncertainty is decomposed into two components: the uncertainty associated with the v n or the v n -scale, and the uncertainty in the shape after adjustment to the same v n or the adjusted v n -shape. The uncertainties are then combined separately for the v n -scale and the adjusted v n -shape. Most shape uncertainties can be attributed to σ vn , such that the remaining uncertainties on the adjusted v n -shape are generally smaller.
To estimate the uncertainty due to the tracking efficiency, the measurement is repeated without applying the efficiency re-weighting. The final distributions are found to have almost identical shape, while the values of v n and σ vn increase by a few percent. This increase can be ascribed mainly to the smaller fraction of low-p T particles, which have smaller v n , so this increase should not be considered as a systematic uncertainty on the v n -scale. Instead, the scale uncertainty is more appropriately estimated from the change in the v EP n when varying the efficiency correction within its uncertainty range. On the other hand, small changes are observed for σ vn / v n and the adjusted v n -shape. Since these changes are small, they are conservatively included in the total systematic uncertainty in the v n -shape.
Additional systematic uncertainties include those from the track selection, dependence on the running periods, and trigger and event selections. These account for the influence of fake tracks, the instability of the detector acceptance and efficiency, and variation of the centrality definition, respectively. All three sources of systematic uncertainty are expected to change only the v n -scale but not the v n -shape, and they are taken directly from the published v EP n measurement [16]. Table 2 summarises the various systematic uncertainties on v n , σ vn and σ vn / v n , for charged particles with p T > 0.5 GeV. The uncertainties on v n and σ vn are strongly correlated, which in many cases leads to a smaller and asymmetric uncertainty on σ vn / v n . In most cases, the uncertainties are dominated by tracking efficiency, as well as the track selection, trigger, and run stability assessed in ref. [16]. The uncertainties associated with -17 -    Table 2. Summary of systematic uncertainties in percent for v n , σ vn and σ vn / v n (n = 2-4) obtained using charged particles with p T > 0.5 GeV. The uncertainties for v n and σ vn are similar so the larger of the two is quoted. The uncertainties associated with track selection, the trigger and stability are taken from ref. [16]. Most uncertainties are asymmetric; the quoted numbers refer to the maximum uncertainty range spanned by various centrality intervals in each group.
the unfolding procedure are usually significant only in peripheral collisions, except for v 4 , where they are important across the full centrality range. The uncertainties in table 2 have also been evaluated separately for charged particles with 0.5 < p T < 1 GeV and p T > 1 GeV. In general, the systematic uncertainties are larger for the group of particles with 0.5 < p T < 1 GeV, due mainly to the increased contributions from residual nonconvergence and the choice of priors. The final v n distributions are shown over a v n range that is chosen such that the statistical uncertainties in all bins are less than 15%, and the results obtained with the default setups between N iter = 32 and N iter = 128 are consistent within 10%. The systematic uncertainties on the adjusted shape from the sources discussed above are then combined to give the final uncertainty for the v n -shape. Within the chosen v n ranges, the statistical uncertainties are found to be always smaller than the systematic uncertainties for the v nshape, and the integrals of the v n distributions outside these ranges are typically < 0.5% for the 5%-wide centrality intervals and < 1% for the 1%-wide centrality intervals. Figure 10 shows the probability density distributions of the EbyE v n in several centrality intervals obtained for charged particles with p T > 0.5 GeV. The shaded bands indicate the systematic uncertainties associated with the shape. These uncertainties are strongly correlated in v n : the data points are allowed to change the shape of the distribution within the band while keeping v n unchanged. The v n distributions are found to broaden from central to peripheral collisions (especially for v 2 ), reflecting the gradual increase of the magnitude of v n for more peripheral collisions [16,33]. The shape of these distributions changes quickly with centrality for v 2 , while it changes more slowly for higher-order harmonics. These distributions are compared with the probability density function obtained from eq. (1.6) (v RP n = 0), which represents a fluctuation-only scenario for v n . These functions, indicated by the solid curves, are calculated directly from the measured v n values via eq. (1.7) for each distribution. The fluctuation-only scenario works reasonably well for v 3 and v 4 over the measured centrality range, but fails for v 2 except for the most central 2% of collisions, i.e. for the centrality interval 0-2%. Hence for v 2 the solid curve representing the fluctuation-only scenario is shown only for the 0-1% centrality interval (the data for the 1-2% interval are not shown). However, there is a small systematic difference between the data and the curve in the tails of the v 3 distributions in mid-central collisions, with a maximum difference of two standard deviations. Using eq. (1.9), this difference is compatible with a non-zero v RP 3 similar to the findings reported in Ref. [44]. Futhermore, since the measured v 4 distribution covers only a limited range (v 4 3δ v 4 ), a non-zero v RP 4 on the order of δ v 4 can not be excluded by this analysis based on eq. (1.9). Figure 11 compares the unfolded v n distributions for charged particles in three p T ranges: p T > 0.5 GeV, 0.5 < p T < 1 GeV and p T > 1 GeV. The v n distributions for p T > 1 GeV are much wider than for 0.5 < p T < 1 GeV, reflecting the fact that the v n values increase strongly with p T in this region [16]. However, once these distributions are rescaled to the same v n as shown in the lower row of figure 11, their shapes are quite similar except in the tails of the distributions for n = 2. This behaviour suggests that the hydrodynamic response of the medium to fluctuations in the initial geometry is nearly independent of p T in the low-p T region; it also demonstrates the robustness of the unfolding performance against the change in the underlying v n distributions and response functions. Figure 12 shows a summary of the quantities derived from the EbyE v n distributions, i.e. v n , σ vn and σ vn / v n , as a function of N part . The shaded bands represent the total systematic uncertainties listed in table 2, which generally are asymmetric. Despite the strong p T dependence of v n and σ vn , the ratio σ vn / v n is relatively stable. For v 2 , the value of σ vn / v n varies strongly with N part , and reaches a minimum of about 0.34 at -19 -  Rescaled to the same mean Figure 11. Top panels: The unfolded distributions for v n in the 20-25% centrality interval for charged particles in the p T > 0.5 GeV, 0.5 < p T < 1 GeV and p T > 1 GeV ranges. Bottom panels: same distributions but rescaled horizontally so the v n values match that for the p T > 0.5 GeV range. The shaded bands represent the systematic uncertainties on the v n -shape.

Results
N part ∼ 200, corresponding to the 20-30% centrality interval. For v 3 and v 4 , the values of σ vn / v n are almost independent of N part , and are consistent with the value expected from the fluctuation-only scenario ( 4/π − 1 via eq. (1.8) as indicated by the dotted lines), except for a small deviation for v 3 in mid-central collisions. This limit is also reached for v 2 in the most central collisions as shown by the top-right panel of figure 12. Figure 13 compares the v n and v 2 n with the v EP n measured using the FCal event -20 - are in between v 2 and v 2 2 . As expected [25,28], and as discussed in section 4.3, they approach v 2 2 only in peripheral collisions, where the resolution factor used in the EP method is small.
The results in figures 10 and 12 imply that the distributions of v 2 in central collisions (centrality interval 0-2%), and of v 3 and v 4 in most of the measured centrality range are described by a 2D Gaussian function of ⇀ v n centred around the origin (eq. 1.6). On the other hand, the deviation of the v 2 distribution from such a description in other centrality intervals suggests that the contribution associated with the average geometry, v RP 2 , becomes important. In order to test this hypothesis, the v 2 distributions have been fitted to the Bessel-Gaussian function eq. (1.4), with v RP 2 not constrained to be zero. The results of the fit are shown in figure 14 for various centrality intervals. The fit works reasonably well up to the 25-30% centrality interval, although systematic deviations in the tails are apparent already in the 15-20% centrality interval. The deviations increase steadily for -21 - more peripheral collisions, which may be due to the fact that the fluctuations of ǫ 2 (eq. 1.2) are no longer Gaussian in peripheral collisions where N part is small [25] (see also figure 18).
The values of v RP 2 and δ v 2 can be estimated from these fits. Since the value of v RP 2 varies rapidly with N part , especially in central collisions, the extracted v RP 2 and δ v 2 values can be affected by their spreads if too broad a centrality interval is used. The effect of the centrality binning has been checked and corrected as follows. Taking the 20-25% interval as an example, the results obtained using the full centrality range within this interval are compared to the results obtained from the average of the five individual 1% intervals: 20-21%,...,24-25%. This procedure has been carried out for each 5% centrality interval, and the difference is found to be significant only for the 0-5% and 5-10% intervals, and negligible for all the others. For the 0-5% interval, results are reported in the individual 1% bins. For the 5-10% bin, results are averaged over the five individual 1% bins.
As a cross-check, the Bessel-Gaussian fits are also performed on the v obs 2 distributions before the unfolding. Systematic devations are also observed between the fit and the v obs 2 data, but the deviations are smaller than those shown in figure 14. The value of v RP 2 from the v obs 2 distribution is found to agree to within a few percent with that from the unfolded v 2 distribution, while the value of δ v 2 from the v obs 2 distribution is significantly larger. This behaviour is expected since the smearing by the response function (eq. 4.7) increases mainly the width, and the value of v RP 2 should be stable. Figure 15 shows the v RP 2 and δ v 2 values extracted from the v 2 distributions as a function of N part . They are compared with values of v 2 and σ v 2 obtained directly from the v 2 distributions. The v RP 2 value is always smaller than the value for v 2 , and it decreases to -22 - zero in the 0-2% centrality interval, consistent with the results shown in figure 10. The value of δ v 2 is close to σ v 2 except in the most central collisions. This behaviour leads to a value of δ v 2 /v RP 2 larger than σ v 2 / v 2 over the full centrality range as shown in the right panel of figure 15. The value of δ v 2 /v RP 2 decreases with N part and reaches a minimum of 0.38 ± 0.02 at N part ≈ 200, but then increases for more central collisions. The two points for the 0-1% and 1-2% centrality intervals are omitted as the corresponding v RP 2 values are consistent with zero.
According to eq. (1.5), when the relative fluctuations are small the value of v n can be approximated by: Similarly, the value of v RP n can be estimated from σ vn and v n without relying on the fit: Both relations are shown in the left panel of figure 15 for v 2 . Good agreement with the data is observed for 100 < N part < 350, corresponding to 5-45% centrality interval in table 1. However, systematic deviations are observed both in central collisions where fluctuations are dominant, and in peripheral collisions where the Bessel-Gaussian function fails to describe the shape of the v 2 distribution.
Multi-particle cumulant methods have been widely used to estimate the values of v RP 2 and δ v 2 , as well as to study the deviation of the v 2 distribution from a Bessel-Gaussian description [21]. The second-order coefficients for the 2k th -order cumulants, v 2 {2k}, are generally calculated via 2k-particle correlations [20]. When the non-flow effects are small, these coefficients can also be calculated analytically from the measured v 2 distributions, labelled as v calc 2 {2k}. The first four terms can be expressed [23] as: The last part of each equation is exact when the v 2 distribution follows the Bessel-Gaussian function eq. (1.4). In this case only the first two cumulants are independent, and the higher-order cumulants do not provide new information. For the same reason, it has been argued that differences between the higher-order cumulants v 2 {2k} (k > 1) are sensitive to non-Gaussian behaviour in the underlying v 2 distribution [23]. Figure 16 shows cumulants agree with the results obtained from the fit, except for N part < 100. However, the calculated coefficients from higher-order cumulants agree with each other over the whole centrality range, despite the poor description of the v 2 distribution by the Bessel-Gaussian function for the centrality interval 25-70%, shown in figure 14. It therefore follows that similar values of v 2 {2k} obtained from higher-order cumulants are no guarantee that the distribution is well described by a Bessel-Gaussian in non-central collisions.
In previous analyses based on cumulants [7,22], several other quantities have been used to study the v 2 fluctuations. The definitions of these quantities are given below and their physical meaning can be understood from eq. (5.3): (5.6) where F 1 , F 2 , and F 3 are calculated from the unfolded distributions, using eq. (5.3). The approximation for F 3 is valid when v RP 2 ≫ δ v 2 . In central collisions where v RP 2 ≪ δ v 2 , the value of F 3 is expected to approach one. Figure 17 compares the calculated values of F 1 , F 2 and F 3 to the rightmost expressions in eqs. (5.4)-(5.6), using δ v 2 , v RP 2 obtained from fits to the Bessel-Gaussian function, and the mean of the unfolded distribution. The value of F 1 is between σ v 2 and δ v 2 . The quantities F 2 and F 3 show similar N part dependence as σ v 2 / v 2 and δ v 2 /v RP 2 , however significant discrepancies are observed, especially in the most central collisions where the flow fluctuation is dominant.
-25 -  million events have been generated and grouped into centrality intervals according to the impact parameter. The ǫ 2 distribution for each centrality interval is rescaled to match the v 2 of the data, and then normalized to form a probability density function. Since v 2 is expected to be proportional to ǫ 2 in most hydrodynamic calculations [6], the deviations between the v 2 distributions and the rescaled ǫ 2 distributions can be used to improve the modeling of the initial geometry. Figure 18 shows that the rescaled ǫ 2 distributions describe the data well for the most central collisions, but not so well for non-central collisions. In peripheral collisions, both the Glauber and MC-KLN models fail to describe the data. A smaller scale factor is generally required for the MC-KLN model, reflecting the fact that the ǫ 2 values from the MC-KLN model are on average larger than those from the Glauber model. Similar comparisons between v n and ǫ n for n = 3 and n = 4 are shown in figure 19 and figure 20, respectively. Agreement with the models is better than in the n = 2 case. However, this could simply reflect the fact that these distributions are dominated by Gaussian fluctuations, which have a universal shape. The shape differences between the v n and ǫ n distributions are also quantified in the right panels of figure 12 by comparing the values of σ vn / v n with the values of σ ǫn / ǫ n . Clearly, both the Glauber and MC-KLN models fail to describe the data consistently, in particular for n = 2, across the most of the measured centrality range. It should be noted that in the centrality intervals that correspond to the more peripheral collisions (e.g centrality interval 55-60%), the ǫ n values have been scaled down by a large factor and the sharp cutoffs of the ǫ n distributions are a natural consequence of the kinematic constraint ǫ n ≤ 1 (see eq. (1.2)). However, this behaviour also implies that v n is not proportional to ǫ n for large ǫ n values.

Summary
Measurements of the event-by-event harmonic flow coefficients v n for n = 2, 3 and 4 for various centrality intervals have been performed using 7 µb −1 of Pb+Pb collision data at √ s NN = 2.76 TeV collected by the ATLAS experiment at the LHC. The v 2 distributions approach that of a radial projection of a 2D Gaussian distribution centred around zero in the 0-2% centrality range, which is consistent with a scenario where fluctuations are the primary contribution to the overall shape (fluctuation-only scenario) for these most central collisions. Starting with the centrality interval 5-10%, the v 2 distributions differ significantly from this scenario, suggesting that they have a significant component associated with the average collision geometry in the reaction plane, v RP 2 . In contrast, the v 3 and v 4 distributions are consistent with a pure 2D Gaussian-fluctuation scenario (i.e., v RP n = 0) over most of the measured centrality range. However, a systematic deviation from this fluctuation-only scenario is observed for v 3 in mid-central collisions, the presence of a non-zero v RP 3 can be allowed. Simiarly, due to the limited range of the measured v 4 distribution, a non-zero v RP 4 on the order of δ v 4 can not be excluded by this analysis.
The v n distributions are also measured separately for charged particles with 0.5 < p T < 1 GeV and p T > 1 GeV. The shape of the unfolded distributions, when rescaled to the same v n , is found to be nearly the same for the two p T ranges. This finding suggests that the hydrodynamic response to the eccentricity of the initial geometry has little variation in this p T region. The ratios of the width to the mean, σ vn / v n , of these distributions are studied as a function of the average number of participating nucleons ( N part ) and p T . The values of σ v 2 / v 2 are observed to reach a minimum of 0.34 for N part ≈ 200, while the values of σ v 3 / v 3 and σ v 4 / v 4 are nearly independent of N part and are close to the value expected for a pure Gaussian-fluctuation scenario.
To understand further the role of average geometry and fluctuations for v 2 , the v 2 distributions have been fitted to a Bessel-Gaussian function to estimate the value of v RP 2 and the width of the fluctuation δ v 2 . In central collisions, where v RP 2 values are small and change rapidly with N part , narrow binning in centrality is necessary in order to obtain reliable estimates for these parameters. The values of δ v 2 /v RP 2 are found to decrease with N part and reach a minimum of 0.38 ± 0.02 at N part ≈ 200, but then increase and are greater than one in central collisions. Furthermore, a systematic deviation of the fit from the data is observed for centralities starting in the 15-20% centrality interval, and becoming more pronounced for the more peripheral collisions, suggesting significant non-Gaussian behaviour in the flow fluctuations for collisions with small N part . Multi-particle cumulant methods have been used to study such non-Gaussian behaviour. However, the v 2 coefficients for the four-, six-and eight-particle cumulants calculated from the v 2 distribution are found to agree with each other over the full centrality range, suggesting that the v 2 values from higher-order cumulants are not sensitive to this scale of deviation.
To elucidate the relation between the azimuthal anisotropy and underlying collision geometry, the measured v n distributions are compared with the eccentricity distributions of the initial geometry from the Glauber model and the MC-KLN model. Both models fail to describe the data consistently over most of the measured centrality range.  [35] ATLAS Collaboration, Measurement of the centrality dependence of the charged particle pseudorapidity distribution in lead-lead collisions at √ s N N = 2.76 TeV with the ATLAS detector, Phys. Lett. B710 (2012) 363, [arXiv:1108].

A Comprehensive performance and data plots
For reference, a more complete set of plots detailing the unfolding performance and crosschecks is presented. Figure 21 shows the distribution of the difference of the flow vectors ⇀ v obs 2 obtained for two half-IDs in peripheral collisions, calculated either without an η gap or with η gap = 2 between the two half-IDs. As mentioned in section. 4.1, due to the small number of charged particles in the events, this distribution is expected to be described by the Student's t-distribution instead of the Gaussian distribution. Figures 22 and 23 show the dependence on the priors for v 2 and v 4 in the 20-25% centrality interval. Despite the very different shapes of the underlying flow distributions, the convergence is robust and independent of the prior. Figure 24 compares the unfolding performance for the singleparticle and 2PC methods for several choices of η gap and demonstrates the consistency of the two methods for all η gap values considered. Figures 25 and 26 are similar to figure 11 -32 -but are obtained for other centrality intervals. They demonstrate that the shapes of the v n distributions are the same for 0.5 < p T < 1 GeV and p T > 1 GeV ranges in all centrality intervals. Figures 27 and 28 show a comparison between v n , v 2 n and v EP n in two different p T ranges; the trends are similar to those seen in figure 13.
Figures 29 compares the v RP 2 from the Bessel-Gaussian fit of the v 2 distribution with v calc 2 {2k} values for 0.5 < p T < 1 GeV and p T > 1 GeV ranges. The trends are similar to those seen in figure 16. However a slightly bigger deviation between v RP 2 and v calc 2 {2k} is observed in peripheral collisions for the p T > 1 GeV range. Figures 30 and 31 show the Bessel-Gaussian fits to the v obs 2 and v 2 distributions in several centrality intervals. Deviations from the fits are observed in both types of distributions with a similar trend, but the deviations are smaller for the v obs 2 distributions. Figure 32 compares the centrality dependence of v RP 2 and δ v 2 from the fits to distributions before and after unfolding. The effects of smearing by the response function increase the values for δ v 2 , but the values of v RP 2 only increase slightly. This is expected since the true v 2 distributions in the bulk region are close to a Bessel-Gaussian shape and the smearing due to the response function is nearly Gaussian (see eqs. 1.4 and 4.7). The slightly larger v RP 2 values for the v obs 2 distribution may reflect non-flow contributions in v obs 2 or deviation of the response function from a Gaussian (figure 21). Figure 33 shows the probability density distributions of v 2 for the full 0-5% centrality interval and the five individual 1% centrality intervals, together with the fits to the Bessel-Gaussian function. The deviation of the data from the Bessel-Gaussian description for the 0-5% centrality interval is due mainly to the rapid change of v RP 2 with centrality within this interval. This supports the need to use small centrality intervals for the most central collisions when calculating quantities such as v 2 , σ v 2 , v RP ) b (arbitrary normalization) for the 65-70% centrality interval with η gap = 0 (left panel) and η gap = 2 (right panel), together with fits to a Gaussian function and a Student's t-distribution. Bottom panels: the ratios of the data to the t-distribution fits. The distributions are wider for η gap = 2 since the average number of tracks used is smaller. 2PC v Figure 24. The unfolded distributions from 2PC method (top row), single-particle method (middle row) and the ratios of the two (bottom row) for different values of η gap in the 20-25% centrality interval and for n = 2 (left panels), n = 3 (middle panels) and n = 4 (right panels).
-35 - Rescaled to the same mean Figure 25. Top panels: The unfolded distributions for v n in the 0-5% centrality interval for charged particles in the p T > 0.5 GeV, 0.5 < p T < 1 GeV and p T > 1 GeV ranges. Bottom panels: same distributions but rescaled horizontally so the v n values match that for the p T > 0.5 GeV range. The shaded bands represent the systematic uncertainties on the v n -shape. Rescaled to the same mean Figure 26. Top panels: The unfolded distributions for v n in the 40-45% centrality interval for charged particles in the p T > 0.5 GeV, 0.5 < p T < 1 GeV and p T > 1 GeV ranges. Bottom panels: same distributions but rescaled horizontally so the v n values match that for the p T > 0.5 GeV range. The shaded bands represent the systematic uncertainties on the v n -shape.
-36 -  to v n . The shaded bands represent the systematic uncertainties. The dotted lines in bottom panels indicate v 2 n / v n = 1.13, the value expected for the radial projection of a 2D Gaussian distribution (eq. (1.8)).
-37 -  Figure 29. Comparison in 0.5 < p T < 1 GeV (left panel) and p T > 1 GeV (right panel) of the v RP 2 obtained from the Bessel-Gaussian fit of the v 2 distributions with the values for four-particle (v calc 2 {4}), six-particle (v calc 2 {6}) and eight-particle (v calc 2 {8}) cumulants calculated directly from the v 2 distributions via eq. (5.3). The bottom part of each panel shows the ratios of the cumulants to the fit results, with the error bars representing the total uncertainties.   Figure 33. The probability density distributions of v 2 for the 0-5% centrality interval and the five individual 1% centrality ranges, together with the fit to the Bessel-Gaussian function.