UvA-DARE (Digital Academic Repository) Measurement of forward-backward multiplicity correlations in lead-lead, proton-lead, and proton-proton collisions with the ATLAS detector

Two-particle pseudorapidity correlations are measured in √ s NN = 2 . 76TeV Pb + Pb, √ s NN = 5 . 02 TeV p + Pb, and √ s = 13 TeV pp collisions at the Large Hadron Collider (LHC), with total integrated luminosities of approximately 7 μ b − 1 , 28 nb − 1 , and 65 nb − 1 , respectively. The correlation function C N ( η 1 ,η 2 ) is measured as a function of event multiplicity using charged particles in the pseudorapidity range | η | < 2 . 4. The correlation function contains a signiﬁcant short-range component, which is estimated and subtracted. After removal of the short-range component, the shape of the correlation function is described approximately by 1 + (cid:4) a 21 (cid:5) 1 / 2 η 1 η 2 in all collision systems over the full multiplicity range. The values of (cid:4) a 21 (cid:5) 1 / 2 are consistent for the opposite-charge pairs and same-charge pairs, and for the three collision systems at similar multiplicity. The values of (cid:4) a 21 (cid:5) 1 / 2 and the magnitude of the short-range component both follow a power-law dependence on the event multiplicity. The short-range component in p + Pb collisions, after symmetrizing the proton and lead directions, is found to be smaller at a given η than in pp collisions with comparable multiplicity. DOI: 10.1103/PhysRevC.95.064914


Measurement of forward-backward multiplicity correlations in lead-lead, proton-lead, and proton-proton collisions with the ATLAS detector
, with total integrated luminosities of approximately 7 μb −1 , 28 nb −1 , and 65 nb −1 , respectively. The correlation function C N (η 1 ,η 2 ) is measured as a function of event multiplicity using charged particles in the pseudorapidity range |η| < 2.4. The correlation function contains a significant short-range component, which is estimated and subtracted. After removal of the short-range component, the shape of the correlation function is described approximately by 1 + a 2 1 1/2 η 1 η 2 in all collision systems over the full multiplicity range. The values of a 2 1 1/2 are consistent for the opposite-charge pairs and same-charge pairs, and for the three collision systems at similar multiplicity. The values of a 2 1 1/2 and the magnitude of the short-range component both follow a power-law dependence on the event multiplicity. The short-range component in p + Pb collisions, after symmetrizing the proton and lead directions, is found to be smaller at a given η than in pp collisions with comparable multiplicity. DOI: 10.1103/PhysRevC.95.064914

I. INTRODUCTION
Heavy-ion collisions at Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) create hot, dense matter whose space-time evolution can be well described by relativistic viscous hydrodynamics [1,2]. Owing to strong event-by-event (EbyE) density fluctuations in the initial state, the space-time evolution of the produced matter in the final state also fluctuates event to event. These fluctuations may lead to correlations of particle multiplicity in momentum space in the transverse and longitudinal directions with respect to the collision axis. Studies of the multiplicity correlation in the transverse plane have revealed strong harmonic modulation of the particle densities in the azimuthal angle, commonly referred to as harmonic flow. The measurements of harmonic flow coefficients v n [3][4][5][6] and their EbyE fluctuations [7][8][9][10] have placed important constraints on the properties of the medium and transverse energy density fluctuations in the initial state.
Two-particle correlations in the transverse plane have also been studied in high-multiplicity pp [11][12][13] and p + Pb [14][15][16][17][18] collisions, and these studies have revealed features that bear considerable similarity to those observed in heavy-ion collisions. These findings have generated many theoretical interpretations [19], and much discussion as to whether the mechanisms that result in the observed correlations are or are not fundamentally the same in the different collision systems.
This paper reports measurements of multiplicity correlations in the longitudinal direction in pp, p + Pb, and Pb+Pb collisions, which are sensitive to the early-time density fluctuations in pseudorapidity (η) [1,2]. These density fluctuations generate long-range correlations (LRC) at the early stages of the collision, well before the onset of any collective behavior, and appear as correlations of the multiplicity densities of produced particles separated in η. For example, the EbyE differences between the partonic flux in the target and the projectile may lead to a long-range asymmetry of the produced system [20][21][22], which manifests itself as a correlation between the multiplicity densities of final-state particles with large η separation.
Longitudinal multiplicity correlations can also be generated during the space-time evolution in the final state as resonance decays, single-jet fragmentation, and Bose-Einstein correlations. These latter correlations are typically localized over a smaller range of η, and are commonly referred to as short-range correlations (SRC). On the other hand, dijet fragmentation may contribute to the LRC if the η separation between the two jets is large.
Many previous studies are based on forward-backward (FB) correlations of particle multiplicity in two η ranges symmetric around the center-of-mass of the collision systems, including e + e − [23], pp [24][25][26][27], and A + A [28,29] collisions where a significant anticorrelation between forward and backward multiplicities has been identified. Recently, the study of multiplicity correlations has been generalized by decomposing the correlation function into orthogonal Legendre polynomial functions, or more generally into principal components, each representing a unique component of the measured FB correlation [21,30].
Particle production in pp collisions is usually described by QCD-inspired models, such as PYTHIA [31] and EPOS [32], implemented in Monte Carlo (MC) event generators with free parameters that are tuned to describe experimental measurements. Previous studies show that these models can generally describe the η and p T dependence of the inclusive charged-particle production [33,34], as well as the underlying event accompanying various hard-scattering processes [35,36]. In many such models, events with large charged-particle multiplicity are produced through multiple parton-parton interactions (MPI), which naturally serve as sources for the FB multiplicity asymmetry describe above. Therefore, a detailed measurement of pseudorapidity correlation in pp collisions also provides new constraints on the longitudinal dynamics of MPI processes in these models.
The two-particle correlation function in pseudorapidity is defined as [37,38] C(η 1 ,η 2 ) = N (η 1 )N(η 2 ) N (η 1 ) N(η 2 ) ≡ ρ(η 1 )ρ(η 2 ) , where N(η) is the multiplicity density distribution in a single event and N(η) is the average distribution for a given event-multiplicity class. The correlation function is directly related to a single-particle quantity ρ(η), which characterizes the fluctuation of multiplicity in a single event relative to the average shape of the event class. Following Refs. [21,38], ρ(η) in the interval [−Y,Y ] is written in terms of Legendre polynomials: ρ(η) ∝ 1+ n a n T n (η), T n (η) ≡ 2n +1 3 Y P n η Y , (2) and the scale factor in Eq. (2) is chosen such that T 1 (η) = η. 1 Using Eqs. (1) and (2), the correlation function C can be expressed in terms of the T n , which involve terms in a 0 a 0 , a 0 a n , and a n a m , with n,m 1. Terms involving a 0 reflect multiplicity fluctuations in the given event class, while the dynamical fluctuations between particles at different pseudorapidities in events of fixed multiplicity are captured by the terms in a n a m , n,m 1. It is the study of these dynamical fluctuations that is the goal of this analysis.
As discussed in more detail in Ref. [38], the terms involving a 0 a n can be removed, provided all deviations from 1 are small, by defining where with a similar expression for C p (η 2 ). The quantities C p (η 1 ) and C p (η 2 ) are referred to as the single-particle modes. The a 0 a 0 term can be removed by renormalizing average value in the η 1 , η 2 phase space to be 1. The final result is and a n,m ≡ a n a m . 1 The T n (η) also satisfy Y −Y T n (η)dη = 0 for n 1, and . From the definition of ρ(η) in Eq. (1), it follows that ∞ n=0 a n T n (η) = 0.
The two-particle Legendre coefficients can be calculated directly from the measured correlation function: The two-particle correlation method measures, in effect, the root-mean-square (rms) values of the EbyE a n , a 2 n 1 /2 , or the cross correlation between a n and a m , a n a m . The correlation functions satisfy the symmetry condition C(η 1 ,η 2 ) = C(η 2 ,η 1 ) and C N (η 1 ,η 2 ) = C N (η 2 ,η 1 ). This paper presents a measurement of the two-dimensional (2-D) correlation function C N (η 1 ,η 2 ) over the pseudorapidity range of |η| < 2.4 in √ s NN = 2.76 TeV Pb+Pb, √ s NN = 5.02 TeV p + Pb, and √ s = 13 TeV pp collisions, using the ATLAS detector. 2 The analysis is performed using events for which the total number of reconstructed charged particles, N rec ch , with |η| < 2.5 and transverse momentum p T > 0.4 GeV, is in the range 10 N rec ch < 300. Both the Pb+Pb and p + Pb data cover this range of N rec ch , but for pp the range extends only to approximately 160. The measured C N (η 1 ,η 2 ) is separated into a short-range component δ SRC (η 1 ,η 2 ) and C sub N (η 1 ,η 2 ), which contains the long-range component. The nature of the FB fluctuation in each collision system is studied by projections as well as Legendre coefficients a n a m of C sub N (η 1 ,η 2 ). The magnitudes of the FB fluctuations are compared for the three systems at similar event multiplicity. A comparison is also made between the pp data and QCD-inspired models.

