Measurement of longitudinal flow decorrelations in Pb+Pb

Measurements of longitudinal flow correlations are presented for charged particles in the pseudorapidity range |η| < 2.4 using 7 μb−1 and 470 μb−1 of Pb+Pb collisions at √ sNN = 2.76 and 5.02 TeV, respectively, recorded by the ATLAS detector at the LHC. It is found that the correlation between the harmonic flow coefficients vn measured in two separated η intervals does not factorise into the product of single-particle coefficients, and this breaking of factorisation, or flow decorrelation, increases linearly with the η separation between the intervals. The flow decorrelation is stronger at 2.76 TeV than at 5.02 TeV. Higher-order moments of the correlations are also measured, and the corresponding linear coefficients for the kth-moment of the vn are found to be proportional to k for v3, but not for v2. The decorrelation effect is separated into contributions from the magnitude of vn and the event-plane orientation, each as a function of η. These two contributions are found to be comparable. The longitudinal flow correlations are also measured between vn of different order in n. The decorrelations of v2 and v3 are found to be independent of each other, while the decorrelations of v4 and v5 are found to be driven by the nonlinear contribution from v2 2 and v2v3, respectively.


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 eventby-event (EbyE) density fluctuations in the initial state, the space-time evolution of the produced matter also fluctuates event by event. These fluctuations lead to correlations of particle multiplicity in momentum space in both the transverse and longitudinal directions with respect to the collision axis. Studies of particle correlations in the transverse plane have revealed strong harmonic modulation of the particle densities in the azimuthal angle: dN /dφ ∝ 1 + 2 ∞ n=1 v n cos n(φ − n ), where v n and n represent the magnitude and evente-mail: atlas.publications@cern.ch plane angle of the n th -order harmonic flow. The measurements of harmonic flow coefficients v n and their EbyE fluctuations, as well as the correlations between n of different order [3][4][5][6][7][8][9], have placed important constraints on the properties of the dense matter and on transverse density fluctuations in the initial state [10][11][12][13][14][15].
Most previous flow studies assumed that the initial condition and space-time evolution of the matter are boostinvariant in the longitudinal direction. Recent model studies of two-particle correlations as a function of pseudorapidity η revealed strong EbyE fluctuations of the flow magnitude and phase between two well-separated pseudorapidities, i.e. v n (η 1 ) = v n (η 2 ) (forward-backward or FB asymmetry) and n (η 1 ) = n (η 2 ) (event-plane twist) [16][17][18]. The CMS Collaboration proposed an observable based on the ratio of two correlations: the correlation between η and η ref and the correlation between −η and η ref . This ratio is sensitive to the correlation between η and −η [19]. The CMS results show that the longitudinal fluctuations lead to a linear decrease of the ratio with η, and the slope of the decrease shows a strong centrality dependence for elliptic flow v 2 but very weak dependences for v 3 and v 4 . This paper extends the CMS result by measuring several new observables based on multi-particle correlations in two or more η intervals [20]. These observables are sensitive to the EbyE fluctuations of the initial condition in the longitudinal direction. They are also sensitive to nonlinear mode-mixing effects, e.g. v 4 contains nonlinear contributions that are proportional to v 2 2 [8,9,[21][22][23]. Furthermore, the measurements are performed at two nucleon-nucleon centre-of-mass collision energies, √ s NN = 2.76 TeV and 5.02 TeV, to evaluate the √ s NN dependence of the longitudinal flow fluctuations. Recent model calculations predict an increase of longitudinal flow fluctuations at lower √ s NN [24]. Therefore, measurements of these observables at two collision energies can provide new insights into the initial condition along the longitudinal direction and should help in the development of full three-dimensional viscous hydrodynamic models. Using these new observables, this paper improves the study of the longitudinal dynamics of collective flow in three ways. Firstly, the CMS measurement, which is effectively the first moment of the correlation between v n in separate η intervals, is extended to the second and the third moments. Secondly, a correlation between four different η intervals is measured to estimate the contributions from the fluctuations of v n amplitudes as well as the contributions from fluctuations of n . Thirdly, correlations between harmonics of different order are also measured, e.g. between v 2 and v 4 in different η intervals, to investigate how mode-mixing effects evolve with rapidity. In this way, this paper presents a measurement of flow decorrelation involving v 2 , v 3 , v 4 and v 5 , using Pb+Pb collisions at √ s NN = 2.76 and 5.02 TeV.

Observables
This section gives a brief summary of the observables measured in this paper, further details can be found in Refs. [19,20,25]. The azimuthal anisotropy of the particle production in an event is conveniently described by harmonic flow vectors V n = v n e in n , 1 where v n and n are the magnitude and phase (or event plane), respectively. The V n are estimated from the observed per-particle normalised flow vector q n [5]: The sums run over all particles in a given η interval of the event, and φ i and w i are the azimuthal angle and the weight assigned to the ith particle, respectively. The weight accounts for detector non-uniformity and tracking inefficiency. The longitudinal flow fluctuations are studied using the correlation between the kth-moment of the nth-order flow vectors in two different η intervals, averaged over events in a given centrality interval, r n|n;k , for k = 1,2,3: where η ref is the reference pseudorapidity common to the numerator and the denominator, the subscript "n|n; k" denotes the kth-moment of the flow vectors of order n at η, combined with the kth moment of the conjugate of the flow vector of order n at η ref . The sine terms vanish in the last expression in Eq.
(2) because any observable must be an even function of n (−η) − n (η ref ). A schematic illustration of the choice of the η (|η| < 2.4) and η ref (4.0 < |η ref | < 4.9) to be discussed in Sect. 5, as well as the relations between different flow vectors, are shown in the left panel of Fig. 1. This observable is effectively a 2k-particle correlator between two subevents as defined in Ref. [28], and the particle multiplets containing duplicated particle indices are removed using the cumulant framework, with particle weights taken into account [20]. The observable measured by the CMS Collaboration [19] corresponds to k = 1, i.e. r n|n;1 . It should be noted that q n = 0 because the event plane changes randomly from event to event. Hence a direct study of the correlation between +η and −η via a quantity such as q n (+η)q * n (−η) /( q n (+η) q * n (−η) ) is not possible. One could also consider a quantity like q n (+η)q * n (−η) / q 2 n (η) q 2 n (−η) 1/2 , but the denominator would be affected by short-range correlations. Hence, it is preferable to work with quantities of the type used in Eq.
(2), which give a correlator sensitive to the flow decorrelation between η and −η through the reference flow vector q k n (η ref ). One important feature of Eq. (2) is that the detector effects at η ref are expected to cancel out to a great extent (see Sect. 5). To ensure a sizeable pseudorapidity gap between the flow vectors in both the numerator and denominator of Eq. (2), η ref is usually chosen to be at large pseudorapidity, e.g. η ref > 4 or η ref < −4, while the pseudorapidity of q n (−η) and q n (η) is usually chosen to be close to midrapidity, |η| < 2.4. If flow harmonics from multi-particle correlations factorise into single-particle flow harmonics, e.g.
, then it is expected that r n|n;k (η) = 1. Therefore, a value of r n|n;k (η) different from 1 implies a factorisation-breaking effect due to longitudinal flow fluctuations, and such an effect is generally referred to as "flow decorrelation".
Based on the CMS measurement [19] and arguments in Ref. [20], the observable r n|n;k (η) is expected to be approximately a linear function of η with a negative slope, and is sensitive to both the asymmetry in the magnitude of v n and the twist of the event-plane angles between η and −η: where F asy n;k and F twi n;k represent the contribution from FB v n asymmetry and event-plane twist, respectively. The r n|n;k results obtained in Ref. [19] were for k = 1 and n = 2, 3, 4. The measured F r n;1 show only a weak dependence on η ref for η ref > 3 or η ref < −3 at the LHC. Measuring r n|n;k for k > 1 provides new information on how the v n asymmetry and event-plane twist fluctuate event by event.
If the amount of decorrelation for the kth-moment of the flow vector is proportional to k, it can be shown that [20]: Deviations from Eq. (4) are sensitive to the detailed EbyE structure of the flow fluctuations in the longitudinal direction.
To estimate the separate contributions of the asymmetry and twist effects, a new observable involving correlations of flow vectors in four η intervals is used [20]: where the notation "2" in the subscript indicates that there are two q n and two q * n in the numerator and denominator. A schematic illustration of the relations between different flow vectors is shown in the right panel of Fig. 1. Since the effect of an asymmetry is the same in both the numerator and the denominator, this correlator is mainly sensitive to the event-plane twist effects: Therefore, the asymmetry and twist contributions can be estimated by combining Eqs.
Measurements of longitudinal flow fluctuations can also be extended to correlations between harmonics of different order: where the comma in the subscripts denotes the combination of q n of different order. If the longitudinal fluctuations for V 2 and V 3 are independent of each other, one would expect r 2,3|2,3 = r 2|2;1 r 3|3;1 [20]. On the other hand, r 2,2|4 and r 2,3|5 are sensitive to the η dependence of the correlations between v n and event planes of different order, for exam- Correlations between different orders have been measured previously at the LHC [8,9,23,29].
It is well established that the V 4 and V 5 in Pb+Pb collisions contain a linear contribution associated with initial geometry and mode-mixing contributions from lower-order harmonics due to nonlinear hydrodynamic response [8,9,14,21,22]: where the linear component V nL is driven by the corresponding eccentricity in the initial geometry [11]. If the linear component of v 4 and v 5 is uncorrelated with lower-order harmonics, i.e. V 2 2 V * 4L ∼ 0 and V 2 V 3 V * 5L ∼ 0, one expects [20]: Furthermore, using Eq. (10) the r n|n;1 correlators involving v 4 and v 5 can be approximated by: Therefore, both the linear and nonlinear components are important for r 4|4;1 and r 5|5;1 .

ATLAS detector and trigger
The ATLAS detector [30] provides nearly full solid-angle coverage of the collision point with tracking detectors, calorimeters, and muon chambers, and is well suited for measurements of multi-particle correlations over a large pseudorapidity range. 2 The measurements were performed using the inner detector (ID), minimum-bias trigger scintillators (MBTS), the forward calorimeters (FCal), and the zerodegree calorimeters (ZDC). The ID detects charged particles within |η| < 2.5 using a combination of silicon pixel detectors, silicon microstrip detectors (SCT), and a strawtube transition-radiation tracker (TRT), all immersed in a 2 T axial magnetic field [31]. An additional pixel layer, the "insertable B-layer" (IBL) [32] installed during the 2013-2015 shutdown between Run 1 and Run 2, is used in the 5.02 TeV measurements. The MBTS system detects charged particles over 2.1 |η| 3.9 using two hodoscopes of counters positioned at z = ±3.6 m. The FCal consists of three sampling layers, longitudinal in shower depth, and covers 3.2 < |η| < 4.9. The ZDC are positioned at ±140 m from the IP, detecting neutrons and photons with |η| > 8.3. This analysis uses approximately 7 and 470 µb −1 of Pb+Pb data at √ s NN = 2.76 and 5.02 TeV, respectively, recorded by the ATLAS experiment at the LHC. The 2.76 TeV data were collected in 2010, while the 5.02 TeV data were collected in 2015.
The ATLAS trigger system [33] consists of a level-1 (L1) trigger implemented using a combination of dedicated electronics and programmable logic, and a high-level trigger (HLT) implemented in general-purpose processors. The trigger requires signals in both ZDC or either of the two MBTS counters. The ZDC trigger thresholds on each side are set below the maximum corresponding to a single neutron. A timing requirement based on signals from each side of the MBTS was imposed to remove beam backgrounds. This trigger selected 7 and 22 µb −1 of minimum-bias Pb+Pb data at √ s NN = 2.76 TeV and √ s NN = 5.02 TeV, respectively. To increase the number of recorded events from very central Pb+Pb collisions, a dedicated L1 trigger was used in 2015 to select events requiring the total transverse energy ( E T ) in the FCal to be more than 4.54 TeV. This ultra-central trigger sampled 470 µb −1 of Pb+Pb collisions at 5.02 TeV and was fully efficient for collisions with centrality 0-0.1% (see Sect. 4).
2 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).

