Two-pion femtoscopy in p-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV

We report the results of the femtoscopic analysis of pairs of identical pions measured in p-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV. Femtoscopic radii are determined as a function of event multiplicity and pair momentum in three spatial dimensions. As in the pp collision system, the analysis is complicated by the presence of sizable background correlation structures in addition to the femtoscopic signal. The radii increase with event multiplicity and decrease with pair transverse momentum. When taken at comparable multiplicity, the radii measured in p-Pb collisions, at high multiplicity and low pair transverse momentum, are 10-20% higher than those observed in pp collisions but below those observed in A-A collisions. The results are compared to hydrodynamic predictions at large event multiplicity as well as discussed in the context of calculations based on gluon saturation.


Introduction
The Large Hadron Collider (LHC) [1] delivered Pb-Pb collisions at √ s NN = 2.76 TeV in 2010 and in 2011. Several signatures of quark-gluon plasma were observed, including a strong suppression of highp T particle production ("jet quenching") [2][3][4], as well as collective behavior at low p T [5,6] which is well described by hydrodynamic models with a low shear viscosity to entropy ratio. A comparison to reference results from pp collisions at √ s = 0.9, 2.76, and 7 TeV shows that these effects cannot be described by an incoherent superposition of nucleon-nucleon collisions. As such, they can be interpreted as final-state phenomena, characteristic of the new state of matter [7][8][9][10] created in heavy-ion collisions. In order to verify the creation of such a state, p-Pb collisions at √ s NN = 5.02 TeV were provided by the LHC. In particular, cold nuclear matter effects, such as gluon saturation, which are expected to influence a number of observables, are being investigated [11].
One of the observables characterizing the bulk collective system is the size of the particle-emitting region at freeze-out which can be extracted with femtoscopic techniques [12,13]. In particular, two-particle correlations of identical pions (referred to as Bose-Einstein, or Hanbury-Brown Twiss "HBT", correlations) provide a detailed picture of the system size and its dependence on the pair transverse momentum and the event multiplicity. Femtoscopy measures the apparent width of the distribution of relative separation of emission points, which is conventionally called the "radius parameter." The radius can be determined independently for three directions: long along the beam axis, out along the pair transverse momentum, and side, perpendicular to the other two. Such measurements were performed at the LHC for central Pb-Pb collisions [14] as well as for pp collisions at √ s = 0.9 and 7 TeV [15][16][17][18], and compared to results from heavy-ion collisions at lower collision energies. Two clear trends were found: (i) In A-A collisions all radii scale approximately linearly with the cube root of the final state charged particle multiplicity density at midrapidity dN ch /dη 1/3 for all three radii separately, consistently with previous findings [13]. For pp collisions, the radii scale linearly as well, however the slope and intercept of the scaling line are clearly different than for A-A. (ii) A significant, universal decrease of the radii with pair momentum has been observed in A-A collisions, while the analogous trend in pp depends on the considered direction (out, side, or long) and event multiplicity. A transverse momentum dependence of the radii similar to A-A was observed for the asymmetric d-Au collision system at RHIC [19,20]. The one-dimensional radii extracted from the ALICE 3-pion cumulant analysis were also investigated in pp, p-Pb, and Pb-Pb collisions at the LHC [21]. For the p-Pb system, at a given multiplicity, the radii were found to be 5-15% larger than those in pp, while the radii in Pb-Pb were 35-55% larger than those in p-Pb.
The A-A pion femtoscopy results are interpreted within the hydrodynamic framework as a signature of collective radial flow. Models including this effect are able to reproduce the ALICE data for central collisions [22,23]. Attempts to describe the pp data in the same framework have not been successful so far and it is speculated that additional effects related to the uncertainty principle may play a role in such small systems [24]. In p-A collisions, hydrodynamic models which assume the creation of a hot and dense system expanding hydrodynamically predict system sizes larger than those observed in pp, and comparable to those observed in lower-energy A-A collisions at the same multiplicity [24,25]. However such models have an inherent uncertainty of the initial state shape and size, which can also differ between pp and peripheral A-A collisions.
Alternatively, a model based on gluon saturation suggests that the initial system size in p-A collisions should be similar to that observed in pp collisions, at least in the transverse direction [26,27]. At that stage both systems are treated in the same manner in the Color Glass Condensate model (CGC), so that their subsequent evolution should lead to comparable radius parameter at freeze-out. Ref. [28] suggests a (small) Yang-Mills evolution in addition. The observation of a larger size in the p-A system with respect to pp would mean that a comparable initial state evolves differently in the two cases, which is not easily explained within the CGC approach alone. The d-Au data measured at RHIC suggest that hydrodynamic evolution may be present in such a system, while the ALICE 3-pion analysis at the LHC [21] leaves room for different interpretations. The pion femtoscopic radii as a function of pair transverse momentum from p-Pb collisions at √ s NN = 5.02 TeV, which are reported in this paper, provide additional constraints on the validity of both approaches.
The paper is organized as follows: In Sec. 2 the data-taking conditions, together with event and track selections are described. The femtoscopic correlation function analysis, as well as the extraction of the radii and associated systematic uncertainties, and the discussion of the fitting procedure are explained in Sec. 3. In Sec. 4 the results for the radii are shown and compared to model predictions. Section 5 concludes the paper.

