Bose-Einstein correlations in pp , pPb , and PbPb collisions at √ sN N = 0 . 9 – 7 TeV

Author(s): CMS Collaboration; Sirunyan, Albert M.; Bachmair, Felix; Bäni, Lukas; Berger, Pirmin; Bianchini, Lorenzo; Casal, Bruno; Dissertori, Günther; Dittmar, Michael; Donegà, Mauro; Grab, Christoph; Heidegger, Constantin; Hits, Dmitry; Hoss, Jan; Kasieczka, Gregor; Klijnsma, Thomas; Lustermann, Werner; Mangano, Boris; Marionneau, Matthieu; Meinhard, Maren T.; Meister, Daniel; Micheli, Francesco; Musella, Pasquale; Nessi-Tedaldi, Francesca; Pandolfi, Francesco; Pata, Joosep; Pauss, Felicitas; Perrin, Gaël; Perrozzi, Luca; Quittnat, Milena; Schönenberger, Myriam; Shchutska, Lesya; Tavolaro, Vittorio R.; Theofilatos, Konstantinos; Vesterbacka Olsson, Minna L.; Wallny, Rainer; Zagozdzinska, Agnieszka; Zhu, De Hua; et al.


I. INTRODUCTION
Studies of quantum-statistical correlations of pairs of identical particles produced in high-energy 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 (the 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 correlation (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 * Full author list given at the end of the article. 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. Second, 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, in pPb collisions at a nucleon-nucleon center-ofmass energy of √ s NN = 5.02 TeV, and in 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 threedimensional 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 Sec. II, while track selections and particle identification are detailed in Sec. III. In Sec. IV, the two analysis methods are described. A compilation of the results is presented in Secs. V and VI. Finally, Sec. VII summarizes and discusses the conclusions of this study.

II. 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 Large Hadron Collider (LHC) bunch structure and timing of the incoming beams.

A. 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 Secs. III A 1 and III A 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 (totaling 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) TABLE I. 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. An ellipsis indicates that there are not enough data to allow for a good-quality measurement.  [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 I. 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 I, two different sets of N rec1 ranges are shown. The upper five (N rec1 = 0-9, 10-14, 15-19, 20-29, and 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 I.
For the BEC analysis the standard CMS highPurity [9] tracks were used. The additional track selection requirements that were 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) 064912-2 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 remains 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, and 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 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 spiraling 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 Sec. III A 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 II. 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  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 (the full theoretical calculation is given in Eq. (33.11) in Ref. [24]).
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.

B. 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 [23]. 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 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 −2lnL 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-hand 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-hand plot. The plots for pp and PbPb interactions are very similar.

IV. THE BOSE-EINSTEIN CORRELATIONS
The correlation functions are investigated in one, two, and three dimensions in terms of q inv = ( p 1 − p 2 ) 2 − (E 1 − E 2 ) 2 , as well as in projections of the relative three-momentum 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-of-mass (c.m.) 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 c.m. 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,

A. 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 chargedparticle multiplicity N trk0.4 , and of the pair average transverse momentum k T . The distributions are determined in the c.m. 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. IV E 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 a 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, which is considered a systematic uncertainty [5,6]. Conversely, 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 III 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].

B. Particle identification and cluster subtraction technique
Similarly to the case of unidentified particles discussed in Sec. IV A, 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 064912-5 TABLE IV. Chosen bin widths for the various variables studied for pions and kaons in the particle identification and cluster subtraction method.

Variable
Range Bin width 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 IV. 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 Sec. IV C 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 • . 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 Sec. IV C. A joint functional form combining all effects is fitted in order to obtain the radius parameter and the correlation function intercept, as discussed in Sec. IV D.

Effect of Coulomb interaction and correction
The BEC method employed for the study of the two-particle correlations not only reflects 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 the case of same-sign and opposite-sign charges, is given by where ζ = αm/q inv is the Landau parameter, α is the finestructure 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 ζ 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 , 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) , the Coulomb correction factor K is well described by the formula (5) for q inv in GeV, and R in femtometers. 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 Coulomb-corrected C 2 q inv [GeV] Contribution of resonance decays to the measured Coulomb-corrected correlation function of π + π − (for 2 N rec2 < 10, top, blue squares) and K + K − (for 2 N rec2 < 60, bottom, blue squares). The φ(1020) peak (bottom panel) is off scale. low-momentum mini-jets, but can be also attributed to the effect of multi-body decays of resonances. In the following we refer to such an effect as being due to fragmentation of clusters, or cluster contribution [35][36][37]. For the purpose of evaluating and later eliminating that contribution, the one-dimensional (1D) OS correlation functions are fitted with the following form: where b and σ b are N rec2 -and k T -dependent parameters, and c OS is a normalization constant. The values of b and σ b are parametrized using  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.28ln( √ 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 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  4. Dependence of the relative amplitude z(N rec2 ) of the cluster contribution (between SS and OS pairs) as a function of N rec2 for two k T bins. Fit results (points with error bars) from all collision types are plotted together with combinatorics-motivated fits (solid curves) and their estimated systematic uncertainties (±0.2 bands). 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 minijet 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 square 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 Sec. VI. 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 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. 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 (top: 40 N rec2 < 50, bottom: 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.

D. Fitting the correlation function
The parametrizations used to fit the correlation functions in one, two, and three dimensions 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 twodimensional (2D) case it has the form and, finally, the form used for three dimensions is . (13) The radius parameters R inv , (R T ,R L ), and (R S ,R L ,R O ), correspond to the lengths of homogeneity in one, two, and three dimensions, 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-hand side, proportional to the fit variables δ inv , δ T , δ L , δ S , and δ 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. (13) T are the corresponding radius parameters in the three-dimensional (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 Sec. V C, the parameters resulting from the fits obtained with the simple and the stretched exponentials are 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  (11), obtained in the double ratio (solid 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. 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 1D 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 1D case. Therefore, the systematic 064912-10 uncertainties found in the 1D case were also applied to the 2D and 3D 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. 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 (left) and λ (right) as a function of k T obtained from an exponential fit to BEC data in three multiplicity bins in pp collisions at 2.76 TeV (open symbols) and 7 TeV (solid symbols). The boxes indicate the systematic uncertainties. Bottom: Radius parameter R inv as a function of N trk0.4 for pp collisions at 2.76 TeV (open squares) and 7 TeV (solid squares), where the systematic uncertainties are shown as error bars and the statistical uncertainties are smaller than the marker size. Previous results at 7 TeV (solid diamonds) [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  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 the case of the double ratio method above. The systematic uncertainties from various sources are added in quadrature.

F. Comparisons of the techniques
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 Sec. IV A, 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 Secs. IV B and IV C 2. Regarding the final-state Coulomb corrections, no significant difference is observed if applying the Gamow factor, Eq. (4), or the full Coulomb correction, Eq. (5), to the correlation functions in pp collisions. In both cases, the 1D correlations are fitted with an exponential function [a = 1 in Eq. (11)], from which the radius and intercept parameters are obtained. Since in the particle identification and cluster subtraction method N trk is fully corrected, as discussed in Sec. III A 2 and summarized in Table II, an additional correction is applied to the values of N trk0.4 in Table I 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 064912-12 The boxes indicate the systematic uncertainties. 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, top) 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 (bottom), 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.

A. One-dimensional results
Single and double ratios for pp collisions at √ s = 2.76 TeV are shown in Fig. 7, integrated over the chargedparticle 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 V, together with those for pp collisions at √ s = 7 TeV, these being compatible within the uncertainties at these two energies.
The top 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 V 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 Relativistic Heavy Ion Collider (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 I. 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 of 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 (Secs. V B and V C).
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)   in Fig. 10. The top 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 bottom 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. 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.   7 TeV,integrated over N trk0.4 and k T . Top: 2D projections over the (q L ,q S ), (q O ,q S ), and (q O ,q L ) planes, for q O < 0.05 GeV (left), q L < 0.05 GeV (middle), and q S < 0.05 GeV (right). Bottom: 1D projections of the same data, with the complementary two variables constrained to be within the first bin (q i,j < 0.05 GeV). The error bars show the corresponding statistical uncertainties. All plots are zoomed along the correlation function axis with the BEC peak cut off above 1.2.

B. 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 c.m. frame and in the LCMS. For clarity, an asterisk (*) is added to the variables defined in the c.m. frame when showing the results in this section. In the case of longitudinally symmetric systems, the LCMS has the advantage (Sec. IV D) 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 top panel shows the 2D plot of the double ratio in terms of q T and q L , with the right-hand plot showing the 2D Lévy fit (with a = 1) defined in Eq. (12), superimposed on the 2D correlation function. The bottom 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 those 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 one dimension, the correlation functions are also investigated in two dimensions in terms of the components (q ( * ) T ,q ( * ) L ), in three intervals of  k T , and in different N trk0.4 bins. The results from the stretched exponential fit to the double ratios are compiled in Fig. 12, performed both in the c.m. frame (top) and in the LCMS (bottom). The behavior is very similar in both frames, clearly showing that all fit parameters R ( * ) L and R ( * ) T increase with increasing multiplicity N trk0.4 . However, the behavior with respect to k T varies, showing low sensitivity to this parameter for the lower N trk0.4 bins. The radius parameters R ( * ) L and R ( * ) T 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 (c.m.) for the same bins of N trk0.4 and k T can be related to the Lorentz contraction of the longitudinal radius in the c.m. frame. Figure 13 shows the results for the correlation strength parameter, λ ( * ) , obtained with the same fit procedure, for both the c.m. 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 VI Table VI 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 and 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.

C. Three-dimensional results
The double ratios can be further investigated in three dimensions as a function of the relative momentum components (q O ,q S ,q L ). 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. For clarity, an asterisk (*) is added to the variables defined in the c.m. frame when showing the results in this section. An illustration of the 2D projections is shown in the top panel of Fig. 15, with data from pp collisions at 7 TeV, integrated in all N trk0.4 and k T ranges. Similarly to Fig. 14, the results are zoomed along the 3D double ratio axis, with the BEC peak cut off around 1.2 for better visualization of the directional behavior of the correlation function. The higher values seen around (q O ,q S ) ∼ (2,2) GeV in the upper middle plot are again artifacts of cutting the double ratio at 1.2, this value being much smaller than the true value of the peak in the low (q O ,q S ) region, and correspond to fluctuations in the number of pairs in those bins.
The bottom 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 3Dstretched 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 064912-19  TeV as a function of q inv or of the combined momentum, in two selected N rec2 bins (left and right), for pion pairs (red triangles) in one (top), two (middle), and three (bottom) dimensions, corrected for Coulomb interaction and cluster contributions (mini-jets and multi-body resonance decays). For better visibility, only a fraction of the points is plotted, and those with statistical uncertainty higher than 10% are in light grey color. The solid curves indicate fits with the stretched exponential parametrization.
the stretched exponential fit in the LCMS and integrating over all N trk0.4 and k T ranges, are summarized in Table VII, 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 064912-21 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 c.m. 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 c.m. 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 c.m. 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 Sec. IV D), which is expected to be higher in the c.m. 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 c.m. 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 22. Track-multiplicity (N trk ) dependence of the 1D pion radius parameter R inv (top) and intercept parameter λ (bottom), obtained from 1D fits (integrated over all k T ) for all collision systems studied. For better visibility, only every second measurement is plotted. Lines are drawn to guide the eye.
with increasing average multiplicity, N trk0. 4 , similar to what is seen in the 1D and 2D cases.
The intercept parameter λ ( * ) is also studied as a function of k T and N trk , both in the c.m. 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 c.m. frame, where the central R L value in three dimensions is roughly 0.3 fm larger than the corresponding R L value in two dimensions, 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. 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 1 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, top  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 threedimensional 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 charged-particle 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 relative amplitude z of the cluster contribution; and low q exclusion) are given by the solid 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).

A. 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 Sec. V A. 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 064912-24  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 Sec. V A. 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 Sec. V A. 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

B. 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 purpose, we have used the following 064912-26  Left-hand column, from top to bottom: R inv from the 1D (q inv ) analysis, R L and R T from the 2D (q L ,q T ) analysis. Right-hand column, from top to bottom: R L , R O , and R S from the 3D (q L ,q O ,q S ) analysis. Fit results are indicated in the figures, see text for details. functional form: R param (N trk ,k T ) = a 2 + bN β trk 2 1/2 (0.2 GeV/k T ) γ , (15) 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,    left-hand 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-hand panels of Figs. 28 and 29. The obtained parameter values are listed in Table VIII. The estimated uncertainty on the parameters is about 10%, based on the joint goodness of fit with the above functional form. Overall the proposed factorization works reasonably well.

VII. SUMMARY AND CONCLUSIONS
Detailed studies were 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 dependencies 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 1), 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 onedimensional 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