Centrality dependence of pion freeze-out radii in Pb-Pb collisions at $\sqrt{\mathbf{s_{NN}}}$=2.76 TeV

We report on the measurement of freeze-out radii for pairs of identical-charge pions measured in Pb--Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV as a function of collision centrality and the average transverse momentum of the pair $k_{\rm T}$. Three-dimensional sizes of the system (femtoscopic radii), as well as direction-averaged one-dimensional radii are extracted. The radii decrease with $k_{\rm T}$, following a power-law behavior. This is qualitatively consistent with expectations from a collectively expanding system, produced in hydrodynamic calculations. The radii also scale linearly with $\left<\mathrm{d}N_{\rm ch}/\mathrm{d}\eta \right>^{1/3}$. This behaviour is compared to world data on femtoscopic radii in heavy-ion collisions. While the dependence is qualitatively similar to results at smaller $\sqrt{s_{\rm NN}}$, a decrease in the $R_{\rm out}/R_{\rm side}$ ratio is seen, which is in qualitative agreement with specific predictions from hydrodynamic models. The results provide further evidence for the production of a collective, strongly coupled system in heavy-ion collisions at the LHC.


Introduction
Femtoscopy in heavy-ion collisions is understood in some detail; for example see the experimental overview in [20] and model calculations in [27][28][29][30]. The dependence of the system size extracted from the data is investigated as a function of collision centrality and average transverse momentum of the pair k T = | p 1,T + p 2,T |/2. As the initial size of the system grows with increasing multiplicity (decreasing centrality), so does the apparent system size at freeze-out, measured by femtoscopy. Such increase is naturally produced in a hydrodynamic calculation. Strong hydrodynamic collective flow in the longitudinal and transverse directions results in the decrease of the apparent size of the system with increasing k T . This is because longitudinal and transverse velocity boosts cause particles emitted from spatially separated parts of the collision region to move away from one another. Such particles cannot have a small momentum difference, and so correlation functions of boosted particles are sensitive to only part of the collision region. This part is referred to as the "homogeneity length" [31]. The decrease of the size with k T is observed in experimental data from heavy-ion collisions at all centralities, various collision energies and colliding system types, and is well described quantitatively in hydrodynamic models [11,30] and qualitatively in hadronic rescattering codes [32].
Taking into account the successful description of the femtoscopic scales at lower energies, the hydrodynamic modelling has been extrapolated to collision energies of the LHC [30,[33][34][35]. The expected increase in initial energy density (temperature) leads to larger evolution times, which in turn produce larger overall system size and stronger transverse and longitudinal flows. At the same time the freeze-out hypersurface evolves to have significant positive space-time correlation. This influences the radii of the system in the plane perpendicular to the beam axis. In particular the radius along the pair transverse momentum (called R out ) is decreased by the correlation with respect to the other transverse radius (called R side ), which decreases the R out /R side ratio All of those effects have been observed in the first measurement for central (0-5%) collisions at the LHC [36]. This work extends this measurement to other centralities and compares the obtained radii to recent hydrodynamic calculations, in order to check their validity in a large range of event multiplicities. A measurement of one-dimensional radii was also performed using the two-pion and three-pion cumulants [37]. This work extends the two-pion measurement to several ranges of pair transverse momentum.
The paper is organized as follows. In Section 2 the data taking conditions and data treatment is described. In Section 3 we give the details of the analysis of the correlation function. In Section 4 the extracted radii are presented and compared to model expectations, while Section 5 summarizes our findings.

Data taking and track reconstruction
This work reports on the analysis of Pb-Pb collisions produced by the LHC during the 2010 datataking period. They were recorded by the ALICE experiment; the detailed description of the detector and the performance of all of its subsystems is given in [38,39]. Here we only briefly describe the specific detectors used in this analysis. The ALICE Time Projection Chamber (TPC) [40] is a large volume gaseous ionization chamber detector, which was used both for tracking at mid-rapidity as well as for particle identification via the measurement of the specific ionization energy loss associated with each track. In addition to the TPC, the information from the ALICE Inner Tracker System (ITS) was used. The ITS consists of six cylindrical layers, two Silicon Pixel Detectors closest to the beam pipe, two Silicon Drift Detector layers in the middle, and two Silicon Strip Detectors on the outside. The information from ITS was used for tracking and primary particle selection, as well as for triggering. However the main triggering detector was the V0. It is a small angle detector consisting of two arrays of 32 scintillating counters. The first (V0A) is located 330 cm away from the vertex and covers 2.8 < η < 5.1, the second (V0C) is fixed at the front of the hadronic absorber of the muon arm and covers −3.7 < η < −1.7. The tracking detectors are located inside the solenoidal ALICE magnet, which provides a uniform magnetic field of 0.5 T along the beam direction. The T0 detector [41] was a main luminometer in the heavy-ion run. It consists of two arrays of 12 Cherenkov counters, covering −3.28 < η − 2.97 and 4.61 < η < 4.92. It has a time resolution of 40 ps.
The minimum-bias trigger required a signal in both V0 detectors, which was consistent with the collision occurring at the center of the ALICE apparatus. A total sample of approximately four million Pb-Pb events was used for this analysis. For the centrality range considered in this work the trigger efficiency was 100%. The centrality was determined by analyzing the signal from the V0 detector with the procedure described in details in [5]. This ensured that the centrality determination was obtained using particles at significantly different rapidities than the ones used for the pion correlation analysis. This work presents results for seven centrality ranges: 0-5%, 5-10%, 10-20%, 20-30%, 30-40%, 40-50%, and 50-60% of the total hadronic cross section. The position V z of the event vertex in the beam direction with respect to the center of the ALICE apparatus was determined for each event. In order to ensure uniform pseudorapidity acceptance, only events with |V z | < 8 cm were used in the analysis.
For each event, a list of tracks identified as primary pions was created, separately for positive and negative particles. Each track was required to leave a signal both in the TPC and the ITS and the two parts of the track had to match. The TPC is divided by the central electrode into two halves, each of them is composed of 18 sectors (covering the full azimuthal angle) with 159 padrows placed radially in each sector. A track signal in the TPC consists of space points (clusters), each of which is reconstructed in one of the padrows. A track was required to be composed out of at least 80 such clusters. The parameters of the track are determined by performing a Kalman fit to a set of clusters. The quality of the fit is judged by the value of χ 2 , which was required to be better than 4 per cluster (each cluster has two d.o.f.). The transverse momentum of each track was determined from its curvature in the uniform magnetic field. Two opposite field polarities were used through the data-taking period, for a check of systematic tracking effects. The momentum from this fit in the TPC was used in the analysis. Tracks which had a kink in the trajectory in the TPC were rejected. Trajectories closer than 3.2 cm in the longitudinal direction and 2.4 cm in the transverse direction to the primary vertex were selected to reduce the number of secondaries. The kinematic range for accepted particles was (0.14, 2.0) GeV/c in transverse momentum and (−0.8, 0.8) in pseudorapidity. Based on the specific ionization energy loss in the TPC gas dE/dx, a probability for each track to be a pion, kaon, proton, and electron was determined after comparing with the corresponding Bethe-Bloch curve. Particles for which the pion probability was the largest were used in this analysis. This resulted in an overall purity above 95%, with small contamination from electrons in the region where the dE/dx curves for the two particle types intersect.
The accepted particles from each event are combined into same-charge pairs. The two-particle detector acceptance effects must be taken into account in this procedure. Two main effects are present: track splitting and track merging. Track splitting occurs when a single trajectory is mistakenly reconstructed as two tracks. The tracking algorithm has been specifically designed to suppress such cases. In a rare event when splitting happens, two tracks are reconstructed mostly from the same clusters in the TPC. Therefore, pairs which share more than 5% of clusters are removed from the sample. Together with the anti-merging cut described below this eliminates the influence of the split pairs. Two-particle correlated efficiency and separation power is affected by track merging. In the TPC, two tracks cannot be distinguished if their trajectories stay close to each other through a significant part of the TPC volume. Although this happens rarely, such pairs by definition have low relative momentum and therefore their absence distorts the correlation function in the signal region. The effect of track merging has been studied in central collisions in the previous work [36]. In this work we have used a similar procedure to correct for the merging effects, through dedicated two-particle selection criteria. More details are given in Sec. 3.3.

Correlation function analysis
The two-particle distribution for same-event pion pairs depends on several factors, including trivial single-particle acceptance effects. To extract only the relevant two-particle correlation effects, the correlation function formalism, described below, is applied.

Correlation function construction
The femtoscopic correlation function C is constructed experimentally as where q = p 1 − p 2 is the pair relative momentum (due to fixed masses of the particles only three components are independent). The magnitude of this vector is referred to as q inv . For a detailed description of the formalism, see e.g. [42]. The signal distribution A is composed of pairs of particles where both come from the same event. The background distribution B is constructed with the "mixing" method, in which the two particles come from different events, which must have similar characteristics, so that their single-particle efficiency and distribution are as close as possible. To form a "mixed" pair, particles must come from two events, for which the centralities differ by no more than 2.5% and vertex positions differ by no more than 4 cm. The correlation function is normalized with the ratio of the number of pairs in the B and A samples in the full q range used (0 − 0.25 GeV/c), so that C at unity means no correlation. The dependence on pair momentum sum is investigated by doing the analysis for various ranges of k T , namely: 0. The momentum difference q is calculated in the Longitudinally Co-Moving System (LCMS), in which the pair total longitudinal momentum vanishes: p 1,L + p 2,L = 0. The three coordinates of q in LCMS are defined as follows: long -along the beam axis, out -along the pair transverse momentum, and sideperpendicular to the other two. In Fig. 1 the projections of three example correlation functions along these axes are shown. A significant, approximately Gaussian enhancement at low relative momentum is seen in all projections. The width of the correlation grows with increasing centrality (lowering multiplicity) as well as with increasing k T .  The pair distributions for identical particles have specific symmetries, which are naturally represented in a spherical harmonic (SH) decomposition [43,44]. In particular for identical particles all odd-l and odd-m components of a correlation function representation vanish. The gross features of the correlation function which are relevant for femtoscopy are fully captured by the low-l components of the decomposition. The l = 0, m = 0 component is sensitive to the overall size of the pion source. The l = 2, m = 0 component is sensitive to the difference between longitudinal and transverse size, and the l = 2, m = 2 component reflects the difference between sideward and outward component of the transverse radius. Therefore, those three components of the SH representation contain the same information as the Cartesian one for the purpose of the femtoscopic analysis. In particular both representations allow for fitting of the correlation function with the same theoretical formula. The next non-vanishing components are for l = 4. Their analysis is beyond the scope of this paper, which focuses on the overall width (variance) of the distribution in three directions.
In Fig. 2 we show the first three non-vanishing components of the spherical harmonics representation of three example correlation functions, the same as in Fig. 1. In the (0, 0) component the enhancement at low-q is clearly visible and its width is increasing with centrality and k T . In this representation the low-q dip coming from Coulomb repulsion is visible better than in the Cartesian one. For large k T , there is a lack of low-q pairs, as a result of track merging (see below). Nevertheless, the correlation clearly extends well beyond the region with low statistics. The other two components: (2, 0) and (2, 2), show a non-trivial correlation structure deviating from zero, indicating that the shape is not spherical, but rather ellipsoidal in LCMS. The width of these structures is consistent with the width of the correlation enhancement in the (0, 0) component, which means that all three reflect the properties of the femtoscopic signal. The analysis of the shape of these structures is the main focus of the next section.
If the available statistics is limited, as is sometimes the case for lower collision energies or particles heavier than pions, the analysis is performed only as a function of magnitude of relative momentum q inv , most naturally calculated in the Pair Rest Frame (PRF). In this work we present results in this variable for completeness. Data was analyzed in the same centrality and pair k T ranges as the ones used in the three-dimensional analysis.

Fitting the correlation function
The space-time characteristics of the source are reflected in the correlation function. They are connected via the Koonin-Pratt equation where r is the pair space-time separation four-vector. S is the source emission function, interpreted as a probability to emit a pair of particles with a given relative momentum and space-time separation. Ψ is the two-particle interaction kernel. In the simplest case of non-interacting particles (e.g. photons) it is the modulus of the pair wave-function. If the Coulomb or strong interaction between the particles (called Final State Interaction -FSI) needs to be taken into account, then Ψ becomes the Bethe-Salpeter amplitude corresponding to the solution of the relevant quantum scattering problem, taken with the inverse time direction [45].
Previous studies at RHIC [16-18, 22, 23, 46, 47] and at the LHC [36] have approximated the source by a Gaussian, treating any difference between the real data and a Gaussian as a correction. This procedure was also universally used in all past pion femtoscopic analyses of heavy-ion collisions. Therefore, we also use it here by writing where r out , r side , and r long are components of the relative separation r. This static form of S is expressed in LCMS, with R out , R side , and R long being the single-particle source sizes of the system later referred to as "femtoscopic radii", or simply "radii". They quantify the lengths of homogeneity of the system in the outwards, sidewards, and longitudinal direction, respectively.
For like-sign pions the strong interaction contribution is small for the source sizes expected here (a few fm) [42], so it is neglected. The remaining Ψ is a convolution of the Coulomb interaction and wave-function symmetrization. As an approximation, the Coulomb part is factorized out and integrated separately in the procedure known as the Bowler-Sinyukov fitting [48,49]. It is well tested and is applicable for pions and for large source sizes expected in this analysis. In this approximation the integration of Eq. (2) with S given by Eq. (3) gives the following fit form for the correlation function where N is the overall normalization factor. The function K C (q inv ) is the Coulomb part of the two-pion wave-function integrated over the spherical Gaussian source with a given radius. For each correlation function it is set to the value from the one-dimensional analysis (see Section 4.3), to reflect the decrease of the source size with multiplicity and k T . Its variation is a source of systematic uncertainty. The dilution parameter λ is introduced to account for the fact that not all measured pion pairs are correlated, and that the real emission function may deviate from a Gaussian form.
The fit is performed with the log-likelihood method for the three dimensional correlation function in Cartesian representation, resulting usually in several thousand degrees of freedom. Examples are shown in Fig. 1. The Gaussian fit is able to reproduce the overall width of the correlation in all cases. Some details of the behavior at low-q may not be perfectly described, which can be attributed to the limitations of the Bowler-Sinyukov formula as well as to the non-Gaussian, long-range tails which may be present in the source. Some deviations from the Gaussian ellipsoid shape for the higher centrality can also be seen for the long direction. We leave the detailed investigation of these effects for future work. Nevertheless, the overall sizes R of the system, which are mostly sensitive to the width of the correlation are well estimated. The deviation of the correlation function from the pure Gaussian shape is smaller than a similar deviation in pp collisions [50].
An equivalent fit is also performed for the SH representation of the correlation. Eq. (4) is numerically integrated on a ϕ − θ sphere for each q LCMS bin, with proper Y m l weights, to produce the 3 components of the SH decomposition. Statistical uncertainties on each component are taken into account, as well as the covariance matrix between components. Examples are shown in Fig. 2. The fit describes the general direction-averaged width of the correlation function, shown in the upper panel. Small deviations can be seen in the shape and the behavior of the fit at small q, as discussed earlier for the Cartesian fit. The deviations from zero in the (2, 0) and (2, 2) components are small but statistically significant for a large number of multiplicity/k T ranges, indicating that the source size is slightly different in all three directions.
For the one-dimensional correlation functions constructed as a function of q inv , the fit was performed with a simplified version of Eq. (4) where the only radius parameter is R inv -the one-dimensional direction-averaged femtoscopic radius in the PRF. Table 1 lists the systematic uncertainty contributions. The range of values is given to provide a general estimate of the importance of each contribution, however the systematic uncertainty is estimated for each point individually. Separate analyses were performed for positive and negative pions, as well as for two datasets collected with opposite polarities of the magnetic field inside ALICE. This results in four independent data samples for which certain systematic effects, most notably the single-track inefficiencies, are different. Correlation functions for all four samples for all centrality and pair momentum ranges are statistically consistent, after all corrections are applied; this is an important systematic cross-check of the methodology. In the following discussion the central values are statistical averages of the fit values obtained for the four samples. The systematic uncertainty arising from differences among the data sample is between 1 and 2% for all radii. The other systematic uncertainties are analyzed for each sample separately; their final value is the convolution of the uncertainties for each sample.

Systematic uncertainties of the radii
Two correlation function representations are used in this work: the Cartesian and spherical harmonics. They are mathematically equivalent, the fitting procedure used the same functional form for both. However the implementation of the fitting procedure is quite different: log-likelihood vs. regular χ 2 fit, three-dimensional Cartesian histogram vs. three one-dimensional histograms, fitting range as threedimensional cube in q out , q side , q long or a three-dimensional sphere with constant q LCMS radius among others. Therefore, the fits to the two representations differ systematically upon variation of the fitting procedure (fit ranges, Bowler-Sinyukov approximation, etc.). The difference between the values for the two fits is taken as a part of the systematic uncertainty. It usually ranges from 1% to 3% and grows with pair k T and multiplicity.
Variation of single particle cuts around the default value results in modifications of single particle acceptances and purities. However the correlation function shape should be to the first order insensitive to those effects. We have checked that for a reasonable variation of the single particle cuts, the resulting radii are consistent within statistical uncertainties.
The measurement of the average event multiplicity for a given centrality range has a known uncertainty of 3-4% for all centrality classes [5]. The femtoscopic radii in heavy-ion collsions were observed to scale linearly with dN ch /dη 1/3 at lower collision energies [20]. Such trend is also predicted by hydrodynamical models and is expected to hold at LHC. Therefore the systematic uncertainty coming from the multiplicity estimation is about 1%.
The dominant systematic effect on the two-particle correlation function is the two-track correlated efficiency. This effect has been studied in previous work [36] for central collisions. The effect of "splitting" is fully removed with the help of dedicated two-track selection criteria, mentioned in Sec. 2, and does not influence the radii. However when two particles' trajectories are located close to each other in the volume of the TPC detector, they may be reconstructed as one track or not be reconstructed at all. This effect, called "merging" in the following, results in a loss of reconstruction efficiency for such pairs of tracks. Pairs of primary pion trajectories close in space correspond to low relative momentum, therefore "merging" will affect the pion correlation function exactly in the femtoscopic signal region. The two-track efficiency was studied in Monte-Carlo simulation of the ALICE detector. For pion pairs at low relative momentum an efficiency loss of up to 20% was observed. A two-track selection was chosen, such that the resulting correlation function was not affected by the inefficiency. "Merging" affects only pairs which are spatially close in the detector. The "closeness" can be quantified by the pseudorapidity difference ∆η for the pair and only pairs with |∆η| < 0.016 are affected. The trajectories must also be close in the transverse plane, where they are curved by the magnetic field. The azimuthal coordinate ϕ * of tracks at a radius of 1.2 m from the collision point (so roughly in the center of the TPC volume) is calculated. Pairs with ∆η 2 + ∆ϕ * 2 < 0.045 are affected by merging. Pairs which simultaneously satisfy this and the previous angular criteria are removed both from the signal and the background sample. The correlation function calculated for pairs surviving this cut is then not affected by "merging". As a systematic check a different procedure to calculate ϕ * was used, where the value was taken not at a fixed radius, but instead at a radius inside the TPC (i.e. between 0.5 and 2.4 m) where ϕ * was the smallest. The differences between radii for the correlation functions calculated with the two methods is taken as systematic uncertainty. It is from 2-13%.
A pair of same charge pions traversing a solenoidal field can have two configurations: a "sailor" where the two trajectories bend away from one another and never cross, or a "cowboy" where the trajectories bend toward each other and can cross [18]. Merging affects the "cowboy" pairs, but only weakly influences the "sailors". The two configurations correspond to different phase space regions in the Cartesian representation of the correlation function. We have performed two fits to the correlations (corrected for merging), restricting the fitting range to either the "cowboy" or the "sailor" region. The consistency of the radii obtained in these fits was used as an estimator of the effectiveness of the cut for removing merged tracks. A less restrictive cut degraded agreement between the "cowboy" and "sailor" radii, while making the cut stricter did not improve the agreement but reduced statistics in the signal region. The comparison allowed optimization of the cut, and provided another estimate of the systematic uncertainty on the two-track correction procedures. The uncertainty estimated in this way is consistent with the un-certainty determined by varying the ϕ * definition. This uncertainty is largest for transverse radii, is most prominent for high multiplicity events (central collisions), and affects a wider region in q for pairs with higher k T .
The separation into the "cowboy" and "sailor" phase-space regions is not feasible for the SH representation of the correlation function. In this case, when the anti-merging cut is not properly applied, one observes significant non-femtoscopic signals, especially in the C 2 2 component of the correlation function and a reliable fit cannot be performed. Therefore, the C 2 2 component serves as a sensitive independent check of the effectiveness of the "anti-merging" cut. This is a good illustration on how the two representations complement each other in the systematic studies.
The fit was performed for several values of the fitting range in q (varying with multiplicity and pair k T , following the changes in the correlation function width). The variation of the fitted radii with the change of the range was taken as another component of the uncertainty, which is less than 5% for all the radii.
In addition to the uncertainties listed above, other systematic effects can influence the extracted radii. The first is the momentum resolution, which was studied in [36]. The correction procedure described there is used in this work, as well. The uncertainty on the radii from this correction is 2%. Another effect is the influence of the Bowler-Sinyukov procedure on the extracted radii and λ parameter (fraction of correlated pairs). The procedure results in an uncertainty of 3% on R out , 1% on the other radii, and lowers the λ value by up to 5% [28].
All the systematic uncertainty components mentioned above are added in quadrature, the range of values of the total systematic uncertainty is given in Table 1.
The one-dimensional analysis is performed in PRF, where the total momentum of the pair vanishes. In the transformation from LCMS to PRF R out is scaled by the γ factor for the pair, depending on k T and as a consequence is larger than the other two components. In such case the direction-averaged onedimensional correlation function becomes non-Gaussian. This produces a dependence of the fit value on the range of the fit, resulting in a systematic uncertainty of up to 10%. Other components of the uncertainty, such as field orientation dependence, momentum resolution and Coulomb correction dependence are comparable to those from the three-dimensional fit.

Three-dimensional radii
The outcome of the fitting procedure are 49 sets of femtoscopic radii, one set for each centrality and k T range. They are shown in panels a)-c) of Fig. 3. The radii in all directions are in the range of 2 to 8.5 fm. The radii universally decrease with increasing k T , in qualitative agreement with a decreasing homogeneity length, as predicted by hydrodynamics. Such behavior is a strong indication of a large degree of collectivity in the created system. The radii are also universally higher for more central collisions, which correspond to growing final-state event multiplicity. For the lowest k T , R long is generally the largest, whereas at large k T there is no universal ordering of the radii.
In the most central collisions the values of the λ parameters are around 0.33 for the lowest k T range, are increasing linearly to 0.43 for k T of 0.65 GeV/c, and falling to 0.36 in the last k T range. For peripheral collisions, they are higher by 0.05 to 0.07. The lowering of λ from a theoretical maximum of 1.0 can be attributed to a number of factors. The pion sample contains daughters of long-lived strongly decaying resonances; their fraction falls with growing p T . As a result the source function contains non-Gaussian tails extending to large relative separation. There may be additional factors influencing the shape of the correlation function, coming from the source dynamics. The exact source shape usually deviates from a Gaussian, in a way that lowers the λ parameter of the Gaussian fit significantly. The detailed study of this shape requires a dedicated methodology and is beyond the scope of this work. In addition, at the largest k T the electron and pion dE/dx become comparable in the TPC, and the pion sample is contaminated by electrons, which also lowers λ. The approximate treatment of the Coulomb interaction in the Bowler-Sinyukov fitting procedure lowers it by 5-10% [28]. This effects is most pronounced for large source sizes (central collisions). Finally, possible coherent emission of pions [51] is expected to lower λ by a few percent. The λ parameter values given here are for reference only and their physics interpretation is not discussed.
In panel (d) of Fig. 3 the R out /R side ratio is shown. Its systematic uncertainty is determined independently from those of R out and R side , to account for the fact that they may be correlated. The ratio is consistent with unity for central collisions. Its value slowly decreases for more peripheral collisions, and reaches 0.85 for peripheral collisions and high k T . Based on hydrodynamic models, the R out /R side ratio has been proposed as a sensitive probe of the shape and space-time correlation present at the freeze-out hypersurface [33,52]. In particular this ratio at the LHC was predicted to be lower than a value of 1.1 measured at top RHIC collision energies.

