Multi-particle azimuthal correlations in p-Pb and Pb-Pb collisions at the CERN Large Hadron Collider

Measurements of multi-particle azimuthal correlations (cumulants) for charged particles in p-Pb and Pb-Pb collisions are presented. They help address the question of whether there is evidence for global, flow-like, azimuthal correlations in the p-Pb system. Comparisons are made to measurements from the larger Pb-Pb system, where such evidence is established. In particular, the second harmonic two-particle cumulants are found to decrease with multiplicity, characteristic of a dominance of few-particle correlations in p-Pb collisions. However, when a $|\Delta \eta|$ gap is placed to suppress such correlations, the two-particle cumulants begin to rise at high-multiplicity, indicating the presence of global azimuthal correlations. The Pb-Pb values are higher than the p-Pb values at similar multiplicities. In both systems, the second harmonic four-particle cumulants exhibit a transition from positive to negative values when the multiplicity increases. The negative values allow for a measurement of $v_{2}\{4\}$ to be made, which is found to be higher in Pb-Pb collisions at similar multiplicities. The second harmonic six-particle cumulants are also found to be higher in Pb-Pb collisions. In Pb-Pb collisions, we generally find $v_{2}\{4\} \simeq v_{2}\{6\}\neq 0$ which is indicative of a Bessel-Gaussian function for the $v_{2}$ distribution. For very high-multiplicity Pb-Pb collisions, we observe that the four- and six-particle cumulants become consistent with 0. Finally, third harmonic two-particle cumulants in p-Pb and Pb-Pb are measured. These are found to be similar for overlapping multiplicities, when a $|\Delta\eta|>1.4$ gap is placed.


Introduction
The primary goal of studies with relativistic heavy-ion collisions is to create the Quark Gluon Plasma (QGP), a unique state of matter where quarks and gluons can move freely over large volumes in comparison to the typical size of a hadron. Studies of azimuthal anisotropy for produced particles have contributed significantly to the characterization of the system created in heavy-ion collisions. These studies are based on a Fourier expansion of the azimuthal distribution given by [1]: where ϕ is the azimuthal angle of produced particles. In heavy-ion collisions, the v n terms generally represent flow coefficients where n is the flow harmonic and Ψ n is the corresponding flow angle. The flow coefficients are believed to reflect the response of the system to spatial anisotropies in the initial state. Measurements of the second harmonic flow coefficient (v 2 , elliptic flow) received keen attention at Relativistic Heavy Ion Collider (RHIC), where the correspondence with hydrodynamic calculations in Au+Au √ s NN = 200 GeV collisions indicated that an almost perfect liquid had been produced in the laboratory [2,3,4,5]. Larger values of integrated v 2 have been observed at the Large Hadron Collider (LHC) in Pb-Pb at √ s NN = 2.76 TeV collisions, indicating that the system created at this new energy regime still behaves as an almost ideal liquid [6]. While the initial state anisotropy is usually dominated by an elliptical overlap area which gives rise to v 2 , measurements of the third harmonic flow (v 3 , triangular flow) demonstrated initial state fluctuations modulate the overlap area, and they provide additional constraints to the transport coefficients of the system (e.g. the value of the shear viscosity over entropy ratio η/s) [7,8,9,10,11]. The combination of the second and higher harmonic flow coefficients manifest themselves in two-particle correlation structures (along ∆η) such as the away-side double hump (∆ϕ ∼ π), and near-side ridge (∆ϕ ∼ 0) observed both at RHIC and the LHC.
The study of p-Pb collisions, which usually provides baseline measurements for the quantification of cold nuclear matter effects, led to a number of unexpected results [12,13,14,15,16,17,18]. The CMS Collaboration reported the development of a near-side ridge-like structure in high-multiplicity p-Pb collisions [12,16]. We discovered a symmetric double ridge structure on both the near-and the away-side after subtracting from the high-multiplicity p-Pb correlation function the dominant jet contribution using the low multiplicity events [13]. The ATLAS Collaboration confirmed the appearance of such structure using a similar subtraction technique [14]. We extended the measurements to identified hadrons and reported a mass ordering in the p T differential v 2 measurements for the different species, with a crossing of p and π v 2 at large p T [17]. Around a similar time, the CMS and ATLAS Collaborations measured finite values of v 2 from four particle correlations. [15,16].
The origin of the ridge structure in p-Pb collisions has been the subject of speculation within the heavyion community [19,20,21,22]. It has been suggested that a high enough energy density is achieved in p-Pb at √ s NN = 5.02 TeV collisions to induce hydrodynamic flow using a Lattice QCD equation of state [19]. Combined with spatial anisotropies in the initial p-Pb state, this mechanism would induce global correlations of soft particles with significant values of v 2 and v 3 . A second proposal is that the ridge arises from collimated (in ∆ϕ) correlated two-gluon production from the Color Glass Condensate (CGC) [20]. This leads to few-particle correlations, rather than a global modulation of soft particles. Finally, the third explanation invokes the CGC initial state with a finite number of sources that form the eccentricity [21]. In contrast to the previous explanation, this approach allows for non-zero values of v 2 from four, six, and eight particle correlations in high multiplicity p-Pb collisions.
Whether the current measurements in high-multiplicity p-Pb events reveal the onset of collective behavior, or can be explained in terms of few-particle correlations (i.e. non-flow), is the main goal of this analysis. We report the multiplicity dependence of the two-, four-, and six-particle correlations (cumulants) for charged particles, that can be used as a tool to investigate multi-particle correlations of various harmonics [23,24]. We present the results in both p-Pb and Pb-Pb collisions at √ s NN = 5.02 TeV and √ s NN = 2.76 TeV respectively. The multiplicity dependence of these measurements will help decipher how flow and non-flow contribute. In Sec. II, we will introduce multi-particle cumulants and discuss their response to non-flow and flow fluctuations. In Sec. III we will describe the analysis details. Sec. IV shows our results, and Sec. V presents our summary.