II. ATLAS DETECTOR AND TRIGGER
The ATLAS detector [39] provides nearly full solidangle coverage of the collision point with tracking detectors, calorimeters, and muon chambers and is well suited for measurement of two-particle correlations over a large pseudorapidity range. The measurements were performed using the inner detector (ID), minimum-bias trigger scintillators (MBTS), the forward calorimeter (FCal), and the zero-degree calorimeters (ZDC). The ID detects charged particles within |η| < 2.5 using a combination of silicon pixel detectors, silicon microstrip detectors (SCT), and a straw-tube transition radiation tracker (TRT), all immersed in a 2-T axial magnetic field [40]. An additional pixel layer, the "insertable B layer" (IBL) [41,42] installed between run 1 and run 2 (2013-2015), is used in the 13-TeV pp measurements. The MBTS system detects charged particles over 2.1 |η| 3.9 using two hodoscopes of counters positioned at z = ± 3.6 m. The FCal consists of three sampling layers, longitudinal in shower depth, 2 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z axis along the beam pipe. The x axis points from the IP to the center of the LHC ring, and the y axis points upward. Cylindrical coordinates (r,φ) are used in the transverse plane, φ being the azimuthal angle around the beam pipe. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). 064914-2 and covers 3.2 < |η| < 4.9. The ZDC, available in the Pb+Pb and p + Pb runs, are positioned at ±140 m from the collision point, detecting neutrons and photons with |η| > 8. 3. This analysis uses approximately 7 μb −1 of Pb+Pb data, 28 nb −1 of p + Pb data, and 65 nb −1 of pp data taken by the ATLAS experiment at the LHC. The Pb+Pb data were collected in 2010 at a nucleon-nucleon center-of-mass energy √ s NN = 2.76 TeV. The p + Pb data were collected in 2013, when 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 results in a rapidity shift of 0.47 of the nucleonnucleon center-of-mass frame towards the proton beam direction relative to the laboratory rest frame. The pp data were collected during a low-luminosity operation of the LHC in June and August 2015 at collision energy √ s = 13 TeV. The ATLAS trigger system [43] consists of a level-1 (L1) trigger implemented using a combination of dedicated electronics and programmable logic, and a high-level trigger (HLT) implemented in processors. The HLT reconstructs charged-particle tracks using methods similar to those applied in the offline analysis, allowing high-multiplicity track (HMT) triggers that select on the number of tracks having p T > 0.4 GeV associated with the vertex with the largest number of associated tracks (primary vertex). The Pb+Pb data used in the analysis are collected by a minimum-bias trigger, while the pp and p + Pb data are collected by a minimum-bias trigger and HMT triggers.
The Pb+Pb trigger requires signals in two ZDCs or either of the two MBTS counters. The ZDC trigger thresholds on each side are set below the peak corresponding to a single neutron. A timing requirement based on signals from each side of the MBTS is imposed to remove beam backgrounds. The minimum-bias trigger for p + Pb is similar, except that only the ZDC on the Pb-fragmentation side is used. For pp collisions, the minimum-bias trigger requires one or more signals in the MBTS.
Two distinct HMT triggers are used for the 13-TeV pp analysis. The first trigger selected events at L1 that have a signal in at least one counter on each side of the MBTS, and at the HLT have at least 900 SCT hits and 60 tracks associated with a primary vertex. The second trigger selects events with a total transverse energy of more than 10 GeV at L1 and at least 1400 SCT hits and 90 tracks associated to a primary vertex at HLT. For the p + Pb data, the HMT triggers were formed from a combination of L1 triggers that applied different thresholds for total transverse energy measured over 3.2 < |η| < 4.9 in the FCal and HLT triggers that placed minimum requirements on the number of reconstructed tracks. Details of the minimumbias and HMT triggers can be found in Refs. [12,33] and Refs. [18,44] for the pp and p + Pb collisions, respectively.