Scaling of the radii
It has been argued in [20] that the femtoscopic volume scales with the final-state event multiplicity, and that each of the three-dimensional radii separately scales with this value taken to the power 1/3. In Fig. 4 we present the dependence of the radii on multiplicity for Pb-Pb collisions. The scaling is evident for all datasets, for all three directions and all analyzed pair momentum ranges.
Similarily, hydrodynamics predicts approximate scaling of the radii with pair transverse mass m T = k 2 T + m 2 π [35]. The slope parameters of the lines shown in Fig. 4 are plotted in Fig. 5  where β and α are free parameters. The slope parameters follow the power law scaling within the current systematic uncertainties, the value of the α parameter is −0.65 ± 0.12 for long, −0.46 ± 0.13 for out, and −0.52 ± 0.11 for side direction. The dependence of the values of femtoscopic radii on centrality and k T factorizes into a linear dependence on dN ch /dη 1/3 and a power law dependence on m T . The results of the one-dimensional fits with Eq. (6) are shown in Fig. 6. Similarly to the three-dimensional case, the radius is increasing with event multiplicity (decreasing centrality). Therefore, the final state shape is reflecting the growth of the initial shape with decreasing centrality. The R inv is also decreasing with pair transverse momentum. This is usually understood as a manifestation of the hydrodynamic collectivity. The one-dimensional radius also serves as a comparison basis with the femtoscopic analysis for heavier particles, where the one-dimensional analysis is standard and the three-dimensional analysis  is challenged by the more complicated description of the pair interaction as well as significantly smaller statistics [54, 55].

