UvA-DARE (Digital Academic Repository) Correlated long-range mixed-harmonic fluctuations measured in pp, p+Pb and low-multiplicity Pb+Pb collisions with the ATLAS detector

Correlations of two ﬂow harmonics v n and v m via three-and four-particle cumulants are measured in 13 TeV pp , 5.02 TeV p +Pb, and 2.76 TeV peripheral Pb+Pb collisions with the ATLAS detector at the LHC. The goal is to understand the multi-particle nature of the long-range collective phenomenon in these collision systems. The large non-ﬂow background from dijet production present in the standard cumulant method is suppressed using a method of subevent cumulants involving two, three and four subevents separated in pseudorapidity. The results show a negative correlation between v 2 and v 3 and a positive correlation between v 2 and v 4 for all collision systems and over the full multiplicity range. However, the magnitudes of the correlations are found to depend on the event multiplicity, the choice of transverse momentum range and collision system. The relative correlation strength, obtained by normalisation of the cumulants with the (cid:2) v 2 n (cid:3) from a two-particle correlation analysis, is similar in the three collision systems and depends weakly on the event multiplicity and transverse momentum. These results based on the subevent methods provide strong evidence of a similar long-range multi-particle collectivity in pp , p +Pb and peripheral Pb+Pb collisions. © 2019 The Author. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP 3 .


Introduction
One of the goals in the studies of azimuthal correlations in high-energy nuclear collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) is to understand the multi-parton dynamics of QCD in the strongly coupled nonperturbative regime [1]. Measurements of azimuthal correlations in small collision systems, such as pp, p+A or d+A collisions, have revealed the ridge phenomenon [2][3][4][5][6]: enhanced production of particle pairs at small azimuthal angle separation, φ, extended over a wide range of pseudorapidity separation, η.
The azimuthal structure has been related to harmonic modulation of particle densities, characterised by a Fourier expansion, where v n and n represent the magnitude and the event-plane angle of the nth-order flow harmonic. They are also conveniently represented by the flow vector: V n = v n e in n . The v n are known to depend on the collision system, but have weak dependence on collision energies [6,7]. The ridge reflects multi-parton dynamics early in the collision and has generated significant interest in the high-energy physics community. A key question is whether the long-range multi-particle collectivity reflects initial momentum correlation from gluon sat-E-mail address: atlas .publications @cern .ch.
uration effects [8], or a final-state hydrodynamic response to the initial transverse collision geometry [9]. Further insight into the ridge phenomenon is obtained via a multi-particle correlation technique, known as cumulants, involving three or more particles [10][11][12]. The multi-particle cumulants probe the event-by-event fluctuation of a single flow harmonic v n , as well as the correlated fluctuations between two flow harmonics, v n and v m . These event-by-event fluctuations are often represented by probability density distributions p(v n ) and p(v n , v m ), respectively. For instance, the four-particle cumulants c n {4} = v 4 n − 2 v 2 n 2 constrain the width of p(v n ) [10], while the four-particle symmetric cumulants sc n,m {4} = v 2 quantify the lowest-order correlation between v n and v m [12]. The three-particle asymmetric cumulants such as ac n {3} = V 2 n v 2n cos 2n( n − 2n ) [5,13] are sensitive to correlations involving both the flow magnitude v n and flow phase n . One of the challenges in the study of azimuthal correlations in small collision systems is how to distinguish the long-range ridge from "non-flow" correlations involving only a few particles, such as resonance decays, jets, or dijets. For two-particle correlations, the non-flow contribution is commonly suppressed by requiring a large using methods similar to those applied in the offline analysis. The HLT enables the high-multiplicity track triggers (HMT) to select events according to the number of tracks having p T > 0.4 GeV matched to the primary vertex, N HLT ch . The different HMT triggers apply additional requirements on either the total transverse energy (E T ) in the calorimeters or the number of hits in the MBTS found by the L1 trigger, as well as on N HLT ch by the HLT trigger. The pp and p+Pb data were collected using combinations of the minimum-bias and HMT triggers. The minimum-bias trigger required either a hit in at least one MBTS counter, or a hit in at least one MBTS counter on each side, or at least one reconstructed track at the HLT seeded by a random trigger at L1. More detailed information about the triggers used for the pp and p+Pb data and their performance can be found in Refs. [7,25] and Refs. [5,26], respectively.

