Femtoscopy of pp collisions at $\sqrt{s}$=0.9 and 7 TeV at the LHC with two-pion Bose-Einstein correlations

We report on the high statistics two-pion correlation functions from pp collisions at $\sqrt{s}=0.9$ TeV and $\sqrt{s}$=7 TeV, measured by the ALICE experiment at the Large Hadron Collider. The correlation functions as well as the extracted source radii scale with event multiplicity and pair momentum. When analyzed in the same multiplicity and pair transverse momentum range, the correlation is similar at the two collision energies. A three-dimensional femtoscopic analysis shows an increase of the emission zone with increasing event multiplicity as well as decreasing homogeneity lengths with increasing transverse momentum. The latter trend gets more pronounced as multiplicity increases. This suggests the development of space-momentum correlations, at least for collisions producing a high multiplicity of particles. We consider these trends in the context of previous femtoscopic studies in high-energy hadron and heavy-ion collisions, and discuss possible underlying physics mechanisms. Detailed analysis of the correlation reveals an exponential shape in the outward and longitudinal directions, while the sideward remains a Gaussian. This is interpreted as a result of a significant contribution of strongly decaying resonances to the emission region shape. Significant non-femtoscopic correlations are observed, and are argued to be the consequence of"mini-jet"-like structures extending to low $p_{\rm T}$. They are well reproduced by the Monte-Carlo generators and seen also in $\pi^+\pi^-$ correlations.

We report on the high statistics two-pion correlation functions from pp collisions at √ s = 0.9 TeV and √ s = 7 TeV, measured by the ALICE experiment at the Large Hadron Collider. The correlation functions as well as the extracted source radii scale with event multiplicity and pair momentum. When analyzed in the same multiplicity and pair transverse momentum range, the correlation is similar at the two collision energies. A threedimensional femtoscopic analysis shows an increase of the emission zone with increasing event multiplicity as well as decreasing homogeneity lengths with increasing transverse momentum. The latter trend gets more pronounced as multiplicity increases. This suggests the development of space-momentum correlations, at least for collisions producing a high multiplicity of particles. We consider these trends in the context of previous femtoscopic studies in high-energy hadron and heavy-ion collisions, and discuss possible underlying physics mechanisms. Detailed analysis of the correlation reveals an exponential shape in the outward and longitudinal directions, while the sideward remains a Gaussian. This is interpreted as a result of a significant contribution of strongly decaying resonances to the emission region shape. Significant non-femtoscopic correlations are observed, and are argued to be the consequence of "mini-jet"-like structures extending to low p T . They are well reproduced by the Monte-Carlo generators and seen also in π + π − correlations.

I. INTRODUCTION
Proton-proton collisions at √ s = 0.9 TeV and √ s = 7 TeV have been recorded by A Large Ion Collider Experiment (AL-ICE) at the Large Hadron Collider (LHC) at CERN in 2010. These collisions provide a unique opportunity to probe Quantum Chromodynamics (QCD) in the new energy regime. The distinguishing feature of QCD is the mechanism of color confinement, the physics of which is not fully understood, due, Goethe in part, to its theoretical intractability [1]. The confinement mechanism has a physical scale of the order of the proton radius and is especially important at low momentum. The study presented in this work aims to measure the space-time extent of the source on this scale. Two-pion correlations at low relative momentum were first shown to be sensitive to the spatial scale of the emitting source inp + p collisions by G. Goldhaber, S. Goldhaber, W. Lee and A. Pais 50 years ago [2]. Since then, they were studied in e + + e − [3], hadron-and lepton-hadron [4], and heavy ion [5] collisions. Especially in the latter case, two-particle femtoscopy has been developed into a precision tool to probe the dynamically-generated geometry structure of the emitting system. In particular, a sharp phase transition between the color-deconfined and confined states was precluded by the observation of short timescales, and femtoscopic measurement of bulk collective flow proved that a strongly self-interacting system was created in the collision [6,7].
Femtoscopy in heavy-ion collisions is believed to be understood in some detail, see e.g. [5]. The spatial scales grow naturally with the multiplicity of the event. Strong hydrodynamical collective flow in the longitudinal and transverse directions is revealed by dynamical dependencies of femtoscopic scales. The main puzzling aspect of the data is the relative energy independence of the results of the measurements.
To some extent, Bose-Einstein correlations in particle physics were initially of interest only as a source of systematic uncertainty in the determination of the W boson mass [8]. But overviews [3,4,9] of femtoscopic measurements in hadron-and lepton-induced collisions reveal systematics surprisingly similar to those mentioned above for heavyion collisions. Moreover, in the first direct comparison of femtoscopy in heavy-ion collisions at Relativistic Heavy-Ion Collider (RHIC) and proton collisions in the same apparatus an essentially identical multiplicity-and momentum-dependence was reported in the two systems [10]. However, the multiplicities at which the femtoscopic measurement in pp collisions at RHIC was made were still significantly smaller than those in even the most peripheral heavy-ion collisions. In this work we are, for the first time, able to compare femtoscopic radii measured in pp and heavy-ion collisions at comparable event multiplicities. At these multiplicities the observed correlations may be influenced by jets [11] while other studies suggest that a system behaving collectively may be created [12].
In our previous work [13] we reported that a multiplicity integrated measurement does not show any pair momentum dependence of the R inv radius measured in the Pair Rest Frame (PRF). Similar analysis from the CMS collaboration [14] also mentions that no momentum dependence was observed. However the analysis in two multiplicity ranges suggested that momentum dependence may change with multiplicity, although any strong conclusions were precluded by limited statistics. In this work we explore this dependence by using high statistics data and more multiplicity ranges. It enabled us to perform the three-dimensional analysis in the Longitudinally Co-Moving System (LCMS), where the pair momentum along the beam vanishes.
The paper is organized as follows: in Section II we describe the ALICE experimental setup and data taking conditions for the sample used in this work. In Section III we present the correlation measurement and characterize the correlation functions themselves. In Section IV A we show the main results of this work: the three-dimensional radii extracted from the data. We discuss various observed features and compare the results to other experiments. In Section V we show, for completeness, the one-dimensional R inv analysis. Finally in Section VI we summarize our results. All the numerical values can be obtained from the Durham Reaction Database [15].