Multi-particle cumulants
The measurements of v n in Eq. 1 can be done using a variety of methods, which have different sensitivities to flow fluctuations (event-wise variations in the flow coefficients) and non-flow. Non-flow refers to correlations not related to the common symmetry plane Ψ n , such as those due to resonances and jets. Multi-particle cumulants are utilized since their response to flow fluctuations and non-flow is considered well understood. For a given harmonic n, the average strength of two-particle correlations is determined by forming the following from all pairs: The ϕ values used in the subtraction will originate from different particles to prevent auto-correlations. The single angular brackets denote averaging of particle pairs within the same event. The two-particle cumulant is obtained by averaging 2 over an event ensemble, and is denoted as: In the absence of non-flow, c n {2} provides a measure of v 2 n without the need to measure Ψ n . Respectively, the average strength of four particle correlations is determined by forming the following from all quadruplets: Consequently, the four-particle cumulant is then: The subtraction removes non-flow contributions present in two particle correlations. In the absence of non-flow, c n {4} provides a measure of v 4 n − 2 v 2 n 2 . Respectively, the average strength of six particle correlations is determined by forming the following from all sextuplets: The six-particle cumulant is then: In this case, the subtraction removes non-flow contributions present in two-and four-particle correlations.
In the absence of non-flow, c n {6} provides a measure of v 6 n −9 v 4 n v 2 n +12 v 2 n 3 . As mentioned earlier, the quantities 2 , 4 , or 6 can be determined by averaging over all particles in a given event. The quantities can also be determined using the Q-cumulants of different harmonics, which offers a highly efficient method of evaluating multi-particle correlations without having to consider all combinations [24]. The flow coefficients from two-, four-, and six-particle cumulants can finally be obtained from:  [25,26]. Therefore, the difference in v n {2} and v n {4} can be used to infer the scale of v n fluctuations σ v n . The presence of non-flow influences the cumulants as follows. Assuming large multiplicity events are a superposition of low multiplicity events, the contribution from non-flow (or few-particle correlations) is expected to be diluted as [25]: where M is the multiplicity of the event. Therefore measuring both c n {2}, c n {4}, and c n {6} as a function of multiplicity will help determine whether the underlying correlations are global or few-particle. One can also suppress non-flow by requiring the particles have a relatively large separation in η, since resonances and jets will produce particles with similar rapidity.

Analysis details
The two data sets analyzed were recorded during the p-Pb (in 2013) and the Pb-Pb (in 2010) runs at a center of mass energy of √ s NN = 5.02 TeV and √ s NN = 2.76 TeV, respectively. The Pb-Pb run had equal beam energies giving a nucleon-nucleon center of mass system with rapidity y NN = 0. However, the p-Pb run had different beam energies per nucleon for the p and Pb beam, and resulted in a center of mass system moving in the laboratory frame with y NN = 0.465. All kinematic variables are reported in the laboratory frame. Charged particles are detected using the Time Projection Chamber (TPC), the primary tracking detector of ALICE. The TPC has an angular acceptance of 0 < ϕ < 2π, |η| < 0.9 for tracks with full radial track length (ϕ is the azimuthal angle and η is the pseudo-rapidity), and |η| < 1.5 for tracks of reduced length. Information from the Inner Tracking System (ITS) is used to improve the spatial resolution of TPC tracks, which helps with the rejection of secondary tracks (i.e. not originating from the primary vertex). Primary vertex information is provided by the TPC and the Silicon Pixel Detector (SPD). Two VZERO counters, each containing two arrays of 32 scintillator tiles and covering 2.8 < η < 5.1 (VZERO-A ) and −3.7 < η < −1.7 (VZERO-C), provide information for triggering and event class determination. A more detailed description of the ALICE detector can be found elsewhere [27].
For Pb-Pb collisions, events are selected using a minimum bias trigger, which requires a coincidence of signals in the two VZERO detectors. We use minimum bias and high-multiplicity triggers for p-Pb collisions. As with Pb-Pb, the p-Pb minimum bias trigger, requires a coincidence of two signals from the VZERO detectors, and accepts 99.2% of the non-single diffractive cross section. The high-multiplicity trigger requires a large number of hits in the SPD. Pile-up events are rejected by removing events with multiple vertices, and ensuring the vertices reconstructed from the TPC and SPD agree within 0.8 cm.
After the pile-up rejection procedures, the results are stable with respect to luminosity. Only events with a reconstructed primary vertex within ±10 cm from the center of the detector along the beam axis are used in the analysis to ensure a uniform acceptance in η. The resulting analyzed event sample consisted of about 110M p-Pb and 12M Pb-Pb minimum bias events. In p-Pb collisions, the high-multiplicity trigger allowed for a factor of 10 increase in high-multiplicity events in the top 0.014% of the cross section, compared to the number of minimum bias events. The p-Pb high multiplicity events are used for the last two data points for n = 2, and the last data point for n = 3. Minimum bias events are used for all other points.
The tracks used to determine the cumulants have kinematic cuts 0.2 < p T < 3 GeV/c and |η| < 1. The tracks use an SPD hit if one exists within the trajectory, if not, they are constrained to the primary vertex. Such a configuration leads to a flat ϕ acceptance. It was found that residual non-uniformities influence the cumulant extraction at a level of less than 0.1%. We therefore do not apply acceptance corrections. Track quality is ensured by requiring tracks to have at least 70 TPC clusters out of a maximum of 159, and a χ 2 per TPC cluster less than 4 for the track fit. In addition, the distances of closest approach to the primary vertex in the xy plane and z direction are required to be less than 2.4 cm and 3.2 cm respectively [28].
The results in this article are reported as a function of the corrected multiplicity, N ch . The multiplicity corresponds to the number of charged tracks with 0.2 < p T < 3 GeV/c and |η| < 1, corrected for tracking efficiencies. The tracking efficiency is calculated from a procedure using HIJNG (Pb-Pb) or DPMJET (p-Pb) events [29,30]. GEANT3 is used for transporting simulated particles, followed by a full calculation of the detector response (including production of secondary particles) and track reconstruction done with the ALICE simulation and reconstruction framework [31,32]. The tracking efficiency is ∼ 70% at p T ∼ 0.2 GeV/c and increases to an approximately constant value of ∼ 80% for p T > 1 GeV/c. There are differences on the order of a few percent when comparing between the two collision systems due to the change in detector performance between each run. The final number of particles ( N ch ) is extracted by correcting the raw transverse momentum spectrum with the p T dependent tracking efficiencies. Tables II  and III show multiplicities for the two systems and the fractional cross section.
To reduce the influence of the tracking efficiency on the cumulants (c n {m}), we flatten the p T dependent efficiencies by randomly rejecting high p T particles. These particles have slightly larger efficiencies compared to the low p T ones, so the procedure effectively re-weights the cumulants in favour of low p T particles. This decreases the integrated value of v n by roughly 3%, since v n generally increases with p T . Regarding the choice of multiplicity bin size, it was previously realized that event by event multiplicity fluctuations within a class having a wide multiplicity range can bias the measurement of c n {4}, particularly in the low multiplicity region [26,16]. We avoid this by first extracting c n {m} in unit multiplicity bins (i.e. N ch = 6, 7, 8....). The number of combinations scheme [24], or simple unit event weights gives the same values of c n {m} for unit multiplicity bins. We then average those values to produce c n {m} for larger bin widths, which have a better statistical precision. The following relation is used for averaging procedure; y = ∑ i w i y i ∑ i w i , where y i is the value of the cumulant in a single multiplicity bin, w i corresponds to a choice of weight, and y is the average value obtained from the number of bins in the sum. Monte Carlo studies with known probability density functions (p.d.f.) show that when using unit weights (i.e. w i = 1), our result lies within < 0.1% from the known input y (from the p.d.f.). Other weighting schemes such as w i = M, where M is the multiplicity of the event, or w i = 1/σ 2 i where σ i is the statistical uncertainty of the bin, gave differences of around 2%.
Additional sources of systematic uncertainties in the calculation of c n {m} were extracted by varying the closest approach to the vertex for the tracks, the cut on the minimum number of TPC clusters, the position of the primary vertex and, finally, by analyzing the event sample separately according to the orientation of the magnetic field.
We also generated events with the AMPT model [33] (which includes flow correlations) that were used as an input to our reconstruction simulations. The cumulants obtained directly from the model were compared to those from reconstructed tracks. We found small differences, which are part of the systematic uncertainties. Table 1 summarizes the systematic uncertainties for each collision system. The final systematic uncertainty is calculated by adding all the individual contributions in quadrature.

The second harmonic two-particle cumulant
The results of c 2 {2} as a function of multiplicity are shown in Figs. 1 and 2 for p-Pb and Pb-Pb respectively. The left column presents the results, using the Q-cumulants methods [24] in the case where no ∆η gap is applied. Charge independent refers to the fact that all available charged tracks are used to determine the cumulants. The left panel of Fig. 1 shows that the filled circles (charge independent | < 1) measurements) in p-Pb collisions exhibit a decrease with increasing multiplicity, qualitatively consistent with the expectation of correlations dominated by non-flow effects. When fitting these data points with the function a/M b at large multiplicity, we find b = 0.3. The value b = 1 is expected if high-multiplicity events are a linear superposition of low multiplicity events [25]. This deviation from 1 might indicate the existence of another mechanism that increases c 2 {2}, or that the relative fraction of few particle correlations is increasing with multiplicity. In the same plot, we present measurements of like-sign correlations, calculated by measuring c 2 {2} for positive and negative tracks separately, and forming the average. The corresponding points, represented by the open circles, are lower than the charge independent results for the majority of the multiplicity ranges. This is expected since few-particle correlations from jets and resonances conserve charge, and thus are more likely to be absent in the like-sign measurements. Conversely, the like-sign measurements are higher for the lowest multiplicity bin. This can be explained by a suppression of unlike sign correlations (e.g. multi-particle jets) induced by the low multiplicity cut. Our results in p-Pb collisions are compared to predictions from the DPMJET model [30]. It includes in a phenomenological way the soft multi-particle production as well as hard scatterings, contains no collective effects and thus can serve as a benchmark to study the effect of non-flow on our measurements. It is seen that the corresponding points for c 2 {2} in DPMJET fall off more rapidly compared to data. When carrying out the a/M b fit to the model, we find b ∼ 0.8. The data is also significantly higher than DPMJET at high multiplicity.
The right panel of Fig. 1 presents the multiplicity dependence of the two-particle cumulants in p-Pb collisions in the case where a ∆η gap is applied. It is seen that for a given multiplicity, increasing the gap decreases c 2 {2}. As mentioned previously, this is expected since tracks from few-particle correlations such as jets and resonances have smaller relative angles, therefore their contribution is suppressed by the applied pseudo-rapidity separation. However for large ∆η values, i.e. for |∆η| > 1, the data points increase with multiplicity which is not expected if non-flow dominates. anisotropies decrease, the values of c 2 {2} decrease as expected. Predictions from the HIJING model are also shown in the same plot. This model, similarly to the DPMJET model, contains only non-flow, and as expected, c 2 {2} attenuates more rapidly than the data. Finally, the right panel of Fig. 2 presents the two-particle results in Pb-Pb collisions after applying a ∆η gap to reduce the contribution from nonflow. It is seen that at multiplicities N ch 1000, the measurements with various ∆η gaps converge, indicating the dominance of anisotropic flow. The measurements at lower multiplicities depend on ∆η gap significantly, indicating non-flow plays a prominent role.
In Fig. 3 flow. Both systems have similar values of c 2 {2} at low multiplicity, however the Pb-Pb data points rise more rapidly for higher multiplicities. This maybe explained by higher eccentricities (therefore higher anisotropies) in Pb-Pb collisions found from a CGC inspired cluster model for the initial conditions at similar multiplicities [22] (not shown). We note other studies are exploring these correlations with the AMPT model [34].

The second harmonic four-particle cumulant
The results of c 2 {4} as a function of multiplicity are shown in Fig. 4

The second harmonic six-particle cumulant
The results of c 2 {6} as a function of multiplicity are shown in Fig. 6 for p-Pb and Pb-Pb collisions. We again use the Q-cumulants methods to obtain c 2 {6}. In p-Pb collisions, these measurements are more limited by finite statistics as we observe fluctuations above and below zero at high multiplicity (within the statistical uncertainties). The solid black line indicates v 2 {6} = 4.5%, which is roughly the value of v 2 {4} in this multiplicity region. The p-Pb measurements will benefit from higher statistics measurements planned for future LHC running. However, it is clear at multiplicities above 100, the values of c 2 {6} are significantly higher for Pb-Pb compared to p-Pb. This again maybe be explained by higher eccentricities in the initial state of the colliding nuclei for the former.

Second harmonic cumulants in very high-multiplicity Pb-Pb collisions
The non-zero values of c 2 {4} in high-multiplicity p-Pb collisions merit a comparison to high-multiplicity Pb-Pb collisions, which have an impact parameter that becomes small. In both cases, initial state fluctuations are expected to dominate the eccentricity since there is no intrinsic eccentricity from the overlapping nuclei. In Fig. 7, cumulants of different orders are compared for high-multiplicity Pb-Pb collisions. At N ch 2800, c 2 {4} becomes consistent with zero, which is in contrast to high-multiplicity p-Pb (where c 2 {4} is negative). The measurements of c 2 {6} also become zero in exactly the same region, which corresponds to the highest ∼ 2.5% of the cross section. Constant fits to c 2 {4} and c 2 {6} for N ch > 2800 give 8.5 × 10 −6 ± 9.3 × 10 −6 and 7.2 × 10 −6 ± 2.2 × 10 −5 respectively (with χ 2 /do f ∼ 1 in each case). An explanation for the difference between p-Pb and Pb-Pb can be found by considering the number of sources which form the eccentricity. When this number is small, eccentricity fluctuations have a powerlaw distribution which will lead to finite values of c 2 {4} and c 2 {6}, assuming v 2 ∝ ε 2 [35]. When the number of sources becomes large enough, the power-law distribution becomes equivalent to the Bessel-Gaussian distribution [36,37]. In the special case of very high multiplicity Pb-Pb collisions where the impact parameter is expected to approach 0, the Bessel-Gaussian distribution gives values of c 2 {4} and c 2 {6} that are zero. Assuming the number of sources are highly correlated with the number participants, the difference between very high multiplicity p-Pb and Pb-Pb can be explained by the larger number of sources in the latter. Finally, these results at the LHC can be compared to those from the STAR Collaboration [38,39]. In Au-Au

Second harmonic flow coefficients in p-Pb and Pb-Pb collisions
A comparison of second harmonic flow coefficients is shown in Fig. 8. and v 2 {6} which are close, but not exactly equal [35]. In addition, unfolded measurements of the v 2 distribution have shown Bessel-Gaussian descriptions work reasonably well for Pb-Pb collisions [40,41]. In the left panel of Fig. 9, we show the measurement of R 2 , defined as: As mentioned in Sec. II, when σ v n ≪ v n , R n = σ v n / v n in case non-flow is negligible. In the overlapping multiplicities, the values for p-Pb appear to be higher than Pb-Pb, demonstrating the greater role of fluctuations in the former. A similar observation is reported by the CMS Collaboration [16]. The trend for R 2 in Pb-Pb is similar to observations for Au-Au √ s NN = 200 GeV collisions [38,42]. The value of R 2 in mid-central (mid-multiplicity) Pb-Pb collisions (∼ 0.35) is between the STAR and PHOBOS results for similar centralities. In the right panel, we show σ v2 / v 2 under the assumption that the v 2 distribution is Bessel-Gaussian. Using this assumption, all the information from distribution can be obtained from just v 2 {2} and v 2 {4}, without the need for the condition σ v n ≪ v n [36]. The dashed lines  denote the σ v2 / v 2 = 4/π − 1 limit, expected when fluctuations dominated the eccentricity [43]. We find the Bessel-Gaussian σ v2 / v 2 is close to this limit for high-multiplicity Pb-Pb collisions.

Two-particle cumulants of the third harmonic
In Fig. 10, we show measurements of the third harmonic two-particle cumulants for p-Pb and Pb-Pb collisions, for different values of the ∆η gap. For p-Pb and low Pb-Pb multiplicities, we generally find a strong dependence on the ∆η. The values with small ∆η gaps decrease with multiplicity in p-Pb, as expected when non-flow is dominant. This behavior was also observed by the STAR Collaboration at lower beam energies [11]. The measurements with larger ∆η gaps show an increase with multiplicity, indicating a contribution from global correlations. For large Pb-Pb multiplicities, measurements with various ∆η gaps converge indicating a dominance of flow. Finally, in Fig. 11 we compare the third harmonic flow coefficients for both systems, again with the largest possible ∆η gap. In contrast to measurements of the second harmonic, we find p-Pb and Pb-Pb are consistent for the same multiplicity. This consistency has also been observed by the CMS Collaboration [16], and points to similar third harmonic eccentricities for p-Pb and Pb-Pb at the same multiplicity. A CGC inspired cluster model for the initial conditions is able to reproduce this observation [22].
Measurements of c 2 {2} using all pairs in the event for p-Pb collisions show a decrease with multiplicity, characteristic of a dominance of few-particle correlations. However, the decrease is shallower than from the expectation high-multiplicity events are a superposition of low multiplicity events.  and v 2 {6} are consistent with 0. A comparison of p-Pb cumulants to those of Pb-Pb at the same multiplicity (for N ch 70) shows stronger correlations in Pb-Pb for all the cumulants. This may be explained by higher eccentricities for similar multiplicities. Finally, we have performed measurements of v 3 {2} for p-Pb and Pb-Pb collisions. They are found to be similar for overlapping multiplicities when a |∆η| > 1.4 gap is placed, indicating initial state third harmonic eccentricities may be similar for both systems. We conclude that our measurements indicate the (double) ridge observed in p-Pb at √ s NN = 5.02 TeV arises from global azimuthal correlations, rather than from few-particle correlations which decrease with multiplicity. These measurements provide key constraints to the initial state and transport properties in p-Pb and Pb-Pb collisions.

Uncorrected Corrected
Fractional of hadronic cross section Fraction of hadronic cross section N ch bin N ch within bin above lower bin edge [6,12] 12.0 0.154 0.826 [12,18] 19.5 0.138 0.673 [18,24] 27.  Table 2: Relation of charged track multiplicity N ch to the fraction of hadronic cross section in p-Pb at √ s NN = 5.02 TeV collisions. There is a 3.5% uncertainty in the cross section values. N ch corresponds to the number of charged tracks with 0.2 < p T < 3 GeV/c and |η| < 1. The corrected values of N ch have a systematic uncertainty of 6.0%

Uncorrected Corrected
Fraction of hadronic cross section Fraction of hadronic cross section N ch bin N ch within bin above lower bin edge [6,26] 19  Table 3: Relation of charged track multiplicity N ch to the fraction of hadronic cross section in Pb-Pb at √ s NN = 2.76 TeV collisions. There is a 1% uncertainty in the cross section values. N ch corresponds to the number of charged tracks with 0.2 < p T < 3 GeV/c and |η| < 1. The corrected values of N ch have a systematic uncertainty of 6.0%. 10 Centro de Investigaciones Energéticas Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain