1 Introduction

One of the primary goals of ultra-relativistic heavy-ion collisions is the study of the hot and dense medium formed there, usually referred to as the quark-gluon plasma (QGP) [1,2,3,4,5]. The existence of the QGP phase of nuclear matter has been confirmed by a wealth of experimental data [5, 6]. In particular, properties related to the collective expansion of the QGP (e.g. the equation of state [7] and shear viscosity [8]) are inferred from measurements of azimuthal anisotropies of produced particles. It is now understood that the azimuthal anisotropy results from large initial pressure gradients in the hot, dense matter created in the collisions [9, 10]. These pressure gradients transform the initial spatial anisotropies of nuclear collisions into momentum anisotropies of the final-state particle production, which are experimentally characterized by Fourier (flow) harmonics of the azimuthal angle distributions of produced particles. The discovery of large flow harmonics at RHIC, and more recently at much higher collision energy at the LHC [11,12,13,14], has significantly deepened the understanding of the QGP, as explored theoretically by the QCD lattice [15]. In particular, the recent measurements of azimuthal anisotropy help to constrain the commonly used modelling of the dynamics of heavy-ion collisions based on relativistic viscous hydrodynamics. Typically, in the hydrodynamic models, a strongly interacting quark–gluon medium is formed shortly after the collision and its evolution is well described by relativistic fluid dynamics [8]. Detailed investigations, based on hydrodynamics, have shown that the produced medium has properties similar to those of an almost ideal fluid characterized by a very low ratio of viscosity to entropy density, \(\eta /s\). Precise azimuthal anisotropy measurements over a wide range in kinematic variables and centrality are key elements to improving our understanding of the strongly coupled QGP because of their unique sensitivity to \(\eta /s\).

The azimuthal angular distribution of single produced particles can be expanded in a Fourier series [16, 17]:

$$\begin{aligned} \frac{\text {d}N}{\text {d}\phi } = \frac{N_{0}}{2\pi }\left( 1+ \sum _{n=1} 2v_{n} \cos \left[ n\left( \phi - \Phi _n\right) \right] \right) , \end{aligned}$$
(1)

where \(N_0\) is the total particle yield, \(\phi \) is the azimuthal angle of the produced particles and the \(v_{n} \) and \(\Phi _n\) are, respectively, the magnitude of the nth-order azimuthal anisotropy and the orientation of the nth-order symmetry plane. The \(v_{n}\) coefficients – also called flow harmonics – are typically measured as a function of particle pseudorapidityFootnote 1 (\(\eta \)), transverse momentum (\(p_{\mathrm{T}}\)), and the degree of overlap between the colliding nuclei (centrality). Event-by-event fluctuations in the number and position of the interacting nucleons give rise to anisotropic flow fluctuations [18].

The first harmonic, \(v_{1}\), is known as directed flow and refers to the sideward motion of participants in ultra-relativistic nuclear collisions, and it carries information from the early stage of the collision. The most extensive studies are related to the second flow harmonic, \(v_{2}\), also known as elliptic flow. Elliptic flow is sensitive to the initial spatial asymmetry of the almond-shaped overlapping zone of the colliding nuclei. The higher-order coefficients, \(n>2\), are also important due to their sensitivity to initial-state geometric fluctuations and viscous effects [16,17,18].

During the first operational period at the LHC (Run 1) lead ions were collided at energy per colliding nucleon–nucleon pair \(\sqrt{s_{_\text {NN}}}\) = 2.76 \(\mathrm{TeV}\), which is about 13 times larger than the highest collision energy attained at RHIC in Au+Au collisions. ATLAS and other LHC experiments collected large samples of heavy-ion data enabling extensive studies of the elliptic flow and higher-order Fourier coefficients. ATLAS measurements of flow harmonics were performed in broad regions of transverse momentum, pseudorapidity and event centrality, using the standard event-plane (EP) method [12], two-particle correlations (2PC) [13] and multi-particle cumulants [19]. Significant (non-zero) flow harmonics up to \(v_{6}\) were measured in Pb+Pb collisions at \(\sqrt{s_{_\text {NN}}}\) = 2.76 \(\mathrm{TeV}\), which provide important constraints on the bulk and shear viscosity of the QGP medium [20]. Additionally, by comparing RHIC (STAR [21] and PHENIX [22]) and LHC (ATLAS [12], ALICE [23] and CMS [24]) results, it was found that for similar centrality of Au+Au and Pb+Pb interactions, \(v_{n}\) as a function of \(p_{\mathrm{T}}\)is approximately independent of collision energy. There is an initial rise of \(v_{n}\) with \(p_{\mathrm{T}}\)up to about 3 \(\mathrm{GeV}\) and then a drop-off at higher values of \(p_{\mathrm{T}}\), and only weak dependence for \(p_{\mathrm{T}}~> ~8{-}9~\mathrm{GeV}\). As a function of centrality, there is similarly little variation with collision energy. The second harmonic, \(v_{2}\), exhibits the most pronounced centrality variation, rising to a maximum for mid-central collisions, and then falling off for the most central collisions, reflecting variations in the shape of the initial collision geometry. The harmonic, \(v_{3}\), referred as triangular flow, which has a value similar to \(v_{2}\) in central collisions, shows a weaker dependence on centrality, as do the higher-order harmonics.

At the start of the second operational period of the LHC (Run 2), in November and December of 2015, lead–lead collisions with higher collision energy per nucleon pair of \(\sqrt{s_{_\text {NN}}}= 5.02\) \(\mathrm{TeV}\) were collected by the LHC experiments. The goal of this paper is to present and discuss the first ATLAS measurements of \(v_{n}\) harmonics at this energy, using the two-particle correlation [17], scalar-product (SP) [25] and event-plane [16, 17] methods. Comparing the 2PC and SP results can quantify the extent to which the two-particle correlations factorize into the product of the flow harmonics corresponding to single-particle angular distributions [26, 27]. While the SP and EP methods are expected to yield similar values of the \(v_{n}\), small variations due to their different sensitivity to initial-state geometric fluctuations can nevertheless occur [28]. To study the energy dependence, the 2PC and EP flow harmonics are compared with previous ATLAS measurements in 2.76 \(\mathrm{TeV}\) Pb+Pb collisions [12, 13]. The results presented in this paper, together with the results on azimuthal anisotropy from other LHC experiments [29, 30], provide further opportunity to study the properties of the QGP, constrain hydrodynamic models, study transport coefficients and extract the temperature dependence of transport coefficients, including \(\eta /s\).

The organization of this paper is as follows. Section 2 gives a brief overview of the ATLAS detector and the subsystems that are used in this analysis. Section 3 describes the datasets, triggers and the offline selection criteria used to select events and charged-particle tracks. Section 4 gives details of the scalar-product, event-plane and two-particle correlation methods, which are used to measure the \(v_{n}\). Section 5 describes the systematic uncertainties associated with the measured \(v_{n}\). Section 6 presents the main results of the analysis, which are the \(p_{\mathrm{T}}\), \(\eta \) and centrality dependence of the \(v_{n}\) and comparisons of results from the different methods. Section 7 gives a summary of the main results and observations.

2 Experimental set-up

The measurements were performed using the ATLAS detector [31] at the LHC. The principal components used in this analysis are the inner detector (ID), minimum-bias trigger scintillators (MBTS), calorimeter, zero-degree calorimeters (ZDC), and the trigger and data acquisition systems. The ID detects charged particles within the pseudorapidity range \(|\eta |~<~2.5\) using a combination of silicon pixel detectors, including the “insertable B-layer” [32, 33] that was installed between Run 1 and Run 2, silicon microstrip detectors (SCT), and a straw-tube transition radiation tracker (TRT), all immersed in a 2 T axial magnetic field [34]. The MBTS system detects charged particles over \(2.07< |\eta | < 3.86\) using two scintillator-based hodoscopes on each side of the detector, positioned at \(z=\pm 3.6~{\mathrm m}\). These hodoscopes were rebuilt between Run 1 and Run 2. The ATLAS calorimeter system consists of a liquid argon (LAr) electromagnetic (EM) calorimeter covering \(|\eta |<3.2\), a steel–scintillator sampling hadronic calorimeter covering \(|\eta |<1.7\), a LAr hadronic calorimeter covering \(1.5<|\eta |<3.2\), and two LAr electromagnetic and hadronic forward calorimeters (FCal) covering \(3.2<|\eta |<4.9\). The ZDC, situated at approximately \(\pm 140\) m from the nominal IP, detect neutral particles, mostly neutrons and photons, with \(|\eta |~>~8.3\). The ZDC use tungsten plates as absorbers, and quartz rods sandwiched between the tungsten plates as the active medium. The ATLAS trigger system [35] consists of a first-level (L1) trigger implemented using a combination of dedicated electronics and programmable logic, and a software-based high-level trigger.

3 Event and track selection