A. Event and track selection
The offline event selection for the p + Pb and pp data requires at least one reconstructed vertex with its z position satisfying |z vtx | < 100 mm. The mean collision rate per crossing μ is around 0.03 for p + Pb data, between 0.002 and 0.04 for the June 2015 pp data, and between 0.05 and 0.6 for the August 2015 pp data. Events containing multiple collisions (pileup) are suppressed by rejecting events with more than one good reconstructed vertex, and results are found to be consistent between the June and August datasets. For the p + Pb events, a time difference of | t| < 10 ns is also required between signals in the MBTS counters on either side of the interaction point to suppress noncollision backgrounds.
The offline event selection for the Pb+Pb data requires a reconstructed vertex with its z position satisfying |z vtx | < 100 mm. The selection also requires a time difference | t| < 3 ns between signals in the MBTS trigger counters on either side of the interaction point to suppress noncollision backgrounds. A coincidence between the ZDC signals at forward and backward pseudorapidity is required to reject a variety of background processes, while maintaining more than 98% efficiency for inelastic processes.
Charged-particle tracks and primary vertices are reconstructed in the ID using algorithms whose implementation was optimized for better performance between LHC runs 1 and 2. In order to compare directly the p + Pb and Pb+Pb systems using event selections based on the multiplicity of the collisions, a subset of data from peripheral Pb+Pb collisions, collected during the 2010 LHC heavy-ion run with a minimum-bias trigger, was reanalyzed using the same track reconstruction algorithm as that used for p + Pb collisions. For the p + Pb and Pb+Pb analyses, tracks are required to have a p T -dependent minimum number of hits in the SCT, and 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. A description of the 2010 Pb+Pb data and 2013 p + Pb data can be found in Refs. [5] and [45], respectively.
For the 13-TeV pp analysis, the track selection criteria were modified slightly to profit from the presence of the IBL in run 2. Furthermore, the requirements of |d z 0 | < 1.5 mm and |z 0 sin θ | < 1.5 mm are applied, where d z 0 is the transverse impact parameter of the track relative to the average beam position. These selection criteria are the same as those in Refs. [12,33].
In this analysis, the correlation functions are constructed using tracks passing the above selection requirements and which have p T > 0.2 GeV and |η| < 2.4. However, slightly different kinematic requirements, p T > 0.4 GeV and |η| < 2.5, are used to count the number of reconstructed charged particles in the event, denoted by N rec ch , to be consistent with the requirements used in the HLT. Figure 1 compares the normalized N rec ch distributions of events in the three colliding systems. The distribution decreases slowly in the Pb+Pb system, but decreases much faster in the p + Pb and pp systems. A major goal of the analysis is to compare the correlation function from the three collisions systems at similar N rec ch values, which can reveal whether the FB multiplicity fluctuation is controlled by the collision geometry or the overall activity of the event.
The efficiency of the track reconstruction and track selection requirements, (η,p T ), is evaluated using simulated p + Pb or Pb+Pb events produced with the HIJING event generator [46] or simulated pp events from the PYTHIA 8 064914-3 [31] event generator using parameter settings according to the so-called A2 tune [47]. The MC sample for Pb+Pb events in the multiplicity region of interest was very small; therefore the reconstruction efficiency for Pb+Pb was taken from the larger p + Pb sample. The p + Pb efficiency was found to be consistent with the efficiency from the Pb+Pb MC simulation, but of much higher precision. The response of the detector to these MC events is simulated using GEANT4 [48,49] and the resulting events are reconstructed with the same algorithms that are applied to the data. The efficiencies for the three datasets are similar for events with similar multiplicity. Small differences are due to changes in the detector conditions in run 1 and changes in the reconstruction algorithm between runs 1 and 2. In the simulated events, the efficiency reduces the measured charged-particle multiplicity relative to the event generator multiplicity for primary charged particles. 3 The reduction factors for N rec ch and the associated efficiency uncertainties are b = 1.29 ± 0.05, 1.29 ± 0.05, and 1.18 ± 0.05 for Pb+Pb, p + Pb, and pp collisions, respectively. The values of these reduction factors are found to be independent of multiplicity over the N rec ch range used in this analysis, 10 N rec ch < 300. Therefore, these factors are used to multiply N rec ch to obtain the efficiency-corrected average number of charged particles with >0.4 GeV and |η| < 2.5, N ch = bN rec ch . The quantity N ch is used when presenting the multiplicity dependence of the SRC and the LRC.

B. Two-particle correlations
The two-particle correlation function defined in Eq. (1) is calculated as the ratio of distributions for same-event pairs S(η 1 ,η 2 ) ∝ N (η 1 )N(η 2 ) and mixed-event pairs B(η 1 ,η 2 ) ∝ 3 For Pb+Pb and p + Pb simulation, the event generator multiplicity includes charged particles that originate directly from the collision or result from decays of particles with cτ < 10 mm. The definition for primary charged particles is somewhat stronger in the pp simulation [33]. N (η 1 ) N (η 2 ) [5]: The mixed-event pair distribution is constructed by combining tracks from one event with those from another event with similar N rec ch (matched within two tracks) and z vtx (matched within 2.5 mm). The events are also required to be close to each other in time to account for possible time-dependent variation of the detector conditions. The mixed-event distribution should account properly for detector inefficiencies and nonuniformity but does not contain physical correlations. The normalization of C(η 1 ,η 2 ) is chosen such that its average value in the (η 1 ,η 2 ) plane is one. The correlation function satisfies the symmetry C(η 1 ,η 2 ) = C(η 2 ,η 1 ) and, for a symmetric collision system, C(η 1 ,η 2 ) = C(−η 1 ,−η 2 ). Therefore, for pp and Pb+Pb collisions, all pairs are entered into one quadrant of the (η 1 ,η 2 ) space defined by η − ≡ η 1 − η 2 > 0 and η + ≡ η 1 + η 2 > 0 and then reflected to the other quadrants. For p + Pb collisions, all pairs are entered into one half of the (η 1 ,η 2 ) space defined by η 1 − η 2 > 0 and then reflected to the other half. To correct S(η 1 ,η 2 ) and B(η 1 ,η 2 ) for the individual inefficiencies of particles in the pair, the pairs are weighted by the inverse product of their tracking efficiencies 1/( 1 2 ). Remaining detector distortions not accounted for by the reconstruction efficiency largely cancel in the same-event to mixed-event ratio.
In a separate analysis, the correlation functions in p + Pb collisions are also symmetrized in the same way as for Pb+Pb and pp collisions such that C(η 1 ,η 2 ) = C(−η 1 ,−η 2 ), and they are compared with correlation functions obtained for symmetric collision systems. This symmetrized p + Pb correlation function is used only at the end of Sec. IV, in relation to Fig. 16. In all other cases, the p + Pb correlation function is unsymmetrized.