Event and track selection
The offline event selection requires a reconstructed vertex with its z position satisfying |Z vtx | < 100 mm. For the √ s NN = 2.76 TeV Pb+Pb data, the selection also requires a time difference | t| < 3 ns between signals in the MBTS trigger counters on either side of the nominal centre of ATLAS to suppress non-collision backgrounds. A coincidence between the ZDC signals at forward and backward pseudorapidity is required to reject a variety of background processes such as elastic collisions and non-collision backgrounds, while maintaining high efficiency for inelastic processes. The fraction of events containing more than one inelastic interaction (pile-up) is estimated to be less than 0.1% at both collision energies. The pile-up contribution is studied by exploiting the correlation between the transverse energy E T measured in the FCal or the 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 pileup, 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 event centrality [34] is characterised by the E T deposited in the FCal over the pseudorapidity range 3.2 < |η| < 4.9 using a calibration employing the electromagnetic calorimeters to set the energy scale [35]. 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. Thus the 0-5% centrality interval, for example, corresponds to the most central 5% of the events. The ultra-central trigger mentioned in Sect. 3 selects events in the 0-0.1% centrality interval with full efficiency. A Monte Carlo Glauber analysis [34,36] is used to estimate the average number of participating nucleons, N part , for each centrality interval. The systematic uncertainty in N part is less than 1% for centrality intervals in the range 0-20% and increases to 6% for centrality intervals in the range 70-80%. The Glauber model also provides a correspondence between the E T distribution and sampling fraction of the total inelastic Pb+Pb cross section, allowing centrality percentiles to be set. For this analysis, a selection of collisions corresponding to 0-70% centrality is used to avoid diffraction or other processes that contribute to very peripheral collisions. Following the convention used in heavy-ion analyses, the centrality dependence of the results in this paper is presented as a function of N part .
Charged-particle tracks and primary vertices [37] are reconstructed from hits in the ID. Tracks are required to have p T > 0.5 GeV and |η| < 2.4. For the 2.76 TeV data, tracks are required to have at least nine hits in the silicon detectors with no missing pixel hits and not more than one missing SCT hit, taking into account the presence of known dead modules. For the 5.02 TeV data, tracks are required to have at least two pixel hits, with the additional requirement of a hit in the first pixel layer when one is expected, at least eight SCT hits, and at most one missing hit in the SCT. In addition, for both datasets, the point of closest approach of the track is required to be within 1 mm of the primary vertex in both the transverse and longitudinal directions [38].
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 [39]. The generated particles in each event were rotated in azimuthal angle according to the procedure described in Ref.
[40] to produce harmonic flow consistent with previous ATLAS measurements [5,41]. The response of the detector was simulated using Geant 4 [42,43] and the resulting events are reconstructed with the same algorithms applied to the data. For the 5.02 TeV Pb+Pb data, the efficiency ranges from 75% at η ≈ 0 to about 50% for |η| > 2 for charged particles with p T > 0.8 GeV, falling by about 5% as p T is reduced to 0.5 GeV. The efficiency varies more strongly with η and event multiplicity. For p T > 0.8 GeV, it ranges from 75% at η ≈ 0 to 50% for |η| > 2 in peripheral collisions, while it ranges from 71% at η ≈ 0 to about 40% for |η| > 2 in central collisions. The tracking efficiency for the 2.76 TeV data has a similar dependence on p T and η. The efficiency is used in the particle weight, as described in Sect. 5. However, because the observables studied are ratios (see Sect. 2), uncertainties in detector and reconstruction efficiencies largely cancel. The rate of falsely reconstructed tracks ("fakes") is also estimated and found to be significant only at p T < 1 GeV in central collisions, where its percentage per-track ranges from 2% at |η| < 1 to 8% at the larger |η|. The fake rate drops rapidly for higher p T and towards more peripheral collisions. The fake rate is accounted for in the tracking efficiency correction following the procedure in Ref. [44].

