Evidence for collectivity in pp collisions at the LHC

Measurements of two- and multi-particle angular correlations in pp collisions at sqrt(s) = 5, 7, and 13 TeV are presented as a function of charged-particle multiplicity. The data, corresponding to integrated luminosities of 1.0 inverse picobarn (5 TeV), 6.2 inverse picobarns (7 TeV), and 0.7 inverse picobarns (13 TeV), were collected using the CMS detector at the LHC. The second-order (v[2]) and third-order (v[3]) azimuthal anisotropy harmonics of unidentified charged particles, as well as v[2] of K0 short and Lambda/anti-Lambda particles, are extracted from long-range two-particle correlations as functions of particle multiplicity and transverse momentum. For high-multiplicity pp events, a mass ordering is observed for the v[2] values of charged hadrons (mostly pions), K0 short, and Lambda/anti-Lambda, with lighter particle species exhibiting a stronger azimuthal anisotropy signal below pt of about 2 GeV/c. For 13 TeV data, the v[2] signals are also extracted from four- and six-particle correlations for the first time in pp collisions, with comparable magnitude to those from two-particle correlations. These observations are similar to those seen in pPb and PbPb collisions, and support the interpretation of a collective origin for the observed long-range correlations in high-multiplicity pp collisions.


1
In systems such as pp and pPb, where the transverse size of the overlap region is comparable to that of a single proton, the formation of a hot and dense fluid-like medium was not expected. The expectations for other small systems like dAu and 3 HeAu were similar. Various theoretical models have been proposed to interpret the origin of the observed long-range correlations in small collision systems [30][31][32][33][34][35][36][37]. These include initial-state gluon correlations without final-state interactions [33,34] or, similar to what is thought to occur in AA systems, hydrodynamic flow that develops in a conjectured high-density medium [35][36][37].
Owing to the magnitude of the correlation signal, significant progress has been made in unraveling the nature of the ridge correlations in high-multiplicity pPb collisions. Measurements of anisotropy Fourier harmonics (v n ), using identified particles [38,39] and multi-particle correlation techniques [40][41][42][43], reveal features that support a collective, hydrodynamic-like origin of the observed correlations [44][45][46].
In high-multiplicity pp collisions, the nature of the observed long-range correlation still remains poorly understood. Long-range correlations in pp collisions were first observed at √ s = 7 TeV [1], and the extraction of anisotropy v 2 harmonics in pp collisions was recently reported using data at √ s = 13 TeV [2]. A better understanding of the underlying particle correlation mechanisms leading to these observations requires more detailed study of the properties of the v 2 and higher-order harmonics in pp collisions. In particular, their dependence on particle species, and other aspects related to their possible collective nature, are the key to scrutinize various theoretical interpretations. Furthermore, a quantitative modeling of longrange correlations in pPb collisions is found to be highly sensitive to the spatial structure of the proton [35]. Fluctuations of the substructure of nucleons are not well understood, but they can be better constrained with studies of anisotropy harmonics in pp collisions.
This paper extends the characterization of long-range correlation phenomena in high-multiplicity pp collisions by presenting a detailed study of two-and multi-particle azimuthal correlations with unidentified charged particles, as well as correlations of reconstructed K 0 S and Λ/Λ particles at various LHC collision energies. The results of v 2 and v 3 harmonics, extracted from two-particle correlations, are studied as functions of particle p T and event multiplicity. The residual contribution to long-range correlations of back-to-back jet correlations is estimated and removed by subtracting correlations obtained from very low multiplicity pp events. The v 2 harmonics are also extracted using the multi-particle cumulant method to shed light on the possible collective nature of the correlations. The pp results are directly compared to those found for pPb and PbPb systems over a broad range of similar multiplicities.

The CMS detector and data sets
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter. Within the solenoid volume, there are a silicon pixel and strip tracker detector, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. The silicon tracker measures charged particles within the pseudorapidity range |η| < 2.5. It consists of 1440 silicon pixel and 15 148 silicon strip detector modules and is located in the 3.8 T field of the solenoid. For non-isolated particles of 1 < p T < 10 GeV/c and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90  µm in the transverse (longitudinal) impact parameter [47]. Iron and quartz-fiber Cherenkov hadron forward (HF) calorimeters cover the range 2.9 < |η| < 5.2 on either side of the interaction region. These HF calorimeters are azimuthally subdivided into 20 • modular wedges and further segmented to form 0.175×0.175 rad 2 (∆η×∆φ) "towers". The beam pickup for timing (BPTX) devices are designed to provide precise information on the LHC bunch structure and timing of the incoming beams. They are located around the beam pipe at a distance of 175 m from the interaction point on either side. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [48]. The detailed Monte Carlo (MC) simulation of the CMS detector response is based on GEANT4 [49].
The data samples of pp collisions used in this analysis were collected by the CMS experiment in 2010 at √ s = 7 TeV, and in 2015 at 5.02 TeV (labeled as 5 TeV for simplicity) and 13 TeV, with integrated luminosities of 6.2, 1.0, and 0.7 pb −1 , respectively.

