Measurements of long-range azimuthal anisotropies and associated Fourier coefficients for $pp$ collisions at $\sqrt{s}=5.02$ and $13$ TeV and $p$+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with the ATLAS detector

ATLAS measurements of two-particle correlations are presented for $\sqrt{s} = 5.02$ and $13$ TeV $pp$ collisions and for $\sqrt{s_{\mathrm{NN}}} = 5.02$ TeV $p$+Pb collisions at the LHC. The correlation functions are measured as a function of relative azimuthal angle $\Delta \phi$, and pseudorapidity separation $\Delta \eta$, using charged particles detected within the pseudorapidity interval $|\eta|{<}2.5$. Azimuthal modulation in the long-range component of the correlation function, with $|\Delta\eta|{>}2$, is studied using a template fitting procedure to remove a"back-to-back"contribution to the correlation function that primarily arises from hard-scattering processes. In addition to the elliptic, $\cos{(2\Delta\phi)}$, modulation observed in a previous measurement, the $pp$ correlation functions exhibit significant $\cos{(3\Delta\phi)}$ and $\cos{(4\Delta\phi)}$ modulation. The Fourier coefficients $v_{n,n}$ associated with the $\cos{(n\Delta\phi)}$ modulation of the correlation functions for $n =$ $2$-$4$ are measured as a function of charged-particle multiplicity and charged-particle transverse momentum. The Fourier coefficients are observed to be compatible with $\cos{(n\phi)}$ modulation of per-event single-particle azimuthal angle distributions. The single-particle Fourier coefficients $v_n$ are measured as a function of charged-particle multiplicity, and charged-particle transverse momentum for $n {=} $ $2$-$4$. The integrated luminosities used in this analysis are, $64$ $\mathrm{nb^{-1}}$ for the $\sqrt{s}=13$ TeV $pp$ data, $170$ $\mathrm{nb^{-1}}$ for the $\sqrt{s}=5.02$ TeV $pp$ data and $28$ $\mathrm{nb^{-1}}$ for the $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV $p$+Pb data.


Introduction
Observations of azimuthal anisotropies in the angular distributions of particles produced in proton-lead (p+Pb) collisions at the LHC [1][2][3][4][5] and in deuteron-gold (d+Au) [6-8] and 3 He+Au [9] collisions at RHIC have garnered much interest due to the remarkable similarities between the phenomena observed in those colliding systems and the effects of collective expansion seen in the Pb+Pb and Au+Au collisions [3, 10-13]. 1 The most intriguing feature of the azimuthal anisotropies is the "ridge": an enhancement in the production of particles with small azimuthal angle (φ) separation which extends over a large range of pseudorapidity (η) separation [1,2,14,15]. In Pb+Pb [3, 10-13] and p+Pb [1][2][3] collisions, the ridge is understood to result from sinusoidal modulation of the single-particle azimuthal angle distributions, and the characteristics of the modulation, for example the p T dependence [16], are remarkably similar in the two systems [4].
While the modulation of the azimuthal angle distributions in Pb+Pb collisions is understood to result from the geometry of the initial state and the imprinting of that geometry on the angular distributions of the particles by the collective expansion (see e.g. [17][18][19] and references therein), there is, as yet, no consensus that the modulation observed in p+Pb collisions results from the same mechanism. Indeed, an alternative explanation for the modulation using perturbative QCD and assuming saturated parton distributions in the lead nucleus is capable of reproducing many features of the p+Pb data [20][21][22][23][24][25][26][27][28][29]. Nonetheless, because of the many similarities between the p+Pb and Pb+Pb observations, extensive theoretical and experimental effort has been devoted to address the question of whether the strong-coupling physics understood to be responsible for the collective dynamics in A+A collisions may persist in smaller systems [30][31][32][33][34][35][36][37][38][39][40].
A recent study by the ATLAS Collaboration of two-particle angular correlations in proton-proton (pp) collisions at center-of-mass energies of √ s = 13 and 2.76 TeV obtained results that are consistent with the presence of an elliptic or cos (2φ) modulation of the per-event single particle azimuthal angle distributions [41]. This result suggests that the ridge previously observed in √ s = 7 TeV pp collisions [14] results from modulation of the single-particle azimuthal angle distributions similar to that seen in Pb+Pb and p+Pb collisions. Indeed, the p T dependence of the modulation was similar to that observed in the other systems. Unexpectedly, the amplitude of the modulation relative to the average differential particle yield dN/dφ , was observed to be constant, within uncertainties, as a function of the charged particle multiplicity of the pp events and to be consistent between the two energies, suggesting that the modulation is an intrinsic feature of high-energy pp collisions. These results provide further urgency to address the question of whether strong coupling and collective dynamics play a significant role in small systems, including the smallest systems accessible at collider energiespp collisions. Since the elliptic modulation observed in the pp data is qualitatively similar to that seen in p+Pb collisions, a direct, quantitative comparison of pp and p+Pb measurements is necessary for evaluating whether the phenomena are related.
The modulation of the single-particle azimuthal angle distributions in A+A, p/d+A, and, most recently, pp collisions is usually characterized using a set of Fourier coefficients v n , that describe the relative amplitudes of the sinusoidal components of the single-particle distributions. More explicitly, the azimuthal angle distributions of the particles are parameterized according to: where the average in the equation indicates an average over azimuthal angle. Here, Ψ n represents one of the n angles at which the nth-order harmonic is maximum; it is frequently referred to as the event-plane angle for the nth harmonic. In Pb+Pb collisions, n = 2 modulation is understood to primarily result from an elliptic anisotropy of the initial state for collisions with non-zero impact parameter; that anisotropy is subsequently imprinted onto the angular distributions of the produced particles by the collective evolution of the medium, producing an elliptic modulation of the produced particle azimuthal angle distributions in each event [17,42,43]. The higher (n > 2) harmonics are understood to result from position-dependent fluctuations in the initial-state energy density which produce higher-order spatial eccentricities that similarly get converted into sinusoidal modulation of the single-particle dN/dφ distribution by the collective dynamics [44][45][46][47][48][49][50][51]. Significant v n values have been observed in Pb+Pb (p+Pb) collisions up to n = 6 [13] (n = 5 [4]). An important, outstanding question is whether n > 2 modulation is present in pp collisions.
The v n,n coefficients can be measured using two-particle angular correlation functions, which, when evaluated as a function of ∆φ ≡ φ a − φ b , where a and b represent the two particles used to construct the correlation function, have an expansion similar to that in Eq. (1): If the modulation of the two-particle correlation function arises solely from the modulation of the singleparticle distributions, then, v n,n = v 2 n . Often, the two-particle correlations are measured using different transverse momentum (p T ) ranges for particles a and b. Since the modulation is observed to vary with p T , then v n,n (p a T , if the modulation of the correlation function results solely from single-particle modulation. 2 This "factorization" hypothesis can be tested experimentally by measuring v n,n (p a T , p b T ) for different ranges of p b T and estimating v n (p a T ) using v n (p a T ) = v n,n (p a T , p b T )/ v n,n (p b T , p b T ) (4) and evaluating whether v n (p a T ) depends on the choice of p b T . In addition to the sinusoidal modulation, the two-particle correlation functions include contributions from hard-scattering processes that produce a jet peak centered at ∆φ = ∆η = 0 and a dijet enhancement at ∆φ = π that extends over a wide range of ∆η. The jet peak can be avoided by studying the long-range part of the correlation function, which is typically chosen to be |∆η| > 2. Because the dijet contribution to the two-particle correlation function is not localized in ∆η, that contribution has to be subtracted from the measured correlation function, typically using the correlation function measured in low-multiplicity ("peripheral") events. Different peripheral subtraction methods have been applied for the p+Pb measurements in the literature [2,4]; all of them relied on the "zero yield at minimum" (ZYAM) [2,4] hypothesis to subtract an assumed flat combinatoric component from the peripheral reference correlation function. These methods were found to be inadequate for pp collisions, where the amplitude of the dijet enhancement at ∆φ = π is much larger than the (absolute) amplitude of the sinusoidal modulation. For the measurements in Ref.
[41], a template fitting method, described below, was developed which is better suited for extracting a small sinusoidal modulation from the data. Application of the template fitting method to the pp data provided an excellent description of the measured correlation functions. It also indicated substantial bias resulting from the application of the ZYAM-subtraction procedure to the peripheral reference correlation function due to the non-zero v 2,2 in low-multiplicity events. As a result, the measurements presented in Ref.
[41] were obtained without using ZYAM subtraction. However, the previously published p+Pb data [4] may be susceptible to an unknown bias due to the use of the ZYAM method. Thus, a reanalysis of the p+Pb data is both warranted and helpful in making comparisons between pp and p+Pb data.
To address the points raised above, this paper extends previous measurements of two-particle correlations in pp collisions at √ s = 13 TeV using additional data acquired by ATLAS subsequent to the measurements in Ref.
[41] and provides new measurements of such correlations in pp collisions at √ s = 5.02 TeV. It also presents a reanalysis of two-particle correlations in 5.02 TeV p+Pb collisions and presents a direct comparison between the pp and p+Pb data at the same per-nucleon center-of-mass energy as well as a comparison between the pp data at the two energies. Two-particle Fourier coefficients v n,n are measured, where statistical precision allows, for n = 2, 3, and 4 as a function of charged-particle multiplicity and transverse energy. Measurements are performed for different p a T and p b T intervals and the factorization of the resulting v n,n values is tested. This paper is organized as follows. Section 2 gives a brief overview of the ATLAS detector subsystems and triggers used in this analysis. Section 3 describes the data sets, and the offline selection criteria used to select events and reconstruct charged-particle tracks. The variables used to characterize the "event activity" of the pp and p+Pb collisions are also described. Section 4 gives details of the two-particle correlation method. Section 5 describes the template fitting of the two-particle correlations, which was originally developed in Ref. [41]. The template fits are used to extract the Fourier harmonics v n,n (Eq. (2)) of the long-range correlation, and the factorization of the v n,n into single-particle harmonics v n (Eq. (3)) is studied. The stability of the v n,n as a function of the pseudorapidity separation between the charged-particle pairs is also checked. Section 6 describes the systematic uncertainties associated with the measured v n,n . Section 7 presents the main results of the analysis, which are the p T and event-activity dependence of the single-particle harmonics, v n . Detailed comparisons of the v n between the three data sets: 13 TeV pp, 5.02 TeV pp, and 5.02 TeV p+Pb are also shown. Section 8 gives a summary of the main results and observations.

Trigger
The ATLAS trigger system [58] consists of a Level-1 (L1) trigger implemented using a combination of dedicated electronics and programmable logic, and a software-based high-level trigger (HLT). Due to the large interaction rates, only a small fraction of minimum-bias events could be recorded for all three data sets. The configuration of the minimum-bias (MB) triggers varied between the different data sets. Minimum-bias p+Pb events were selected by requiring a hit in at least one MBTS counter on each side (MBTS_1_1) or a signal in the ZDC on the Pb-fragmentation side with the trigger threshold set just below the peak corresponding to a single neutron. In the 13 TeV pp data, MB events were selected by a L1 trigger that requires a signal in at least one MBTS counter (MBTS_1). In the 5.02 TeV pp data, MB events were selected using the logical OR of the MBTS_1, MBTS_1_1, and a third trigger that required at least one reconstructed track at the HLT. In order to increase the number of events having high charged-particle multiplicity, several high-multiplicity (HMT) triggers were implemented. These apply a L1 requirement on either the transverse energy (E T ) in the calorimeters or on the number of hits in the MBTS, and an HLT requirement on the multiplicity of HLT-reconstructed charged-particle tracks. That multiplicity, N HLT trk , is evaluated for tracks having p T >0.4 GeV that are associated with the reconstructed vertex with the highest multiplicity in the event. This last requirement suppresses the selection of events with multiple collisions (pileup), as long as the collision vertices are not so close as to be indistinguishable. The HMT trigger configurations used in this analysis are summarized in Table 1.

Data sets
The √ s = 13 and 5.02 TeV pp data were collected during Run 2 of the LHC. The 13 TeV pp data were recorded over two periods: a set of low-luminosity runs in June 2015 (used in Ref. [41]) for which the number of collisions per bunch crossing, µ, varied between 0.002 and 0.04, and a set of intermediateluminosity runs in August 2015 where µ varied between 0.05 and 0.6. The 5.02 TeV pp data were recorded during November 2015 in a set of intermediate-luminosity runs with µ of ∼1.5. The p+Pb data were recorded in Run 1 during p+Pb operation of the LHC in January 2013. During that period, the LHC was configured with a 4 TeV proton beam and a 1.57 TeV per-nucleon Pb beam that together produced collisions at √ s NN =5.02 TeV. The higher energy of the proton beam produces a net rapidity shift of the nucleon-nucleon center-of-mass frame by 0.47 units in the proton-going direction, relative to the ATLAS reference system. The p+Pb data were collected in two periods between which the directions of the proton and lead beams were reversed.

Event and track selection
In the offline analysis, additional requirements are imposed on the events selected by the MB and HMT triggers. The events are required to have a reconstructed vertex with the z-position of the vertex restricted to ±150 mm. In the p+Pb data, non-collision backgrounds are suppressed by requiring at least one hit in a MBTS counter on each side of the interaction point, and the time-difference measured between the two sides of the MBTS to be less than 10 ns. In the 2013 p+Pb run, the luminosity conditions provided by the LHC resulted in an average probability of 3% for pileup events. The pileup events are suppressed by rejecting events containing more than one good reconstructed vertex. The remaining pileup events are further suppressed using the number of detected neutrons, N n , measured in the ZDC on the Pbfragmentation side. The distribution of N n in events with pileup is broader than that for the events without pileup. Hence, rejecting events at the high tail end of the ZDC signal distribution further suppresses the pileup, while retaining more than 98% of the events without pileup. In the pp data, pileup is suppressed by only using tracks associated with the vertex having the largest p 2 T , where the sum is over all tracks associated with the vertex. Systematic uncertainties in the measured v n associated with the residual pileup are estimated in Section 6.
In the p+Pb analysis, charged-particle tracks are reconstructed in the ID using an algorithm optimized for pp minimum-bias measurements [59]. The tracks are required to have p T >0.4 GeV and |η|<2.5, at least one pixel hit, with the additional requirement of a hit in the first pixel layer when one is expected, 4 and at least six SCT hits. In addition, 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. They are also required to satisfy |d 0 |/σ d 0 <3 and |z 0 sin(θ)|/σ z 0 sin(θ) <3, where σ d 0 and σ z 0 sin(θ) are uncertainties in d 0 and z 0 sin(θ), respectively.
In the pp analysis, charged-particle tracks and primary vertices are reconstructed in the ID using an algorithm similar to that used in Run 1, but substantially modified to improve performance [60,61]. The reconstructed tracks are required to satisfy the following selection criteria: p T >0.4 GeV and |η|<2.5; at least one pixel hit, with the additional requirement of a hit in the IBL if one is expected (if a hit is not expected in the IBL, a hit in the next pixel layer is required if such a hit is expected); a minimum of six hits in the SCT; |d 0 |<1.5 mm and |z 0 sin(θ)|<1.5 mm. 5 Finally, in order to remove tracks with mismeasured p T due to interactions with the material or other effects, the track-fit χ 2 probability is required to be larger than 0.01 for tracks having p T > 10 GeV.
The efficiencies (p T , η) of track reconstruction for the above track selection cuts are obtained using Monte Carlo (MC) generated events that are passed through a GEANT4 [62] simulation [63] of the ATLAS detector response and reconstructed using the algorithms applied to the data. For determining the p+Pb efficiencies, the events are generated with version 1.38b of the HIJING event generator [64] with a center-of-mass boost matching the beam conditions. For determining the pp efficiencies, nondiffractive 13 TeV pp events obtained from the Pythia 8 [65] event generator (with the A2 set of tuned parameters [66] and the MSTW2008LO PDFs [67]) are used. Both the pp and p+Pb efficiencies increase by ∼3% from 0.4 GeV to 0.6 GeV and vary only weakly with p T for p T >0.6 GeV. In the p+Pb case, the efficiency at p T ∼0.6 GeV ranges from 81% at η=0 to 73% at |η|=1.5 and 65% at |η|>2.0. The efficiency is also found to vary by less than 2% over the multiplicity range used in the analysis. In the pp case, the efficiency at p T ∼0.6 GeV ranges from 87% at η=0 to 76% at |η|=1.5 and 69% for |η|>2.0.

Event-activity classes
As in previous ATLAS analyses of long-range correlations in p+Pb [2,4] and pp [41] collisions, the event activity is quantified by N rec ch : the total number of reconstructed charged-particle tracks with p T >0.4 GeV, passing the track selections discussed in Section 3.1. From the simulated events (Section 3.1), it is determined that the tracking efficiency reduces the measured N rec ch relative to the event generator multiplicity for p T >0.4 GeV primary charged particles 6 by approximately multiplicity-independent factors. The reduction factors and their uncertainties are 1.29 ± 0.05 and 1.18 ± 0.05 for the p+Pb and pp collisions, respectively.
For p+Pb collisions there is a direct correlation between N rec ch and the number of participating nucleons in the Pb nucleus: events with larger N rec ch values have, on average, a larger number of participating nucleons in the Pb nucleus and a smaller impact parameter. In this case, the concept of centrality used in A+A collisions is applicable, and in this paper the terms "central" and "peripheral" are used to refer to events with large and small values of N rec ch , respectively. For pp collisions there may not be a correlation between N rec ch and impact parameter. However, for convenience, the pp events with large and small N rec ch are also termed as "central" and "peripheral", respectively. Figure 1 shows the N rec ch distributions for the three data sets used in this paper. The discontinuities in the distributions result from the different HMT triggers, for which an offline requirement of N rec ch >N HLT trk is applied. This requirement ensures that the HMT triggered events are used only where the HLT trigger is almost fully efficient. The pp event activity can also be quantified using the total transverse energy deposited in the FCal (E FCal T ). This quantity has been used to determine the centrality in all ATLAS heavy-ion analyses. Using the E FCal T to characterize the event activity has the advantage that independent sets of particles are used to determine the event activity and to measure the long-range correlations. Similarly in the p+Pb case, the event activity can be characterized by the sum of transverse energy measured on the Pb-fragmentation side of the FCal (E FCal,Pb T ) [2,4]. Results presented in this paper use both N rec ch and the E FCal T (or E FCal,Pb T ) to quantify the event activity.

Two-particle correlation analysis
The study of two-particle correlations in this paper follows previous ATLAS measurements in Pb+Pb [13,69,70], p+Pb [2,4] and pp [41] collisions. For a given event class, the two-particle correlations are measured as a function of the relative azimuthal angle ∆φ ≡ φ a − φ b and pseudorapidity ∆η ≡ η a − η b separation. The labels a and b denote the two particles in the pair, which may be selected from different p T intervals. The particles a and b are conventionally referred to as the "trigger" and "associated" particles, respectively. The correlation function is defined as: where S and B represent pair distributions constructed from the same event and from "mixed events" [71], respectively. The same-event distribution S is constructed using all particle pairs that can be formed in each event from tracks that have passed the selections described in Section 3.1. The S distribution contains both the physical correlations between particle pairs and correlations arising from detector acceptance effects. The mixed-event distribution B(∆η, ∆φ) is similarly constructed by choosing the two particles in the pair from different events. The B distribution does not contain physical correlations, but has detector acceptance effects similar to those in S . In taking the ratio, S /B in Eq. (5), the detector acceptance effects largely cancel, and the resulting C(∆η, ∆φ) contains physical correlations only. The pair of events used in the mixing are required to have similar N rec ch (|∆N rec ch |<10) and similar z vtx (|∆z vtx |<10 mm), so that acceptance effects in S(∆η, ∆φ) are properly reflected in, and compensated by, corresponding variations in B(∆η, ∆φ). To correct S(∆η, ∆φ) and B(∆η, ∆φ) for the individual φ-averaged inefficiencies of particles a and b, the pairs are weighted by the inverse product of their tracking efficiencies 1/( a b ). Statistical uncertainties are calculated for C(∆η, ∆φ) using standard error-propagation procedures assuming no correlation between S and B, and with the statistical variance of S and B in each ∆η and ∆φ bin taken to be 1/( a b ) 2 where the sum runs over all of the pairs included in the bin. Typically, the two-particle correlations are used only to study the shape of the correlations in ∆φ, and are conveniently normalized. In this paper, the normalization of C(∆η, ∆φ) is chosen such that the ∆φ-averaged value of C(∆η, ∆φ) is unity for |∆η| > 2.
Examples of correlation functions are shown in Figure 2 for 0.5<p a,b T <5 GeV and for two different N rec ch ranges for each of the three data sets: 13 TeV pp (top), 5.02 TeV pp (middle), and 5.02 TeV p+Pb (bottom). The left panels show results for 0≤N rec ch <20 while the right panels show representative highmultiplicity ranges of N rec ch ≥120 for the 13 TeV pp data, 90≤N rec ch <100 for the 5.02 TeV pp data and N rec ch ≥220 for the 5.02 TeV p+Pb data. The correlation functions are plotted over the range −π/2<∆φ<3π/2; the periodicity of the measurement requires that C(∆η, 3π/2)=C(∆η, −π/2). The low-multiplicity correlation functions exhibit features that are understood to result primarily from hard-scattering processes: a peak centered at ∆η=∆φ=0 that arises primarily from jets and an enhancement centered at ∆φ=π and extending over the full ∆η range which results from dijets. These features also dominate the highmultiplicity correlation functions. Additionally, in the high-multiplicity correlation functions, each of the three systems exhibit a ridge -an enhancement centered at ∆φ=0 that extends over the entire measured ∆η range.
One-dimensional correlation functions C(∆φ) are obtained by integrating the numerator and denominator of Eq. (5) over 2<|∆η|<5 prior to taking the ratio: This |∆η| range is chosen to focus on the long-range features of the correlation functions. From the one-dimensional correlation functions, "per-trigger-particle yields," Y(∆φ) are calculated [2,4,71]: where N a denotes the total number of trigger particles, corrected to account for the tracking efficiency. The Y(∆φ) distribution is identical in shape to C(∆φ), but has a physically relevant normalization: it represents the average number of associated particles per trigger particle in a given ∆φ interval. This allows operations, such as subtraction of the Y(∆φ) distribution in one event-activity class from the Y(∆φ) distribution in another, which have been used in studying the p+Pb ridge [2,4].

Template fitting
In order to separate the ridge from other sources of angular correlation, such as dijets, the ATLAS Collaboration developed a template fitting procedure described in Ref. [41]. In this procedure, the measured Y(∆φ) distributions are assumed to result from a superposition of a "peripheral" Y(∆φ) distribution, Y periph (∆φ), scaled up by a multiplicative factor and a constant modulated by cos(n∆φ) for n ≥2. The resulting template fit function,  Figure 2: Two-particle correlation functions C(∆η, ∆φ) in 13 TeV pp collisions (top panels), 5.02 TeV pp collisions (middle panels) and in 5.02 TeV p+Pb collisions (bottom panels). The left panels correspond to a lowermultiplicity range of 0≤N rec ch <20. The right panels correspond to higher multiplicity ranges of N rec ch ≥120 for 13 TeV pp, 90≤N rec ch <100 for the 5.02 TeV pp and N rec ch ≥220 for the 5.02 TeV p+Pb. The plots are for charged particles having 0.5<p a,b T <5 GeV. The distributions have been truncated to suppress the peak at ∆η=∆φ=0 and are plotted over |∆η|<4.6 (|∆η|<4.0 for middle row) to avoid statistical fluctuations at larger |∆η|. For the middle-right panel, the peak at ∆φ=π has also been truncated. where has free parameters F and v n,n . The parameter F is the multiplicative factor by which the Y periph (∆φ) is scaled. The coefficient G, which represents the magnitude of the combinatoric component of Y ridge (∆φ), is fixed by requiring that the integral of Y templ (∆φ) be equal to the integral of the measured Y(∆φ): In this paper, when studying the N rec ch dependence of the long-range correlation, the 0≤N rec ch <20 multiplicity interval is used to produce Y periph (∆φ). When studying the E FCal The template fitting procedure is similar to the peripheral subtraction procedure used in previous ATLAS p+Pb ridge analyses [4]. In those analyses, the scale factor for the peripheral reference, analogous to F in Eq. (8), was determined by matching the near-side jet peaks between the peripheral and central samples. A more important difference, however, lies in the treatment of the peripheral bin. In the earlier analyses, a ZYAM procedure was performed on the peripheral reference, and only the modulated part of Y periph (∆φ), Y periph (∆φ) − Y periph (0), was used in the peripheral subtraction. 7 The ZYAM procedure makes several assumptions, the most relevant of which for the present analysis is that there is no longrange correlation in the peripheral bin. As pointed out in Ref.
[41], neglecting the non-zero modulation present in Y periph (∆φ) significantly biases the measured v n,n values. Results from an alternative version of the template fitting, where a ZYAM procedure is performed on the peripheral reference, by using , are also presented in this paper. This ZYAMbased template fit is similar to the p+Pb peripheral subtraction procedure. These results are included mainly to compare with previous measurements and to demonstrate the improvements obtained using the present method.
In Ref.
[41] the template fitting procedure only included the second-order harmonic v 2,2 , but was able to reproduce the N rec ch -dependent evolution of Y(∆φ) on both the near and away sides. The left panel of Figure 3 shows such a template fit, in the 13 TeV pp data, that only includes v 2,2 . The right panel shows the difference between the Y(∆φ) and the Y templ (∆φ) distributions demonstrating the presence of small (compared to v 2,2 ), but significant residual v 3,3 and v 4,4 components. While it is possible that cos 3∆φ and cos 4∆φ contributions could arise in the template fitting method due to small multiplicity-dependent changes in the shape of the dijet component of the correlation function, such effects would not produce the excess at ∆φ∼0 observed in the right-hand panel in Figure 3. That excess and the fact that its magnitude is compatible with the remainder of the distribution indicates that there is real cos 3∆φ and cos 4∆φ modulation in the two-particle correlation functions. Thus this paper extends the v 2,2 results in Ref.
[41] by including v 3,3 and v 4,4 as well. A study of these higher-order harmonics, including their N rec ch and p T dependence and factorization (Eq. (4)), can help in better understanding the origin of the long-range correlations. Figure 4 shows template fits to the 13 TeV (left panels) and 5.02 TeV pp data (right panels), for 0.5 < p a,b T < 5 GeV. From top to bottom, each panel represents a different N rec ch range. The template fits (Eq. (9)) include harmonics 2-4. Visually, a ridge, i.e. a peak on the near-side, cannot be seen in the top two rows, which correspond to low and intermediate N rec ch intervals, respectively. However, the template fits indicate the presence of a large modulated component of Y ridge (∆φ) even in these N rec ch intervals. Several prior pp ridge measurements rely on the ZYAM method [71, 72] to extract yields on the near-side [14,15]. In these analyses, the yield of excess pairs in the ridge above the minimum of the Y(∆φ) distribution is considered to be the strength of the ridge. Figure 4 shows that such a procedure would give zero yields in low-and intermediate-multiplicity collisions where the minimum of Y(∆φ) occurs at ∆φ∼0. In high-multiplicity events the ZYAM-based yields, while non-zero, are still underestimated. Figure 5 shows the template fits to the p+Pb data in a format similar to Figure 4. The template fits describe the data well across the entire N rec ch range used in this paper. Previous p+Pb ridge analyses used a peripheral subtraction procedure to remove the jet component from Y(∆φ) [1]. That procedure is similar to the ZYAM-based template fitting procedure, in that it assumes absence of any long-range correlations in the peripheral events. In the following sections, comparisons between the v n,n obtained from these two methods are shown.     Figure 4: Template fits to the per-trigger particle yields Y(∆φ), in 13 TeV (left panels) and in 5.02 TeV (right panels) pp collisions for charged-particle pairs with 0.5<p a,b T <5 GeV and 2<|∆η|<5. The template fitting includes secondorder, third-order and fourth-order harmonics. From top to bottom, each panel represents a different N rec ch range. The solid points indicate the measured Y(∆φ), the open points and curves show different components of the template (see legend) that are shifted along the y-axis by G or by FY periph (0), where necessary, for presentation.    Figure 5: Template fits to the per-trigger particle yields Y(∆φ) in 5.02 TeV p+Pb collisions for charged-particle pairs with 0.5<p a,b T <5 GeV and 2<|∆η|<5. The template fitting includes second-order, third-order and fourth-order harmonics. Each panel represents a different N rec ch range. The solid points indicate the measured Y(∆φ), the open points and curves show different components of the template (see legend) that are shifted along the y-axis by G or by FY periph (0), where necessary, for presentation. Figure 6 shows the v n,n obtained from the template fits in the 13 TeV pp data, as a function of N rec ch and E FCal T . The v n,n from the ZYAM-based template fits as well as the coefficients obtained from a direct Fourier transform of Y(∆φ):

Fourier coefficients
are also shown for comparison. While the template-v n,n are the most physically meaningful quantities, the Fourier-v n,n are also included to demonstrate how the template fitting removes the hard contribution. Similarly, the ZYAM-based template-v n,n are also included, as the ZYAM-based fitting is similar to the peripheral subtraction procedure used in prior p+Pb analyses [2,4], and comparing with the ZYAM-based results illustrates the improvement brought about in the template fitting procedure.
The v 2,2 values are nearly independent of N rec ch throughout the measured range. As concluded in Ref.
[41], this implies that the long-range correlation is not unique to high-multiplicity events, but is in fact present even at very low multiplicities. In the E FCal T dependence, however, v 2,2 shows a systematic decrease at low E FCal T . Further, the asymptotic value of the template-v 2,2 at large N rec ch is also observed to be ∼10% larger than the asymptotic value at large E FCal T . This might indicate that the v 2,2 at a given rapidity is more correlated with the local multiplicity than the global multiplicity.
The removal of the hard-process contribution to v 2,2 in the template fitting can be seen by comparing to the Fourier-v 2,2 values. The Fourier-v 2,2 values are always larger than the template-v 2,2 and show a systematic increase at small N rec ch (E FCal T ). This indicates the presence of a relatively large contribution from backto-back dijets over this range. Asymptotically, at large N rec ch the Fourier-v 2,2 values become stable, but show a small decreasing trend in the E FCal T dependence. The ZYAM-based v 2,2 values are smaller than the template-v 2,2 values for all N rec ch (E FCal T ), and by construction systematically decrease to zero for the lower N rec ch (E FCal T ) intervals. However, at larger N rec ch (E FCal T ) they also show only a weak dependence on N rec ch (E FCal T ). Asymptotically, at large N rec ch the v 2,2 values from the Fourier transform and the default template fits match to within ∼10% (relative). In general, the v 2,2 values from all three methods agree within ±15% at large N rec ch or E FCal T . This implies that at very high multiplicities, N rec ch ∼ 120, the ridge signal is sufficiently strong that the assumptions made in removing the hard contributions to Y(∆φ) do not make a large difference. However, for the highest p T values used in this analysis, p a T >7 GeV, it is observed that the width of the dijet peak in the pp correlation functions broadens with increasing multiplicity. This change is opposite to that seen at lower p T where v 2,2 causes the dijet peak to become narrower. As a result, the measured v 2,2 values become negative. This bias from the multiplicity dependence of the hardscattering contribution likely affects the correlation functions at lower p a,b T values and its potential impact is discussed below.
The v 2,2 component is dominant, with a magnitude approximately 30 times larger than v 3,3 and v 4,4 , which are comparable to each other. This is in stark contrast to Pb+Pb collisions where in the most central events, where the average geometry has less influence, the v n,n have comparable magnitudes [13]. The Fourier-v 3,3 shows considerable N rec ch (E FCal T ) dependence and is negative almost everywhere. However, the v 3,3 values from the template fits are mostly positive. As the factorization of the v n,n requires that the v n,n be positive (Eq. (3)), the negative Fourier-v 3,3 clearly does not arise from single-particle modulation. However, because the template-v 3,3 is positive, its origin from single-particle modulation cannot be ruled out. Within statistical uncertainties, the v 4,4 values from all three methods are positive throughout the   larger compared to the systematic uncertainties in the v 2,2 values (discussed later in Section 6). This is possibly indicative of a systematic change in the average collision geometry which is present in p+Pb but not in pp collisions. A similar increase of the v 2,2 values is also observed in the E FCal,Pb T dependence. The higher-order harmonics v 3,3 and v 4,4 show a stronger relative increase with increasing N rec ch and E FCal,Pb T . This also implies that the assumption made in the template-fitting, regarding the independence or weak dependence of the v n,n on N rec ch , is not strictly correct for v 3,3 and v 4,4 . Figure 8 also compares the Fourier and ZYAM-based template-v n,n values. The v n,n from the peripheral subtraction procedure used in a previous ATLAS p+Pb long-range correlation analysis [4] are also shown. The peripheral-subtracted v n,n values are nearly identical to the values obtained from the ZYAM-based template fits. This is expected, as the treatment of the peripheral bin is identical in both cases: both use the ZYAM-subtracted Y periph (∆φ) as the peripheral reference. What differs procedurally between the two methods is determination of the scale factor by which Y periph (∆φ) is scaled up when subtracting it from Y(∆φ). In the peripheral subtraction case, this scale factor, analogous to the parameter F in Eq. (8), is determined by matching the near-side jet peaks over the region |∆η|<1 and |∆φ|<1. In the template-fitting case, the parameter F is determined by the jet contribution to the away-side peak. The similarity of the v 2,2 values from the two procedures implies that whether the matching is done in the near-side jet peak, or over the away-side peak, identical values of the scale factor are obtained. The Fourier-v 2,2 and . This is unlike the pp case (Figures 6 and 7), where the values differed by ∼15% (relative) at large N rec ch . This similarity does not hold for v 3,3 where the values from the template fit are systematically larger than the values obtained from Fourier decomposition. For all harmonics, the relative difference in the v n,n decreases with increasing event activity. Like in the pp case (Figure 6), this implies that at large enough event activity, the v n,n are less sensitive to the assumptions made in removing the hard contributions.

Test of factorization in template fits
If the v n,n obtained from the template fits are the result of single-particle modulations, then the v n,n should factorize as in Eq. (3), and the v n (p a T ) obtained by correlating trigger particles at a given p a T with associated particles in several different intervals of p b T (Eq. (4)) should be independent of the choice of the p b T interval. Figure 9 demonstrates the factorization of the v 2,2 in the 13 TeV pp data, as a function of N rec ch . The left panel shows the v 2,2 values for 0.5<p a T <5 GeV and for four different choices of the associated particle p T : 0.5-5, 0.5-1, 1-2 and 2-3 GeV. The right panel shows the corresponding v 2 (p a T ) obtained using Eq. (4). While the v 2,2 (p a T , p b T ) values vary by a factor of ∼2 between the different choices of the p b T interval, the corresponding v 2 (p a T ) values agree quite well. Similar plots for the p+Pb data are shown in Figure 10. Here due to higher statistical precision in the data, the factorization is tested for both v 2,2 and v 3,3 . The variation of v 2,2 (p a T , p b T ) between the four p b T intervals is a factor of ∼2 while the variation of v 3,3 (p a T , p b T ) is more than a factor of 3. However, the corresponding v n (p a T ) values are in good agreement with each other, with the only exception being the v 2,2 values for 2<p b T <3 GeV where some deviation from this behavior is seen for N rec ch 60. Figure 11 studies the p a T dependence of the factorization in the 13 TeV pp data for v 2,2 (top panels) and v 3,3 (bottom panels). The results are shown for the N rec ch ≥90 multiplicity range. The left panels show the v n,n as a function of p a T for four different choices of the associated particle p T : 0.5-5, 0.5-1, 1-2 and 2-3 GeV. The right panels show the corresponding v n (p a T ) obtained using Eq. (4). In the v 2,2 case, factorization holds reasonably well for p a T ≤3 GeV, and becomes worse at higher p T . This breakdown at higher p T is likely caused by the above-discussed multiplicity-dependent distortions of the dijet component of the correlation function which are not accounted for in the template fitting procedure. For v 3,3 , the factorization holds reasonably well for p b T >1 GeV. The 0.5<p b T <1 GeV case shows a larger deviation in the factorization, but has much larger associated statistical uncertainties. Similar plots for the p+Pb case are shown in Figure 12

Dependence of u n,n on ∆η gap
A systematic study of the ∆η dependence of the v n,n can also help in determining the origin of the longrange correlation. If it arises from mechanisms that only correlate a few particles in an event, such as jets, then a strong dependence of the correlation on the ∆η gap between particle pairs is expected. Figure 13 shows the measured v n,n (left panels) and v n = √ v n,n (right panels), as a function of |∆η| for |∆η|>1 in the 13 TeV pp data. Also shown for comparison are the Fourier and ZYAM-based template-v n,n . The template-v 2,2 (top left panel) and v 2 (top right panel) are quite stable, especially for |∆η|>1.5, where the influence of the near-side jet is diminished. In contrast, the Fourier-v 2,2 show a strong |∆η| dependence. The ∆η dependence is largest at small |∆η| because of the presence of the sharply peaked near-side jet, but is considerable even for |∆η|>2. Similarly, the Fourier-v 3,3 shows large |∆η| dependence, going from positive values at |∆η|∼1 to negative values at large |∆η|, while the template-v 3,3 change only weakly in comparison. The Fourier-v 3,3 is often negative, ruling out the possibility of it being generated by singleparticle anisotropies, which require that v n,n = v 2 n be positive. For points where v 3,3 is negative, v 3 is not defined and hence not plotted. The template-v 3,3 is, however, positive and, therefore, consistent with a single-particle anisotropy as its origin, except for the highest |∆η| interval where it is consistent with zero. The v 4,4 values, like the v 2,2 and v 3,3 values, vary only weakly with |∆η|. These observations further support the conclusion that the template-v n,n are coefficients of genuine long-range correlations.    Figure 13: The |∆η| dependence of the v n,n (left panels) and v n (right panels) in the 13 TeV pp data. From top to bottom the rows correspond to n=2, 3 and 4, respectively. The ZYAM-template and Fourier-v n,n values are also shown for comparison. Only the range |∆η|>1 is shown to suppress the large Fourier-v n,n at |∆η|∼0 that arise due to the near-side jet peak. Plots are for the N rec ch ≥90 multiplicity range and for 0.5<p a,b T <5 GeV. The error bars indicate statistical uncertainties only. For points where v 3,3 is negative, v 3 is not defined and hence not plotted.

Systematic uncertainties and cross-checks
The systematic uncertainties in this analysis arise from choosing the peripheral bin used in the template fits, pileup, tracking efficiency, pair-acceptance and Monte Carlo consistency. Each source is discussed separately below.
Peripheral interval: As explained in Section 5, the template fitting procedure makes two assumptions. First it assumes that the contributions to Y(∆φ) from hard processes have identical shape across all event activity ranges, and only change in overall scale. Second, it assumes that the v n,n are only weakly dependent on the event activity. The assumptions are self-consistent for the N rec ch dependence of the v n,n in the 5.02 and 13 TeV pp data (Figures 6-7), where the measured template-v n,n values do turn out to be nearly independent of N rec ch . However, for the E FCal T dependence in the pp data, and for both the N rec ch and E FCal,Pb T dependence in the p+Pb data, a systematic increase of the template-v 2,2 with event activity is seen at small event activity. This indicates the breakdown of one of the above two assumptions. To test the sensitivity of the measured v n,n to any residual changes in the width of the away-side jet peak and to the v n,n present in the peripheral reference, the analysis is repeated using 0≤N rec ch <10 and 10≤N rec ch <20 intervals to form Y periph (∆φ). The variations in the v n,n for the different chosen peripheral intervals are taken to be a systematic uncertainty. For a given dataset, this uncertainty is strongly correlated across all multiplicity intervals. Choosing a peripheral interval with larger mean multiplicity typically decreases the measured v n,n .
The sensitivity of the template-v 2 to which peripheral interval is chosen is demonstrated in the left panels of Figure 14,  <20. In both the 13 TeV and 5.02 TeV pp data, except at very low N rec ch , the v 2 values are nearly independent of the chosen peripheral reference. In the 13 TeV pp case, the variation is ∼6% at N rec ch ∼30 and decreases to ∼1% for N rec ch ≥60. Even in the p+Pb case, where the measured template-v 2,2 exhibits some dependence on N rec ch , the dependence of the template-v 2 on the choice of peripheral bin is quite small: ∼6% at N rec ch ∼30 and decreases to ∼2% for N rec ch ∼60. Also shown for comparison are the corresponding v 2 values obtained from the ZYAM-based template fitting method (right panels of Figure 14). These exhibit considerable dependence on the peripheral reference. For the 13 TeV pp case, the variation in the ZYAM-based v 2 is ∼40% at N rec ch ∼30, and decreases to ∼12% at N rec ch ∼60 and asymptotically at large N rec ch is ∼7%. For the p+Pb case, the variation is even larger: ∼35% at N rec ch ∼30 and ∼14% for N rec ch ∼60. These results show that the template-v 2 is quite stable as the peripheral interval is varied, while the ZYAM-based result is very sensitive. This is one of the advantages of the new method. For the ZYAM-based results, as the upper edge of the peripheral interval is moved to lower multiplicities, the measured v 2 becomes less and less dependent on N rec ch . Qualitatively, it seems that in the limit of N rec,periph ch → 0 the ZYAM-based pp-v 2 would be nearly independent of N rec ch , thus contradicting the assumption of zero v 2 made in the ZYAM method, and supporting the flat-v 2 assumption made in the new method.
Pileup: Pileup events, when included in the two-particle correlation measurement, dilute the v n,n signal since they produce pairs where the trigger and associated particle are from different collisions and thus have no physical correlations. The maximal fractional dilution in the v n,n is equal to the pileup rate. In the p+Pb data, nearly all of the events containing pileup are removed by the procedure described in Section 3. The influence of the residual pileup is evaluated by relaxing the pileup rejection criteria and then calculating the change in the Y(∆φ) and v n values. The differences are taken as an estimate of the  uncertainty for the v n,n , and are found to be negligible in low event activity classes, and increase to 4% for events with N rec ch ∼300. In the pp data, for events containing multiple vertices, only tracks associated with the vertex having the largest p 2 T , where the sum is over all tracks associated with the vertex, are used in the analysis. Events with multiple unresolved vertices affect the results by increasing the combinatoric pedestal in Y(∆φ). The fraction of events with merged vertices is estimated and taken as the relative uncertainty associated with pileup in the pp analysis. The merged-vertex rate in the 13 TeV pp data is 0-3% over the 0-150 N rec ch range. In the 5.02 TeV pp data, it is 0-4% over the 0-120 N rec ch range. Track reconstruction efficiency: In evaluating Y(∆φ), each particle is weighted by 1/ (p T , η) to account for the tracking efficiency. The systematic uncertainties in the efficiency (p T , η) thus need to be propagated into Y(∆φ) and the final v n,n measurements. Unlike Y(∆φ), which is strongly affected by the efficiency, the v n,n are mostly insensitive to the tracking efficiency. This is because the v n,n measure the relative variation of the yields in ∆φ; an overall increase or decrease in the efficiency changes the yields but does not affect the v n,n . However, as the tracking efficiency and its uncertainties have p T and η dependence, there is some residual effect on the v n,n . The corresponding uncertainty in the v n,n is estimated by repeating the analysis while varying the efficiency to its upper and lower extremes. In the pp analysis, this uncertainty is estimated to be 0.5% for v 2,2 and 2.5% for v 3,3 and v 4,4 . The corresponding uncertainties in the p + Pb data are 0.8%, 1.6% and 2.4% for v 2,2 , v 3,3 and v 4,4 , respectively.
Pair-acceptance: As described in Section 4, this analysis uses the mixed-event distributions B(∆η, ∆φ) and B(∆φ) to estimate and correct for the pair-acceptance of the detector. The mixed-event distributions are in general quite flat in ∆φ. The Fourier coefficients of the mixed-event distributions, v det n,n , which quantify the magnitude of the corrections, are ∼10 −4 in the p+Pb data, and ∼2 × 10 −5 in the pp data. In the p+Pb analysis, potential systematic uncertainties in the v n,n due to residual pair-acceptance effects not corrected by the mixed-events are evaluated following Ref. [13]. This uncertainty is found to be smaller than ∼10 −5 . In the pp analysis, since the mixed-event corrections are themselves quite small, the entire correction is conservatively taken as the systematic uncertainty.
MC closure: The analysis procedure is validated by measuring the v n,n of reconstructed particles in fully simulated Pythia 8 and HIJING events and comparing them to those obtained using the generated particles. The difference between the generated and reconstructed v n,n varies between 10 −5 and 10 −4 (absolute) in the pp case and between 2% and 8% (relative) in the p+Pb case, for the different harmonics. This difference is an estimate of possible systematic effects that are not accounted for in the measurement, such as a mismatch between the true and reconstructed momentum for charged particles, and is included as a systematic uncertainty.
As a cross-check, the dependence of the long-range correlations on the relative charge of the two particles used in the correlation is studied. If the long-range correlations arise from phenomena that correlate only a few particles in an event, such as jets or decays, then a dependence of the correlation on the relative sign of the particles making up the pair is expected. Figure 15 shows the measured v 2 from the template fits for both the same-charge and opposite-charge pairs. No systematic difference between the two is observed. Tables 2 and 3 list the systematic uncertainties in the v n,n for the 13 TeV and 5.02 TeV pp data, respectively. Most uncertainties are listed as relative uncertainties (in percentages of the v n,n ), while some are listed as absolute uncertainties. Uncertainties for the p+Pb data are listed in Table 4. The corresponding uncertainties in the v n are obtained by propagating the uncertainties in the v n,n when using Eq. (3) to obtain the v n . In some cases the systematic uncertainties in the v n,n are larger than 100%. In these cases 0.5 2.5 2.5 Pair acceptance (absolute) 2 × 10 −5 2 × 10 −5 2 × 10 −5 MC closure (absolute) 1 × 10 −4 2 × 10 −5 2 × 10 −5 Table 2: Systematic uncertainties for the v n,n obtained from the template analysis in the 13 TeV pp data. Where ranges are provided for both multiplicity and the uncertainty, the uncertainty varies from the first value to the second value as the multiplicity varies from the lower to upper limits of the range. Where no multiplicity range is provided the uncertainty is multiplicity-independent.
Source   Table 4: Systematic uncertainties for the v n,n obtained from the template analysis in the 5.02 TeV p+Pb data. Where ranges are provided for both multiplicity and the uncertainty, the uncertainty varies from the first value to the second value as the multiplicity varies from the lower to upper limits of the range. Where no multiplicity range is provided the uncertainty is multiplicity-independent.

Results
Figure 16 provides a summary of the main results of this paper in the inclusive p T interval 0.5<p T <5 GeV. It compares the v n obtained from the 5.02 TeV, 13 TeV pp and 5.02 TeV p+Pb template fits. The left panels show v 2 , v 3 and v 4 as a function of N rec ch while the right panels show the results as a function of p a T for the N rec ch ≥60 multiplicity range. The measured v 3 and v 4 in the 5.02 TeV pp data for 0.5<p a,b T <5 GeV have large systematic uncertainties associated with the choice of peripheral reference and are not shown in Figure 16. They are shown in Figure 18 for a different p T interval of 1<p a,b T <5 GeV. Figure 16 shows that the p+Pb v 2 increases with increasing N rec ch as previously observed [4] while the pp v 2 is N rec ch -independent within uncertainties. The p+Pb v 3 is significantly larger than the pp v 3 and also shows a systematic increase with N rec ch , while the pp v 3 is consistent with being N rec ch -independent. The pp and p+Pb v 4 are consistent within large uncertainties, and the p+Pb v 4 increases weakly with increasing N rec ch . The difference between the pp and p+Pb results for the N rec ch dependence of the v n is expected. Studies of the centrality dependence of the multiplicity distributions in p+Pb collisions show a strong correlation between the multiplicity and the number of participants, or equivalently, the number of scatterings of the proton in the nucleus [73]. Regardless of the interpretation of the results, a dependence of the v n on the geometry of the p+Pb collisions is expected [74]. In contrast, the relationship between multiplicity and geometry in pp collisions is poorly understood and necessarily different as there are, by definition, only two colliding nucleons. However, an early study of this problem accounting for perturbative evolution did predict a weak dependence of v 2 on multiplicity, as observed in this measurement [75]. A more recent study that models the proton substructure and fluctuations in the multiplicity of the final particles, showed that the eccentricities 2 and 3 of the initial entropy-density distributions in pp collisions have no correlation with the final particle multiplicity [76]. If the v n in pp collisions are directly related to the n , then the calculations in Ref.
[76] are consistent with the trends observed in the measured v n .
The pp and p+Pb v 2 (p T ) shown in Figure 16 display similar trends with both increasing with p T at low p T , reaching a maximum near 3 GeV and decreasing at higher p T . The v 2 (p T ) values for the 5.02 and 13 TeV pp data agree within uncertainties. The p T dependence of the v 3 and v 4 values is similar to that of v 2 at low p T , where the p+Pb results increase more rapidly with increasing p T . However, unlike for v 2 , the values of v 3 and v 4 are similar at high p T for the pp and p+Pb data. A direct test of the similarity of the p T dependence of the Fourier coefficients in pp and p+Pb collisions is provided in Figure 17 for n = 2. The pp v 2 values have been multiplied by 1.51, the ratio (p+Pb to pp) of the maximum v 2 in the top right panel in Figure 16. The resulting v 2 (p a T ) values for (scaled) pp and p+Pb agree well for p a T up to 5 GeV. At higher p a T the pp v 2 decreases more rapidly due to the above-described multiplicity-dependent change in the shape of the dijet peak in the two-particle correlation function at high p T . After the scaling, the pp v 2 (p a T ) are slightly higher than the p+Pb at low p a T , but the similarity of the shapes of the p T dependence is, nonetheless, striking.
A separate evaluation of the N rec ch -dependence of the v 2 , v 3 and v 4 values is shown in Figure 18 for the 1<p a,b T <5 GeV interval, where the 5.02 TeV pp measurements yield meaningful v 3 and v 4 results. The figure shows agreement between the 5.02 and 13 TeV pp data for all three Fourier coefficients. It also shows that the p+Pb v 2 , v 3 and v 4 rise monotonically with increasing N rec ch while the pp results are generally N rec ch -independent. One possible exception to this statement is that the 13 TeV data indicate a small (∼15%) decrease in v 2 in the two lowest N rec ch intervals. The pp and p+Pb v 3 and v 4 agree at low N rec ch while v 2 still differs significantly, although by a smaller amount than at larger N rec ch . This behavior is different from that observed in the inclusive p T interval, which may, in turn, reflect the convergence of the v 2 (p T ) between the pp and p+Pb data shown in the top, right panel of Figure 16.
Measurements [70,77] and theoretical analyses [78][79][80][81][82] of the correlations between the Fourier coefficients and event-plane angles of different flow harmonics in Pb+Pb collisions have indicated significant "non-linearity" resulting from collective expansion such that the response of the medium to an initial elliptic eccentricity can contribute to cos (4φ) modulation of the produced particles. In Pb+Pb collisions, the non-linear contribution to v 4 is found to dominate over the geometric contribution except for the most central collisions where the initial-state fluctuations have the greatest impact. The non-linear contribution to v 4 is expected to be proportional to v 2 2 so a comparison of the measured v 4 to v 2 2 in pp and p+Pb collisions may be of interest. The results are presented in Figure 19, which shows v 4 /v 2 2 versus N rec ch for the 13 TeV pp and the p+Pb data. In the ratio, the correlated systematic uncertainties between the measured v 4 and v 2 2 cancel. The ratio is observed to be constant as a function of N rec ch for both data sets even  Figure 16: Left panels: comparison of the v n obtained from the template fitting procedure in the 13 TeV pp, 5.02 TeV pp, and 5.02 TeV p+Pb data, as a function of N rec ch . The results are for 0.5<p a,b T <5 GeV. Right panels: the p T dependence of the v n for the N rec ch ≥60 multiplicity range. From top to bottom the rows correspond to n=2, 3 and 4, respectively. The error bars and shaded bands indicate statistical and systematic uncertainties, respectively. though the p+Pb v 2 and v 4 increase with N rec ch . The v 4 /v 2 2 ratio is observed to be 50% larger in the pp data than in the p+Pb data. Naively, this would indicate a larger non-linear contribution to v 4 in pp collisions compared to p+Pb collisions.

Conclusion
In summary, this paper presents results of two-charged-particle correlation measurements made by AT-LAS in √ s = 13 and 5.02 TeV pp collisions and in 5.02 TeV p+Pb collisions at the LHC. This measurement uses integrated luminosities of 64 nb −1 for the √ s=13 TeV pp data, 170 nb −1 for the √ s=5.02 TeV pp data and 28 nb −1 for the √ s NN =5.02 TeV p+Pb data. The 13 TeV measurements represent an extension of results presented in Ref.
[41] using a larger data sample. The p+Pb results are obtained from a reanalysis of Run 1 data presented in Ref.
[4] using a template fitting procedure developed for pp collisions and applied in Ref.
[41]. The correlation functions are measured for different intervals of measured chargedparticle multiplicity and FCal transverse energy and for different intervals of charged-particle transverse momentum; many of the results are presented for an "inclusive" p T interval 0.5 < p T < 5 GeV.
One-dimensional distributions of per-trigger-particle yields as a function of azimuthal angle separation, Y(∆φ), are obtained from the long-range (|∆η| > 2) component of the correlation functions. A template fitting procedure is applied to the Y(∆φ) distributions to remove the contributions from hard-scattering processes and to measure the relative amplitudes v n,n of the sinusoidal modulation of the soft underlying event. Results for v 2,2 , v 3,3 , and v 4,4 are obtained for all three colliding systems. An analysis of the factorizability of the v n,n shows good factorization for most of the measured N rec ch and p T intervals although factorization is observed to break down for the most extreme combinations of p a T and p b T in the lowest and highest multiplicity or transverse energy intervals. Since the v n,n results are observed to be consistent with the presence of single-particle modulation of the per-event dN/dφ distributions, single-particle v n values are extracted and plotted versus N rec ch and p T . Comparisons of the v 2 , v 3 and v 4 values between 13 and 5.02 TeV pp collisions show no significant variation in these quantities with center-of-mass energy. As observed in Ref.
[41], the v 2 values obtained in pp collisions at both energies are observed to be independent of N rec ch within uncertainties for the inclusive p T interval. However, for the 1<p T <5 GeV interval a ∼15% decrease in v 2 is seen in the lowest N rec ch intervals. The p+Pb v 2 values are larger than the pp v 2 values for all multiplicities and are observed to increase slowly with N rec ch . However, the p+Pb trend appears to converge with the pp values for the lowest multiplicities, at least in the inclusive p T interval. For the 1<p T <5 GeV interval, the v 2 (p T ) trends do not show the same convergence between pp and p+Pb results. Similar to the results for v 2 , the pp v 3 and v 4 values are consistent with being independent of N rec ch within uncertainties and the p+Pb values are observed to increase with N rec ch . The pp and p+Pb v 3 and v 4 values are consistent within uncertainties in the lowest measured N rec ch intervals. The p T dependence of the pp and p+Pb v 2 values is similar: both rise approximately linearly with p T and reach a maximum near 3 GeV. The maximum p+Pb v 2 value is approximately 50% larger than the maximum v 2 values for the 13 and 5.02 TeV pp data, which are consistent within uncertainties. The p+Pb v 3 and v 4 values also increase more rapidly with increasing p T than the corresponding pp values for p T < 2 GeV, but the p+Pb v 3 values saturate above 3 GeV while the measured 13 TeV pp v 3 values continue to increase with increasing p T over the full range of the measurement. A test of the similarity of the p T dependence of the pp and p+Pb v 2 values rescaling pp v 2 values shows that the pp and p+Pb v 2 (p a T ) distributions are remarkably similar in shape for p a T <5 GeV. An evaluation of the v 4 /v 2 2 ratio in the inclusive p T interval shows results that are N rec ch -independent for both the 13 TeV pp data and the p+Pb data. This ratio is observed to be 50% larger for the pp data than for the p+Pb data.
The similarities between the pp and p+Pb results presented here suggest a common physical origin for the azimuthal anisotropies. The difference in the observed multiplicity dependence of the Fourier coefficients likely arises from the different geometry of the pp and p+Pb collisions.
[3] CMS Collaboration, Multiplicity and transverse momentum dependence of two-and four-particle correlations in pPb and PbPb collisions, Phys. [13] ATLAS Collaboration, Measurement of the azimuthal anisotropy for charged particle production in √ s NN = 2.76 TeV lead-lead collisions with the ATLAS detector,