Femtoscopy with identified charged pions in proton-lead collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with ATLAS

Bose-Einstein correlations between identified charged pions are measured for $p$+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV using data recorded by the ATLAS detector at the LHC corresponding to a total integrated luminosity of $28$ $\mathrm{nb}^{-1}$. Pions are identified using ionization energy loss measured in the pixel detector. Two-particle correlation functions and the extracted source radii are presented as a function of collision centrality as well as the average transverse momentum ($k_{\mathrm{T}}$) and rapidity ($y^{\star}_{\pi\pi}$) of the pair. Pairs are selected with a rapidity $-2<y^{\star}_{\pi\pi}<1$ and with an average transverse momentum $0.1<k_{\mathrm{T}}<0.8$ GeV. The effect of jet fragmentation on the two-particle correlation function is studied, and a method using opposite-charge pair data to constrain its contributions to the measured correlations is described. The measured source sizes are substantially larger in more central collisions and are observed to decrease with increasing pair $k_{\mathrm{T}}$. A correlation of the radii with the local charged-particle density is demonstrated. The scaling of the extracted radii with the mean number of participating nucleons is also used to compare a selection of initial-geometry models. The cross-term $R_\mathrm{ol}$ is measured as a function of rapidity, and a nonzero value is observed with $5.1\sigma$ combined significance for $-1<y^{\star}_{\pi\pi}<1$ in the most central events.

to a total integrated luminosity of 28 nb −1 . Pions are identified using ionization energy loss measured in the pixel detector. Two-particle correlation functions and the extracted source radii are presented as a function of collision centrality as well as the average transverse momentum (k T ) and rapidity (y ππ ) of the pair. Pairs are selected with a rapidity −2 < y ππ < 1 and with an average transverse momentum 0.1 < k T < 0.8 GeV. The effect of jet fragmentation on the two-particle correlation function is studied, and a method using oppositecharge pair data to constrain its contributions to the measured correlations is described. The measured source sizes are substantially larger in more central collisions and are observed to decrease with increasing pair k T . A correlation of the radii with the local chargedparticle density is demonstrated. The scaling of the extracted radii with the mean number of participating nucleons is also used to compare a selection of initial-geometry models. The cross-term R ol is measured as a function of rapidity, and a nonzero value is observed with 5.1σ combined significance for −1 < y ππ < 1 in the most central events.

Introduction
Studies of multi-particle correlations in proton-lead (p+Pb) [1][2][3][4][5] and proton-proton (pp) [6] collisions at the LHC and in deuteron-gold (d+Au) [7-9] and helium-3-gold ( 3 He+Au) [10] collisions at RHIC have shown that these correlation functions exhibit features similar to those observed in nucleus-nucleus collisions [11][12][13][14][15][16] that are attributed to collective dynamics of the strongly coupled quark-gluon plasma. In particular, two-particle angular correlations studied in high multiplicity p+Pb [1,4,17] and pp [6,18] collisions at the LHC show a "ridge" -an enhancement in the correlation function at small relative azimuthal angle (∆φ) that extends over a range of relative pseudorapidity (∆η). The ridge in both systems is generally understood to result from a combination of sinusoidal modulations of the two-particle correlation function of different harmonics [1][2][3][4][5]. Hydrodynamic calculations show that such a modulation can arise from initial-state spatial anisotropies that, through the collective expansion of the medium, are imprinted on the azimuthal-angle distributions of the final-state particles (see e.g. Refs. [19][20][21] and references therein). Such hydrodynamic models can also reproduce the modulation observed in p+Pb, d+Au, and 3 He+Au collisions [22][23][24][25], but the suitability of hydrodynamics in "small" systems remains a topic of active debate. Alternatively, the observed modulation in these collisions has been explained by so-called "glasma" models that invoke saturation of the nuclear parton distributions [26][27][28][29][30]. To disentangle these competing explanations and, more specifically, to test whether collective phenomena are present in p+Pb collisions at the LHC, additional measurements are required to constrain the source geometry.
Hanbury Brown and Twiss (HBT) correlations, which probe the space-time extent of a particle-emitting source (see Ref. [31] and references therein), may provide valuable insight into the problems described above. The HBT method originated in astronomy [32,33], where space-time correlations of photons due to wave function symmetrization are used to measure the size of distant stars. The procedure can be adapted to the extremely small sources encountered in hadronic collisions if identical-particle Bose-Einstein correlations are instead studied in relative momentum space [34]. The two-particle correlation function C(q), parameterized as a function of relative momentum, is sensitive to the two-particle source density function S(r) through the two-particle final-state wavefunction [31]: where q and k are, respectively, the relative and average momentum of a pair of particles, r is the distance between the origin points of the two particles, and the two-particle source function S k is normalized so that ∫ d 3 r S k (r) = 1. In the case of a non-interacting identical boson wave-function, the term within the parentheses of Eq. 1 is a cosine and the correlation function is enhanced by the Fourier transform of the source function. Thus, the Bose-Einstein modification of the relative momentum distributions produces an enhancement at small q whose range in q is inversely related to the size of the source.
In a typical HBT analysis, the correlation functions are fit to a function of relative momentum that is often a Gaussian or exponential function, or a stretched exponential function that can interpolate between these two. The parameters of the fits that relate to the space-time extent of the source function are referred to as the "HBT radii". Measurements of Bose-Einstein correlations in pp collisions at center-of-mass energies It is also observed that the extracted radii increase with particle multiplicity but saturate at the highest multiplicities.
Although Bose-Einstein correlations are the most straightforward to measure experimentally, any finalstate interaction can in principle be used to image the source density. The term "femtoscopy" is often used to refer to any measurement that provides spatio-temporal information about a hadronic source [38]. The measured source radii are interpreted as the dimensions of the region of homogeneity of the source at freeze-out, after all interactions between final-state particles and the bulk have ceased; thus, they are sensitive to the space-time evolution of the event. In particular, an increase in radii at low average transverse momentum k T indicates radial expansion since higher-momentum particles are more likely to be produced earlier in the event [39]. The k T -scaling of HBT radii in p+Pb systems is of significant interest when studied as a function of centrality, an experimental proxy for the impact parameter. Thus, these measurements can provide insight into the conditions necessary for hydrodynamic behavior in small systems.
In many HBT measurements, the correlation functions are evaluated in one dimension using the invariant relative momentum q inv ≡ −q µ q µ , where q = p a − p b for a pair of particles a and b with four-momenta p a and p b . In three dimensions, HBT correlations are studied using the "out-side-long" convention [40][41][42][43].
In this system, q out , the outwards component, is the projection along k T ; q side , the sideways component, is the projection alongẑ × k T (with the z-axis along the beamline); and q long is the longitudinal component. The relative momentum of the pair is evaluated in the longitudinally comoving frame (LCMF), i.e., the frame boosted such that k z = 0. This formulation of the HBT analysis has the advantage that it decomposes the correlation function into components that emphasize distinct physical effects. In particular, the spatial extent of the source in the longitudinal and transverse directions is likely to be different. The out and side radii are also expected to differ due to the effects of the Lorentz boost in the out direction and, if the system exhibits collectivity, due to space-momentum correlations. In a fully boost-invariant system, observables evaluated in the LCMF should be independent of k z (or rapidity). The inherent asymmetry of p+Pb collisions seen, for example, in the charged-particle pseudorapidity distributions [44,45], provides a unique opportunity to study the correlations between source sizes and the pair's rapidity, collision centrality, or the local (in rapidity) charged-particle density. The results of such a study may provide insight into or constrain theoretical models of the underlying dynamics responsible for producing the final-state particles.
To address the topics and questions discussed above, this paper presents measurements of correlations between identified charged pions in 5.02 TeV p+Pb collisions which were performed by the ATLAS experiment at the LHC. While femtoscopic methods have already been applied to p+Pb systems at the LHC [46,47], this paper presents a new data-driven technique to constrain the significant background contribution from jet fragmentation, referred to in this paper as the "hard process" background. It also provides new measurements of the dependence of the source radii on the pair's rapidity y ππ , calculated assuming both particles have the mass of the pion, over the range −2 < y ππ < 1. Results are presented for one-and three-dimensional source radii as a function of the pair's average transverse momentum, k T , over the range 0.1 < k T < 0.8 GeV and for several p+Pb centrality intervals with the most central case being 0-1%. The p+Pb collision centrality is characterized using ΣE Pb T , the total transverse energy measured in the Pb-going forward calorimeter (FCal) [45]. It is defined such that central events, with large ΣE Pb T , have a low centrality percentage, and peripheral events, with a small ΣE Pb T , have a high centrality percentage. Using the measured centrality dependence of the source radii, the scaling of the system size with the number of nucleon participants N part is also investigated, using a generalization of the Glauber model [48].

ATLAS detector
The ATLAS detector is described in detail in Ref. [49]. The measurements presented in this paper have been performed using the inner detector, minimum-bias trigger scintillators (MBTS), FCal, zerodegree calorimeter (ZDC), and the trigger and data acquisition systems. The inner detector [50], which is immersed in a 2 T axial magnetic field, is used to reconstruct charged particles within |η| < 2.5.1 It consists of a silicon pixel detector, a semiconductor tracker (SCT) made of double-sided silicon microstrips, and a transition radiation tracker made of straw tubes. All three detectors consist of a barrel and two symmetrically placed endcap sections. A particle traveling from the interaction point (IP) with |η| < 2 crosses at least 3 pixel layers, 4 double-sided microstrip layers and typically 36 straw tubes. In addition to hit information, the pixel detector provides time-over-threshold for each hit pixel which is proportional 1 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).
to the deposited energy and which is used to provide measurements of specific energy loss (dE/dx) for particle identification.
The FCal covers a pseudorapidity region of 3.1 < |η| < 4.9 and is used to estimate the centrality of each collision. The FCal uses liquid argon as the active medium with tungsten and copper absorbers. The MBTS, consisting of two arrays of scintillation counters, are positioned at z = ±3.6 m and cover 2.1 < |η| < 3.9. The ZDCs, situated approximately 140 m from the nominal IP, detect neutral particles, mostly neutrons and photons, that have |η| > 8. 3. They are used to distinguish pileup events (bunch crossings involving more than one collision) from central collisions by detecting spectator nucleons that did not participate in the interaction. The calorimeters use tungsten plates as absorbers and quartz rods sandwiched between the tungsten plates as the active medium.
Events used for the analysis presented in this paper were primarily obtained from a combination of minimum-bias (MinBias) triggers that required either at least two hit scintillators in the MBTS or at least one hit on each side of the MBTS. An additional requirement on the number of hits in the SCT was imposed on both of these minimum-bias triggers to remove false triggers. To increase the number of events available in the highest ΣE Pb T interval, the analysis includes a separate sample of events selected by a trigger (HighET) that required a total transverse energy in both sides of the FCal of at least 65 GeV.

LHC data
This analysis uses data from the LHC 2013 p+Pb run at √ s NN = 5.02 TeV with an integrated luminosity of 28 nb −1 . The Pb ions had an energy per nucleon of 1.57 TeV and collided with the 4 TeV proton beam to yield a center-of-mass energy √ s NN = 5.02 TeV with a longitudinal boost of y CM = 0.465 in the proton direction relative to the ATLAS laboratory frame. The p+Pb run was divided into two periods between which the directions of the proton and lead beams were reversed. The data in this paper are presented using the convention that the proton beam travels in the forward (+z) direction and the lead beam travels in the backward (−z) direction. When the data from these two periods are combined, the MinBias triggers sampled a total luminosity, after prescale, of 24.5 µb −1 and yielded a total of 44 million events; the HighET trigger sampled a total luminosity of 41.4 µb −1 after prescale and yielded 700 thousand events.

Monte Carlo event generators
The effects of charged-particle reconstruction and selection are studied in a p+Pb sample generated using H [51] and simulated with the GEANT4 package [52]. Five million events are generated at a centerof-mass energy per nucleon-nucleon pair of √ s NN = 5.02 TeVwith a longitudinal boost of y CM = 0.465 in the proton direction. The sample is fully reconstructed with the same conditions as the data [53].
Four additional Monte Carlo generator samples are used to study the background from hard processes, as described in Section 4.2. No detector simulation is performed on these samples, as the net effects of the simulation and reconstruction were studied using the fully reconstructed p+Pb simulation events and found to be negligible. The two-particle reconstruction effects occur only at very low q (as discussed in Section 5.1), but these generated samples are used only to study correlations from jet fragmentation which span a much broader range of q. In each of the following samples, 50 million (250 million for

Event selection and centrality
In the offline analysis, charged-particle tracks and collision vertices are reconstructed using the same algorithms and methods applied in previous minimum-bias pp and p+Pb measurements [45,60]. Events included in this analysis are required to pass either of the two MinBias triggers or the HighET trigger, to have a hit on each side of the MBTS with a difference in average particle arrival times measured on the two sides of the MBTS which is less than 10 ns, a reconstructed primary vertex (PV), and at least two tracks satisfying the selection criteria listed in Section 3.4. Events that have more than one reconstructed vertex (including secondary vertices) with either more than ten tracks or a sum of track transverse momentum (p T ) greater than 6 GeV are rejected. An upper limit is placed on the activity measured in the Pb-going ZDC to further reject pileup events.    [45] for each centrality interval in the Glauber model as well as the two choices for the Glauber-Gribov model with color fluctuations (GGCF) [61] (and references therein), along with the average multiplicity with p T > 100MeV and |η| < 1.5 also obtained from Ref. [45]. The parameter ω σ represents the size of fluctuations in the nucleon-nucleon cross section. Asymmetric systematic uncertainties are shown for N part . The uncertainties in dN ch /dη are given in the order of statistical followed by systematic.
The centralities of the p+Pb events are characterized following the procedures described in Ref. [45], using ΣE Pb T , the total transverse energy in the Pb-going side of the FCal. The use of the FCal for measuring centrality has the advantage that it is not sensitive to multiplicity fluctuations in the kinematic region covered by the inner detector, where the measurements are performed. Measurements are presented in this paper for the centrality intervals listed in Table 1. The events selected using the HighET trigger are used only in the 0-1% centrality interval. Figure 1 shows the distribution of ΣE Pb T values obtained from events included in this measurement. The discontinuity in the spectrum occurs at the low edge of the 0-1% centrality interval, above which the HighET events are included.
For each centrality interval, the average multiplicity of charged particles with p T > 100 MeV and |η| < 1.5, dN ch /dη , and the corresponding average number of participating nucleons, N part , are obtained from a previous publication [45]. Since this analysis uses finer centrality intervals (no wider than 10% of the total centrality range) than those used in Ref. [45], a linear interpolation over the Glauber N part is used to construct additional values for dN ch /dη based on the published results. This interpolation is justified by the result in Ref. [45] that charged-particle multiplicity is proportional to N part in the peripheral region. The values and uncertainties from this procedure are listed in Table 1.

Charged-particle selection and pion identification
Reconstructed tracks used in the HBT analysis are required to have |η| < 2.5 and p T > 0.1 GeV and to satisfy a standard set of selection criteria [60]: a minimum of one pixel hit is required, and if the track crosses an active module in the innermost layer, a hit in that layer is required; for a track with p T greater than 0.1, 0.2, or 0.3 GeV there must be at least two, four, or six hits respectively in the SCT; the transverse impact parameter with respect to the primary vertex, d PV 0 , must be such that |d PV 0 | < 1.5 mm; and the p+Pb events, as a function of the pair's average transverse momentum k T and rapidity y ππ . The rapidity is calculated using the pion mass for both reconstructed charged particles. corresponding longitudinal impact parameter must satisfy |z PV 0 sin θ| < 1.5 mm. To reduce contributions from secondary decays, a stronger constraint on the pointing of the track to the primary vertex is applied. Namely, neither |d PV 0 | nor |z PV 0 sin θ| can be larger than three times its uncertainty as derived from the covariance matrix of the track fit.
Particle identification (PID) is performed through measurements of the specific energy loss dE/dx derived from the ionization charge deposited in the pixel clusters associated with a track. The dE/dx of a track is calculated as a truncated mean of the dE/dx in individual pixel clusters as described in Ref. [62], since the truncated mean gives a better resolution than the mean. Relative likelihoods that the track is a π, K, and p are formed by fitting the dE/dx distributions to √ s = 7 TeV pp data in several momentum intervals as explained in Ref. [63]. Three PID selection levels are defined: one designed to have a high efficiency for pions, one designed to result in high pion purity, and one in between that was chosen as the nominal selection level and is used throughout the analysis if other PID selections are not explicitly mentioned. The efficiency and purity of these selections are studied in the fully reconstructed simulated sample. The resulting purity of track pairs in the nominal selection is shown in Figure 2 as a function of pair's k T and y ππ . The results are also evaluated at the looser and tighter PID definitions (also in Figure 2), and the differences are incorporated into the systematic uncertainty (see Section 5).

Pair selection
Track pairs are required to have |∆φ| < π/2 to avoid an enhancement in the correlation function arising primarily from dijets. This enhancement does not directly affect the signal region but can influence the results by affecting the overall normalization factor in the fits. The pair's rapidity y ππ , measured with respect to the nucleon-nucleon center of mass, must lie in the range −2 < y ππ < 1. This requirement is more stringent than the single-track requirement |η| < 2.5. When analyzing track pairs of opposite charge, common particle resonances are removed via requirements on the invariant mass so that m ππ − m ρ 0 > The q inv , q long , |q side |, and q out distributions of the pairs obtained through these procedures are shown in Figure 3 for the 0-1% and 60-80% centrality intervals. The one-dimensional q inv distribution necessarily decreases to zero at q inv = 0 due to the scaling of the phase-space volume element d 3 q ∝ q 2 dq. In contrast the three-dimensional quantities remain finite at zero relative momentum. The distributions are nearly identical for the two centrality intervals, although differences can be seen at small relative momentum in all four distributions.   Figure 3: Pair-normalized distributions of the invariant relative momentum q inv (a) and the three-dimensional relative momentum components q long (b), q side (c) and q out (d) for identified same-charge pion pairs with 0.2 < k T < 0.3 GeV obtained from two centrality intervals: 60-80% and 0-1%.

Correlation function analysis
The two-particle correlation function is defined as the ratio of two-particle to single-particle momentum spectra: for pairs of particles with four-momenta p a and p b . This definition has the useful feature that most single-particle efficiency, acceptance, and resolution effects cancel in the ratio. The correlation function is expressed as a function of the relative momentum2 q ≡ p a − p b in intervals of average momentum k ≡ p a + p b /2.
The relative momentum distribution A(q) ≡ dN/dq| same (Figure 3) is formed by selecting like-charge (or unlike-charge) pairs of particles from each event in an event class, which is defined by the collision centrality and z position of the primary vertex (z PV ). The combinatorial background B(q) ≡ dN/dq| mix is constructed by event mixing, that is, by selecting one particle from each of two events in the same event class as A(q). Each particle in the background fulfills the same selection requirements as those used in the same-event distribution. Event classes are categorized by centrality so that events are only compared to others with similar multiplicities and momentum distributions. Events are sorted by z PV so that the background distribution is constructed with pairs of tracks originating from nearby space points, which is necessary for B(q) to accurately represent the as-installed detector. The A(q) and B(q) distributions are combined over z PV intervals in such a way that each of them samples the same z PV distribution. The ratio of the distributions defines the correlation function:

Parameterization of the correlation function
Assuming that all particles are identical pions created in a fully chaotic source and that they have no final-state interaction, the enhancement in the correlation function is the Fourier transform of the source density.
The Bowler-Sinyukov formalism [64,65] is used to account for final-state corrections: where K is a correction factor for final-state Coulomb interactions, and C BE (q) = 1 + F [S k ] (q) with F [S k ] (q) denoting the Fourier transform of the two-particle source density function S k (r). Several factors influence the value of the parameter λ. Including non-identical particles decreases this parameter, as does coherent particle emission. Products of weak decays or long-lived resonances also lead to a decrease in λ, as they are emitted at a length scale greater than can be resolved by femtoscopic methods given the momentum resolution of the detector. These additional contributions to the source density are not Coulomb-corrected within the Bowler-Sinyukov formalism. When describing pion pairs of opposite charge, there is no Bose-Einstein enhancement and C BE → 1.
Coulomb interactions suppress the correlation at small relative momentum for identical charged particles. The particular choice of correction factor K(q) is determined using the formalism in Ref. [66]. This uses the approximation that the Coulomb correction is effectively applied not from a point source, but over a Gaussian source density of radius R eff : where a = 388 fm is the Bohr radius [67] of a two-pion state, 2 F 2 is a generalized hypergeometric function, and G(q inv ) is the Gamow factor [68,69] For opposite-charge pairs, a is taken to be negative, since its definition includes a product of the two charges.
The Bose-Einstein enhancement in the invariant correlation functions is fit to an exponential form: where R inv is the Lorentz-invariant HBT radius. This function corresponds to an underlying Breit-Wigner source density.
The Bose-Einstein component of the three-dimensional correlation functions is fit to a function of the form where R is a symmetric matrix of the form The off-diagonal entries other than R ol can be argued to vanish by the average azimuthal symmetry of the source. In hydrodynamic models the out-long term R ol is sensitive to spatio-temporal correlations, and therefore, to the lifetime of the source [70,71]. It couples radial and transverse expansion, and is expected to vanish in the absence of either. If the source is fully boost invariant then this term vanishes, so an observation of a nonzero value demonstrates that the homogeneity region is not boost invariant.
In order to reduce computational demands, a few symmetry arguments are considered. The order of the pairs is chosen such that q out is always positive, which can be done so long as C(−q) = C(q). The average azimuthal symmetry of the source is invoked in order to allow only the absolute value of q side to be considered. The sign of q long cannot be similarly discarded if a nonzero R ol is allowed.
A Gaussian form for the Bose-Einstein enhancement is often used in the three-dimensional correlation function. However, this form was found to give a poor description of ATLAS data, relative to an exponential form. This was also observed in the ATLAS pp results in Ref. [35]. The chosen form of F [S k ](q) must be taken into account when interpreting source radii, and there is no simple correspondence between parameters estimated using one form and those from another. An ad hoc factor of √ π is often invoked to relate Gaussian radii to exponential radii by assuming that the first q-moment of the invariant correlation function should be preserved, but this assumption is not rigorously justified and the argument fails in general for three-dimensional correlation functions.

Hard-process contribution
Additional non-femtoscopic enhancements to the correlation functions at q inv 0.5 − 1 GeV are observed in both the opposite-charge (+−) and the same-charge (±±) pairs. As discussed later in this section, the enhancement is more prominent in +− pairs than it is in ±± pairs. Monte Carlo (MC) generators do not simulate the final-state two-particle interactions used for femtoscopy, but they do describe the background correlations. The MC generators used in this section to constrain the background description are described in Section 3.2.
The non-femtoscopic enhancement is more prominent for higher k T and lower multiplicities. This suggests that the correlation is primarily due to jet fragmentation. This hypothesis is verified by studying correlation functions in H , by increasing the minimum hard-scattering p T (p HS,min T ) from 2 to 20 GeV.
Increasing p HS,min T has the effect of suppressing most hard processes in typical events. Without the resulting jet fragmentation, the non-femtoscopic enhancement is removed from the correlation function, as demonstrated in Figure 4 by comparing the panels on the left and right of each row.
The amplitude of the hard-process contribution tends to be larger in the Monte Carlo events than it is in the data. Thus, attempting to account for it by studying the double-ratio C data (q)/C MC (q) leads to a depletion that is apparent in the region where the Bose-Einstein enhancement disappears [35]. Another commonly used method is to parameterize the mini-jet contribution using simulation and to allow one or more parameters of the description to vary in the fit [46,47].
To avoid too much reliance on either a full MC description or arbitrary additional free parameters, a data-driven method is derived here to constrain the correlations from jet fragmentation. Opposite-charge correlation functions are used to predict the jet contribution in the same-charge correlation function. This poses two challenges. Firstly, resonance decays appear prominently in the opposite-charge correlations.   Figure 4: Correlation functions of charged particles from H for opposite-(a,b) and same-charge (c,d) pairs with transverse momentum 0.7 < k T < 0.8 GeV, using events with a generated multiplicity 26 ≤ N tr ch ≤ 36. The generator is run with the minimum hard-scattering p HS, min T at the default setting of 2 GeV (a,c) and increased to 20 GeV (b,d) to remove the contribution from hard processes. The gaps in the opposite-charge correlation functions are a result of the requirements described in Section 3.5, which remove the largest resonance contributions.
The most prominent of these are removed by requirements on the invariant mass of the opposite-charge pairs (as described in Section 3.5), and the fits to the opposite-charge correlation functions are restricted to q inv > 0.1 GeV. The lower bound on the domain of the fit reduces sensitivity to effects such as three-body decays that are unrelated to jet fragmentation, which is significant over a broader range of q. Secondly, jet fragmentation does not affect opposite-charge and same-charge correlations in an identical manner. This is in part because opposite-charge pairs are more likely to have a closer common ancestor in a jet's fragmentation into hadrons.
To account for the remaining differences between +− and ±± pairs, a study of both classes of correlation functions is performed in P 8. In order to isolate the effect of jet fragmentation, decays from the relatively longer-lived particles η, η , and ω are excluded. Pairs of particles from two-body resonance decays are also neglected, in order to remove mass peaks in the correlation function. The same-pair mass cut around the ρ resonance that is used in the data is also applied in P 8 events, since the removal of the corresponding region of phase space has a significant effect on the shape of the correlation function.

Jet fragmentation in q inv
To describe the jet fragmentation in the invariant correlation functions, fits are performed in P 8 pp to a stretched exponential function of the form where N is a normalization factor and the other parameters depend on the charge combination and on k T .
The Ω(q inv ) function above is applied as a multiplicative factor to the femtoscopic correlation function. The strategy employed is to estimate these parameters for same-charge correlation functions based on values determined using opposite-charge correlations. First, the shape parameter α inv bkgd is determined with fits to same-charge correlation functions, with all parameters allowed to be free. It is only weakly dependent on multiplicity, so a function is fit to parameterize α inv bkgd in P 8 as a function of k T (with k T in GeV): α inv bkgd (k T ) = 2 − 0.050 ln 1 + e 50.9(k T −0.49) .
The fits are well described by a Gaussian form (α inv bkgd = 2) at k T 0.4 GeV, and α inv bkgd decreases to a value around 1.3 in the highest k T interval.
The fits are performed again to the P 8 correlation functions, with α inv bkgd now fixed to the same value in same-and opposite-charged pairs, and a comparison is made between the width parameters R inv bkgd +− and R inv bkgd ±± . The width of the jet fragmentation correlation for same-charge pairs is found to be correlated to that for opposite-charge pairs, as shown in the right plot of Figure 5. Four intervals of charged particle multiplicity, N ch , calculated for particles with p T > 100 MeV and |η| < 2.5 are shown: 26 ≤ N ch ≤ 36, 37 ≤ N ch ≤ 48, 49 ≤ N ch ≤ 64, and 65 ≤ N ch . The relationship between the invariant background widths is modeled as a direct proportionality, with a value of ρ = 1.3 extracted from P 8. This proportionality begins to break down at low k T , but the model becomes increasingly accurate at larger k T , where hard processes give a larger contribution to the correlation function.  Next, R inv bkgd ±± is fixed from R inv bkgd +− using the value of ρ, and the fits are performed again to parameterize the relationship between the amplitudes, As shown in the left-hand plot of Figure 5, µ and ν are fit in each k T interval to describe four multiplicity intervals. The power-law scaling of Eq. 9 is found to provide a good description of the relation between the same-and opposite-charge amplitudes across all four multiplicity intervals studied. The multiplicityindependence of µ and ν is important in justifying the use of these parameters in p+Pb.
The correspondence between opposite-and same-charge pairs in both pp and p+Pb systems is studied in H , since the study described in this section is performed with P 8 in a pp system. While the mapping is mostly consistent between the two systems, it is found that µ is larger in p+Pb than in pp by 8.5% on average. When the mapping is applied to the data, this attenuation factor (along with a corresponding systematic uncertainty described in Section 5) is also taken into account.
With α inv bkgd (k T ), µ(k T ), ν(k T ), and ρ determined from Monte Carlo generator samples, the mapping can be applied to the p+Pb data. As illustrated in Figure 6, the +− correlation function is fit to Eq. 7 for q inv > 0.1 GeV, with α bkgd fixed from P 8 and λ inv bkgd +− and R inv bkgd +− as free parameters. The µ, ν, and ρ parameters are used to infer λ inv bkgd ±± and R inv bkgd ±± , which are fixed before the femtoscopic part of the correlation function is fit to ±± data.

Jet fragmentation in three dimensions
In the longitudinally comoving frame of a particle pair produced in a jet, the axis of the jet is aligned on average with the "out" direction and the plane transverse to the jet's momentum is spanned by the Figure 6: Correlation functions in p+Pb data for opposite-charge (teal circles) and same-charge (red squares) pairs. The opposite-charge correlation function, with the most prominent resonances removed, is fit to a function of the form in Eq. 7 (blue dashed line). The violet dotted line is the estimated jet contribution in the same-charge correlation function, also of the form of Eq. 7, and the dark red line is the full fit of Eq. 19 to the same-charge data.
"side" and "long" directions. In three dimensions the correlation from jet fragmentation is factorized into components which separately describe the "out" direction and both the "side" and "long" directions: where q sl = q 2 side + q 2 long , λ osl bkgd is the background amplitude, and α out bkgd and α sl bkgd parameterize the shape of the fragmentation contribution along and transverse to the jet axis, respectively. The shape parameters α out bkgd and α sl bkgd are taken from P 8 and fixed to 1.5 and 1.7 respectively. Fits of these parameters to P 8 correlation functions are not fully consistent with these chosen numerical constants at all k T and multiplicities. However, the impact of the somewhat arbitrary choice of fixing these parameters is tested by varying them both by 0.1, and the changes in the results are less than 1%.
Similarly to the procedure used for the q inv correlation functions, the width parameters are compared between opposite-and same-charge correlation functions (bottom panels of Figure 7); however, in three dimensions the relationships are parameterized as a function of k T . Next, as for q inv , the amplitudes for three-dimensional jet correlations are compared between opposite-and same-charge pairs (top panel of Figure 7). While the relationships between opposite-and same-charge correlations are not well described everywhere by the fitted lines, the model becomes increasingly accurate at larger k T , where hard processes give a larger contribution to the correlation function. The functional forms of the mappings from oppositeto same-charge three-dimensional parameters are: The widths used in the three-dimensional background description are related by a k T -dependent additive factor [ Eqs. 12 and 13 ] because a simple proportionality [ Eq. 8 ] is not as successful in describing the behavior.  The invariant and three-dimensional (3D) fragmentation amplitudes are strongly correlated and the mappings from opposite-to same-charge correlation functions are quantitatively similar. Thus, the same 8.5% attenuation factor for µ derived from H for the invariant mapping is used for the three-dimensional fits as well.
The numerical values used for mapping the amplitude λ osl bkgd , fragmentation width colinear with the jet axis R out bkgd , and fragmentation width transverse to the jet axis R sl bkgd are given by the following parameterizations (k T in GeV and R bkgd in GeV −1 ): The jet fragmentation parameters of the p+Pb data depend on centrality and k T . The same-charge amplitude of the background ranges from being negligible at low k T up to a maximum of roughly 0.25 at the largest measured k T of 0.8 GeV for the most peripheral events. The widths of the same-charge fragmentation correlation have length scales which are typically in a range of 0.3-0.5 fm at the largest k T where they are most relevant. The p+Pb femtoscopic measurement is most challenging at high k T in peripheral events, where the fragmentation background amplitude is a significant fraction of the Bose-Einstein amplitude and the HBT radii are smaller and closer in magnitude to the length scale of the jet correlation.

Fitting procedure
The bin contents of the histogram representations of A(q) and B(q) are assumed to be Poisson distributed.
The correlation function C(q) is assumed to be fit best by the ratio of their means, and a flat Bayesian prior is assumed for both means. A corresponding χ 2 analog [72], the negative log-likelihood ratio L [73], is minimized using the Minuit package [74]: Here A and B are the signal and background relative momentum distributions in Eq. 2 when represented as histograms, such that A i and B i are the contents in bin i, C i is shorthand for C(q i ) where q i is the bin center, and C(q) is the fitting function describing the correlation. The multiplicative factor of −2 causes this statistic to approach χ 2 as the sample size increases. The 1σ statistical uncertainties in the fit parameters of C(q) are evaluated using the MINOS routine, and are selected from the points in the parameter space where −2 ln L = min (−2 ln L) + 1.
The full form of the invariant-correlation-function fit to like-charge track pair data including the hardprocess background description is where C BE (q) is given by Eq. 4 or Eq. 5, K(q inv ) is given by Eq. 3, and Ω(q) is given by Eq. 7 or Eq. 10.
As discussed in Section 4.2, the opposite-charge correlation functions are fit in the regions where q inv (or |q| in 3D) is greater than 100 MeV. The opposite-charge parameters are highly insensitive to the choice of cutoff, as the q distributions contribute more statistical weight at larger q. The same-charge correlation functions are fit in the regions q inv > 30 MeV for the invariant fits and |q| > 25 MeV in three dimensions.

Sources of systematic uncertainty
The systematic uncertainties in the extracted parameter values have contributions from several sources: the jet fragmentation description, PID, the effective Coulomb-correction size R eff , charge asymmetry, and particle reconstruction effects.
One of the largest sources of uncertainty originates from the description of the background correlations Ω(q) from jet fragmentation. For the uncertainty in the hard-process contribution, three effects are considered. First, the extrapolation from a pp to a p+Pb system is represented with an uncertainty in the background amplitude. Also, to investigate the uncertainty in the Monte Carlo description of jet fragmentation, the amplitude of C +− (q inv ) C ±± (q inv ) is studied in both P and Herwig. Herwig does not predict enough difference between +− and ±± correlations to describe the data. Thus, instead of using the ratio of the predicted scalings of the two generators, the standard deviation of the ratio amplitude (across a selection of k T and multiplicity intervals) is used as a variation reflecting this systematic uncertainty. The hard-process amplitude λ bkgd is scaled up and down by 12.3%, the quadrature sum of the relative variation from the difference between the pp and p+Pb systems (4.1%) and from the generator difference (11.6%). The widths of the background description are highly correlated with the amplitude in the P fit results, so varying the widths in addition to the amplitude would overstate the uncertainty. The choice of varying the amplitude instead of the width is found to provide a larger and more consistent variation in the radii, so only the amplitude of the background is varied. The variation from the combination of the generator and the collision system are indicated by a label of "Gen ⊕ Sys" in the figures of Section 5.2. Additionally, the procedure described in Section 4.2 to control the jet fragmentation correlations is repeated in both the central (|y ππ | < 1) and forward (−2 < y ππ < −1) rapidity intervals. While the relationship between the fragmentation widths (R inv bkgd , R sl bkgd , and R out bkgd ) is fairly robust, the mappings of the amplitudes (λ inv bkgd and λ osl bkgd ) from opposite-to same-charge correlations vary between the two rapidity intervals. This variation represents the breakdown of the assumptions used to describe the jet fragmentation. The mapping procedure is repeated with the results from each rapidity interval, and the variation is used as an additional systematic uncertainty in the amplitude. The HBT radii and amplitudes are both highly correlated with the amplitude of the jet background, so varying λ bkgd is a robust method of evaluating the uncertainties from the background description procedure. This systematic variation is represented by the "Jet y ππ " label in the figures of Section 5.2.
The analysis is repeated at both a looser and a tighter PID selection than the nominal definition, and the variations are included as a systematic uncertainty. The effect on the radii is at the 1-2% level for the lower k T intervals, but becomes more significant at higher momentum, where there are relatively more kaons and protons and the dE/dx separation is not as large. In the highest k T intervals studied, variations are typically in the range of 5-30%. The PID systematic variation is labeled by "PID" in the figures of Section 5.2.
The nonzero effective size of the Coulomb correction R eff should only cause a bin-by-bin change of a few percent in the correlation function, even with a value up to several femtometers, since the Bohr radius of pion pairs is nearly 400 fm. However, since this parameter changes the width in q inv over which the Coulomb correction is applied, varying this parameter can affect the source radii measurably. The effective size is assumed to scale with the size of the source itself, so a scaling constant ξ is chosen such that R eff = ξ R inv . The nominal value of ξ is taken to be equal to 1 and the associated systematic uncertainty is evaluated by varying this between 1/2 and 2. The Coulomb size systematic variation is indicated by a label of "R eff " in the figures of Section 5.2.
A small difference between positive and negative charge pairs is observed, attributable to detector effects such as the orientation of the azimuthal overlap of the inner detector's component staves. The nominal results use all of the same-charge pairs, and a systematic variation accounting for this charge asymmetry is assigned which covers the results for both of the separate charge states. The variation from this effect is labeled by "+ + /−−" in the figures of Section 5.2.
Single-particle correction factors for track reconstruction efficiency cancel in the ratio A(q)/B(q). However, two-particle effects in the track reconstruction can affect the correlation function at small relative    momentum. Single-and multi-track reconstruction effects are both studied with the fully simulated H sample. The generator-level and reconstructed correlation functions are compared, and a deficit in the latter, due to the impact of the two-particle reconstruction efficiency, is observed at q inv below approximately 50 MeV. At larger q inv the two-particle reconstruction efficiency is found to not depend on q inv within statistical uncertainties. A minimum q cutoff is applied in the fits to minimize the impact of these detector effects. The sensitivity of the results to this cutoff is checked by taking q min inv = 30 ± 10 MeV in the one-dimensional fits, and symmetrizing the effect of the variation from |q| min = 25 to 50 MeV in the 3D fits. Because this variation has only a small effect on the radii, this procedure is taken to be sufficient to account for two-particle reconstruction effects. The effects of this variation have the label "Min q" in the figures of Section 5.2.

Magnitude of systematic effects
In this section the contributions of each source of systematic uncertainty are illustrated. Examples of the systematic uncertainties in the invariant parameters R inv and λ inv are shown as a function of k T and centrality in Figure 8 and 9.
Systematic uncertainties are also shown for the 3D radii R out (Figure 10), R side (Figure 11), R long (Figure 12), and R ol (Figure 13), as well as the ratio R out /R side (Figure 14) and the amplitude (Figure 15). These are all shown for typical choices of centrality, k T , and y ππ so that they represent standard, rather than exceptional, values of the uncertainties.
The uncertainties in the HBT radii (Figures 8, 10, 11 and 12) are dominated by the jet background    description. At larger k T the generator (P vs. Herwig) and system (pp vs. p+Pb) contributions constitute the larger portion of this, and at lower k T the variation of the mapping over y ππ is more significant.
The Bose-Einstein amplitudes λ inv (Figure 9) and λ osl (Figure 15) are also affected strongly by the jet fragmentation description, but at sufficiently large k T ( 0.4 GeV) pion identification contributes a comparable systematic uncertainty. This is expected because other particles misidentified as pions do not exhibit Bose-Einstein interference with real pions, and the most significant effect of their inclusion in the correlation function is to decrease the amplitude of the Bose-Einstein enhancement.    The systematic uncertainties in the ratio R out /R side (Figure 14) are estimated by evaluating the ratio after each variation reflecting a systematic uncertainty and taking the difference from the nominal value. Thus, the uncertainties that are correlated between R out and R side cancel properly in the ratio. The uncertainties in the ratio are not universally dominated by any single effect. At sufficiently large k T in central events, the effective Coulomb size becomes the largest contributor. This is understandable because the Coulomb correction is applied as a function of q inv , and as a result it is applied over a wider range in q out than in q side at larger k T .
Similarly, systematic effects in the cross-term R ol (Figure 13) are not dominated by any one source, and in fact the uncertainties in this quantity are predominantly statistical. At large k T the systematic uncertainties from PID appear large. However, these are mostly a result of statistical fluctuations that arise because the   fit in the tight PID selection is performed on a sample even smaller than the nominal dataset. Therefore, the reported systematic uncertainties are overly conservative for this quantity. For a similar reason, the systematic uncertainty for charge-asymmetric detector effects is not included in the error bars for R ol .
No systematic dependence of R ol is observed when measuring ++ and −− correlations independently, so they are excluded from the total systematic uncertainties in order not to include additional statistical fluctuations as systematic effects.

Results
This section shows examples of one-and three-dimensional fits to correlation functions, then presents results for extracted invariant and 3D source radii. The results are shown as a function of k T , which can illustrate the time-dependence of the source size. They are also shown as a function of y ππ , showing any variations in source size along the collision axis, and against several quantities related to multiplicity and centrality. These results show the freeze-out density and the evolution of the source with the size of the initial geometry.

Performance of fit procedure
An example of a one-dimensional fit to C(q inv ) using the functional form of Eq. 19 is included in Figure 6. Additional examples of one-dimensional fits for different k T intervals are shown in Figures 16-18   Slices of a three-dimensional fit of C(q) to the three-dimensional variant of Eq. 19 are shown in Figure 19. The apparently imperfect fit along the q out axis is characteristic of q side ≈ q long ≈ 0, and away from this slice the fit agrees better with the data (the test statistic per degree of freedom is 1.03 for the fit shown).

One-dimensional results
The results from fits of C(q inv ) to Eq. 19 for the invariant radius, R inv , are shown in Figure 20 in four selected centrality intervals. Only an intermediate rapidity interval −1 < y ππ < 0 is shown for these  Figure 19: Results of the 3D fit to the correlation function in the 0.4 < k T < 0.5 GeV, −1 < y ππ < 0 kinematic intervals and for the 10-20% centrality interval. The left, middle, and right panels show the distributions versus q out , q side , and q long , respectively, with limits on the other two components of q such that |q i | < 40 MeV. The dashed blue line indicates the description of the contribution from hard processes and the red line shows the full correlation function fit. The dotted red line indicates the extrapolation of the fit function beyond the interval over which the fit is performed. and similar results as a function of k T , as the qualitative behavior is consistent in forward and backward rapidities. The clear decrease in size with increasing k T that is observed in central events is not observed in peripheral events. This is consistent with the interpretation that central events undergo transverse expansion, since in hydrodynamic models higher-p T particles are more likely to freeze out earlier in the event. Another way of understanding this trend as evidence for transverse expansion is that there is a smaller homogeneity region for particles with higher p T [39]. At low k T , ultra-central (0-1%) events have an invariant radius significantly greater than peripheral (70-80%) events by a factor of about 2.6. This difference becomes less prominent at high k T . In central events R inv is larger on the lead-going side than on the proton-going side, while in peripheral events the rapidity dependence of the radius becomes constant.
Invariant radii are shown for several centralities in Figure 21 (left) as a function of the cube root of average dN ch /dη. For both k T intervals shown, the scaling of R inv with dN ch /dη 1/3 is close to linear but with a slightly increasing slope at higher multiplicities. The invariant radius, R inv , has a steeper trend versus multiplicity at lower k T . Figure 21 (right) shows R inv in several centrality and rapidity intervals as a function of the local particle density, dN ch /dy , which is evaluated by taking the average over the same interval used for the pair's rapidity. The extracted radius and the local particle density are seen to be tightly correlated, such that the radius can be predicted, within uncertainties, by the local density alone.
The Bose-Einstein amplitude of the invariant fits, λ inv , is shown in Figure 22 as a function of k T and y ππ . At low k T , λ inv has values near unity, and it decreases with rising k T . In the lower k T intervals a systematic difference is observed between centrality intervals, with λ inv having larger values in central events. In contrast, at larger k T the amplitudes are indistinguishable between different centralities. The amplitude exhibits no significant variation over rapidity.  Figure 20: The exponential invariant radii, R inv , obtained from one-dimensional fits to the q inv correlation functions shown as a function of pair transverse momentum, k T , (a) and rapidity, y ππ (b). Four non-adjacent centrality intervals are shown. The vertical size of each box represents the quadrature sum of the systematic uncertainties described in Section 5, and statistical uncertainties are shown with vertical lines. The horizontal positions of the points are the average k T or y ππ in each interval, and the horizontal lines indicate the standard deviation of k T or y ππ . The widths of the boxes differ among centrality intervals only for visual clarity.  Figure 21: Exponential fit results for R inv as a function of the cube root of average charged-particle multiplicity dN ch /dη 1/3 (a), where the average is taken over |η| < 1.5, and as a function of the local density, dN ch /dy , over several centrality and rapidity intervals (b). In the left plot the systematic uncertainties from pion identification and from the generator and collision system components of the background amplitude are treated as correlated and shown as error bands, and the systematic uncertainties from charge asymmetry, R eff , the rapidity variation of the jet fragmentation description, and two-particle reconstruction are treated as uncorrelated and indicated by the height of the boxes. The horizontal error bars indicate the systematic uncertainty from dN ch /dη or dN ch /dy .  Figure 22: The Bose-Einstein amplitude, λ inv , obtained from one-dimensional fits to the q inv correlation functions shown as a function of pair transverse momentum, k T , (a) and rapidity, y ππ (b). Four non-adjacent centrality intervals are shown. The vertical size of each box represents the quadrature sum of the systematic uncertainties described in Section 5, and statistical uncertainties are shown with vertical lines. The horizontal positions of the points are the average k T or y ππ in each interval, and the horizontal lines indicate the standard deviation of k T or y ππ . The widths of the boxes differ among centrality intervals only for visual clarity.

Three-dimensional results
The three-dimensional radii R out , R side , and R long are shown as a function of k T and y ππ in four selected centrality intervals in Figures 23-25. In central collisions, the 3D radii exhibit an even steeper decrease with increasing k T relative to that observed for the invariant radii in Figure 20. A similar, but weaker trend is present in peripheral events. Central collisions exhibit larger radii on the backward (Pb-going) side of the event, while peripheral events show no distinguishable variation of the radii with rapidity.
The 3D radii are also shown as a function of the cube root of both average event multiplicity and local density in Figures 26-28. These plots demonstrate the relationship between the size and the density of the source at freeze-out. All of the radii are seen to be very strongly correlated with the local density. The scaling of the radii is not far from being linear with the cube root of multiplicity. This behavior is qualitatively similar to the scaling of R inv with dN ch /dη in Figure 21.
The Bose-Einstein amplitude in the 3D fits, λ osl , is shown in Figure 29 as a function of k T and y ππ . Like the invariant amplitude, at low k T it is larger for central events than for peripheral ones. The three-dimensional amplitude does not decrease significantly with rising k T as the invariant amplitude does, except in the most peripheral events. The 3D amplitude also exhibits no significant variation over rapidity.
The ratio R out /R side ( Figure 30) is often studied because in models with radial flow, R out includes components of the source's lifetime but R side does not (see, for instance, the discussion in Ref. [31]). A value of R out /R side less than one is observed and it decreases with increasing k T . The ratio is observed to be the same in different centrality intervals within uncertainties. As explained in Ref.
[75], several improvements to naive hydrodynamic models-primarily pre-thermal acceleration, a stiffer equation of state, and shear viscosity-all result in more sudden emission. This implies that a value of R out /R side 1 does not necessarily rule out collective behavior.  Figure 23: Exponential fit results for the 3D source radius, R out , as a function of pair transverse momentum, k T , (a) and rapidity, y ππ (b). Four non-adjacent centrality intervals are shown. The vertical size of each box represents the quadrature sum of the systematic uncertainties described in Section 5, and statistical uncertainties are shown with vertical lines. The horizontal positions of the points are the average k T or y ππ in each interval, and the horizontal lines indicate the standard deviation of k T or y ππ . The widths of the boxes differ among centrality intervals only for visual clarity. The transverse area scale R out R side is shown in Figure 31 as a function of both event and local density. At lower k T , the transverse area scales linearly with multiplicity over all centralities and rapidities. This result is consistent with a picture in which the longitudinal dynamics can be separated from the transverse particle production, and low-k T particles freeze out at a constant transverse area density.
The determinant of the 3D radius matrix, det (R) (Eq. 6), is shown in Figure 32 as a function of both the average and local density. While the transverse area scales linearly with multiplicity at low k T , the volume scale grows linearly at higher k T , implying a constant freeze-out volume density for particles with higher momentum. Figure 33 compares the volume scaling with N part for the standard Glauber model as well as for two choices of the Glauber-Gribov color fluctuation (GGCF) model [61]. The parameter ω σ controls the size of the fluctuations in the nucleon-nucleon cross section within the Glauber-Gribov model. With the Glauber model, the scaling of the volume element with N part has a significant upwards curvature. Including Glauber-Gribov fluctuations in the N part calculation results in a more modest curvature in the scaling of det (R). This result suggests that the fluctuations in the nucleon-nucleon cross section are a crucial component of the initial geometry description in p+Pb systems. The values and systematic uncertainties of N part in each model are listed in Table 1.
The cross-term, R ol (Eq. 6), which couples to the lifetime of the source [70], is shown in Figures 34 and  35. A significant departure from zero is observed in this parameter in central events, but only for rapidities y ππ −1. For the 0-1% centrality interval, in 0.2 < k T < 0.4 and −1 < y ππ < 1, R ol is measured to be nonzero with a significance of 7.1/7.3/5.1 σ (statistical/systematic/combined). The next most central interval, 1-5%, has a nonzero R ol with a significance of 5.2/5.8/3.9 σ (statistical/systematic/combined). This suggests that the particle production at middle and forward rapidities is sensitive to the local zasymmetry of the system.
The argument from Section 4.1 for why the order of the particles in a pair can be chosen so that q out is greater than zero relies on the assumption that both particles in the pair are the same species, or at least that   Figure 26: Exponential fit results for R out as a function of (a) the cube root of average charged-particle multiplicity, dN ch /dη 1/3 , where the average is taken over |η| < 1.5 and (b) the local density, dN ch /dy , in intervals of y ππ . In the left plot the systematic uncertainties from pion identification and from the generator and collision system components of the background amplitude are treated as correlated and shown as error bands, while the systematic uncertainties from charge asymmetry, R eff , the rapidity variation of the jet fragmentation description, and twoparticle reconstruction are treated as uncorrelated and indicated by the height of the boxes. The horizontal error bars indicate the systematic uncertainty from dN ch /dη or dN ch /dy . they are characterised by the same momentum distributions. In principle, final-state interactions between different particle species could break this symmetry of the correlation function and lead to a nonzero R ol term. However, the systematic uncertainties shown in Figure 13 demonstrate that R ol is not sensitive to particle identification, particularly at low k T . At larger k T the systematic effect from PID looks larger, but the variations are likely driven by statistical fluctuations.  Figure 27: Exponential fit results for R side as a function of (a) the cube root of average charged-particle multiplicity, dN ch /dη 1/3 , where the average is taken over |η| < 1.5 and (b) the local density, dN ch /dy , in intervals of y ππ . In the left plot the systematic uncertainties from pion identification and from the generator and collision system components of the background amplitude are treated as correlated and shown as error bands, while the systematic uncertainties from charge asymmetry, R eff , the rapidity variation of the jet fragmentation description, and twoparticle reconstruction are treated as uncorrelated and indicated by the height of the boxes. The horizontal error bars indicate the systematic uncertainty from dN ch /dη or dN ch /dy .  Figure 28: Exponential fit results for R long as a function of (a) the cube root of average charged-particle multiplicity, dN ch /dη 1/3 , where the average is taken over |η| < 1.5 and (b) the local density, dN ch /dy , in intervals of y ππ . In the left plot the systematic uncertainties from pion identification and from the generator and collision system components of the background amplitude are treated as correlated and shown as error bands, while the systematic uncertainties from charge asymmetry, R eff , the rapidity variation of the jet fragmentation description, and twoparticle reconstruction are treated as uncorrelated and indicated by the height of the boxes. The horizontal error bars indicate the systematic uncertainty from dN ch /dη or dN ch /dy .  Figure 29: The Bose-Einstein amplitude, λ osl , as a function of pair transverse momentum, k T , (la) and rapidity, y ππ (b). Four non-adjacent centrality intervals are shown. The vertical size of each box represents the quadrature sum of the systematic uncertainties described in Section 5, and statistical uncertainties are shown with vertical lines. The horizontal positions of the points are the average k T or y ππ in each interval, and the horizontal lines indicate the standard deviation of k T or y ππ . The widths of the boxes differ among centrality intervals only for visual clarity.  Figure 30: The ratio of exponential radii R out /R side as a function of pair transverse momentum, k T , (a) and rapidity, y ππ (b). Four non-adjacent centrality intervals are shown. The vertical size of each box represents the quadrature sum of the systematic uncertainties described in Section 5, and statistical uncertainties are shown with vertical lines. The horizontal positions of the points are the average k T or y ππ in each interval, and the horizontal lines indicate the standard deviation of k T or y ππ . The widths of the boxes differ among centrality intervals only for visual clarity.  Figure 31: The transverse area scale, R out R side , plotted against the average multiplicity, dN ch /dη , (a) and the local density, dN ch /dy , as a function of rapidity (b). The systematic uncertainties from pion identification and the generator and collision system components of the jet background description are treated as correlated and shown as error bands. The systematic uncertainties from charge asymmetry, R eff , rapidity variation of the jet fragmentation, and two-particle track reconstruction effects are treated as uncorrelated and indicated by the height of the boxes. The horizontal error bars indicate the systematic uncertainty from dN ch /dη or dN ch /dy . The slope and intercept of the best fit to the right-hand plot are shown with combined statistical and systematic uncertainties. * y -2 < < -1 π π * y -1.5 < < -0.5 π π * y -1 < < 0 π π * y -0.5 < < 0.5 π π * y 0 < < 1 π π * y 0.5 <  Figure 32: The volume scale, det (R), plotted against the average multiplicity, dN ch /dη , (a) and the local density, dN ch /dy , as a function of rapidity (b). In the left plot, the systematic uncertainties from pion identification and the generator and collision system components of the jet background description are treated as correlated and shown as error bands, while those from charge asymmetry, R eff , rapidity variation of the jet fragmentation, and two-particle track reconstruction effects are treated as uncorrelated and indicated by the height of the boxes. The widths of the boxes indicate the systematic uncertainty in dN ch /dη or dN ch /dy .  Figure 33: The scaling of the volume element, det (R), with N part calculated with three initial geometry models: standard Glauber as well as Glauber-Gribov (GGCF) for two choices of the color fluctuation parameter, ω σ . Each of the panels shows a different k T interval. The systematic uncertainties from the pion identification and the generator and collision system components of the background description are treated as correlated and shown as error bands. The systematic uncertainties from charge asymmetry, R eff , rapidity variation of the jet background, and two-particle reconstruction are treated as uncorrelated and indicated by the height of the boxes. The horizontal error bars indicate the systematic uncertainties in N part .   Figure 35: The cross-term, R ol , as a function of average event multiplicity, dN ch /dη , (a) and the local density, dN ch /dy (b). The vertical size of each box represents the quadrature sum of the systematic uncertainties described in Section 5, and statistical uncertainties are shown with vertical lines. The widths of the boxes indicate the systematic uncertainties in the corresponding quantities.

Summary and conclusions
This paper presents ATLAS measurements of two-identified-pion HBT correlations in p+Pb collisions at the LHC at √ s NN = 5.02 TeV using a total integrated luminosity of 28 nb −1 . Two-particle correlation functions were measured in one dimension as a function of q inv and in three dimensions using the out-sidelong decomposition as a function of the pair's average transverse momentum and the pair's rapidity. The measurements were performed for several intervals of p+Pb centrality characterized by ΣE Pb T , the total transverse energy measured in the Pb-going forward calorimeter. A data-driven technique was developed for constraining the contribution of jet fragmentation to the correlation function. The correlation functions were fit to a Bowler-Sinyukov parameterization of the Bose-Einstein and Coulomb form with a contribution that describes the jet background. The HBT correlation is described by an exponential form that provides a good description of the data. For the out-side-long fits, the parameterization includes a non-diagonal coupling between q out and q long in the correlation function.
The radii extracted from the one-dimensional and three-dimensional fits show a significant variation with transverse momentum k T that is strongest for the most central events and weakest or not present in the most peripheral centrality interval. For the three-dimensional fits, the k T dependence is found to be the largest for the out and long directions. A small but significant dependence of the three-dimensional source radii on the pair's rapidity is observed in the more central collisions while in the most peripheral collisions the radii do not depend on rapidity.
The one-dimensional and three-dimensional source radii increase monotonically between peripheral and central collisions with a slope that decreases with rising k T . The dependence of the radii on centrality was studied as a function of both the cube root of the rapidity-averaged charged-particle multiplicity, dN ch /dη 1/3 , and the local charged-particle density, dN ch /dy . When evaluated in intervals of both centrality and rapidity, the radii as a function of dN ch /dy fall on a single curve. At low k T the rapidityaveraged radii are observed to increase approximately linearly with dN ch /dη 1/3 . A nonzero out-long cross-term R ol is observed in central (0-5%) collisions for rapidities greater than −1, with a combined significance of 5.1σ in the most central (0-1%) events.
The transverse area, R side R long , is observed to vary linearly with dN ch /dy at low k T . The volume scale represented by det (R) increases faster than linearly with dN ch /dη 1/3 or dN ch /dy at low k T , but increases approximately linearly with dN ch /dη 1/3 at higher momentum (k T > 0.5 GeV). When plotted versus the mean number of nucleon participants N part obtained from three different geometric models, the volume shows a steady increase with N part for the GGCF-derived N part values, but a sudden increase with N part for N part 12 when N part is obtained from the Glauber model. While the freeze-out volume scale det (R) should only be strictly linear with the initial size, represented by N part , if the expansion is independent of centrality, an extreme deviation from a naive linear scaling is not expected. This observation supports the conclusion drawn from previous studies that the Glauber-Gribov color fluctuation model provides a better description of p+Pb collision geometry.
The R out to R side ratio is found to be less than unity for all centrality and kinematic selections studied in this analysis, and it is observed to decrease approximately linearly with increasing k T . This result, combined with the k T dependence of the radii, suggests a collective, explosive expansion of the source. The nonzero out-long cross-term indicates that the freeze-out behavior of the source is sensitive to the local z-asymmetry of the particle production away from the Pb-going region.
The results presented in this paper provide detailed measurements of the space-time extent of the particle source in p+Pb systems. In particular, the rapidity sensitivity of the results demonstrate the asymmetry of the p+Pb system and show that k T and local charged-particle density are sufficient to predict the radii. These conclusions present a significant opportunity for theoretical models.
[4] ATLAS Collaboration, Measurement of long-range pseudorapidity correlations and azimuthal harmonics in √ s NN = 5.02 TeV proton-lead collisions with the ATLAS detector,