The Pb+Pb dataset used in this paper corresponds to an integrated luminosity of 0.49 \({\mathrm {nb}}^{-1}\). Minimum-bias events were selected by two mutually exclusive triggers :

  • Events with smaller impact parameter (semi-central and central collisions) were selected by a trigger that required the total transverse energy (\(E_{\mathrm{T}}\)) deposited in the calorimeters at L1 to be above 50 \(\mathrm{GeV}\).

  • Collisions with large impact parameter (peripheral events) were selected by a trigger that required the total transverse energy at L1 to be less than 50 \(\mathrm{GeV}\), one neutron on either side in the ZDC (\(|\eta | > 8.3\)), and at least one reconstructed track in the ID.

The minimum-bias triggers sampled a total luminosity of 22 \(\upmu {{\mathrm {b}}}^{-1}\). To enhance the statistics of ultra-central collisions, additional data samples were recorded by two dedicated triggers – UCC-1 and UCC-2 – that selected events in which the total \(E_{\mathrm{T}}\) in the FCal at L1 was more than 4.21 \(\mathrm{TeV}\) and 4.54 \(\mathrm{TeV}\), respectively. The UCC-1 trigger sampled a luminosity of 45 \(\upmu {\mathrm {b}}^{-1}\) while the UCC-2 trigger sampled the entire luminosity of 0.49 \({\mathrm {nb}}^{-1}\). The luminosities sampled by the different triggers are listed in Table 1.

Table 1 The luminosities sampled by the triggers used in the analysis

In the offline analysis the z coordinate of the primary vertex [36] is required to be within 10 cm of the nominal interaction point. The fraction of events containing more than one inelastic interaction (pile-up) is estimated to be at the level of 0.2%. The fraction varies with \(\Sigma E_\mathrm {T}^\mathrm {FCal}\), and for ultra-central collisions it amounts to a few percent. Pile-up events were removed by exploiting the correlations between the transverse energy measured in the FCal and in the ZDC as well as the number of tracks associated with the primary vertex, \(N^{\mathrm {rec}}_{\mathrm {ch}}\). As the pile-up is very small, in a typical pile-up event the track multiplicity associated with the primary vertex belongs to a single Pb+Pb collision, while the energy deposited in calorimeters contains contributions from multiple, mostly two, collisions. Therefore, events with small values of \(N^{\mathrm {rec}}_{\mathrm {ch}}\) and large \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) that differ markedly from those of the majority of Pb+Pb collisions are removed from the analysis [19]. In addition, the anti-correlation between the \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) and the number of neutrons detected in ZDC is also used to suppress pile-up events. Events with the number of neutrons (as recorded in the ZDC) much higher than the number expected from the bulk of events for a given value \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) are rejected.

The heavy-ion collision geometry is defined by its impact parameter, b. As the actual event-by-event impact parameter is not accessible experimentally, the centrality classification is based on the transverse energy measured in the forward calorimeter, \(\Sigma E_\mathrm {T}^\mathrm {FCal}\), which exhibits a strong monotonic correlation with b. A model based on the Monte Carlo (MC) Glauber approach [37, 38] is used to obtain the mapping from the observed \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) to the primary properties, such as the number of binary nucleon–nucleon interactions, \(N_{\mathrm {coll}}\), or the number of nucleons participating in the nuclear collision, \(N_{\mathrm {part}}\), for each centrality interval. The Glauber model also provides a correspondence between the \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) distribution and the sampling fraction of the total inelastic Pb+Pb cross-section, allowing the setting of the centrality percentiles [12]. For this analysis a selection of the 80% most central collisions (i.e. centrality 0–80%) is used to avoid any diffractive, photonuclear, and other inelastic processes that contribute significantly to very peripheral collisions (centrality 80–100%). Additionally, the events selected by UCC-1 and UCC-2 are used only over the 0–1% and 0–0.1% centrality intervals, respectively. Figure 1 shows the distribution of \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) in the data, and thresholds for the selection of several centrality intervals. The correspondence of centrality intervals to \(\langle N_{\mathrm {part}}\rangle \) values is provided in Table 2.

Fig. 1
figure 1

The \(\Sigma E_{\mathrm{T}}^{\mathrm {FCal}}\) distribution in \(\sqrt{s_{_\text {NN}}}\) = 5.02 \(\mathrm{TeV}\) Pb+Pb data for events selected by the minimum-bias trigger. The \(\Sigma E_{\mathrm{T}}^\mathrm {FCal}\) thresholds for several centrality intervals are marked with vertical lines and labelled on the plot. Also shown are the number of events over the 0–1% and 0–0.1% centrality intervals selected by the ultra-central triggers

Table 2 The correspondence between centrality intervals used in the analysis and \(\langle N_{\mathrm {part}} \rangle \) values

In order to study the performance of the ATLAS detector, a minimum-bias sample of 4M Pb+Pb MC events was generated using version 1.38b of HIJING [39]. The effect of flow was added after the generation using an “afterburner” [40] procedure in which the \(p_{\mathrm{T}}\), \(\eta \) and centrality dependence of the \(v_{n}\), as measured in the \(\sqrt{s_{_\text {NN}}}=2.76~\mathrm{TeV}\) Pb+Pb data [13], is implemented by artificially rearranging the \(\phi \) positions of the generated particles. The generated sample was passed through a full simulation of the ATLAS detector using \({\textsc {geant}}\)4 [41], and the simulated events are reconstructed using the same algorithms as used for real data. Charged-particle tracks are reconstructed from the signals in the ID. A reconstruction procedure developed for tracking in dense environments in \(pp\)collisions, and optimized for heavy-ion collisions, was used for this purpose [42]. In the analysis the set of reconstructed tracks is filtered using several selection criteria. The tracks are required to have \(p_{\mathrm{T}}>0.5~\mathrm{GeV}\), \(|\eta |<2.5\), 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. A hit is expected if the extrapolated track crosses an active region of a pixel module that has not been disabled, and a hit is said to be missing when it is expected but not found. In addition, the transverse (\(d_0\)) and longitudinal (\(z_0\sin \theta \)) impact parameters of the track relative to the vertex are required to be less than 1 mm. The track-fit quality parameter \(\chi ^2/\)ndof is required to be less than 6.

The MC sample is used to determine the track-reconstruction efficiency as a function of \(p_{\mathrm{T}}\), \(\eta \) and centrality, \(\epsilon (p_{\mathrm{T}},\eta , {\mathrm {centrality}})\). The efficiency is defined as the fraction of primary [36] charged particles matched to reconstructed tracks. The matching criterion is that the weighted fraction of hits in a reconstructed track originating from a given generated particle is above 30%. Different weights are assigned to pixel, SCT and TRT signals to be more robust against fake tracks, which are defined below. At mid-rapidity (\(|\eta |<1\)) and for events with centrality \(<5\%\), the reconstruction efficiency is \(\sim \) 60% at low \(p_{\mathrm{T}}\) and increases to \(\sim \) 75% at higher \(p_{\mathrm{T}}\). For \(|\eta |>1\) the efficiency decreases to about 40–60% depending on the \(p_{\mathrm{T}}\)and centrality. The reconstruction efficiency depends weakly on the centrality for low-\(p_{\mathrm{T}}\) tracks, for which it is smaller in the most central events by about 5% as compared to mid-central and peripheral collisions. For tracks with \(p_{\mathrm{T}}> 1~\mathrm{GeV}\) the dependence on centrality is less than 1%.

The fraction of tracks that are not matched to primary, generated MC particles or are produced from random combinations of hits in the ID, both referred to as “fake tracks”, is found to depend significantly on \(\eta \). For \(|\eta |<1\), it is \(\sim \)10% for low-\(p_{\mathrm{T}}\) tracks in the most central 5% Pb+Pb events, and about 5% for more peripheral collisions. In the forward part of the detector, especially for \(1<|\eta |<2\) where detector services reside, the fake rate can reach 18% for low \(p_{\mathrm{T}}\)tracks in the most central collisions. The fake rate drops rapidly for higher \(p_{\mathrm{T}}\) and also decreases gradually towards more peripheral collisions. For \(p_{\mathrm{T}}\ > 10~\mathrm{GeV}\) and 0–5% centrality it rises to about 5%.

4 Analysis procedure

Three analysis techniques are used to determine the flow harmonics: the two-particle correlation method, which uses only the information from the tracking detectors, and the scalar-product and event-plane methods, which also use information from the FCal.

In all approaches the differential flow harmonics are first obtained in narrow intervals of \(p_{\mathrm{T}}\), \(\eta \) and centrality. Integrated quantities are obtained by taking into account the track reconstruction efficiency, \(\epsilon \), and fake rate, f. A \(p_{\mathrm{T}}\)-, \(\eta \)- and centrality-dependent weight factor \(w = (1-f)/\epsilon \) is applied to each track in the 2PC measurement and to scale each bin of the differential \(v_{n} \) distributions in the SP and EP methods.