Datasets and Monte Carlo simulations
This analysis is based on ATLAS datasets corresponding to integrated luminosities of 0.9 pb −1 of pp data recorded at  Table 1. The track reconstruction efficiency was determined using simulated Monte Carlo (MC) event samples (Section 4). The pp events were simulated with the Pythia8 MC event generator [27] using the A2 set of tuned parameters with MSTW2008LO parton distribution functions [28]. The HIJING event generator [29] was used to produce Pb+Pb and p+Pb collisions with the same energy and the same boost of the centre-of-mass system as in the data. The detector response was simulated using Geant4 [30,31] with detector conditions matching those during the data-taking. The simulated events and data events are reconstructed with the same algorithms. The MC sample for Pb+Pb events in the multiplicity region of interest is very small, and so the track reconstruction efficiency for Pb+Pb was taken from the larger p+Pb sample reconstructed with the same reconstruction algorithm. The efficiency in p+Pb events was found to be consistent with the efficiency from the Pb+Pb MC simulation [17].

Event and track selection
The offline event selection for the pp and p+Pb data requires at least one reconstructed vertex with its longitudinal position satisfying |z vtx | < 100 mm relative to the nominal interaction point.
The vertex is required to have at least two associated tracks with p T > 0.4 GeV. The mean number of collisions per bunch crossing, μ, was 0.002-0.8 for the 13 TeV pp data, 0.03 for the 2013 p+Pb data, and 0.001-0.006 for the 2016 p+Pb data. In order to suppress additional interactions in the same bunch crossing (referred to as pile-up) in pp collisions, events containing additional vertices with at least four associated tracks are rejected. In p+Pb collisions, events with more than one good vertex, defined as any vertex for which the scalar sum of the p T of the associated tracks is greater than 5 GeV, are rejected. The remaining pile-up events are further suppressed by using the signal in the ZDC in the direction of the Pb beam. This signal is calibrated to the number of detected neutrons, N n , by using the location of the peak corresponding to a single neutron. The distribution of N n in events with pile-up is broader than that for the events without pile-up. Hence a simple requirement on the ZDC signal distribution is used to further suppress events with pile-up, while retaining more than 98% of events without pile-up. The impact of residual pile-up, at the level of 10 −3 , is studied by comparing the results obtained from data with different μ values.
The offline event selection for the Pb+Pb data requires |z vtx | < 100 mm. The selection also requires a time difference | t| < 3 ns between signals in the MBTS trigger counters on either side of the interaction point 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, while maintaining high efficiency for inelastic processes. The fraction of events with more than one interaction after applying these selection criteria is less than 10 −4 .
Charged-particle tracks and collision vertices are reconstructed using algorithms optimised for improved performance for Run-2. In order to compare directly with the pp and p+Pb systems using event selections based on the multiplicity of the collisions, a subset of data from low-multiplicity Pb+Pb collisions, collected during the 2010 LHC heavy-ion run with a minimum-bias trigger, was analysed using the same track reconstruction algorithm as that used for p+Pb collisions. For the Pb+Pb and 2013 p+Pb analyses, tracks are required to have a p T -dependent minimum number of hits in the SCT. The transverse (d 0 ) and longitudinal (z 0 sin θ ) impact parameters of the track relative to the vertex are required to be less than 1.5 mm. Additional requirements |d 0 |/σ d 0 < 3 and |z 0 sin θ|/σ z 0 < 3 are imposed, where σ d 0 and σ z 0 are the uncertainties of the transverse and longitudinal impact parameter values, respectively. A more detailed description of the track selection for the 2010 Pb+Pb data and 2013 p+Pb data can be found in Refs. [5,17].
For all the data taken since the start of Run-2, the track selection criteria make use of the IBL, as described in Refs. [14,25]. For the pp and 2016 p+Pb analyses, the tracks are required to satisfy |d BL 0 | < 1.5 mm and |z 0 sin θ| < 1.5 mm, where d BL 0 is the transverse impact parameter of the track relative to the beam line (BL).
The cumulants are calculated using tracks passing the above selection requirements, and having |η| < 2.5 and 0.3 < p T < 3 GeV or 0.5 < p T < 5 GeV. These two p T ranges are chosen because they were often used in the previous ridge measurements at the LHC [6,7,14,15,17]. However, to count the number of reconstructed charged particles for event-class definition (denoted by N rec ch ), tracks with p T > 0.4 GeV and |η| < 2.5 are used for compatibility with the requirements in the HLT selections described above. Due to different trigger requirements, most of the p+Pb events with N rec ch > 150 are provided by the 2013 dataset, while the 2016 dataset provides most of the events at lower N rec ch . The efficiency of the combined track reconstruction and track selection requirements is estimated using MC samples reconstructed with the same algorithms and selection requirements as in data. Efficiencies, (η, p T ), are evaluated as a function of track η, p T and the number of reconstructed charged-particle tracks, but averaged over the full range in azimuth. The efficiencies are simi-lar for events with the same multiplicity. For all collision systems, the efficiency increases by about 4% as track p T increases from 0.3 GeV to 0.6 GeV. Above 0.6 GeV, the efficiency is independent of p T and reaches 86% (72%) for Run-1 pp and p+Pb, and 83% (70%) for Pb+Pb and Run-2 p+Pb collisions, at η ≈ 0 (|η| > 2). The efficiency is independent of the event multiplicity for N rec ch > 40.
For lower-multiplicity events the efficiency is smaller by up to 3% due to broader d 0 and z 0 sin θ distributions [17].
The fraction of falsely reconstructed charged-particle tracks is also estimated and found to be negligibly small in all datasets. This fraction decreases with increasing track p T , and even at the lowest transverse momenta of 0.3 GeV it is below 1% of the total number of tracks. Therefore, there is no correction for the presence of such tracks in the analysis.
In the simulated events, the reconstruction efficiency reduces the measured charged-particle multiplicity relative to the generated multiplicity for primary charged particles. A correction factor b is used to correct N rec ch to obtain the efficiency-corrected average number of charged particles per event, N ch = b N rec ch . The value of the correction factor is obtained from the MC samples described above, and is found to be nearly independent of N rec ch in the range used in this analysis, N rec ch < 400. Its value and the associated uncertainties are b = 1.29 ± 0.05 for the Pb+Pb and 2013 p+Pb collisions and b = 1.18± 0.05 for Run-2 p+Pb and pp collisions [32]. Both sc n,m {4} and ac 2 {3} are then studied as a function of N ch .

