Measurement of long-range pseudorapidity correlations and azimuthal harmonics in √sNN = 5.02 TeV proton-lead collisions with the ATLAS detector

Measurements of two-particle correlation functions and the ﬁrst ﬁve azimuthal harmonics, v 1 to v 5 , are presented, using 28 nb − 1 of p + Pb collisions at a nucleon-nucleon center-of-mass energy of √ s NN = 5 . 02 TeV measured with the ATLAS detector at the LHC. Signiﬁcant long-range “ridgelike” correlations are observed for pairs with small relative azimuthal angle ( | (cid:2)φ | < π/ 3) and back-to-back pairs ( | (cid:2)φ | > 2 π/ 3) over the transverse momentum range 0 . 4 < p T < 12 GeV and in different intervals of event activity. The event activity is deﬁned by either the number of reconstructed tracks or the total transverse energy on the Pb-fragmentation side. The azimuthal structure of such long-range correlations is Fourier decomposed to obtain the harmonics v n as a function of p T and event activity. The extracted v n values for n = 2 to 5 decrease with n . The v 2 and v 3 values are found to be positive in the measured p T range. The v 1 is also measured as a function of p T and is observed to change sign around p T ≈ 1 . 5–2.0 GeV and then increase to about 0.1 for p T > 4 GeV. The v 2 ( p T ), v 3 ( p T ), and v 4 ( p T ) are compared to the v n coefﬁcients in Pb + Pb collisions at √ s NN = 2 . 76 TeV with similar event multiplicities. Reasonable agreement is observed after accounting for the difference in the average p T of particles produced in the two collision systems.


I. INTRODUCTION
One striking observation in high-energy nucleus-nucleus (A + A) collisions is the large anisotropy of particle production in the azimuthal angle φ [1,2]. This anisotropy is often studied via a two-particle correlation of particle pairs in relative pseudorapidity ( η) and azimuthal angle ( φ) [3,4]. The anisotropy manifests itself as a strong excess of pairs at φ ∼ 0 and π , and the magnitude of the excess is relatively constant out to large | η| [5][6][7][8][9]. The azimuthal structure of this "ridgelike" correlation is commonly characterized by its Fourier harmonics, dN pairs /d φ ∼ 1 + 2 n v 2 n cos n φ. While the elliptic flow, v 2 , and triangular flow, v 3 , are the dominant harmonics in A + A collisions, significant v 1 , v 4 , v 5 , and v 6 harmonics have also been measured [8][9][10][11][12][13]. These v n values are commonly interpreted as the collective hydrodynamic response of the created matter to the collision geometry and its fluctuations in the initial state [14]. The success of hydrodynamic models in describing the anisotropy of particle production in heavy-ion collisions at BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC) places important constraints on the transport properties of the produced matter [15][16][17][18][19][20].
For a small collision system, such as proton-proton (p + p) or proton-nucleus (p + A) collisions, it was assumed that the transverse size of the produced system is too small for the hydrodynamic flow description to be applicable. Thus, it came * Full author list given at the end of the article. as a surprise that ridge-like structures were also observed in two-particle correlations in high-multiplicity p + p [21] and proton-lead (p + Pb) [22][23][24] collisions at the LHC and later in deuteron-gold collisions [25] at RHIC. A Fourier decomposition technique has been employed to study the azimuthal distribution of the ridge in p + Pb collisions. The transverse momentum (p T ) dependence of the extracted v 2 and v 3 [23,24], and the particle-mass dependence of v 2 [26] are found to be similar to those measured in A + A collisions.
The interpretation of the long-range correlations in highmultiplicity p + p and p + Pb collisions is currently a subject of intense study. References [30][31][32][33] argue that the produced matter in these collisions is sufficiently voluminous and dense that the hydrodynamic model framework may still apply. However, models based on gluon saturation and color connections suggest that the long-range correlations are an initialstate effect, intrinsic to QCD at high gluon density [34][35][36][37][38].
Recently, a hybrid model that takes into account both the initial-and the final-state effects has been proposed [39]. All these approaches can describe, qualitatively and even quantitatively, the v 2 and v 3 data in the p + Pb collisions.
To provide more insights into the nature of the ridge correlation and to discriminate between different theoretical interpretations, this paper provides a detailed measurement of the two-particle correlation and v n coefficients in p + Pb collisions at a nucleon-nucleon center-of-mass energy of √ s NN = 5.02 TeV. The data correspond to an integrated luminosity of approximately 28 nb −1 , recorded in 2013 by the ATLAS experiment at the LHC. This measurement benefits from a dedicated high-multiplicity trigger (see Sec. II B) that provides a large sample of high-multiplicity events, not only extending the previous v 2 and v 3 results to higher p T , each side, 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. A timing requirement based on signals from each side of the MBTS is imposed to suppress beam backgrounds. Owing to the high event rate during the run, only a small fraction of the MB events (∼1/1000) were recorded.
To enhance the number of events with high multiplicity, a dedicated high-multiplicity trigger (HMT) was implemented, which uses the ATLAS L1 and high-level trigger (HLT) systems [42]. At L1, the total transverse energy E L1 T in the FCal rapidity interval is required to be above a certain threshold. In the HLT, the charged-particle tracks are reconstructed by requiring at least two hits in the pixel detector and three hits in the SCT. For each event, the collision vertex reconstructed with the highest number of online tracks is selected, and the number of tracks (N HLT trk ) associated with this vertex with p T > 0.4 GeV and a distance of closest approach of less than 4 mm is calculated.
The HMT triggers are implemented by requiring different thresholds on the values of E L1 T and N HLT trk with prescale factors adjusted to the instantaneous luminosity provided by the LHC [42]. This analysis uses the six pairs of thresholds on E L1 T and N HLT trk listed in Table I. The N HLT trk 225 trigger was not prescaled throughout the entire running period.