Data taking and track reconstruction
The LHC delivered p-Pb collisions at the beginning of 2013 at √ s NN = 5.02 TeV (4 TeV and 1.58 TeV per nucleon for the p and Pb beams, respectively). The nucleon-nucleon center-of-mass system is shifted with respect to the ALICE laboratory system by 0.465 unit of rapidity in the direction of the proton beam.
The ALICE detector and its performance are described in Refs. [29,30]. The main triggering detector is the V0, consisting of two arrays of 32 scintillator counters, which are installed on each side of the interaction point and cover 2.8 < η lab < 5.1 (V0A, located on the Pb-remnant side), and −3.7 < η lab < −1.7 (V0C). The minimum-bias trigger requires a signal in both V0 detectors within a time window that is consistent with the collision occurring at the center of the ALICE detector. Additionally, specific selection criteria to remove pile-up collisions are applied [30]. Approximately 80 million minimum-bias events were analyzed.
The analysis was performed in multiplicity classes, which were determined based on the signal from the V0A detector, located along the Pb-going beam. This ensures that the multiplicity determination procedure uses particles at rapidities significantly different from the ones used for the pion correlation analysis, avoiding potential auto-correlation effects. Events are grouped in four multiplicity classes: 0-20%, 20-40%, 40-60%, and 60-90%, defined as fractions of the analyzed event sample sorted by decreasing V0A signal, which is proportional to the multiplicity within the acceptance of this detector. Table 1 shows the multiplicity class definitions and the corresponding mean charged-particle multiplicity densities dN ch /dη averaged over |η lab | < 0.5 as obtained using the method presented in Ref. [31]. The dN ch /dη values are not corrected for trigger and vertex-reconstruction inefficiency, which is about 4% for non-single diffractive events [31].