II. ALICE DATA TAKING
In this study we report on the analysis of pp collisions recorded by the ALICE experiment during the 2010 run of the LHC. Approximately 8 million events, triggered by a minimum bias trigger at the injection energy of √ s = 0.9 TeV, and 100 million events with similar trigger at the maximum LHC energy to date, √ s = 7 TeV, were analyzed in this work. The ALICE Time Projection Chamber (TPC) [16] was used to record charged particle tracks as they left ionization trails in the Ne − CO 2 gas. The ionization drifts up to 2.5 m from the central electrode to the end-caps to be measured on 159 padrows, which are grouped into 18 sectors; the position at which the track crossed the padrow was determined with resolutions of 2 mm and 3 mm in the drift and transverse directions, respectively. The momentum resolution is ∼ 1% for pions with p T = 0.5 GeV/c. The ALICE Inner Tracking System (ITS) was also used for tracking. It consists of six silicon layers, two innermost Silicon Pixel Detector (SPD) layers, two Silicon Drift Detector (SDD) layers, and two outer Silicon Strip Detector (SSD) layers, which provide up to six space points for each track. The tracks used in this analysis were reconstructed using the information from both the TPC and the ITS, such tracks were also used to reconstruct the primary vertex of the collision. For details of this procedure and its efficiency see [17]. The forward scintillator detectors VZERO are placed along the beam line at +3 m and −0.9 m from the nominal interaction point. They cover a region 2.8 < η < 5.1 and −3.7 < η < −1.7 respectively. They were used in the minimum bias trigger and their timing signal was used to reject the beam-gas and beam-halo collisions.
The minimum bias trigger required a signal in either of the two VZERO counters or one of the two inner layers of the Silicon Pixel Detector (SPD). Within this sample, we selected events based on the measured charged-particle multiplicity within the pseudorapidity range |η| < 1.2. Events were required to have a primary vertex within 1 mm of the beam line, and 10 cm of the center of the 5 m-long TPC. This provides almost uniform acceptance for particles with |η| < 1 for all events in the sample. It decreases for 1.0 < |η| < 1.2. In addition, we require events to have at least one charged particle reconstructed within |η| < 1.2.
The minimum number of clusters associated to the track in the TPC is 70 (out of the maximum of 159) and 2 in the ITS (out of the maximum of 6). The quality of the track is determined by the χ 2 /N value for the Kalman fit to the re-constructed position of the TPC clusters (N is the number of clusters attached to the track); the track is rejected if the value is larger than 4.0 (2 degrees of freedom per cluster). Tracks with |η| < 1.2 are taken for the analysis. The p T of accepted particles has a lower limit of 0.13 GeV/c, because tracks with lower p T do not cross enough padrows in the TPC. The efficiency of particle reconstruction is about 50% at this lowest limit and then quickly increases and reaches a stable value of approximately 80% for p T > 0.2 GeV/c. In order to reduce the number of secondary particles in our sample, we require the track to project back to the primary interaction vertex within 0.018 + 0.035p −1.01 T cm in the transverse plane and 0.3 cm in the longitudinal direction (so-called Distance of Closest Approach or DCA selection). ALICE provides an excellent particle identification capability, through the combination of the measurement of the specific ionization (dE/dx) in the TPC and the ITS and the timing signals in the ALICE Time Of Flight (TOF). In the momentum range covered here (0.13 GeV/c to 0.7 GeV/c) pions constitute the majority of particles. We use only the TPC measurement for Particle IDentification (PID) in this work, as the other detectors offer significant improvement at higher p T than used here. This PID procedure results in a small contamination of the pion sample by electrons at p T < 0.2 GeV/c and kaons at p T > 0.65 GeV/c. Allowing other particles into our sample has only a minor effect of lowering the strength of the correlation (the λ parameter), while it does not affect the femtoscopic radius, so we do not correct for it explicitly. The amount of electron contamination is less than 5%.

III. CORRELATION FUNCTION MEASUREMENT
Experimentally, the two-particle correlation function is defined as the ratio C (q) = A (q) /B (q), where A (q) is the measured two-pion distribution of pair momentum difference q = p 2 − p 1 , and B (q) is a similar distribution formed by using pairs of particles from different events [18].
The size of the data sample used for this analysis allowed for a highly differential measurement. In order to address the physics topics mentioned in the introduction, the analysis was performed simultaneously as a function of the total event multiplicity N ch and pair transverse momentum k T = | p T,1 + p T,2 | /2. For the multiplicity determination we counted the tracks reconstructed simultaneously in the ITS and the TPC, plus the tracks reconstructed only in the ITS in case the track was outside of the TPC η acceptance. The total number of events accepted after applying the selection criteria in the √ s = 7 TeV sample was 60 × 10 6 and in the √ s = 0.9 TeV sample it was 4.42 × 10 6 . We divided the full multiplicity range into eight and four ranges for the two energies respectively in such a way that the like-charge pion pair multiplicity in each of them was comparable. Table I gives (a) values for the range of raw charged particle multiplicity that was used to categorize the event, (b) the corresponding mean charge particle density dN ch /dη as well as (c) number of events and (d) the number of identical pion pairs in each range. The femtoscopic measurement requires the events to have at  least one charged pion identified 1 and its momentum determined. We give the dN ch /dη values in Tab. I for this event sample. We denote this value as dN ch /dη | N ch ≥1 ; its typical uncertainty is 10%. We note that for the the lowest multiplicity this charged particle density is biased towards higher values with respect to the full sample of inelastic events.