Event and track selection
Minimum bias (MB) pp events were triggered by requiring the coincidence of signals from both BPTX devices, indicating the presence of two proton bunches crossing at the interaction point (zero bias condition). The data used in this study were recorded with an average number of pp interactions per bunch crossing ranging from 0.1 to 1.4. Because of hardware limits on the data acquisition rate, only a small fraction (∼10 −3 ) of all MB triggered events were recorded. In order to collect a large sample of high-multiplicity pp collisions, a dedicated trigger was implemented using the CMS level-1 (L1) and high-level trigger (HLT) systems. At L1, the total transverse energy summed over ECAL and HCAL was required to be greater than a given threshold (40, 45 and 55 GeV thresholds are used). Online track reconstruction for the HLT was based on the three layers of pixel detectors, requiring the track origin to be located within a cylindrical region centered on the nominal interaction point with a length of 30 cm along the beam and a radius of 0.2 cm perpendicular to the beam. For each event, the vertex reconstructed with the highest number of pixel tracks was selected. The vertex reconstruction efficiency with high track multiplicities is 100%. The number of pixel tracks (N online trk ) with |η| < 2.4, p T > 0.4 GeV/c, and a distance of closest approach less than 0.12 cm to this vertex, was determined for each event. Data were taken with HLT thresholds of N online trk ≥ 60, 85, 110, seeded with L1 total transverse energy thresholds of 40, 45, and 55 GeV, respectively.
In the offline analysis, hadronic collisions are selected by requiring at least one HF calorimeter tower with more than 3 GeV of total energy in each of the two HF detectors. Events are also required to contain at least one reconstructed primary vertex within 15 cm of the nominal interaction point along the beam axis and within 0.15 cm in the direction transverse to the beam trajectory. Beam related background is suppressed by rejecting events for which less than 25% of all reconstructed tracks pass the high-purity selection (as defined in Ref. [47]). With these selection criteria, 94-96% of the pp interactions simulated with PYTHIA 6 tune Z2 [50] and PYTHIA 8 tune CUETP8M1 [51] event generators that have at least one particle with energy E > 3 GeV in both ranges −5 < η < −3 and 3 < η < 5 are selected.
In this analysis, primary tracks, i.e. tracks that emanate from the primary vertex and that satisfy the high-purity criteria, are used to define the event charged-particle multiplicity and perform correlation measurements. Additional requirements are also applied to enhance the purity of primary tracks. The significance of the separation along the beam axis (z) between the track and the primary vertex, d z /σ(d z ), and the significance of the impact parameter relative to the primary vertex transverse to the beam, d T /σ(d T ), must be smaller than 3. The relative uncertainty of the transverse-momentum measurement, σ(p T )/p T , must be smaller than 10%. To ensure high tracking efficiency and to reduce the rate of misreconstructed tracks, only tracks in the region |η| < 2.4 and p T > 0.3 GeV/c are included. The p T threshold is raised to 0.4 GeV/c for purposes of the event multiplicity determination, to match the requirement in the HLT, while tracks down to 0.3 GeV/c are used in the correlation analysis. Based on simulation studies using GEANT4 to propagate particles from the PYTHIA 8 event generator, the combined geometrical acceptance and efficiency for primary track reconstruction exceeds 60% for p T ≈ 0.3 GeV/c and |η| < 2.4. The efficiency is greater than 90% in the |η| < 1 region for p T > 0.6 GeV/c. For the event multiplicity range studied in this paper, no dependence of the tracking efficiency on the multiplicity is found and the rate of misreconstructed tracks is 1-2%.
Additionally, the CMS loose [47] tracks are also used to incorporate secondary-track candidates with larger track impact parameters, for reconstructing K 0 S and Λ/Λ candidates (also called V 0 candidates). The reconstruction of V 0 candidates in this analysis is identical to that in Refs. [39,52], where more details can be found. Oppositely charged tracks that are detached (having large d z /σ(d z ) and d T /σ(d T ) values) from the primary vertex are selected to determine if they point to a common secondary vertex. The pair of tracks are assumed to be π + π − in K 0 S reconstruction, while the assumption of π − p(π + p) is used in Λ (Λ) reconstruction. The angle θ point between the V 0 momentum vector and the vector connecting the primary and V 0 vertices is required to satisfy cos θ point > 0.999. The K 0 S (Λ/Λ) reconstruction efficiency is about 6% (1%) for p T ≈ 1 GeV/c and 17% (7%) for p T > 3 GeV/c within |η| < 2.4. This efficiency includes the effects of acceptance and the branching ratio for V 0 particle decays into neutral particles. The relatively low reconstruction efficiency of the V 0 candidates is primarily due to the low efficiency for reconstructing daughter tracks with p T < 0.3 GeV/c or larger impact parameters.
Following the same procedure as used in previous analyses [1,4,43], the pp data sets at different energies are divided into classes of reconstructed track multiplicity, N offline trk , where primary tracks with |η| < 2.4 and p T > 0.4 GeV/c are counted. Details of the multiplicity classification in this analysis, including the fractional inelastic cross section and the average number of primary tracks before and after correcting for detector effects in each multiplicity range, are provided in Table 1. Within a given N offline trk range, the average event multiplicity is expected to be larger for higher √ s data, as suggested by the average uncorrected N offline trk values. However, due to a slightly higher tracking efficiency, and hence a smaller efficiency correction, the corrected av-erage multiplicity, N corrected trk , for the same N offline trk range happens to be smaller for 7 TeV data than for 5 TeV data. Uncertainties on N corrected trk values come from systematic uncertainties on the tracking efficiency correction factor estimated from MC simulations.