All analysis methods utilize the minimum-bias sample of 22 \(\upmu {\mathrm {b}} ^{-1}\). In addition, the SP and EP analyses use the ultra-central samples of 45 \(\upmu {\mathrm {b}} ^{-1}\) and 0.49 \({\mathrm {nb}}^{-1}\).

4.1 Two-particle correlation analysis

The 2PC method has been used extensively by ATLAS for correlation measurements [13, 43,44,45,46,47,48]. In the 2PC method, the distribution of particle pairs in relative azimuthal angle \(\Delta \phi =\phi ^a-\phi ^b \) and pseudorapidity separation \(\Delta \eta =\eta ^a-\eta ^b \) is measured. Here the labels a and b denote the two particles used to make the pair. They are conventionally called the “trigger” and “associated” particles, respectively. In this analysis, the two particles are charged particles reconstructed by the ATLAS tracking system, over the full azimuth and \(|\eta |<2.5\), resulting in a pair-acceptance coverage of \(\pm 5.0\) units in \(\Delta \eta \).

In order to account for the detector acceptance effects, the correlation is constructed from the ratio of the distribution in which the trigger and associated particles are taken from the same event to the distribution in which the trigger and associated particles are taken from two different events. These two distributions are referred to as the “same-event” (S) or “foreground” distribution and the “mixed-event” or “background” (B) distribution, respectively, and the ratio is written as:

$$\begin{aligned} C(\Delta \eta ,\Delta \phi ) =\frac{S(\Delta \phi ,\Delta \eta )}{B(\Delta \phi ,\Delta \eta )}. \end{aligned}$$

The same-event distribution includes both the physical correlations and correlations arising from detector acceptance effects. On the other hand, the mixed-event distribution reflects only the effects of detector inefficiencies and non-uniformity, but contains no physical correlations. To ensure that the acceptance effects in the B distribution match closely those in the S distribution, the B distribution is constructed from particles from two different events that have similar multiplicity and z-vertex. Furthermore, in order to account for the effects of tracking efficiency \(\epsilon (p_{\mathrm{T}},\eta )\), and fakes \(f(p_{\mathrm{T}},\eta )\), each pair is weighted by

$$\begin{aligned} w_{a,b} = \frac{ (1-f(p_{\mathrm {T}}^{a},\eta ^a))(1-f(p_{\mathrm {T}}^{b},\eta ^b))}{\epsilon (p_{\mathrm {T}}^{a},\eta ^a)\epsilon (p_{\mathrm {T}}^{b},\eta ^b)} \end{aligned}$$

for S and B. In the ratio C, the acceptance effects largely cancel out and only the physical correlations remain [49]. Typically, the two-particle correlations are used only to study the shape of the correlations in \(\Delta \phi \), and are conveniently normalized. In this paper, the normalization of \(C(\Delta \eta ,\Delta \phi )\) is chosen such that the \(\Delta \phi \)-averaged value of \(C(\Delta \eta ,\Delta \phi )\) is unity for \(|\Delta \eta |>2\).

Figure 2 shows \(C(\Delta \eta ,\Delta \phi )\) for several centrality intervals for \(2~<~p_{\mathrm {T}}^{a,b} <~3~\mathrm{GeV}\). In all cases a peak is seen in the correlation at \((\Delta \eta ,\Delta \phi )\sim (0,0)\). This “near-side” peak arises from a range of sources including resonance decays, Hanbury Brown and Twiss (HBT) correlations [50] and jet fragmentation [51]. The long-range (large \(|\Delta \eta |\)) correlations are the result of the global anisotropy of the event and are the focus of the study in this paper.

Fig. 2
figure 2

Two-particle correlation functions \(C(\Delta \eta ,\Delta \phi )\) in 5.02 \(\mathrm{TeV}\) Pb+Pb collisions for \(2~<~p_{\mathrm {T}}^{a,b} <~3~\mathrm{GeV}\). The left, middle and right panels correspond to the 0–5%, 30–40% and 60–70% centrality classes, respectively. The distributions are truncated to suppress the peak at \(\Delta \eta = \Delta \phi = 0\) to show the long-range correlations in greater detail

To investigate the \(\Delta \phi \) dependence of the long-range (\(|\Delta \eta |~>~2\)) correlation in more detail, the projection on to the \(\Delta \phi \) axis is constructed as follows:

$$\begin{aligned} C(\Delta \phi ) = \frac{\int _{2}^{5} \text {d}|\Delta \eta | \; S(\Delta \phi ,|\Delta \eta |)}{\int _{2}^{5} \text {d}|\Delta \eta | \; B(\Delta \phi ,|\Delta \eta |)}\equiv \frac{S(\Delta \phi )}{B(\Delta \phi )}. \end{aligned}$$

The \(|\Delta \eta |~>~2\) requirement is imposed to reject the near-side jet peak and focus on the long-range features of the correlation functions.

In a similar fashion to the single-particle distribution (Eq. (1)), the 2PC can be expanded as a Fourier series:

$$\begin{aligned} C(\Delta \phi ) =C_{0}\left( 1+\Sigma _{n=1}^{\infty }v_{n,n} (p_{\mathrm {T}}^{a},p_{\mathrm {T}}^{b}) \cos (n\Delta \phi )\right) , \end{aligned}$$
(2)

where the \(v_{n,n}\) are the Fourier coefficients of the 2PC, and \(C_0\) is its average value. If the two-particle distribution is simply the product of two single-particle distributions, then it can be shown that the Fourier coefficients of the 2PC factorize as [49]:

$$\begin{aligned} v_{n,n} (p_{\mathrm {T}}^{a},p_{\mathrm {T}}^{b})=v_{n} (p_{\mathrm {T}}^{a})v_{n} (p_{\mathrm {T}}^{b}).\ \end{aligned}$$
(3)

In Ref. [13] it was demonstrated that the factorization of \(v_{n,n}\), given by Eq. (3), is valid in central and mid-central Pb+Pb collisions at \(\sqrt{s_{NN}}\) = 2.76 GeV as long as one of the correlated particles is from a low \(p_{\mathrm{T}}\)range. A breakdown of the factorization is expected when the anisotropy does not arise from flow, e.g. in peripheral collisions at high \(p_{\mathrm{T}}\). The factorization is also expected to break when the \(\eta \) separation between the particles is small, and short-range correlations dominate [13]. However, the \(|\Delta \eta |~>~2\) requirement removes most such short-range correlations. In the phase-space region where Eq. (3) holds, the \(v_{n}\) (\(p_{\mathrm {T}}^{b}\)) can be evaluated from the measured \(v_{n,n}\) as:

$$\begin{aligned} v_{n} (p_{\mathrm {T}}^{b})=\frac{v_{n,n} (p_{\mathrm {T}}^{a},p_{\mathrm {T}}^{b})}{v_{n} (p_{\mathrm {T}}^{a})}= \frac{v_{n,n} (p_{\mathrm {T}}^{a},p_{\mathrm {T}}^{b})}{\sqrt{v_{n,n} (p_{\mathrm {T}}^{a},p_{\mathrm {T}}^{a})}}, \end{aligned}$$
(4)

where \(v_{n,n} (p_{\mathrm {T}}^{a},p_{\mathrm {T}}^{a})=v_{n} ^2(p_{\mathrm {T}}^{a})\) is used in the denominator. In this analysis, for most of the 2PC results the \(v_{n}\) (\(p_{\mathrm {T}}^{b}\)) will be evaluated using Eq. (4) with \(0.5~<~p_{\mathrm {T}}^{a} <~5.0~\mathrm{GeV}\). The lower cut-off of 0.5 \(\mathrm{GeV}\) on \(p_{\mathrm {T}}^{a}\) is the lower limit of \(p_{\mathrm{T}}\)measurements in this paper. The upper cut-off on \(p_{\mathrm {T}}^{a}\) is chosen to exclude high-\(p_{\mathrm{T}}\) particles, which predominantly come from jets and are not expected to obey Eq. (4).

Figure 3 shows one-dimensional 2PCs as a function of \(\Delta \phi \) for \(2<p_{\mathrm {T}}^{a,b} <3~\mathrm{GeV}\) and for several different centrality intervals. The correlations are normalized to have a mean value (\(C_{0}\) in Eq. (2)) of 1.0. The continuous line in Fig. 3 is a Fourier fit to the correlation (Eq. (2)) that includes harmonics up to \(n=6\). The contribution of the individual \(v_{n,n}\) are also shown. The modulation in the correlation about its mean value is smallest in the most central events (top left panel) and increases towards mid-central events, reaching a maximum in the 40–50% centrality interval and then decreases. In central collisions, the \(v_{2,2}\)\(v_{4,4}\) are of comparable magnitude. But for other centralities, where the average collision geometry is elongated, the \(v_{2,2}\) is significantly larger than the other \(v_{n,n}\) for \(n\ge 3\). In the central events the “away-side” peak at \(\Delta \phi \sim \pi \) is also much broader because all the significant harmonics are of similar magnitude, while in mid-central events the near-side and away-side peaks are quite symmetric as the \(v_{2,2}\) dominates. In central and mid-central events, the near-side peak is larger than the away-side peak. However, for the 60–70% and more peripheral centralities, the away-side peak becomes larger due to the presence of a large negative \(v_{1,1}\) component. This negative \(v_{1,1}\) component in the peripheral 2PCs arises largely from dijets: while the near-side jet peak is rejected by the \(|\Delta \eta |~>~2\) requirement, the “away-side jet” correlation that arises from back-to-back jets and contributes at \(\Delta \phi =\pi \), cannot be rejected entirely as its position varies in \(|\Delta \eta |\) from event to event. In the peripheral multiplicity intervals, the away-side jet significantly affects the 2PC. It produces a large negative \(v_{1,1}\) and also affects the other harmonics by adding alternately positive and negative contributions to \(v_{n,n}\) harmonics of even and odd order, respectively. In peripheral events the \(v_{n,n}\) are strongly biased by dijets especially at higher \(p_{\mathrm{T}}\). The presence of the jets also results in the breakdown of the factorization relation (Eq. (3)).