Comparison to previous measurements
In Fig. 7 the heavy-ion data from Pb-Pb collisions at the LHC reported in this work are compared to the previous measurements, including results obtained at lower collision energies. It has been argued [20] that three-dimensional femtoscopic radii scale with the cube root of measured charged particle multiplicity not only for a single energy and collision system, but universally, across all collision energies and initial system sizes. The dashed lines in the figure are linear fits to heavy-ion data available before the startup of the LHC (the dotted lines show one sigma contours of these fits). At lower energies the linear scaling was followed well in long and side directions and only approximately in out, with some outliers such as the most peripheral collisions reported by STAR. Our data at higher collision energy show that the scaling in the long direction is preserved. The data for the side direction fall below the scaling trend, although still within the statistical uncertainty. A clear departure from the linear scaling is seen in the out direction; data from the LHC lie clearly below the trend from lower energies. Such behavior was predicted by hydrodynamic calculations [33] and is the result of the modification of the freeze-out shape. Larger initial deposited energy produces larger temperature gradients and longer evolution time at LHC. This results in a change from outside-in to inside-out freeze-out and this modification of the space-time correlation drives the R out /R side ratio to values lower than at RHIC. Therefore, already for heavy-ion data the simple universal linear scaling is broken in the transverse direction. As already reported in [50], the femtoscopic radii for pp collisions also exhibit linear scaling in the same variable, albeit with significantly different parameters. In this case the scaling between different colliding systems is broken again in the longitudinal direction.
Other scaling variables were proposed for the femtoscopic radii, based on Monte-Carlo simulations of the initial state, such as average number of participants N part [46] or the characteristic initial transverse sizē R [56, 57]. If data from different centralities but the same collision energy are considered, a linear scaling is indeed observed for radii in all directions, if they are plotted as a function of either of these variables. However, no such scaling is observed when data from two collision energies (e.g. top RHIC and top LHC) is compared. In that sense these variables are less adequate than dN ch /dη 1/3 , for which the linear scaling is preserved, across many collision energies, for at least one direction (R long ). This observation is consistent with an expectation that the final freeze-out volume, reflected in the femtoscopic radii, should scale with the final-state observable (such as e.g. dN ch /dη 1/3 ), while the simple "geometric" initial state variables do not contain enough information. They must be replaced with modeling of the full evolution process, such as the one given by hydrodynamics, which depends on additional parameters apart from the initial size, such as initial energy density.

Model comparisons
The hydrodynamic models predict both the centrality and pair momentum dependence of the femtoscopic radii. Usually the parameters of the model (initial energy density profile, the equation of state and the freeze-out condition) are adjusted to reproduce the shape of the single particle inclusive transverse momentum spectra as well as the elliptic flow. The charged particle multiplicity must also be reproduced, preferably as a function of pseudorapidity if the model employs full three-dimensional modeling. The total particle multiplicity is usually determined on the freeze-out hypersurface by employing statistical hadronization. After converting the continuous medium to hadrons, final state interactions are taken into account either in the simplified form, with propagation and decay of hadronic resonances or with the full rescattering simulation. We compare our results to predictions from the Therminator model coupled to (3+1)D viscous hydrodynamics [34,35,58]. Similar results have been obtained in the HKM model [30].
In Fig. 8 we show the comparison of our data to the calculations from the (3+1)D hydrodynamic model coupled to the Therminator statistical hadronization code. The model is fully three-dimensional and it is able to reproduce values of R long for all centralities, with some overprediction of the overall magnitude and the slope of the k T dependence, especially at low momentum. This is an indication that the longitudinal dynamics is reasonably described in the model, both in momentum and space-time sectors. The R out is well described for all centralities, both in magnitude as well in the slope of the k T dependence. The slope of the k T dependence is also well described for R side , but the magnitude is lower than in data, although within the systematic uncertainty. The intercept of R side at low k T is usually associated with the overall geometrical size of the system, while the slope of the k T dependence of both transverse radii depends on the amount of flow in the system. Both are well reproduced, so the hydrodynamic approach is in good agreement with our data.
Another model based on hydrodynamic formalism, the HKM [30] is also shown in Fig. 8 for 0-5% most central collisions. It differs from the previous in the implementation of the freeze-out process. It also directly treats hadron rescattering with the UrQMD simulation. Nevertheless, the pion femtoscopy in central Pb-Pb collisions at the LHC is reproduced in the calculation, with some underprediction of R side and R long magnitude. Therefore, the approximate agreement with data is a universal feature of such models, not of a particular implementation. The particular choice of initial conditions and the equation of state for these models was motivated by the analysis of RHIC femtoscopic data. This choice is essential for the correct description of the data, due to the collective flow mechanism in the hydrodynamic evolution at the LHC energies. It also indicates that the details of the freeze-out process have limited influence on pion femtoscopy. Some studies suggest that femtoscopy of heavier particles might be a more sensitive probe in this case [30].
The data to model comparison of the R out /R side ratio is plotted in Fig. 9. It shows values consistent with unity for central collisions, both for models and data. Such low values are associated with the change to outside-in freeze-out scenario at LHC collision energies. For most peripheral collisions the ratio decreases with k T even more, to values smaller than unity. This decrease is qualitatively reproduced in the model, although the calculations are at the upper edge of the experimental systematic uncertainty.  Fig. 9: Comparison of the R out /R side ratio, as a function of pair transverse momentum, with the calculation from the Therminator and (3+1)D hydro model [35], for 0-5% centrality in panel (a) and 50-60% centrality in panel (b). The comparison to the HKM model [30] is also shown for the central data. Closed symbols are experimental data (with statistical and systematic uncertainty), bands and dashed lines are the model calculations.

Conclusions
We reported on the centrality and pair k T dependence of the three-dimensional and direction-averaged one-dimensional pion femtoscopic radii in Pb-Pb collisions at √ s NN =2.76 TeV. The behavior of the femtoscopic radii can be factorized into a linear dependence on cube root of charged particle density (separately for each of the three-dimensional radii) and a power-law dependence on pair transverse momentum, with slightly different exponents for each direction. The dependence for the longitudinal radius is steeper than for the transverse ones. The dependence was also compared to the one observed in heavy-ion collisions at lower energies and to other collisions systems. The radii at the LHC follow the "universal" dN ch /dη 1/3 scaling in the long direction, but the radii in the transverse directions are below the "universal" curve. Simple linear scaling predictions are not valid when the collision energy is increased by an order of magnitude. The details of the dynamic evolution of the system influence the results significantly. This is in qualitative agreement with predictions from hydrodynamic models. In particular, when moving from RHIC to LHC collision energies, they produce a change in freeze-out shape, larger transverse radial flow and longer system evolution time. Comparison of the full dataset to the calculations from the recent hydrodynamic models, including three-dimensional evolution as well as hadronic stage, generally show a good agreement, which is complementary to similar agreement observed for momentum-only observables, such as momentum spectra and elliptic flow. The existence of such agreement both in the space-time as well as in momentum sectors, provides strong arguments for the validity of hydrodynamic models for the description of flowing bulk matter created in heavy-ion collisions at the LHC.