Cumulant method
The multi-particle cumulant method [10] has the advantage of directly reducing non-flow correlations from jets and dijets. The mathematical framework for the standard cumulant is based on the Q-cumulants discussed in Refs. [11,12,33]. It was extended recently to the case of subevent cumulants in Refs. [13,16]. These methods are briefly summarised below.

Cumulants in the standard method
The standard cumulant method calculates k-particle azimuthal correlations, {k} , in one event using a complex number notation [11,12]: where " " denotes a single-event average over all pairs, triplets or quadruplets, respectively. The averages from Eq. (1) can be expressed in terms of per-particle normalised flow vectors q n;l with l = 1, 2... in each event [11]: where the sum runs over all tracks in the event and w j is a weight assigned to the jth track. This weight is constructed to correct for both detector non-uniformity and tracking inefficiency as explained in Section 6. The multi-particle asymmetric and symmetric cumulants are obtained from {k} as: where the averages are taken over the events. This analysis measures three types of cumulants defined in Eq. (3): sc 2,3 {4}, sc 2,4 {4} and ac 2 {3}.

Cumulants in the subevent method
In the standard cumulant method described above, all k-particle multiplets involved in {k} n and {k} n,m are selected using tracks in the entire ID acceptance of |η| < η max = 2.5. To suppress further the non-flow correlations that typically involve a few particles within a localised region in η, the tracks are divided into several subevents, each covering a unique η interval. The multi-particle correlations are then constructed by only correlating tracks between different subevents.
In the two-subevent cumulant method, the tracks are divided into two subevents, labelled by a and b, according to −η max < η a < 0 and 0 ≤ η b < η max . The per-event k-particle azimuthal correlations are evaluated as: where the superscript or subscript a (b) indicates tracks chosen from the subevent a (b). Here the three-and four-particle cumulants are defined as: The two-subevent method suppresses correlations within a single jet (intra-jet correlations), since particles from one jet usually fall in one subevent.
In the three-subevent cumulant method, tracks in each event are divided into three subevents a, b and c, each covering one third of the η range, −η max < η a < −η max /3, |η b | ≤ η max /3 and η max /3 < η c < η max . The multi-particle azimuthal correlations and cumulants are then evaluated as: Since a dijet event usually produces particles in at most two subevents, the three-subevent method efficiently suppresses the non-flow contribution from inter-jet correlations associated with dijets. To maximise the statistical precision, the η range for subevent a is swapped with that for subevent b or c, and the results are averaged to obtain the final values.
The four-subevent cumulant method is only relevant for the symmetric cumulants sc n,m {4}. Tracks in each event are divided into four subevents a, b, c, and d, each covering one quarter of the η range: −η max < η a < −η max /2, −η max /2 ≤ η b < 0, 0 ≤ η c < η max /2, and η max /2 ≤ η d < η max . The multi-particle azimuthal correlations and cumulants are then evaluated as: sc a,b|c,d The four-subevent method based on Eqs. (9) and (10) should further suppress the residual non-flow contributions, for instance when each of the two jets from the dijet falls across the boundary between two neighbouring subevents. To maximise the statistical precision, the η ranges for the four subevents are swapped with each other, and the results are averaged to obtain the final values.