Analysis technique
The analysis techniques for two-and multi-particle correlations are identical to those used in Refs. [3,4,13,14,39,40,43]. They are briefly summarized in this section for the analysis of the new pp data samples.

Two-particle correlations and Fourier harmonics
For each track multiplicity class, "trigger" particles are defined as charged particles or V 0 candidates with |η| < 2.4 and originating from the primary vertex within a given p trig T range. The number of trigger particles in the event is denoted by N trig . Particle pairs are then formed by associating each trigger particle with the remaining charged primary tracks with |η| < 2.4 and from a specified p assoc T interval (which can be either the same as or different from the p trig T range). A pair is removed if the associated particle is the daughter of any trigger V 0 candidate (this contribution is negligible since associated particles are mostly primary tracks). The twodimensional (2D) per-trigger-particle associated yield is defined in the same way as done in previous analyses, 1 N trig where ∆η and ∆φ are the differences in η and φ of the pair. The same-event pair distribution, S(∆η, ∆φ), represents the yield of particle pairs normalized by N trig from the same event, The mixed-event pair distribution, is constructed by pairing the trigger particles in each event with the associated charged particles from 20 different randomly selected events in the same 0.5 cm wide z vtx range and from the same track multiplicity class. The same-event and mixed-event pair distributions are first calculated for each event, and then averaged over all the events within the track multiplicity class. The ratio B(0, 0)/B(∆η, ∆φ) mainly accounts for pair acceptance effects, with B(0, 0) representing the mixed-event associated yield for both particles of the pair going in approximately the same direction and thus having maximum pair acceptance. Following the procedure described in Refs. [4,13,14,39,43], each reconstructed primary track and V 0 candidate is weighted by a correction factor derived from MC simulations, which accounts for detector effects including the reconstruction efficiency, the detector acceptance, and the fraction of misreconstructed tracks.
The azimuthal anisotropy harmonics of charged particles, K 0 S and Λ/Λ particles can be extracted via a Fourier decomposition of long-range two-particle ∆φ correlation functions, obtained by averaging the 2D two-particle correlation function over |∆η| > 2 (to remove shortrange correlations, such as jet fragmentation), where V n∆ are the Fourier coefficients and N assoc represents the average number of pairs per trigger particle for a given (p trig T , p assoc T ) bin. The first three Fourier terms are included in the fits to the dihadron correlation functions. Including additional terms has a negligible effect on the results of the Fourier fit.
Assuming V n∆ coefficients can be factorized into the product of single-particle anisotropies [43], the elliptic and triangular anisotropy harmonics, v 2 {2, |∆η| > 2} and v 3 {2, |∆η| > 2}, of trigger particles can be extracted as a function of p T from the fitted Fourier coefficients from the twoparticle correlation method, Here, a fixed p ref T range for the "reference" charged primary particles is chosen to be 0.3 < p T < 3.0 GeV/c to minimize correlations from jets at higher p T .
To extract v 2 signal for V 0 candidates, the invariant mass distributions are fitted by a sum of two Gaussian functions with a common mean to describe the true V 0 signal peak, and a fourthorder polynomial function to model the background, as done in Ref. [39]. The v 2 values are first extracted for V 0 candidates from the peak region (which contains contributions from genuine V 0 , as well as background V 0 candidates from random combinatorics) and from a sideband region (which contains only background V 0 s from random combinatorics), denoted as v obs 2 and v bkg 2 . Here the peak region is defined as the mass window of ±2σ around the center of the V 0 candidate mass peak, where σ is found from the addition in quadrature of the standard deviations of the two Gaussian functions weighted by their yields. The sideband region is defined as the mass window outside ±3σ mass range around the V 0 candidate mass peak to upper limit of 0.565 (1.135) GeV and lower limit of 0.430 (1.155) GeV for K 0 S (Λ/Λ) particles. The v 2 signal for V 0 candidates can then be calculated as where f sig represents the signal yield fraction in the peak region determined from the fits to the mass distribution.
Although a requirement of |∆η| > 2 can largely exclude near-side jet-like correlations for v n {2} extraction, contributions from back-to-back (i.e. dijet) correlations are still present in the longrange, away-side (∆φ ≈ π) region, especially for pp collisions. By assuming that the shape of the jet-induced correlations is invariant with event multiplicity, a procedure of removing jet-like correlations in pPb collisions was proposed in Refs. [15,16]. The method consists of subtracting the results for low-multiplicity events, where the ridge signal is not present, from those for high-multiplicity events. For this analysis, a very similar low-multiplicity subtraction method developed for pPb collisions [43] is employed. The Fourier coefficients, V n∆ , extracted from Eq. (4) for 10 ≤ N offline trk < 20 are subtracted from the V n∆ coefficients extracted in the higher-multiplicity region, with . (7) Here, Y jet represents the near-side jet yield obtained by integrating the difference of the shortand long-range event-normalized associated yields for each multiplicity class as shown for 105 ≤ N offline trk < 150 in Fig. 2 (to be described in Section 5.1) over |∆φ| < 1.2. The ratio, Y jet /Y jet (10 ≤ N offline trk < 20), is introduced to account for the enhanced jet correlations resulting from the selection of higher-multiplicity events. This jet subtraction procedure is verified using PYTHIA 6 (Z2) and PYTHIA 8 tune CUETP8M1 pp simulations, where no jet modification from initial-or final-state effects is present. The residual V n∆ after subtraction is found to be consistent with zero. The azimuthal anisotropy harmonics v n after correcting for back-to-back jet correlations estimated from low-multiplicity data (denoted as v sub n ) can be extracted from V sub n∆ using Eq. (5) and (6). In this paper, both the v n and v sub n results are presented.