Event class
23.2 ± 0.5 0-20% 35.5 ± 0.8 Table 1: Definition of the V0A multiplicity classes as fractions of the analyzed event sample and their corresponding dN ch /dη(|η lab | < 0.5, p T > 0) . The given uncertainties are systematic only since the statistical uncertainties are negligible.
Charged track reconstruction is performed using the Time Projection Chamber (TPC) and the Inner Tracking System (ITS). The TPC is a large volume cylindrical gaseous tracking chamber, providing information of particle trajectories and their specific energy loss. The read-out chambers mounted on the endcaps are arranged in 18 sectors on each side (covering the full azimuthal angle) measuring up to 159 samples per track. The TPC covers an acceptance of |η lab | < 0.8 for tracks which reach the outer radius of the TPC and |η lab | < 1.5 for shorter tracks. The ITS is composed of position-sensitive silicon detectors. It consists of six cylindrical layers: two layers of Silicon Pixel Detector (SPD) closest to the beam pipe covering |η lab | < 2.0 and |η lab | < 1.4 for inner and outer layers respectively, two layers of Silicon Drift Detector in the middle covering |η lab | < 0.9, and two layers of Silicon Strip Detector on the outside covering |η lab | < 1.0. The information from the ITS is used for tracking and primary particle selection. The momentum of each track is determined from its curvature in the uniform magnetic field of 0.5 T oriented along the beam axis, provided by the ALICE solenoidal magnet.
The primary-vertex position is determined with tracks reconstructed in the ITS and TPC as described in Ref. [32]. Events are selected if the vertex position along the beam direction is within ±10 cm of the center of the detector. This ensures a uniform acceptance in η lab .
Each track is required to exploit signals in both TPC and ITS. The track segments from both detectors have to match. Additionally, each track is required to have at least one hit in the SPD. A TPC track segment is reconstructed from space points (clusters). Each track is required to be composed of at least 50 out of the 159 such clusters. The parameters of the track are determined by a Kalman fit to the set of TPC+ITS clusters. The quality of the fit χ 2 was required to be better than 4 per cluster in the TPC and better than 36 in ITS. Tracks that show a kink topology in the TPC are rejected. To ensure that dominantly primary-particle tracks are selected, the distance of closest approach to the primary vertex is required to be closer than 2.0 cm in the longitudinal direction and (0.0105 + 0.0350 · p −1.1 T ) cm, with p T in GeV/c, in the transverse direction. The kinematic range of particles selected for this analysis is 0.12 < p T < 4.0 GeV/c and |η lab | < 0.8.
The Time-Of-Flight (TOF) detector is used together with the TPC for pion identification. The TOF is a cylindrical detector of modular structure, consisting of 18 azimuthal sectors divided into 5 modules along the beam axis at a radius r ≃ 380 cm. The active elements are Multigap Resistive Chambers (MRPC). For both TPC and TOF, the signal (specific energy loss dE/dx for the TPC and the time-of-flight for TOF) for each reconstructed particle is compared with the one expected for a pion. The difference is confronted with the detector resolution. The allowed deviations vary between 2 to 5σ for the TPC and 2 to 3σ for TOF depending on the momentum of the particle. The selection criteria are optimized to obtain a high-purity sample while maximizing efficiency, especially in the regions where the expected signal for other particles (electrons, kaons, and protons for the TPC, kaons for TOF) approaches the pion value. The purity of the pion sample is above 98%.
The accepted particles from each event are combined to pairs. The two-particle detector acceptance effects, track splitting and track merging, are present. Track splitting occurs when a single trajectory is mistakenly reconstructed as two tracks. ALICE tracking algorithm has been specifically designed to suppress such cases. In a rare event when splitting happens, two tracks are reconstructed mostly from the same clusters in the ALICE TPC. Therefore pairs which share more than 5% of clusters are removed from the sample. Together with the anti-merging cut described below this eliminates the influence of the split pairs. Track merging can be understood as two-particle correlated efficiency and separation power. In the ALICE TPC, two tracks cannot be distinguished if their trajectories stay close to each other through a significant part of the TPC volume. Although this happens rarely such pairs by definition have low relative momentum and therefore their absence distorts the correlation function in the signal region. Track splitting and track merging are taken into account and corrected with the procedure described in Ref. [16]. The effect of the two-particle detector acceptance on the final results is similar to what was observed in pp and is limited to low pair relative momentum, where it slightly affects the shape of the correlation function. However, in p-Pb collisions the femtoscopic effect is an order of magnitude wider than any region affected by this inefficiency and, as a consequence, the extracted radii are not affected by the two-track acceptance.   Fig. 1: Projections of the three-dimensional π + π + correlation functions for three selected multiplicity and k T ranges along the out (top), side (middle), and long (bottom) direction. The other components are integrated over the four bins closest to zero in their respective q directions.