C. Outline of the procedure for separating SRC and LRC
As explained in the introduction, the aim of this analysis is to measure and parametrize the long-range correlation, which requires the separation and subtraction of the short-range component. The separation of SRC and LRC is quite involved and so is briefly summarized here, with details left to the relevant later sections.
The core of the separation method is to exploit the difference between the correlations for opposite-charge and same-charge pairs, C +− (η 1 ,η 2 ) and C ±± (η 1 ,η 2 ), respectively. The SRC component centered around η − (≡η 1 − η 2 ) ∼ 0 is found to be much stronger for opposite-charge pairs, primarily due to local charge conservation, while the LRC and single-particle modes are expected to be independent of the charge combination. With this assumption, the ratio is given approximately by This analysis assumes further that the dependence of δ SRC on η − and η + (≡η 1 + η 2 ) factorizes and that the dependence , and the ratio R(η 1 ,η 2 ) = C +− (η 1 ,η 2 )/C ±± (η 1 ,η 2 ) (top right panel) for Pb+Pb collisions with 200 N rec ch < 220. The width and magnitude of the short-range peak of the ratio are shown, as a function of η + , in the lower middle panel and lower right panels, respectively. The error bars represent the statistical uncertainties, and the solid lines indicate a quadratic fit. The dotted line in the bottom right panel serves to indicate better the deviation of f (η + ) from 1.
where g +− (η − ) and g ±± (η − ) are allowed to differ in both shape and magnitude. With these assumptions, 4 f (η + ) can be determined from R by suitable integration over η − , as described in Sec. III D.
To complete the determination of δ ±± SRC , the quantity g ±± is determined and parameterized from suitable projections of C ±± N (η + ,η − ) in the η − direction, as described in Sec. III E. The use of C ±± N rather than C ±± is because the former does not contain the single-particle modes. The procedure to obtain a correlation function with the SRC subtracted is also described in Sec. III E. With δ ±± SRC determined, δ +− SRC is obtained directly from Eq. (9). The δ ±± SRC and δ +− SRC are then averaged to obtain the SRC for all charge combinations, δ SRC . Figure 2 shows separately the correlation functions for same-charge pairs and opposite-charge pairs from Pb+Pb 4 The validity of the various assumptions is confirmed in the data from the extracted δ +− SRC (η + ,η − ) and δ ±± SRC (η + ,η − ) after applying the separation procedure. collisions with 200 N rec ch < 220. The ratio of the two, R(η 1 ,η 2 ) via Eq. (8), is shown in the top right panel. The correlation functions show a narrow ridgelike shape along η 1 ≈ η 2 or η − ≈ 0, and a falloff towards the corners at η 1 = −η 2 ≈ ±2.4. The magnitude of the ridge for the oppositecharge pairs is stronger than that for the same-charge pairs, which is characteristic of the influence from SRC from jet fragmentation or resonance decays. In regions away from the SRC, i.e., large values of |η − |, the ratio approaches unity, suggesting that the magnitude of the LRC is independent of the charge combinations. To quantify the shape of the SRC in the ratio along η + , R is expressed in terms of η + and η − , R(η + ,η − ), and the following quantity is calculated:

D. Probing the SRC via the same-charge and opposite-charge correlations
As shown in Fig. 2, the quantity f (η + ) is nearly constant in Pb+Pb collisions, implying that the SRC is consistent with being independent of η + . To quantify the shape of the SRC along the η − direction, R(η + ,η − ) is fit to a Gaussian function in slices of η + . The width, as shown in the bottom middle panel of Fig. 2, is constant, which may suggest that the shape of the SRC in η − is the same for different η + slices.  Figure 3 shows the correlation function in p + Pb collisions with multiplicity similar to the Pb+Pb data in Fig. 2. The correlation function shows a significant asymmetry between the proton-going side (positive η + ) and lead-going side (negative η + ). However, much of this asymmetry appears to be confined to a small |η − | region where the SRC dominates. The magnitude of the SRC, estimated by f (η + ) shown in the bottom-right panel, increases by about 50% from the lead-going side (negative η + ) to the proton-going side (positive η + ), but the width of the SRC in η − is independent of η + as shown in the bottom-middle panel. In contrast, the LRC has no dependence on the charge combinations, since the value of R approaches unity at large |η − |. Figure 4 shows the width in η − of the short-range component as a function of N ch in the three collision systems. The width is obtained as the Gaussian width of R(η + ,η − ) in the η − direction, and then averaged over η + as the width is observed to be independent of η + , as shown in Figs. 2 and 3. This width reflects the extent of the short-range correlation in η, and it is observed to decrease with increasing N ch in all collision systems. At the same N ch value, the width is smallest in pp collisions and largest in Pb+Pb collisions. In Fig. 5, the width of the short-range component from pp data is compared with PYTHIA 8 based on the A2 tune [50] and EPOS based on the LHC tune [32]. The width is underestimated by PYTHIA 8 A2 and overestimated by EPOS LHC.
TeV, compared between data and two models. The y axis is zero suppressed to demonstrate better the difference between data and models.