A. Correlation function representations
The correlations are measured as a function of pair relative momentum four-vector q. We deal with pions, so the masses of the particles are fixed -in this case q reduces to a vector: q. The one-dimensional analysis is performed versus the magnitude of the invariant momentum difference q inv = | q|, in PRF. The large available statistics for this work allowed us to perform a detailed analysis also for the 3D functions. In forming them, we calculate the momentum difference in LCMS and decompose this q LCMS according to the Bertsch-Pratt [19,20] "out-side-long" (sometimes indicated by o, s, and l subscripts) parametrization. Here, q long is parallel to the beam, q out is parallel to the pair transverse momentum, and q side is perpendicular to q long and q out . If one wishes to compare the radii measured in LCMS to R inv one needs to multiply one of the transverse radii in LCMS (the one along the pair transverse momentum) by the Lorentz γ corresponding to the pair transverse velocity, and then average the three radii. Therefore an R inv constant with momentum is consistent with the radii in LCMS decreasing with momentum. Figure 1 shows one-dimensional projections of the 3-dimensional correlation function C(q out , q side , q long ) onto the q out , q side , and q long axes, for π + pairs from one of the multiplicity/k T ranges from the √ s = 7 TeV sample. The function is normalized with a factor that is a result of the fit (the details of the procedure are described in Sec. III D); unity means no correlation.
The 1-dimensional projections, shown in Fig. 1, present a limited view of the 3-dimensional structure of the correlation function. It is increasingly common to represent correlation functions in a harmonic analysis [21][22][23]; this provides a more complete representation of the 3-dimensional structure of the correlation, a better diagnostic of non-femtoscopic correlations [22], and a more direct relation to the shape of the source [24]. The moments of the Spherical Harmonic (SH) decomposition are given by Here, the out-side-long space is mapped onto Euler angles in which q long = | q| cos θ and q out = | q| sin θ cos φ. For pairs of identical particles in collider experiments done with symmetrical beams, including the analysis in this work, the odd l and the imaginary and odd m components for even l vanish. The first three non-vanishing moments, which capture essentially all of the 3-dimensional structure, are then C 0 0 , C 0 2 , and C 2 2 . These are shown in Fig. 2. The components for l ≥ 4 represent the fine details of the correlation structure and are not analyzed in this work. The C 0 0 is the angle-averaged component. It captures the general shape of the correlation. The width of the peak near q = 0 is inversely proportional to the overall femtoscopic size of the system. The C 0 2 component is the correlation weighed with the cos 2 (θ). If it differs from 0, it signifies that the longitudinal and transverse sizes of the emission region differ. The C 2 2 is weighed with cos 2 (φ). If it differs from 0, it signals that the outward and sideward sizes differ. The correlation function is normalized to the number of pairs in the background divided by the number of pairs in the signal.

