Fluctuations of anisotropic flow in Pb + Pb collisions at √ s NN = 5 . 02 TeV with the ATLAS detector The ATLAS

Multi-particle azimuthal cumulants are measured as a function of centrality and transverse momentum using 470 μb−1 of Pb+Pb collisions atsNN = 5.02 TeV with the ATLAS detector at the LHC. These cumulants provide information on the event-by-event fluctuations of harmonic flow coefficients vn and correlated fluctuations between two harmonics vn and vm. For the first time, a non-zero four-particle cumulant is observed for dipolar flow, v1. The four-particle cumulants for elliptic flow, v2, and triangular flow, v3, exhibit a strong centrality dependence and change sign in ultra-central collisions. This sign change is consistent with significant non-Gaussian fluctuations in v2 and v3. The four-particle cumulant for quadrangular flow, v4, is found to change sign in mid-central collisions. Correlations between two harmonics are studied with threeand four-particle mixed-harmonic cumulants, which indicate an anticorrelation between v2 and v3, and a positive correlation between v2 and v4. These correlations decrease in strength towards central collisions and either approach zero or change sign in ultra-central collisions. To investigate the possible flow fluctuations arising from intrinsic centrality or volume fluctuations, the results are compared between two different event classes used for centrality definitions. In peripheral and mid-central collisions where the cumulant signals are large, only small differences are observed. In ultra-central collisions, the differences are much larger and transverse momentum dependent. These results provide new information to disentangle flow fluctuations from the initial and final states, as well as new insights on the influence of centrality fluctuations.


Fluctuations of anisotropic flow in Pb+Pb collisions at √ s NN = 5.02 TeV with the ATLAS detector
The ATLAS Collaboration Multi-particle azimuthal cumulants are measured as a function of centrality and transverse momentum using 470 µb −1 of Pb+Pb collisions at √ s NN = 5.02 TeV with the ATLAS detector at the LHC. These cumulants provide information on the event-by-event fluctuations of harmonic flow coefficients v n and correlated fluctuations between two harmonics v n and v m . For the first time, a non-zero four-particle cumulant is observed for dipolar flow, v 1 . The four-particle cumulants for elliptic flow, v 2 , and triangular flow, v 3 , exhibit a strong centrality dependence and change sign in ultra-central collisions. This sign change is consistent with significant non-Gaussian fluctuations in v 2 and v 3 . The four-particle cumulant for quadrangular flow, v 4 , is found to change sign in mid-central collisions. Correlations between two harmonics are studied with three-and four-particle mixed-harmonic cumulants, which indicate an anticorrelation between v 2 and v 3 , and a positive correlation between v 2 and v 4 . These correlations decrease in strength towards central collisions and either approach zero or change sign in ultra-central collisions. To investigate the possible flow fluctuations arising from intrinsic centrality or volume fluctuations, the results are compared between two different event classes used for centrality definitions. In peripheral and mid-central collisions where the cumulant signals are large, only small differences are observed. In ultra-central collisions, the differences are much larger and transverse momentum dependent. These results provide new information to disentangle flow fluctuations from the initial and final states, as well as new insights on the influence of centrality fluctuations.