Normalised cumulants
Although the cumulants reflect the nature of the correlation between v n and v m , their magnitudes also depend on the square of single flow harmonics v 2 n and v 2 m , see Eq. (4). The dependence on the single flow harmonics can be scaled out via the normalised cumulants [34,35]: where the v n {2} 2 = v 2 n are flow harmonics obtained using a twoparticle correlation method based on a peripheral subtraction technique [7,14], and c 2 {4} = v 4 2 − 2 v 2 2 2 are four-particle cumulant results from Refs. [17,18]. This definition for nac 2 {3} is motivated by Ref. [36].

Analysis procedure
The measurement of the sc n,m {4} and ac 2 {3} follows the same analysis procedure as for the four-particle cumulants c n {4} in Ref. [18]. The multi-particle cumulants are calculated in three steps using charged particles with |η| < 2.5. In the first step, {2} n , {3} n and {4} n,m from Eqs. (1), (6), (7) and (9) are calculated for each event from particles in one of two different p T ranges, 0.3 < p T < 3 GeV and 0.5 < p T < 5 GeV. The numbers of reconstructed charged particles in these p T ranges are denoted by N sel1 ch and N sel2 ch , respectively.
In the second step, the correlators {k} for 0.3 < p T < 3 GeV (0.5 < p T < 5 GeV) are averaged over events with the same N sel1 The mapping procedure is necessary so that sc 2,3 {4}, sc 2,4 {4} and ac 2 {3} obtained for the two different p T ranges can be compared using a common x-axis defined by N rec ch . The N rec ch value is then converted to N ch , the efficiency-corrected average number of charged particles with p T > 0.4 GeV, as discussed in Section 4.
In order to account for detector inefficiencies and non-uniformity, particle weights used in Eq. (2) are defined as: The additional weight factor d(φ, η) accounts for non-uniformities in the azimuthal acceptance of the detector as a function of η.
All reconstructed charged particles with p T > 0.3 GeV are entered into a two-dimensional histogram N(φ, η), and the weight factor is then obtained as is the track density averaged over φ in the given η bin. This procedure removes most of the φ-dependent non-uniformity in the detector acceptance [17].
In order to calculate the normalised cumulants from Eqs.
where superscripts "peri" and "tmp" indicate quantities for the N rec ch < 20 event class and quantities after the template fit for the event class of interest, respectively. The scale factor F and pedestal

