Study of Bose-Einstein correlations in pp, pPb, and PbPb collisions at $\sqrt{s_\mathrm{NN}} =$ 0.9-7 TeV

Quantum statistical (Bose-Einstein) two-particle correlations are measured in pp collisions at $\sqrt{s}=$ 0.9, 2.76, and 7 TeV, as well as in pPb and peripheral PbPb collisions at nucleon-nucleon center-of-mass energies of 5.02 and 2.76 TeV, respectively, using the CMS detector at the LHC. Separate analyses are performed for same-sign unidentified charged particles as well as for same-sign pions and kaons identified via their energy loss in the silicon tracker. The characteristics of the one-, two-, and three-dimensional correlation functions are studied as functions of the pair average transverse momentum ($k_\mathrm{T}$) and the charged-particle multiplicity in the event. For all systems, the extracted correlation radii steadily increase with the event multiplicity, and decrease with increasing $k_\mathrm{T}$. The radii are in the range 1-5 fm, the largest values corresponding to very high multiplicity pPb interactions and to peripheral PbPb collisions with multiplicities similar to those seen in pPb data. It is also observed that the dependencies of the radii on multiplicity and $k_\mathrm{T}$ largely factorize. At the same multiplicity, the radii are relatively independent of the colliding system and center-of-mass energy.


Introduction
Studies of quantum-statistical correlations of pairs of identical particles produced in highenergy collisions provide valuable information about the size and shape of the underlying emitting system. The general technique is similar to the intensity interference mechanism used by Hanbury-Brown and Twiss (HBT effect) [1][2][3] for estimating angular dimensions of stars. In high-energy collisions, the equivalent of this astronomical effect in the femtoscopic realm was discovered in antiproton-proton collisions at √ s = 1.05 GeV [4]. The effect is known as Bose-Einstein correlations (BEC) when dealing with bosonic pairs, or femtoscopy since the characteristic probed lengths are in the femtometer range. It relates the joint probability of observing two identical particles to the product of the isolated probabilities of observing each one independently. The result is a correlation function in terms of the relative momentum of the particles in the pair, reflecting the so-called length of homogeneity of the particle-emitting source.
A broad investigation of correlations of like-sign charged particles as a function of the invariant four-momentum difference of the particles was performed by CMS using pp collisions at √ s = 0.9 TeV [5,6], 2.36 TeV [5], and 7 TeV [6]. The present paper extends the investigation of such femtoscopic correlations using two different analysis methods. First, correlations in pp collisions at √ s = 2.76 and 7 TeV are extracted using the "double ratio" technique [5,6] for unidentified pairs of same-sign charged particles with respect to different components of the relative three-momentum of the pair, thereby allowing the exploration of the source extent in various spatial directions. This procedure has the advantage of suppressing non-BEC effects coming from multi-body resonance decays, mini-jets, and energy-momentum conservation with the help of collision events simulated without Bose-Einstein correlations. Secondly, BEC effects are also studied using pairs of identical charged pions and kaons, identified via their energy loss in the CMS silicon tracker, in pp collisions at √ s = 0.9, 2.76, and 7 TeV, pPb collisions at a nucleon-nucleon center-of-mass energy of √ s NN = 5.02 TeV, and peripheral PbPb collisions at √ s NN = 2.76 TeV. The suppression of non-BEC contributions is less model-dependent with this "particle identification and cluster subtraction" approach, which also has different systematic uncertainties than the double ratio method. In both cases, the characteristics of the one-, two-, and three-dimensional correlation functions are investigated as functions of pair average transverse momentum, k T , and charged-particle multiplicity in the pseudorapidity range |η| < 2.4.
The paper is organized as follows. The CMS detector is introduced in Section 2, while track selections and particle identification are detailed in Section 3. In Section 4, the two analysis methods are described. A compilation of the results is presented in Sections 5 and 6. Finally, Section 7 summarizes and discusses the conclusions of this study.

The CMS detector
A detailed description of the CMS detector can be found in Ref. [7]. The CMS experiment uses a right-handed coordinate system, with the origin at the nominal interaction point (IP) and the z axis along the counterclockwise-beam direction. The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter. Within the 3.8 T field volume are the silicon pixel and strip tracker, the crystal electromagnetic calorimeter, and the brass and scintillator hadronic calorimeter. The tracker measures charged particles within the |η| < 2.4 range. It has 1440 silicon pixel and 15 148 silicon strip detector modules, ordered in 13 tracking layers. In addition to the electromagnetic and hadron calorimeters composed of a barrel and two endcap sections, CMS has extensive forward calorimetry. Steel and quartz fiber hadron forward (HF) calorimeters cover 3 < |η| < 5. Beam Pick-up Timing for the eXperiments (BPTX) devices are used to trigger the detector readout. They are located around the beam pipe at a distance of 175 m from the IP on either side, and are designed to provide precise information on the LHC bunch structure and timing of the incoming beams.

Event and track selections
At the trigger level, the coincidence of signals from both BPTX devices is required, indicating the presence of two bunches of colliding protons and/or ions crossing the IP. In addition, at least one track with p T > 0.4 GeV within |η| < 2.4 is required in the pixel detector. In the offline selection, the presence of at least one tower with energy above 3 GeV in each of the HF calorimeters, and at least one interaction vertex reconstructed in the tracking detectors are required. Beam halo and beam-induced background events, which usually produce an anomalously large number of pixel hits [8], are suppressed. The two analysis methods employ slightly differing additional event and particle selection criteria, which are detailed in Sections 3.1.1 and 3.1.2.

Double ratio method
In the case of the double ratio method, minimum bias events are selected in a similar manner to that described in Refs. [5,6]. Events are required to have at least one reconstructed primary vertex within 15 cm of the nominal IP along the beam axis and within 0.15 cm transverse to the beam trajectory. At least two reconstructed tracks are required to be associated with the primary vertex. Beam-induced background is suppressed by rejecting events containing more than 10 tracks for which it is also found that less than 25% of all the reconstructed tracks in the event pass the highPurity track selection defined in Ref. [9].
A combined event sample produced in pp collisions at √ s = 7 TeV is considered, which uses data from different periods of CMS data taking, i.e. from the commissioning run (23 million events), as well as from later runs (totalling 20 million events). The first data set contains almost exclusively events with a single interaction per bunch crossing, whereas the later ones have a non-negligible fraction of events with multiple pp collision events (pileup). In such cases, the reconstructed vertex with the largest number of associated tracks is selected. An alternative event selection for reducing pileup contamination is also investigated by considering only events with a single reconstructed vertex. This study is then used to assess the related systematic uncertainty. Minimum bias events in pp collisions at 7 TeV (22 million) and at 2.76 TeV (2 million) simulated with the Monte Carlo (MC) generator PYTHIA6.426 [10] with the tune Z2 [11] are also used to construct the single ratios employed in the double ratio technique. Two other tunes (PYTHIA6 D6T [12] and PYTHIA6 Z2* [13,14]) are used for estimating systematic uncertainties related to the choice of the MC tune. The samples for each of these tunes contain 2 million events. The selected tunes describe reasonably well the particle spectra and their multiplicity dependence [8,15,16].
The selected events are then categorized by the multiplicity of reconstructed tracks, N rec1 , which is obtained in the region |η| < 2.4, after imposing additional conditions: p T > 0.4 GeV, |d z /σ(d z )| < 3.0 (impact parameter significance of the track with respect to the primary vertex along the beam axis), |d xy /σ(d xy )| < 3.0 (impact parameter significance in the transverse plane), and σ(p T )/p T < 0.1 (relative p T uncertainty). The resulting multiplicity range is then divided in three (for pp collisions at 2.76 TeV) or four (for pp collisions at 7 TeV) bins. The corrected average charged-particle multiplicity, N trk0. 4 , in the same acceptance range |η| < 2.4 and p T > 0.4 GeV, is determined using the efficiency estimated with PYTHIA, as shown in Table 1. While N rec1 is a measured quantity used to bin the data, for a given bin in this variable, the calculated N trk0.4 value (which is only known on average) facilitates comparisons of the data with models.
In Table 1, two different sets of N rec1 ranges are shown. The upper five (N rec1 = 0-9, 10-14, 15-19, 20-29, 30-79) are used for comparisons with previous CMS results [6]. The larger sample recorded in pp collisions at 7 TeV allowed the analysis to be extended to include the largest bin in multiplicity shown in Table 1.
For the BEC analysis the standard CMS highPurity [9] tracks were used. The additional track selection requirements applied to all the samples mentioned above follow the same criteria as in Refs. [5,6], i.e. primary tracks falling in the kinematic range of |η| < 2.4 with full azimuthal coverage and p T > 200 MeV are used. To remove spurious tracks, primary tracks with a loose selection on the χ 2 of the track fit (χ 2 /ndf < 5, where ndf is the number of degrees of freedom) are used. After fulfilling these requirements, tracks are further constrained to have an impact parameter with respect to the primary vertex of |d xy | < 0.15 cm. Furthermore, the innermost measured point of the track must be within 20 cm of the primary vertex in the plane transverse to the beam direction. This requirement reduces the number of electrons and positrons from photon conversions in the detector material, as well as secondary particles from the decays of long-lived hadrons (K 0 S , Λ, etc.). After applying these requirements, a small residual amount of duplicated tracks remain in the sample. In order to eliminate them, track pairs with an opening angle smaller than 9 mrad are rejected, if the difference of their p T is smaller than 0.04 GeV, a requirement that is found not to bias the BEC results [5,6].