Introduction
Heavy-ion collisions at RHIC and the LHC create hot, dense matter whose space-time evolution is well described by relativistic viscous hydrodynamics [1,2]. Owing to strong event-by-event density fluctuations in the initial state, the distributions of the final-state particles also fluctuate event by event. These fluctuations produce an effect in the azimuthal angle φ distribution of the final-state particles, characterized by a Fourier expansion dN dφ ∝ 1 + 2 ∑ ∞ n=1 v n cos n(φ − Φ n ), where v n and Φ n represent the magnitude and event-plane angle of the n th -order harmonic flow. These quantities also are conveniently represented by the 'flow vector' V n = v n e inΦn in each event. The V n value reflects the hydrodynamic response of the produced medium to the n th -order initial-state eccentricity vector [3,4], denoted by E n = n e inΨn . Model calculations show that V n is approximately proportional to E n in general for n = 2 and 3, and for n = 4 in the case of central collisions [3,5,6]. The measurements of v n and Φ n [7][8][9][10][11][12][13][14] place important constraints on the properties of the medium and on the density fluctuations in the initial state [4][5][6][15][16][17].
In order to disentangle the initial-and final-state effects, one needs detailed knowledge of the probability density distribution (or the event-by-event fluctuation) for single harmonics, p(v n ), and two harmonics, p(v n , v m ). These distributions are often studied through multi-particle azimuthal correlations within the cumulant framework [18][19][20][21]. In this framework, the moments of the p(v n ) distributions are measured by the 2k-particle cumulants, c n {2k}, for instance, c n {2} = ⟨v 2 n ⟩ and c n {4} = ⟨v 4 n ⟩ − 2 ⟨v 2 n ⟩ 2 which are then used to define flow harmonics v n {2k} such as v n {2} = (c n {2}) 1 2 and v n {4} = (−c n {4}) 1 4 . The four-particle cumulants c 2 {4} and c 3 {4} have been measured at RHIC and the LHC [22][23][24][25][26][27][28][29]. Most models of the initial state of A+A collisions predict a p(v n ) with shape that is close to Gaussian, and these models predict zero or negative values for c n {4} [30,31]. The values of c 2 {4} and c 3 {4} are found to be negative, except that c 2 {4} in very central Au+Au collisions at RHIC is positive [25]. Six-and eight-particle cumulants for v 2 have also been measured [22,26,32].
Assuming that the scaling between V n and E n is exactly linear, then p(v n ) and p(v n , v m ) should be the same as p( n ) and p( n , m ) up to a global rescaling factor. In order to isolate the initial eccentricity fluctuations, it was proposed in Ref. [35] to measure the ratios of two cumulants of different order, for instance nc n {4} ≡ c n {4} (c n {2}) 2 4 . Similar cumulant ratios can be constructed for symmetric and asymmetric cumulants such as nsc n,m {4} ≡ sc n,m {4} (⟨v 2 n ⟩ ⟨v 2 m ⟩) and nac n {3} = ac n {3} (⟨v 4 n ⟩ ⟨v 2 2n ⟩) 1 2 . In addition, hydrodynamic model calculations suggest strong p Tdependent fluctuations of v n and Φ n even in a single event [36,37]. Such final-state intra-event flow fluctuations may change the shape of p(v n ) or p(v n , v m ) in a p T -dependent way and can be quantified by comparing cumulant ratios using particles from different p T ranges.
In heavy-ion collisions, v n coefficients are often calculated for events with similar centrality, defined by the particle multiplicity in a fixed pseudorapidity range, which is also referred to as the reference multiplicity. The event ensemble, selected using a given reference multiplicity, is referred to as a reference event class. Due to fluctuations in the particle production process, the true centrality for events with the same reference multiplicity still fluctuates from event to event. Since the v n values vary with centrality, the fluctuations of centrality can lead to additional fluctuations of v n and change the underlying p(v n ) and p(v n , v m ) distributions [38]. Consequently, the cumulants c n {2k}, sc n,m {4}, and ac n {3} could be affected by the centrality resolution effects that are associated with the definition of the reference event class. Such centrality fluctuations, also known as volume fluctuations, have been shown to contribute significantly to event-by-event fluctuations of conserved quantities, especially in ultra-central collisions [39][40][41]. Recently, the centrality fluctuations were found to affect flow fluctuations as indicated by the sign change of c 2 {4} measured in ultra-central collisions [38]. A detailed study of c n {2k}, sc n,m {4} and ac n {3} for different choices of the reference event class helps clarify the meaning of centrality and provides insight into the sources of particle production in heavy-ion collisions. In this paper, two reference event-class definitions are used to study the influence of centrality fluctuations on flow cumulants. The total transverse energy in the forward pseudorapidity range 3.2 < η < 4.9 is taken as the default definition and a second definition uses the number of reconstructed charged particles in the mid-rapidity range η < 2.5.  [27,33] in order to quantify the influence of non-flow correlations such as resonance decays and jets. Results using the two reference event-class definitions are compared in order to understand the role of centrality fluctuations and to probe the particle production mechanism which directly influences the size of centrality fluctuations.
The paper is organized as follows. Sections 2 and 3 describe the detector, trigger and datasets, as well as event and track selections. The mathematical framework for the multi-particle cumulants and the list of cumulant observables are provided in Section 4. The correlation analysis and systematic uncertainties are described in Sections 5 and 6, respectively. Section 7 first presents the results for various cumulant observables and then investigates the role of centrality fluctuations by making a detailed comparison of the cumulants calculated using two reference event classes. A summary is given in Section 8.

ATLAS detector and trigger
The ATLAS detector [42] provides nearly full solid-angle coverage with tracking detectors, calorimeters, and muon chambers, and is well suited for measurements of multi-particle correlations over a large pseudorapidity range.1 The measurements are performed using the inner detector (ID), the forward calorimeters (FCal), and the zero-degree calorimeters (ZDC). The ID detects charged particles within η < 2.5 using a combination of silicon pixel detectors, silicon microstrip detectors (SCT), and a straw-tube transition-radiation tracker, all immersed in a 2 T axial magnetic field [43]. An additional pixel layer, the 'insertable B-layer' [44,45], was installed during the 2013-2015 shutdown between Run 1 and Run 2, and is used in the present analysis. The FCal consists of three sampling layers, longitudinal in shower depth, and covers 3.2 < η < 4.9. The ZDC, positioned at ±140 m from the IP, detects neutrons and photons with η > 8.3.
The ATLAS trigger system [46] consists of a level-1 (L1) trigger implemented using a combination of dedicated electronics and programmable logic, and a high-level trigger (HLT), which uses software algorithms similar to those applied in the offline event reconstruction. Events for this analysis were selected by two types of trigger. The minimum-bias trigger required either a scalar sum, over the whole calorimeter system, of transverse energy ΣE tot T greater than 0.05 TeV or the presence of at least one neutron on both sides of the ZDC in coincidence with a track identified by the HLT. This trigger selected 22 µb −1 of Pb+Pb data. The number of recorded events from very central Pb+Pb collisions was increased by using a dedicated trigger selecting on the ΣE tot T at L1 and ΣE T , the total transverse energy in the FCal, at HLT. The combined trigger selects events with ΣE T larger than one of the three threshold values: 4.21 TeV, 4.37 TeV and 4.54 TeV. This ultra-central trigger has a very sharp turn-on as a function of ΣE T and for these thresholds was fully efficient for the 1.3%, 0.5% and 0.1% of events with the highest transverse 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). energy in the FCal. The trigger collected 52 µb −1 , 140 µb −1 and 470 µb −1 of Pb+Pb collisions for the three thresholds, respectively.
In the offline data analysis, events from the minimum-bias and ultra-central triggers are combined as a function of ΣE T by applying an event-by-event weight calculated as the ratio of the number of minimum-bias events to the total number of events. This procedure ensures that the weighted distribution as a function of ΣE T for the combined dataset follows the distribution of the minimum-bias events, and the results measured as a function of ΣE T or centrality (see Section 3) are not biased in their ΣE T or centrality values.

Event and track selection
The analysis uses approximately 470 µb −1 of √ s NN = 5.02 TeV Pb+Pb data collected in 2015. The offline event selection requires a reconstructed primary vertex with a z position satisfying z vtx < 100 mm. A coincidence between the ZDC signals at forward and backward pseudorapidity rejects a variety of background processes such as elastic collisions and non-collision backgrounds, while maintaining high efficiency for inelastic processes. The contribution from events containing more than one inelastic interaction (pile-up) is studied by exploiting the correlation between the transverse energy, ΣE T , measured in the FCal or the estimated number of neutrons N n in the ZDC and the number of tracks associated with a primary vertex N rec ch . Since the distribution of ΣE T or N n in events with pile-up is broader than that for the events without pile-up, pile-up events are suppressed by rejecting events with an abnormally large ΣE T or N n as a function of N rec ch . The remaining pile-up contribution after this procedure is estimated to be less than 0.1%.
The Pb+Pb event centrality [47] is characterized by the ΣE T deposited in the FCal over the pseudorapidity range 3.2 < η < 4.9. The FCal ΣE T distribution is divided into a set of centrality intervals. A centrality interval refers to a percentile range, starting at 0% relative to the most central collisions at the largest ΣE T value. Thus the 0-5% centrality interval, for example, corresponds to the most central 5% of the events. The ultra-central trigger mentioned in Section 2 enhances the number of events in the 0-1.3%, 0-0.5% and 0-0.1% centrality intervals with full efficiency for the three L1 ΣE T thresholds, respectively. Centrality percentiles are set by using a Monte Carlo Glauber analysis [47,48] to provide a correspondence between the ΣE T distribution and the sampling fraction of the total inelastic Pb+Pb cross section.
Charged-particle tracks [49] are reconstructed from hits in the ID and are then used to construct the primary vertices. Tracks are required to have p T > 0.5 GeV and η < 2.5. They are required to have at least one pixel hit, with the additional requirement of a hit in the first pixel layer when one is expected, and at least six SCT hits. In addition, each track must have transverse and longitudinal impact parameters relative to the primary vertex which satisfy d 0 < 1.5 mm and z 0 sin θ < 1.5 mm, respectively [50].
The efficiency (p T , η) of the track reconstruction and track selection criteria is evaluated using Pb+Pb Monte Carlo events produced with the HIJING event generator [51]. The generated particles in each event are rotated in azimuthal angle according to the procedure described in Ref. [52] in order to produce a harmonic flow that is consistent with the previous ATLAS measurements [9,50]. The response of the detector is simulated using G 4 [53,54] and the resulting events are reconstructed with the same algorithms as are applied to the data. For peripheral collisions, the efficiency ranges from 75% at η ≈ 0 to about 50% for η > 2 for charged particles with p T > 0.8 GeV. The efficiency falls by about 5% for a p T of 0.5 GeV. The efficiency in central collisions ranges from 71% at η ≈ 0 to about 40% for η > 2 for charged particles with p T > 0.8 GeV, falling by about 8% for a p T of 0.5 GeV. The rate of falsely reconstructed tracks ('fake' tracks) is also estimated and found to be significant only at p T < 1 GeV in central collisions where it ranges from 2% for η < 1 to 8% at larger η . The fake-track rate drops rapidly for higher p T and for more peripheral collisions. The fake-track rate is accounted for in the tracking efficiency correction following the procedure in Ref. [22].

Observables
Both the standard cumulant method [18] and the three-subevent cumulant method [27,33,55] are used to calculate the cumulants c n {4}, sc n,m {4} and ac n {3}. However, only the standard method is used to calculate the six-particle cumulants c n {6}.

Cumulants in the standard method
The standard cumulant method calculates the 2k-particle (k = 1,2...) cumulants c n {2k} from the 2k-particle azimuthal correlations ⟨{2k} n ⟩, which are calculated for each event as [19,20] where '⟨⟩' denotes a single-event average over all pairs, quadruplets or sextuplets, respectively. The averages from Eq. (1) can be expressed in terms of per-particle normalized flow vectors q n;l with l = 1, 2... in each event [19] where the sum runs over all particles in the event and w j is a weight assigned to the j th particle. This weight is constructed to correct for both detector non-uniformity and tracking inefficiency as explained in Section 5.
The multi-particle cumulants are obtained from the azimuthal correlations using where '⟪⟫' represents a weighted average of ⟨{2k} n ⟩ over an event ensemble with similar ΣE T or N rec ch . In the absence of non-flow correlations, the c n {2k} values are related to the moments of the p(v n ) distribution by the expression given in the last part of each equation chain. In particular, the higher moments of p(v n ) can be obtained by combining the cumulants of different order, for example ⟨v which extends the standard definition [20] of v n {2k} to regions where c n {4} > 0 and c n {6} < 0.
If the fluctuation of the event-by-event flow-vector V n = v n e inΦn is described in the plane transverse to the beam by a two-dimensional Gaussian function 2 given by [11]. The parameter δ n is the width of the Gaussian function and v 0 n is related to the average geometry of the overlap region. However, if the shape of p(v n ) has significant non-Gaussian fluctuations at large v n , both c n {4} and c n {6} may change sign, giving negative values for v n {4} and v n {6} [56].
The four-particle symmetric cumulants sc n,m {4} and three-particle asymmetric cumulants ac n {3} are related to multi-particle azimuthal correlations for two flow harmonics of different order by [20,55] The first average is over all distinct quadruplets, triplets or pairs in one event to obtain ⟨{4} n,m ⟩, ⟨{3} n ⟩, ⟨{2} n ⟩ and ⟨{2} m ⟩, and the second average is over an event ensemble with the same ΣE T or N rec ch to obtain sc n,m {4} and ac n {3}. In the absence of non-flow correlations, sc n,m {4} and ac n {3} are related to the correlation between v n and v m or between v n and v 2n , respectively: This analysis measures three types of cumulants defined by Eq.