A. Event and track selections
In the off-line analysis, p + Pb events are required to have a reconstructed vertex containing at least two associated offline tracks, with its z position satisfying |z vtx | < 150 mm. Noncollision backgrounds and photonuclear interactions are suppressed by requiring at least one hit in a MBTS counter on each side of the IP and the difference between times measured on the two sides to be less than 10 ns. In the 2013 p + Pb run, the luminosity conditions provided by the LHC result in an average probability of 3% that an event contains two or more p + Pb collisions (pileup). The pileup events are suppressed by rejecting events containing more than one good reconstructed vertex. The remaining pileup events are further suppressed based on the signal in the ZDC on the Pb-fragmentation side. This signal is calibrated to the number of detected neutrons (N n ) based on the location of the peak corresponding to a single neutron. The distribution of N n in events with pileup is broader than that for the events without pileup. Hence, a simple cut on the high tail end of the ZDC signal distribution further suppresses the pileup, while retaining more than 98% of the events without pileup. After this pileup rejection procedure, the residual pileup fraction is estimated to be 10 −2 in the event class with the highest track multiplicity studied in this analysis. About 57 × 10 6 MB-selected events and 15 × 10 6 HMT-selected events are included in this analysis. 044906-2 Charged-particle tracks are reconstructed in the ID using an algorithm optimized for p + p MB measurements [43]: The tracks are required to have p T > 0.3 GeV and |η| < 2.5, at least seven hits in the pixel detector and the SCT, and a hit in the first pixel layer when one is expected. 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 | < 3, respectively, where σ d 0 and σ z are uncertainties on d 0 and z 0 sin θ obtained from the track-fit covariance matrix.
The efficiency, (p T ,η), for track reconstruction and track selection cuts is obtained using p + Pb Monte Carlo events produced with version 1.38b of the HIJING event generator [44] with a center-of-mass boost matching the beam conditions. The response of the detector is simulated using GEANT4 [45,46] and the resulting events are reconstructed with the same algorithms as applied to the data. The efficiency increases with p T by 6% between 0.3 and 0.5 GeV and varies only weakly for p T > 0.5 GeV, where it ranges from 82% at η = 0 to 70% at |η| = 2 and 60% at |η| > 2.4. The efficiency is also found to vary by less than 2% over the multiplicity range used in the analysis. The extracted efficiency function (p T ,η) is used in the correlation analysis, as well as to estimate the average efficiencycorrected charged-particle multiplicity in the collisions.