Fig. 3
figure 3

One-dimensional two-particle correlation functions \(C(\Delta \phi )\) in 5.02 \(\mathrm{TeV}\) Pb+Pb collisions for \(2~<~p_{\mathrm {T}}^{a,b} ~<~3~\mathrm{GeV}\) (points). The solid-black line indicates a fit to Eq. (2) containing harmonics \(v_{n,n}\) up to \(n=6\). The dashed grey line shows the contribution of the \(v_{1,1}\). The contributions of the \(v_{2,2}\)\(v_{6,6}\) are indicated by the coloured lines (\(v_{2,2}\)-red, \(v_{3,3}\)-blue, \(v_{4,4}\)-magenta, \(v_{5,5}\)-orange, \(v_{6,6}\)-green), and can be identified by the number of peaks that they have. Each panel corresponds to a different centrality class. The y-axis range for the different panels is different

4.2 Scalar product and event plane analysis

The SP method was introduced by the STAR Collaboration [25] and is further discussed in Ref. [17]. The SP method is very similar to the Event Plane method (EP) widely used in earlier analyses [12, 13]. It is superior to the EP as \(v_{n}\{\mathrm {SP}\}\) is an estimator of \(\sqrt{\langle v_{n} ^{2}\rangle }\), independent of the detector resolution and acceptance, whereas \(v_{n}\{\mathrm {EP}\}\) produces a detector-dependent estimate of \(v_{n}\) that lies between \(\langle v_{n} \rangle \) and \(\sqrt{\langle v_{n} ^{2}\rangle }\) [28].

Both the SP and EP method use flow vectors \(Q_{n}\) and \(q_{n,j}\) defined as:

$$\begin{aligned} Q_{n} = |Q_{n}|\text {e}^{\text {i}n\Psi _{n}} = \frac{1}{M}\sum _{j=1,M}q_{n,j} = \frac{1}{M}\sum _{j=1,M}w_j \text {e}^{\text {i}n\phi _{j}}, \end{aligned}$$
(5)

where the sum runs over M particles in a single event. The \(\phi _{j}\) is the particle azimuthal angle and n is the harmonic order. In this analysis the flow vectors are established separately for the two sides of the FCal and are denoted \(Q_n^{N|P}\), where the N and P correspond to \(\eta < 0\) and \(\eta >0\) sides, respectively. The sum in Eq. (5) in this case runs over the calorimeter towers of approximate granularity \(\Delta \eta \times \Delta \phi = 0.1 \times 0.1\) and the weights \(w_i\) are the transverse energies \(E_{\mathrm{T}}\) measured in the FCal towers. The flow vectors are also calculated using charged-particle tracks. In this case the sum in Eq. (5) is over tracks and \(w_{j}\) is obtained as the MC tracking weight (\((1-f)/\epsilon \)) multiplied by a factor that depends on azimuthal angle to correct for non-uniformity in the azimuthal-angle distribution of reconstructed tracks. This latter factor is obtained run-by-run from the data as the average track multiplicity in one \(\eta \) slice of 0.1 divided by the multiplicity in the narrow \(\Delta \eta \times \Delta \phi = 0.1 \times 0.1\) interval.

The main idea of the SP method is to correlate single-track unit flow vectors with the flow vector of all particles measured in the FCal region (\(3.2< |\eta | < 4.9\)). Therefore, the SP method differs from the two-particle correlation method, in which each single track is correlated with all tracks of \(|\Delta \eta |>2\) in the event. The values of \(v_{n} \) in this analysis are obtained as:

$$\begin{aligned} v_{n}\{\mathrm {SP}\}= & {} Re \frac{\langle q_{n,j}Q_{n}^{N|P*} \rangle }{\sqrt{\langle Q_{n}^{N}Q_{n}^{P*}\rangle }}\nonumber \\= & {} \frac{\langle |q_{n,j}||Q_{n}^{N|P}|\cos [n(\phi _{j} - \Psi _{n}^{N|P})] \rangle }{\sqrt{\langle |Q_{n}^{N}||Q_{n}^{P}|\cos [n(\Psi _{n}^{N} - \Psi _{n}^{P})] \rangle }}, \end{aligned}$$
(6)

where \(q_{n,j}\) is the flow vector obtained for a small (\(\eta , p_{\mathrm{T}}\)) interval (typically 0.1 in \(\eta \) and 0.1 \(\mathrm{GeV}\) in \(p_{\mathrm{T}}\)below 5 \(\mathrm{GeV}\) and 1 \(\mathrm{GeV}\) at higher \(p_{\mathrm{T}}\)) using tracks, \(Q_{n}^{N|P}\) is the flow vector obtained using either the N or P side of the FCal, chosen so that the \(\eta \) gap between the \(q_{n,j}\) and \(Q_{n}\) is maximized, the * denotes complex conjugation, the \(\Psi _{n}\) are estimates of the \(n^{\text {th}}\)-order reaction-plane angles (Eq. (1)) and the angular brackets indicate an average over all events. In the last line of Eq. (6) it is assumed that the sine terms disappear, as required from symmetry. The correction factor, \(1/\sqrt{\langle Q_{n}^{N}Q_{n}^{P*}\rangle }\), (Eq. (6)) depends on the harmonic order and \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) as shown in the left panel of Fig. 4. The event-plane angles, \(\Psi _{n}\), and the \(Q_{n}\) vectors, both measured in the FCal, may be biased due to non-uniform detector response. As \(\Psi _{n}\) varies randomly from event to event, its distribution should be uniform, and the components of the \(Q_n\) vector, \(Q_{n,x} = |Q_{n,}|\mathrm {cos}(\Psi _{n})\) and \(Q_{n,y} = |Q_{n}|\mathrm {sin}(\Psi _{n})\), should be zero when averaged over many events. This is achieved by the following procedure. In its first step, non-zero offsets of the mean of raw flow vector coordinates are removed for each run: \(Q_{n,i} =Q^{\mathrm {raw}}_{n,i} - \langle Q^{\mathrm {raw}}_{n,i}\rangle \) where \(i=x,y\) and \(\langle Q^{\mathrm {raw}}_{n}\rangle \) is the mean calculated for each run. However, even after this correction, residual higher-order non-uniformities persist, indicated by non-zero values of \(\langle Q_{n,x}Q_{n,y}\rangle \). These are removed by rotating the \(Q_{n}\) vector so that the corrected \(Q_{n}\) vector has no skew (\(\langle Q_{n,x}^{2}\rangle = \langle Q_{n,y}^{2}\rangle \); \(\langle Q_{n,x} Q_{n,y}\rangle = 0\)) and the distributions of the resulting EP angles, \(\Psi _{n}\), are uniform [52].

Fig. 4
figure 4

The dependence of the correction factor in the SP method, \(\sqrt{\langle Q_{n}^{N}Q_{n}^{P*}\rangle }\) (left panel), and EP method, \(\sqrt{\Big \langle \frac{Q_{n}^{N} Q_{n}^{P*}}{|Q_{n}^{N} ||Q_{n}^{P}|}\Big \rangle }\) (right panel), for all measured harmonics as a function of \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) binned according to the centrality bins definition

In the Event Plane analysis the reference \(Q_{n}\) vectors are normalized to unity, \(Q_n^{N|P} \rightarrow Q_n^{N|P}/|Q_n^{N|P}|\), before using them in Eq. (6). So the \(v_{n}\) estimate is obtained as:

$$\begin{aligned} v_{n} \{\mathrm {EP}\} = Re \frac{\Big \langle q_{n,j}\frac{Q_{n}^{N|P*}}{|Q_{n}^{N|P}|}\Big \rangle }{\sqrt{\Big \langle \frac{Q_{n}^{N}}{|Q_{n}^{N}|} \frac{Q_{n}^{P*}}{|Q_{n}^{P}|} \Big \rangle }} = \frac{\langle \cos [n(\phi _{j} - \Psi _{n}^{N|P})] \rangle }{\sqrt{\langle \cos [n(\Psi _{n}^{N} - \Psi _{n}^{P})] \rangle }}. \end{aligned}$$
(7)