Systematic uncertainties
The evaluation of the systematic uncertainties follows closely the procedure established for the four-particle cumulants c n {4} and described in Ref. [18]. The main sources of systematic uncertainties are related to the detector azimuthal non-uniformity, track selection, track reconstruction efficiency, trigger efficiency and pile-up. Due to the relatively poor statistics and larger non-flow effects, the systematic uncertainties are typically larger in pp collisions. The systematic uncertainties are also generally larger, in percentage, for four-particle cumulants sc n,m {4} than for the three-particle cumulants ac 2 {3}, since the |sc n,m {4}| values are much smaller than those for ac 2 {3}. The systematic uncertainties are generally similar among the two-and three-and four-subevent methods, but are different from those for the standard method, which is strongly influenced by non-flow correlations. The following discussion focuses on the three-subevent method, which is the default method used to present the final results. The effect of detector azimuthal non-uniformity is accounted for using the weight factor d(φ, η). The impact of the weighting procedure is studied by fixing the weight to unity and repeating the analysis. The results are mostly consistent with the nominal results. The corresponding uncertainties for sc n,m {4} vary in the range of 0-4%, 0-2% and 1-2% in pp, p+Pb and Pb+Pb collisions, respectively. The uncertainties for ac 2 {3} vary in the range of 0-2% in pp collisions, and 0-1% in p+Pb and Pb+Pb collisions, respectively.
The systematic uncertainty associated with the track selection is estimated by tightening the |d 0 | and |z 0 sin θ| requirements.
They are each varied from the default requirement of less than 1.5 mm to less than 1 mm. In p+Pb and Pb+Pb collisions, the requirement on the significance of impact parameters, |d 0 |/σ d 0 and |z 0 sin θ|/σ z 0 are also varied from less than 3 to less than 2. For each variation, the tracking efficiency is re-evaluated and the analysis is repeated. For ac 2 {3}, which has a large flow signal, the differences from the nominal results are observed to be less than 2% for all collision systems. For sc n,m {4}, for which the signal is small, the differences from the nominal results are found to be in the range of 2-10% in pp collisions, 2-7% in p+Pb collisions and 2-4% in Pb+Pb collisions. The differences are smaller for results obtained for 0.5 < p T < 5 GeV than those obtained for 0.3 < p T < 3 GeV.
Previous measurements indicate that the azimuthal correlations (both the flow and non-flow components) have a strong dependence on p T , but a relatively weak dependence on η [5,7]. Therefore, p T -dependent systematic effects in the track reconstruction efficiency could affect the cumulant values. The uncertainty in the track reconstruction efficiency is mainly due to differences in the detector conditions and material description between the simulation and the data. The efficiency uncertainty varies between 1% and 4%, depending on track η and p T [7,17]. Its impact on multiparticle cumulants is evaluated by repeating the analysis with the tracking efficiency varied up and down by its corresponding uncertainty as a function of track p T . For the standard cumulant method, which is more sensitive to jets and dijets, the evaluated uncertainty amounts to 2-6% in pp collisions and less than 2% in p+Pb collisions for N ch > 100. For the subevent methods, the evaluated uncertainty is typically less than 3% for most of the N ch ranges.
Most events in pp and p+Pb collisions are collected with the HMT triggers with several online N rec ch thresholds. In order to estimate the possible bias due to trigger inefficiency as a function of N ch , the offline N rec ch requirements are changed such that the HMT trigger efficiency is at least 50% or 80%. The results are obtained independently for each variation. These results are found to be consistent with each other for the subevent methods, and show some differences for the standard cumulant method in the low N ch region. The nominal analysis is performed using the 50% efficiency selection and the differences between the nominal results and those from the 80% efficiency selection are included in the systematic uncertainty. The changes for pp collisions are in the range of 5-15% for sc 2,3 {4}, 2-8% for sc 2,4 {4} and 1-5% for ac 2 {3}.
The ranges for p+Pb collisions are much smaller due to the much sharper turn-on of the trigger efficiency and larger signal: they are estimated to be 1-3% for sc 2,3 {4}, 2-4% for sc 2,4 {4} and 1-2% for ac 2 {3}.
In this analysis, a pile-up rejection criterion is applied to reject events containing additional vertices in pp and p+Pb collisions. In order to check the impact of residual pile-up, the analysis is repeated without the pile-up rejection criterion. No differences are observed in p+Pb collisions, as is expected since the μ values in p+Pb are modest. For the 13 TeV pp dataset, the differences with and without pile-up rejection are in the range of 0-7% for sc 2,3 {4}, 2-15% for sc 2,4 {4} and 2-3% for ac 2 {3}. As a cross-check, the pp data are divided into two samples with approximately equal number of events based on the μ value: μ > 0.4 and μ < 0.4, and the results are compared. No systematic differences are observed between the two independent datasets.
The systematic uncertainties from different sources are added in quadrature to determine the total systematic uncertainty. In p+Pb and Pb+Pb collisions, the total uncertainties are in the range of 3-8% for sc 2,3 {4}, 1-5% for sc 2,4 {4} and 1-4% for ac 2 {3}. In pp collisions, the total uncertainties are larger, mainly due to larger non-flow contribution, larger pile-up and the less sharp turn-on of the HMT triggers. They are in the ranges of 10-20% for sc 2,3 {4}, 10-20% for sc 2,4 {4} and 2-5% for ac 2 {3}. The total systematic uncertainties are generally smaller than the statistical uncertainties.
The v n {2} values used to obtain normalised cumulants from Eqs. (11)-(13) are measured following the prescription of the previous ATLAS publications [7,14], resulting in very similar systematic uncertainties. The correction for the bias of the template fit procedure, as described in Section 6, reduces the sensitivity to the choice of the peripheral N rec ch bin. The uncertainties of normalised cumulants are obtained by propagation of the uncertainties from the original cumulants and v n {2}, taking into account that the correlated systematic uncertainties partially cancel out.

Results
The results are presented in two parts. Section 8.1 presents a detailed comparison between the standard method and subevent methods to demonstrate the ability of the subevent methods to suppress non-flow correlations. Section 8.2 compares the cumulants among pp, p+Pb and Pb+Pb collisions to provide insight into the common nature of collectivity in these systems.

Comparison between standard and subevent methods
The top row of Fig. 1 compares the sc 2,3 {4} values obtained from the standard, two-, three-and four-subevent methods from pp collisions in 0.3 < p T < 3 GeV (left panel) and 0.5 < p T < 5 GeV (right panel). The values from the standard method are positive over the full N ch range, and are larger at lower N ch or in the higher p T range. This behaviour suggests that the sc 2,3 {4} values from the standard method in pp collisions, including those from Ref. [19], are strongly influenced by non-flow effects in all N ch and p T ranges [16]. In contrast, the values from the subevent methods are negative over the full N ch range, and they are slightly more negative at lowest N ch and also more negative at higher p T . The results are consistent among the various subevent methods for 0.3 < p T < 3 GeV. For the high p T region of 0.5 < p T < 5 GeV, results from the two-subevent method are systematically lower than those from the three-and four-subevent methods, suggesting that the two-subevent method may be affected by negative non-flow contributions. Such negative non-flow correlation has been observed in a Pythia8 calculation [16].
The middle row of Fig. 1 shows sc 2 The results are consistent among all four methods across most of the N ch range. In the low N ch region, where the non-flow contribution is expected to be significant, the uncertainties of the results are too large to distinguish between different methods.
The results for the symmetric cumulant sc 2,4 {4} are presented in Fig. 2. The top row shows the sc 2,4 {4} obtained from the standard, two-subevent, three-subevent and four-subevent methods from pp collisions in 0.3 < p T < 3 GeV (left panel) and 0.5 < p T < 5 GeV (right panel). The values of sc 2,4 {4} are positive for all four methods. However, the results from the standard method are much larger than those from the subevent methods and also exhibit a much stronger increase towards the lower N ch region. This behaviour is consistent with the expectation that the standard method is more affected by dijets. Significant differences are also observed between the two-subevent and three-or four-subevent methods at low N ch , but these differences decrease and disappear for N ch > 100. Within the statistical uncertainties of the measurement, no differences are observed between the three-and four-subevent methods. This comparison suggests that the twosubevent method may not be sufficient to reject non-flow correlations from dijets in pp collisions, and methods with three or more subevents are required to suppress the non-flow contribution over the measured N ch range.
Significant differences are observed between the standard method and the subevent methods over the full N ch range. However, no differences are observed among the various subevent methods. These results suggest that the standard method is contaminated by large contributions from non-flow correlations at low N ch , and these contributions may not vanish even at large N ch values. All subevent methods suggest an increase of sc 2,4 {4} toward lower N ch for N ch < 40, which may reflect some residual nonflow correlations in this region.
The bottom row of Fig. 2 shows sc 2,4 {4} from Pb+Pb collisions. The sc 2,4 {4} values increase gradually with N ch for all four methods. This increase reflects the known fact that the v 2 increases with N ch in Pb+Pb collisions [37]. The values from the standard method are systematically larger than those from the subevent methods, and this difference varies slowly with N ch , similar to the behaviour observed in p+Pb collisions in the high N ch region.
The results for the asymmetric cumulant ac 2 {3} are presented in Fig. 3. The top row shows the results obtained from the standard, two-subevent, and three-subevent methods from pp collisions in 0.3 < p T < 3 GeV (left panel) and 0.5 < p T < 5 GeV (right panel). The results are positive for all methods. The results from the standard method are much larger than those from the subevent methods, consistent with the expectation that the standard method is more affected by non-flow correlations from dijets. Significant differences are also observed between the two-subevent and three-subevent methods at low N ch , but these differences decrease and disappear at N ch > 100. The ac 2 {3} values from the three-subevent method show a slight increase for N ch < 40 but are nearly constant for N ch > 40. This behaviour suggests that in the three-subevent method, the non-flow contribution may play some role at N ch < 40, but is negligible for N ch > 40. Therefore, the ac 2 {3} from the three-subevent method supports the existence of a three-particle long-range collective flow that is nearly independent of N ch in pp collisions, consistent with the N ch -independent behaviour of v 2 and v 4 observed previously in the two-particle correlation analysis   In the subevent methods, the influence of non-flow contributions is very small for N ch > 60 in both collision systems, and therefore the N ch dependence of ac 2 {3} reflects the N ch dependence of the v 2 and v 4 . The ac 2 {3} values from the subevent methods increase with N ch , and the increase is stronger in Pb+Pb collisions. This is consistent with previous observations that v 2 and v 4 increase with N ch more strongly in Pb+Pb than in p+Pb collisions [17].
The values of sc 2,4 {4} and ac 2 {3}, which are both measures of correlations between v 2 and v 4 , show significant differences between the standard method and the subevent methods, as shown in Figs. 2 and 3. The N ch dependence of these differences decreases gradually with N ch , and is consistent with an influence of non-flow that is expected to scale as 1/ N ch . However, these differences seem to persist for N ch > 200 in p+Pb collisions and for N ch > 150 in Pb+Pb collisions, which is not compatible with the predicted behaviour of non-flow correlations. The differences at with v 2 , and therefore they are expected to reduce the sc 2,4 {4} and ac 2 {3} in the subevent method. Therefore, the observed differences between the standard method and subevent method reflect the combined contribution from non-flow correlations, which dominates in the low N ch region, and decorrelation, which is more important at large N ch (see further discussion in the Appendix B).
The results presented above suggest that the three-subevent method is sufficient to suppress most of the non-flow effects. It is therefore used as the default method for the discussion below. These results support the existence of a negative correlation between v 2 and v 3 and a positive correlation between v 2 and v 4 . In the multiplicity range covered by the pp collisions, N ch < 150, the results for symmetric cumulants sc 2,3 {4} and sc 2,4 {4} are similar among the three systems. In the range N ch > 150, |sc 2,3 {4}| and sc 2,4 {4} are larger in Pb+Pb than in p+Pb collisions. The results for ac 2 {3} are similar among the three systems at N ch < 100, but they deviate from each other at higher N ch . The pp data are approximately constant or decrease slightly with N ch , while the p+Pb and Pb+Pb data show significant increases as a function of N ch . The bottom row shows the results for the higher p T range of 0.5 < p T < 5 GeV, where similar trends are observed.  This suggests that the v 2 3 values from the template fit method [7] may be significantly underestimated. As pointed out in Ref.

Comparison between collision systems
[7] and emphasised in Appendix A, the template fit method, and other methods based on peripheral subtraction in general [5,15], tend to underestimate the odd flow harmonics, due to the presence of a large away-side peak at φ ∼ π in the two-particle correlation function. The comparison of sc 2,3 {4} and nsc 2,3 {4} among different collision systems provides indirect evidence of this underestimation of v 2 3 . Fig. 5 shows that the normalised cumulants are consistent between 0.3 < p T < 3 GeV and 0.5 < p T < 5 GeV. On the other hand, the magnitudes of the cumulants in Fig. 4 differ by a large factor between the two p T ranges: about a factor of three for sc 2,3 {4} and sc 2,4 {4}, and a factor of two for ac 2 {3}. These results suggest that the p T dependence of sc 2,3 {4}, sc 2,4 {4} and ac 2 {3} largely reflects the p T dependence of the v n at the single-particle level.

Discussion
Three-and four-particle cumulants involving correlations between two harmonics of different order v n and v m are measured in √ s = 13 TeV pp, √ s NN = 5.02 TeV p+Pb, and low-multiplicity √ s NN = 2.76 TeV Pb+Pb collisions with the ATLAS detector at the LHC, with total integrated luminosities of 0.9 pb −1 , 28 nb −1 , and 7 μb −1 , respectively. The correlation between v n and v m is studied using four-particle symmetric cumulants, sc 2,3 {4} and Significant differences are observed between the standard method and the subevent methods over the full N ch range in pp collisions, as well as over the low N ch range in p+Pb and Pb+Pb collisions. The differences are larger for particles at higher p T or at smaller N ch . When analysed with the standard method in pp collisions, this behaviour is compatible with the dominance of the non-flow correlations rather than the long-range collective flow correlations. Systematic, but much smaller, differences are also observed in the low N ch region between the two-subevent method and three-or four-subevent methods, which indicate that the two-subevent method may still be affected by correlations arising from jets. On the other hand no differences are observed between the three-subevent and four-subevent methods, within experimental uncertainties, suggesting that methods with three or more subevents are sufficient to reject non-flow correlations from jets. Therefore, the three-subevent method is used to present the main results in this analysis.
The three-subevent method provides a measurement of negative sc 2,3 {4} and positive sc 2,4 {4} and ac 2 {3} over nearly the full N ch range and in all three collision systems. These results indicate a negative correlation between v 2 and v 3 and a positive are also similar to each other for 0.5 < p T < 5 GeV as well as 0.3 < p T < 3 GeV. This suggests that the N ch , p T and system dependence of the sc 2,3 {4}, sc 2,4 {4} and ac 2 {3} reflect mostly the N ch , p T and system dependence of v 2 n , but the relative strengths of the correlations are similar for the three collision systems.
The new results obtained with the subevent cumulant technique provide further evidence that the ridge is indeed a longrange collective phenomenon involving many particles distributed across a broad rapidity interval. The similarity between different collision systems for nsc 2,3 {4}, nsc 2,4 {4} and nac 2 {3}, and the weak dependence of these observables on the p T range and N ch , largely free from non-flow effects, provide an important input towards understanding the space-time dynamics and the properties of the medium created in small collision systems. These results provide inputs to distinguish between models based on initial-state momentum correlations and models based on final-state hydrodynamics.  [45].

Appendix A. Improvement to the template fit procedure
In order to separate the long-range ridge from other nonflow sources, especially dijets, the ATLAS Collaboration developed a template fitting procedure described in Refs. [7,14]. The first step is to construct a φ distribution of particle pairs with large pseudorapidity separation | η| > 2, the so-called "per-trigger" particle yield, Y ( φ), for a given N rec ch range. The | η| > 2 requirement suppresses the intra-jet and other short-range correlations, and in small collision systems the resulting Y ( φ) distributions are known to be dominated by away-side jet correlations [4,5,14]. This away-side non-flow component is peaked at φ ∼ π , and leads to a significant bias in the flow coefficients v n , especially for the odd harmonics.
To subtract the away-side jet correlations, the measured Y ( φ) distribution in a given N rec ch interval is assumed to be a sum of a scaled "peripheral" distribution Y ( φ) peri , obtained for lowmultiplicity events N rec ch < 20, and a constant pedestal modulated by cos(n φ) for n ≥ 2 [7,14]: The scale factor F and pedestal G tmp are fixed by the fit, and With the assumption that the shape of the dijet component is independent of N rec ch , and the magnitudes of the dijet components are related by the scale factor F : peri jet , Eq. (14) can be written as: Comparing with Eqs. (15) and (16), one obtains G cent = G tmp + F G peri and the following relation: One important feature of the template fit analysis is the assumption that the dijet component Y ( φ) jet is independent of N ch . In Ref.
[7], the uncertainty associated with this assumption is studied by changing the default peripheral interval from N rec ch < 20 to N rec ch < 10 and 10 ≤ N rec ch < 20. It was found that the v n {2, tmp} values are relatively insensitive to the choice of peripheral interval for n = 2 and n = 4, but the sensitivity is much larger for n = 3. This finding is reproduced in Fig. 6 for pp collisions, which shows that the v 3 {2, tmp} 2 values obtained via Eq. (14) differ substantially for the different N rec ch ranges. In addition to the template fit with and without the above mentioned correction procedure, the ATLAS and CMS collaborations also calculated directly the v n {2} values via a Fourier transform of the Y ( φ) distribution without dijet subtraction [7,19]. The differences between the direct Fourier transform and template fit reflect mainly the away-side jet contribution subtracted by the template fit procedure, and therefore give a sense of the magnitude of unknown systematic uncertainties associated with the template fit procedure. If these differences are too large, the v n {2, tmp} values may be sensitive to the systematic effects associated with the assumption that the shape of Y ( φ) jet is independent of N rec ch .  Fig. 7 shows that the changes introduced by the correction procedure described above are small in all cases and for all harmonics. The values of the even-order harmonics, v 2 and v 4 , are also quite similar to those obtained from the direct Fourier transformation, reflecting the fact that the dijet correlations have very little influence on the even-order harmonics. On the other hand, significant differences are observed between the direct Fourier transform and template fit for v 3 , especially in the pp collisions, due to the influence of Y ( φ) jet , a trend observed and discussed previously in Refs. [7,15]. The template fit procedure is able to subtract the dijet correlations and change the sign of v 3 , but also introduces a large uncertainty associated with the procedure. As discussed in Section 8.2, the behaviour of the symmetric cumulants sc 2,3 {4} in Fig. 4 and normalised cumulants nsc 2,3 {4} in Fig. 5 in pp collisions, suggest that the v 3 values from the template fit procedure are significantly underestimated due to the presence  of a large residual non-flow bias. In contrast, the differences of v 3 between the direct Fourier transform and the template fit are much smaller in the p+Pb and the Pb+Pb collisions, except in the very low N ch region. Therefore, the v 3 values in p+Pb and Pb+Pb systems extracted from the template fit procedure are expected to be less affected by the dijets.

Appendix B. Effects of flow decorrelations in the subevent cumulant methods
As discussed in Section 8, the differences between the standard method and subevent methods for ac 2 {3} can be partially at-  {3} at larger N ch region. This is because the subevent for V 4 is in between the two subevents used to calculate V 2 . This configuration has much smaller decorrelation effects [41]. For ac a,b|c or b,c|a 2 {3}, the two subevents for V 2 are on the same side of the subevent for V 4 , leading to larger decorrelation effects. Interestingly, such configuration gives results that are very similar to those from the two-subevent method. Figs. 9 and 10 show the results for the p+Pb and pp collisions, respectively. Similar observations as in Pb+Pb collisions can be made, although in pp collisions the results from the two-subevent method are larger than those obtained with the three-subevent method, due to significant non-flow contribution even in the large N ch region.  [40] CMS Collaboration, Evidence for transverse momentum and pseudorapidity dependent event plane fluctuations in PbPb and pPb collisions, Phys.