Particle identification and cluster subtraction method
The analysis methods (event selection, reconstruction of charged particles in the silicon tracker, finding interaction vertices, treatment of pileup) used for this method are identical to the ones used in the previous CMS papers on the spectra of identified charged hadrons produced in √ s = 0.9, 2.76, and 7 TeV pp [16], and in √ s NN = 5.02 TeV pPb [17] collisions. In the case of pPb collisions, due to the asymmetric beam energies, the nucleon-nucleon center-of-mass is not at Table 1: Ranges in charged-particle multiplicity over |η| < 2.4, N rec1 , and corresponding average corrected number of tracks with p T > 0.4 GeV, N trk0.4 , in pp collisions at 2.76 and 7 TeV considered in the double ratio method. The N trk0.4 values are rounded off to the nearest integer. A dash indicates that there is not enough data to allow for a good-quality measurement. TeV 7 TeV  0-9  7  6  10-14  14  14  15-19  20  20  20-29  28  28  30-79  47  52   0-9  7  6  10-24  18  19  25-79  42  47  80-110  -105 rest with respect to the laboratory frame, but moves with a velocity |β| = 0.434. Data were taken with both directions of the colliding proton and Pb beams, and are combined in this analysis by reversing the z axis for one of the beam directions.
For this method, 9.0, 9.6, and 6.2 million minimum bias events from pp collisions at √ s = 0.9, 2.76, and 7 TeV, respectively, as well as 9.0 million minimum bias events from pPb collisions at √ s NN = 5.02 TeV are used. The data sets have small fractions of events with pileup. The data samples are complemented with 3.1 million peripheral (60-100% centrality) PbPb events at √ s NN = 2.76 TeV, where 100% corresponds to fully peripheral and 0% to fully central (headon) collisions. The centrality percentages of the total inelastic hadronic cross section for PbPb collisions are determined by measuring the sum of the energies in the HF calorimeters.
This analysis extends charged-particle reconstruction down to p T ≈ 0.1 GeV by exploiting special tracking algorithms used in previous CMS studies [8,[15][16][17] that provide high reconstruction efficiency and low background rates. The acceptance of the tracker is flat within 96-98% in the range −2 < η < 2 and p T > 0.4 GeV. The loss of acceptance for p T < 0.4 GeV is caused by energy loss and multiple scattering of particles, both depending on the particle mass. Likewise, the reconstruction efficiency is about 75-90% for p T > 0.4 GeV, degrading at low p T , also in a mass-dependent way. The misreconstructed-track rate, defined as the fraction of reconstructed primary tracks without a corresponding genuine primary charged particle, is very small, reaching 1% only for p T < 0.2 GeV. The probability of reconstructing multiple tracks from a single true track is about 0.1%, mostly due to particles spiralling in the strong magnetic field of the CMS solenoid. For the range of event multiplicities included in the current study, the efficiencies and background rates do not depend on the multiplicity.
An agglomerative vertex reconstruction algorithm [18] is used, having as input the z coordinates (and their uncertainties) of the tracks at the point of closest approach to the beam axis. The vertex reconstruction resolution in the z direction depends strongly on the number of reconstructed tracks, but is always smaller than 0.1 cm. Only tracks associated with a primary vertex are used in the analysis. If multiple vertices are present, the tracks from the highest multiplicity vertex are used.
The multiplicity of reconstructed tracks, N rec2 , is obtained in the region |η| < 2.4. It should be noted that N rec2 differs from N rec1 . As defined in Section 3.1.1, N rec1 has a threshold of p T > 0.4 GeV applied to the reconstructed tracks, while N rec2 has no such threshold and also includes a correction for the extrapolation down to p T = 0. Over the range 0 < N rec2 < 240, the events are divided into 24 classes, defined in Table 2. This range is a good match to the multiplicity in the 60-100% centrality in PbPb collisions. To facilitate comparisons with models, the corresponding corrected average charged-particle multiplicity N trk in the same region defined by |η| < 2.4 and down to p T = 0 is also determined. For each multiplicity class, the correction from N rec2 to N trk uses the efficiency estimated with MC event generators, followed by the CMS detector response simulation based on GEANT4 [19]. The employed event generators are PYTHIA6.426 (with the tunes D6T and Z2) for pp collisions, and HIJING 2.1 [20,21] for pPb collisions. The corrected data are then integrated over p T , down to zero yield at p T = 0 with a linear extrapolation below p T = 0.1 GeV. The yield in the extrapolated region is about 6% of the total yield. The systematic uncertainty due to the extrapolation is small, well below 1%. Finally, the integrals for each η slice are summed up. In the case of PbPb collisions, events generated with the HYDJET 1.8 [22] MC event generator are simulated and reconstructed. The N rec2 − N trk relationship is parametrized using a fourth order polynomial. Table 2: Relation between the number of reconstructed tracks (N rec2 , p T > 0.1 GeV) and the average number of corrected tracks ( N trk , p T > 0) in the region |η| < 2.4 for the 24 multiplicity classes considered in the particle identification and cluster subtraction method. The values are rounded off to the nearest integer. The corrected N trk values listed for pp collisions are common to all three measured energies. A dash indicates that there is not enough data to allow for a good-quality measurement.

Particle identification
The reconstruction of charged particles in CMS is limited by the acceptance of the tracker and by the decreasing tracking efficiency at low momentum. Particle-by-particle identification using specific ionization losses in the tracker is possible in the momentum range p < 0.15 GeV for electrons, p < 1.15 GeV for pions and kaons, and p < 2.00 GeV for protons (note that protons were not used for correlation studies in this analysis). In view of the (η, p T ) regions where pions, kaons, and protons can all be identified, only particles in the range −1 < η < 1 (in the laboratory frame) are used for this measurement.
For the identification of charged particles, the estimated most probable energy loss rate ε at a reference path length (l 0 = 450 µm) is used [24]. For an accurate determination of ε and its variance σ 2 for each individual track, the responses of all readout chips in the tracker are calibrated with multiplicative gain correction factors. The procedures for gain calibration and track-by-track determination of ε values are the same as in previous CMS analyses [16,17].
Charged particles in the region −1 < η < 1 and 0 < p T < 2 GeV are sorted into bins with  Figure 1: Logarithm of the most probable energy loss rate ε normalized at a reference path length l 0 = 450 µm in pPb collisions at √ s NN = 5.02 TeV. The distribution of ln ε is shown as a function of total momentum p for positively charged particles (left), as well as for identified charged particles (both charges) with high purity selection (> 99.5%, right). The linear z axis scale is shown in arbitrary units. The curves show the expected ln ε for electrons, pions, kaons, and protons (full theoretical calculation, Eq. (33.11) in Ref. [23]). a width of 0.1 units in η, and 0.05 GeV in p T . Since the ratios of particle yields do not change significantly with the charged-particle multiplicity of the event [16,17], the data are not subdivided into bins of N rec2 . The relative abundances of different particle species in a given (η, p T ) bin are extracted by minimizing the joint log-likelihood where µ k are the expected means of ln ε for the different particle species, and P k are their relative probabilities. The value of −2 ln L is minimized by varying P k and µ k starting from reasonable initial values. In the above formula ε i and σ 2 i are the estimated most probable value and its estimated variance, respectively, for each individual track. Since the estimated variance can slightly differ from its true value, a scale factor ξ is introduced. The µ k values are used to determine a unique ln ε(p/m) function. A third-order polynomial closely approximates the collected (p/m, ln ε) pairs, within the corresponding uncertainties. This information is reused by fixing µ k values according to the polynomial, and then re-determining P k .
In this analysis a very high purity particle identification is achieved by requiring P k /(∑ k=π,K,p P k ) > 0.995. If none of the particle type assumptions yields a P k value in this range, the particle is regarded as unidentified. This requirement ensures that less than 1% of the examined particle pairs are misidentified. The degree of purity achieved in particle identification is indicated by distributions of ln ε as a function of total momentum p for pPb collisions at √ s NN = 5.02 TeV in Fig. 1. The left plot shows the initial distribution for positive particles. The plot for negatively charged particles is very similar. This distribution is compared to the theoretically expected energy loss for electrons, pions, kaons, and protons. In Fig. 1, the distribution of ln ε as a function of total momentum for identified charged particles with high purity is shown in the right plot. The plots for pp and PbPb interactions are very similar.

The Bose-Einstein correlations
The correlation functions are investigated in one (1D), two (2D), and three (3D) dimensions in terms of q inv = ( p 1 − p 2 ) 2 − (E 1 − E 2 ) 2 , as well as in projections of the relative threemomentum of the pair, q = p 1 − p 2 , in two (in terms of longitudinal and transverse momenta, q L , q T ), and three (in terms of out-short-long momenta, q L , q S , q O ) directions. In the center-ofmass (CM) frame the variables are defined as: where e k T ≡ k T /k T is a unit vector along the direction of the pair average transverse momentum, k T . In the case of pp collisions, and when dealing with the double ratio technique, the multidimensional investigations are carried out both in the CM frame and in the local co-moving system (LCMS). The latter is the frame in which the longitudinal component of the pair average momentum, k L = (k L,1 + k L,2 )/2, is zero. In the remainder of this paper, a common notation is used to refer to the magnitudes of the average and relative momentum vectors,

Double ratio technique
The analysis procedure using the double ratio technique is the same as that described in Refs. [5,6], where no particle identification is considered. However, the contamination by non-pions is expected to be small, since pions are the dominant type of hadrons in the sample (the ratio of kaons to pions is about 12%, and of protons to pions is roughly 6% [16]). The unidentified kaons and protons would contaminate mainly the low relative momentum region of the correlation functions. The impact of this level of contamination is covered by the systematic uncertainties.
For each event, the signal containing the BEC is formed by pairing same-sign tracks in the same event originating from the primary vertex, with |η| < 2.4 and p T > 0.2 GeV. The distributions in terms of the relative momentum of the pair (the invariant q inv , or the components q T , q L or q O , q S , q L ) are divided into bins of the reconstructed charged-particle multiplicity N trk0.4 , and of the pair average transverse momentum k T . The distributions are determined in the CM and LCMS reference frames.
The background distribution or the reference sample can be constructed in several ways, most commonly formed by mixing tracks from different events (mixed-event technique), which can also be selected in several ways. In the first method employed in this analysis, the reference sample is constructed by pairing particles (all charge combinations) from mixed events and within the same η range, as in Refs. [5,6]. In this case, the full |η| < 2.4 interval is divided in three subranges: −2.4 < η < −0.8, −0.8 ≤ η < 0.8, and 0.8 ≤ η < 2.4. Alternative techniques considered for estimating the systematic uncertainties associated with this choice of reference sample are discussed in Sec. 4.5.1.
The correlation function is initially defined as a single ratio (SR), having the signal pair distribution as numerator and the reference sample as denominator, with the appropriate normalization, where C 2 (q) is two-particle correlation function, N sig is the integral of the signal content, whereas N ref is the equivalent for the reference sample, both obtained by summing up the pair distributions for all the events in the sample. For obtaining the parameters of the BEC effect in this method, a double ratio (DR) is constructed [5,6] as where SR(q) MC is the single ratio as in Eq. (2), but computed with simulated events generated without BEC effects.
Since the reference samples in single ratios for data and simulation are constructed with exactly the same procedure, the bias due to such background construction could be, in principle, reduced when the DR is taken. However, even in this case, the correlation functions are still sensitive to the different choices of reference sample, leading to a spread in the parameters fitted to the double ratios, that is considered as a systematic uncertainty [5,6]. On the other hand, by selecting MC tunes that closely describe the behavior seen in the data, this DR technique should remove the contamination due to non-Bose-Einstein correlations. Table 3 shows the ranges and bin widths of relative momentum components used in fits to the double ratios. The bins in k T are also listed, which coincide with those in Refs. [5,6].

Particle identification and cluster subtraction technique
Similarly to the case of unidentified particles discussed in Section 4.1, the identified hadron pair distributions are binned in the number of reconstructed charged particles N rec2 of the event, in the pair average transverse momentum k T , and also in the relative momentum variables in the LCMS of the pair in terms of q inv , (q L , q T ) or (q L , q O , q S ). The chosen bin widths for pions and kaons are shown in Table 4. Table 4: Chosen bin widths for the various variables studied for pions and kaons in the particle identification and cluster subtraction method. The construction of the signal distribution is analogous to that described above: all valid particle pairs from the same event are taken and the corresponding distributions are obtained. Several types of pair combinations are used to study the BEC: π + π + , π − π − , K + K + , and K − K − ; while π + π − and K + K − are employed to correct for non-BEC contributions, as described in Section 4.3.2. For the reference sample an event mixing procedure is adopted, in which particles from the same event are paired with particles from 25 preceding events. Only events belonging to the same multiplicity (N rec2 ) class are mixed. Two additional cases are also investigated. In one of them particles from the same event are paired, but the laboratory momentum vector of the second particle is rotated around the beam axis by 90 degrees. Another case considers pairs of particles from the same event, but with an opposite laboratory momentum vector of the second particle (mirrored). Based on the goodness-of-fit distributions of correlation function fits, the first event mixing prescription is used, while the rotated and mirrored versions, which give considerably worse χ 2 /ndf values, are employed in the estimation of the systematic uncertainty.
The measured two-particle correlation function C 2 (q) is itself the single ratio of signal and background distributions, as written in Eq. (2). This single ratio contains the BEC from quantum statistics, while non-BEC effects also contribute. Such undesired contributions are taken into account as explained in Section 4.3. A joint functional form combining all effects is fitted in order to obtain the radius parameter and the correlation function intercept, as discussed in Section 4.4.

Effect of Coulomb interaction and correction
The BEC method employed for the study of the two-particle correlations reflects not only the quantum statistics of the pair of identical particles, but is also sensitive to final-state interactions. Indeed, the charged-hadron correlation function may be distorted by strong, as well as by Coulomb effects. For pions, the strong interactions can usually be neglected in femtoscopic measurements [25].
The Coulomb interaction affects the two-particle relative momentum distribution in a different way in the case of pairs with same-sign (SS) and opposite-sign (OS) electric charge. This leads to a depletion (enhancement) in the correlation function caused by repulsion (attraction) of SS (OS) charges. The effect of the mutual Coulomb interaction is taken into account by the factor K, the squared average of the relative wave function Ψ, as K(q inv ) = d 3 r f ( r) |Ψ( k, r)| 2 , where f ( r) is the source intensity of emitted particles and k is the relative momentum of the pair. Pointlike sources, f ( r) = δ( r), result in an expression for K(q inv ) coincident with the Gamow factor [26] which, in case of same sign and opposite sign charges, is given by where ζ = αm/q inv , is the Landau parameter, α is the fine-structure constant, m is the particle mass, and q inv is the invariant relative momentum of the pair [27].
For extended sources, a more elaborate treatment is needed [28,29], and the Bowler-Sinyukov formula [25,30] is most commonly used. Although this does not differ significantly from the Gamow factor correction in the case of pions, this full estimate is required for kaons. The possible screening of Coulomb interaction, sometimes seen in data of heavy ion collisions [31,32], is not considered. The value of zeta is typically ζ 1 in the q range studied in this analysis. The absolute square of the confluent hypergeometric function of the first kind F, present in CMS pp, √ s = 0.9 TeV Figure 2: Contribution of resonance decays to the measured Coulomb-corrected correlation function of π + π − (for 2 ≤ N rec2 < 10, left, blue squares) and K + K − (for 2 ≤ N rec2 < 60, right, blue squares). The φ(1020) peak (right panel) is off-scale.
Ψ, can be approximated as |F| 2 ≈ 1 + 2ζ Si(x), where Si is the sine integral function [33]. Furthermore, for Cauchy-type source functions ( f (r) = R/{2π 2 r 2 + (R/2) 2 2 }, where f is normalized, such that f (r) d 3 r = 1), the Coulomb-correction factor K is well described by the formula for q inv in GeV, and R in fm. The factor π in the approximation comes from the fact that, for large kr arguments, Si(kr) → π/2. Otherwise it is a simple but accurate approximation of the result of a numerical calculation, with deviations below 0.5%.

Clusters: mini-jets, multi-body decays of resonances
The measured opposite-sign correlation functions show contributions from various hadronic resonances [34]. Selected Coulomb-corrected correlation functions for low N rec2 bins are shown in Fig. 2. The observed resonances include K 0 S , ρ(770), f 0 (980), and f 2 (1270) decaying to π + π − , and φ(1020) decaying to K + K − . Photon conversions into e + e − pairs, when misidentified as pion pairs, can also appear as a peak at very low q inv in the π + π − spectrum. With increasing N rec2 values the effect of resonances diminishes, since their contribution is quickly exceeded by the combinatorics of unrelated particles.
The Coulomb-corrected OS correlation functions are not always close to unity even at low q inv , but show a Gaussian-like hump (Fig. 3). That structure has a varying amplitude with a stable scale (σ of the corresponding Gaussian) of about 0.4 GeV. This feature is often related to particles emitted inside low-momentum mini-jets, but can be also attributed to the effect of multi-body decays of resonances. In the following we will refer to such an effect as due to fragmentation of clusters, or cluster contribution [35][36][37]. For the purpose of evaluating and later eliminating that contribution, the one-dimensional OS correlation functions are fitted with the following form: The amplitude b is inversely proportional to a power of N rec2 (with the exponent n b in the range 0.80-0.96). The parameter b 0 is much bigger for 5 TeV pPb than for 2.76 TeV PbPb data; for pp collisions the dependence is described by 0.28 ln( √ s/0.48 TeV). For pPb collisions b 0 is about 1.6-1.7 times higher than would be predicted from the pp curve at the same √ s NN . At the same time b increases exponentially with k T , the parameter k 0 being in the range of 1.5-2.5 GeV. The Gaussian width σ b of the hump at low q inv first decreases with increasing N rec2 , but for N rec2 > 70 remains constant at about 0.35-0.55 GeV. The width increases with k T as a power law, with the exponent n T in the range 0.19-0.33. The expression bσ 2 b N rec2 is proportional to the fraction of particles that have a cluster-related partner. Our data show that this fraction does not change substantially with N rec2 , but increases with k T and √ s for pp collisions. In summary, the object that generates this structure is assumed to always emit particles in a similar way, with a relative abundance that increases with the center-of-mass energy.
The cluster contribution can also be extracted in the case of the SS correlation function, if the momentum scale of the BEC effects and that of the cluster contribution (≈0.4 GeV) are different enough. An important constraint in both mini-jet and multi-body resonance decays is the conservation of electric charge that results in a stronger correlation for OS pairs than for SS pairs. Consequently, the cluster contribution is expected to be present as well for SS pairs, with similar shape but a somewhat smaller amplitude. The form of the cluster-related contribution obtained from OS pairs, multiplied by a relative amplitude z(N rec2 , k T ) (ratio of measured OS to SS contributions), is used to fit the SS correlations. The dependence of this relative amplitude on N rec2 in k T bins is shown in Fig. 4. Results (points with error bars) from all collision systems are plotted together with a parametrization ((aN rec2 + b)/(1 + aN rec2 + b)) based on the combinatorics of SS and OS particle pairs (solid curves). Instead of the points, these fitted curves are used to describe the N rec2 dependence of z in the following. The ±0.2 grey bands are chosen to cover most of the measurements and are used in determining a systematic uncertainty due to   the subtraction of the cluster contribution.
Thus, the fit function for SS pairs is where K ±± describes the Coulomb effect, the squared brackets encompass the cluster contribution, and C 2,BE (q inv ) contains the quantum correlation to be extracted. Only the normalization c SS and the parameters of C 2,BE (q inv ) are fitted, all the other variables (those in K ±± , z(N rec2 ), b, and σ b ) were fixed using Eq. (7).
In the case of two and three dimensions, the measured OS correlation functions are constructed as a function of a length given by a weighted sum of q components, instead of q inv . The Coulomb-corrected OS correlation functions of pions as a function of a combination of q L and q T , q h = q 2 L + (aq T ) 2 , in selected N rec2 bins for all k T , are shown in Fig. 5. The dependence of the a factor on N rec2 and k T is described by: This particular form is motivated by the physics results pointing to an elongated source, shown later in Section 6. The small peak at q h ≈ 0 is from photon conversions where both the electron and the positron are misidentified as a pion. The very low q h region is not included in the fits. Instead of the measurements, the fitted curves shown in Fig. 5 are used to describe the N rec2 and k T dependence of these parameters in the following. Hence the fit function for SS pairs in the multidimensional case is:  Figure 5: Coulomb-corrected OS correlation function of pions (blue squares) as a function of a combination of q L and q T , in two selected N rec2 bins for all k T values (left: 40 ≤ N rec2 < 50, right: 160 ≤ N rec2 < 170). For better visibility, only a randomly selected small fraction of the points is plotted. In addition, those with statistical uncertainty higher than 10% are plotted with light grey color. The solid curves indicate the prediction of the parametrized cluster contribution.
In the formula above, only the normalization c and the parameters of C 2,BE ( q) are fitted, all the other variables (those in K ±± , z(N rec2 ), b, σ b , and a) are fixed using the parametrization described above.
There are about two orders of magnitude less kaon pairs than pion pairs, and therefore a detailed study of the cluster contribution is not possible. We assume that their cluster contribution is similar to that of pions, and that the parametrization deduced for pions with identical parameters can be used (for b and σ b in Eq. (7), and for a in Eq. (9)). This choice is partly justified by OS kaon correlation functions shown later, but based on the observed difference, the systematic uncertainty on z(N rec2 ) for kaons is increased to 0.3.

Fitting the correlation function
The parametrizations used to fit the correlation functions in 1D, 2D, and 3D in terms of q inv , (q L , q T ), and (q S , q L , q O ), respectively, are the following [5,6,38]: is the fit function employed in the 1D case, while in the 2D case it has the form: and, finally, the form used for 3D is: The radius parameters R inv , (R T , R L ), (R S , R L , R O ), correspond to the lengths of homogeneity in 1D, 2D and 3D, respectively. In the above expressions, λ is the intercept parameter (intensity of the correlation at q = 0) and C is the overall normalization. The polynomial factors on the right side, proportional to the fit variables δ inv , δ T , δ L , δ S , δ O , are introduced for accommodating possible deviations of the baseline from unity at large values of these variables, corresponding to long-range correlations. They are only used in the double ratio method, their values are uniformly set to zero in the case of the particle identification and cluster subtraction method that contains no such factors.
In Eq. (12), T,L are the geometrical transverse and longitudinal radii, respectively, whereas β T = k T /k 0 and β L = k L /k 0 are the transverse and longitudinal velocities, respectively; φ is the angle between the directions of q T and k T , and τ is the source lifetime. Similarly, in Eq.
T are the corresponding radius parameters in the 3D case. When the analysis is performed in the LCMS, the cross-terms q T q L and q O q L do not contribute for sources symmetric along the longitudinal direction.
Theoretical studies show that for the class of Lévy stable distributions (with index of stability 0 < a ≤ 2) [38] the BEC function has a stretched exponential shape, corresponding to a = 1 in Eqs. (12) and (13). This is to be distinguished from the simple exponential parametrization, in which the coefficients of the exponential would be given by ∑ i R i q i . In Section 5.3, the parameters resulting from the fits obtained with the simple and the stretched exponentials will be discussed in the 3D case. If a is treated as a free parameter, it is usually between a = 1 and a = 2 (which corresponds to the Gaussian distribution [38]). In particular, for a = 1, the exponential term in expressions given by Eqs. (12) and (13) coincides with the Fourier transform of the source function, characterized by a Cauchy distribution.
The measurements follow Poisson statistics, hence the correlation functions are fitted by minimizing the corresponding binned negative log-likelihood.

Double ratio method
Various sources of systematic uncertainties are investigated with respect to the double ratio results, similar to those discussed in Refs. [5,6]: the MC tune used, the reference sample employed to estimate the single ratios, the Coulomb correction function, and the selection requirements applied to the tracks used in the analysis. From these, the major sources of systematic uncertainty come from the choice of the MC tune, followed by the type of reference sample adopted. In the latter case, alternative techniques to those considered in this analysis for constructing the reference sample are studied: mixing tracks from different events selected at random, taking tracks from different events with similar invariant mass to that in the original event (within 20%, calculated using all reconstructed tracks), as well as combining oppositely charged particles from the same event. Each of these procedures produces different shapes for the single ratios, and also leads to significant differences in the fit parameters. The exponential function in Eq. (11) (Lévy type with a = 1) is adopted as the fit function in the systematic studies.
The results of the analysis do not show any significant dependence on the hadron charges, as obtained in a study separating positive and negative charges in the single ratios. Similarly, the effect of pileup, investigated by comparing the default results to those where only single-vertex events are considered, does not introduce an extra source of systematic uncertainty.
The sources of systematic uncertainties discussed above are considered to be common to pp collisions at both 2.76 and 7 TeV. In the 1-D case the overall systematic uncertainty associated with R inv and that associated with λ are estimated using the change in the fitted values found when performing the variations in procedure mentioned above for the four major sources of uncertainties at both energies. This leads to 15.5% and 10% for the systematic uncertainties in the fit radius and intercept parameters, respectively, which are considered to be common to all N trk and k T bins. The estimated values in the multidimensional cases were of similar order or smaller than those found in the 1-D case. Therefore, the systematic uncertainties found in the 1-D case were also applied to the 2-D and 3-D results.

Particle identification and cluster subtraction method
The systematic uncertainties are dominated by two sources: the construction of the background distribution, and the amplitude z of the cluster contribution for SS pairs with respect to that for OS ones. Several choices for the construction of the background distribution are studied. The goodness-of-fit distributions clearly favor the event mixing prescription, while the rotated version gives poorer results, and the mirrored background yields the worst performance. The associated systematic uncertainty is calculated by performing the full analysis using the mixed event and the rotated variant, and calculating the root-mean-square of the final results (radii and λ).
Although the dependence of the ratio of the mini-jet contribution z (between SS and OS pairs), as a function of N rec2 could be well described by theory-motivated fits in all k T regions, the extracted ratios show some deviations within a given colliding system, and also between systems. In order to cover those variations, the analysis is performed by moving the fitted ratio up and down by 0.2 for pions, and by 0.3 for kaons. The associated systematic uncertainties are calculated by taking the average deviation of the final results (radii and λ) from their central values.
In order to suppress the effect of multiply-reconstructed particles and misidentified photon conversions, the low-q region (q < 0.02 GeV) is excluded from the fits. To study the sensitivity of the results to the size of the excluded range, its upper limit is doubled and tripled. Both the radii and λ decreased slightly with increasing exclusion region. As a result, the contribution of this effect to the systematic uncertainty is estimated to be 0.1 fm for the radii and 0.05 for λ.
The systematic uncertainties from all the sources related to the MC tune, Coulomb corrections, track selection, and the differences between positively and negatively charged particles are similar to those found in case of the double ratio method above. The systematic uncertainties from various sources are added in quadrature.

Comparisons of the techniques
In order to illustrate the level of consistency between both BEC analyses, the fit results for 1D R inv and λ obtained employing the double ratio and the particle identification and cluster subtraction methods are compared. In the double ratio method, mixed events having similar multiplicities within the same η ranges are used as the reference sample. The non-BEC effects are corrected with the double ratio technique, as discussed in Section 4.1, which also reduces the bias due to the choice of reference sample. In the particle identification and cluster subtraction method, single ratios are measured considering mixed pairs from 25 events chosen at random, whose contamination from resonance decays and mini-jet contributions are corrected as discussed in  11), obtained in the double ratio (filled markers) and the particle identification (for pions, empty squares) methods in pp collisions at 7 TeV, together with the results from Ref. [6]. The inner (outer) error bars correspond to statistical (systematic) uncertainties.
parameters are obtained. Since in the particle identification and cluster subtraction method N trk is fully corrected, as discussed in Section 3.1.2 and summarized in Table 2, an additional correction is applied to the values of N trk0.4 in Table 1 in the double ratio method, thus allowing the comparison of both sets of results plotted at the corresponding value of N trk . Such a correction is applied with a multiplicative factor obtained from the ratio of the total charged hadron multiplicities in both analyses, N ch (p T → 0)/N ch (p T > 0.4 GeV) ∼ 1.7. The measurements are plotted in Fig. 6 for charged hadrons in the double ratio method and from identified charged pions from the particle identification and cluster subtraction method, as well as those for charged hadrons from Ref. [6]. In the latter, the correction factor is obtained in a similar fashion as N ch (p T → 0)/N ch (p T > 0.2 GeV) ∼ 1.1, reflecting the p T > 0.2 GeV requirement applied in that case.
The radius parameter R inv (Fig. 6, left) steadily increases with the track multiplicity. The results using the two different correlation techniques agree within the uncertainties, and they also agree well with our previous results [6]. In Fig. 6 (right), the corresponding results for the intercept parameter, λ, are shown integrated over all k T . The overall behavior of the correlation strength is again similar in both cases, showing first a decrease with N trk and then a flattening, approaching a constant for large values of the track multiplicity. The differences in the absolute values of the fit parameter λ in the two approaches are related to the differences in p T (and k T ) acceptances for the tracking methods used in the two analyses, and the availability of particle identification, with both factors lowering the fit value of λ for the double ratio method.