Data analysis
Measurement of the longitudinal flow dynamics requires the calculation of the flow vector q n via Eq. (1) in the ID and the FCal. The flow vector from the FCal serves as the reference q n (η ref ), while the ID provides the flow vector as a function of pseudorapidity q n (η).
In order to account for detector inefficiencies and nonuniformity, a particle weight for the ith-particle in the ID for the flow vector from Eq. (1) is defined as: similar to the procedure in Ref.
[44]. The determination of track efficiency (η, p T ) is described in Sect. 4. The additional weight factor d ID (η, φ) corrects for variation of tracking efficiency or non-uniformity of detector acceptance as a function of η and φ. For a given η interval of 0.1, the distribu-tion in azimuthal bins, N (φ, η), is built up from reconstructed charged particles summed over all events. The weight factor is then obtained as . This "flattening" procedure removes most φ-dependent non-uniformity from track reconstruction, which is important for any azimuthal correlation analysis. Similarly, the weight in the FCal for the flow vector from Eq. (1) is defined as: where E T,i is the transverse energy measured in the ith tower in the FCal at η and φ. The azimuthal weight d FCal (η, φ) is calculated in narrow η intervals in a similar way to what is done for the ID. It ensures that the E T -weighted distribution, averaged over all events in a given centrality interval, is uniform in φ. The flow vectors q n (η) and q n (η ref ) are further corrected by an event-averaged offset: q n − q n evts [8].
The flow vectors obtained after these reweighting and offset procedures are used in the correlation analysis. The correlation quantities used in r n|n;k are calculated as: where subscripts "s" and "b" represent the correlator constructed from the same event ("signal") and from the mixedevent ("background"), respectively. The mixed-event quantity is constructed by combining q k n (η) from each event with q * k n (η ref ) obtained in other events with similar centrality (within 1%) and similar Z vtx (| Z vtx | < 5 mm). The q k n (η)q * k n (η ref ) b , which is typically more than two orders of magnitude smaller than the corresponding signal term, is subtracted to account for any residual detector non-uniformity effects that result from a correlation between different η ranges.
For correlators involving flow vectors in two different η ranges, mixed events are constructed from two different events. For example, the correlation for r 2,3|5 is calculated as: The mixed-event correlator is constructed by combining q 2 (η)q 3 (η) from one event with q * 5 (η ref ) obtained in another event with similar centrality (within 1%) and similar Z vtx (| Z vtx | < 5 mm). On the other hand, for correlators involving more than two different η ranges, mixed events are constructed from more than two different events, one for each unique η range. One such example is R n|n;2 , for which each mixed event is constructed from four different events with similar centrality and Z vtx . Table 1 The list of observables measured in this analysis