The denominator of Eq. (7), shown in the right panel of Fig. 4, can be thought of as a resolution. It is distinct for each harmonic and depends on \(\Sigma E_\mathrm {T}^\mathrm {FCal}\).

In this analysis the EP method is used only for the purpose of a direct comparison with the results obtained in Run 1 [13], in which only the EP method was used.

The analysis is performed in intervals of centrality. The \(v_{n}\) values are obtained in narrow bins of \(p_{\mathrm{T}}\) and \(\eta \), which are summed, taking into account tracking efficiency and fake rate, to obtain the integrated results.

A detailed study based on a HIJING [39] Monte Carlo sample showed a difference for the most central events between the \(v_{n}\) obtained with generated particles and the \(v_{n}\) obtained with reconstructed tracks (the “MC closure test”). The discrepancies are due to the presence of fake tracks, which at low \(p_{\mathrm{T}}\) distort the \(v_{n}\) measurements, and a tracking inefficiency in the event-plane direction due to increased detector occupancy resulting from the flow phenomenon itself, which lowers the measured \(v_{n}\) values. The study based on the \(d_0\) distribution also showed that the fake-track rates are overestimated in MC simulation as compared to the data. This disagreement was removed by weighting MC tracks so that the \(d_0\)-distribution tails (\(2< |d_0| < 10 \, \hbox {mm}\)) match those in data, following the procedure described in Ref. [53]. It was observed that the contribution of fakes to the “MC non-closure” is significant only for events with centrality < 30% and at low \(p_{\mathrm{T}}\), which is the region where the fake rate is the largest. In this modified MC sample, the relative differences between values of the \(v_{n}\) measured with generated particles and reconstructed tracks are used as corrections to account for both effects; the fakes and the \(\Psi _{n}\)-dependent inefficiency. Corrections are at most 5–10%. For example, for \(v_{2}\) in the 0–5% centrality interval, the correction is as large as 7% at low \(p_{\mathrm{T}}\)and becomes negligible above \(p_{\mathrm{T}}> 2~\mathrm{GeV}\). Corrections of a similar magnitude are obtained for higher-order harmonics.

5 Systematic uncertainties

The systematic uncertainties of the measured \(v_{n}\) are evaluated by varying several aspects of the analysis. As the EP and SP results are subject to the same uncertainty sources, the uncertainty values are of the similar magnitude and are not discussed separately. Similarly, some of the uncertainty sources are common to the SP/EP and the 2PC methods and are discussed together. The uncertainties for two representative \(p_{\mathrm{T}}\) intervals are summarized in the Tables 3 and 4 for the 2PC and SP/EP methods, respectively. In the discussion below, other \(p_{\mathrm{T}}\) ranges are referred to if uncertainties are significantly higher than in the \(p_{\mathrm{T}}\) ranges shown in the tables. The following sources of uncertainty are considered:

  • Track selection: The tracking selection requirements control the relative contribution of genuine charged particles and fake tracks entering the analysis. The stability of the results to the track selection is evaluated by varying the requirements imposed on the reconstructed tracks. Two sets of variations are used. In the first case the required number of pixel and SCT hits on the reconstructed track are relaxed to one and six, respectively. Additionally, the requirements on the transverse and longitudinal impact parameters of the track are relaxed to 1.5 mm. In the second case, the track selection is based on requirements used for the baseline measurement, but the transverse and longitudinal impact parameters of the track are restricted to 0.5 mm. For each variation, the entire analysis is repeated including the evaluation of the corresponding efficiencies and fake rates. The fake rate is largest at the lowest \(p_{\mathrm{T}}\)(0.5 GeV) and for the most central events, and consequently the variation in the \(v_{n}\) values obtained from this procedure is largest, typically 10%, in this region of phase space.

  • Tracking efficiency: As mentioned above, the tracks are weighted by a factor \((1-f)/\epsilon (p_{\mathrm{T}},\eta )\) when calculating the \(v_{n}\) to account for the effects of the tracking efficiency. Uncertainties in the efficiency, resulting e.g. from an uncertainty in the amount of detector material, need to be propagated into the measured \(v_{n}\) [54]. This uncertainty is evaluated by varying the efficiency up and down within its uncertainties in a \(p_{\mathrm{T}}\)-dependent manner and re-evaluating the \(v_{n}\). This contribution to the overall uncertainty is very small and amounts to less than 1% on average. This is because the change of efficiency largely cancels out in the differential \(v_{n} (p_{\mathrm{T}})\) measurement, and for \(v_{n}\) integrated over \(p_{\mathrm{T}}\)the low-\(p_{\mathrm{T}}\)particles dominate the measurement. It does not change significantly with centrality nor with the order of harmonics.

  • Centrality determination: An uncertainty in the flow harmonics comes from the uncertainty in the fraction of the total inelastic cross-section accepted by the trigger and the event selection criteria, which was estimated to be at a level of 1%. The \(v_{n}\) uncertainty is evaluated by repeating the analysis with the modified centrality selections on the \(\Sigma E_\mathrm {T}^\mathrm {FCal}\) distribution shown in Fig. 1, that give the ± 1% uncertainty in the sampled fraction of the cross-section [12]. The changes in the \(v_{n}\) are largest in the peripheral-centrality intervals, for which the bin definitions are significantly changed when remapping the centralities. For \(v_{2}\), a change of \(\sim \)0.8% (2PC) and \(\sim \)1% (SP) is also observed in the most central events. This is because the \(v_{2}\) changes rapidly with centrality in central events, so slight variations in the centrality definition result in significant change in \(v_{2}\). For \(v_{3}\) this uncertainty varies from less than 0.5% over the 0–50% centrality range to \(\sim \)5% in the 70–80% centrality interval. For the higher-order harmonics, \(n>3\), the uncertainty is less than 0.5% over the 0–50% centrality range and increases to about 2% for more peripheral bins. The variation in the \(v_{n}\) when using these alternative centrality definitions is taken as a systematic uncertainty. To limit the statistical instability of \(v_{6}\) and \(v_{7}\) in uncertainty estimation, the variations for this measurement were determined over a wide range of \(p_{\mathrm{T}}\)= 0.5–60 \(\mathrm{GeV}\).

  • MC corrections: To assess the uncertainty related to the MC corrections the closure test is repeated with the two selections of tracks described in the “track selection” paragraph. Differences between the correction factors obtained with loose, nominal and tight tracking selections are compared. The difference between them is largest at low \(p_{\mathrm{T}}\)and central events and amounts typically to a few percent. It is negligibly small above 2 \(\mathrm{GeV}\). The larger of the two differences (between the nominal and loose tracking selections) is used as an uncertainty estimate.

  • Residual sine term: The ability of the detector to measure small \(v_{n}\) signals can be quantified by comparing the value of the \(v_{n}\) calculated as the real part of the flow vector product (SP) in Eq. (6) with its imaginary part. The ratio \(Im(SP)/v_{n} \) is taken as a contribution to the systematic uncertainty. The contribution from this source is \(\sim \)1% in most of the phase space, while for the higher harmonics (\(n=6, 7\)) it can reach 20% in the most central collisions. This uncertainty is only relevant for the \(v_{n}\) values measured by the EP and SP methods.

  • Variation of FCal acceptance in the \(Q_n^{N|P}\) estimation: In order to quantify an uncertainty arising from the FCal acceptance in the \(Q_n^{N|P}\) estimation, \(v_{n}\) harmonics are compared for two distinct FCal regions \(3.2<|\eta |<4\) and \(4<|\eta |<4.8\) used for the determination of the reference flow vector, \(Q_n\). The differences between the \(v_{n}\) values are treated as the systematic uncertainty, which, similarly to the \(\eta \) symmetry (next paragraph), quantifies the ability of the detector to measure small signals. Accordingly, this contribution is small (\(\sim \) 1%) for \(v_{2}\) and \(v_{3}\) and starts growing for higher-order harmonics, reaching about 27% for \(v_{7}\). This uncertainty is only relevant to the \(v_{n}\) values measured by the EP and SP methods.

  • Detector non-uniformity: Due to the symmetry of the Pb+Pb collision system the \(v_{n} (\eta )\) are expected to be on average symmetric in rapidity. Any difference between the event-averaged \(v_{n}\) at \(\pm \eta \) arises from residual detector non-uniformity. The difference between the \(v_{n}\) values measured in opposite hemispheres is treated as the systematic uncertainty quantifying non-perfect detector performance. This uncertainty is in general very low (less than 1%) except for high-order harmonics: \(v_{5}\) and \(v_{6}\) at high \(p_{\mathrm{T}}\)and \(v_{7}\) at all \(p_{\mathrm{T}}\). This uncertainty only contributes to the \(v_{n}\) values measured by the EP and SP methods. For the 2PC method, the residual non-uniformity is estimated by varying the event-mixing procedure.

  • Event-mixing: As explained in Sect. 4.1, the 2PC analysis uses the event-mixing technique to estimate and correct for the detector-acceptance effects. Potential systematic uncertainties in the \(v_{n}\) due to the residual pair-acceptance effects, which were not removed by the mixed events, are evaluated by varying the multiplicity and z-vertex matching criteria used to make the mixed-event distributions, following Ref. [13]. The resulting uncertainty for \(v_{2}\)\(v_{5}\) is between 1–3%, and for \(v_{6}\) is between 4–8% for most of the centrality and \(p_{\mathrm{T}}\) ranges measured in this paper. However, the uncertainties for \(v_{4}\)\(v_{6}\) are significantly larger for \(p_{\mathrm{T}}\ {<} ~0.7~\mathrm{GeV}\), where the \(v_{n}\) signals are quite small and very susceptible to acceptance effects, and for \(v_{6}\) are correlated with statistical uncertainties. The uncertainties are also significantly larger for \(p_{\mathrm{T}} {>} ~10~\mathrm{GeV}\), where they are difficult to determine due to large statistical uncertainties in the measurements.