One-dimensional results
Single and double ratios for pp collisions at √ s = 2.76 TeV are shown in Fig. 7, integrated over the charged-particle multiplicity, N trk0.4 , and the pair momentum, k T . The fit to the double   ratio is performed using the exponential function in Eq. (11), with a = 1. The corresponding fit values of the invariant radii, R inv , and the intercept parameter, λ, are shown in Table 5, together with those for pp collisions at √ s = 7 TeV; these being compatible within the uncertainties at these two energies. The upper panel of Fig. 8 shows the results for the invariant radius R inv (left) and the intercept parameter λ (right), obtained with the exponential fit function in three bins of both N trk0.4 and k T , for pp collisions at √ s = 2.76 and 7 TeV. The results corresponding to the two energies again show good agreement in the various (N trk0.4 , k T ) bins considered, extending the observation made with respect to the integrated values in Table 5 to different bin combinations of these variables. The parameter λ decreases with increasing k T , and its dependence on the charged multiplicity N trk0.4 follows a similar trend. The fit invariant radius, R inv , decreases with k T , a behavior previously observed in several different collision systems and energy ranges [5,6,34,[39][40][41][42], and expected for expanding emitting sources. In the case of heavy ion collisions, as well as for deuteron-nucleus and proton-nucleus collisions, both at the BNL RHIC and at the CERN LHC, this observed trend of the fit radii is compatible with a collective behavior of the source. The radius fit parameter is also investigated in Fig. 8 (bottom) for results integrated over k T bins and in five intervals of the charged-particle multiplicity, as in Table 1. The value of R inv steadily increases with multiplicity for pp collisions at 2.76 and 7 TeV, as previously observed at 7 TeV in Ref. [6]. The various curves are fits to the data, and are proportional to (N trk0.4 ) 1/3 , showing also that the results are approximately independent on the center-of-mass energy.
As discussed in Ref. [6], the exponential fit alone in Eq. (11) is not able to describe the overall behavior of the correlation function in the range 0.2 < q inv < 2 GeV, for which an additional long-range term proportional to q inv is required. Although the general trend can be described by this parametrization, it misses the turnover of the baseline at large values of q inv . The resulting poor fit quality is partially due to very small statistical uncertainties which restrict the variation of the parameters, but mainly due to the presence of an anticorrelation, i.e. correlation function values below the expected asymptotic minimum at unity, which is observed in the intermediate q inv range (0.4 < q inv < 1.1 GeV). It is unlikely that this anticorrelation is due to the contamination of non-pion pairs in the sample, which would mainly give some contribution in the lower q inv region.
The anticorrelation depth depends on the range of N trk0.4 considered, as reported in Ref. [6]. This is investigated further in this analysis by considering different k T bins for pp collisions at 2.76 and 7 TeV. The results for the lower energy are shown in Fig. 9, where the data and the exponential fit are shown together with those based on Eq. (14). As previously observed in e + e − collisions [43] and also in Ref. [6], this anticorrelation feature is compatible with effects as suggested by the τ model [44], where particle production has a broad distribution in proper time and the phase space distribution of the emitted particles is dominated by strong correlations of the space-time coordinate and momentum components. In this model the time evolution of the source is parameterized by means of a one-sided asymmetric Lévy distribution, leading to the fit function: where the parameter r 0 is related to the proper time of the onset of particle emission, r ν is a scale parameter entering in both the exponential and the oscillating factors, and ν corresponds to the Lévy index of stability.
The appearance of such anticorrelations may be due to the fact that hadrons are composite, extended objects. The authors of Ref. [45] propose a simple model in which the final-state hadrons are not allowed to interpenetrate each other, giving rise to correlated hadronic positions. This assumption results in a distortion in the correlation function, leading to values below unity.
Fitting the double ratios with Eq. (14) in the 1D case considerably improves the values of the χ 2 /ndf as compared to those fitted with the exponential function in Eq. (11), for pp collisions at both 2.76 and 7 TeV [6]. From Fig. 9, it can be seen that the fit based on Eq. (14) describes the overall behavior of the measurements more closely in all the bins of N trk0.4 and k T investigated. In all but one of those bins, the fit based on the τ model results in 1 < χ 2 /ndf < 2. In addition, a significant improvement is also observed when 2D and 3D global fits to the correlation functions are employed, mainly when the Lévy fit with a = 1 is adopted (Sections 5.2 and 5.3).
The polynomial form C(1 + δ q) in Eq. (14) is used as a baseline to quantify the depth of the dip [6]. This is estimated by the difference, ∆, of the baseline function and the fit based on Eq. (14) at the minimum, leading to the results shown in Fig. 10. The left plot shows the present results for pp collisions at 2.76 and at 7 TeV as a function of the multiplicity N trk0.4 integrated over k T bins, exhibiting their consistency at different energies, as well as close similarity with the previous results [6]. The right plot in Fig. 10 shows the current results corresponding to the 2.76 and 7 TeV data further scrutinized in different N trk0.4 and k T bins. From this plot, it can be seen that the depth of the dip decreases with increasing N trk0.4 for both energies; however, the dependence on k T seems to decrease and then flatten out. In any case, the dip structure observed in the correlation function is a small effect, ranging from ≈3% in the lower N trk0. 4 and k T range investigated, to ≈1% in the largest one.
In summary, the BEC results obtained in pp collisions and discussed in this section are shown to have a similar behavior to those observed in such systems in previous measurements at the LHC, by CMS [5,6], as well as by the ALICE [46] and ATLAS [47] experiments. This behavior is reflected in the decreasing radii for increasing pair average transverse momentum k T , and is also observed in e + e − collisions at LEP [48]. Such behavior is compatible with nonstatic, expanding emitting sources. Correlations between the particle production points and their momenta, (x, p) correlations, appear as a consequence of this expansion, generating a dependence of the correlation function on the average pair momentum, and no longer only on their relative momentum; (x, p) correlations were modeled in Refs. [44,45] as well. Similar behavior is also observed in various collision systems and energy ranges [5,6,34,[39][40][41][42]. In the case of heavy ion, deuteron-nucleus, and proton-nucleus collisions, both at RHIC and at the LHC, these observations are compatible with a collective behavior of the source.

Two-dimensional results
The 1D analysis is extended to the 2D case in terms of the components (q T , q L ) of the pair relative momentum, for data samples at √ s = 2.76 and 7 TeV. In the 2D and 3D cases the measurement is conducted both in the CM frame and in the LCMS. For clarity, the asterisk symbol (*) is added to the variables defined in the CM frame when showing the results in this section. In the case of longitudinally symmetric systems, the LCMS has the advantage (Section 4.4) of omitting the cross-term proportional to q T q L in the 2D, or q L q O in the 3D fit functions.
The 2D correlation function provides information on the behavior of the emitting source along and transverse to the beam direction. Typical examples are illustrated in Fig. 11 for pp collisions at √ s = 2.76 TeV, in the LCMS. The upper panel shows the 2D plot of the double ratio in terms of q T and q L , with the right plot showing the 2D Lévy fit (with a = 1) defined in Eq. (12), superimposed on the 2D correlation function. The lower plots show 1D projections in terms of q L (left) and q T (right) of the 2D double ratios, and the corresponding 1D projections of the 2D Lévy (with a = 1) fit function. When plotting in terms of q L , only the first bin in q T (for q T < 0.05 GeV) is considered, and vice-versa. This choice is made to exhibit the largest values of the 2D correlation function in each direction, since the BEC decreases with increasing q T , q L . The functions shown in the 1D figures are not fits to that data, but are projections of the global 2D-stretched exponential fit to the 2D correlation function. This explains the poor description of the data by the Lévy fit with a = 1, in Fig. 11, left. The statistical uncertainties are also larger in the q L case, as well as the fluctuations in the results along this direction.
Analogously to the studies performed in 1D, the correlation functions are also investigated in 2D in terms of the components (q decrease with k T for larger multiplicities, a behavior similar to that observed in the 1D case, and expected for expanding sources. The observation that R L (LCMS) > R * L (CM) for the same bins of N trk0.4 and k T can be related to the Lorentz-contraction of the longitudinal radius in the CM frame. Figure 13 shows the results for the correlation strength parameter, λ ( * ) , obtained with the same fit procedure, for both the CM and LCMS frames. No significant sensitivity of the intercept is seen as a function of k T . However, within each k T range, λ ( * ) slowly decreases with increasing track multiplicity, in a similar way in both frames. The values obtained with the stretched exponential fit in the LCMS for the radius parameters R T , R L , together with the intercept λ, corresponding to data integrated both in N trk0.4 and k T , are collected in Table 6 for pp collisions at 2.76 and 7 TeV. The values fitted to the cross-term R LT in the CM frame are of order 0.5 fm or less, and in the LCMS such terms are negligible (O[10 −7 ]), as expected. From Table 6 it can be seen that, in the LCMS and at both energies, R L ≈ 4R T /3, suggesting that the source is longitudinally elongated.
In Fig. 14 the 2D results for the double ratios versus (q L , q T ) in the LCMS are shown zoomed along the correlation function axis, with the BEC peak cut off around 1.2, for four chargedmultiplicity bins, N trk0.4 . These plots illustrate the directional dependence of the k T -integrated correlation functions along the beam and transverse to its direction, as well as the variation of their shapes when increasing the multiplicity range (from the upper left to the lower right plot) of the events considered. The edge effects seen for q T > 1.5 GeV in the upper left plot (for N trk0.4 = 6) are an artifact of the cutoff at 1.2, being much smaller than the true value of the peak in the low (q L , q T ) region, and correspond to fluctuations in the number of pairs in those bins, due to a lower event count than in the other N trk0.4 intervals. The N trk0.4 ranges are the same as in Figs. 12-13, although the fit parameters shown in those figures are obtained considering differential bins in both N trk0.4 and k T . Figure 14 also shows the 1D projections of the 2D double ratios in terms of q T and q L (for values q L,T < 0.05 GeV). The shallow anticorrelation seen in the 1D double ratios in terms of q inv seems to be also present in the 2D plots, suggesting slightly different depths and shapes along the q L and q T directions.    [6] are also shown, with their corresponding statistical (inner error bars) and statistical plus systematic uncertainties added in quadrature (outer error bars). The curves indicate the respective fits to a N 1/3 trk0.4 functional form.    Top: 2D correlation functions as a function of (q T , q L ) without (left) and with (right) the 2D Lévy (with a = 1) fit function superimposed. Bottom: Corresponding 1D projections in terms of q L (for q T < 0.05 GeV, left) and q T (for q L < 0.05 GeV, right). The error bars indicate statistical uncertainties.      Figure 14: Results obtained in the LCMS for the 2D double ratios in the (q L , q T ) plane, zoomed along the correlation function axis, for four charged multiplicity bins, N trk0.4 (integrated over all k T ), increasing from upper left to lower right. The BEC peak is cut off at ∼1.2 for a better visualization of the directional behavior. The 1D projections in q L (for q T < 0.05 GeV) and q T (for q L < 0.05 GeV) are shown to the right of the bidimensional plots, with error bars corresponding to statistical uncertainties.

Three-dimensional results
The double ratios can be further investigated in 3D as a function of the relative momentum components (q O , q S , q L ). In order to obtain the corresponding fit parameters, (R O , R S , R L ), as well as the correlation function strength, λ, the fits are performed to the full 3D double ratios and, if desired for illustration purposes, projected onto the directions of the variables q O , q S , q L , similarly to the projections done for the measurements. The 3D correlation functions can be visualized through 2D projections in terms of the combinations (q S , q L ), (q L , q O ), and (q O , q S ), with the complementary components constrained to q O < 0.05 GeV, q S < 0.05 GeV, and q L < 0.05 GeV, respectively, corresponding to the first bin in each of these variables. The lower panel in Fig. 15 shows the corresponding 1D projections of the full 3D correlation function and fits. Similarly to what was stressed regarding Fig. 11, the functions shown in Fig. 15 are projections of the global 3D-stretched exponential fit to the 3D correlation function. This explains the poor description of the 1D data by the Lévy fit with a = 1.  The values of the radii in the Bertsch-Pratt parametrization [49,50] for pp collisions at both 2.76 and 7 TeV, obtained with the stretched exponential fit in the LCMS and integrating over all N trk0.4 and k T ranges, are summarized in Table 7, together with the corresponding intercept fit parameter. Comparing the LCMS radii in the 3D case for pp collisions at 7 TeV, the hierarchy R L > R O and R L R S can be seen (although the values of R L and R S are also compatible with each other within their statistical and systematic uncertainties). Therefore, as in the 2D case, the 3D source also seems to be more elongated along the longitudinal direction in the LCMS, with the relationship between the radii being given by R L R S > R O . A similar inequality holds for pp collisions at 2.76 TeV.
The data sample obtained in pp collisions at 2.76 TeV in 2013 is considerably smaller than at 7 TeV and would lead to large statistical uncertainties if the double ratios in the 3D case are measured differentially in N trk0.4 and k T , so this study is not considered at 2.76 TeV.
For the data from pp collisions at 7 TeV, the fits to the double ratios are investigated in three bins of the pair average transverse momentum, k T (integrating over charged multiplicity, N trk0.4 ). The fit parameters are obtained with exponential and Lévy (with a = 1) type functions. The corresponding results are compiled in Fig. 16, both in the CM frame and in the LCMS, showing that the values are very sensitive to the expression adopted for the fit function, as also suggested in Fig. 15 by the 1D projections of the global 3D fit performed to the 3D double ratios.
In the LCMS, the behavior of R S with increasing pair average transverse momentum is very similar to that of R * S in the CM frame, and the fit values are consistent within statistical and systematic uncertainties, as would be expected, since the boost to the LCMS is in the longitudinal direction. It can also be seen in both frames that the sideways radius is practically independent of k T , within the uncertainties. The outwards fit parameter R ( * ) O shows a similar decrease with respect to k T in both frames. However, in the CM frame the R * O fit values are slightly larger than those in the LCMS (R O ) for the same (N trk0.4 , k T ) bins. This is compatible with the dependence of R ( * ) O on the source lifetime (as seen in Section 4.4), which is expected to be higher in the CM frame than in the LCMS, because of the Lorentz time dilation. In the LCMS, R L decreases with increasing k T . It has larger values than in the CM frame for both energies, suggesting an effect related to the Lorentz boost, similar to the 2D case. Also in this frame, its decrease with increasing k T is more pronounced for the same reason.
The fits to the double ratios are also studied in four bins of N trk0.4 integrating over all k T , and the corresponding results are shown in Fig. 17, both in the CM frame (upper) and in the LCMS (lower). A clear trend can be seen, common to both fit functions and in all directions of the relative momentum components: the radii R The intercept parameter λ ( * ) is also studied as a function of k T and N trk , both in the CM frame and in the LCMS. The corresponding results are shown in Fig. 18. A moderate decrease with increasing k T is observed. As a function of N trk0.4 , λ ( * ) first decreases with increasing charged- particle multiplicity and then it seems to saturate.
Since the longitudinal radius parameter represents the same length of homogeneity along the beam direction, it would be expected that the corresponding fit values in the 2D and 3D cases are comparable. In fact, they are consistent within the experimental uncertainties, as can be seen in Fig. 19, except in the first bin of N trk0.4 investigated in the CM frame, where the central R L value in 3D is roughly 0.3 fm larger than the corresponding R L value in 2D, for pp collisions at both 2.76 and 7 TeV. This difference, however, may reflect a better description of the source in the 3D case, where the transverse component is decomposed into two independent orthogonal directions, whereas it may not be fully accommodated by considering just one transverse parameter as in the 2D case. In addition, it can be seen in Fig. 19 that there is an approximate energy independence of R * L (R L ) as a function of N trk0. 4 , for results at √ s = 2.76 and 7 TeV, similar to what is observed in the 1D case.  Figure 17: 3D radii parameters as a function of N trk0.4 (integrated over all k T ), obtained from double ratios fitted to the exponential (left) and Lévy (with a = 1, right) functions for pp collisions at √ s = 7 TeV in the CM frame (top) and in the LCMS (bottom). The boxes indicate the systematic uncertainties.

Results with particle identification and cluster subtraction in pp, pPb, and PbPb collisions
A selection of correlation functions obtained using the mixed-event prescription, as described in Section 4.2, and the corresponding fits are shown in Figs. 20 (pions) and 21 (kaons). The SS correlation functions, corrected for Coulomb interaction and cluster contribution (mini-jets and multi-body resonance decays), as a function of q inv , (q L , q T ), and (q L , q O , q S ) in selected N rec2 bins are plotted for all k T . The solid curves indicate fits with the stretched exponential parametrization. Although the probability of reconstructing multiple tracks from a single true track is very small (about 0.1%), the bin with the smallest | q| is not used in the fits in order to avoid regions potentially containing pairs of multiply-reconstructed tracks.
The fits are mostly of good quality, the corresponding χ 2 /ndf values being close to one in the 2D and 3D cases, for pions and kaons alike. The stretched exponential parametrization provides an excellent match to the data. The 1D correlation function shows a wide minimum around 0.5 GeV (Fig. 20, upper panel). The amplitude of that very shallow dip is of the order of a few percent only, and is compatible with the previously discussed results.
The characteristics of the extracted one-, two-, and three-dimensional correlation functions of identified pions and kaons in pp, pPb, and peripheral PbPb collisions are presented in Figs. 22-30 as functions of the pair average transverse momentum k T and of the corrected chargedparticle multiplicity N trk of the event (in the range |η| < 2.4 in the laboratory frame, with integrated results extrapolated to p T = 0). In all plots, the results of positively and negatively charged pions and kaons are averaged. The central values for the radii and intercept parameter λ are given by markers. The statistical uncertainties are indicated by vertical error bars, while the combined systematic uncertainties (choice of background method; uncertainty in the rela- tive amplitude z of the cluster contribution; low q exclusion) are given by the filled boxes. Unless indicated, the lines are drawn to guide the eye (cubic splines whose coefficients are found by weighting the measurements with the inverse of their squared statistical uncertainty).

Dependencies on N trk and k T
The N trk dependence of the 1D radius R inv and intercept λ parameters for pions (Eq. (8)), averaged over several k T bins, for all studied systems are shown in Fig. 22. The radii found for pions from pp collisions span a similar range to those found previously [5,6], as well as those obtained for charged hadrons via the double ratio technique, as reported in Section 5.1. Although illustrated in Fig. 22 for the k T -integrated case only, the dependence on the number of tracks of the 1D pion radius parameter is remarkably similar for pp and pPb collisions in all k T bins. This also applies to PbPb collisions when k T > 0.4 GeV. The values of the pion intercept parameter λ are usually below unity, in the range 0.7-1, and seem to remain approximately constant as a function of N trk for all systems and k T bins. Similar results obtained for kaons in the 1D analysis are displayed in Fig. 23. The 1D kaon radius parameter increases with N trk , but with a smaller slope than seen for pions.
The k T dependencies of the 1D radius parameters for pions, in several N trk bins for the pp and pPb systems are plotted in Fig. 24. They show a decrease with increasing k T for k T > 0.4 GeV.
The N trk dependence of the 2D radii and the corresponding intercept parameters for pions (Eq. (10)), in several k T bins, for all studied systems are shown in Fig. 25. The N trk dependence of 2D pion radius parameters, R L and R T , is similar for pp and pPb collisions in all k T bins. In general R L > R T , indicating that the source is elongated in the beam direction, a behavior also seen in pp collisions with the double ratio technique, as discussed in Section 5.1. In the case of PbPb collisions R L ≈ R T , indicating that the source is approximately symmetric. In the case of kaons, the results corresponding to the 2D case are shown in Fig. 26.
The N trk dependence of the 3D radii and the corresponding intercept parameters for pions (Eq. (10)), in several k T bins for all studied systems are shown in Fig. 27. The 3D pion radius parameters also show a similar pattern. The values of R L , R O , and R T are similar for pp and pPb collisions in all k T bins. In general R L > R S > R O , indicating once again that the source is elongated along the beam, which once more coincides with the behavior seen in pp collisions with the double ratio technique in Section 5.1. In the case of PbPb collisions, R L ≈ R O ≈ R S , indicating that the source is approximately symmetric. In addition, PbPb data show a slightly different N trk dependence. The R L radii extracted from the two-and three-dimensional fits differ slightly. While they are the same within statistical uncertainties for pp collisions at √ s = 0.9 TeV, the 3D radii are on average smaller by 0. 15   data is observed in R O , which could point to a different lifetime, τ, of the created systems in those collisions, since the outward radius is related to the emitting source lifetime by T , as discussed in Section 4.4.

Scaling for pions
The extracted radius parameters for pions are in the 1-5 fm range, reaching their largest values for very high multiplicity pPb collisions, as well as for similar multiplicity PbPb collisions, and decrease with increasing k T . By fitting the radii with a product of two independent functions of N trk and k T , the dependencies on multiplicity and pair momentum can be factorized. For that  Figure 24: k T -dependence of the 1D pion radius R inv , for several N trk bins, extracted for the pp and pPb systems. Lines are drawn to guide the eye.
purpose, we have used the following functional form: where the minimal radius a and the exponent γ of k T are kept fixed for a given radius component, for all collision types. This choice of parametrization is based on previous results [51], namely the minimal radius can be connected to the size of the proton, while the power law dependence on N trk is often attributed to the freeze-out density of hadrons.
We demonstrate the performance of the functional form of Eq. (15) by scaling one of the two variables (N trk or k T ) to a common value and showing the dependence on the other. Radius parameters scaled to k T = 0.45 GeV, that is R(k T /0.45 GeV) γ , as a function of N trk are shown in the left panels of Figs. 28 and 29. The ratio of the radius parameter over the value of the above parametrization at k T = 0.45 GeV, that is R/R param (N trk , k T = 0.45 GeV), as a function k T is shown in the right panels of Figs. 28 and 29. The obtained parameter values are listed in Table 8. The estimated uncertainty on the parameters is about 10%, based on the joint goodnessof-fits with the above functional form. Overall the proposed factorization works reasonably well.

Summary and conclusions
Detailed studies have been presented on femtoscopic correlations between pairs of same-sign hadrons produced in pp collisions at √ s = 0.9, 2.76, and 7 TeV, as well as in pPb collisions at √ s NN = 5.02 TeV and peripheral PbPb collisions at √ s NN = 2.76 TeV. The characteristics of two-particle Bose-Einstein correlations are investigated in one (R inv ), two (R T and R L ), and three (R L , R S , and R O ) dimensions, as functions of charged multiplicity N trk and pair average transverse momentum k T .
Two different analysis techniques are employed. The first "double ratio" technique is used to study Bose-Einstein correlations of pairs of charged hadrons emitted in pp collisions at √ s = 2.76 and 7 TeV, both in the collision center-of-mass frame and in the local co-moving system. The second "particle identification and cluster subtraction" method is used to study the characteristics of two-particle correlation functions of identical charged pions and kaons, identified via their energy loss in the CMS silicon tracker, in collisions of pp, pPb, and peripheral PbPb collisions at various center-of-mass energies. The similarities and differences between the three colliding systems have been investigated.
The quantum correlations are well described by an exponential parametrization as a function of the relative momentum of the particle pair, both in one and in multiple dimensions, consistent with a Cauchy-Lorentz spatial source distribution. The fitted radius parameters of the emitting source, obtained for inclusive charged hadrons as well as for identified pions, increase along with charged-particle multiplicity for all colliding systems and center-of-mass energies, for one, two, and three dimensions alike. The radii are in the range 1-5 fm, reaching the largest values for very high multiplicity pPb collisions, close to those observed in peripheral PbPb collisions with similar multiplicity. In the one-dimensional case, R inv in pp collisions steadily increases with the charged multiplicity following a N 1/3 trk dependence, as well as showing an approximate independence of the center-of-mass energy. This behavior is also observed for the longitudinal radius parameter (R * L in the center-of-mass frame, and R L in the local co-moving system) in the two-dimensional case.
The multiplicity dependence of R L and R T is similar for pp and pPb collisions in all k T bins, and that similarity also applies to peripheral PbPb collisions for k T > 0.4 GeV. In general, the observed orderings, R L > R T and R L R S > R O , indicate that the pp and pPb sources are elongated in the beam direction. For peripheral PbPb collisions, the source is quite symmetric, and shows a slightly different N trk dependence, with the largest differences between the systems for R T and R O , while R L and R S are approximately equal. The most visible differences between pp, pPb, and PbPb collisions are seen in R O , which could point to a different lifetime of the emitting systems in the three cases. The kaon radius parameters show a smaller increase with N trk than observed for pions. The radius parameters of charged hadrons and pions decrease with increasing k T , an observation compatible with expanding emitting sources. The k T -and N trk -dependences of the radius parameters factorize, and are largely independent of beam energy and colliding system for lower multiplicities, although some system dependence appears at higher values of N trk .
A shallow anticorrelation (a region where the correlation function falls below one) observed using the double ratio technique in a previous analysis of pp collisions, is also seen in the present data. The depth of this dip in the one-dimensional correlation functions decreases with increasing particle multiplicity and also decreases slightly with increasing k T .
Finally, the similar multiplicity dependencies of the radii extracted for pp, pPb, and peripheral PbPb collisions suggest a common freeze-out density of hadrons for all collision systems, since the correlation techniques measure the characteristic size of the system near the time of the last interactions. In the case of pp collisions, extending the investigation of these similarities to higher multiplicities than those accessible with the current data can provide new information on final state collective effects in such events.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: the Austrian Federal Ministry of Science, Research and Economy and the Austrian Science Fund; the Belgian