Fourier harmonics from multi-particle correlations
To avoid a model-dependent correction of jet-like correlations in extracting v n harmonics from two-particle correlations, a multi-particle cumulant analysis using the Q-cumulant method [53] is employed to determine the second-order elliptic harmonic, similar to what was done in pPb and PbPb collisions [40,43]. By simultaneously correlating several (no less than four) particles, the multi-particle cumulant technique has the advantage of suppressing short-range twoparticle correlations such as jets and resonance decays. It also serves as a powerful tool to directly probe the collective nature of the observed azimuthal correlations.
The two-, four-, and six-particle azimuthal correlations [53] are evaluated as: Here φ i (i = 1, . . . , 6) are the azimuthal angles of one unique combination of multiple particles in an event, n is the harmonic number, and · · · represents the average over all combinations from all events within a given multiplicity range. The corresponding cumulants, c n {4} and c n {6}, are calculated as follows [53]: The Fourier harmonics v n that characterize the global azimuthal behavior are related to the multi-particle cumulants [53] using Note that a non-imaginary v n {4} coefficient, which would suggest a bulk medium collective behavior, requires having a negative c n {4} value.

Systematic uncertainties
Systematic uncertainties in this analysis include those from the experimental procedure for obtaining the two-particle v n harmonics, as well as from the jet subtraction procedure.
The experimental systematic effects are evaluated by varying conditions in extracting v sub 2 {2} and v 2 {4} values. The systematic uncertainties are found to have no significant dependence on p T and √ s so they are quoted to be constant percentages over the entire p T range for all collision energies. Experimental systematic uncertainties due to track quality requirements are examined by varying the track selection thresholds for d z /σ(d z ) and d xy /σ(d xy ) from 2 to 5. A comparison of high-multiplicity pp data for a given multiplicity range but collected by two different HLT triggers with different trigger efficiencies is made to study potential trigger biases. The possible contamination of residual pileup events, especially for 5 and 7 TeV data, is also investigated by varying the pileup selection of events in performing the analysis, from no pileup rejection at all to selecting events with only one reconstructed vertex. The sensitivity of the results to the primary vertex position (z vtx ) is quantified by comparing results at different z vtx locations over a 30 cm wide range.
In the jet subtraction procedure for v n {2} measurements, while the factor Y jet /Y jet (10 ≤ N offline trk < 20) accounts for any bias in the magnitude of jet-like associated yield due to multiplicity selection, a change in the ∆φ width of away-side yields could lead to residual jet effects in v n {2} results. This systematic uncertainty is evaluated by integrating the associated yields of the long-range region over a fixed ∆φ window of |∆φ| < π/3 and |∆φ − π| < π/3 on the near and away side. By taking the difference of the yields in two ∆φ windows symmetric around ∆φ = π/2, contributions from the second and fourth Fourier components are cancelled. By choosing the ∆φ window size to be 2π/3, any contribution from the third Fourier component to the near-and away-side associated yields is also cancelled. Any dependence of this yield difference on the event multiplicity (beyond that induced by the Y jet /Y jet (10 ≤ N offline trk < 20) factor) would indicate a modification of jet correlation width in ∆φ. The systematic uncertainty of v n due to this effect is estimated to be 16%, 9%, and 6% for N offline trk < 40, 40 ≤ N offline trk < 85, and N offline trk > 85, respectively. In the same sense, any multiplicity dependence of the ∆η distribution of the away-side would indicate a modification of the jet correlation. The ∆η distribution is investigated in a fixed window |∆φ − π| < π/16 for different N offline trk ranges, resulting in systematic uncertainties of 8%, 3%, and 2.5% for N offline trk < 40, 40 ≤ N offline trk < 85, and N offline trk > 85, respectively. In addition, by separating events in a given multiplicity range into two groups corresponding to the top and bottom 30% in the leading particle p T distribution, jet correlations are either strongly enhanced or suppressed in a controlled manner. After applying the subtraction procedure, the v n results for the two event groups are consistent within 5%. Table 2 summarizes various sources of systematic uncertainties in this analysis for multiplicitydependent results. The same sources apply to p T differential results, leading to total experimental systematic uncertainty of 5% and method uncertainties of 9%, 13%, 23%, and 37% for p   Table 2: Summary of systematic uncertainties for multiplicity-dependent v sub n {2} from twoparticle correlations (after correcting for jet correlations), and v 2 {4}, v 2 {6} from multi-particle correlations in pp collisions. Different multiplicity ranges are represented as [m, n). Systematic uncertainties originating from different sources are added in quadrature to obtain the overall systematic uncertainty shown as boxes in the figures. No energy dependence has been observed for the systematic uncertainties, therefore, they are only shown for 13 TeV results. Because of insufficient statistical precision, the uncertainties in v 3 are assumed to be the same as those in v 2 , as was done in Refs. [39,43]. Systematic uncertainties of v 2 {2} results for V 0 particles are taken from Ref. [39] in pPb collisions except for the HLT trigger bias and jet subtraction method, which are quoted to be the same as for the inclusive charged particles. Different particle species have different η distributions, which can affect the comparison of results if there is a strong η dependence. This effect is found to be negligible by comparing v 2 {2} results for V 0 particles with different reconstruction efficiency corrections for the η distribution. The relative systematic uncertainties for the two-particle V n∆ coefficients as a function of N offline trk in Fig. 4 (described in Section 5.2) are exactly twice those for the corresponding v n harmonics, since V n∆ = v 2 n when trigger and associated particles are selected from the same p T range. In the same way, relative systematic uncertainties for multi-particle c 2 {m} measurements as a function of N offline trk in Fig. 9 (described in Section 5.3) are exactly m times those for the corresponding v 2 {m} harmonics, where m = 4 or 6.  Figure 1 shows the 2D ∆η-∆φ correlation functions, for pairs of a charged (top), a K 0 S (middle), or a Λ/Λ (bottom) trigger particle with a charged associated particle, in low-multiplicity (10 ≤ N offline trk < 20, left) and high-multiplicity (105 ≤ N offline trk < 150, right) pp collisions at √ s = 13 TeV. Both trigger and associated particles are selected from the p T range of 1-3 GeV/c. For all three types of particles at high multiplicity, in addition to the correlation peak near (∆η, ∆φ) = (0, 0) that results from jet fragmentation, a long-range ridge structure is seen at ∆φ ≈ 0 extending at least 4 units in |∆η|, while such a structure is not observed in low multiplicity events. On the away-side (∆φ ≈ π) of the correlation functions, a long-range structure is also seen and found to exhibit a much larger magnitude compared to that on the near-side for this p T range. This away-side correlation structure contains contributions from back-to-back jets, which need to be accounted for before extracting any other source of correlations.