B. Measured correlations
In Figs. 1 and 2 we show selected correlations to illustrate how they depend on multiplicity. This is done for k T of (0.2, 0.3) GeV/c; the behavior in other k T ranges and at the lower collision energy is qualitatively the same. The narrowing of the correlation peak with increasing multiplicity is apparent, corresponding to the increase of the size of the emitting region. The behavior of the correlation function at large q is also changing, the low multiplicity baseline is not flat, goes below 1.0 around q = 1 GeV/c and then rises again at larger q, for higher multiplicities the background becomes flatter at large q. In Cartesian representation shown in Fig. 1, areas with no data points (acceptance holes) are seen in q out projection near q = 0.5 GeV/c and in q long above 0.6 GeV/c. Since q long is proportional to the difference of longitudinal momenta, its value is limited due to η acceptance. In the out direction the hole appears due to a combination of lower p T cut-off and the selected k T range. It can be simply understood as follows: For the projection in the upper panel of Fig. 1, we take the value of q side and q long small. The value of q side is proportional to the azimuthal angle difference, while q long is proportional to polar angle difference. For q side , q long = 0, q out is simply p T,2 − p T,1 and k T is (p T,1 + p T,2 )/2, where p T is no longer a 2-vector, but just a scalar. The particles are either fully aligned (both p T 's are positive or both are negative) or back-to-back (one p T is positive, the other negative). When we combine the lower p T cut-off |p T | > 0.13 GeV/c and the k T selection 0.2 ≤ k T ≤ 0.3, it can be shown that some range of the q out values is excluded. This range will depend on the k T selection.
The k T dependence of the correlation function is shown in Figs. 3 and 4, for multiplicity 17 ≤ N ch ≤ 22. The behavior in other multiplicity ranges and at lower energy is qualitatively similar (except the lowest multiplicity bin where the behavior is more complicated -see the discussion of the extracted radii in Sec. III D for details). We see a strong change of the correlation with k T , with two apparent effects. At low k T the correlation appears to be dominated by the femtoscopic effect at q < 0.3 GeV/c, and is flat at larger q. As k T grows, the femtoscopic peak broadens (corresponding to a decrease in size of the emitting region). In addition, a wide structure, extending up to 1.0 GeV/c in q for the highest k T range, appears. We analyze this structure in further detail later in this work. We also see that, according to expectations, the acceptance holes in the out and long region move as we change the k T range. Figure 5 shows the example of the correlation function, for the same multiplicity/k T range, for pp collisions at two collision energies. We note a similarity between the two functions; the same is seen for other k T 's and overlapping mul- pp @ 0.9 TeV pp @ 7 TeV tiplicity ranges. The similarity is not trivial: changing the multiplicity by 50%, as seen in Fig. 2 or k T by 30% as seen in Fig. 3 has a stronger influence on the correlation function than changing the collision energy by an order of magnitude. We conclude that the main scaling variables for the correlation function are global event multiplicity and transverse momentum of the pair; the dependence on collision energy is small. The energy independence of the emission region size is the first important physics result of this work. We emphasize that it can be already drawn from the analysis of the correlation functions themselves, but we will also perform more qualitative checks and discussions when we report the fitted emission region sizes in Section IV.

C. Non-femtoscopic correlation structures
In Fig. 3 we noted the appearance of long-range structures in the correlation functions for large k T . If these were of femtoscopic origin, they would correspond to an unusually small emission region size of 0.2 fm. We reported the observation of these structures in our previous analysis [13] at √ s = 0.9 TeV, where they were interpreted as non-femtoscopic correlations coming from "mini-jet" like structures at p T < 1 GeV/c. Here we further analyze this hypothesis. In Fig. 6 we show the comparison of the correlation function at multiplicity 12 ≤ N ch ≤ 16 in an intermediate k T range, where the long-range correlations are apparent, to the Monte-Carlo (MC) calculation. The simulation used the PYTHIA generator [25], Perugia-0 tune [26] as input and was propagated through the full simulation of the ALICE detector [16]. Then it was reconstructed and analyzed in exactly the same way as our real data, us-ing the same multiplicity and k T ranges. The MC calculation does not include the wave-function symmetrization for identical particles; hence, the absence of the femtoscopic peak at low q is expected. In the angle-averaged C 0 0 component a significant correlation structure is seen, up to 1 GeV/c, with a slope similar to the data outside of the peak at low q. Similarly, in the C 0 2 component a weak and wide correlation dip is seen around q = 0.5 GeV/c, which is also seen in the data. In MC, the correlation in C 0 2 disappears at lower q, while for the data it extends to much lower q, exactly where the femtoscopic peak is expected and seen in C 0 0 . Our hypothesis is that both the long-range peak in C 0 0 and the dip in C 0 2 are of a "mini-jet" origin. They need to be taken into account when fitting the correlation function from data, so that the femtoscopic peak can be properly extracted and characterized. The calculations was also carried out with a second Monte-Carlo, the PHOJET model [27,28], and gave similar results. The differences between the two models are reflected in the systematic error.
In order to characterize the non-femtoscopic background we study in detail the correlation structure in the MC generators, in exactly the same multiplicity/k T ranges as used for data analysis. We see trends that are consistent with the "mini-jet" hypothesis. The correlation is small or non-existent for low p T (first k T range) and it grows strongly with p T . In Fig. 7 we show this structure for selected multiplicity/k T at both energies. At the highest k T the effect has the magnitude of 0.3 at low q, comparable to the height of the femtoscopic peak. The appearance of these correlations is the main limiting factor in the analysis of the k T dependence. We tried to analyze the correlations at k T higher than 0.7 GeV/c but we were unable to obtain a meaningful femtoscopic result, because the "minijet" structure was dominating the correlation. The strength of the correlation decreases with growing multiplicity (as expected), slower than 1/M, so that it is still significant at the highest multiplicity. We studied other tunes of the PYTHIA model and found that the Perugia-0 tune reproduces the "minijet" correlation structures best, which is why it is our choice. Its limitation though is a relatively small multiplicity reach, smaller than the one observed in data. As a result the MC calculation for our highest multiplicity range is less reliablethis is reflected in the systematic error.
Analyzing the shape of the underlying event correlation for identical particle pairs in MC is important; however, it does not ensure that the behavior of the correlation at very low q is reproduced well in MC. We compared the identical particle MC and data in the large q region, where the femtoscopic effect is expected to disappear, and found them to be very similar in all multiplicity/k T . However, if the "mini-jet" hypothesis is correct, the same phenomenon causes similar correlations to appear for non-identical pions. The magnitude is expected to be higher than for identical pions, because it is easier to produce an oppositely-charged pair from a fragmenting "mini-jet" than it is to create an identically-charged pair, due to local charge conservation. Moreover, the femtoscopic effect for non-identical pions comes from the Coulomb interaction only. It is limited to very low q, below 0.1 GeV/c. It is therefore possible to test the low-q behavior of the "mini-jet" correlation with such correlations. In Fig. 8 we show the measured π + π − correlation functions, in selected multiplicity/k T ranges, compared to the corresponding correlations from the same MC sample which was used to produce correlations in Fig. 7. The underlying event long-range correlation is well reproduced in the MC. We see some deviation in the lowest multiplicity range, which is taken into account in the systematic error estimation. At larger multiplicities the strength of the correlation is well reproduced. By comparing the 3D function in SH we checked that the shape in 3D q space is also in agreement between data and MC. The magnitude for non-identical pions is slightly bigger than for identical pions, as expected. The femtoscopic Coulomb effect at q < 0.1 GeV/c is also visible. Another strong effect, even dominating at low multiplicity, are the peaks produced by the correlated pairs of pions coming from strong resonance decays. They do appear in the MC as well, but they are shifted and have different magnitude. This is the effect of the simplified treatment of resonance decays in PYTHIA, where phase space and final state rescattering are not taken into account. By analyzing some of the correlation functions in Fig 8 we were able to identify signals from at least the following decays: two-body ρ, f 0 , and f 2 mesons decays, three-body ω meson decay, and also possibly η meson two-body decay. Some residual K 0 S weak decay pairs, which are not removed by our DCA selection, can also be seen. All of these contribute through the full q range (0.0, 1.2) GeV/c. This fact, in addition to the stronger "mini-jet" contribution to non-identical (as compared to identical) correlations, makes the non-identical correlation not suitable for the background estimation for identical pion pairs. We also note that there appears to be very rich physics content in the analysis of resonances decaying strongly in the π + π − channel; however, we leave the investigation of this topic for separate studies.
The study of the π + π − correlations confirms that the MC generator of choice reproduces the underlying event structures also at low q. We found that they are adequately described by a Gaussian in LCMS for the C 0 0 component. The dashed lines in Fig. 7 show the fit of this form to the correlation in MC. The results of this fit, taken bin-by-bin for all multiplicity/k T ranges, are the input to the fitting procedure described in Section III D. Similarly, the observed C 0 2 correlation can be characterized well by a Gaussian, with the magnitude of −0.01 or less and a peak around q = 0.5 GeV/c with a width of 0.25 − 0.5 GeV/c. We proceed in the same way as for C 0 0 ; we fit the MC correlation structures with this functional form and take the results as fixed input parameters in the fitting of the measured correlations.

D. Fitting the correlation function
Having qualitatively analyzed the correlation functions themselves we move to the quantitative analysis. The femtoscopic part of the correlation function is defined theoretically via the Koonin-Pratt equation [29,30]: where q is the pair 3-momentum difference (the fourth component is not independent for pairs of identical pions since masses of particles are fixed), k is the pair total momentum, r is the pair space-time separation at the time when the second particle undergoes its last interaction, Ψ is the wave function of the pair, and S is the pair separation distribution. The aim in the quantitative analysis of the correlation function is to learn as much as possible about S from the analysis of the measured C. The correlation function C is, in the most general form, a 6-dimensional object. We reduce the dimensionality to 3 by factorizing out the pair momentum k. We do not study the dependence on the longitudinal component of k in this work. The dependence on the transverse component of k is studied via the k T binning, introduced in Section II. We assume that S is independent of k inside each of the k T ranges. We also note that for identical pions the emission function S is a convolution of two identical single particle emission functions S 1 .
In order to perform the integral in Eq. (2) we must postulate the functional form of S or S 1 . We assume that S does not depend on q. The first analysis is performed with S 1 as a 3-dimensional ellipsoid with Gaussian density profile. This produces S which is also a Gaussian (with σ larger by a factor of √ 2): where R G out , R G side and R G long are pion femtoscopic radii, also known as "HBT radii" or "homogeneity lengths", and r o , r s and r l are components of the pair separation vector. For identical charged pions Ψ should take into account the proper symmetrization, as well as Coulomb and strong interaction in the final state. In the case of the analysis shown in this work, with pions emitted from a region with the expected size not larger than 2 − 3 fm, the strong interaction contribution is relatively small and can be neglected [31]. The influence of the Coulomb interaction is approximated with the Bowler-Sinyukov method. It assumes that the Coulomb part can be factorized out from Ψ and integrated independently. There are well-known limitations to this approximation but they have minor influence for the analysis shown in this work. With these assumptions Ψ is a sum of two plane waves modified by a proper symmetrization. By putting Eq. (3) into Eq. (2) the integration can be done analytically and yields the quantum statistics-only correlation C qs : where λ is the fraction of correlated pairs for which both pions were correctly identified. The 3-dimensional correlation function is then modified with the Bowler-Sinyukov formula to obtain the complete femtoscopic component of the correlation C f : where K is the Coulomb like-sign pion pair wave function squared averaged over the Gaussian source with a radius of 1 fm. Changing this radius within the range of values measured in this work has negligible effect on the extracted radii. Eq. (5) describes properly the femtoscopic part of the two pion correlation function. However, in the previous section we have shown that our experimental functions also contain other, non-femtoscopic correlations. We studied them in all multiplicity/k T ranges and found that they can be generally described by a combination of an angle-averaged Gaussian in LCMS plus a small Gaussian deviation in the C 0 2 component:

B( q LCMS
where A h , A w , B h , B m and B w are parameters. They are obtained, bin-by-bin, from the fit to the MC simulated correlation functions shown in Fig. 7. They are fixed in the procedure of fitting the data. The final functional form that is used for fitting is then: where N is the overall normalization. Projections of the Cartesian representation of the correlation functions, shown in Figs. 1 and 4, are normalized with this factor. Function (7) is used to fit both the SH and Cartesian representation of the 3D correlations.
In Fig. 9 an example of the fit to one of our correlation functions is shown. The SH representation of the data is shown as points; the result of the fit is a black dashed line. The femtoscopic component is shown as a blue dotted line, the non-femtoscopic background as green dash-dotted line. The correlation function in this range has significant contribution from the background and is reasonably reproduced by the fit. At q < 0.1 GeV/c the fit misses the data points in C 0 0 and C 2 2 ; we will discuss what can be done to improve the agreement in Section IV C. In Fig. 10 the same correlation is shown as projections of the 3D Cartesian representation. The other q components are integrated over the range of 0 − 0.16 GeV/c. The fit, shown as lines, is similarly projected. In this plot the fit does not describe the shape of the correlation perfectly; however, the width is reasonably reproduced.

A. Results of the 3D Gaussian fits
We fitted all 72 correlation functions (4+8 multiplicity ranges for two energies times 6 k T ranges) with Eq. (7). We show the resulting femtoscopic radii in Fig. 11 as a function of k T . The strength of the correlation λ is relatively independent of k T , is 0.55 for the lowest multiplicity, decreases monotonically with multiplicity and reaches the value of 0.42 for the highest multiplicity range. The radii shown in the Fig. 11 are the main results of this work. Let us now discuss many aspects of the data visible in this figure.
Firstly, the comparison between the radii for two energies, in the same multiplicity/k T ranges reveals that they are universally similar, at all multiplicities, all k T 's and all directions. This confirms what we have already seen directly in the measured correlation functions. The comparison to √ s = 200 GeV pp collisions at RHIC is complicated by the fact that these data are not available in multiplicity ranges. The multiplicity reach at RHIC corresponds to a combination  11. Parameters of the 3D Gaussian fits to the complete set of the correlation functions in 8 ranges in multiplicity and 6 in k T for pp collisions at √ s = 7 TeV, and 4 ranges in multiplicity and 6 in k T for pp collisions at √ s = 0.9 TeV. All points at given k T bin should be at the same value of k T , but we shifted them to improve visibility. Open black squares show values for pp collisions at √ s = 200 GeV from STAR [10]. Lines connecting the points for lowest and highest multiplicity range were added to highlight the trends. of the first three multiplicity ranges in our study. No strong change is seen between the RHIC and LHC energies. It shows that the space-time characteristics of the soft particle production in pp collisions are only weakly dependent on collision energy in the range between 0.9 TeV to 7 TeV, if viewed in narrow multiplicity/k T ranges. Obviously the √ s = 7 TeV data have a higher multiplicity reach, so the minimum-bias (multiplicity/k T integrated) correlation function for the two energies is different.
Secondly, we analyze the slope of the k T dependence. R G long falls with k T at all multiplicities and both energies. R G out and R G side show an interesting behavior -at low multiplicity the k T dependence is flat for R G side and for R G out it rises at low k T and then falls again. For higher multiplicities both transverse radii develop a negative slope as multiplicity increases. At high multiplicity the slope is bigger for R G out , while R G side grows universally at all k T 's while developing a smaller negative slope. The difference in the evolution of shapes of R G out and R G side is best seen in their ratio, shown in panel d) of Fig. 11. At low multiplicities the ratio is close to 1.0, then it decreases monotonically with multiplicity. We note that a negative slope in R G out and R G side was universally observed in all heavy-ion measurements at RHIC energies and sometimes also at lower energies. It is interpreted as a signature of a strong collective behavior of matter created in such collisions. The observation of the development, with increasing multiplicity, of such slope in pp collisions is consistent with the hypothesis, that the larger the produced multiplicity, the more self-interacting and collective is the produced system.
Thirdly, all the measured radii grow with event multiplicity, in each k T range separately. This is shown more clearly in Fig. 12, where we plot the radii as a function of dN ch /dη (1/3) (For our pp data we use the dN ch /dη (1/3) | N ch ≥1 given in Tab. I). R G side and R G long grow linearly with the cube root of √ s = 0.9 TeV and 7 TeV, compared to the results from (heavy-)ion collisions at RHIC [32,33] and SPS [34]. Panel a) shows R G out , b) shows R G side , c) shows R G long . All results are for k T = 0.4 GeV/c, except the values from the PHENIX experiment, which are at k T = 0.45 GeV/c. charged particle multiplicity, for all k T ranges. Data, at both energies, follow the same scaling. For R G out the situation is similar for medium k T ranges. The lowest k T points show the strongest growth with multiplicity, while the highest hardly grows at all. That is the result of the strong change of the slope of k T dependence with multiplicity, noted in the discussion of Fig. 11.
Similar multiplicity scaling was observed in heavy-ion collisions at RHIC energies and below. In Fig. 13 we compare our results to heavy-ion results from collision energies above 15 AGeV. This is the first time that one can directly compare pp and heavy-ion radii at the same dN ch /dη , as we measure dN ch /dη comparable to the one in peripheral AuAu and CuCu collisions at RHIC. Since the value of the radius strongly depends on k T , we carefully selected the results to have the same average k T = 0.4 GeV/c. The picture at other k T 's is qualitatively similar. While both the heavy-ion and pp data scale linearly with dN ch /dη (1/3) , the slope of the dependence is clearly different, for all directions. The pp results are systematically below the heavy-ion ones at similar multiplicity; therefore, the "universal" multiplicity scaling [5] observed in heavy-ion collisions does not hold for pp collisions at √ s = 0.9 and 7 TeV. The pp radii do scale linearly with multiplicity, but with a different slope.
We speculate that the difference comes from a different way that the two types of collisions arrive at similar multiplicity.
To produce a large number of particles in pp collision one needs a particularly energetic elementary collision that produces a lot of soft particles. The region where they are created is on the order of the incoming proton size and the growth of the size with multiplicity comes from further reinteraction between particles after they are born. In contrast, in heavy-ion collision we have a combination of many elementary nucleon scatterings, each of them producing a relatively low multiplicity. But each scattering happens at a different space-time point and it is the distribution of these points that mostly determines the final observed size. In this picture, one would expect the heavy-ion sizes to be larger than the ones observed in pp at the same multiplicity.