Table 3 The systematic uncertainties associated with the 2PC \(v_{n}\) measurements for selected intervals of \(p_{\mathrm{T}}\)and for 5–10% and 40–50% centrality bins. The contributions are expressed in %. The total systematic uncertainty is obtained by adding the contribution of the individual sources in quadrature
Table 4 The systematic uncertainties associated with the SP and EP (in parentheses) \(v_{n}\) measurements for \(v_{n}\) in 5–10% and 40–50% centrality bins. The contributions are expressed in %. The total systematic uncertainty is obtained by adding the contribution of the individual sources in quadrature

6 Results

Fig. 5
figure 5

The \(v_{n}\) obtained with the SP method as a function of transverse momentum \(p_{\mathrm{T}}\), integrated over \(|\eta | < 2.5\) in 11 centrality intervals, from the most central at the top left panel to the most peripheral at the bottom right panel. Results are averaged over the intervals indicated by horizontal error bars. The vertical error bars indicate statistical uncertainties; the shaded boxes indicate systematic uncertainties

Fig. 6
figure 6

The \(v_{n}\) values obtained with the 2PC method as a function of \(p_{\mathrm {T}}^{b}\) for \(0.5~<~p_{\mathrm {T}}^{a} <~5~\mathrm{GeV}\). Each panel represents a different centrality interval. The vertical error bars indicate statistical uncertainties. The shaded bands indicate systematic uncertainties

6.1 The \(p_{\mathrm{T}}\) dependence of \(v_{n}\)

Figures 5 and 6 show the \(v_{n}\) obtained from the SP and 2PC methods, respectively, as a function of \(p_{\mathrm{T}}\) for several centrality intervals. For the SP method the \(v_{2}\)\(v_{5}\) harmonics are also shown for the 0–0.1% and 0–1% ultra-central collisions. The SP results are integrated over the pseudorapidity \(|\eta |<2.5\) and the 2PC results are obtained with \(0.5~<~p_{\mathrm {T}}^{a} <~5~\mathrm{GeV}\) and for \(2~<~|\Delta \eta |~<~5\). The \(v_{n}\) values show a similar \(p_{\mathrm{T}}\)dependence across all centralities: a nearly linear rise to about 2 \(\mathrm{GeV}\), followed by a gradual increase to reach a maximum around 2–4 \(\mathrm{GeV}\) and a gradual fall at higher \(p_{\mathrm{T}}\). However, significant \(v_{n}\) values persist at high \(p_{\mathrm{T}}\)(\(\sim \)20 \(\mathrm{GeV}\)). The \(v_{2}\) is positive even at the highest measured \(p_{\mathrm{T}}\)of 60 \(\mathrm{GeV}\) (Fig. 5). This indicates the parton energy loss in the created medium [30]. Such elliptic flow is expected due to path-length dependence of the energy loss of high-\(p_{\mathrm{T}}\) partons traversing the hot and dense medium. In peripheral events, at the highest \(p_{\mathrm{T}}\), the 2PC and SP \(v_{2}\) values again show an increasing trend due to the increasing influence of the away-side jet. The increased \(v_{2}\) in peripheral collisions at high-\(p_{\mathrm{T}}\) is accompanied by reduced values of \(v_{3}\) and increased values of \(v_{4}\), which is characteristic of a large away-side peak, as described in Sect. 4.1. This is most clearly seen in the 70–80% centrality interval.

The \(v_{2}\) varies significantly with centrality, reflecting a change in the shape of the average initial collision geometry, from nearly circular in ultra-central collisions to an almond shape in peripheral events. The higher harmonics do not show similar behaviour, as neither higher-order eccentricities nor the fluctuations vary so significantly with the centrality. The \(v_{2}\) is dominant at all centralities, except for the most central collisions interval where, at intermediate \(p_{\mathrm{T}}\), \(v_{3}\) and \(v_{4}\) become larger than \(v_{2}\), indicating that the dominant source of observed flow comes from the initial geometry fluctuations. This change in the \(v_{n}\) ordering is even more pronounced in the 1% and 0.1% ultra-central collisions measured using the SP method, which shows that, in the \(p_{\mathrm{T}}\)region around the \(v_{n}\) peak, \(v_{3}> v_{4} > v_{2} \approx v_{5} \). The \(v_{4}\), similarly to \(v_{2}\), exhibits an increase beyond \(p_{\mathrm{T}}\sim 10\ \mathrm{GeV}\), which can be attributed to the presence of the events with dijets in the data. In the SP measurement the \(v_{7}\) results are also presented. The characteristics of \(v_{7}\) are similar to the other high-order harmonics, but the values are smaller and significant, given the uncertainties, only in central and mid-central collisions and for the \(p_{\mathrm{T}}\) range of 2–6 \(\mathrm{GeV}\).

6.1.1 The scalar product and event plane methods comparison

Figure 7 compares the \(v_{n}\) values measured with the EP and SP methods as a function of \(p_{\mathrm{T}}\)and \(N_{\mathrm {part}}\) for the integrated \(p_{\mathrm{T}}\) range of \(0.5<p_{\mathrm{T}}<60~\mathrm{GeV}\). A small difference is seen between the \(v_{2}\) values measured with the two methods. The difference is largest in mid-central events: about 3% in the 20–30% and 40–50% centrality intervals, about 1% in the 0–5% most central collisions and negligible in peripheral collisions. This difference is expected according to Ref. [28] as the SP method measures \(\sqrt{\langle v_{n} ^2\rangle }\) while the EP method measures values between \(\langle v_{n} \rangle \) and \(\sqrt{\langle v_{n} ^2\rangle }\), with the former value attained in the limit of a small correction factor (the inverse of the denominator in Eq. (7)) and approaching the latter when the correction factor is large. In the most central and peripheral events, where the correction is large for the second-order harmonic, the EP \(v_{2}\) values are closer to the SP ones, while for the mid-central events where the correction is small, the EP \(v_{2}\) values are systematically lower than the SP \(v_{2}\) values. Higher-order EP and SP \(v_{n}\) harmonics are consistent with each other.

Fig. 7
figure 7

Comparison of the \(v_{n}\) obtained with EP and SP methods as a function of \(p_{\mathrm{T}}\) in three centrality bins: 0–5%, 20–30% and 40–50%. The right bottom panel shows the \(v_{n}\) as a function of \(N_{\mathrm {part}}\), integrated over \(0.5~<~p_{\mathrm{T}}~<~60~\mathrm{GeV}\). The correspondence of \(N_{\mathrm {part}}\) to centrality intervals is provided in Table 2. In the inset the \(v_{6}\) and \(v_{7}\) integrated over \(0.5~<~p_{\mathrm{T}}~<~60~\mathrm{GeV}\) are shown with adjusted scale. For the \(v_{n} (p_{\mathrm{T}})\) comparisons, the results are averaged over the intervals indicated by horizontal error bars. The vertical error bars indicate the quadrature sum of statistical and systematical uncertainties

Fig. 8
figure 8

Comparison of the \(v_{n}\) obtained with 2PC and SP methods as a function of \(p_{\mathrm{T}}\). Each panel shows the comparison for a different order harmonic. The comparisons are shown for three different centrality intervals: 0–5%, 20–30% and 40–50%. The vertical error bars indicate statistical uncertainties only

6.1.2 The scalar product and two-particle correlation methods comparison