Construction of the correlation function
The correlation function C(p 1 , p 2 ) of two particles with momenta p 1 and p 2 is defined as The signal distribution A is constructed from pairs of particles from the same event. The background distribution B should be constructed from uncorrelated particles measured with the same single-particle acceptance. It is built using the event mixing method with the two particles coming from two different events for which the vertex positions in beam direction agree within 2 cm and the multiplicities differ by no more than 1/4 of the width of the given event class. The denominator is normalized to the number of entries in the numerator, so that the absence of correlation gives a correlation function at unity.
The femtoscopic correlation is measured as a function of the momentum difference of the pair q = p 2 −p 1 as The momentum difference q is evaluated in the Longitudinally Co-Moving System (LCMS) frame in which the total longitudinal pair momentum vanishes: p 1,L + p 2,L = 0, similarly to previous measurements in small systems [16]. In Fig. 1 correlation functions are shown, projected over 128 MeV/c-wide slices along the q out , q side , and q long axes. An enhancement at low relative momentum is seen in all projections. The width of this correlation peak grows with decreasing multiplicity or with increasing k T . The femtoscopic effect is expected to disappear at large q = |q|, with the correlation function approaching unity. We observe, especially for large k T and small multiplicities, that the correlation function is not flat in this region and has different values in different projections 1 . The cause may be non-femtoscopic correlations, which are presumably also affecting the shape of the correlation function in the femtoscopic (low q) region. This issue is a major source of systematic uncertainty on the extracted radii and is discussed The pair distributions and the correlation function can be represented in spherical harmonics (SH) [33,34] alternatively to the traditionally-used Cartesian coordinates. All odd-l and odd-m components of such a representation vanish for symmetry reasons. The important features of the correlation function are fully captured by the following ones: l = 0, m = 0 is sensitive to the overall size of the pion source, l = 2, m = 0 is sensitive to the difference between the longitudinal and transverse sizes, and l = 2, m = 2 reflects the difference between the sidewards and outwards transverse radii. Therefore, three independent sizes of the source can also be extracted from these three SH components.
In Fig. 2 we show the first three non-vanishing components of the spherical harmonics representation corresponding to the correlation functions shown in Fig. 1. In the (0, 0) component the enhancement at low-q is clearly visible, decreasing (increasing) in width with multiplicity (k T ). The other two components, (2, 0) and (2, 2), also show structures in this region, indicating that the source shape is not spherically symmetric in the LCMS frame.

Non-femtoscopic structures
As mentioned in the discussion of Figs. 1 and 2, a significant non-femtoscopic correlation is observed in the range in q that is much larger than the characteristic width of the femtoscopic effect. As an example, in Fig. 3 we show the correlation in the SH representation up to 2.0 GeV/c in q. For the lowest Fig. 4: First three non-vanishing components of the SH representation of the π + π + correlation functions for a selected multiplicity and k T range, compared to a calculation from EPOS 3.076 [35,36] (generator level only) and PYTHIA 6.4 Perugia-0 tune [37,38] for pp at √ s = 7 TeV (full simulation with detector response).
multiplicity, and to a smaller degree at higher multiplicities, a significant slope in low q region is seen in the (0, 0) component and a deviation from zero in the (2, 0) component up to approximately 1 GeV/c. Similar correlations have been observed by ALICE in pp collisions [16]. They were interpreted, based on Monte-Carlo model simulations, to be a manifestation of mini-jets, the collimated fragmentation of partons scattered with modest momentum transfer. The lowest multiplicities observed in p-Pb collisions are comparable to those in pp collisions at √ s = 7 TeV. Therefore a similar interpretation of the non-femtoscopic correlations in this analysis is natural. Similar structures have been observed in d-Au collisions by STAR [19]. This picture is corroborated by the analysis using the 3-pion cumulants, where expectedly the mini-jet contribution is suppressed [21].
Two important features of the non-femtoscopic correlation affect the interpretation of our results. First, it is a broad structure, extending up to 1 GeV/c and we have to assume that it also extends to 0 GeV/c in q. Therefore it affects the extracted femtoscopic radii and has to be taken into account in the fitting procedure. It can be quantified in the large q region and then extrapolated, with some assumptions, to the low q region, under the femtoscopic peak. The procedure leads to a systematic uncertainty. Secondly, it becomes visibly larger as multiplicity decreases and also as k T increases, which is consistent with the mini-jet-like correlation.
The background was studied in Monte Carlo models, such as, for p-Pb collisions, AMPT [39], HIJING [40], DPMJET [41], EPOS 3.076 [35,36], and PYTHIA 6.4 (Perugia-0 tune) [37,38] for pp collisions at similar multiplicities. In all cases the Monte Carlo correlation functions exhibit significant structures similar to the long-range effects observed in data, which is another argument for their non-femtoscopic origin. However, quantitative differences in the magnitude and shape of these structures when compared to those observed in data are seen for AMPT, HIJING and DPMJET. These models are therefore unsuitable for a precise characterization of the background, which is needed for the fitting procedure. The only models that qualitatively describe the features of the background (enhancement at low q, growing with k T , and falling with multiplicity) are EPOS 3.076 and PYTHIA 6.4 (Perugia-0 tune), which was also used in the pp analysis [16]. We note that PYTHIA simulation included full detector response modelling, while the EPOS 3.076 one did not. The comparison with data is shown in Fig. 4. The behavior of the correlation is well reproduced above 0.5 GeV/c in q, where non-femtoscopic correlations are expected to dominate. At low q, below 0.3 GeV/c, the data and models diverge, which is expected, as the femtoscopic correlations are not included in the model calculation. Above 0.3 GeV/c, EPOS reproduces the (0, 0) component well, PYTHIA slightly overestimates the data. For the (2, 0) and (2, 2) components, which describe the three-dimensional shape of the non-femtoscopic correlations, PYTHIA is closer to the experimental data. Overall, for like-sign pairs, both models are reasonable approximations of the non-femtoscopic background. We use these models to fix the background parameters in the fitting procedure.
Similarly to the pp analysis [16], the unlike-sign pairs have also been studied. We found that correlations in the (0, 0) component of PYTHIA are slightly larger than in data in the femtoscopic region for all k T ranges, and similar to data in the (2, 0) and (2, 2) components. EPOS was found to reasonably describe the unlike-sign pairs for low k T ranges and has smaller correlation than data in the (0, 0) component in higher k T ranges.

Fitting the correlation functions
The space-time characteristics of the source are reflected in the correlation function where r is the pair space-time separation four-vector. S is the source emission function, interpreted as a probability density describing the emission of a pair of particles with a given relative momentum and space-time separation. Ψ is the two-particle interaction kernel.
The Gaussian source provides a commonly used approximation of the source size and was used to compare to other experimental results, especially the ones coming from A-A collisions, where the source shape is more Gaussian than in small systems. While pursuing the standard procedure with the Gaussian assumption, we also carefully look for any deviations between the fit function and data which might suggest a significantly non-Gaussian shape of the source, which would be an important similarity to the pp case.
In the analysis of pp collisions by ALICE [16], a Gaussian is used together with other source shapes, exponential and Lorentzian [16]. A Lorentzian parametrization in the out and long directions and a Gaussian parametrization in the side direction were found to fit the data best according to χ 2 /ndf. Therefore, we use this source parametrization also in the analysis of p-Pb collisions The corresponding source sizes in out and long are R E out and R E long , while for the side direction the size parameter R G side is the same as in the Gaussian case. For identical pions, which are bosons, Ψ must be symmetrized. Since charged pions also interact via the Coulomb and strong Final State Interactions (FSI), |Ψ| 2 corresponds to the Bethe-Salpeter amplitude [50]. For like-sign pion pairs the contribution of the strong interaction is small for the expected source sizes [50], and is neglected here. The used Ψ therefore is a symmetrized Coulomb wave function. It is approximated by separating the Coulomb part and integrating it separately, following the procedure known as Bowler-Sinyukov fitting [51,52], which was used previously for larger sizes observed in Pb-Pb [14] as well as smaller sizes observed in pp collisions [16]. In this approximation the integration of Eq. (3) with S given by Eq. (4) results in the following functional form for the correlation function which is used to fit the data The function K C (q) is the Coulomb part of the two-pion wave function integrated over the spherical Gaussian source with a fixed radius. The value of this radius is chosen to be 2 fm. Its uncertainty has systematic effects on the final results (see Sec. 3.4). This form of the correlation function from Eq. (6) is denoted in the following as GGG. Similarly for the source shape given by Eq. (5) the correlation function is It has an exponential shape in out and long and a Gaussian shape in the side direction. Therefore it is referred to as EGE form of the correlation function. Parameter λ in Eqs. (6) and (7) represents the correlation strength.
Additionally, a component Ω describing non-femtoscopic correlations needs to be introduced. There is no a priori functional form which can be used for this component. Several of its features can be deduced from the correlations shown in Figs. 1, 2, and 3: it has to allow for different shapes in the out, side, and long directions; in the (0, 0) component it has to extrapolate smoothly to low q and have a vanishing slope at q = 0. Since this structure is not known, maximum information about its shape and magnitude should be gained from an observation of the raw correlation functions and the corresponding effects in Monte Carlo, in as many representations as possible. It is therefore crucial to simultaneously use the Cartesian and SH representations as they provide complementary ways to study the correlation shape.
The fit is performed with the log-likelihood method in three dimensions for the Cartesian representation. The Gaussian fit reproduces the overall width of the femtoscopic correlation in all cases. The background component describes the behavior of the correlation at large q, but can also have non-zero correlation at 0 in q.
A corresponding fit is also performed for the SH representation of the correlation, which is shown in Fig. 5. The formula from Eq. (6) or Eq. (7) (for the GGG fit or the EGE fit, respectively) is numerically integrated on a ϕ − θ sphere for each q bin, with proper Y m l weights, to produce the three components of the SH decomposition. Statistical uncertainties on each component as well as the covariance matrix between them are taken into account in this simultaneous fit to the three histograms. The results are shown in Fig. 5. The fit describes the general direction-averaged width of the correlation function, shown in the upper panel. The background component Ω describes the behavior at large q but also contributes to the correlation at low q. The shape in three-dimensional space, captured by the (2, 0) and (2, 2) components, is also a combination of the femtoscopic and non-femtoscopic correlations.
Overall the GGG fit describes the width of the correlation but the data at low q are not perfectly reproduced, which can be attributed to the limitations of the Bowler-Sinyukov formula as well as to the non-Gaussian, long-range tails which are possibly present in the source. Some deviations from the pure Gaussian shape can also be seen for the long direction for the higher multiplicities. The EGE fit (Eq. (7)) better reproduces the correlation peak in the (0,0) component, as shown in Fig. 6. The (2,0) and (2,2) components show similar quality of the fit. The χ 2 values for both fits are comparable.