E. Separation of the SRC and the LRC
As discussed in Sec. III C, the ratio of the correlation function between opposite-charge and same-charge pairs R(η + ,η − ) is the key to the separation of the SRC and LRC. Following Eqs. (8) and (9), this ratio can be approximated by where f (η + ) describes the shape along η + and can be calculated via Eq. (10). The functions g +− and g ±± describe the SRC along the η − direction for the two charge combinations, which differ in both magnitude and shape. In order to estimate the g ±± (η − ) function for same-charged pairs, the C N (η + ,η − ) distributions for same-charge pairs are projected into one-dimensional (1-D) η − distributions over a narrow slice |η + | < 0.4. The distributions are denoted by C N (η − ). They are shown, after a small iterative correction discussed below, in the second column of Fig. 6 for the same-charge pairs in Pb+Pb and p + Pb collisions. The SRC appears as a narrow peak on top of a distribution that has an approximately quadratic shape. Therefore, a quadratic fit is applied to the data in the region of |η − | > 1.5, and the difference between the data and fit in the |η − | < 2 region is taken as the estimated SRC component or the g ±± (η − ) function, which is assumed to be zero for |η − | > 2. This range (|η − | > 1.5) is about twice the width of the short-range peak in the R(η + ,η − ) distribution along the η − direction (examples are given in the bottom middle panel of Figs. 2 and 3). This width is observed to decrease from 1.0 to 0.7 as a function of N rec ch in the p + Pb collisions, and is slightly broader in Pb+Pb collisions and slightly narrower in pp collisions at the same N rec ch . The range of the fit is varied from |η − | > 1.0 to |η − | > 2.0 to check the sensitivity of the SRC estimation, and the variation is included in the final systematic uncertainties. Furthermore, this study is also repeated for C N (η − ) obtained in several other η + slices within |η + | < 1.2, and consistent results are obtained. Once the distribution g ±± (η − ) for same-charge pairs is obtained from the fit, it is multiplied by the f (η + ) function calculated from R(η 1 ,η 2 ) using Eq. (10), to obtain the δ SRC (η 1 ,η 2 ) from Eq. (11) in the full phase space. Subtracting this distribution from the C N (η 1 ,η 2 ) distribution, one obtains the initial estimate of the correlation function containing mostly the LRC component.
The LRC obtained via this procedure is still affected by a small bias from the SRC via the normalization procedure of Eq. (3). This bias appears because the δ SRC (η 1 ,η 2 ) contribution is removed from the numerator but is still included in the denominator via C p (η). This contribution is not uniform in η: If the first particle is near midrapidity η 1 ≈ 0 then all pairs in δ SRC (η 1 ,η 2 ) contribute to C p (η 1 ), whereas if the first particle is near the edge of the acceptance η 1 ≈ ±Y then only half of the pairs in δ SRC (η 1 ,η 2 ) contribute to C p (η 1 ). The acceptance bias in C p is removed via a simple iterative procedure: First, the δ SRC contribution determined from the above procedure is used to eliminate the SRC contribution to the single-particle mode: with a similar expression for C sub p (η 2 ). The C sub p (η 1 ), C sub p (η 2 ) are then used to redefine the C N function: This distribution, which is very close to the distribution before correction, is shown in the second column of Fig. 6 for projection over a narrow slice |η + | < 0.4. The estimation of δ SRC (η 1 ,η 2 ) is repeated using the previously described procedure for the C N (η 1 ,η 2 ), and the extracted distribution is shown in the third column of Fig. 6. Subtracting this distribution from C N (η 1 ,η 2 ), one obtains the correlation function containing only the LRC component. The resulting correlation function, denoted C sub N (η 1 ,η 2 ), is shown in the last column of Fig. 6. The results presented in this paper are obtained using the iterative procedure discussed above. In most cases, the results obtained from the iterative procedure are consistent with the one obtained without iteration. In p + Pb and Pb+Pb collisions, where the SRC component is small, the difference between the two methods is found to be less than 2%. In pp collisions with N rec ch > 100, the difference between the two methods reaches 4% where the SRC is large and therefore the bias correction is more important.
In principle, the same analysis procedure can be applied to opposite-charge and all-charge pairs. However, due to the much larger SRC, the extracted LRC for opposite-charge pairs has larger uncertainties. Instead, the SRC for opposite-charge pairs is obtained directly by rearranging the terms in Eq. (9) as The SRC for all-charge pairs is calculated as the average of δ ±± SRC and δ +− SRC weighted by the number of same-charge and opposite-charge pairs. The LRC is then obtained by subtracting the SRC from the modified C N (η 1 ,η 2 ) using the same procedure as that for the same-charge pairs.
For pp collisions, the pseudorapidity correlations are also compared with the PYTHIA 8 A2 and EPOS LHC event generators mentioned above. The analysis procedure used on the data is repeated for the two models in order to extract the SRC and LRC components. The correlation is carried out on the generated, as opposed to the reconstructed, charged particles.

F. Quantifying the magnitude of the forward-backward multiplicity fluctuations
In the azimuthal correlation analysis, the azimuthal structure of the correlation function is characterized by harmonic coefficients v n obtained via a Fourier decomposition [5,51]. A similar approach can be applied for pseudorapidity correlations [21,38]. Following Eq. (5), the correlation functions are expanded into Legendre polynomial functions, and the twoparticle Legendre coefficients a n a m are calculated directly from the correlation function according to Eq. (6). The two-particle correlation method measures, in effect, the rms values of the EbyE a n , and the final results for the coefficients are presented in terms of √ | a n a m |. As a consequence of the condition for a symmetric collision system, the odd and even coefficients should be uncorrelated in pp and Pb+Pb collisions: a n,n+1 = a n a n+1 = 0.
However, even in p + Pb collisions, the correlation function after SRC removal, C sub N (η 1 ,η 2 ), is observed to be nearly symmetric between η and −η (right column of Fig. 6), and hence the a n a n+1 values are very small and considered to be negligible in this paper.
The shape of the first two Legendre bases in 2-D are shown in Fig. 7. The first basis function has the shape of η 1 × η 2 and is directly sensitive to the FB asymmetry of the EbyE fluctuation. The second basis function has a quadratic shape in the η 1 and η 2 directions and is sensitive to the EbyE fluctuation in the width of the N (η) distribution. It is shown in Sec. IV that the data require only the first term, in which case the shape of the correlation function can be approximated by Therefore, a quadratic shape is expected along the two diagonal directions η + and η − of the correlation function, and the a 2 1 1 /2 064914-8 coefficient can be calculated by a simple quadratic fit of C sub N in narrow slices of η − or η + . Alternatively, a 2 1 1 /2 can also be estimated from a correlator constructed from a simple ratio: where η ref is a narrow interval of 0.2. This correlator has the advantage that most of the single-particle modes are even functions in η, so they cancel in the ratios. Therefore, this correlator provides a robust consistency check of any potential bias induced by the renormalization procedure of Eq. (3). A similar quantity can also be calculated for C N (η 1 ,η 2 ), denoted by r N (η,η ref ). In summary, this paper uses the following four different methods to estimate a 2 The three fitting methods (2, 3, and 4) use the correlation function in limited and largely nonoverlapping regions of the η 1 and η 2 phase space, and therefore are independent of each other and largely independent of the Legendre decomposition method. Moreover, if the correlation function is dominated by the a 2 1 term, the results from all four methods should be consistent.