A comparison of the SP and 2PC results is presented in Fig. 8. In general, results from the two methods are quite consistent. There is a significant difference in \(v_{2}\) from the two methods in the phase-space region \(p_{\mathrm{T}}<5~\mathrm{GeV}\), 0–5% centrality. This difference decreases considerably for 20–30% mid-central events, where the \(v_{2}\) values match within 2–5% up to \(p_{\mathrm{T}}\sim \)10 \(\mathrm{GeV}\). For \(v_{3}\)\(v_{5}\), where there are enough events for a clear comparison, the \(v_{n}\) values match within \(\sim \)4% for \(p_{\mathrm{T}}<4~\mathrm{GeV}\) for the three centrality intervals shown in Fig. 8. In principle, both the SP and 2PC methods measure \(\sqrt{\langle v_{n} ^2\rangle }\) and the flow harmonics measured by the two methods should be identical. However, a breakdown of factorization (Eq. (3)) results in systematic differences in the flow harmonics measurement. Such factorization breakdown has been observed to be significant for \(v_{2}\) in central events [55], and in general for all \(v_{n}\) at higher \(p_{\mathrm{T}}\), and is the leading source of disagreement between the 2PC and SP results. Furthermore, in the 2PC method the \(\Delta \eta \) gap between the reference and associated particles is chosen to be \(|\Delta \eta |~>~2\), while in the SP method, where the reference flow is measured in the FCal, the minimum gap between the tracks and the FCal is 3.2 units in \(\eta \). The presence of longitudinal-flow fluctuations, in which the event-plane angle can change with \(\eta \), can result in different \(v_{n}\) values depending on the \(\eta \) range where the reference flow is measured [27, 56]. This effect is also found to be larger in central events and relatively smaller in mid-central events [56]. These effects can further contribute to the observed difference between the SP and 2PC \(v_{n}\) values.

6.1.3 Comparison to Pb+Pb results at \(\sqrt{s_{_\text {NN}}}\) = 2.76 \(\mathrm{TeV}\)

Figure 9 shows a comparison of the \(v_{n}\) measured in the present analysis at \(\sqrt{s_{_\text {NN}}}=5.02\) \(\mathrm{TeV}\) with the corresponding measurements at \(\sqrt{s_{_\text {NN}}}=2.76\) \(\mathrm{TeV}\) for harmonics \(v_{2}\) to \(v_{6}\) obtained using the 2PC method [13]. The comparisons are shown for two centralities: a central interval of 0–5% and a mid-central interval of 20–30%. Figure 10 shows a similar comparison of results obtained using the EP method for 0–5%, 20–30% and 40–50% centrality bins. The \(v_{n}\) at the two energies are quite similar and almost consistent throughout within systematic and statistical uncertainties, even though the MC non-closure correction was applied only in the \(\sqrt{s_{_\text {NN}}}=5.02\) \(\mathrm{TeV}\) measurement. These results are consistent with the recent ALICE measurements comparing the measurement of \(v_{n}\) at the two collision energies [29].

Fig. 9
figure 9

Comparisons of the 2PC \(v_{n}\) harmonics measured at \(\sqrt{s_{_\text {NN}}}\) = 2.76 \(\mathrm{TeV}\) (Run 1) and at \(\sqrt{s_{_\text {NN}}}\) = 5.02 \(\mathrm{TeV}\) (Run 2). The results are plotted as a function of \(p_{\mathrm {T}}^{b}\) for \(1~<~p_{\mathrm {T}}^{a} ~<~2~\mathrm{GeV}\) for two centralities: 0–5% and 20–30%. Each panel corresponds to a different harmonic. Results are averaged over the intervals indicated by horizontal error bars. The vertical error bars indicate statistical uncertainties

Fig. 10
figure 10

Comparison of the \(v_{n}\) obtained with EP method using Run 1 and Run 2 data as a function of \(p_{\mathrm{T}}\). The results are shown in three centrality intervals: 0–5%, 20–30% and 40–50%. Results are averaged over the intervals indicated by horizontal error bars. The vertical error bars indicate statistical uncertainties. The shaded areas indicate systematic uncertainties

Fig. 11
figure 11

The \(v_{n}\) as a function of pseudorapidity obtained with the SP method, for transverse momentum ranges: \(0.8<p_{\mathrm{T}}<1~\mathrm{GeV}\) (left column), \(2<p_{\mathrm{T}}<3~\mathrm{GeV}\) (middle column) and \(7< p_{\mathrm{T}}< 60~\mathrm{GeV}\) (right column) and for centrality intervals: 0–0.1% (top row), 0–5%, 10–20%, 30–40% and 60–70% (bottom row). The vertical error bars indicate statistical uncertainties. The shaded boxes indicate systematic uncertainties

6.2 The \(\eta \) dependence of \(v_{n}\)

The pseudorapidity dependence of the \(v_{2}\)\(v_{7}\) obtained from the SP method is shown in Fig. 11 as a function of \(|\eta |\). Benefiting from the symmetry of \(v_{n} (\eta )\) with respect to \(\eta =0\), the \(v_{n}\) pseudorapidity dependence over the full range of pseudorapidity was folded into the \(\eta \) range 0–2.5. The \(\eta \)-dependence is shown for three ranges of transverse momenta \(0.8<p_{\mathrm{T}}<1~\mathrm{GeV}\), \(2<p_{\mathrm{T}}<3~\mathrm{GeV}\) and \(7<p_{\mathrm{T}}<60~\mathrm{GeV}\) and for five centrality intervals 0–0.1%, 0–5%, 10–20%, 30–40% and 60–70%. The strong dependence of flow harmonics on \(p_{\mathrm{T}}\)and centrality shown across different panels (all vertical axes in Fig. 11 have the same range) is discussed in the previous section. On the other hand, no strong pseudorapidity dependence of \(v_{n}\) harmonics is observed. The \(v_{2}\) harmonic in central and mid-central collisions for \(p_{\mathrm{T}}<3~\mathrm{GeV}\) drops only by about 2–4% over the pseudorapidity range \(|\eta |\) = 0–2.5. For peripheral collisions and for high \(p_{\mathrm{T}}>7~\mathrm{GeV}\) a larger decrease of about 10% is observed. The \(v_{3}\) harmonic in central and mid-central collisions over the \(p_{\mathrm{T}}\) range from 2 to 3 \(\mathrm{GeV}\) decreases by about 10% with a larger drop of \(\sim \)15% for peripheral collisions. Similar pseudorapidity dependence is measured for the \(v_{4}\) harmonic in central and mid-central collisions over the \(p_{\mathrm{T}}\) range from 2 to 3 \(\mathrm{GeV}\) where it changes by about 10%, but a larger drop of 25% is observed in peripheral collisions. In other cases, \(v_{n}\) harmonics are almost consistent with a uniform distribution within the statistical and systematic uncertainties. As observed in the earlier measurement of \(v_{n}\) harmonics in 2.76 \(\mathrm{TeV}\) Pb+Pb collisions [19], such a weak \(\eta \) dependence of \(v_{2}\) may be partially attributed to a contribution of “non-flow” short-range two-particle correlations.

6.3 Centrality dependence of \(v_{n}\)

Figure 12 shows the \(N_{\mathrm {part}}\) dependence of \(v_{n}\) integrated over \(|\eta |<2.5\) and for various ranges of \(p_{\mathrm{T}}\)using the SP method. The elliptic flow is the dominant anisotropy, except at the largest \(N_{\mathrm {part}}\) (\(N_{\mathrm {part}}\gtrsim 350\)), which corresponds to the 0-5% most central collisions. For \(p_{\mathrm{T}}<8~\mathrm{GeV}\), a clear dependence on initial geometry can be observed as \(v_{2}\) is highest in mid-central collisions, where this asymmetry is most significant. For \(p_{\mathrm{T}}>8 ~\mathrm{GeV}\), \(v_{2}\) is still the dominant harmonic, and it is non-zero even in peripheral collisions as non-flow effects start to contribute in this region. A hierarchy \(v_{n+1} < v_n\) is observed for \(n=2\)–7 for all ranges of \(p_{\mathrm{T}}\) and all centralities, except for the two bins 0–0.1% and 0–1% of the ultra-central collisions at intermediate \(p_{\mathrm{T}}\), which are characterised by the flow harmonics ordering \(v_{3}> v_{4} > v_{2} \approx v_{5} \).

Fig. 12
figure 12

Integrated \(v_{n}\{\mathrm {SP}\}\) vs. \(N_{\mathrm {part}}\) for six \(p_{\mathrm{T}}\)ranges shown in the panels from lowest \(p_{\mathrm{T}}\)range at the top left to the highest at the bottom middle, measured using the scalar-product method. In the inset in the bottom-right panel the \(v_{6}\) and \(v_{7}\) integrated over \(0.5<p_{\mathrm{T}}<60~\mathrm{GeV}\) are shown with adjusted scale. The correspondence of \(N_{\mathrm {part}}\) to centrality intervals is provided in Table 2. Results are averaged over the intervals indicated by horizontal error bars. The vertical error bars indicate statistical uncertainties. The shaded areas indicate systematic uncertainties

6.4 \(v_{n} (p_{\mathrm{T}})\) scaling