Systematic uncertainties on the radii
The analysis was performed separately for positively and negatively charged pions. For the practically zero-net-baryon-density system produced at the LHC they are expected to give consistent results. Both datasets are statistically consistent at the correlation function level.
The main contributions to the systematic uncertainty are given in Table 2 for GGG radii and in Table 3 for EGE radii.
We used two alternative representations (Cartesian and spherical harmonics) of the correlation function.   The same functional form for both of them was used for the fitting procedure. However, the implementation of the fitting procedure is quite different: log-likelihood for Cartesian vs. regular χ 2 fit for SH, 3D Cartesian histogram vs. three 1D histograms, cubic or spherical fitting range in the (q out , q side , q long ) space. Therefore the fits to the two representations may react in a systematically different way to the variation of the fitting procedure (fit ranges, Bowler-Sinyukov approximation, etc.).
The fitting procedure requires the knowledge of the non-femtoscopic background shape and magnitude. Two models were used to estimate it, EPOS 3.076 [36] and PYTHIA 6.4 (Perugia-0 tune) [37,38], as described in Sec. 3.2.
In addition the correlation function shape is not ideally described by a Gaussian form. The EGE form is better (lower χ 2 values for the fit), but still not exactly accurate. As a result the fit values depend on the fitting range used in the procedure of radius extraction. We have performed fits with an upper limit of the fit range varied between 0.3 GeV/c and 1.1 GeV/c.
The three effects mentioned above are the main sources of systematic uncertainty on the radii. Their influence, averaged over the event multiplicity and pair k T , is given in Table 2 and Table 3. The background parametrization and the CF representation effects lead to systematic uncertainties less than 10% at low k T and up to 35% for large k T and low multiplicities. In particular, the radius could not be reliably extracted for the two highest k T ranges in the lowest multiplicity range, therefore these two sets of radii are not shown. Moreover, radii obtained with the background parametrization from PYTHIA are always larger than the ones obtained with the EPOS parametrization. These uncertainties are correlated between k T ranges. Similarly, the radii from the narrow fit range are always on average 10% higher than the ones from the wide fit range. This also gives a correlated contribution to the systematic uncertainty. The final radii are calculated as an average of four sets of radii -the two representations with both EPOS and PYTHIA background parametrization. The systematic uncertainties are symmetric and equal to the largest difference between the radius and one of the four sets of radii.
The effect of the momentum resolution on the correlation function was studied using a Monte Carlo simulation. For tracks with a low p T , below 1 GeV/c, the momentum resolution in the TPC is better than 1%. Smearing of the single particle momenta reduces the height and increases the width of the correlation function. It was estimated that this effect changes the reconstructed radius by 2% for a system size of 2 fm and 3% for a size of 3 fm. Therefore, the 3% correlated contribution from momentum resolution is always added to the final systematic uncertainty estimation.
Smaller sources of systematic uncertainties include those originating from the difference between positively and negatively charged pion pairs (around 3%), track selection variation (less than 1%) and from the Coulomb factor (less than 1%). All the uncertainties are added in quadrature.

Three-dimensional radii
We have extracted R out , R side , and R long in intervals of multiplicity and k T , which results in 26 radii in each direction. The fit procedure did not allow to reliably extract values of the radii for the two highest k T ranges in the 60-90% multiplicity class. For the GGG fit, they are shown in Fig. 7. The radii are in the range of 0.6 to 2.4 fm in all directions and universally decrease with k T . The magnitude of this decrease is similar for all multiplicities in the out and long directions, and is visibly increasing with multiplicity in the side direction. The radii rise with event multiplicity. The plot also shows data from pp collisions at √ s = 7 TeV [16] at the highest multiplicities measured by ALICE, which is slightly higher than the multiplicity measured for the 20-40% V0A signal range in the p-Pb analysis. At small k T the pp radii are lower by 10% (for side) to 20% (for out) than the p-Pb radii at the same dN ch /dη 1/3 . At high k T the difference in radius grows for R out , while for R long the radii for both systems become comparable. The distinct decrease of radii with k T is observed both in pp and p-Pb.
The correlation strength λ increases with k T from 0.44 to 0.58 for the collisions with highest multiplicities. It is also higher for low multiplicity collisions, with a difference of 0.1 between collisions with highest multiplicities and lowest multiplicities. A non-constant λ parameter as a function of k T is an indication of a non-Gaussian shape of the correlation function. The correlation functions are normalized to the ratio of the number of pairs in the signal and background histograms. The positive correlation at low q has to be then compensated by the normalization parameter N, which is in the range of 0.9-1.0. The χ 2 /ndf for the three-dimensional fit is on the order of 1.2.
The extracted background parameters indicate that this contribution increases with k T and decreases with multiplicity, which is consistent with qualitative expectations for the mini-jet effect. The shape of the background is not spherical, leading to finite contributions to the (2, 0) and (2, 2) components. The constant shift in these components, given by β 0 2 and β 2 2 respectively, is only significant for the (2, 0) component in lower multiplicities.
The corresponding fit results for the EGE fit are shown in Fig. 8. In the side direction the radii are consistent with the GGG results. The radii in the out and long directions are not Gaussian widths in this case and cannot be directly compared to previous fits. However, all the trends are qualitatively the same in both cases: radii increase with event multiplicity and decrease with pair transverse momentum. The values are 10% (for side and long) to 20% (for out) higher than those measured in pp collisions at similar event multiplicity [16]. The λ parameter for the EGE fit is on the order of 0.7 for the SH fits, and growing from 0.7 at low k T to approximately 0.9 at the highest k T for the Cartesian fit, therefore  significantly higher than in the Gaussian case. The observation that λ is closer to unity when moving from GGG to EGE fits is expected, as the EGE fit describes the shape of the correlation much better at low q and therefore better accounts for the non-Gaussian tails in the source function.

Model comparisons
Hydrodynamic model calculations for p-Pb collisions [24,25], shown as lines in Fig. 7, predict the existence of a collectively expanding system. Both models employ two initial transverse size assumptions: R init = 1.5 fm and R init = 0.9 fm, which correspond to two different scenarios of the energy deposition in the wounded nucleon model [25]. The resulting charged particle multiplicity densities dN ch /dη of 45 [25] and 35 [24] are equal or higher than the one in the ALICE 0-20% multiplicity class. The calculations for R out overestimate the measured radii, while the ones with large initial size strongly overpredict the radii. The scenarios with lower initial size are closer to the data. For R side , the calculations are in good agreement with the data in the highest multiplicity class, both in magnitude and in the slope of the k T dependence. Only the Shapoval et al. [24] calculation for large initial size shows higher values than data. For R long , calculations by Bożek and Broniowski [25] overshoot the measurement by at least 30% for the most central data, while those by Shapoval et al. are consistent within systematics. Again, the slope of the k T dependence is comparable. The study shows that the calculation with large initial size is disfavored by data. The calculations with lower initial size are closer to the experimental results, but are still overpredicting the overall magnitude of the radii by 10-30%. Further refinement of the initial conditions may lead to a better agreement of the models with the data, especially at large multiplicities. The slope of the k T dependence is usually interpreted as a signature of collectivity. Interestingly, it is very similar in data and the models in all directions, which suggests that the system dynamics might be correctly modeled by hydrodynamics.
Also in the data the source shape is distinctly non-Gaussian. Further studies would require examination of the source shape in p-Pb collision models to see if similar deviations from a Gaussian form are observed.
The CGC approach has provided a qualitative statement on the initial size of the system in p-Pb collisions, suggesting that it is similar to that in pp collisions [27,28]. The measured radii, at high multiplicities and low k T , are 10-20% larger than those observed at similar multiplicities in pp data. For lower multiplicities the differences are smaller. These differences could still be accommodated in CGC calculations. Furthermore, the evolution of the slope of the k T dependence is similar between pp and p-Pb collisions in the side direction. Another similarity is the distinctly non-Gaussian shape of the source, which in pp and p-Pb is better described by an exponential-Gaussian-exponential form. It appears that data in p-Pb collisions still exhibit strong similarities to results from pp collisions. However some deviations, which make the p-Pb more similar to A-A collisions, are also observed, especially at high multiplicity. The differences between small systems such as pp and p-Pb and peripheral A-A data are most naturally explained by the significantly different initial states in the two scenarios. Dedicated theoretical investigation of this issue is needed for a more definite answer, which may be able to accommodate both CGC and hydrodynamic picture.

Comparison to the world data
In Fig. 9 the results from this analysis of the p-Pb data from the LHC (red filled circles) are compared to the world heavy-ion data, including results obtained at lower collision energies, as well as to results from pp collisions from ALICE and STAR. It has been observed [13] that the 3-dimensional femtoscopic radii scale roughly with the cube root of the measured charged-particle multiplicity density not only for a single energy and collision system, but also across many collision energies and initial system sizes. The pp and A-A datasets show significantly different scaling behavior, although both are linear in dN ch /dη 1/3 .
The p-Pb radii agree with those in pp collisions at low multiplicities. With increasing multiplicity, the radii for the two systems start to diverge. An analysis of one-dimensional averaged radii in pp, p-Pb and Pb-Pb collisions using the 3-pion cumulant correlations technique reveals that the multiplicity scaling for p-Pb lies between pp and Pb-Pb trends [21], which is consistent with results presented here. On the other hand, the deviation of the correlation shape from Gaussian is similar to that observed in pp collisions, and unlike the shapes observed in A-A collisions.

Conclusions
We reported on the three-dimensional pion femtoscopic radii in p-Pb collisions at √ s NN = 5.02 TeV, measured in four multiplicity and seven pair momentum intervals. The radii are found to decrease with k T in all cases, similar to measurements in A-A and high-multiplicity pp collisions. The radii increase with event multiplicity. At low multiplicities they are comparable to the pp values, while at higher multiplicities and low pair transverse momentum they are larger by 10-20%. However, they do not reach the values observed in A-A collisions at lower energies. The high multiplicity data are compared to predictions from two models, both of them incorporating a fast hydrodynamic expansion of the created medium. They overpredict the values of the R out and R long parameters, however the introduction of a smaller initial size results in a better description. The values of the R side parameter and the slope of the k T dependence of the radii are in reasonable agreement. The models based on the CGC formalism suggest sizes similar to those obtained in pp data. The observed differences of about 10-20% for high multiplicity p-Pb collisions might not exclude this scenario. The observed non-Gaussian shape of the correlation is also similar in the pp and p-Pb collision systems.  Fig. 9: Comparison of femtoscopic radii (Gaussian), as a function of the measured charged-particle multiplicity density, measured for various collision systems and energies by CERES [42], STAR [45,46,53], PHENIX [54], and ALICE [16].