G. Systematic uncertainties
The systematic uncertainties in this analysis arise from the event mixing, track reconstruction and selection efficiency, pair acceptance, and use of simulated events to test the analysis process by comparing results from the generated charged particles with those from reconstructed tracks. These uncertainties apply to C N (η 1 ,η 2 ) or C sub N (η 1 ,η 2 ) and the associated Legendre coefficients. However, the systematic uncertainty for C sub N (η 1 ,η 2 ) also depends on the procedure for separating the SRC from the LRC.
A natural way of quantifying these systematic uncertainties, used in this analysis, is to calculate C N (η 1 ,η 2 ) or C sub N (η 1 ,η 2 ) under a different condition and then construct the ratio to the default analysis: D(η 1 ,η 2 ). The average deviation of D(η 1 ,η 2 ) from unity can be compared with the correlation signal to estimate the systematic uncertainties in the correlation function. The same D(η 1 ,η 2 ) function can also be expanded into a Legendre series (Eq. (5)), and the resulting coefficients a d n,m can be used to estimate the systematic uncertainties for the a n,m coefficients. For the three fitting methods discussed in Sec. III F, the fits are repeated for each check to estimate the uncertainties in the resulting a 2 1 1 /2 values. These uncertainties are not always the same for C N and C sub N because C sub N is not sensitive to the variation in the short-range region, η − ≈ 0. In the following, the uncertainty from each source is discussed.
The main source of uncertainty for C sub N (η 1 ,η 2 ) arises from the procedure to separate the SRC and the LRC. Since the estimated SRC component for the opposite-charge pairs is more than a factor of two larger than that for the same-charge pairs (e.g., Figs. 2 and 3), the difference between C sub,+− N and C sub,±± N is a conservative check of the robustness of the subtraction procedure. This difference is typically small for events with large N rec ch , and it is found to be within 0.2-2.2% of the correlation signal and 1-6% for a 2 1 1 /2 in the three collision systems. The stability of LRC is also checked by varying the fit range and varying the η + slice used to obtain the δ SRC (η − ) distribution for same-charge pairs. This uncertainty amounts to 1-2% in the correlation signal and 1-5% for a 2 1 1 /2 in Pb+Pb collisions, and is larger in p + Pb and pp collisions due to a stronger SRC for events with the same N rec ch . Uncertainties due to the event mixing are evaluated by varying the criteria for matching events in N rec ch and z vtx . The a d n,m values are calculated for each case. The uncertainty from variation of the matching range in z vtx is less than 0.5% of the correlation signal for both C N and C sub N . The bin size in N rec ch for event matching is varied such that the number of events in each bin varies by a factor of three. Most of the changes appear as modulations of the projections of the correlation function in η 1 or η 2 as defined in Eq. (4), and the renormalized correlation functions C N (η 1 ,η 2 ) and C sub N (η 1 ,η 2 ) are very stable. The difference between different variations amounts to at most 2% of the correlation signal or a 2 1 1 /2 . The analysis is also repeated separately for events with |z vtx | < 50 mm and 50 < |z vtx | < 100 mm. Good agreement is seen between the two. To evaluate the stability of the correlation function, the entire dataset is divided into several groups of runs, and the correlation functions and a n coefficients are calculated for each group. The results are found to be consistent within 2% for a 2 1 1 /2 . The 13-TeV pp results are obtained from the June 2015 and August 2015 datasets with different μ values. The influence of the residual pileup is evaluated by comparing the results obtained separately from these two running periods, and no systematic difference is observed between the results.
The shape of the correlation function is not very sensitive to the uncertainty in the tracking efficiency correction, since this correction is applied in both the numerator and denominator. On the other hand, both the correlation signal and reconstruction efficiency are observed to increase with p T , and hence the correlation signal and associated a n a m coefficients are expected to be smaller when corrected for reconstruction efficiency. Indeed, a 1-2% decrease in a 2 n 1 /2 is observed after applying this correction. This change is conservatively included in the systematic uncertainty.
The correlation function C N (η 1 ,η 2 ) has some small localized structures that are not compatible with statistical fluctuations. These structures are due to residual detector effects in the pair acceptance that are not removed by the eventmixing procedure, which can be important for extraction of the higher-order coefficients. Indeed, the Legendre coefficients for 064914-9 n 8 show significant nonstatistical fluctuations around zero. Therefore, the spread of a 2 n 1 /2 for n 10 and √ | a n a n+2 | for n 8 are quoted as uncertainties for the Legendre coefficients. These uncertainties are less than 0.5 × 10 −5 for a n a m calculated from C sub N (η 1 ,η 2 ) in all collision systems, and are larger for those calculated from C N (η 1 ,η 2 ). The corresponding relative uncertainty for a 2 1 is negligible. The HIJING and PYTHIA 8 events used for evaluating the reconstruction efficiency have a significant correlation signal and sizable a n,m coefficients for C N . The correlation functions obtained using the reconstructed tracks are compared with those obtained using the generated charged particles. The ratio of the two is then used to vary the measured C N (η 1 ,η 2 ), the procedure for removal of the SRC is repeated, and the variations of C sub N and a n,m are calculated. The differences in the correlation function reflect mainly the uncertainty in the efficiency correction, but also the influence of secondary decays and fake tracks. These differences are found to be mostly concentrated in a region around η − ≈ 0; hence, they affect mostly the estimation of the SRC component and have very little impact on C sub N and associated a n,m . The differences in Legendre coefficients are found to be up to 5% for a n calculated from C N and are 0.2-3.5% for a 2 1 1 /2 calculated from C sub N . The systematic uncertainties from the different sources described above are added in quadrature to give the total systematic uncertainties for the correlation functions and a 2 1 1 /2 values for both C N and C sub N . The systematic uncertainties associated with C sub N (η 1 ,η 2 ) and a 2 1 1 /2 are given in Tables I  and II, respectively. Since there are four methods for extracting a 2 1 1 /2 , they are given separately in Table II. The systematic uncertainty quoted for each source in both tables covers the maximum uncertainty in the specified collision system.