The left panels of Fig. 13 compare the \(p_{\mathrm{T}}\) dependence of \(v_{2}\) and \(v_{3}\) for different centrality intervals obtained from the 2PC method. Two distinct features are observed to change in the shape of the \(v_{n} (p_{\mathrm{T}})\) (for the same n) between the different centrality intervals:

  1. 1.

    Change in the \(v_{n}\)-scale: The overall magnitude of the \(v_{n} \) changes from one centrality to another. This effect is particularly large for \(v_{2}\), and has to do with the changing of the average collision geometry from one centrality to another.

  2. 2.

    Change in the \(p_{\mathrm{T}}\) scale: The \(p_{\mathrm{T}}\) value at which the \(v_{n} \) reaches its maximum also changes systematically from one centrality to another.

In a recent ATLAS paper [45] it was observed that for a given order n, the \(v_{n} (p_{\mathrm{T}})\) in p+Pb and Pb+Pb collisions have a very similar \(p_{\mathrm{T}}\) dependence. In fact, after a scaling along the \(p_{\mathrm{T}}\) axis to account for the difference in \(\langle p_{\mathrm{T}}\rangle \) between p+Pb and Pb+Pb collisions, and an empirical scaling along the \(v_{n} \) (i.e. y-axis) to account for the difference in the collision geometry between the two systems, the scaled \(v_{n} (p_{\mathrm{T}})\) in p+Pb and Pb+Pb collisions were found to be nearly identical. In this section, a check is done to see if such a scaling along the \(p_{\mathrm{T}}\) and \(v_{n} \) axes can yield universal shapes for the \(v_{n} \) across the different centrality classes. Accordingly, the \(v_{n}(p_{\mathrm{T}})\) are scaled along the x- and y-axes to match their shapes across the different centrality intervals. For the matching, the 0–60% centrality interval is chosen as the reference, and the \(v_{n} (p_{\mathrm{T}})\) for the individual 5%-wide centralities are scaled to match best the \(v_{n} (p_{\mathrm{T}})\) in the 0–60% centrality interval over the 0.5–5 \(\mathrm{GeV}\) \(p_{\mathrm{T}}\) range, with the scales along the x- and y-axes treated as fit parameters. The right panels of Fig. 13 show the scaled-\(v_{n}\), for \(n=2\) and 3. Overall, the \(v_{n} (p_{\mathrm{T}})\) shapes match well after the \(p_{\mathrm{T}}\) and \(v_{n} \) scalings.

Fig. 13
figure 13

Left panels: the \(v_{2}\) (\(p_{\mathrm{T}}\)) (top) and \(v_{3}\) (\(p_{\mathrm{T}}\)) (bottom) for different centrality intervals. Right panels: the corresponding scaled-\(v_{n}(p_{\mathrm{T}})\). The error bars indicate statistical and systematic uncertainties added in quadrature and are typically too small to be seen

Figure 14 shows the x (or \(p_{\mathrm{T}}\)) and y-scale factors obtained for \(v_{2}\) and \(v_{3}\) as a function of centrality. It is interesting to note that the \(p_{\mathrm{T}}\)-scale factors are quite comparable between \(v_{2}\) and \(v_{3}\) across most of the centrality ranges. However, in the 0–10% most central events, some significant difference is observed between the two scale factors. This difference could be due to larger jet-bias and factorization-breaking effects in the \(v_{2}\) as compared to \(v_{3}\). The right panel of Fig. 14 shows the y-scale factors for \(v_{2}\) and \(v_{3}\) as well. Their centrality dependence is very different for the two harmonics. This is to be expected as the y-scale factors are mostly indicative of the changing collision geometry, which becomes more and more elliptic from central to mid-central events resulting in a large increase in \(v_{2}\), while \(v_{3}\), which is driven by fluctuations, changes only gradually.

In order to check if the change in \(\langle p_{\mathrm{T}}\rangle \) between the different centralities accounts for the change in the x-scale of the \(v_{n} \), the \(1/\langle p_{\mathrm{T}}\rangle \) for protons and pions, as measured by the ALICE Collaboration [57], are also shown for comparison. While the centrality dependence of these \(1/\langle p_{\mathrm{T}}\rangle \) factors is qualitatively similar to the scale factors for the \(v_{n} \), the relative variation in \(1/\langle p_{\mathrm{T}}\rangle \) is significantly smaller than the variation in the x-scale of the \(v_{n}\), indicating that there are additional effects at play besides the change in the mean \(p_{\mathrm{T}}\). However, whatever the origin of these effects may be, they result in nearly identical scaling factors for \(v_{2}\) and \(v_{3}\).

Fig. 14
figure 14

Left panel: the x-scale factors for the \(v_{2}(p_{\mathrm{T}})\) and \(v_{3}(p_{\mathrm{T}})\) (see text) as a function of the collision centrality. Also shown for comparison are the ALICE 1/\(\langle p_{\mathrm{T}}\rangle \) for positively charged protons and pions (scaled by constant factors of 1.4 \(\mathrm{GeV}\) and 0.57 \(\mathrm{GeV}\), respectively for plotting purposes). Right panel: the y-scale factors. The error bars on the scale factors and on the ALICE data indicate systematic and statistical uncertainties added in quadrature

7 Summary

This paper presents the first ATLAS measurements of azimuthal anisotropy of charged particles in Pb+Pb collisions at 5.02 \(\mathrm{TeV}\) collected during LHC running in 2015. The measurements are based on the Pb+Pb sample of 0.49 \({\mathrm {nb}}^{-1}\) integrated luminosity and are performed with the two-particle correlation, scalar-product and event-plane methods. The azimuthal anisotropy is quantified by the flow harmonics \(v_{2} \)\(v_{6} \) and \(v_{2} \)\(v_{7} \) for measurements based on the 2PC and SP/EP methods, respectively. The flow harmonics are obtained in wide transverse momentum (\(0.5< p_{\mathrm{T}}< ~ 60~\mathrm{GeV}\)), pseudorapidity (\(|\eta |< 2.5\)) and centrality 0–80% ranges. All harmonics show a similar trend in the \(p_{\mathrm{T}}\)dependence; first increasing with \(p_{\mathrm{T}}\)up to a maximum around 3–4 GeV and then decreasing for higher \(p_{\mathrm{T}}\). Significant values of the second-order harmonic \(v_{2}\) persist up to 60 \(\mathrm{GeV}\). The \(v_{2}\) results at high \(p_{\mathrm{T}}\)provide a useful handle for the study of partonic energy loss in the dense medium, and so can improve our understanding of the QGP properties. The values of the flow harmonics decrease strongly with increasing harmonic order for all centralities, except for 0–5% central collisions where a different ordering is seen: \(v_{3}> v_{4} > v_{2} \approx v_{5} \) for \(p_{\mathrm{T}}\)above 1 \(\mathrm{GeV}\), which is indicative of the presence of significant \(v_{2}\) fluctuations in these ultra-central collisions. The elliptic flow signal is strongly dependent on event centrality and is largest in mid-central collisions of 30–50%. The higher-order harmonics show a weak centrality dependence, which is consistent with an anisotropy associated with fluctuations in the initial geometry. After scaling the \(p_{\mathrm{T}}\)and magnitude of the differential elliptic and triangular, \(v_{3}\), flow harmonics at different centralities to best match the reference \(v_{2}\) and \(v_{3}\) for the 0–60% centrality interval over the \(p_{\mathrm{T}}\)range \(0.5~<p_{\mathrm{T}}<~5~\hbox {GeV}\), the shapes of the rescaled harmonics agree well, indicating similarity among \(p_{\mathrm{T}}\)-shapes of flow harmonics evolving from different initial conditions. The \(v_{n}\) coefficients are shown to exhibit a weak \(\eta \)-dependence, irrespective of the harmonic order, centrality and \(p_{\mathrm{T}}\). A weaker pseudorapidity dependence of \(v_{n} (\eta )\) in more central collisions is suggested by data. The results obtained using the EP and SP methods are consistent for harmonics of order \(n\ge 3\). A small, systematic difference is observed for \(v_{2}\), where the values obtained from the SP method are up to 3% larger than the values obtained using the EP method. The 2PC and SP methods give values for \(v_{n}\) that are quite consistent up to \(\sim \)10 \(\mathrm{GeV}\). However, in the most central events the SP method gives systematically larger values for \(v_{2}\) for \(p_{\mathrm{T}}> 2~\mathrm{GeV}\). Comparisons with measurements in Pb+Pb collisions at \(\sqrt{s_{_\text {NN}}}=2.76~\mathrm{TeV}\) show that the \(p_{\mathrm{T}}\)dependence of the \(v_{n}\) shows no change from \(\sqrt{s_{_\text {NN}}}=2.76~\mathrm{TeV}\) to \(\sqrt{s_{_\text {NN}}}=5.02~\mathrm{TeV}\).

The set of results on flow harmonics presented in this paper provides a tool for studies of the underlying mechanism leading to the large azimuthal anisotropy observed in the Pb+Pb collision system, and can be used to constrain the theoretical modelling of the initial state of the QGP, its subsequent hydrodynamic evolution as well as partonic energy loss in the dense and hot medium.