Two-particle correlation functions
To investigate the observed correlations in finer detail, the 2D distributions shown in Fig. 1 are reduced to one-dimensional (1D) distributions in ∆φ by averaging over |∆η| < 1 (defined as the "short-range region") and |∆η| > 2 (defined as the "long-range region"), respectively, as done in Refs. [1,4,13,14]. Figure 2 shows examples of 1D ∆φ correlation functions for trigger particles composed of inclusive charged particles (left), K 0 S particles (middle), and Λ/Λ particles (right), in the multiplicity range 10 ≤ N offline trk < 20 (open symbols) and 105 ≤ N offline trk < 150 (filled symbols). The curves show the Fourier fits from Eq. (4) to the long-range region, which will be discussed in detail in Section 5.2. To represent the correlated portion of the associated yield, each distribution is shifted to have zero associated yield at its minimum following the standard zero-yield-at-minimum (ZYAM) procedure [43]. An enhanced correlation at ∆φ ≈ 0 in the long-range region is observed for 105 ≤ N offline trk < 150, while such a structure is not presented for 10 ≤ N offline trk < 20. As illustrated in Fig. 1 (right), the near-side long-range ridge structure remains nearly constant in ∆η. Therefore, the near-side jet correlation can be extracted by taking a difference of 1D ∆φ projections between the short-and long-range regions, as shown in the bottom panels in Fig. 2.  After subtracting the results, with the ZYAM procedure applied, for low-multiplicity 10 ≤ N offline trk < 20 scaled by Y jet /Y jet (10 ≤ N offline trk < 20) as in Eq. (7), the long-range 1D ∆φ correlation functions in the high-multiplicity range 105 ≤ N offline trk < 150 for pp collisions at √ s = 13 TeV are shown in Fig. 3, for trigger particles composed of inclusive charged particles (left), K 0 S (middle), and Λ/Λ (right) particles. A "double-ridge" structure on the near and away side is observed after subtraction of jet correlations. The shape of this structure, which is dominated by a second-order Fourier component, is similar to what has been observed in pPb [4][5][6][7] and PbPb [13][14][15][16][17]