Cumulants in the subevent method
In the 'standard' cumulant method described so far, all the k-particle multiplets involved in ⟨{k} n ⟩ and ⟨{k} n,m ⟩ are selected using charged tracks that are in the entire ID acceptance of η < 2.5. In order to further suppress the non-flow correlations that typically involve particles emitted within a localized region in η, the charged tracks are grouped into three subevents, labelled a, b and c, that each cover a unique η range [33]: . Various subevent cumulants are then constructed by correlating particles between different subevents: where The statistical precision is enhanced by interchanging the η range for subevent a with that for subevent b or c which results in three independent measurements for each of c n {4}, sc n,m {4} and ac n {3}. They are averaged to obtain the final result.
It is well known that the values of c n {2} and v n {2} calculated using the standard cumulant method have a significant contribution from non-flow effects [57]. Therefore, in this analysis, they are measured using the two-subevent method following the expressions used in previous publications [58]: This definition ensures that the non-flow correlations in v n {2} are greatly reduced by requiring a minimum pseudorapidity gap of 1.67 between subevents a and c. For k-particle cumulants with k > 2, the standard method is used as the default since they are less influenced by non-flow correlations, and this assumption is additionally verified with the three-subevent method.

Normalized cumulants and cumulant ratios
Any quantity which is linearly proportional to v n has the same cumulants, up to a global factor. Therefore the shapes of p(v n ) and p(v n , v m ) can be more directly probed using the ratio of the cumulants [59,60]: where the two-particle cumulants c n {2} in the denominator of these equations are calculated from subevents a and c using Eq. (7). If v n is exactly proportional to n , the normalized cumulants defined above would be the same as the normalized cumulants calculated from eccentricities in the initial state, i.e. nc n {2k} = nc n {2k, }, nsc n,m {4} = nsc n,m {4, } and nac n {3} = nac n {3, }. In practice, final-state effects, such as p T -dependent fluctuations of v n and Φ n [36,37], hydrodynamic noise [61] and non-linear mode-mixing between harmonics of different order [3,62] can break this equality. Therefore, studying the p T dependence of these normalized cumulants can help in understanding the influence of dynamical effects from the final state.
The nc n {4} and nc n {6} cumulants defined above contain the same information as the previously proposed The

Data analysis
The cumulants are calculated in three steps following examples from Refs. [27,55] using the standard and subevent methods. Since these steps are the same for c n {2k}, sc n,m {4} and ac n {3}, they are explained using c n {2k} as an example.
In the first step, the multi-particle correlators ⟨{2k} n ⟩ are calculated for each event from particles in one of four p T ranges: 0.5 < p T < 5 GeV, 1.0 < p T < 5 GeV, 1.5 < p T < 5 GeV, and 2 < p T < 5 GeV.
In the second step, the correlators ⟨{2k} n ⟩ are averaged over an event ensemble, defined as events in either a narrow interval of ΣE T (0.002 TeV) or a narrow interval of N rec ch (track bin width is 1) taken as the number of reconstructed charged particles in the range 0.5 < p T < 5 GeV. The c n {2k} values are then calculated separately for these two types of reference event classes, denoted by c n {2k, ΣE T } and c n {2k, N rec ch }, respectively. In order to obtain statistically significant results, in the final step the c n {2k} values from several neighbouring ΣE T or N rec ch intervals are combined, weighted by the number of events in each interval. The p T dependence of the cumulants is studied by simultaneously varying the p T range for all particles in each 2k-multiplet in the cumulant analysis. This approach is different from previous studies where the p T range of only one particle in the multiplet is varied [20,22,26,63,64].
The left panel of Figure 1 shows the correlation between ΣE T and N rec ch . The two quantities have an approximately linear correlation, but events with the same ΣE T have significant fluctuations in N rec ch and vice versa. Due to these relative fluctuations, the reference event class based on N rec ch may have centrality fluctuations that differ from those of the reference event class based on ΣE T , even if both are matched to have the same ⟨ΣE T ⟩ or the same ⟨N rec ch ⟩. The correlation between ΣE T and N rec ch is studied using events divided into narrow intervals in either ΣE T or N rec ch . The mean and root-mean-square values of the N rec ch (ΣE T ) distributions are calculated for each ΣE T (N rec ch ) interval, and the results are shown in the middle and right panels of Figure 1, respectively. A linear relation is observed between ⟨N rec ch ⟩ and ΣE T over the full ΣE T range, while a significant non-linear relation is observed between ⟨ΣE T ⟩ and N rec ch at large N rec ch . This latter behaviour suggests that, in ultra-central collisions, ΣE T retains sensitivity to the ⟨N rec ch ⟩ of the events, while N rec ch has relatively poorer sensitivity to the ⟨ΣE T ⟩ of the events. This implies that the true centrality is more smeared for events with the same N rec ch than for events with the same ΣE T .
Since v n changes with centrality, any centrality fluctuations could lead to additional fluctuation of v n , and subsequently to a change in the flow cumulants. Indeed, previous ATLAS studies [27,55,58] have shown that the c n {2k} values depend on the definition of the reference event class used for averaging. A comparison of the results based on these two reference event classes can shed light on the details of flow fluctuations and how they are affected by centrality fluctuations.   Figure 1. The inserted panels show the local first-order derivatives of the one-dimensional ΣE T or N rec ch distributions in the most central collisions. The derivative for the ΣE T distribution is relatively independent of ΣE T up to 4.1 TeV and then decreases and reaches a local minimum at around 4.4 TeV. The derivative for the N rec ch distribution is mostly flat up to 2800 and then decreases and reaches a local minimum at around 3100. The locations where the derivative starts to depart from a constant are defined as the knee of the ΣE T or N rec ch distribution and is given by (ΣE T ) knee = 4.1 TeV and (N rec ch ) knee = 2800. Events with ΣE T > (ΣE T ) knee correspond to the top 1.9% centrality and events with N rec ch > (N rec ch ) knee correspond to top 2.7% centrality when mapped to the equivalent ⟨ΣE T ⟩. The knees mark the locations where multiplicity distributions start to decrease sharply and the underlying centrality fluctuations are expected to deviate significantly from a Gaussian distribution [38,41]. The knee values are important in discussing the trends of cumulants in ultra-central collisions in Section 7.3.
The particle weights used in Eq. (2) that account for detector inefficiencies and non-uniformity are defined as [58] where (η, p T ) is the efficiency for reconstructing charged particles. The additional weight factor d(φ, η) accounts for non-uniformities in the efficiency as a function of φ in each η range. All reconstructed charged particles with p T > 0.5 GeV are entered into a two-dimensional histogram N(φ, η), and the weight factor is then obtained as where ⟨N(η)⟩ is the track density averaged over φ in the given η interval. This procedure corrects most of the φ-dependent non-uniformity that results from track reconstruction [58].

Systematic uncertainties
The systematic uncertainties of the measurements presented in this paper are evaluated by varying different aspects of the analysis and comparing c n {2k}, sc 2,3 {4}, sc 2,4 {4} and ac 2 {3} with their baseline values. The main sources of systematic uncertainty are track selection, the track reconstruction efficiency, the pile-up contribution, and differences between data and Monte Carlo simulation. The uncertainties are generally small when the absolute values of the cumulants are large. The relative uncertainties are larger in central or very peripheral collisions where the signal is small. The uncertainties also decrease rapidly with increasing p T , due to a larger flow signal at higher p T and are typically less than a few percent for p T > 1 GeV. Therefore, the following discussion focuses mainly on the results obtained for charged particles in the 0.5 < p T < 5 GeV range. The systematic uncertainties are also found to be similar between the standard method and the three-subevent method.
The systematic uncertainty associated with track selection is evaluated by applying more restrictive requirements. The requirement on d 0 and z 0 sin θ is changed to be less than 1.0 mm instead of the nominal value of 1.5 mm. The numbers of pixel and SCT hits required are also increased, to two and eight respectively, to further reduce the fake-track rates. The uncertainties are less than 2% for c n {2}, less than 3% for c 2 {4}, c 2 {6} and c 3 {4}, less than 5% for c 1 {4} and c 4 {4}, and are in the range of 1-5% for sc 2,3 {4}, sc 2,4 {4} and ac 2 {3}.
Previous measurements [9] show that the v n signal has a strong dependence on p T but a relatively weak dependence on η. Therefore, a p T -dependent uncertainty in the track reconstruction efficiency (η, p T ) could affect the measured cumulants through the particle weights in Eqs. (2) and (13). The uncertainty of (η, p T ) arises from differences in the detector conditions and known differences in the material between data and simulations. This uncertainty varies between 1% and 4%, depending on η and p T [22]. Its impact on cumulants is evaluated by repeating the analysis with the tracking efficiency varied up and down by its corresponding uncertainty. The impact on cumulants is in the range of 1-5% for c n {2}, 0.5-12% for c n {4} and c n {6}, and in the range of 2-8% for sc n,m {4} and ac 2 {3}.
Pile-up events are suppressed by exploiting the correlation, discussed in Section 3, between ΣE T measured in the FCal and the number of neutrons N n in the ZDC. In the ultra-central collisions, where the pile-up fraction is the largest, the residual pile-up is estimated to be less than 0.1%. The impact of the pile-up is evaluated by tightening and relaxing pile-up rejection criteria, and the resulting variation is included in the systematic uncertainty. The uncertainty is in the range of 0.1-1% for all cumulants.
The analysis procedure is also validated through Monte Carlo studies by comparing the observables calculated with generated particles with those obtained from reconstructed particles, using the same analysis chain and correction procedure as for data. In the low p T region, where tracking performance suffers from low efficiency and high fake-track rates, systematic differences are observed between the cumulants calculated at the generator level and at the reconstruction level. These differences are included as part of the systematic uncertainty. They amount to 0.1-3% in mid-central and peripheral collisions and up to 10% in the most central collisions.
The systematic uncertainties from different sources are added in quadrature to determine the total systematic uncertainties. These uncertainties for two-particle cumulants are in the range of 1-5% for c 2 {2}, 2-7% for c 3 {2} and 4-9% for c 4 {2}. The uncertainties for normalized cumulants, nc n {4}, nc 2 {6}, nsc 2,3 {4}, nsc 2,4 {4} and nac 2 {3}, are calculated separately for each source of systematic uncertainty discussed above, and are similar to the baseline results. Most of the systematic uncertainties cancel out in these ratios. In mid-central and peripheral collisions, the total uncertainties are in the range of 1-5% depending on the observables. However, the total uncertainties are larger in ultra-central collisions, reaching as high as 10% for nc 2 {6} and nac 2 {3}.

