1 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 (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: d\(N/{\text {d}}\phi \propto 1+2\sum _{n=1}^{\infty }v_{n}\cos n(\phi -\Phi _{n})\), where \(v_n\) and \(\Phi _n\) represent the magnitude and event-plane angle of the \(n^{\mathrm {th}}\)-order harmonic flow. The measurements of harmonic flow coefficients \(v_n\) and their EbyE fluctuations, as well as the correlations between \(\Phi _{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 boost-invariant in the longitudinal direction. Recent model studies of two-particle correlations as a function of pseudorapidity \(\eta \) revealed strong EbyE fluctuations of the flow magnitude and phase between two well-separated pseudorapidities, i.e. \(v_n(\eta _1)\ne v_n(\eta _2)\) (forward-backward or FB asymmetry) and \(\Phi _{n}(\eta _1)\ne \Phi _{n}(\eta _2)\) (event-plane twist) [16,17,18]. The CMS Collaboration proposed an observable based on the ratio of two correlations: the correlation between \(\eta \) and \(\eta _{\mathrm {ref}}\) and the correlation between \(-\eta \) and \(\eta _{\mathrm {ref}}\). This ratio is sensitive to the correlation between \(\eta \) and \(-\eta \) [19]. The CMS results show that the longitudinal fluctuations lead to a linear decrease of the ratio with \(\eta \), 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 \(\eta \) 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, \(\sqrt{s_{\text {NN}}}=2.76\) TeV and 5.02 TeV, to evaluate the \(\sqrt{s_{\text {NN}}}\) dependence of the longitudinal flow fluctuations. Recent model calculations predict an increase of longitudinal flow fluctuations at lower \(\sqrt{s_{\text {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 \(\eta \) intervals, is extended to the second and the third moments. Secondly, a correlation between four different \(\eta \) intervals is measured to estimate the contributions from the fluctuations of \(v_n\) amplitudes as well as the contributions from fluctuations of \(\Phi _{n}\). Thirdly, correlations between harmonics of different order are also measured, e.g. between \(v_2\) and \(v_4\) in different \(\eta \) 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 \(\sqrt{s_{\text {NN}}}=2.76\) and 5.02 TeV.

2 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 \( {\varvec{V}}_n = v_n e^{{\text {i}}n\Phi _n}\) Footnote 1, where \(v_n\) and \(\Phi _n\) are the magnitude and phase (or event plane), respectively. The \( {\varvec{V}}_n\) are estimated from the observed per-particle normalised flow vector \({\varvec{q}}_n\) [5]:

$$\begin{aligned} {\varvec{q}}_{n} \equiv \frac{{\sum _{i}} {w_{i}}e^{{\text{i}}n{\phi_{i}}}}{{\sum _i}{w_i}}. \end{aligned}$$
(1)

The sums run over all particles in a given \(\eta \) interval of the event, and \(\phi _i\) and \(w_i\) are the azimuthal angle and the weight assigned to the \(i{\text {th}}\) particle, respectively. The weight accounts for detector non-uniformity and tracking inefficiency.

The longitudinal flow fluctuations are studied using the correlation between the \(k{\text {th}}\)-moment of the \(n{\text {th}}\)-order flow vectors in two different \(\eta \) intervals, averaged over events in a given centrality interval, \(r_{n|n;k}\), for \(k=1\),2,3:

$$\begin{aligned} r_{n|n;k}(\eta )= & {} \frac{\left\langle {\varvec{q}}_n^k (-\eta ) {\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle }{\left\langle {\varvec{q}}_n^k (\eta ){\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle }\nonumber \\= & {} \frac{\left\langle \left[ v_n(-\eta ) v_n(\eta _{\mathrm {ref}})\right] ^k \cos kn(\Phi _n(-\eta )-\Phi _n(\eta _{\mathrm {ref}}))\right\rangle }{\left\langle \left[ v_n(\eta ) v_n(\eta _{\mathrm {ref}})\right] ^k \cos kn(\Phi _n(\eta )-\Phi _n(\eta _{\mathrm {ref}}))\right\rangle }\;, \end{aligned}$$
(2)

where \(\eta _{\mathrm {ref}}\) is the reference pseudorapidity common to the numerator and the denominator, the subscript “n|nk” denotes the \(k{\text {th}}\)-moment of the flow vectors of order n at \(\eta \), combined with the \(k{\text {th}}\) moment of the conjugate of the flow vector of order n at \(\eta _{\mathrm {ref}}\). The sine terms vanish in the last expression in Eq. (2) because any observable must be an even function of \(\Phi _n(-\eta )-\Phi _n(\eta _{\mathrm {ref}})\). A schematic illustration of the choice of the \(\eta \) (\(|\eta |<2.4\)) and \(\eta _{\mathrm {ref}}\) (\(4.0<|\eta _{\mathrm {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].

Fig. 1
figure 1

Schematic illustration of the procedure for constructing the corrlators \(r_{n|n;k}(\eta )\) Eq. (2) (left panel) and \(R_{n|n;2}(\eta )\) Eq. (5) (right panel). The acceptance coverages for the ATLAS tracker used for \(\eta \) and reference detector used for \(\eta _{\mathrm {ref}}\) are discussed in Sect. 5

The observable measured by the CMS Collaboration [19] corresponds to \(k=1\), i.e. \(r_{n|n;1}\) . It should be noted that \(\left\langle {\varvec{q}}_n\right\rangle =0\) because the event plane changes randomly from event to event. Hence a direct study of the correlation between \(+\eta \) and \(-\eta \) via a quantity such as \(\left\langle {\varvec{q}}_n (+\eta ) {\varvec{q}}_n^{*}(-\eta )\right\rangle /(\left\langle {\varvec{q}}_n (+\eta )\right\rangle \left\langle {\varvec{q}}_n^{*} (-\eta )\right\rangle )\) is not possible. One could also consider a quantity like \(\left\langle {\varvec{q}}_n (+\eta ) {\varvec{q}}_n^{*}(-\eta )\right\rangle /\left( \left\langle q_n^2 (\eta )\right\rangle \left\langle q_n^2 (-\eta )\right\rangle \right) ^{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 \(\eta \) and \(-\eta \) through the reference flow vector \({\varvec{q}}_n^{k}(\eta _{\mathrm {ref}})\).

One important feature of Eq. (2) is that the detector effects at \(\eta _{\mathrm {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), \(\eta _{\mathrm {ref}}\) is usually chosen to be at large pseudorapidity, e.g. \(\eta _{\mathrm {ref}}>4\) or \(\eta _{\mathrm {ref}}<-4\), while the pseudorapidity of \({\varvec{q}}_n(-\eta )\) and \({\varvec{q}}_n(\eta )\) is usually chosen to be close to mid-rapidity, \(|\eta |<2.4\). If flow harmonics from multi-particle correlations factorise into single-particle flow harmonics, e.g. \(\left\langle {\varvec{V}}_n^k (\eta ) {\varvec{V}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle ^2=\left\langle v_n^{2k} (\eta )\right\rangle \left\langle v_n^{2k}(\eta _{\mathrm {ref}})\right\rangle \), then it is expected that \(r_{n|n;k}(\eta )=1\). Therefore, a value of \(r_{n|n;k}(\eta )\) 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}(\eta )\) is expected to be approximately a linear function of \(\eta \) 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 \(\eta \) and \(-\eta \):

$$\begin{aligned} r_{n|n;k}(\eta ) \approx 1-2F_{n;k}^{\text {r}}\eta ,\;\; F_{n;k}^{\text {r}}= F_{n;k}^{\text {asy}}+F_{n;k}^{\text {twi}}, \end{aligned}$$
(3)

where \(F_{n;k}^{\text {asy}}\) and \(F_{n;k}^{\text {twi}}\) 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_{n;1}^{\text {r}}\) show only a weak dependence on \(\eta _{\mathrm {ref}}\) for \(\eta _{\mathrm {ref}}>3\) or \(\eta _{\mathrm {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 \(k{\text {th}}\)-moment of the flow vector is proportional to k, it can be shown that [20]:

$$\begin{aligned} r_{n|n;k}\approx r_{n|n;1}^k,\;\;F_{n;k}^{\text {r}}\approx k F_{n;1}^{\text {r}}. \end{aligned}$$
(4)

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 \(\eta \) intervals is used [20]:

$$\begin{aligned} R_{n|n;2}(\eta )= & {} \frac{\left\langle {\varvec{q}}_n(-\eta _{\mathrm {ref}}){\varvec{q}}_n^{*} (\eta ) {\varvec{q}}_n (-\eta ) {\varvec{q}}_n^*(\eta _{\mathrm {ref}})\right\rangle }{\left\langle {\varvec{q}}_n(-\eta _{\mathrm {ref}}){\varvec{q}}_n^* (-\eta ){\varvec{q}}_n (\eta ) {\varvec{q}}_n^*(\eta _{\mathrm {ref}})\right\rangle }\nonumber \\= & {} \frac{\left\langle v_n(-\eta _{\mathrm {ref}})v_n(-\eta )v_n(\eta )v_n(\eta _{\mathrm {ref}})\cos n\left[ \Phi _n(-\eta _{\mathrm {ref}})-\Phi _n(\eta _{\mathrm {ref}})+(\Phi _n(-\eta )-\Phi _n(\eta ))\right] \right\rangle }{\left\langle v_n(-\eta _{\mathrm {ref}})v_n(-\eta )v_n(\eta )v_n(\eta _{\mathrm {ref}})\cos n\left[ \Phi _n(-\eta _{\mathrm {ref}})-\Phi _n(\eta _{\mathrm {ref}})-(\Phi _n(-\eta )-\Phi _n(\eta ))\right] \right\rangle }, \end{aligned}$$
(5)

where the notation “2” in the subscript indicates that there are two \({\varvec{q}}_n\) and two \({\varvec{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:

$$\begin{aligned} R_{n|n;2} (\eta ) \approx 1-2F_{n;2}^{\text {R}}\eta , F_{n;2}^{\text {R}}= F_{n;2}^{\text {twi}}. \end{aligned}$$
(6)

Therefore, the asymmetry and twist contributions can be estimated by combining Eqs. (3) and (6).

Measurements of longitudinal flow fluctuations can also be extended to correlations between harmonics of different order:

$$\begin{aligned} r_{2,3|2,3}(\eta )= & {} \frac{\left\langle {\varvec{q}}_2 (-\eta ){\varvec{q}}_2^{*}(\eta _{\mathrm {ref}}){\varvec{q}}_3 (-\eta ){\varvec{q}}_3^{*}(\eta _{\mathrm {ref}})\right\rangle }{\left\langle {\varvec{q}}_2 (\eta ){\varvec{q}}_2^{*}(\eta _{\mathrm {ref}}){\varvec{q}}_3 (\eta ){\varvec{q}}_3^{*}(\eta _{\mathrm {ref}})\right\rangle }, \end{aligned}$$
(7)
$$\begin{aligned} r_{2,2|4}(\eta )= & {} \frac{\left\langle {\varvec{q}}_2^2 (-\eta ) {\varvec{q}}_4^{*}(\eta _{\mathrm {ref}})\right\rangle +\left\langle {\varvec{q}}_2^2 (\eta _{\mathrm {ref}}) {\varvec{q}}_4^{*}(-\eta )\right\rangle }{\left\langle {\varvec{q}}_2^2 (\eta ){\varvec{q}}_4^{*}(\eta _{\mathrm {ref}})\right\rangle +\left\langle {\varvec{q}}_2^2 (\eta _{\mathrm {ref}}){\varvec{q}}_4^{*}(\eta )\right\rangle }, \end{aligned}$$
(8)
$$\begin{aligned} r_{2,3|5}(\eta )= & {} \frac{\left\langle {\varvec{q}}_2(-\eta ){\varvec{q}}_3(-\eta ) {\varvec{q}}_5^{*}(\eta _{\mathrm {ref}}) \right\rangle +\left\langle {\varvec{q}}_2(\eta _{\mathrm {ref}}){\varvec{q}}_3(\eta _{\mathrm {ref}}) {\varvec{q}}_5^{*}(-\eta ) \right\rangle }{\left\langle {\varvec{q}}_2 (\eta ){\varvec{q}}_3 (\eta ){\varvec{q}}_5^{*}(\eta _{\mathrm {ref}})\right\rangle +\left\langle {\varvec{q}}_2 (\eta _{\mathrm {ref}}){\varvec{q}}_3 (\eta _{\mathrm {ref}}){\varvec{q}}_5^{*}(\eta )\right\rangle }, \end{aligned}$$
(9)

where the comma in the subscripts denotes the combination of \({\varvec{q}}_n\) of different order. If the longitudinal fluctuations for \({\varvec{V}}_2\) and \({\varvec{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 \(\eta \) dependence of the correlations between \(v_n\) and event planes of different order, for example \(\left\langle {\varvec{q}}_2^2 (-\eta ) {\varvec{q}}_4^{*}(\eta _{\mathrm {ref}})\right\rangle = \left\langle v_2^2(-\eta ) v_4(\eta _{\mathrm {ref}}) \cos 4(\Phi _2(-\eta )-\Phi _4(\eta _{\mathrm {ref}}))\right\rangle \). Correlations between different orders have been measured previously at the LHC [8, 9, 23, 29].

It is well established that the \({\varvec{V}}_4\) and \({\varvec{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]:

$$\begin{aligned} {\varvec{V}}_4 = {\varvec{V}}_{4\text {L}} +\chi _{4} {\varvec{V}}_2^2,\;\;{\varvec{V}}_5 = {\varvec{V}}_{5\text {L}}+\chi _{5} {\varvec{V}}_2{\varvec{V}}_3, \end{aligned}$$
(10)

where the linear component \({\varvec{V}}_{n\text {L}}\) 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. \({\varvec{V}}_2^2{\varvec{V}}_{4\text {L}}^{*}\sim 0\) and \({\varvec{V}}_2{\varvec{V}}_3{\varvec{V}}_{5\text {L}}^{*}\sim 0\), one expects [20]:

$$\begin{aligned} r_{2,2|4}\approx r_{2|2;2},\;\;r_{2,3|5}\approx r_{2,3|2,3}. \end{aligned}$$
(11)

Furthermore, using Eq. (10) the \(r_{n|n;1}\) correlators involving \(v_4\) and \(v_5\) can be approximated by:

$$\begin{aligned} r_{4|4;1}(\eta )\approx & {} \frac{\left\langle {\varvec{V}}_{4\text {L}}(-\eta ){\varvec{V}}_{4\text {L}}^{*}(\eta _{\mathrm {ref}})\right\rangle +\chi _{4}^2\left\langle {\varvec{V}}_2^2 (-\eta ) {\varvec{V}}_2^{*2}(\eta _{\mathrm {ref}})\right\rangle }{\left\langle {\varvec{V}}_{4\text {L}}(\eta ){\varvec{V}}_{4\text {L}}^{*}(\eta _{\mathrm {ref}})\right\rangle +\chi _{4}^2\left\langle {\varvec{V}}_2^2 (\eta ) {\varvec{V}}_2^{*2}(\eta _{\mathrm {ref}})\right\rangle }, \end{aligned}$$
(12)
$$\begin{aligned} r_{5|5;1}(\eta )\approx & {} \frac{\left\langle {\varvec{V}}_{5\text {L}}(-\eta ){\varvec{V}}_{5\text {L}}^{*}(\eta _{\mathrm {ref}})\right\rangle +\chi _{5}^2\left\langle {\varvec{V}}_2 (-\eta ){\varvec{V}}_2^{*}(\eta _{\mathrm {ref}}){\varvec{V}}_3 (-\eta ){\varvec{V}}_3^{*}(\eta _{\mathrm {ref}})\right\rangle }{\left\langle {\varvec{V}}_{5\text {L}}(\eta ){\varvec{V}}_{5\text {L}}^{*}(\eta _{\mathrm {ref}})\right\rangle +\chi _{5}^2\left\langle {\varvec{V}}_2 (\eta ){\varvec{V}}_2^{*}(\eta _{\mathrm {ref}}){\varvec{V}}_3 (\eta ){\varvec{V}}_3^{*}(\eta _{\mathrm {ref}})\right\rangle }. \end{aligned}$$
(13)

Therefore, both the linear and nonlinear components are important for \(r_{4|4;1}\) and \(r_{5|5;1}\).

3 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.Footnote 2 The measurements were performed using the inner detector (ID), minimum-bias trigger scintillators (MBTS), the forward calorimeters (FCal), and the zero-degree calorimeters (ZDC). The ID detects charged particles within \(|\eta | < 2.5\) using a combination of silicon pixel detectors, silicon microstrip detectors (SCT), and a straw-tube 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 \(\mathrm{TeV}\) measurements. The MBTS system detects charged particles over \(2.1\lesssim |\eta |\lesssim 3.9\) using two hodoscopes of counters positioned at \(z = \pm \)3.6 m. The FCal consists of three sampling layers, longitudinal in shower depth, and covers \(3.2<|\eta |< 4.9\). The ZDC are positioned at ±140 m from the IP, detecting neutrons and photons with \(|\eta |>8.3\).

This analysis uses approximately 7 and 470 \(\upmu \text {b}^{-1}\) of Pb+Pb data at \(\sqrt{s_{\mathrm {NN}}}= 2.76\) and 5.02 \(\mathrm{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 \(\upmu \text {b}^{-1}\) and 22 \(\upmu \text {b}^{-1}\) of minimum-bias Pb+Pb data at \(\sqrt{s_{\mathrm {NN}}}= 2.76\) \(\mathrm{TeV}\) and \(\sqrt{s_{\mathrm {NN}}}= 5.02\) \(\mathrm{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 (\(\Sigma E_{\text {T}}\,\)) in the FCal to be more than 4.54 TeV. This ultra-central trigger sampled 470 \(\upmu \hbox {b}^{-1}\) of Pb+Pb collisions at 5.02 TeV and was fully efficient for collisions with centrality 0–0.1% (see Sect. 4).

4 Event and track selection

The offline event selection requires a reconstructed vertex with its z position satisfying \(|Z_{\text {vtx}}\,|< 100\) mm. For the \(\sqrt{s_{\mathrm {NN}}}= 2.76\) \(\mathrm{TeV}\) Pb+Pb data, the selection also requires a time difference \(|\Delta 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 \(\Sigma E_{\text {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_{\mathrm {ch}}^{\mathrm {rec}}\). Since the distribution of \(\Sigma E_{\text {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 \(\Sigma E_{\text {T}}\,\) or \(N_n\) as a function of \(N_{\mathrm {ch}}^{\mathrm {rec}}\).

The event centrality [34] is characterised by the \(\Sigma E_{\text {T}}\,\) deposited in the FCal over the pseudorapidity range \(3.2< |\eta | < 4.9\) using a calibration employing the electromagnetic calorimeters to set the energy scale [35]. The FCal \(\Sigma E_{\text {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_{\mathrm {part}}\), for each centrality interval. The systematic uncertainty in \(N_{\mathrm {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 \(\Sigma E_{\text {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_{\mathrm {part}}\).

Charged-particle tracks and primary vertices [37] are reconstructed from hits in the ID. Tracks are required to have \(p_{\text {T}}\,>0.5\) GeV and \(|\eta |<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, \(\epsilon (p_{\text {T}}\,,\eta )\), 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 \(\textsc {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 \(\eta \approx 0\) to about 50% for \(|\eta | > 2\) for charged particles with \(p_{\text {T}}\,> 0.8\) GeV, falling by about 5% as \(p_{\text {T}}\,\) is reduced to 0.5 GeV. The efficiency varies more strongly with \(\eta \) and event multiplicity. For \(p_{\text {T}}\,> 0.8\) GeV, it ranges from 75% at \(\eta \approx 0\) to 50% for \(|\eta | > 2\) in peripheral collisions, while it ranges from 71% at \(\eta \approx 0\) to about 40% for \(|\eta | > 2\) in central collisions. The tracking efficiency for the 2.76 TeV data has a similar dependence on \(p_{\text {T}}\,\) and \(\eta \). 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_{\text {T}}\,<1\) GeV in central collisions, where its percentage per-track ranges from 2% at \(|\eta |<1\) to 8% at the larger \(|\eta |\). The fake rate drops rapidly for higher \(p_{\text {T}}\,\) and towards more peripheral collisions. The fake rate is accounted for in the tracking efficiency correction following the procedure in Ref. [44].

5 Data analysis

Measurement of the longitudinal flow dynamics requires the calculation of the flow vector \({\varvec{q}}_n\) via Eq. (1) in the ID and the FCal. The flow vector from the FCal serves as the reference \({\varvec{q}}_n(\eta _{\mathrm {ref}})\), while the ID provides the flow vector as a function of pseudorapidity \({\varvec{q}}_n(\eta )\).

In order to account for detector inefficiencies and non-uniformity, a particle weight for the \(i{\text {th}}\)-particle in the ID for the flow vector from Eq. (1) is defined as:

$$\begin{aligned} w_i^{\text {ID}}(\eta ,\phi ,p_{\text {T}}\,) = d_{\text {ID}} (\eta ,\phi )/\epsilon (\eta ,p_{\text {T}}\,), \end{aligned}$$
(14)

similar to the procedure in Ref. [44]. The determination of track efficiency \(\epsilon (\eta ,p_{\text {T}}\,)\) is described in Sect. 4. The additional weight factor \(d_{\text {ID}}(\eta ,\phi )\) corrects for variation of tracking efficiency or non-uniformity of detector acceptance as a function of \(\eta \) and \(\phi \). For a given \(\eta \) interval of 0.1, the distribution in azimuthal bins, \(N(\phi ,\eta )\), is built up from reconstructed charged particles summed over all events. The weight factor is then obtained as \(d_{\text {ID}}(\eta ,\phi ) \equiv \left\langle N(\eta )\right\rangle /N(\phi ,\eta )\), where \(\left\langle N(\eta )\right\rangle \) is the average of \(N(\phi ,\eta )\). This “flattening” procedure removes most \(\phi \)-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:

$$\begin{aligned} w_i^{\text {FCal}}(\eta ,\phi ) = d_{\text {FCal}}(\eta ,\phi ) E_{{\text {T}},i}, \end{aligned}$$
(15)

where \(E_{{\text {T}},i}\) is the transverse energy measured in the \(\text {i}{\text {th}}\) tower in the FCal at \(\eta \) and \(\phi \). The azimuthal weight \(d_{\text {FCal}}(\eta ,\phi )\) is calculated in narrow \(\eta \) intervals in a similar way to what is done for the ID. It ensures that the \(E_{\text {T}}\)-weighted distribution, averaged over all events in a given centrality interval, is uniform in \(\phi \). The flow vectors \({\varvec{q}}_n(\eta )\) and \({\varvec{q}}_n(\eta _{\mathrm {ref}})\) are further corrected by an event-averaged offset: \({\varvec{q}}_n -\left\langle {\varvec{q}}_n\right\rangle _{\text {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:

$$\begin{aligned} \left\langle {\varvec{q}}_n^k (\eta ) {\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle\equiv & {} \left\langle {\varvec{q}}_n^k (\eta ) {\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle _{\text {s}}\nonumber \\&-\left\langle {\varvec{q}}_n^k (\eta ) {\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle _{\text {b}}, \end{aligned}$$
(16)

where subscripts “s” and “b” represent the correlator constructed from the same event (“signal”) and from the mixed-event (“background”), respectively. The mixed-event quantity is constructed by combining \({\varvec{q}}_n^k (\eta )\) from each event with \({\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\) obtained in other events with similar centrality (within 1%) and similar \(Z_{\text {vtx}}\,\) (\(|\Delta Z_{\text {vtx}}\,|<5\) mm). The \(\left\langle {\varvec{q}}_n^k (\eta ) {\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle _{\text {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 \(\eta \) ranges.

For correlators involving flow vectors in two different \(\eta \) ranges, mixed events are constructed from two different events. For example, the correlation for \(r_{2,3|5}\) is calculated as:

$$\begin{aligned} \left\langle {\varvec{q}}_2 (\eta ) {\varvec{q}}_3 (\eta ) {\varvec{q}}_5^{*}(\eta _{\mathrm {ref}})\right\rangle\equiv & {} \left\langle {\varvec{q}}_2 (\eta ) {\varvec{q}}_3 (\eta ) {\varvec{q}}_5^{*}(\eta _{\mathrm {ref}})\right\rangle _{\text {s}}\nonumber \\&-\left\langle {\varvec{q}}_2 (\eta ) {\varvec{q}}_3 (\eta ) {\varvec{q}}_5^{*}(\eta _{\mathrm {ref}})\right\rangle _{\text {b}}. \end{aligned}$$
(17)

The mixed-event correlator is constructed by combining \({\varvec{q}}_2 (\eta ) {\varvec{q}}_3 (\eta )\) from one event with \({\varvec{q}}_5^{*}(\eta _{\mathrm {ref}})\) obtained in another event with similar centrality (within 1%) and similar \(Z_{\text {vtx}}\,\) (\(|\Delta Z_{\text {vtx}}\,|<5\) mm). On the other hand, for correlators involving more than two different \(\eta \) ranges, mixed events are constructed from more than two different events, one for each unique \(\eta \) 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_{\text {vtx}}\,\).

Most correlators can be symmetrised. For example, in a symmetric system such as Pb+Pb collisions, the condition \(\left\langle {\varvec{q}}_n^k (-\eta ) {\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})\right\rangle = \left\langle {\varvec{q}}_n^k (\eta ) {\varvec{q}}_n^{*k}(-\eta _{\mathrm {ref}})\right\rangle \) holds. So instead of Eq. (2), the actual measured observable is:

$$\begin{aligned} r_{n|n;k}(\eta ) = \frac{\left\langle {\varvec{q}}_n^k (-\eta ) {\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})+{\varvec{q}}_n^k (\eta ) {\varvec{q}}_n^{*k}(-\eta _{\mathrm {ref}})\right\rangle }{\left\langle {\varvec{q}}_n^k (\eta ){\varvec{q}}_n^{*k}(\eta _{\mathrm {ref}})+{\varvec{q}}_n^k (-\eta ){\varvec{q}}_n^{*k}(-\eta _{\mathrm {ref}})\right\rangle }. \end{aligned}$$
(18)

The symmetrisation procedure also allows further cancellation of possible differences between \(\eta \) and \(-\eta \) 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 \(\eta \) for \(|\eta |<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 compare with results obtained from the 5.02 TeV data and study the dependence on collision energy.

Table 1 The list of observables measured in this analysis

Figures 2, 3 show the sensitivity of \(r_{2|2;1}\) and \(r_{3|3;1}\), respectively, to the choice of the range of \(\eta _{\mathrm {ref}}\). A smaller \(\eta _{\mathrm {ref}}\) value implies a smaller pseudorapidity gap between \(\eta \) and \(\eta _{\mathrm {ref}}\). The values of \(r_{n|n;1}\) generally decrease with decreasing \(\eta _{\mathrm {ref}}\), possibly reflecting the contributions from the dijet correlations [5]. However, such contributions should be reduced in the most central collisions due to large charged-particle multiplicity and jet-quenching [45] effects. Therefore, the decrease of \(r_{n|n;1}\) in the most central collisions may also reflect the \(\eta _{\mathrm {ref}}\) dependence of \(F_{n;1}^{\text {r}}\), as defined in Eq. (3). In this analysis, the reference flow vector is calculated from \(4.0<\eta _{\mathrm {ref}}<4.9\), which reduces the effect of dijets and provides good statistical precision. For this choice of \(\eta _{\mathrm {ref}}\) range, \(r_{2|2;1}\) and \(r_{3|3;1}\) show a linear decrease as a function of \(\eta \) in most centrality intervals, indicating a significant breakdown of factorisation. A similar comparison for \(r_{4|4;1}\) can be found in the “Appendix”.

Fig. 2
figure 2

The \(r_{2|2;1}(\eta )\) measured for several \(\eta _{\mathrm {ref}}\) ranges. Each panel shows the results for one centrality range. The error bars are statistical only

Fig. 3
figure 3

The \(r_{3|3;1}(\eta )\) measured for several \(\eta _{\mathrm {ref}}\) ranges. Each panel shows the results for one centrality range. The error bars are statistical only

Figures 4, 5 show \(r_{2|2;1}\) and \(r_{3|3;1}\) calculated for several \(p_{\text {T}}\,\) ranges of the charged particles in the ID. A similar comparison for \(r_{4|4;1}\) can be found in the “Appendix”. If the longitudinal-flow asymmetry and twist reflect global properties of the event, the values of \(r_{n|n;1}\) should not depend strongly on \(p_{\text {T}}\,\). Indeed no dependence is observed, except for \(r_{2|2;1}\) in the most central collisions and very peripheral collisions. The behaviour in central collisions may be related to the factorisation breaking of the \(v_2\) as a function of \(p_{\text {T}}\,\) and \(\eta \) [5, 19]. The behaviour in peripheral collisions is presumably due to increasing relative contributions from jets and dijets at higher \(p_{\text {T}}\,\) and for peripheral collisions. Based on this, the measurements are performed using charged particles with \(0.5<p_{\text {T}}\,<3\) GeV.

Fig. 4
figure 4

The \(r_{2|2;1}(\eta )\) measured in several \(p_{\text {T}}\,\) ranges. Each panel shows the results for one centrality range. The error bars are statistical only

Fig. 5
figure 5

The \(r_{3|3;1}(\eta )\) measured in several \(p_{\text {T}}\,\) ranges. Each panel shows the results for one centrality range. The error bars are statistical only

6 Systematic uncertainties

Since all observables are found to follow an approximately linear decrease with \(\eta \), i.e. \(D(\eta ) \approx 1-c\eta \) for a given observable \(D(\eta )\) where c is a constant, the systematic uncertainty is presented as the relative uncertainty for \(1-D(\eta )\) at \(\eta =1.2\), the mid-point of the \(\eta \) 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(\eta ,\phi )\) 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 statistical 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_{\mathrm {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 \theta |\) 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_{\text {T}}\,\) but relatively weak dependence on \(\eta \). Therefore, a \(p_{\text {T}}\,\)-dependent uncertainty in the track reconstruction efficiency \(\epsilon (\eta ,p_{\text {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 \(\eta \) and \(p_{\text {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 \({\varvec{q}}_n(\eta _{\mathrm {ref}})\) calculated from the azimuthal distribution of the \(E_{\text {T}}\,\) via Eqs. (1) and (15) differs from the true azimuthal distribution. However, since \({\varvec{q}}_n(\eta _{\mathrm {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_{\text {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_{\text {T}}\,\) above the \(50^{\mathrm {th}}\) percentile are used to calculate the \({\varvec{q}}_n(\eta _{\mathrm {ref}})\). The \(|{\varvec{q}}_n(\eta _{\mathrm {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 \({\varvec{q}}_n(\eta _{\mathrm {ref}})\) is calculated using both the generated \(E_{\text {T}}\,\) and the reconstructed \(E_{\text {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.

Table 2 Systematic uncertainties in percent for \(1-r_{2|2;k}\) and \(1-r_{3|3;k}\) at \(\eta =1.2\) in selected centrality intervals
Table 3 Systematic uncertainties in percent for \(1-R_{2|2;2}\), \(1-R_{3|3;2}\), \(1-r_{4|4;1}\) and \(1-r_{5|5;1}\) at \(\eta =1.2\) in selected centrality intervals
Table 4 Systematic uncertainties in percent for \(1-r_{2,3|2,3}\), \(1-r_{2,2|4}\) and \(1-r_{2,3|5}\) at \(\eta =1.2\) in selected centrality intervals

7 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 mixed-harmonic 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_{\text {T}}\,<3\) GeV, and the reference flow vector is calculated with \(4.0<|\eta _{\mathrm {ref}}|<4.9\). Most results are shown for the \(\sqrt{s_{\mathrm {NN}}}=5.02\) TeV Pb+Pb dataset, which has better statistical precision. The results for the \(\sqrt{s_{\mathrm {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

Fig. 6
figure 6

The \(r_{2|2;1}(\eta )\) compared between the two collision energies. Each panel shows results from one centrality interval. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

Figure 6 shows \(r_{2|2;1}\) in various centrality intervals at the two collision energies. The correlator shows a linear decrease with \(\eta \), 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 \(\sqrt{s_{\mathrm {NN}}}=2.76\) TeV is slightly stronger than that at \(\sqrt{s_{\mathrm {NN}}}=5.02\) TeV, which is expected as the collision system becomes less boost-invariant 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 \(\eta \) 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.

Fig. 7
figure 7

The \(r_{3|3;1}(\eta )\) compared between the two collision energies. Each panel shows results from one centrality interval. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

Fig. 8
figure 8

The \(r_{4|4;1}(\eta )\) compared between the two collision energies. Each panel shows results from one centrality interval. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

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 \(\eta \) 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 \(\eta _{\mathrm {ref}}\) chosen to be \(4.4<|\eta _{\mathrm {ref}}|<4.9\) to match the CMS choice of \(\eta _{\mathrm {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 \(\eta \) and the decreasing trends are stronger at lower collision energy.

Fig. 9
figure 9

The \(R_{2|2;2}(\eta )\) compared between the two collision energies. Each panel shows results from one centrality interval. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

Fig. 10
figure 10

The \(R_{3|3;2}(\eta )\) compared between the two collision energies. Each panel shows results from one centrality interval. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

The measured \(r_{n|n;k}\) and \(R_{n|n;2}\) are parameterised with linear functions,

$$\begin{aligned} r_{n|n;k} = 1-2F_{n;k}^{\text {r}}\;\eta ,\;\;R_{n|n;2} = 1-2F_{n;2}^{\text {R}}\;\eta , \end{aligned}$$
(19)

where the slope parameters are calculated as linear-regression coefficients,

$$\begin{aligned} F_{n;k}^{\text {r}}= & {} \frac{\sum _i(1-r_{n|n;k}(\eta _i))\eta _i}{2\sum _i\eta _i^2},\nonumber \\ F_{n;2}^{\text {R}}= & {} \frac{\sum _i(1-R_{n|n;2}(\eta _i))\eta _i}{2\sum _i\eta _i^2}, \end{aligned}$$
(20)

which characterise the average \(\eta \)-weighted deviation of \(r_{n|n;1}(\eta )\) and \(R_{n|n;2}(\eta )\) from unity. The sum runs over all data points. If \(r_{n|n;k}\) and \(R_{n|n;2}\) are a linear function in \(\eta \), 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_{n;1}^{\text {r}}\) and \(F_{n;2}^{\text {R}}\) are plotted as a function of centrality in terms of \(N_{\mathrm {part}}\), in Figs. 11 and 12, respectively. The values of \(F_{2;1}^{\text {r}}\) and \(F_{2;2}^{\text {R}}\) first decrease and then increase as a function of increasing \(N_{\mathrm {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_{\mathrm {part}}\), a slight decrease in \(F_{3;1}^{\text {r}}\) and \(F_{3;2}^{\text {R}}\) is observed for \(N_{\mathrm {part}}>200\), as well as an increase in \(F_{4;1}^{\text {r}}\) for \(N_{\mathrm {part}}<100\). The values of \(F_{n;1}^{\text {r}}\) and \(F_{n;2}^{\text {R}}\) are larger with decreasing \(\sqrt{s_{\mathrm {NN}}}\), as the rapidity profile of the initial state is more compressed due to smaller beam rapidity \(y_{\mathrm {beam}}\) at lower \(\sqrt{s_{\mathrm {NN}}}\). This energy dependence has been predicted for \(F_{n;1}^{\text {r}}\) in hydrodynamic model calculations [24], and it is quantified in Fig. 13 via the ratio of \(F_{2;1}^{\text {r}}\) values and of \(F_{2;2}^{\text {R}}\) values at the two energies. The weighted averages of the ratios calculated in the range \(30<N_{\mathrm {part}}<400\) are given in Table 5. Compared to \(\sqrt{s_{\mathrm {NN}}}=5.02\) TeV, the values of \(F_{2;1}^{\text {r}}\) and \(F_{2;2}^{\text {R}}\) at \(\sqrt{s_{\mathrm {NN}}}=2.76\) TeV are about 10% higher, and the values of \(F_{3;1}^{\text {r}}\) and \(F_{4;1}^{\text {r}}\) are about 16% higher.

If the change of correlators with \(\sqrt{s_{\mathrm {NN}}}\) were entirely due to the change of \(y_{\mathrm {beam}}\), then the correlators would be expected to follow a universal curve when they are rescaled by \(y_{\mathrm {beam}}\), i. e. \(r_{n|n;k}(\eta /y_{\mathrm {beam}})\) and \(R_{n|n;2}(\eta /y_{\mathrm {beam}})\) should not depend on \(\sqrt{s_{\mathrm {NN}}}\). In this case, the slopes parameters multiplified by the beam rapidity, \(\hat{F}_{n;1}^{\text {r}}\equiv F_{n;1}^{\text {r}}y_{\mathrm {beam}}\) and \(\hat{F}_{n;2}^{\text {R}}\equiv F_{n;2}^{\text {R}}y_{\mathrm {beam}}\), should not depend on \(\sqrt{s_{\mathrm {NN}}}\). The beam rapidity is \(y_{\mathrm {beam}}=7.92\) and 8.52 for \(\sqrt{s_{\mathrm {NN}}}=2.76\) and 5.02 TeV, respectively, which leads to a 7.5% reduction in the ratio. Figure 14 shows the ratio of \(\hat{F}_{2;1}^{\text {r}}\) values and of \(\hat{F}_{2;2}^{\text {R}}\) values at the two energies, and the weighted averages of the ratios calculated in the range \(30<N_{\mathrm {part}}<400\) are given in Table 5. The \(y_{\mathrm {beam}}\)-scaling accounts for a large part of the \(\sqrt{s_{\mathrm {NN}}}\) dependence. Compared to \(\sqrt{s_{\mathrm {NN}}}=5.02\) TeV, the values of \(\hat{F}_{2;1}^{\text {r}}\) and \(\hat{F}_{2;2}^{\text {R}}\) at \(\sqrt{s_{\mathrm {NN}}}=2.76\) TeV are about 3% higher, and the values of \(\hat{F}_{3;1}^{\text {r}}\) and \(\hat{F}_{4;1}^{\text {r}}\) are about 8% higher, so this level of difference remains after accounting for the change in the beam rapidity.

Fig. 11
figure 11

Centrality dependence of \(F_{2;1}^{\text {r}}\) (left panel), \(F_{3;1}^{\text {r}}\) (middle panel) and \(F_{4;1}^{\text {r}}\) (right panel) for Pb+Pb at 2.76 TeV (circles) and 5.02 TeV (squares). 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

Fig. 12
figure 12

Centrality dependence of \(F_{2;2}^{\text {R}}\) (left panel), \(F_{3;2}^{\text {R}}\) (middle panel) and \(F_{4;2}^{\text {R}}\) (right panel) for Pb+Pb at 2.76 TeV (circles) and 5.02 TeV (squares). 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

Fig. 13
figure 13

Centrality dependence of ratio of \(F_{n;1}^{\text {r}}\) values (left panel) and \(F_{n;2}^{\text {R}}\) values (right panel) at 2.76 and 5.02 TeV. The lines indicate the average values in the range \(30<N_{\mathrm {part}}<400\), with the results and fit uncertainties given by Table 5. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

Fig. 14
figure 14

Centrality dependence of ratio of \(\hat{F}_{n;1}^{\text {r}} \equiv F_{n;1}^{\text {r}}y_{\mathrm {beam}}\) values (left panel) and \(\hat{F}_{n;2}^{\text {R}}\equiv F_{n;2}^{\text {R}}y_{\mathrm {beam}}\) values (right panel) at 2.76 and 5.02 TeV. The lines indicate the average values in the range \(30<N_{\mathrm {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_{n;1}^{\text {r}}\), \(F_{n;2}^{\text {R}}\), \(\hat{F}_{n;1}^{\text {r}}\equiv F_{n;1}^{\text {r}}y_{\mathrm {beam}}\) and \(\hat{F}_{n;2}^{\text {R}}\equiv F_{n;2}^{\text {R}}y_{\mathrm {beam}}\) values at the two energies in the range \(30<N_{\mathrm {part}}<400\) shown in Figs. 13 and 14. The uncertainties include both statistical and systematic uncertainties

7.2 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).

Figure 15 compares the results for \(r_{2|2;k}\) for \(k=1\)–3 (solid symbols) with \(r_{2|2;1}^k\) for \(k=2\)–3 (open symbols). The data follow the scaling relation from Eq. (4) in the most central collisions (0–5% centrality) 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 \(\eta \) than \(r_{2|2;1}^k\).

Fig. 15
figure 15

The \(r_{2|2;k}\) for \(k=1\)–3 compared with \(r_{2|2;1}^k\) 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

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.

Fig. 16
figure 16

The \(r_{3|3;k}\) for \(k=1\)–3 compared with \(r_{3|3;1}^k\) 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

To quantify the difference between \(r_{n|n;k}\) and \(r_{n|n;1}^k\), the slopes (\(F_{n;k}^{\text {r}}\)) of \(r_{n|n;k}\) are calculated via Eqs. (19) and (20). The scaled quantities, \(F_{n;k}^{\text {r}}/k\), are then compared with each other as a function of centrality in Fig. 17. For second-order harmonics, the data show clearly that over most of the centrality range \(F_{2;3}^{\text {r}}/3>F_{2;2}^{\text {r}}/2>F_{2;1}^{\text {r}}\), implying \(F_{2;k}^{\text {r}}>kF_{2;1}^{\text {r}}\). However, for the most central and most peripheral collisions the quantities approach each other. On the other hand, a slightly opposite trend for the third-order harmonics, \(F_{3;3}^{\text {r}}/3\lesssim F_{3;2}^{\text {r}}/2\lesssim F_{3;1}^{\text {r}}\), i.e. \(F_{3;k}^{\text {r}}\lesssim kF_{3;1}^{\text {r}}\), is observed in mid-central collisions (\(150<N_{\mathrm {part}}<350\)).

Fig. 17
figure 17

The values of \(F_{n;k}^{\text {r}}/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

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].

Fig. 18
figure 18

The \(r_{2|2;2}(\eta )\) and \(R_{2|2;2}(\eta )\) in various centrality intervals for Pb+Pb collisions at 5.02 TeV. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

Fig. 19
figure 19

The \(r_{3|3;2}(\eta )\) and \(R_{3|3;2}(\eta )\) 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 in 40–50% centrality interval are rebinned to reduce the uncertainty

Following the discussion in Sect. 2, Eqs. (3) and (6), the measured \(F_{n;2}^{\text {r}}\) and \(F_{n;2}^{\text {R}}\) values can be used to estimate the separate contributions from FB asymmetry and event-plane twist, \(F_{n;2}^{\text {asy}}\) and \(F_{n;2}^{\text {twi}}\), respectively, via the relation:

$$\begin{aligned} F_{n;2}^{\text {twi}} = F_{n;2}^{\text {R}},\;\;\; F_{n;2}^{\text {asy}} = F_{n;2}^{\text {r}}- F_{n;2}^{\text {R}}. \end{aligned}$$
(21)

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.

Fig. 20
figure 20

The estimated event-plane twist component \(F_{n;2}^{\text {twi}}\) and FB asymmetry component \(F_{n;2}^{\text {asy}}\) as a function of \(N_{\mathrm {part}}\) for \(n=2\) and 3 for Pb+Pb collisions at 5.02 TeV. The error bars and shaded boxes are statistical and systematic uncertainties, respectively

7.3 Mixed-harmonics correlation

Fig. 21
figure 21

The \(r_{2,3|2,3}\) (circles) and \(r_{2|2;1}r_{3|3;1}\) (squares) as a function of \(\eta \) for several centrality intervals. The error bars and shaded boxes are statistical and systematic uncertainties, respectively. The \(r_{2,3|2,3}\) data in the 50–60% centrality interval are rebinned to reduce the uncertainty

Fig. 22
figure 22

Comparison of \(r_{2|2;2}\), \(r_{2,2|4}\) and \(r_{4|4;1}\) for several centrality intervals. The error bars and shaded boxes are statistical and systematic uncertainties, respectively. The data points in some centrality intervals are rebinned to reduce the uncertainty

Figure 21 compares the \(r_{2,3|2,3}\) with the product of \(r_{2|2;1}\) and \(r_{3|3;1}\). The data show that they are consistent with each other, suggesting the previously observed anticorrelation beween \(v_2\) and \(v_3\) is a property of the entire event [9, 48], and that longitudinal fluctuations of \(v_2\) and \(v_3\) are uncorrelated. Figure 22 compares \(r_{2|2;2}\) with the mixed-harmonic correlator \(r_{2,2|4}\), as well as \(r_{4|4;1}\). As discussed in Sect. 2 in the context of the first relation in Eq. (10), if the linear and non-linear components of \(v_4\) in Eq. (10) are uncorrelated, then \(r_{2,2|4}\) would be expected to be similar to \(r_{2|2;24}\). This is indeed confirmed by the comparisons of the \(\eta \) and centrality dependence of \(r_{2|2;2}\) and \(r_{2,2|4}\) in Fig. 22. Figure 22 also shows that the \(\eta \) 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)).

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.

Fig. 23
figure 23

Comparison of \(r_{2,3|2,3}\), \(r_{2,3|5}\) and \(r_{5|5;1}\) for several centrality intervals. The error bars and shaded boxes are statistical and systematic uncertainties, respectively. The \(r_{5|5;1}\) data in some centrality intervals are rebinned to reduce the uncertainty

The decorrelations shown in Figs. 21, 22 and 23 can be quantified by calculating the slopes of the distributions in each centrality interval and presenting the results as a function of centrality. Following the example for \(r_{n|n;k}\), the slopes for the mixed-harmonic correlators are obtained via the linear regression procedure of Eqs. (19) and (20):

$$\begin{aligned} r_{2,3|2,3} = 1-2F_{2,3|2,3}^{\text {r}}\;\eta ,\;\;r_{2,2|4} = 1-2F_{2,2|4}^{\text {r}}\;\eta ,\;\;r_{2,3|5} = 1-2F_{2,3|5}^{\text {r}}\;\eta . \end{aligned}$$
(22)

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}\).

Fig. 24
figure 24

Comparison of the slopes of the correlators as a function of \(N_{\mathrm {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

8 Summary

Measurements of longitudinal flow correlations for charged particles are presented in the pseudorapidity range \(|\eta |<2.4\) using 7 and 470 \(\upmu \hbox {b}^{-1}\) of Pb+Pb data at \(\sqrt{s_{\mathrm {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 \(\eta \) separation between the two particles. The slope of this dependence is nearly independent of centrality and \(p_{\text {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_{\text {T}}\,\). Furthermore, the effect is found to be larger 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 \(\eta \)-dependent flow correlations are also measured and the corresponding linear coefficients of the \(\eta \) dependence are extracted. The coefficient for the \(k{\text {th}}\)-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_2v_3\) between two \(\eta \) 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\propto v_2^2\). Similarly, the correlations between \(v_5\) and \(v_2v_3\) suggest that the longitudinal fluctuations of \(v_5\) are driven by the nonlinear contribution from \(v_2v_3\), i.e. \(v_5\propto v_2v_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.