B. Systematic uncertainty
The correlation function is, to the first order, independent of the single particle acceptance and efficiency. We performed the analysis independently for many samples of data, that naturally had single particle efficiencies different by up to 5%. We analyzed positive and negative pions separately, data at two magnetic field polarities, data from three different monthlong "LHC periods", each of them having a slightly different detector setup. Two-particle correlations from all these analyzes were consistent within statistical errors.
We studied the effect of momentum resolution on the correlation peak with the MC simulation of our detector. At this low p T , below 1 GeV/c, the momentum resolution for tracks reconstructed in the TPC is below 1%. This was confirmed by several methods, including the reconstruction of tracks from cosmic rays, and comparison of the reconstructed K 0 S mass peak position with the expected value. The smearing of single particle momenta does result in the smearing of the correlation peak: it makes it appear smaller and wider. We estimated that this changes the reconstructed radius by 1% for the femtoscopic size of 1 fm; the effect grows to 4% for the size of 2 fm, as it corresponds to a narrower correlation peak.
In contrast to single particle acceptance, the femtoscopic correlation function is sensitive to the two-track reconstruction effects, usually called "splitting" and "merging". The "splitting" occurs when one track is mistakenly reconstructed as two. Both tracks have then very close momenta. This results in a sharp correlation peak at low relative momentum. We have seen such effects in the data and we took several steps to remove them. Firstly, the requirement that the track is simultaneously reconstructed in the TPC and ITS decreases splitting significantly. In addition, each cluster in the TPC is flagged as "shared" if it is used in the reconstruction of more than one track. The split tracks tend to produce pairs which share most of their clusters; therefore, we removed pairs that share more than 5% of the TPC clusters. We also look for configurations where a single track is split in two segments in the TPC, e.g. by the TPC central membrane or a TPC sector boundary. Such segments should be correctly connected in the tracking procedure to form a single track if the detector calibration is perfect. However, in a few rare cases this does TABLE II. Systematic uncertainty coming from varying "mini-jet" background height/width by 10%/5% up/down. not happen and a split track can appear. Such pairs would consist of two tracks that have a relatively small number of TPC clusters and they would rarely both have a cluster in the same TPC padrow. Therefore, we count, for each pair, the number of times that both tracks have a separate (non-shared) cluster in a TPC padrow. Pairs for which this number is low are removed. Both selections are applied in the same way to the signal and background distributions. As a consequence, the fake low-momentum pairs from splitting are almost completely removed, and the remaining ones are concentrated in a very narrow relative momentum q range, corresponding essentially to the first correlation function bin. The inclusion of this bin has a negligible effect on the fitting result; hence, we do not assign any systematic error on the fitting values from these procedures. Another two track effect is merging, where two distinct tracks are reconstructed as one, due to finite detector spacepoint resolution. The ALICE detector was specifically designed to cope with the track densities expected in heavyion Pb+Pb collisions, which are expected to be orders of magnitude higher than the ones measured in pp collisions. More specifically, the ITS detector granularity as well as TPC tracking procedure, which allows for cluster sharing between tracks, make merging unlikely. We confirmed with detailed MC simulation of our detector setup that merging, if it appears at all, would only affect the correlation function in the lowest q bin, which means that it would not affect the measured radii.
In summary, the systematic uncertainty on the raw measurement, the correlation functions itself, is small.
The most significant systematic uncertainty on the extracted radii comes from the fact that we rely on the MC simulation of the "mini-jet" underlying event correlations. We fix the parameters of the B function in Eq. (7) by fitting it to the correlations obtained from the MC generated events. We confirmed with the analysis of the non-identical π + π − pairs that our Monte-Carlos of choice, the Perugia-0 tune of the PYTHIA 6 model, and the PHOJET model reproduce the height and the width of the "mini-jet" peak with an accuracy better than 10%, except the first multiplicity range where the differences go up to 20% for the highest k T range. We performed the fits to the correlation function varying the parameters A h and B h of the B function by ±10%, and A w by 5%. The fit values for the case when A h , B h are decreased and A w is increased (corresponding to smaller "mini-jet" correlations) are systematically below the standard values. For larger "mini-jet" correlations  Fig. 3. In an ideal case both procedures should produce identical fit results; therefore, we take the difference between the radii obtained from the two procedures as an estimate of the systematic uncertainty incurred by the fitting procedure itself. The error is shown in Table V as a function of k T . The large error at low k T is coming from the fact that the two procedures are sensitive to the holes in the acceptance in a different way. It reflects the experimental fact that, in these k T ranges, pairs in certain kinematic regions are not measured; therefore, the femtoscopic radius cannot be obtained with better accuracy. In the highest k T range the "mini-jet" underlying correlation is highest and broadest. If our simple phenomenological parametrization of it does not perfectly describe its behavior in full 3D space, it can affect differently a fit in Cartesian and SH representations.
In summary, the combined systematic error is 10% for all k T and multiplicity ranges except the ones at the lower and upper edge. It is 20% for the lowest and highest k T and for the lowest and highest multiplicity range at each collision energy. It is also never smaller than 0.1 fm.