Results
The results for various cumulant observables are presented in Sections 7.1 and 7.2. The cumulants are calculated using the reference event class based on ΣE T and with the procedure discussed in Section 5. The results are presented as a function of centrality calculated from ΣE T . Section 7.1 discusses the cumulants related to single harmonics: c n {2k, ΣE T }, v n {2k, ΣE T }, and nc n {2k, ΣE T }. Section 7.2 presents correlations between two flow harmonics: nsc 2,3 {4, ΣE T }, nsc 2,4 {4, ΣE T } and nac 2 {3, ΣE T }. The results are shown for four p T ranges: 0.5 < p T < 5 GeV, 1.0 < p T < 5 GeV, 1.5 < p T < 5 GeV, and 2 < p T < 5 GeV. The default results are obtained using the standard cumulant method and are compared with those obtained using the three-subevent cumulant method. The comparisons are shown only if significant differences are observed; otherwise, they are included in Appendix B.  Figure 3 shows the v n {2} values for n = 2, 3, 4 for charged particles in several p T ranges, calculated for the event class based on FCal ΣE T and then plotted as a function of centrality. The v n {2} values are obtained from two-particle cumulants with a pseudorapidity gap according to Eq. (7). For all p T ranges, v 2 {2} first increases and then decreases toward central collisions, reflecting the typical centrality dependence behaviour of the eccentricity 2 [57]. The magnitude of v 2 {2} also increases strongly with p T . The centrality and p T dependences of v 3 {2} and v 4 {2} are similar, but the tendency to decrease from mid-central toward central collisions is less pronounced.   Figure 4 shows the centrality dependence of normalized four-particle cumulants nc 2 {4}, nc 3 {4}, and nc 4 {4} in four p T ranges using the standard method (top row) and the three-subevent method (bottom row). The advantage of using nc n {4} instead of c n {4} is that the p T dependence of v n , seen in Figure 3, is largely cancelled out and that nc n {4} directly reflects the shape of the p(v n ) distributions [11]. Overall, the results based on the three-subevent method behave similarly to those from the standard cumulant method, implying that the influence of non-flow correlations is small. Therefore, the remaining discussion is focused on the standard method in the top row.  [22,23,65] to be driven by the centrality dependence of the four-particle cumulants for 2 and 3 , respectively. The normalized cumulants still show some residual dependence on p T . Namely, the nc 2 {4} values are smaller for the higher-p T particles, while the values of nc 3 {4} are larger for the higher p T range. Furthermore, the values of nc 2 {4} are also observed to change sign in ultra-central collisions and the pattern of these sign changes also has significant p T dependence. The observed behaviour of nc n {4} in ultra-central collisions is closely related to centrality fluctuations and is discussed further in Section 7.3.  The nc 4 {4} values, as shown in the right panels of Figure 4, are negative in central collisions but change sign around a centrality range of 25-30%. The centrality value at which the sign change occurs shifts towards more peripheral collisions as the p T of the particles increases. It is well established that V 4 in Pb+Pb collisions contains a linear contribution associated with the initial geometry and a mode-mixing contribution from lower-order harmonics due to the non-linear hydrodynamic response [3,12,13,16,62],

Flow cumulants for p(v n )
where the linear component V 4L is driven by the corresponding eccentricity 4 in the initial geometry [5], and χ 2 is a constant. Previous measurements [12,13] show that the V 4L term dominates in central collisions, while the V 2 2 term dominates in more peripheral collisions. Therefore, the sign change of nc 4 {4} could reflect an interplay between these two contributions [66]. In central collisions, nc 4 {4} is dominated by a negative contribution from p(v 4L ), while in peripheral collisions nc 4 {4} is dominated by a positive contribution from p(v 2 2 ). The change of the crossing point with p T suggests that the relative contribution from these two sources is also a function of p T .
If the v n value is driven only by n , then p(v n ) should have the same shape as p( n ). On the other hand, the significant p T dependence of nc n {4} in Figure 4 suggests that the shape of p(v n ) also changes with p T . Such p T -dependent behaviour implies that the eccentricity fluctuations in the initial state are not the only source for flow fluctuations. Dynamical fluctuations in the momentum space in the initial or final state may also change p(v n ).   From the measured nc 2 {6} and nc 2 {4}, the ratio of the six-particle cumulant to the fourth-particle cumulant, v 2 {6} v 2 {4}, can be obtained. The results are shown in Figure 7. For the Gaussian fluctuation model in Eq. (5), this ratio is expected to be one. The apparent deviation of the ratio from one suggests non-Gaussianity of p(v 2 ) over a broad centrality range. The results for different p T ranges are close to each other, but nevertheless show systematic-and centrality-dependent differences. In general, the results from higher p T are larger in central collisions and smaller in peripheral collisions than those from lower p T .
The multi-particle correlations are also calculated to obtain cumulants for the dipolar flow, v 1 . Figure 8 shows the centrality dependence of c 1 {4} in several p T ranges, which is obtained from the reference event class based on ΣE T . In the hydrodynamic picture, c 1 {4} is sensitive to event-by-event fluctuations  of the dipolar eccentricity 1 associated with initial-state geometry [5]. This measurement has a large uncertainty because both ⟪{4} 1 ⟫ and ⟪{2} 1 ⟫ in Eq. (3) contain a significant contribution from global momentum-conservation effects [9,67]. This contribution cancels out for c 1 {4} but leads to a large statistical uncertainty. A negative c 1 {4} for p T > 1.5 GeV is observed in both the standard and three-subevent cumulant methods, which reflects the event-by-event fluctuations of the dipolar eccentricity.
Previously, ATLAS measured v 1 using the two-particle correlation method in Pb+Pb collisions at √ s NN = 2.76 TeV where an explicit procedure was employed to subtract the global momentum-conservation   effects [9]. The v 1 {2} values was observed to be negative at low p T , change sign at p T ≈ 1.2 GeV and increase quickly for higher p T . Therefore, a c 1 {4} signal is expected to be larger and easier to measure at higher p T . Figure 9 shows the v 1 {4} values calculated from c 1 {4} for the two highest p T ranges: 1.5 < p T < 5 GeV and 2 < p T < 5 GeV. The v 1 {4} values increase both with p T and in more peripheral collisions, and are in the range of 0.02-0.04 for 2 < p T < 5 GeV.

Flow cumulants for p(v n , v m )
The correlation between flow harmonics of different order is studied using the four-particle normalized symmetric cumulant nsc 2,3 {4} and nsc 2,4 {4}, and the three-particle normalized asymmetric cumulant nac 3 {3}. Figure 10 shows the centrality dependence of nsc 2,3 {4} in several p T ranges which probes the correlation between the v 2 and v 3 . The nsc 2,3 {4} is negative in most of the centrality range, indicating an anti-correlation between the v 2 and v 3 . This anti-correlation has been observed in previous studies based on the same observable [14] and using an event-shape engineering technique [13]. The strength of the anti-correlation has significant p T dependence. For higher-p T particles, the anti-correlation is stronger in peripheral collisions and weaker in central collisions. In the ultra-central collisions, nsc 2,3 {4} changes sign and becomes positive. This positive correlation is related to centrality fluctuations and is discussed further in Section 7.3. The behaviour of the overall centrality and p T dependence is also found to be similar between the standard cumulant method and the three-subevent cumulant method. This suggests that these features are not caused by non-flow correlations.  Figure 11 shows the centrality dependence of nsc 2,4 {4} in several p T ranges which probes the correlation between v 2 and v 4 . The nsc 2,4 {4} value is positive over the entire centrality range, indicating a positive correlation between v 2 and v 4 . The signal is very small in central collisions but increases rapidly towards peripheral collisions. The correlations are similar among different p T ranges in central collisions but are slightly weaker for higher-p T particles in mid-central collisions. This behaviour is also predicted by hydrodynamic models [6,68]. Compared with the three-subevent method, the nsc 2,4 {4} values from the standard method have better statistical precision but slightly higher values in peripheral collisions, indicating that the non-flow effects may become significant for events beyond 60% centrality. Figure 12 shows the centrality dependence of ac 2 {3} in several p T ranges which also probes the correlation between v 2 and v 4 . The ac 2 {3} value is positive over the entire centrality range. The correlation is weak in the central collisions, increases rapidly as the centrality approaches about 20-30% and then increases slowly toward more peripheral collisions. The correlation patterns among different p T ranges are similar in central collisions but are slightly weaker for higher-p T particles in mid-central collisions. Compared with results obtained from the three-subevent method, the results from the standard method are slightly larger in peripheral collisions, indicating that non-flow fluctuations may contribute for events with more than 60% centrality. The similar p T and centrality dependences for nsc 2,4 {4} and nac 2 {3} are related to the non-linear mode-mixing effects between v 2 and v 4 described by Eq. (14) [59].

Dependence on reference event class and the role of centrality fluctuations
This section presents the ⟨ΣE T ⟩ or ⟨N rec ch ⟩ dependence of various cumulants for the two reference event classes. Section 5 describes how the role of centrality fluctuations associated with the reference event class used in the calculation of the cumulants can be understood by extracting the results for each observable in narrow ranges of ΣE T and N rec ch . These results are presented as a function of ⟨ΣE T ⟩ (ΣE T ) knee and ⟨N rec ch ⟩ (N rec ch ) knee , where (ΣE T ) knee = 4.1 TeV and (N rec ch ) knee = 2800 are the knee values of the ΣE T and N rec ch distributions shown in Figure 2. It should be noted that c n {2k, ΣE T } (and other observables as well) as a function of ⟨ΣE T ⟩ (ΣE T ) knee contains the same information as the centrality dependence of c n {2k, ΣE T } shown in two previous sections. However, x-axes based on ⟨ΣE T ⟩ (ΣE T ) knee and ⟨N rec ch ⟩ (N rec ch ) knee more naturally characterize the size of the overlap region in Pb+Pb collisions and allow a more detailed visualization of the ultra-central region, where the impacts of centrality fluctuations is strongest.

Two-particle cumulants
The top panels of Figure 13 show v n {2, ΣE T } as a function of ⟨ΣE T ⟩. The v n {2, ΣE T } values tend to reflect the same centrality and p T dependence behaviour already shown in Figure 3. In ultra-central collisions, the v n {2, ΣE T } values are nearly constant. Similar trends are also observed for v n {2, N rec ch } which are shown in the bottom panels of Figure 13 as a function of ⟨N rec ch ⟩. These results suggest that the underlying initial geometry, in terms of ⟨ 2 n ⟩, is quite similar between the two reference event classes.
In order to quantify differences between the two reference event classes, v n {2, N rec ch } is mapped to a ⟨ΣE T ⟩ dependence and v n {2, ΣE T } is mapped to a ⟨N rec ch ⟩ dependence. The ratio v n {2, N rec ch } v n {2, ΣE T } is then calculated at a given ⟨ΣE T ⟩ or at a given ⟨N rec ch ⟩. The top row of Figure 14 shows The bottom row of Figure 14 shows the same ratio, v n {2, N rec ch } v n {2, ΣE T }, but instead as a function of ⟨N rec ch ⟩. Compared with the upper row of Figure 14, the ratio for v 2 shows a larger deviation from unity which reaches 7% in ultra-central collisions. Smaller, but significant differences are also observed for v 3 and v 4 in ultra-central collisions. This is probably because v n {2, N rec ch } has even more contributions from less central events than v n {2, ΣE T } when both are matched to the same ⟨N rec ch ⟩ instead of the same ⟨ΣE T ⟩. This is consistent with the hypothesis in which N rec ch has poorer centrality resolution and therefore larger centrality fluctuations than ΣE T , when mapped to the same average event activity in the final state.
Due to the steep decrease of the ΣE T and N rec ch distributions in the ultra-central region, the centrality fluctuations and the shapes of the p( n ) and p(v n ) distributions are expected to exhibit a significant departure from a Gaussian shape [38,39]. The flow cumulants with four or more particles are more sensitive to a non-Gaussian shape of p(v n ) than the two-particle cumulants. Therefore, they are expected to exhibit larger differences between the two reference event classes. This is the topic of the next section.

