Correlated long-range mixed-harmonic fluctuations measured in $pp$, $p$+Pb and low-multiplicity Pb+Pb collisions with the ATLAS detector

Correlations of two flow 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-flow 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 strongly 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 $\langle v_n^2\rangle$ 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.


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 non-perturbative 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, dN dφ ∝ 1 + 2 ∑ ∞ n=1 v n cos n(φ − Φ n ), where v n and Φ n represent the magnitude and the event-plane angle of the n th -order flow harmonic.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 saturation 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 eventby-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  n v 2 m ⟩ − ⟨v 2 n ⟩ ⟨v 2 m ⟩ 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 ∆η gap between the two particles in each pair and a peripheral subtraction procedure [3-5, 7, 14, 15].For multi-particle cumulants, the non-flow contributions can be suppressed by requiring correlation between particles from different subevents separated in η, while preserving the genuine long-range multi-particle correlations associated with the ridge.Here each subevent is a collection of particles in a given η range.This so-called "subevent method" has been demonstrated to measure reliably c n {4} and sc n,m {4} [13,16].In contrast, c n {4} and sc n,m {4} based on the standard cumulant method are contaminated by non-flow correlations over the full multiplicity range in pp collisions and the low multiplicity region in p+A collisions [16].In small collision systems, measurements have been performed for c n {4} with both the standard [15,17] and subevent methods [18], and for sc n,m {4} with the standard method [19].The subevent method has not yet been used to measure sc n,m {4}, and no measurements of ac n {3} have ever been attempted in small collision systems.This Letter presents measurements of sc 2,3 {4}, sc 2,4 {4} and ac 2 {3} in pp collisions at √ s = 13 TeV, p+Pb collisions at √ s NN = 5.02 TeV and low-multiplicity Pb+Pb collisions at √ s NN = 2.76 TeV.They are obtained using two-, three-and four-subevent cumulant methods and are compared with results from the standard cumulant method.The cumulants are normalised by the ⟨v 2 n ⟩ obtained from a two-particle correlation analysis [7] to quantify their relative correlation strength.The measurements suggest that the results obtained with the standard method are strongly contaminated by correlations from non-flow data at √ s NN = 2.76 TeV.The 2.76 TeV Pb+Pb data were collected in 2010.The p+Pb data were mainly collected in 2013, but also include 0.3 nb −1 of data collected in 2016, which increase the number of events at moderate multiplicity (see Section 4).During both p+Pb runs, the LHC was configured to provide a 4 TeV proton beam and a 1.57TeV per-nucleon Pb beam, which produced collisions at √ s NN = 5.02 TeV, with a rapidity shift of 0.465 of the nucleon-nucleon centre-of-mass frame towards the proton beam direction relative to the ATLAS rest frame.The direction of the Pb beam is always defined to have negative pseudorapidity.The 13 TeV pp data were collected during several special runs of the LHC with low pile-up in 2015 and 2016.A summary of the datasets used in this analysis is shown in Table 1.
Table 1: The list of datasets used in this analysis.

Pb+Pb p+Pb pp
Integrated luminosity (year) 7 µb −1 (2010) 28 nb −1 (2013) 0.07 pb −1 (2015) 0.3 nb −1 (2016) 0.84 pb −1 (2016) The track reconstruction efficiency was determined using simulated Monte Carlo (MC) event samples (Section 4).The pp events were simulated with the P 8 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 G 4 [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.8for the 13 TeV pp data, 0.03 for the 2013 p+Pb data, and 0.001-0.006for 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 similar 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 j th 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 "⟪⟫" represents a weighted average of ⟨{k}⟩ over an event ensemble with similar N rec ch .One averages first over all distinct pairs, triplets or quadruplets in one event to obtain ⟨{2} n ⟩, ⟨{2} m ⟩, ⟨{3} n ⟩ and ⟨{4} n,m ⟩.Then the obtained values are averaged over an event ensemble with similar N rec ch to obtain sc n,m {4} and ac n {3}.In the absence of non-flow correlations, sc n,m {4} and ac n {3} measure the correlation between v n and v m or between v n and v 2n : 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 threeand 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: and 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: , and η max 2 ≤ η d < η max .The multi-particle azimuthal correlations and cumulants are then evaluated as: The four-subevent method based on Eqs. ( 7) and ( 8) 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 two-particle correlation method based on a peripheral subtraction technique [7,14].In this method, the non-flow effects due to dijets are estimated using low-multiplicity events and then subtracted via a "template fit" procedure [7,14].

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), ( 5), ( 6) and ( 7) 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 ch (N sel2 ch ) to obtain ⟪{k}⟫, and then sc 2,3 {4}, sc 2,4 {4} and ac 2 {3}.The sc 2,3 {4}, sc 2,4 {4} and ac 2 {3} values are then averaged in broader multiplicity ranges of the event ensemble, weighted by number of events, to obtain statistically significant results.
In the third step, the sc 2,3 {4}, sc 2,4 {4} and ac 2 {3} values obtained for a given N sel1 ch or N sel2 ch are mapped to ⟨N rec ch ⟩, the average number of reconstructed charged particles with p T > 0.4 GeV.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 d(φ, η) ≡ ⟨N(η)⟩ N(φ, η), where ⟨N(η)⟩ 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. ( 9)-( 11), the flow harmonics v n {2} are obtained from a "template fit" of two-particle ∆φ correlation as described in Refs.[7,14].The v n {2} values are calculated identically to the procedure used in the previous ATLAS publications [7,14], but are further corrected for a bias, which exists only if v n {2} changes with N rec ch .The details of the correction procedure are given in the Appendix A and are discussed briefly below.
The standard procedure of Refs.[7,14] first constructs a ∆φ distribution for pairs of tracks with ∆η > 2: the per-trigger-particle yield Y (∆φ) for a given N rec ch range.The dominating non-flow jet peak at ∆φ ∼ π is estimated using low-multiplicity events with N rec ch < 20 and separated via a template fit procedure, and the harmonic modulation of the remaining component is taken as the v n {2} 2 [7]: 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 G tmp are fixed by the fit, and v n {2, tmp} are calculated from a Fourier transform.This procedure implicitly assumes that v n {2} is independent of N rec ch , and requires a small correction if v n {2} does change with N rec ch (Appendix A).In p+Pb and Pb+Pb collisions, this correction in the N rec ch > 100 region amounts to a 2-6% reduction for v 2 {2, tmp} and a 4-9% reduction for v 3 {2, tmp} and v 4 {2, tmp}.The correction is smaller for v 2 {2, tmp} in pp collisions as it is nearly independent of N rec ch [7].

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 multi-particle 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. ( 9)- (11) 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 Figure 1 compares the sc 2,3 {4} values obtained from the standard, two-, three-and foursubevent 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 P 8 calculation [16].
The middle row of Figure 1 shows sc 2,3 {4} from p+Pb collisions.At ⟨N ch ⟩ > 140, the values are negative and consistent among all four methods, reflecting genuine long-range collective correlations.At ⟨N ch ⟩ < 140, the values are different between the standard method and the subevent methods.The sc 2,3 {4} from the standard method changes sign around ⟨N ch ⟩ ∼ 80 and remains positive at lower ⟨N ch ⟩, reflecting the contribution from non-flow correlations.In contrast, the sc 2,3 {4} from various subevent methods are negative and consistent with each other at ⟨N ch ⟩ < 140, suggesting that they mainly reflect the genuine long-range correlations.
The bottom row of Figure 1 shows sc 2,3 {4} from Pb+Pb collisions.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 Figure 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 two-subevent 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.
The middle row of Figure 2 shows sc 2,4 {4} from p+Pb collisions.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 non-flow correlations in this region.
The bottom row of Figure 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 [36].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 Figure 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 [7].
The middle and bottom rows of Figure 3 show ac 2 {3} from p+Pb and Pb+Pb collisions, respectively.The ac 2 {3} values from the standard method have a significant non-flow contribution up to ⟨N ch ⟩ ∼ 200 in p+Pb collisions and ⟨N ch ⟩ ∼ 80 in Pb+Pb collisions.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 Figures 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 large ⟨N ch ⟩ may arise from longitudinal flow decorrelations [37,38], which have been measured by CMS [39] and ATLAS [40].Decorrelation effects are found to be large for v 4 and strongly correlated 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   which dominates in the low ⟨N ch ⟩ region, and decorrelation, which is more important at large ⟨N ch ⟩.
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.

Comparison between collision systems
Figure 4 shows a direct comparison of cumulants for the three collision systems.The three panels in the top row show the results for sc 2,3 {4}, sc 2,4 {4} and ac 2 {3}, respectively, for 0.3 < p T < 3 GeV.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 .Such correlation patterns have previously been observed in large collision systems [41][42][43], but are now confirmed also in the small collision systems, once non-flow effects are adequately suppressed.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.Figure 5 shows the results for normalised cumulants, nsc 2,3 {4}, nsc 2,4 {4} and nac 2 {3}, compared among the three systems.The normalised cumulants generally show a much weaker ⟨N ch ⟩ dependence at ⟨N ch ⟩ > 100, where the statistical uncertainties are small.This behaviour implies that the strong ⟨N ch ⟩ dependence of the sc n,m {4} and ac 2 {3} values reflects the ⟨N ch ⟩ dependence of the v n values, and these dependences are removed in the normalised cumulants.The normalised cumulants are also similar among different collision systems at large ⟨N ch ⟩, although some differences at the relative level of 20-30% are observed for smaller ⟨N ch ⟩.The only exception is nsc 2,3 {4}, whose values in the pp collisions are very different from those in p+Pb and Pb+Pb collisions.In contrast, the sc 2,3 {4} values in Figure 4 are close among different systems.This suggests that the ⟨v 2  3 ⟩ values from the template fit method [7] may be significantly underestimated.As pointed out in Ref. [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 ⟩. Figure 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 Figure 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 sc 2,4 {4}, and the three-particle asymmetric cumulant ac 2 {3}.The symmetric cumulants sc n,m {4} = ⟨v 2 n v 2 m ⟩ − ⟨v 2 n ⟩ ⟨v 2 m ⟩ probe the correlation of the flow magnitudes, while the asymmetric cumulant ac 2 {3} = ⟨v 2  2 v 4 cos 4(Φ 2 − Φ 4 )⟩ is sensitive to correlations involving both the flow magnitude v n and flow phase Φ n .They are calculated using the standard cumulant method, as well as the two-, three-and four-subevent methods to suppress non-flow effects.The final results are presented as a function of the average number of charged particles with p T > 0.4 GeV, ⟨N ch ⟩.
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 correlation between v 2 and v 4 .Such correlation patterns have previously been observed in large collision systems [41][42][43], but are now confirmed in small collision systems, once non-flow effects are adequately suppressed.The values of sc 2,3 {4} and sc 2,4 {4} are consistent in pp and p+Pb collisions over the same ⟨N ch ⟩ range, but their magnitudes at large ⟨N ch ⟩ are much smaller than those for Pb+Pb collisions.The values of ac 2 {3} are similar at very low ⟨N ch ⟩ among the three systems, but are very different at large ⟨N ch ⟩.On the other hand, after scaling by the ⟨v 2 n ⟩ estimated from a two-particle analysis [7,14], the resulting normalised cumulants nsc 2,3 {4}, nsc 2,4 {4} and nac 2 {3} show a much weaker dependence on ⟨N ch ⟩, and their values are much closer to each other among the three systems.The magnitudes of the normalised cumulants 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 long-range 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.

Appendix A Improvement to the template fit procedure
In order to separate the long-range ridge from other non-flow 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 intrajet 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 low-multiplicity 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 v n {2, tmp} are calculated from a Fourier transform.On the other hand, both Y (∆φ) and Y (∆φ) peri contain a dijet component and flow component: 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: Y (∆φ) cent jet = FY (∆φ) peri jet , Eq. ( 12) can be written as: Comparing with Eqs. ( 13) and ( 14), one obtains G cent = G tmp + FG peri and the following relation: which shows that v n {2, tmp} from the template fit differs from the true v n {2} by a correction term that vanishes if and only if v n {2} is independent of N rec ch .Since the true flow harmonics in the peripheral interval v n {2, peri} are unknown in principle, the correction is applied starting from the third-lowest N rec ch interval (40 ≤ N rec ch < 60) in this analysis, by using v n {2, tmp} of the second N rec ch interval (20 ≤ N rec ch < 40) as an estimate of the true flow harmonics.Since the non-flow contribution primarily affects the odd harmonics, the v 3 {2, tmp} 2 may become negative in the first few N rec ch intervals in pp collisions.In such cases, the correction starts from the second N rec ch interval with positive v 3 {2, tmp} 2 (60 ≤ N rec ch < 80) by using v 3 {2, tmp} from the previous N rec ch interval (40 ≤ N rec ch < 60).
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 Figure 6 for pp collisions, which shows that the v 3 {2, tmp} 2 values obtained via Eq.( 12) 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 .Figure 7 compares the v n {2} in 0.3 < p T < 3 GeV obtained from Y (∆φ) using three methods: a direct Fourier transform (solid circles), a template fit (open circles) and a template fit corrected for the bias (open squares), as described above.The systematic uncertainties for the template fit results are nearly the same as those from Ref. [7]. Figure 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 Figure 4 and normalised cumulants nsc 2,3 {4} in Figure 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.

Switzerland.
f Also at Departament de Fisica de la Universitat Autonoma de Barcelona, Barcelona; Spain.g Also at Departamento de Física Teorica y del Cosmos, Universidad de Granada, Granada (Spain); p Also at Department of Physics, University of Fribourg, Fribourg; Switzerland.q Also at Department of Physics, University of Michigan, Ann Arbor MI; United States of America.r Also at Dipartimento di Fisica E. Fermi, Università di Pisa, Pisa; Italy.s Also at Giresun University, Faculty of Engineering, Giresun; Turkey.t Also at Graduate School of Science, Osaka University, Osaka; Japan.u Also at Hellenic Open University, Patras; Greece.v Also at Horia Hulubei National Institute of Physics and Nuclear Engineering, Bucharest; Romania.w Also at II.Physikalisches Institut, Georg-August-Universität Göttingen, Göttingen; Germany.

Figure 1 :
Figure 1: The symmetric cumulant sc 2,3 {4} as a function of ⟨N ch ⟩ for 0.3 < p T < 3 GeV (left panels) and 0.5 < p T < 5 GeV (right panels) obtained for pp collisions (top row), p+Pb collisions (middle row) and lowmultiplicity Pb+Pb collisions (bottom row).In each panel, the sc 2,3 {4} is obtained from the standard method (filled symbol), the two-subevent method (open circles), three-subevent method (open squares) and four-subevent method (open diamonds).The error bars and shaded boxes represent the statistical and systematic uncertainties, respectively.

Figure 3 :
Figure 3: The asymmetric cumulant ac 2 {3} as a function of ⟨N ch ⟩ for 0.3 < p T < 3 GeV (left panels) and 0.5 < p T < 5 GeV (right panels) obtained for pp collisions (top row), p+Pb collisions (middle row) and lowmultiplicity Pb+Pb collisions (bottom row).In each panel, the ac 2 {3} is obtained from the standard method (filled symbol), two-subevent method (open circles), and three-subevent method (open squares).The error bars and shaded boxes represent the statistical and systematic uncertainties, respectively.

Figure 6 :
Figure 6: The values of v n {2, tmp} 2 obtained following the template fit procedure given in Eq. (12) [7] in pp collisions for n = 2 (left panel), n = 3 (middle panel) and n = 4 (right panel).In each panel, the values are calculated for three peripheral N rec ch intervals: N rec ch < 20, N rec ch < 10 and 10 ≤ N rec ch < 20.Only statistical uncertainties are shown.

Figure 7 :
Figure 7: The v 2 (left column), v 3 (middle column) and v 4 (right column) obtained from two-particle correlations in 0.3 < p T < 3 GeV in pp (top row), p+Pb (middle row) and Pb+Pb (bottom row) collisions.In each panel, they are compared between three methods: direct Fourier transformation (solid circles), template fit (open circles) and the improved template fit (open squares).The error bars and shaded boxes represent the statistical and systematic uncertainties, respectively.
Also at Department of Applied Physics and Astronomy, University of Sharjah, Sharjah; United Arab Emirates.i Also at Department of Financial and Management Engineering, University of the Aegean, Chios; Also at Department of Physics and Astronomy, University of Louisville, Louisville, KY; United States of America.k Also at Department of Physics and Astronomy, University of Sheffield, Sheffield; United Kingdom.l Also at Department of Physics, California State University, Fresno CA; United States of America.m Also at Department of Physics, California State University, Sacramento CA; United States of America.n Also at Department of Physics, King's College London, London; United Kingdom.o Also at Department of Physics, St. Petersburg State Polytechnical University, St. Petersburg; Russia.
Spain.h j Also at Institut für Experimentalphysik, Universität Hamburg, Hamburg; Germany.z Also at Institute for Mathematics, Astrophysics and Particle Physics, Radboud University Nijmegen/Nikhef, Nijmegen; Netherlands.aa Also at Institute for Particle and Nuclear Physics, Wigner Research Centre for Physics, Budapest; Hungary.ab Also at Institute of Particle Physics (IPP); Canada.ac Also at Institute of Physics, Academia Sinica, Taipei; Taiwan.ad Also at Institute of Physics, Azerbaijan Academy of Sciences, Baku; Azerbaijan.ae Also at Institute of Theoretical Physics, Ilia State University, Tbilisi; Georgia.a f Also at LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, Orsay; France.ag Also at Louisiana Tech University, Ruston LA; United States of America.ah Also at Manhattan College, New York NY; United States of America.ai Also at Moscow Institute of Physics and Technology State University, Dolgoprudny; Russia.a j Also at National Research Nuclear University MEPhI, Moscow; Russia.ak Also at Physikalisches Institut, Albert-Ludwigs-Universität Freiburg, Freiburg; Germany.al Also at School of Physics, Sun Yat-sen University, Guangzhou; China.am Also at The City College of New York, New York NY; United States of America.an Also at The Collaborative Innovation Center of Quantum Matter (CICQM), Beijing; China.ao Also at Tomsk State University, Tomsk, and Moscow Institute of Physics and Technology State University, Dolgoprudny; Russia.ap Also at TRIUMF, Vancouver BC; Canada.aq Also at Universita di Napoli Parthenope, Napoli; Italy.
y * Deceased