C. Non-Gaussian fits
In the discussion of Fig. 9 we note that the measured correlation function is not perfectly reproduced by a 3D Gaussian fit. In our previous work [13] and in the work of the CMS collaboration [14] it was noted that the shape of the 1dimensional correlation in the Pair Rest Frame is better described by an exponential shape. Also, model studies [35] suggest that pion production at these energies has large contribution from strongly decaying resonances. This is confirmed by the observation of significant resonance peaks in the π + π − correlation functions, seen e.g. in Fig. 8. Resonances decay after random time governed by the exponential decay law, which transforms into an exponential shape in space via the pair velocity. By definition pair velocity exists in the out and long direction, and vanishes in side. It is then reasonable to attempt to fit the correlation with a functional form other than a simple Gaussian, at least for the out and long components.
If we keep the assumption that the emission function factorizes into the out, side and long directions, we can write a general form of the pair emission function: We can independently change each component. We stress, however, that only for a Gaussian there is an analytically known correspondence between the pair separation distribution S and single particle emission function S 1 . Two commonly used forms of S are exponential and Lorentzian. They have the desired feature that the integration in Eq. (2) can be analytically carried out and produce a Lorentzian and exponential in C respectively. In order to select the proper combination of functional forms we seek guidance from models. They suggest that at least in the out and long direction the emission function is not Gaussian, and in some cases seems to be well described by a Lorentzian. We performed a study of all 27 combinations of the fitting functions for selected multiplicity/k T ranges. We found that universally the out correlation function was best described by an exponential, corresponding to Lorentzian emission function, which agrees with model expectations. In contrast, the side direction is equally well described by a Gaussian or a Lorentzian: we chose the former because the Lorentzian correlation function would correspond to exponential pair emission function with a sharp peak at 0. We deem this unlikely, given that the models do not produce such shapes. In long, the correlation function is not Gaussian; hence, we chose the exponential shape in C for the fit. In conclusion, we postulate that the source has the following shape: which corresponds to the following form of the femtoscopic part of the correlation function formula: In Figs. 14 and 15 we show an example of the exponential-Gaussian-exponential fit to the correlation functions at multiplicity 23 ≤ N ch ≤ 29 and k T in (0.3, 0.4) GeV/c. In the SH representation we see improvements over the Gaussian fit from Fig. 9. The behavior in C 0 0 at low q is now well described. In C 2 2 the "wiggle" in the correlation is also reproduced -this is possible because the functional forms for the out and side directions are now different. In the Cartesian projections the improvement is also seen, however it is not illustrated as clearly as in the SH.
We then proceed with the fitting of the full set of 72 correlation functions. The resulting fit parameters are summarized in Fig. 16. The quality of the fit (judged by the value of χ 2 /N do f ) is better than for the 3D Gaussian fit. The λ parameter is higher by up to 0.2, as compared to the pure Gaussian fit, reflecting the fact that the new functional form accounts for the pairs contributing to the narrow correlation peak at small q. The resulting exponential radii cannot be directly compared in magnitude to the Gaussian radii from other experiments. However all the features seen in dependencies of the Gaussian radii on multiplicity and k T are also visible here. This confirms that with a functional form that fits our correlation   (10)) as a function of pair momentum k T for all multiplicity ranges and for two collision energies. Panel a) shows R E out , panel b) shows R G side , panel c) shows R E long , panel d) shows R E out /R G side ratio. All points at given k T bin should be at the same value of k T , but we shifted them to improve visibility.
valid. The study of the fit functional form shows that the correlation does not have a Gaussian shape in out and long.
The R G side from this fit should be equal to the R G side from the 3D Gaussian fit with two caveats. The first is the assumption that the emission function fully factorizes into separate functions for out, side, and long directions. In the fitting of the 3D correlation functions the residual correlation between the value of the λ parameter and the values of the radii is often observed. We noted already that the non-gaussian fit produces larger values of λ, so R G side could be affected. Nevertheless we observe very good agreement (within statistical errors for multiplicities above 16 radius for all multiplicity and k T ranges for the √ s = 0.9 TeV data. The points for different multiplicities were slightly shifted in k T for clarity. The systematic error, typically on the order of 10% is not shown [15]. Closed and open stars show the previously published result from [13] for two ranges of the multiplicity M. panel d) of Fig. 16. Again, the picture seen for the Gaussian radii is confirmed; the higher the multiplicity of the collision and the collision energy, the lower the value of the ratio.