Two-particle fourier harmonics v n
Fourier coefficients, V n∆ , extracted from 1D ∆φ two-particle correlation functions for the longrange ∆η region using Eq. (4), are first studied for inclusive charged hadrons. Figure 4 shows the V 2∆ and V 3∆ values for pairs of inclusive charged particles averaged over 0.3 < p T < 3.0 GeV/c as a function of multiplicity in pp collisions at √ s = 13 TeV, before and after correcting for back-to-back jet correlations estimated from low-multiplicity data (10 ≤ N offline trk < 20). The V n∆ results for 5 and 7 TeV are equal to the 13 TeV results within the uncertainties.
Before corrections, the V 2∆ coefficients are found to remain relatively constant as a function of multiplicity. This behavior is very different from the PYTHIA 8 tune CUETP8M1 MC simulation, where the only source of long-range correlations is back-to-back jets and the V 2∆ coefficients decrease with N offline trk . The V 3∆ coefficients found using the PYTHIA 8 simulation are always negative because of dominant contributions at ∆φ ≈ π from back-to-back jets [15], with their magnitudes decreasing as a function of N offline trk . A similar trend is seen in the data for the low multiplicity region, N offline trk < 90. However, for N offline trk ≥ 90, the V 3∆ coefficients in pp data change to positive values. This transition directly indicates a new phenomena that is not present in the PYTHIA 8 simulation. After applying the jet correction procedure detailed in Section 4.1, V 2∆ exhibits an increase with multiplicity for N offline trk 100, and reaches a relatively constant value for the higher N offline trk region. The V 3∆ values after subtraction of jet correlations become positive over the entire multiplicity range and increase with multiplicity.
The elliptic (v 2 ) and triangular (v 3 ) flow harmonics for charged particles with 0.3 < p T < 3.0 GeV/c, after applying the jet correction procedure, are then extracted from the two-particle Fourier coefficients obtained using Eq. (5), and are shown in Fig. 5   The v 2 values reported by the ATLAS collaboration for pp collisions at √ s = 13 TeV in Ref. [2] have a different multiplicity dependence than the results presented in this paper. A nearly constant v 2 value is observed over the entire multiplicity range. This distinct difference, especially in the low multiplicity region, is rooted in the different approaches employed in identifying the v 2 signal from jet-like correlations. In the method from CMS and also previous ATLAS and ALICE analyses [5,6,54], v 2 is always extracted with respect to all of the particles in each event. As seen in Eq. (4), the Fourier coefficients in the current analysis represent an oscillation multiplied by the full N assoc . In the new approach by the ATLAS collaboration [2], v 2 is extracted with respect to a subset of particles in each event. In Ref. [2], their equivalent of our Eq. (4) uses a smaller number than the full N assoc in the events, thereby assuming that some of the particles do not participate in the underlying processes producing the observed azimuthal correlations. As a result, using the method of Ref. [2], a cos(2∆φ) modulation with exactly the same amplitude would yield a bigger V n∆ parameter compared to that found using Eq. (4). This, in turn, leads to larger v 2 values comparing to results obtained with respect to all of the particles. The difference between the two methods becomes larger as N offline trk decreases. It was checked that when applying exactly the same kinematic selections and analysis methods, no discrepancy is found between the two experiments. In the study of v 2 from multiparticle correlations, as will be discussed in Section 5.3, the v 2 is always considered with respect to all the particles in the event for each multiplicity class, which is consistent with the method used in this paper to extract v 2 from two-particle correlations.
The v 2 results as a function of p T for high-multiplicity pp events at √ s = 5, 7, and 13 TeV are shown in Fig. 6 before (left) and after (right) correcting for jet correlations. To compare results with similar average N offline trk , 105 ≤ N offline trk < 150 is chosen for 13 TeV while 110 ≤ N offline trk < 150 is chosen for 5 and 7 TeV. Little energy dependence is observed for the p Tdifferential v 2 results, especially before correcting for jet correlations, as shown in Fig. 6 (left). This conclusion also holds after jet correction procedure for v sub 2 results (Fig. 6, right) within systematic uncertainties, although systematic uncertainties for v sub 2 are significantly higher at high p T because of the large magnitude of the subtracted term. This observation is consistent with the energy independence of associated long-range yields on the near-side reported in Ref. [3]. The observed p T dependence of v sub 2 , in high-multiplicity pp events with peak values at 2-3 GeV/c at various energies, is similar to that in pPb [38,43,54] and PbPb [14,55,56] collisions.   < 150, after correcting for back-to-back jet correlations estimated from low-multiplicity data. Bottom: the n q -scaled v sub 2 results for K 0 S and Λ/Λ particles as a function of KE T /n q . Ratios of v sub 2 /n q for K 0 S and Λ/Λ particles to a smooth fit function of data for K 0 S particles are also shown. The error bars correspond to the statistical uncertainties, while the shaded areas denote the systematic uncertainties.
The dependence of the elliptic flow harmonic on particle species can shed further light on the nature of the correlations. The v 2 data as a function of p T for identified K 0 S and Λ/Λ particles are extracted for pp collisions at √ s = 13 TeV. Figure 7 shows the results for a low (10 ≤ N offline trk < 20) and a high (105 ≤ N offline trk < 150) multiplicity range before applying the jet correction procedure.
At low multiplicity (Fig. 7 left), the v 2 values are found to be similar for charged particles, K 0 S and Λ/Λ hadrons across most of the p T range within statistical uncertainties, similar to the observation in pPb collisions at √ s NN = 5 TeV [39]. This would be consistent with the expectation that back-to-back jets are the dominant source of long-range correlations on the away side in low-multiplicity pp events. Moving to high-multiplicity pp events (105 ≤ N offline trk < 150, Fig. 7 right), a clear deviation of v 2 among various particle species is observed. In the lower p T region of 2.5 GeV/c, the v 2 value of K 0 S is greater than that of Λ/Λ at a given p T value. Both are consistently below the inclusive charged particle v 2 values. Since most charged particles are pions in this p T range, this indicates that lighter particle species exhibit a stronger azimuthal anisotropy signal. A similar trend was first observed in AA collisions at RHIC [57,58], and later also seen in pPb collisions at the LHC [38,39]. This behavior is found to be qualitatively consistent with hydrodynamic models [44,45]. At p T > 2.5 GeV/c, the v 2 values of Λ/Λ particles become greater than those of K 0 S particles. This reversed ordering of K 0 S and Λ/Λ at high p T is similar to what was previously observed in pPb and PbPb collisions [39].
After applying the correction for jet correlations, the v sub 2 results as a function of p T for 105 ≤ N offline trk < 150 are shown in Fig. 8 (top) for the identified particles and charged hadrons. The v sub 2 values for all three types of particles are found to increase with p T , reaching 0.08-0.10 at 2 < p T < 3 GeV/c, and then to decrease for higher p T values. The particle mass ordering of v 2 values in the lower p T region is also observed after applying jet correction procedure, while at higher p T the ordering is reversed. As done in Ref. [39], the scaling behavior of v sub 2 divided by the number of constituent quarks, n q , as a function of transverse kinetic energy per quark, KE T /n q , is investigated for high-multiplicity pp events in Fig. 8 (bottom). The dashed curve corresponds to a polynomial fit to the K 0 S data. The ratio of n q -scaled v sub 2 results for K 0 S and Λ/Λ particles divided by this polynomial function fit is also shown in Fig. 8 (bottom). An approximate scaling is seen for KE T /n q 0.2 GeV within about ±10%. To further reduce the residual jet correlations on the away side and explore the possible collective nature of the long-range correlations, a four-and six-particle cumulant analysis is used to extract the elliptic flow harmonics, v 2 {4} and v 2 {6}. The four-particle cumulant c 2 {4} values for charged particles with 0.3 < p T < 3.0 GeV/c are shown in Fig. 9 (left), as a function of N offline trk for pp collisions at √ s = 5, 7, and 13 TeV. The pPb data at √ s NN = 5 TeV [43] are also plotted for comparison. The six-particle cumulant c 2 {6} values for pp collisions at √ s = 13 TeV are shown in Fig. 9 (right), compared with pPb data at √ s NN = 5 TeV [43]. Due to statistical limitations,

Multi-particle correlations and collectivity
values are only derived for high multiplicities (i.e., N offline trk ≈ 100) for 13 TeV pp data.
The c 2 {4} values for pp data at all energies show a decreasing trend with increasing multiplicity, similar to that found for pPb collisions. An indication of energy dependence of c 2 {4} values is seen in Fig. 9 (left), where c 2 {4} tends to be more positive for a given N offline trk range at lower √ s energies. As average p T values are slightly smaller at lower collision energies, the observed energy dependence may be related to smaller negative contribution to c 2 {4} from smaller p T -averaged v 2 {4} signals. In addition, when selecting from a fixed multiplicity range, a larger positive contribution to c 2 {4} from larger jet-like correlations in the much rarer highmultiplicity events in lower energy pp collisions can also result in an energy dependence. At The  Fig. 10, as a function of event multiplicity. The v 2 data obtained from long-range two-particle correlations after correcting for jet correlations (v sub 2 ) are also shown for comparison. Within experimental uncertainties, the multi-particle cumulant v 2 {4} and v 2 {6} values in highmultiplicity pp collisions are consistent with each other, similar to what was observed previously in pPb and PbPb collisions [40]. This provides strong evidence for the collective nature of the long-range correlations observed in pp collisions. However, unlike for pPb and PbPb colli-sions where v 2 {2} values show a larger magnitude than multi-particle cumulant v 2 results, the v 2 values obtained from two-, four-, and six-particle correlations are comparable in pp collisions at √ s = 13 TeV within uncertainties. In the context of hydrodynamic models, the relative ratios of v 2 among two-and various orders of multi-particle correlations provide insights to the details of initial-state geometry fluctuations in pp and pPb systems. As shown in Ref. [46], the ratio of v 2 {4} to v 2 {2} is related to the total number of fluctuating sources in the initial stage of a collision. The comparable magnitudes of v 2 {2} and v 2 {4} signals observed in pp collisions, compared to pPb collisions at similar multiplicities, may indicate a smaller number of initial fluctuating sources that drive the long-range correlations seen in the final state. Meanwhile, it remains to be seen whether other proposed mechanisms [32][33][34] in interpreting the long-range correlations in pPb and PbPb collisions can also describe the features of multi-particle correlations seen in pp collisions.