Observables
Pb+Pb datasets (Tev) r n|n;k for n = 2, 3, 4 and k = 1 2.76 and 5.02 R n|n;2 for n = 2, 3 2.76 and 5.02 r n|n;k for n = 5 and k = 1 5 . 0 2 r n|n;k for n = 2, 3 and k = 2,3 5.02 R n|n;2 for n = 4 5 . 0 2 r 2,2|4 , r 2,3|5 , r 2,3|2,3 5.02 Most correlators can be symmetrised. For example, in a symmetric system such as Pb+Pb collisions, the condition (2), the actual measured observable is: The symmetrisation procedure also allows further cancellation of possible differences between η and −η in the tracking efficiency or detector acceptance. Table 1 gives a summary of the set of correlators measured in this analysis. The analysis is performed in intervals of centrality and the results are presented as a function of η for |η| < 2.4. The main results are obtained using 5.02 TeV Pb+Pb data. The 2010 2.76 TeV Pb+Pb data are statistically limited, and are used only to obtain r n|n;1 and R n|n;2 to com-pare with results obtained from the 5.02 TeV data and study the dependence on collision energy.  19]. The behaviour in peripheral collisions is presumably due to increasing relative contributions from jets and dijets at higher p T and for peripheral collisions. Based on this, the measurements are performed using charged particles with 0.5 < p T < 3 GeV.