Multi-particle cumulants
The top panels of Figure 15 show nc n {4, ΣE T } as a function of ⟨ΣE T ⟩. This figure contains the same information as the results shown in Figure 4, except for a change in the scale of the x-axis which shows the central region in more detail. The nc 2 {4, ΣE T } value changes sign for ⟨ΣE T ⟩ ≳ (ΣE T ) knee , where it first increases, reaches a maximum and then decreases to close to zero. The value of the maximum also increases with the p T of the particles. The nc 3 {4, ΣE T } value is negative and approaches zero in ultra-central collisions and only changes sign for the highest p T range used in this analysis. The nc 4 {4, ΣE T } value changes from positive in peripheral collisions to negative in mid-central collisions, reaches a minimum and then turns back and approaches zero in the ultra-central collisions.
The bottom panels of Figure 15 show nc n {4, N rec ch } as a function of ⟨N rec ch ⟩. The overall ⟨N rec ch ⟩ and p T -dependent trends are similar to those in the top panels. However, the maximum of nc 2 {4, N rec ch } is more than a factor of two larger, and nc 3 {4, N rec ch } shows a clear sign change for the two highest p T ranges used in this analysis. Furthermore, nc 4 {4, N rec ch } shows a local maximum in ultra-central collisions, a feature absent for nc 4 {4, ΣE T }.
If V n ∝ E n is valid, then the shape of p(v n ) should be the same as the shape of p( n ) and nc n {4} = nc n {4, } [35,38]. The c n {4, } values can be estimated from a simple Glauber model framework using participating nucleons in the overlap region. The c n {4, } value is found to be always negative when the reference event class is defined using the number of participating nucleons N part or the impact parameter of the collisions [65]. However, a positive nc n {4, } is observed in ultra-central collisions when the reference event class is defined using the final-state particle multiplicity [38,69]. Due to multiplicity smearing, events with the same final-state multiplicity can have different N part , and therefore different n . The positive nc n {4, } reflects the non-Gaussian shape of p( n ) due to the smearing in N part for events with the same final-state multiplicity. The larger values of nc n {4, N rec ch } in comparison with nc n {4, ΣE T } in ultra-central collisions could be due to stronger multiplicity smearing for nc n {4, N rec ch }. Figure 16 compares nc n {4, ΣE T } and nc n {4, N rec ch } as a function of ⟨ΣE T ⟩ obtained for 1.5 < p T < 5 GeV. In both cases, the normalized cumulants for v 2 and v 3 show significant differences between the two reference event classes, while the difference is smaller for v 4 . The values of nc n {4, N rec ch } for n = 2 and 3 are significantly larger than those for nc n {4, ΣE T } over a broad centrality range, not only limited to the ultra-central collisions. This implies that the influence of centrality fluctuations on flow fluctuations is potentially important even in mid-central collisions.
The left two panels of Figure 17 show the six-particle normalized cumulants for v 2 obtained using the two reference event classes, nc 2 {6, ΣE T } and nc 2 {6, N rec ch }, respectively. The nc 2 {6} values are positive in most of the centrality range but decrease to zero at around ⟨ΣE T ⟩ = (ΣE T ) knee or ⟨N rec ch ⟩ = (N rec ch ) knee and stay close to zero above that. The right panel of Figure 17 Figure 18 shows v 2 {6} v 2 {4} obtained using the event class based on N rec ch but then mapped onto ⟨ΣE T ⟩. The differences between the results for the various p T ranges are larger for most of the centrality range, which again implies that the centrality fluctuations influence the ratios between multi-particle cumulants over a broad centrality range.

Multi-particle mixed-harmonic cumulants
{4} exhibit a small but significant p T dependence, suggesting flow fluctuations may also arise directly in the momentum space through the initial-state correlations or final-state interactions.
This paper also presents a detailed measurement of the four-particle symmetric cumulants nsc 2,3 {4} and nsc 2,4 {4} and the three-particle asymmetric cumulant nac 2 {3}. The symmetric cumulants probe the correlation between the magnitudes of two flow harmonics, while the asymmetric cumulant is sensitive to correlations involving both the magnitude and phase of flow. Over most of the centrality range, nsc 2,3 {4} is found to be negative, reflecting an anti-correlation between v 2 and v 3 . The nsc 2,4 {4} and nac 2 {3} values are found to be positive, and their dependence on centrality is consistent with non-linear mode-mixing effects between v 2 and v 4 .
In experimental measurements, the flow cumulants are always calculated for events with similar activity. However, for a given activity measure, fluctuations in the particle production process lead to irreducible centrality fluctuations, also known as volume fluctuations. Since v n changes with centrality, centrality fluctuations lead to an additional fluctuation of v n , and consequently a change in the flow cumulants. In order to study the influence of centrality fluctuations, cumulant observables are calculated for two reference event classes with different centrality resolution: the total transverse energy in 3.2 < η < 4.9, and number of reconstructed charged particles with η < 2.5 and 0.5 < p T < 5 GeV. In ultra-central collisions, the cumulants nc 2 {4}, nc 3 {4}, and nsc 2,3 {4} are observed to change sign, indicating a significant influence of centrality fluctuations on the multi-particle cumulants of p(v 2 ), p(v 3 ) and p(v 2 , v 3 ). The sign change patterns are more pronounced for the event class based on ⟨N rec ch ⟩, consistent with larger centrality fluctuations. The differences between the two event classes are found to extend, with decreasing magnitude, to mid-central collisions, which may suggest that the centrality fluctuations influence the flow fluctuations over a broad centrality range. The sign-change patterns are found to be more pronounced at higher p T , which may indicate that the flow fluctuations have significant p T dependence. Such p T dependence cannot be explained by considering only fluctuations in the initial geometry.
These results provide comprehensive information about the nature of flow fluctuations and the contributions coming from both the initial state and the final state. They also shed light on the influence of centrality fluctuations on flow fluctuations, especially in the ultra-central collisions, which can help to clarify the meaning of centrality and provide insights into the sources of particle production in heavy-ion collisions.

Appendix
A Flow harmonics v n {2k} from 2k-particle correlations  show the flow coefficients from the four-particle cumulants v 2 {4}, v 3 {4} and v 4 {4}, respectively. Figure 24 shows the elliptic flow coefficient from the six-particle cumulant, v 2 {6}. They are all obtained from Eq. (4) and are shown as a function of centrality, ΣE T and N rec ch . The apparent discontinuities correspond to the locations where the corresponding c n {2k} changes sign.

B Comparison between standard method and three-subevent method
This appendix shows a comparison between the standard cumulant method and the three-subevent method for various cumulant observables.       The ATLAS Collaboration