V. FITTING 1D CORRELATIONS
For completeness, we also repeated the 1-dimensional study in Pair Rest Frame, using all the methods and fitting functions described in the previous work of ALICE [13]. The 1-dimensional correlation functions are fit with the standard Gaussian form, modified with the approximate Bowler-Sinyukov formula to account for the Coulomb interaction between charged pions: where K is the Coulomb function averaged over a spherical source of the size 1.0 fm, R inv is the femtoscopic radius and B is the function describing the non-femtoscopic background. In Fig. 17 we plot the gaussian 1-dimensional invariant radius as a function of multiplicity and k T . The closed and open stars are the results from our earlier work, which are consistent with the more precise results from this analysis. The systematic error is on the order of 10% and is now dominating the precision of the measurement. At √ s = 0.9 TeV we see that, for the lowest multiplicity, the radius is not falling with k T , while it develops a slope as one goes to higher multiplicity. The 1-dimensional analysis is consistent with the 3D measurement -one needs to take into account that by going from LCMS (3D measurement) to PRF (1D measurement) because it is necessary to boost the out radius by pair velocity, which is 1-dimensional R inv radius versus multiplicity and k T for the √ s = 7 TeV data. The points for different multiplicities were slightly shifted in k T for clarity. The systematic error, typically on the order of 10% is not shown [15]. defined by k T . Then, one averages the radii in three directions to obtain the 1D R inv .
In Fig. 18 we show the same analysis performed for the √ s = 7 TeV data. The radii are again comparable at the same multiplicity/k T range. In addition, as one goes to higher multiplicities, the k T dependence of R inv is getting more pronounced. The results are again consistent with the 3D analysis.