Systematic uncertainties
Since all observables are found to follow an approximately linear decrease with η, i.e. D(η) ≈ 1−cη for a given observable D(η) where c is a constant, the systematic uncertainty is presented as the relative uncertainty for 1 − D(η) at η = 1.2, the mid-point of the η range. The systematic uncertainties in this analysis arise from event mixing, track selection, and reconstruction efficiency. Most of the systematic uncertainties enter the analysis through the particle weights in Eqs. (14) and (15). In general, the uncertainties for r n|n;k increase with n and k, the uncertainties for R n|n;2 increase with n, and all uncertainties are larger in the most central and more peripheral collisions. For r 2,3|2,3 , r 2,2|4 and r 2,3|5 , the uncertainties are significantly larger than for the other correlators. Each source is discussed separately below. The effect of detector azimuthal non-uniformity is accounted for by the weight factor d(η, φ) in Eqs. (14) and (15). The effect of reweighting is studied by setting the weight to unity and repeating the analysis. The results are consistent with the default (weighted) results within statisti-cal uncertainties, so no additional systematic uncertainty is included. Possible residual detector effects for each observable are further removed by subtracting those obtained from mixed events as described in Sect. 5. Uncertainties due to the event-mixing procedure are estimated by varying the criteria for matching events in centrality and z vtx . The resulting uncertainty is in general found to be smaller than the statistical uncertainties. The event-mixing uncertainty for r 2|2;k and r 3|3;k is less than 1% for k = 1 and changes to about 0.4-8% for k = 2 and 0.6-10% for k = 3, while the uncertainty for r 4|4;1 and r 5|5;1 is in the range 1.5-3% and 5-13%, respectively. The uncertainty for R n|n;2 is 1.5-6% for n = 2 and 3-14% for n = 3. The uncertainties for r 2,3|2,3 , r 2,2|4 and r 2,2|5 are typically larger: 1-4%, 1.5-16% and 3-15%.
The systematic uncertainty associated with the track quality selections is estimated by tightening or loosening the requirements on transverse impact parameter |d 0 | and longitudinal impact parameter |z 0 sin θ | used to select tracks. In each case, the tracking efficiency is re-evaluated and the analysis is repeated. The difference is observed to be larger in the most central collisions where the flow signal is smaller and the influence of falsely reconstructed tracks is higher. The difference is observed to be in the range 0.2-12% for r 2|2;k and r 3|3;k , 1.1-2% for r 4|4;1 , 3-6% for r 5|5;1 , 0.5-13% for R n|n;2 , and 1-14% for r 2,3|2,3 , r 2,2|4 and r 2,2|5 .
From previous measurements [5,6,46], the v n signal has been shown to have a strong dependence on p T but relatively weak dependence on η. Therefore, a p T -dependent uncer-  tainty in the track reconstruction efficiency (η, p T ) could affect the measured longitudinal flow correlation, through the particle weights. The uncertainty in the track reconstruction efficiency is due to differences in the detector conditions and known differences in the material between data and simulations. The uncertainty in the efficiency varies between 1% and 4%, depending on η and p T [44]. The systematic uncertainty for each observable in Table 1 is evaluated by repeating the analysis with the tracking efficiency varied up and down by its corresponding uncertainty. For r n|n;k the uncertainties are in the range 0.1-2%, depending on n and k. For R n|n;2 the uncertainties are in the range 0.1-1%. For r 2,3|2,3 , r 2,2|4 and r 2,3|5 , the uncertainties are in the range 0.1-2%. Due to the finite energy resolution and energy scale uncertainty of the FCal, the q n (η ref ) calculated from the azimuthal distribution of the E T via Eqs. (1) and (15) differs from the true azimuthal distribution. However, since q n (η ref ) appears in both the numerator and the denominator of the correlators studied in this paper, most of the effects associated with the FCal E T response are expected to cancel out. Two cross-checks are also performed to study the influence of the FCal response. In the first cross-check, only the FCal towers with E T above the 50th percentile are used to calculate the q n (η ref ). The |q n (η ref )| value is different from the default analysis, but the values of the correlators are found to be consistent. In the second cross-check, HIJING events with imposed flow (see Sect. 4) are used to study the FCal response. The q n (η ref ) is calculated using both the generated E T and the reconstructed E T , and the resulting correlators are compared with each other. The results are found to be consistent. Accordingly, no additional systematic uncertainty is added for the FCal response.
The systematic uncertainties from the different sources described above are added in quadrature to give the total systematic uncertainty for each observable. They are listed in Tables 2, 3 and 4.

Results
The presentation of the results is structured as follows. Section 7.1 presents the results for r n|n;1 and R n|n;2 and the comparison between the two collision energies. Section 7.2 shows the results for r n|n;k for k > 1. The scaling relation from Eq. (4) is tested and the contributions from v n FB asymmetry and event-plane twist are estimated. Results for the mixedharmonic correlators, Eqs. (7)-(9), are presented in Sect. 7.3 and checked for compatibility with the hydrodynamical picture. The measurements are performed using charged particles with 0.5 < p T < 3 GeV, and the reference flow vector is calculated with 4.0 < |η ref | < 4.9. Most results are shown for the √ s NN = 5.02 TeV Pb+Pb dataset, which has better statistical precision. The results for the √ s NN = 2.76 TeV Pb+Pb dataset are shown only for r n|n;1 and R n|n;2 . 7.1 r n|n;1 and R n|n;2 at two collision energies Figure 6 shows r 2|2;1 in various centrality intervals at the two collision energies. The correlator shows a linear decrease with η, except in the most central collisions. The decreasing trend is weakest around the 20-30% centrality range, and is more pronounced in both more central and more peripheral collisions. This centrality dependence is the result of a strong centrality dependence of the v 2 associated with the average elliptic geometry [47]. The decreasing trend at √ s NN = 2.76 TeV is slightly stronger than that at √ s NN = 5.02 TeV, which is expected as the collision system becomes less boostinvariant at lower collision energy [24]. Figures 7 and 8 show the results for r 3|3;1 and r 4|4;1 , respectively, at the two collision energies. A linear decrease as a function of η is observed for both correlators, and the rate of the decrease is approximately independent of centrality. This centrality independence could be due to the fact that v 3 and v 4 are driven mainly by fluctuations in the initial state. The rate of the decrease is also observed to be slightly stronger at lower collision energy. The decreasing trend of r n|n;1 for n = 2-4 in Figs. 6, 7 and 8 indicates significant breakdown of the factorisation of two-particle flow harmonics into those between different η ranges. However, the size of the factorisation breakdown depends on the harmonic order n, collision centrality, and collision energy. The results have also been compared with those from the CMS Collaboration [19], with the η ref chosen to be 4.4 < |η ref | < 4.9 to match the CMS choice of η ref .
The two results agree very well with each other, and details are shown in the "Appendix". Figures 9 and 10 show R 2|2;2 and R 3|3;2 in several centrality intervals. Both observables follow a linear decrease with where the slope parameters are calculated as linear-regression coefficients, which characterise the average η-weighted deviation of r n|n;1 (η) and R n|n;2 (η) from unity. The sum runs over all data points. If r n|n;k and R n|n;2 are a linear function in η, the linear-regression coefficients are equivalent to a fit to Eq. (19). However, these coefficients are well defined even if the observables have significant nonlinear behaviour, which is the case for r 2|2;k and R 2|2;2 in the 0-20% centrality range. The extracted slope parameters F r n;1 and F R n;2 are plotted as a function of centrality in terms of N part , in Figs. 11 and 12, respectively. The values of F r 2;1 and F R 2;2 first decrease and then increase as a function of increasing N part . The larger values in central and peripheral collisions are related to the fact that v 2 is more dominated by the initial geometry fluctuations. The slopes for higher-order harmonics are significantly larger. As a function of N part , a slight decrease in F r 3;1 and F R 3;2 is observed for N part > 200, as well as an increase in √ s NN , as the rapidity profile of the initial state is more compressed due to smaller beam rapidity y beam at lower √ s NN . This energy dependence has been predicted for F r n;1 in hydrodynamic model calculations [24], and it is quantified in Fig. 13  √ s NN were entirely due to the change of y beam , then the correlators would be expected to follow a universal curve when they are rescaled by y beam , i. e. r n|n;k (η/y beam ) and R n|n;2 (η/y beam ) should not depend on √ s NN . In this case, the slopes parameters multiplified by the beam rapidity,F r n;1 ≡ F r n;1 y beam andF R n;2 ≡ F R n;2 y beam , should not depend on √ s NN . The beam rapidity is y beam = 7.92 and 8.52 for √ s NN = 2.76 and 5.02 TeV, respectively, which leads to a 7.5% reduction in the ratio. Figure 14 shows the ratio ofF r 2;1 values and ofF R 2;2 values at the two energies, and the weighted averages of the ratios calculated in the range 30 < N part < 400 are given in

Higher-order moments
The longitudinal correlations of higher-order moments of harmonic flow carry information about the EbyE flow fluctuations in pseudorapidity. In the simple model described in Ref. [20], the decrease in r n|n;k is expected to scale with k as given by Eq. (4).   Fig. 14 Centrality dependence of ratio ofF r n;1 ≡ F r n;1 y beam values (left panel) andF R n;2 ≡ F R n;2 y beam values (right panel) at 2.76 and 5.02 TeV. The lines indicate the average values in the range 30 < N part < 400, with the results and fit uncertainties given by Table 5. The error bars and shaded boxes are statistical and systematic uncertainties, respectively Table 5 Results of the fits to the ratio of F r n;1 , F R n;2 ,F r n;1 ≡ F r n;1 y beam andF R n;2 ≡ F R n;2 y beam values at the two energies in the range 30 < N part < 400 shown in Figs. 13  where v 2 is driven by the initial-state fluctuations. In other centrality intervals, where the average geometry is more important for v 2 , the r 2|2;k (k = 2 and 3) data show stronger decreases with η than r k 2|2;1 . A similar study is performed for third-order harmonics, and the results are shown in Fig. 16. The data follow approximately the scaling relation Eq. (4) in all centrality intervals.
To quantify the difference between r n|n;k and r k n|n;1 , the slopes (F r n;k ) of r n|n;k are calculated via Eqs. (19) and (20). The scaled quantities, F r n;k /k, are then compared with each other as a function of centrality in Fig. 17. For secondorder harmonics, the data show clearly that over most of the centrality range F r 2;3 /3 > F r 2;2 /2 > F r 2;1 , implying F r 2;k > k F r 2;1 . However, for the most central and most peripheral collisions the quantities approach each other. On the  Fig. 16 The r 3|3;k for k = 1-3 compared with r k 3|3;1 for k = 2-3 in various centrality intervals for Pb+Pb collisions at 5.02 TeV. The error bars and shaded boxes are statistical and systematic uncertainties, respectively. The data points for k = 2 or 3 in some centrality intervals are rebinned to reduce the uncertainty Fig. 17 The values of F r n;k /k for k = 1,2 and 3 for n = 2 (left panel) and n = 3 (right panel), respectively. The error bars and shaded boxes are statistical and systematic uncertainties, respectively. The widths of the centrality intervals are not fixed but are optimised to reduce the uncertainty other hand, a slightly opposite trend for the third-order harmonics, F r 3;3 /3 F r 3;2 /2 F r 3;1 , i.e. F r 3;k k F r 3;1 , is observed in mid-central collisions (150 < N part < 350). Figures 18 and 19 compare the r n|n;2 with R n|n;2 for n = 2 and n = 3, respectively. The decorrelation of R n|n;2 is significantly weaker than that for the r n|n;2 . This is because the R n|n;2 is mainly affected by the event-plane twist effects, while the r n|n;2 receives contributions from both FB asymmetry and event-plane twist [20].
Following the discussion in Sect. 2, Eqs.
(3) and (6), the measured F r n;2 and F R n;2 values can be used to estimate the separate contributions from FB asymmetry and event-plane twist, F asy n;2 and F twi n;2 , respectively, via the relation: F twi n;2 = F R n;2 , F asy n;2 = F r n;2 − F R n;2 .
The results are shown in Fig. 20. The contributions from the two components are similar to each other for n = 2, for which the harmonic flow arises primarily from the average collision shape, as well as for n = 3, for which the harmonic flow is driven mainly by fluctuations in the initial geometry.  Figure 22 also shows that the η dependence for r 4|4;1 is stronger than for r 2|2;2 in all centrality intervals, suggesting that the decorrelation effects are stronger for the linear component of v 4 than for the nonlinear component (see Eq. (12)).

Mixed-harmonics correlation
A similar study of the influence of the linear and nonlinear effects for v 5 was also performed, and results are shown in Fig. 23. The three observables r 2,3|2,3 , r 2,3|5 , and r 5|5;1 show similar values in all centrality intervals, albeit with large statistical uncertainties.
The results are summarised in Fig. 24, with each panel corresponding to the slopes of distributions in Figs. 21, 22, and 23, respectively. The only significant difference is seen between F 4|4;1 and F 2|2;2 or F 2,2|4 .

Summary
Measurements of longitudinal flow correlations for charged particles are presented in the pseudorapidity range |η| < 2.4 using 7 and 470 µb −1 of Pb+Pb data at √ s NN = 2.76 and 5.02 TeV, respectively, recorded by the ATLAS detector at the LHC. The factorisation of two-particle azimuthal correlations into single-particle flow harmonics v n is found to be broken, and the amount of factorisation breakdown increases approximately linearly as a function of the η separation between the two particles. The slope of this dependence is nearly independent of centrality and p T for n > 2. However, for n = 2 the effect is smallest in mid-central collisions and increases toward more central or more peripheral collisions, and in central collisions the effect also depends strongly on p T . Furthermore, the effect is found to be larger    Fig. 24 Comparison of the slopes of the correlators as a function of N part for three groups of correlators: r 2,3|2,3 and r 2|2;1 r 3|3;1 (for which the slope is F 2|2;1 + F 3|3;1 ) in Fig. 21 (left panel), r 2|2;2 , r 2,2|4 and r 4|4;1 in Fig. 22 (middle panel), and r 2,3|2,3 , r 2,3|5 and r 5|5;1 in Fig. 23 (right panel). The error bars and shaded boxes are statistical and systematic uncertainties, respectively at 2.76 than 5.02 TeV for all harmonics, which cannot be explained entirely by the change in the beam rapidity. The higher moments of the η-dependent flow correlations are also measured and the corresponding linear coefficients of the η dependence are extracted. The coefficient for the kth-moment of v n scales with k for n > 2, but scales faster than k for n = 2. The factorisation breakdown is separated into contributions from forward-backward asymmetry of the flow magnitude and event-plane twist, which are found to be comparable to each other.
The longitudinal flow correlations are also measured between harmonic flows of different order. The correlation of v 2 v 3 between two η ranges is found to factorise into the product of the correlation for v 2 and the correlation for v 3 , suggesting that the longitudinal fluctuations of v 2 and v 3 are independent of each other. The correlations between v 4 and v 2 2 suggest that the longitudinal fluctuations of v 4 have a significant nonlinear contribution from v 2 , i.e. v 4 ∝ v 2 2 . Similarly, the correlations between v 5 and v 2 v 3 suggest that the longitudinal fluctuations of v 5 are driven by the nonlinear contribution from v 2 v 3 , i.e. v 5 ∝ v 2 v 3 . The results presented in this paper provide new insights into the fluctuations and correlations of harmonic flow in the longitudinal direction, which can be used to improve full three-dimensional viscous hydrodynamic models.  Figure 29 compiles the results of r n|n;1 and R n|n;2 for 0-0.1% ultra-central collisions.  Fig. 25 The values of r 2|2;1 measured by ATLAS and by CMS [19]