IV. RESULTS
The top row of Fig. 8 shows the correlation functions C N (η 1 ,η 2 ) in the three collision systems for events with similar multiplicity 100 N rec ch < 120. The corresponding estimated SRC component δ SRC (η 1 ,η 2 ) and long-range component C sub N (η 1 ,η 2 ) are shown in the middle and bottom rows, respectively. The magnitude of the SRC in p + Pb is observed to be larger in the proton-going direction than in the lead-going direction, reflecting the fact that the particle multiplicity is smaller in the proton-going direction. However, this forwardbackward asymmetry in p + Pb collisions is mainly associated TABLE II. Summary of systematic uncertainties for a 2 1 1/2 with p T > 0.2 GeV, calculated with four different methods: Legendre expansion of C sub N (η 1 ,η 2 ), quadratic fit of the η − dependence of C sub N (η 1 ,η 2 ) for |η + | < 0.1, quadratic fit of the η + dependence of C sub N (η 1 ,η 2 ) for 0.9 < |η − | < 1.  with the SRC component, and the C sub N (η 1 ,η 2 ) distribution shows very little asymmetry. The C N (η 1 ,η 2 ) distributions show significant differences between the three systems, which is mainly due to their differences in δ SRC (η 1 ,η 2 ). In fact, the estimated long-range component C sub N (η 1 ,η 2 ) shows similar shape and similar overall magnitude for the three systems.
To characterize the shape of the correlation functions, the Legendre coefficients a n a m for the distributions C N and C sub N shown in Fig. 8 are calculated via Eq. (6) and plotted in Fig. 9. The a n a m values are shown for the first six diagonal terms a 2 n and the first five mixed terms a n a n+2 , and they are also compared with coefficients calculated for opposite-charge pairs and same-charge pairs for the same event class. The magnitudes of the a n a m coefficients calculated for C N differ significantly for the different charge combinations, and they also increase as the size of the collision system decreases, i.e., | a n a m | p+p > | a n a m | p+Pb > | a n a m | Pb+Pb . This is consistent with a large contribution from SRC to all a n a m coefficients obtained from C N . After removal of the SRC, the a 2 1 coefficient is quite consistent between different charge combinations and different collision systems. All higher-order coefficients are much smaller, and they are very close to zero within the systematic uncertainties. Therefore, the rest of the paper focuses on the a 2 1 1 /2 results. To quantify further the shape of the LRC in C sub N (η 1 ,η 2 ), the a 2 1 1 /2 coefficients are also calculated by fitting the 1-D distributions from the three projection methods as outlined in Sec. III F: (1) quadratic fit of C sub N (η − ) in a narrow range of η + , (2) quadratic fit of C sub N (η + ) in a narrow range of η − , and in p + Pb collisions and pp collisions, respectively. Results are quite similar to those in Pb+Pb collisions, albeit with larger systematic uncertainties arising from the subtraction of a larger short-range component. For p + Pb (Fig. 11), the small FB asymmetry in the C sub N distribution along the η + direction is responsible for the difference in a 2 1 1 /2 between η + and −η + in the bottom left panel and between η ref and −η ref in the bottom right panel, but they still agree within their respective systematic uncertainties. Figure 13 shows a comparison of the a 2 1 1 /2 values extracted by the four methods as a function of N ch in the three collision systems. Good agreement between the different methods is observed.
On the other hand, the SRC is expected to have strong dependence on the charge combinations and collision systems, as shown by Figs. 8 and 9. The magnitude of the SRC is quantified by δ SRC (η 1 ,η 2 ) averaged over the two-particle pseudorapidity phase space: The corresponding contribution of the SRC at the singleparticle level is √ SRC , which can be directly compared with the strength of the LRC characterized by a 2 1 1 /2 . Figure 14 shows the values of √ SRC as a function of N ch for different charge combinations in the three collision systems. The strength of the SRC always decreases with N ch , and it is larger for smaller collision systems and opposite-charge pairs. Figure 15 compares the strength of the SRC in terms of √ SRC and the LRC in terms of a 2 1 1 /2 for the three collision systems. The values of √ SRC are observed to differ significantly while the values of a 2 1 1 /2 agree within ±10% between the three collision systems.
The strength of the SRC and LRC can be related to the number of clusters n contributing to the final multiplicity N ch , where n is the sum of clusters from the projectile and target nucleon or nucleus, n = n F + n B . The LRC is expected to be related to the asymmetry between n F and n B : The clusters could include the participating nucleons, subnucleonic degrees of freedom such as the fragmentation of scattered partons, or resonance decays. In an independent cluster model [37], of clusters, and hence, assuming n and N ch are proportional, the √ SRC and a 2 1 1 /2 values in Fig. 15 are expected to follow a simple power-law function in N ch : A power index that is less than one half, α < 0.5, would suggest that n grows more slowly than N rec ch , and vice versa. To test this idea, the √ SRC and a 2 1 1 /2 data in Fig. 15 are fit to a power-law function: c/N α ch . The function describes the N ch dependence in all three collision systems, with a reduced χ 2 values ranging between 0.2 and 0.9. The extracted power index values are summarized in Table III. The values of α for the SRC are found to be smaller for smaller collision systems, they are close to 0.5 in the Pb+Pb collisions and are significantly smaller than 0.5 in the pp collisions. In contrast, the values of α for a 2 1 1 /2 agree within uncertainties between the three systems and are slightly below 0.5.
One striking feature of the correlation function in p + Pb collisions, for example in Fig. 8  SRC, δ SRC (η 1 ,η 2 ) along the η + direction. Even in pp collisions, the δ SRC distribution is not uniform, but instead shows a quadratic increase towards large |η + | values. According to the discussion in Sec. III B, the shape of the δ SRC distribution in η + is described by the f (η + ) defined in Eq. (10). Examples of the f (η + ) are shown in Fig. 16 for p + Pb, symmetrized-p + Pb, pp, and Pb+Pb collisions with 100 N rec ch < 120. As described in Sec. III B, symmetrized-p + Pb results are obtained by averaging the proton-going and lead-going directions such that C(η 1 ,η 2 ) = C(−η 1 ,−η 2 ).
The independent cluster picture discussed above offers a simple interpretation of the shape of f (η + ). Assuming the population of clusters is a function of η, n c (η), and on average each cluster produces m charged particles according to a Poisson distribution, then the number of the SRC pairs scales as n c m(m − 1) = n c m 2 and the number of the combinatorial pairs scales as (n c m ) 2 . Therefore the strength of the SRC at given η is expected to scale as where n c (η) is assumed to be proportional to the local charge-particle multiplicity density dN ch / dη. Hence, the fact that f (η + ) is larger in the proton-going direction than in the Pb-going direction in p + Pb collisions simply reflects the asymmetric shape of the dN ch / dη distribution in each event [52]. The quadratic shape of f (η + ) for pp and symmetrizedp + Pb system therefore reflects a large, intrinsic FB asymmetry of dN ch /dη on an event-by-event level. The FB asymmetry in pp collisions is slightly larger than p + Pb collisions at comparable N ch , but is significantly less in Pb+Pb collisions. This observation suggests that the FB asymmetry for particle production in pp collisions could be as large as that in p + Pb collisions at comparable event activity, whereas the FB asymmetry for particle production is smaller in Pb+Pb collisions.

V. COMPARISON TO MODELS
QCD-inspired models such as PYTHIA and EPOS are often used to describe the particle production in pp collisions. ATLAS has previously compared the predictions of the PTYHIA8 A2 and EPOS LHC tunes with various single-particle distributions, such as the p T , η and the  event-by-event N ch distributions, fully unfolded for detector effects [33,34]. Reasonable agreement has been observed for these single-particle observables. In order to perform a data model comparison, the multiplicity correlation procedure used on the data is repeated for the two models to extract the SRC and LRC components. The extracted LRC in these models is then decomposed into Legendre coefficients of different order. The coefficients are found to be dominated by a 2 1 1 /2 , consistent with the observation that the shapes of the LRC are similar to those in the pp data in Fig. 8. However, the values of a 2 1 1 /2 predicted by the models are found to be much smaller than the pp data at the same N ch .
For a more direct comparison, Fig. 17 shows the N ch dependence of SRC and LRC from the data and the two models in pp collisions. The systematic uncertainties on the model predictions are dominated by the uncertainty in separating the SRC and LRC, as discussed in Sec. III G. However, at large N ch , they are also limited by the available MC statistics. There is some indication that the values of √ SRC from data are larger than the EPOS predictions and smaller than those from PYTHIA 8. Furthermore, the values from PYTHIA 8 increase for N ch > 120, a trend not supported by the data. On the other hand, both models underestimate significantly the values of a 2 1 1 /2 , suggesting that the FB multiplicity fluctuations in both models are significantly weaker than in the pp data. Therefore, these two models, which were tuned to describe many single-particle observables, fail to describe the longitudinal correlations between the produced charged particles.

VI. SUMMARY
Two-particle pseudorapidity correlations are measured with the ATLAS detector in √ s NN = 2.76 TeV Pb + Pb, √ s NN = 5.02 TeV p + Pb, and √ s = 13 TeV pp collisions at the LHC, with total integrated luminosities of approximately 7 μb −1 , 28 nb −1 , and 65 nb −1 , respectively. The correlation function C N (η 1 ,η 2 ) is measured using charged particles in the pseudorapidity range |η| < 2.4 with transverse momentum p T > 0.2 GeV, and it is measured as a function of event multiplicity N ch defined by the total number of charged particles with |η| < 2.5 and p T > 0.4 GeV. The correlation function shows an enhancement along the η 1 ≈ η 2 direction and suppression at η 1 ≈ −η 2 ∼ ±2.4, consistent with the expectation from an event-by-event forward-backward asymmetry in the multiplicity fluctuation (the long-range correlations or LRC). However, the correlation function also has a large narrow ridge along the η 1 ≈ η 2 direction associated with short-range correlations (SRC). The magnitudes of the SRC in p + Pb is found to be larger in the proton-going direction than the lead-going direction, reflecting the fact that the particle multiplicity is smaller in the proton-going direction. This is consistent with the observation that the SRC strength increases for smaller N ch . The SRC is observed to be much stronger for opposite-charge pairs than for the same-charge pairs, while the LRC is found to be similar for the two charge combinations. Based on this, a data-driven subtraction method was developed to separate the SRC and the LRC. The magnitudes of the SRC and the LRC are then compared for the three collision systems at similar values of N ch . After subtracting out the SRC δ SRC (η 1 ,η 2 ), the correlation function C sub N (η 1 ,η 2 ) is decomposed into a sum of products of Legendre polynomials that describe the different shape components, and the coefficients a n a m are calculated. Significant values are observed for a 2 1 in all N ch ranges and higher-order coefficients are consistent with zero, and suggesting that C sub N has an approximate functional form C sub N ≈ 1 + a 2 1 η 1 η 2 . The quantity a 2 1 is also estimated by parametrization of the shape of the correlation function in narrow ranges of η − = η 1 − η 2 and η + = η 1 + η 2 , or from a ratio C sub N (η 1 ,η 2 )/C sub N (−η 1 ,η 2 ), and consistent results are obtained. The magnitude of the SRC and a 2 1 1 /2 are compared for the three collision systems as a function of N ch . Large differences are observed for the SRC, but the values of a 2 1 1 /2 agree within ±10% at the same N ch . The N ch dependences of both the SRC and a 2 is approximately the same for the three collision systems. In contrast, the power-law index for the SRC is smaller for smaller collision systems. The SRC distribution shows strong dependence on η + in p + Pb and pp, but much weaker dependence in Pb+Pb collisions. The δ SRC (η + ) distribution, after symmetrizing the proton and lead directions, is found to be similar to the SRC in pp collisions with comparable N ch , suggesting that the event-by-event FB asymmetry for particle production is similar in pp and p + Pb collisions with comparable event activity. The PYTHIA 8 A2 and EPOS LHC models, which were tuned to describe many single-particle observables in pp collisions, fail to describe the SRC and the LRC observed in the pp data.