VI. SUMMARY
In summary, ALICE measured two-pion correlation functions in pp collisions at √ s = 0.9 TeV and at √ s = 7 TeV at the LHC. The analysis was performed in multiplicity and pair transverse momentum ranges. When viewed in the same multiplicity and pair momentum range, correlation functions at the two collision energies are similar.
The correlations are analyzed quantitatively by extracting the emission source sizes in three dimensions: outward, sideward and longitudinal. The longitudinal size shows expected behavior. It decreases with pair momentum and increases with event multiplicity, consistent with all previous measurements in elementary and heavy-ion collisions. The transverse sizes show more complicated behavior. The sideward radius grows with multiplicity and has a negative correlation with pair momentum. The outward radius at the lowest multiplicity is small for the lowest k T , increases for larger k T and then decreases. As the multiplicity grows the shape of the k T dependence gradually changes to the one monotonically falling with k T . The resulting ratio of outward to sideward radii gets smaller as multiplicity grows. Similar dependencies in heavy-ion collisions were interpreted as signatures of the collective behavior of matter. One possible interpretation of the results in this work is that as one moves towards pp collisions producing high multiplicity of particles, similar collectivity develops. More experimental and theoretical information is needed to address this intriguing possibility.
The upper range of multiplicities produced in pp collisions at √ s = 7 TeV is comparable to the multiplicities measured in peripheral heavy-ion collisions at RHIC. When plotted versus dN ch /dη (1/3) the radii in pp show linear scaling, but with different slope and offset than those observed in heavy-ion collisions. Therefore our observations violate the "universal" dN ch /dη (1/3) scaling. This proves that the initial geometry of the collision does influence the final measured radii and must be taken into account in any scaling arguments.
The analysis is complicated by the existence of the longrange underlying event correlations. We assume these are the "mini-jet" structures which are visible at values of p T as low as 0.5 GeV/c. The Monte-Carlo studies are consistent with such a hypothesis and are used to parametrize and take into account the influence of mini-jets on the fitted femtoscopic radii. Studies of the π + π − correlations are also consistent with such hypothesis. Nevertheless, the need to account for this effect remains the main source of the systematic error.
Finally, the detailed analysis of the correlation reveals that the three-dimensional Gaussian describes the measurement only approximately. A better shape, exponential-Gaussianexponential, is postulated, based on Monte-Carlo studies, and is found to better agree with the data. The resulting radii and their behavior versus event multiplicity and pair momentum are fully consistent with the one obtained with the Gaussian approximation.

VII. ACKNOWLEDGEMENTS
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the