B. Characterization of the event activity
The two-particle correlation (2PC) analyses are performed in event classes with different overall activity. The event activity is characterized by either E Pb T , the sum of transverse energy measured on the Pb-fragmentation side of the FCal with −4.9 < η < −3.2, or N rec ch , the off-line-reconstructed track multiplicity in the ID with |η| < 2.5 and p T > 0.4 GeV. These event-activity definitions have been used in previous p + Pb analyses [21,22,24,27,28]. Events with larger activity have on average a larger number of participating nucleons in the Pb nucleus and a smaller impact parameter. Hence, the term "centrality," familiar in A + A collisions, is used to refer to the event activity. The terms "central" and "peripheral" are used to refer to events with large activity and small activity, respectively.
Owing to the wide range of trigger thresholds and the prescale values required by the HMT triggers, the E Pb T and N rec ch distributions are very different for the HMT events and the MB events. To properly include the HMT events in the event-activity classification, an event-by-event weight, w = 1/P , is utilized. The combined probability, P , for a given event to be accepted by the MB trigger or any of the HMT triggers is calculated via the inclusion-exclusion principle as where N is the total number of triggers and p i is the probability for the ith trigger to accept the event, defined as zero if the event does not fire the trigger and otherwise as the inverse of the prescale factor of the trigger. The higher-order terms in Eq. (1) account for the probability of more than one trigger being fired. The weight factor, w, is calculated and applied event by event.
The distribution for all events after reweighting has the same shape as the distribution for MB events, as should be the case if the reweighting is done correctly. Figure 1 shows the distribution of N rec ch (left panels) and E Pb T (right panels) for the MB and MB + HMT events before (top panels) and after (bottom panels) the reweighting procedure. For MB-selected events, the reweighted distribution differs from the original distribution by a constant factor, reflecting the average prescale. The multiple steps in the N rec ch distribution (top-left panel) reflect the rapid turn-on behavior of individual HMT triggers in N rec ch . The broad shoulder in the E Pb T distribution (top-right panel) is attributable to the finite width of the N rec ch vs E Pb T correlation, which smears the contributions from different HMT triggers in E Pb T . All these structures disappear after the reweighting procedure. The results of this analysis are obtained using the MB + HMT combined dataset with event reweighting.
Owing to the relatively slow turn on of the HMT triggers as a function of E Pb T [ Fig. 1(b)], the events selected in a given E Pb T range typically receive contributions from several HMT triggers with very different weights. Hence, the effective increase in the number of events from the HMT triggers in the large E Pb T region is much smaller than the increase in the large N rec ch region. Figure 2(a) shows the correlation between E Pb T and N rec ch from MB + HMT p + Pb events. This distribution is similar to that obtained for the MB events, except that the HMT triggers greatly extend the reach in both quantities. The E Pb T value grows with increasing N rec ch , suggesting that, on average, E Pb T in the nucleus direction correlates well with the particle production at midrapidity. However, the broad distribution of E Pb T at fixed N rec ch also implies significant fluctuations. To study the relation between E Pb T and N rec ch , events are divided into narrow bins in N rec ch , and the mean and root-mean-square values of the E Pb T distribution are calculated for each bin. The results are shown in Fig. 2(b). A nearly linear relation between E Pb T and N rec ch is observed. This relationship is used to match a given N rec ch event class to the corresponding E Pb T event class. This approximately linear relation can also be parameterized [indicated by the solid line in Fig. 2 The 2PC analysis is performed in different intervals of the event activity defined by either E Pb T or N rec ch . Table II gives a list of representative event-activity classes, together with the fraction of MB + HMT events [after reweighting as shown in Fig. 2(a)] contained in each event class. The table also provides the average N rec ch and E Pb T values for each event-activity class, as well as the efficiency-corrected number of charged particles within |η| < 2.5 and p T > 0.4 GeV, N ch . The event classes defined in narrow E Pb T or N rec ch ranges are used for detailed studies of the centrality dependence of the 2PC, while the event classes in broad E Pb T or N rec ch ranges are optimized for the studies of the p T dependence. As the number of events at large E Pb T is smaller than at large N rec ch , the main results in this paper are presented for event classes defined in N rec ch .
044906-3 ch (left panels) and E Pb T (right panels) for MB and MB + HMT events before (top panels) and after (bottom panels) applying an event-by-event weight (see text). The smaller symbols in the top panels represent the distributions from the six HMT triggers listed in Table I.

C. Two-particle correlation
For a given event class, the 2PCs are measured as functions of relative azimuthal angle, φ = φ a − φ b , and relative pseudorapidity, η = η a − η b , with | η| η max = 5. 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 strength, expressed in terms of the number of pairs per trigger particle, is defined as [4][5][6] where S and B represent pair distributions constructed from the same event and from "mixed events" [4], respectively, which are then normalized by the number of trigger particles in the event. These distributions are also referred to as per-trigger yield distributions. The mixed-event distribution, B( φ, η), measures the distribution of uncorrelated pairs. The B( φ, η) distribution is constructed by choosing the two particles in the pair from different events of similar N rec ch (match to | N rec ch | < 10 tracks), E Pb T (match to | E Pb T | < 10 GeV), and z vtx (match to | z vtx | < 10 mm), so that B( φ, η) properly reflects the known detector effects in S( φ, η). The onedimensional (1D) distributions S( φ) and B( φ) are obtained by integrating S( φ, η) and B( φ, η), respectively, over a η range. The region | η| < 1 is chosen to focus on the short-range features of the correlation functions, while the region | η| > 2 is chosen to focus on the long-range features of the correlation functions. These two regions are hence referred to as the "short-range region" and the "long-range region," respectively. The normalization factors in front of the S/B ratio are chosen such that the ( φ, η)-averaged 044906-4   GeV in low-activity events, E Pb T < 10 GeV or N rec ch < 20, in the top panels, and high-activity events, E Pb T > 100 GeV or N rec ch > 220, in the bottom panels. The correlation for low-activity events shows a sharp peak centered at ( φ, η) = (0,0) owing to short-range correlations for pairs resulting from jets, high-p T resonance decays, and Bose-Einstein correlations. The correlation function also shows a broad structure at φ ∼ π from low-p T resonances, dijets, and momentum conservation that is collectively referred to as "recoil" [24] in the remainder of this paper. In the high-activity events, the correlation reveals a flat ridgelike structure at φ ∼ 0 (the "near side") that extends over the full measured η range. This η independence is quantified by integrating the 2D correlation functions over The yield associated with the near-side short-range correlation peak centered at ( φ, η) = (0,0) can then be estimated as where the second term accounts for the contribution of uncorrelated pairs and the ridge component under the near-side peak. The default value of Y n-peak is obtained with a lower-end of the integration range of η min = 2, but the value of η min is varied from 2 to 4 to check the stability of Y n-peak . The distribution at φ ∼ π (the "away side") is also broadened in high-activity events, consistent with the presence of a longrange component in addition to the recoil component [24]. This recoil component can be estimated from the low-activity events and subtracted from the high-activity events using the procedure detailed in the next section.

D. Recoil subtraction
The correlated structure above a flat pedestal in the correlation functions is calculated using a zero-yield-atminimum (ZYAM) method [4,47] following previous measurements [22][23][24], where the parameter b ZYAM represents the pedestal formed by uncorrelated pairs. A second-order polynomial fit to the 1D Y ( φ) distribution in the long-range region is used to find the location of the minimum point, φ ZYAM , and from this the value 044906-6 of b ZYAM is determined and subtracted from the 2D correlation function. The Y corr ( φ, η) functions differ, therefore, by a constant from the Y ( φ, η) functions, such as those in Fig. 3.
In low-activity events, Y corr ( φ, η) contains mainly the short-range correlation component and the recoil component. In high-activity events, the contribution from the long-range "ridge" correlation also becomes important. This long-range component of the correlation function in a given event class is obtained by estimating the short-range correlation component using the peripheral events and is then subtracted, where the Y corr in a low-activity or peripheral event class, denoted by Y corr peri , is used to estimate and subtract [hence, the superscript "sub" in Eq. (6)] the short-range correlation at the near side and the recoil at the away side. The parameter α is chosen to adjust the near-side short-range correlation yield in the peripheral events to match that in the given event class for each p a T and p b T combination, α = Y n-peak /Y n-peak peri . This scaling procedure is necessary to account for enhanced short-range correlations and away-side recoil in higher-activity events, under the assumption that the relative contribution of the near-side short-range correlation and away-side recoil is independent of the event activity. A similar rescaling procedure has also been used by the CMS Collaboration [28]. The default peripheral event class is chosen to be E Pb T < E 0 T = 10 GeV. However, the results have also been checked with other E 0 T values, as well as with a peripheral event class defined by N rec ch < 20. In the events with the highest multiplicity, the value of α determined with the default peripheral event class varies from ∼2 at p T ≈ 0.5 GeV to ∼1 for p T > 3 GeV, with a p T -dependent uncertainty of 3%-5%.
The uncertainty on b ZYAM only affects the recoil-subtracted correlation functions through the Y corr peri term in Eq. (6). This uncertainty is usually very small in high-activity p + Pb collisions, owing to their much larger pedestal level than for the peripheral event class.
Figures 4(a) and 4(b) show, respectively, the 2-D correlation functions before and after the subtraction procedure given by Eq. (6). Most of the short-range peak and away-side recoil structures are removed by the subtraction, and the remaining distributions exhibit a φ-symmetric double ridge that is almost independent of η. Figure 4(c) shows the corresponding 1D correlation functions before and after recoil subtraction in the long-range region of | η| > 2. The distribution at the near-side is not affected because the near-side short-range peak Per-trigger yield 11 11.2 11.4 is narrow in η [ Fig. 4(a)], while the away-side distribution is reduced owing to the removal of the recoil component.

E. Extraction of the azimuthal harmonics associated with long-range correlation
The azimuthal structure of the long-range correlation is studied via a Fourier decomposition similar to the approach used in the analysis of Pb + Pb collisions [7,9], where v n,n are the Fourier coefficients calculated via a discrete Fourier transformation, where N = 24 is the number of φ bins from 0 to π . The first five Fourier coefficients are calculated as functions of p a T and p b T for each event-activity class. The azimuthal anisotropy coefficients for single particles, v n , can be obtained via the factorization relation commonly used for heavy-ion collisions [7,9,48], From this the p T dependence of v n for n = 2-5 are calculated as v n p a T = v n,n p a where the default transverse momentum range for the associated particle (b) is chosen to be 1 < p b T < 3 GeV, and the Fourier coefficient as a function of the transverse momentum of the trigger particle is denoted by v n (p a T ) or simply v n (p T ) where appropriate. The extraction of v 1 requires a slight modification and is discussed separately in Sec. IV C. The factorization behavior is checked by comparing the v n (p a T ) obtained for different p b T ranges, as discussed in Sec. IV B. A similar Fourier decomposition procedure is also carried out for correlation functions without peripheral subtraction, i.e., Y ( φ). The harmonics obtained in this way are denoted by v unsub n,n and v unsub n , respectively. for | η| > 1, reflecting the removal of the away-side recoil contribution.

F. Systematic uncertainties
The systematic uncertainties in this analysis arise from pair acceptance, the ZYAM procedure, tracking efficiency, Monte Carlo consistency, residual pileup, and the recoil subtraction. Each source is discussed separately below.
The correlation functions rely on the pair acceptance functions, B( φ, η) and B( φ) in Eq. (3), to reproduce detector acceptance effects in the signal distribution. A natural way of quantifying the influence of detector effects on v n,n and v n is to express the single-particle and pair acceptance functions as Fourier series, similar to Eq. (7). The resulting coefficients for pair acceptance v det n,n are the product of those for the two single-particle acceptances v det,a n and v det,b n . In general, the pair acceptance function in φ is quite flat: The maximum fractional variation from its average value is observed to be less than 0.001 for pairs integrated in 2 < | η| < 5, and the corresponding |v det n,n | values are found to be less than 2 × 10 −4 . These v det n,n values are expected to mostly cancel in the correlation function, and only a small fraction contributes to the uncertainties in the pair acceptance function. Possible residual effects on the pair acceptance are evaluated following Ref. [9] by varying the criteria for matching in N rec ch , E Pb T , and z vtx . In each case, the residual v det n,n values are evaluated by a Fourier expansion of the ratio of the pair acceptances before and after the variation. This uncertainty varies in the range of (5-8) × 10 −6 . It is negligible for v 2 and v 3 , but becomes sizable for higher-order harmonics, particularly at low p T , where the v n values are small.
As discussed in Sec. III D, the value of b ZYAM is determined by a second-order polynomial fit of the Y ( φ) distribution. The stability of the fit is studied by varying the φ range in the fit. The uncertainty in b ZYAM depends on the local curvature around φ ZYAM and is estimated to be 0.0003-0.001 of the minimum value of Y ( φ). This uncertainty contributes directly to Y corr ( φ), but contributes to Y sub ( φ) and v n indirectly through the peripheral subtraction [see Eq. (6)]. The resulting uncertainty on v n is found to be less than 2%, for all n.
The values of per-trigger yields, Y ( φ), Y corr ( φ), and Y sub ( φ), are sensitive to the uncertainty on the tracking efficiency correction for the associated particles. This uncertainty is estimated by varying the track quality cuts and the detector material in the simulation, re-analyzing the data using corresponding Monte Carlo efficiencies, and evaluating the change in the extracted yields. The resulting uncertainty is estimated to be 2.5% owing to the track selection and 2%-3% related to our limited knowledge of the detector material. The v n,n and v n values depend only on the shape of the Y sub ( φ) distribution and hence are not sensitive to the tracking efficiency.
The analysis procedure is also validated by measuring v n values in fully simulated HIJING events [45,46] and comparing them to those measured using the generated particles. A small but systematic difference between the two results are included in the systematic uncertainties.
Nearly all of the events containing pileup are removed by the procedure described in Sec. III A. The influence of the residual pileup is evaluated by relaxing the pileup rejection criteria and then calculating the change in the per-trigger yields and v n values. The differences are taken as an estimate of the uncertainty and are found to be negligible in low event-activity classes and increase to 2% for events with E Pb T > 200 GeV or N rec ch > 300. 044906-8 According to Table II, the low-activity events used in the peripheral subtraction (E Pb T < E 0 T = 10 GeV) correspond to 28% of the MB-triggered events. The pair distributions for these events may contain a small genuine long-range component, leading to a reduction of the long-range correlation signal in a high-activity class via the peripheral subtraction procedure. The influence of this oversubtraction is evaluated by varying the definition of the low-activity events in the range of E 0 T = 5 GeV to E 0 T = 20 GeV. The Y sub ( φ) and v n values are calculated for each variation. The v n values are found to decrease approximately linearly with increasing E 0 T . The amount of oversubtraction can be estimated by extrapolating E 0 T to zero. The estimated changes of v n and Y sub ( φ) vary from less than 1% for E Pb T > 100 GeV or N rec ch > 150 and increase for lower event-activity classes approximately as 1.5/N rec ch . The relative change of v n is also found to be independent of p T . As a cross check, the analysis is repeated by defining peripheral events as N rec ch < 20. The variation of v n values is found to be consistent with the variation from varying E 0 T . The stability of the scale factor, α, is evaluated by varying the η window of the long-range region in Eq. (4). A 3%-5% uncertainty is quoted for α from these variations. The resulting uncertainty on v n for n = 2-5 is within 1% at low p T (<4 GeV) and increases to ∼10% at the highest p T . However, the v 1 extraction is directly affected by the subtraction of the recoil component, and hence the v 1 value is very sensitive to the uncertainty in α. The estimated uncertainty is 8%-12% for p T < 1 GeV and is about 20%-30% for p T > 3 GeV.
The different sources of the systematic uncertainties described above are added in quadrature to give the total systematic uncertainties for per-trigger yields and v n , which are summarized in Tables III and IV, respectively. The systematic uncertainty quoted for each source usually covers the maxmium uncertainty over the measured p T range and event-activity range. However, because v 1 (p T ) changes sign within 1.5-2.0 GeV (see Fig. 15), the relative uncertainties are quoted for p T < 1 GeV and p T > 3 GeV. The uncertainty of pair acceptance, which is less than 8 × 10 −6 for v n,n , was converted to percent uncertainties. This uncertainty can be significant at high p T . Figure 5 shows the 1D correlation functions after the ZYAM procedure, Y corr ( φ), in various ranges of p a T for a fixed p b T range of 1-3 GeV. The correlation functions are obtained in the long-range region (| η| > 2) and are shown for events selected by N rec ch 220. This event class contains a small fraction (3 × 10 −5 ) of the MB p + Pb events with highest multiplicity. The correlation functions are compared to the distributions of the recoil component, αY corr peri ( φ) in Eq. (6), estimated from the peripheral event class defined by E Pb T < 10 GeV. The scale factor α is chosen such that the near-side short-range yield matches between the two event classes [see Eq. (6) and discussion around it]. Figure 5 shows a clear near-side excess in the full p a T range studied in this analysis. An excess above the estimated recoil contribution is also observed on the away side over the same p T range.

A. Correlation functions and integrated yields
To further quantify the properties of the long-range components, the Y corr ( φ) distributions are integrated over | φ| < π/3 and | φ| > 2π/3, similar to the procedure used in previous analyses [23,24]. The integrated yields, Y int , are obtained in several event classes and are plotted as functions of p a T in Fig. 6. The near-side yields increase with trigger p T , reach a maximum at p T ∼ 3 GeV, and then decrease to a value close to zero at p T > 10 GeV. This trend is characteristic of the p T dependence of the Fourier harmonics in A + A collisions. In contrast, the away-side yields show a continuous increase across the full p T range, owing to the contribution of the recoil component that mostly results from dijets. Figure 7 shows the centrality dependence of the long-range integrated yields for the event-activity based on N rec ch (left) and E Pb T (right) for particles in 1 < p a,b T < 3 GeV range. The near-side yield is close to zero in low-activity events and increases with E Pb T or N rec ch . The away-side yield shows a similar increase as a function of E Pb T or N rec ch , but it starts at a value significantly above zero. The yield difference between these two regions is found to vary slowly with E Pb T or N rec ch , indicating  that the growth in the integrated yield with increasing event activity is similar on the near side and the away side. This behavior suggests the existence of an away-side long-range component that has a magnitude similar to the near-side long-range component. Figure 7 also shows (solid lines) the recoil component estimated from the low event-activity class (E Pb T < 10 GeV) via the rescaling procedure discussed in Sec. III D. The yield difference between the away side and the near side in this p T range is reproduced by this estimate of the recoil component. In other p T ranges, a systematic difference between the recoil component and the yield difference is observed and is attributed to the contribution of a genuine dipolar flow, v 1,1 , to the correlation function (see discussion in Sec. IV C).
[GeV] a T p 0  To quantify the φ dependence of the measured longrange correlations, the first five harmonics of the correlation functions, v 1 to v 5 , are extracted via the procedure described in Sec. III E. The following section summarizes the results for v 2 -v 5 , and the results for v 1 are discussed in Sec. IV C. Figure 8 shows the v 2 , v 3 , and v 4 obtained using the 2PC method described in Sec. III E for 1 < p b T < 3 GeV. The results are shown both before (denoted by v unsub n ) and after the subtraction of the recoil component [Eq. (6)]. The recoil contribution affects slightly the v n values for trigger p T < 3 GeV, but becomes increasingly important for higher trigger p T and higher-order harmonics. This behavior is expected as the dijet contributions, the dominant contribution to the recoil component, increase rapidly with p T (for example, see Fig. 5 or Ref. [9]). At high p T , the contribution of dijets appears as a narrow peak at the away side, leading to v unsub n coefficients with alternating sign: (−1) n [9]. In contrast, the v n values after recoil subtraction are positive across the full measured p T range. Hence, the recoil subtraction is necessary for the reliable extraction of the long-range correlations, especially at high p T . Figure 9 shows the trigger p T dependence of the v 2 -v 5 in several N rec ch event classes. The v 5 measurement is available only for three event-activity classes in a limited p T range. All flow harmonics show similar trends; i.e., they increase with The v n (p a T ) with n = 2 to 5 for six N rec ch event-activity classes obtained for | η| > 2 and the p b T range of 1-3 GeV. The error bars and shaded boxes represent the statistical and systematic uncertainties, respectively. Results in 220 N rec ch < 260 are compared to the CMS data [28] obtained by subtracting the peripheral events (the number of off-line tracks N off trk < 20), shown by the solid and dashed lines.

B. Fourier coefficients v 2 -v 5
p T up to 3-5 GeV and then decrease, but remain positive at higher p T . For all event classes, the magnitude of the v n is largest for n = 2, and decreases quickly with increasing n. The ATLAS data are compared to the measurement by the CMS experiment [28] for an event-activity class in which the number of off-line reconstructed tracks, N off trk , within |η| < 2.4 and p T > 0.4 GeV is 220 N off trk < 260. This is comparable to the 220 N rec ch < 260 event class used in the ATLAS analysis. A similar recoil removal procedure, with N off trk < 20 as the peripheral events, has been used for the CMS data. Excellent agreement is observed between the two results.
The extraction of the v n from v n,n relies on the factorization relation in Eq. (9). This factorization is checked by calculating v n using different ranges of p b T for events with N rec ch 220 as shown in Fig. 10. The factorization behavior can also be studied via the ratio [49,50] with r n = 1 for perfect factorization. The results with recoil subtraction (r n ) and without subtraction (r unsub n ) are summarized in Fig. 11, and they are shown as functions of p b T − p a T , because by construction the ratios equal 1 for p b T = p a T . This second method is limited to p a,b T 4 GeV, because requiring both particles to be at high p T reduces the number of the available pairs for v n,n (p a T ,p a T ) or v n,n (p b T ,p b T ). In contrast, for the results shown in Fig. 10, using Eqs. (9) and (10), the restriction applies to only one of the particles, i.e., p b T 4 GeV.
Results in Figs. 10 and 11 show that, in the region where the statistical uncertainty is small, the factorization holds to within a few percent for v 2 over 0.5 < p a,b T < 4 GeV, within 10% for v 3 over 0.5 < p a,b T < 3 GeV, and within 20%-30% for v 4 over 0.5 < p a,b T < 4 GeV (Fig. 10 only). Furthermore, in this p T region, the differences between r n and r unsub n are very small (<10%) as shown by Fig. 11, consistent with the observation in Fig. 8. This level of factorization is similar to what was observed in peripheral Pb + Pb collisions [9]. Figure 11 also compares the r n data with a theoretical calculation from a viscous hydrodynamic model [51]. The model predicts at most a few percent deviation of r n from 1, which is attributed to p T -dependent decorrelation effects associated with event-by-event flow fluctuations [49]. In most cases, the data are consistent with the prediction within uncertainties. Figure 12 shows the centrality dependence of v 2 , v 3 , and v 4 as functions of N rec ch and E Pb T . The results are obtained for 0.4 < p a,b T < 3 GeV, both before and after subtraction of the recoil contribution. The difference between v unsub n and v n is very small in central collisions, up to 3%-4% for both event-activity definitions. For more peripheral collisions, the difference is 044906-12  larger and reaches 20%-30% for N rec ch ∼ 40 or E Pb T ∼ 30 GeV. The sign of the difference also alternates in n (already seen in Fig. 8): i.e., v unsub n > v n for even n and v unsub n < v n for odd n. This behavior is characteristic of the influence of the away-side dijet contribution to v unsub n .
The v n values in Fig. 12 exhibit modest centrality dependence. The change of v 2 is less than 8% over 140 < N rec ch < 300 (top 0.5% of MB-triggered events) or 130 < E Pb T < 240 GeV (top 0.05% of MB-triggered events), covering about half of the full dynamic range. The centrality dependence of v 3 is stronger and exhibits a nearly linear increase with N rec ch and E Pb T . Figure 12 shows that the overall centrality dependence is similar for N rec ch and E Pb T . The correlation data [not the fit, Eq. (2)] in Fig. 2 are used to map the N rec ch dependence in the top row of Fig. 12 to a corresponding E Pb T dependence. The E Pb T dependence of v n mapped from N rec ch dependence is then compared to the directly measured E Pb T dependence in Fig. 13. Good agreement is seen for v 2 and v 3 .

C. First-order Fourier coefficient v 1
A similar analysis is performed to extract the dipolar flow v 1 . Figure 14 shows the v 1,1 values as functions of p a T in various ranges of p b T before and after the recoil subtraction.
Before the recoil subtraction, v unsub 1,1 values are always negative and decrease nearly linearly with p a T and p b T , except for the p T region around 3-4 GeV where a shoulderlike structure is seen. This shoulder is very similar to that observed in A + A collisions, which is understood as a combined contribution from the negative recoil and positive dipolar flow in this p T range [9,52] according to the form [53,54]: where M and p 2 T are the multiplicity and the average squared transverse momentum of the particles in the whole event, respectively. The negative correction term reflects the global momentum conservation contribution, which is important in low-multiplicity events and at high p T . The shoulderlike structure in Fig. 14 reflects the contribution of the dipolar flow term v 1 (p a T )v 1 (p b T ). After the recoil subtraction, the magnitude of v 1,1 is greatly reduced, suggesting that most of the momentum conservation contribution has been removed. The resulting v 1,1 values cross each other at around p a T ∼ 1. 5-2.0 GeV. This behavior is consistent with the expectation that the v 1 (p T ) function crosses zero at p T ∼ 1-2 GeV, a feature that is also observed in A + A 044906-14 collisions [9,52]. The trigger p T dependence of v 1 is obtained via a factorization procedure very similar to that discussed in where the dipolar flow in the associated p T bin, v 1 (p b T ), is defined as where sgn(p b T − p 0 T ) is the sign of the v 1 , defined to be negative for p b T < p 0 T = 1.5 GeV and positive otherwise. This function is necessary to account for the sign change of v 1 at low p T .
To obtain the v 1 (p a T ), three reference p b T ranges, 0.5-1, 3-4, and 4-5 GeV, are used to first calculate v 1 (p b T ). These values are then inserted into Eq. (13) to obtain three v 1 (p a T ) functions. The uncertainties on the v 1 (p a T ) values are calculated via an error propagation through Eqs. (13) and (14). The calculation is not possible for p b T in the range of 1-3 GeV, where the v 1,1 values are close to zero and, hence, the resulting v 1 (p b T ) have large uncertainties.
The results for v 1 (p a T ) are shown in Fig. 15 for these three reference p b T bins. They are consistent with each other. The v 1 value is negative at low p T , crosses zero at around p T ∼ 1.5 GeV, and increases to 0.1 at 4-6 GeV. This p T dependence is similar to the v 1 (p T ) measured by ATLAS experiment in Pb + Pb collisions at √ s NN = 2.76 TeV [9], except that the v 1 value in Pb + Pb collisions crosses zero at lower p T (∼1.1 GeV), which reflects the fact that the p T in Pb + Pb at √ s NN = 2.76 TeV is smaller than that in p + Pb at √ s NN = 5.02 TeV.
[GeV]  (13) and (14) in three reference p b T ranges for events with N rec ch 220. The error bars and shaded boxes represent the statistical and systematic uncertainties, respectively.

D. Comparison of v n results between high-multiplicity p + Pb and peripheral Pb + Pb collisions
In the highest multiplicity p + Pb collisions, the chargedparticle multiplicity, N rec ch , can reach more than 350 in |η| < 2.5 and E Pb T close to 300 GeV on the Pb-fragmentation side. This activity is comparable to Pb + Pb collisions at √ s NN = 2.76 TeV in the 45%-50% centrality interval, where the longrange correlation is known to be dominated by collective flow. Hence, a comparison of the v n coefficients in similar event activity for the two collision systems can improve our current understanding of the origin of the long-range correlations. The left column of Fig. 16 compares the v n values from p + Pb collisions with 220 N rec ch < 260 to the v n values for Pb + Pb collisions in the 55%-60% centrality interval from Ref. [9]. These two event classes are chosen to have similar efficiency-corrected multiplicity of charged particles with p T > 0.5 GeV and |η| < 2.5, characterized by its average value ( N ch ) and its standard deviation (σ ): N ch ± σ ≈ 259 ± 13 for p + Pb collisions and N ch ± σ ≈ 241 ± 43 for Pb + Pb collisions.
[GeV] , v 3 (middle row), and v 4 (bottom row) as functions of p T compared between p + Pb collisions with 220 N rec ch < 260 in this analysis and Pb + Pb collisions in 55%-60% centrality from Ref. [9]. The left column shows the original data with their statistical (error bars) and systematic uncertainties (shaded boxes). In the right column, the same Pb + Pb data are rescaled horizontally by a constant factor of 1.25, and the v 2 and v 4 are also downscaled by an empirical factor of 0.66 to match the p + Pb data.

044906-16
The Pb + Pb results on v n [9] were obtained via an eventplane method by correlating tracks in η > 0 (η < 0) with the event plane determined in the FCal in the opposite hemisphere. The larger v 2 values in Pb + Pb collisions can be attributed to the elliptic collision geometry of the Pb + Pb system, while the larger v 4 values are attributable to the nonlinear coupling between v 2 and v 4 in the collective expansion [55].
The v 3 data for Pb + Pb collisions are similar in magnitude to those in p + Pb collisions. However, the p T dependence of v n is different for the two systems. These observations are consistent with similar comparisons performed by the CMS experiment [28].
Recently, Basar and Teaney [56] have proposed a method to rescale the Pb + Pb data for a proper comparison to the p + Pb data. They argue that the v n (p T ) shape in the two collision systems are related to each other by a constant scale factor of K = 1.25, accounting for the difference in their p T , and that one should observe a similar v n (p T ) dependence shape after rescaling the p T measured in Pb + Pb collisions to get v n (p T /K). The difference in the overall magnitude of v 2 after the p T rescaling is entirely attributable to the elliptic geometry of Pb + Pb collisions.
To test this idea, the p T for Pb + Pb collisions are rescaled by the constant factor of 1. 25 and v n values with rescaled p T are displayed in the right column of Fig. 16. Furthermore, the magnitudes of v 2 and v 4 are also rescaled by a common empirical value of 0.66 to approximately match the magnitude of the corresponding p + Pb v n data. The rescaled v n results are shown in the right column and compared to the p + Pb v n data. They agree well with each other, in particular, in the lowp T region (p T < 2-4 GeV), where the statistical uncertainties are small.

V. SUMMARY
This paper presents measurements of 2PC functions and the first five azimuthal harmonics v 1 -v 5 in √ s NN = 5.02 TeV p + Pb collisions with a total integrated luminosity of approximately 28 nb −1 recorded by the ATLAS detector at the LHC. The 2PCs and v n coefficients are obtained as functions of p T for pairs with 2 < | η| < 5 in different intervals of event activity, defined by either N rec ch , the number of reconstructed tracks with p T > 0.4 GeV and |η| < 2.5, or E Pb T , the total transverse energy over −4.9 < η < −3.2 on the Pb-fragmentation side.
Significant long-range correlations (extending to | η| = 5) are observed for pairs at the near side (| φ| < π/3) over a wide range of transverse momentum (p T < 12 GeV) and broad ranges of N rec ch and E Pb T . A similar long-range correlation is also observed on the away side (| φ| > 2π/3), after subtracting the recoil contribution estimated using the 2PC in low-activity events.
The azimuthal structure of these long-range correlations is quantified using the Fourier coefficients v 2 -v 5 as functions of p T . The v n values increase with p T to 3-4 GeV and then decrease for higher p T , but remain positive in the measured p T range. The overall magnitude of v n (p T ) is observed to decrease with n. The magnitudes of v n also increase with both N rec ch and E Pb T . The v 2 values seem to saturate at large N rec ch or E Pb T values, while the v 3 values show a linear increase over the measured N rec ch or E Pb T range. The first-order harmonic v 1 is also extracted from the 2PC. The v 1 (p T ) function is observed to change sign at p T ≈ 1.5-2.0 GeV and to increase to about 0.1 at p T > 4 GeV.
The extracted v 2 (p T ), v 3 (p T ), and v 4 (p T ) are compared to the v n coefficients in Pb + Pb collisions at √ s NN = 2.76 TeV with similar N rec ch . After applying a scale factor of K = 1.25 that accounts for the difference of mean p T in the two collision systems as suggested in Ref. [56], the shape of the v n (p T /K) distribution in Pb + Pb collision is found to be similar to the shape of v n (p T ) distribution in p + Pb collisions. This suggests that the long-range ridge correlations in high-multiplicity p + Pb collisions and peripheral Pb + Pb collisions are driven by similar dynamics.

ACKNOWLEDGMENTS
We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions, without whom ATLAS could not be operated efficiently. We acknowledge the support of ANPCyT