Summary
The CMS detector has been used to measure two-and multi-particle azimuthal correlations with K 0 S , Λ/Λ and inclusive charged particles over a broad pseudorapidity and transverse momentum range in pp collisions at √ s = 5, 7, and 13 TeV. With the implementation of highmultiplicity triggers during the LHC 2010 and 2015 pp runs, the correlation data are explored over a broad particle multiplicity range. The observed long-range (|∆η| > 2) correlations are quantified in terms of azimuthal anisotropy Fourier harmonics (v n ). The elliptic (v 2 ) and triangular (v 3 ) flow Fourier harmonics are extracted from long-range two-particle correlations. After subtracting contributions from back-to-back jet correlations estimated using low-multiplicity data, the v 2 and v 3 values are found to increase with multiplicity for N offline trk 100, and reach a relatively constant value at higher values of N offline trk . The p T dependence of the v 2 harmonics in high-multiplicity pp events is found to have no or very weak dependence on the collision energy. In low-multiplicity events, similar v 2 values as a function of p T are observed for inclusive charged particles, K 0 S and Λ/Λ, possibly reflecting a common back-to-back jet origin of the correlations for all particle species. Moving to the higher-multiplicity region, a particle species dependence of v 2 is observed with and without correcting for jet correlations. For p T 2 GeV/c, the v 2 of K 0 S is found to be larger than that of Λ/Λ. This behavior, which is consistent with predictions of hydrodynamic models, is similar to what was previously observed for identified particles produced in pPb and AA collisions at RHIC and the LHC. This mass ordering is reversed at higher p T values. Finally, v 2 signals based on four-and six-particle correlations are observed for the first time in pp collisions. The v 2 values obtained with two-, four-, and six-particle correlations at √ s = 13 TeV are found to be comparable within uncertainties. These observations provide strong evidence supporting the interpretation of a collective origin for the observed long-range correlations in high-multiplicity pp collisions.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW