Two-particle Bose-Einstein correlations and their L\'evy parameters in PbPb collisions at $\sqrt{s_\mathrm{NN}}$ = 5.02 TeV

Two-particle Bose-Einstein momentum correlation functions are studied for charged-hadron pairs in lead-lead collisions at a center-of-mass energy per nucleon pair of $\sqrt{s_\mathrm{NN}}$ = 5.02 TeV. The data sample, containing 4.27 $\times$ $10^{9}$ minimum bias events corresponding to an integrated luminosity of 0.607 nb$^{-1}$, was collected by the CMS experiment in 2018. The experimental results are discussed in terms of a L\'evy-type source distribution. The parameters of this distribution are extracted as functions of particle pair average transverse mass and collision centrality. These parameters include the L\'evy index or shape parameter ($\alpha$), the L\'evy scale parameter ($R$), and the correlation strength parameter ($\lambda$). The source shape, characterized by $\alpha$, is found to be neither Cauchy nor Gaussian, implying the need for a full L\'evy analysis. Similarly to what was previously found for systems characterized by Gaussian source radii, a hydrodynamical scaling is observed for the L\'evy $R$ parameter. The $\lambda$ parameter is studied in terms of the core-halo model.


Introduction
The effect of quantum statistics is evident in the momentum correlations observed for identical particles emitted in high energy hadronic and nuclear collisions.These correlations were first studied using pion pairs produced in proton-antiproton collisions [1], where they were explained by the Bose-Einstein nature of pions [2].Subsequent to the first Bose-Einstein momentum correlation measurements, it was shown [3][4][5] that the particle correlations are closely related to earlier intensity correlation measurements used to determine the angular diameters of stars, as performed by R. Hanbury Brown and R.Q.Twiss [6], and are now referred to as HBT correlations.The observed momentum correlations can be related to the spatio-temporal particle emission probability, the so-called "source distribution", through a Fourier transform.This allows for a determination of the spatial structure of the particle-emitting source at the femtometer scale.
Measuring and interpreting such momentum correlations is called "femtoscopy" and these studies have significantly advanced our understanding of the processes governing heavy ion collisions [7].The fluid nature of the strongly coupled quark-gluon plasma created in heavy ion collisions was established, in part, by Bose-Einstein correlation measurements [8].In particular, the pair transverse momentum dependence of the Gaussian source radii [8][9][10] can be explained as a consequence of the hydrodynamical expansion of the medium [11,12].In addition, the measured source radii provide information on the transition from a quark-gluon plasma to a hadronic state of matter [13,14] and have also been important for extending our understanding of quantum chromodynamics to new regions of phase space [15].
In past femtoscopic measurements, either Gaussian [8,[16][17][18] or Cauchy [19,20] source distributions have generally been assumed.However, recent high precision correlation measurements in gold-gold (AuAu) collisions at a center-of-mass energy per nucleon pair of √ s NN = 200 GeV at BNL RHIC [21] and in beryllium-beryllium collisions at 150A GeV/c beam momentum at CERN SPS [22] have shown that neither the Gaussian nor the Cauchy approximation can adequately reproduce the measured results.Rather, a generalization of these distributions, the so-called Lévy alpha-stable distribution [23], is required to describe the data.The Lévy exponent α quantifies the deviation of the source shape from a Gaussian distribution and may be influenced by various physical processes [24][25][26][27][28].
In this paper, the source geometry is investigated for lead-lead (PbPb) collisions at √ s NN = 5.02 TeV by measuring quantum-statistical two-particle Bose-Einstein correlations of unidentified charged hadrons.The measured correlation functions are fitted with a formula based on a Lévy source distribution, allowing for the possibility of a non-Gaussian source shape.The extracted source parameters are studied as functions of pair transverse momentum and collision geometry.
The paper is structured as follows: The theory of Bose-Einstein correlations and the connection between the two-particle correlation function and the geometry of the particle emitting source are discussed in Section 2. The experimental setup is described in Section 3. The analysis details and systematic uncertainties are discussed in Section 4. The results are presented and discussed in Section 5, with a summary presented in Section 6. Tabulated results are provided in the HEPData record for this analysis [29].

Bose-Einstein correlations
The general definition of the two-particle momentum correlation function is given by where p 1 and p 2 are the particle four-momenta, N 1 (p 1 ) and N 2 (p 2 ) are the single particle momentum distributions, and N 12 (p 1 , p 2 ) is the two-particle momentum distribution.For identical bosons in high-multiplicity heavy ion collisions, the main source of correlations is the Bose-Einstein effect [21].This can be understood if one expresses the momentum distributions in terms of the source distribution, S(x, p), which is defined as the probability density distribution of creating a particle with four-momentum p at spacetime point x.As shown in Refs.[30,31], in the approximation that p 1 ≈ p 2 ≈ K, with K = (p 1 + p 2 )/2 being the average four-momentum, the relationship between the correlation function and the source distribution becomes where S(Q, K) is the Fourier transform of the source distribution, with Q = p 1 − p 2 denoting the relative four-momentum of the two particles.In Eq. ( 2), the superscript (0) denotes that final-state interactions of hadrons, such as the Coulomb interaction, are ignored.The effect of the Coulomb interaction will be discussed later in this section.
The above expressions imply C (0) This cannot be directly verified since the two-particle resolution limits the measurement of C 2 at small values of relative momentum Q.However, based on extrapolations of observed correlation functions, it is found that, in general, This observation has led to the development of a "core-halo" picture [32][33][34][35][36], where the particle emitting source is divided into two parts, a core of primordial hadrons (including also the decay products of short lived resonances) and a halo of long-lived resonances.The latter is experimentally unresolvable due to its large size that corresponds to very small momentum difference in Fourier space.Taking S to represent the core part of the source and λ as the square of the core fraction, one obtains ( Here N core is the number of hadrons in the core (for a given species) and N halo is the number of hadrons in the halo.The λ parameter is called the correlation strength, because it qualitatively describes how strong the two-particle correlation is.Other approaches to account for the effect of resonances are discussed in Refs.[34][35][36].
Using Eq. ( 4), the correlation function can be calculated for an assumed source distribution, with the parameters of the source distribution then determined by fits to the experimental results.The Gaussian, Cauchy, and a generalization of these two distribution functions, the so-called Lévy alpha-stable distribution, have been studied in Ref. [21].The motivation behind these studies is that the mean free path increases over time in an expanding hadron gas.Hence, if one describes the evolution of hadrons by generating random steps, the distribution of the size of these steps may not have a finite variance.In this case, the central limit theorem cannot be applied, thus the limiting distribution of the sum of step sizes for a given hadron is not a Gaussian.This effect has been referred to as "anomalous diffusion."However, a more generalized version of the central limit theorem is applicable and the limiting distribution in this case falls in the class of symmetric Lévy alpha-stable distributions [37,38], defined as with r being the variable of the distribution and q being an arbitrary integration variable.In this paper, L(r; α, R) is assumed for the shape of the source distribution, i.e., S mentioned above (and used in Eq. ( 3)) is assumed to be a Lévy distribution L, with parameters α and R depending on average momentum K.In this case, R is the Lévy scale parameter and α is the Lévy index of stability.Stability means that for α ≤ 2 this distribution represents the limit of properly normalized sums of independent random variables, and the sum of Lévy distributed random variables (with the same stability index) is also Lévy distributed.The Gaussian distribution is obtained with α = 2 and the Cauchy distribution is obtained with α = 1.A value of α > 2 is not physically allowed for a stable Lévy distribution (as discussed in Ref. [23]).The α parameter describes the shape of the source and its deviation from the Gaussian distribution.
In case of a Gaussian source, R is simply the quadratic mean of the distribution.However, for α < 2, the stable distributions do not have a second moment.In those cases, R is still proportional to the full width at half maximum (with the proportionality constant depending on α), and therefore represents the spatial scale of the source.The above considerations lead to the following correlation function for a general value of α [38]: Here q is a one-dimensional (1-D) variable, the magnitude of the spatial part of Q in a given frame: q = |Q|.The dependence on the average momentum K is carried by the parameters: λ, α, and R [21].
The above formulas are based on a plane wave analysis and are only valid for interactionfree particles.Final-state interactions distort this simple picture.In the so-called Bowler-Sinyukov method [39], the most significant final-state interaction in case of charged particles, the Coulomb interaction, is handled by including a correction term K C (q; R, α), with This correction can be calculated based on the Coulomb wave function, as determined by solving the two-particle Schr ödinger equation [21,40].Such a calculation requires numerical integration of the Lévy distribution multiplied by the square of the absolute value of the Coulomb wave function.A parametrization of the results of this calculation, as a function of R and α, is given in Ref. [40].The measured correlation functions are fitted with the formula described in Eq. (8).

The CMS detector
The CMS detector is built around a superconducting solenoid of 6 m internal diameter.Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromag-netic calorimeter (|η| < 3), and a brass and scintillator hadron calorimeter (|η| < 3), each composed of a barrel and two endcap sections, where η is the pseudorapidity.In addition to the barrel and endcap detectors, quartz-fiber Cherenkov hadron forward (HF) calorimeters (3 < |η| < 5) complement the coverage provided by the barrel and endcap detectors on either side of the interaction point.These HF calorimeters are azimuthally subdivided in ϕ into 20 • modular wedges and further segmented to form 0.175×0.175(∆η×∆ϕ) "towers", where ϕ is the azimuthal angle.A muon system located outside the solenoid and embedded in the steel flux-return yoke is used for the reconstruction and identification of muons with |η| < 2.4.
The silicon tracker measures the properties of charged particles with |η| < 2.5.During the 2018 LHC running period corresponding to the data used in this paper, the silicon tracker consisted of 1856 silicon pixel and 15 148 silicon strip detector modules.Details on the pixel detector can be found in Ref. [41].For nonisolated particles with a transverse momentum of 1 < p T < 10 GeV/c and |η| < 2.5, the track resolutions are typically 1.5% in p T and 20-75 µm in the transverse impact parameter (d xy ) [42].Events of interest are selected using a two-tiered trigger system [43].The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz.The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.A detailed description of the CMS detector, together with a definition of the coordinate system and kinematic variables, can be found in Ref. [44].
4 Measurement details

Event and track selection
The analysis presented in this paper used 4.27×10 9 minimum bias triggered events, corresponding to an integrated luminosity of 0.607 nb −1 [45,46], from PbPb collisions collected by the CMS experiment during the 2018 LHC run at √ s NN = 5.02 TeV.The minimum bias events were triggered by requiring signals above thresholds in the range of ∼6-12 GeV energy in at least one channel of each of the HF calorimeters [43].Further selections were applied to reject events from beam-gas interactions and nonhadronic collisions [47].Events were also required to have precisely one interaction vertex, reconstructed based on two or more tracks, and within a distance of less than 15 cm from the center of the nominal interaction point along the beam axis.The shapes of the clusters in the pixel detector had to be compatible with those expected from particles produced at the interaction vertex location.After these selections, 2.65×10 9 usable events were obtained.The event centrality was calculated from the transverse energy deposited in both HF calorimeters, using the methodology detailed in Ref [48].The centrality describes the fraction of the total inelastic hadronic cross section, and is connected to the degree of overlap of the two incident nuclei, with 0% corresponding to the most central (total overlap) and 100% to the least central (least overlap) collisions.
Only the standard CMS highPurity tracks defined in Ref. [49] were used for the analysis.Furthermore, individual particles were required to have p T > 0.5 GeV/c with a relative uncertainty δp T < 10%.To focus on mid-rapidity Bose-Einstein correlations, the analysis was restricted to tracks with |η| < 0.95.As the correlation function parameters depend on the average pseudorapidity of the considered pair [50], this limited pseudorapidity interval provides a cleaner data sample where the correlation function shape is not influenced by the η dependence.Tracks were required to have at least 2 hits (N pixel-hit ) in the silicon pixel detector, to have at least 11 hits (N hit ) in the strip tracking detector, and to have a relative χ 2 of the track fit (χ 2 N −1 dof N −1 layer ) of less than 0.18, where N dof is the number of degrees of freedom during the track fitting and N layer is the number of tracking layers used out of the 10 barrel and 4 pixel layers.A reconstructed track is only considered as a candidate track from the vertex if the significance of the separation along the beam axis (z) between the track and the best vertex, |d z /σ(d z )|, and the significance of the track impact parameter measured in a plane transverse to the beam, |d xy /σ(d xy )|, are each less than 3.
In CMS, particle identification in central to mid central PbPb collisions is hindered by the overlapping strip tracker clusters in a high track density environment.Therefore, in this measurement, no particle identification was done, and all charged-particle tracks were used.The approach presented in Section 2 is only applicable for a single particle species, hence it was assumed that all detected charged particles are pions.Depending on the transverse mass and the centrality range, 60-90% of the investigated particles are pions [51].It was found that this assumption does not modify the physical results, apart from the λ parameter.This is because the HBT peak in the correlation function at low q is not modified by nonidentical particle pairs or by identical particle pairs of species that have been misidentified as pions and, therefore, calculated to have the wrong q value.The effect of the lack of particle identification on the λ parameter is discussed in Section 5.3.

Event mixing and pair selection
For determining the correlation functions, two distributions are formed: a signal (actual) distribution A(q) of pairs containing the Bose-Einstein correlation, and a background distribution B(q), using the event mixing method to correct for correlations not stemming from quantum statistics or final state interactions.The signal pair distribution is formed with pairs of particles belonging to the same event.The q variable is taken as the absolute value of the three-dimensional (3-D) momentum difference (|q| LCMS ) in the longitudinally comoving frame (LCMS), which is the reference system where the longitudinal component of the average pair momentum is zero.This distribution is affected by various effects resulting from the event kinematics and the detector acceptance.To correct for these effects a background distribution is created with particle pairs taken from different events.For this purpose, a pool of events is formed based on centrality (within a width of 10%) and collision vertex location (within 3 cm in z-vertex).Then, for each data event, a mixed event was formed using the pool associated with the data event class, making sure that each particle of this mixed event belongs to a different data event.Pairs within the mixed events were used to create the above mentioned B(q) distribution.More details about the event mixing can be found in Refs.[21,52].The A(q) and B(q) distributions were divided and normalized to obtain the correlation function The integrals are performed over a range [q 1 , q 2 ] = [4.8,6.4] GeV/c where the correlation function is not expected to exhibit quantum-statistical features, which are only present for q ≲ 0.2 GeV/c [21,53], for source sizes of several femtometers.
Tracking efficiency correction factors were used as weights when measuring pair distributions.Each reconstructed track was weighted by the inverse of the efficiency factor, ϵ trk (η, p T , centrality).The efficiency weighting factor accounts for the reconstruction efficiency E r , and the fraction of misidentified tracks, with f fake (η, p T , centrality), and ϵ trk = E r /(1 − f fake ).
The obtained correlation functions exhibit effects stemming from the finite segmentation of the detectors and the imperfect tracking algorithm, allowing for tracks to split into two recon-structed particles, or separate tracks to merge into a single reconstructed particle.To remove these effects from our sample, a two-dimensional (2-D) histogram of pseudorapidity difference (∆η) vs. azimuthal angle difference (∆ϕ) was studied for the track pairs.A uniform distribution outside of a small ∆η and ∆ϕ region was found, with the detector artifacts limited to this small region.The pairs were then required to satisfy the condition with ∆η min = 0.014 and ∆ϕ min = 0.022, to select pairs free of these artifacts.
While in this analysis the correlation function is determined in terms of |q| LCMS , it can also be expressed in terms of the absolute value of the invariant momentum difference, q inv [21,54], which is the absolute value of the momentum difference in the pair center-of-mass system.Correlation functions depending on the transverse and longitudinal components of the 3-D momentum difference, as well as on the energy difference, were investigated.It was found that the correlation functions had their maximal values for small |q| LCMS , rather than small q inv values, indicating that the most sensitive 1-D variable of the measured correlations is |q| LCMS , hence the choice q ≡ |q| LCMS to investigate the physics parameters.
In the LCMS, the longitudinal component of the average momentum K vanishes.In this case, the Fourier transform of the source distribution S can be expressed in terms of the transverse component K T of the average momentum.Alternatively, the transverse mass, m T = √ m 2 + K 2 T , may be used, if the mass m of the investigated particle species is known.The current analysis is done in ranges of the K T variable, with the results expressed in terms of m T .Here, the pion mass is assumed, to facilitate comparisons with other experimental results and with theory.

Correlation functions
The A(q) and B(q) distributions were measured up to q = 8 GeV/c in 24 ranges of K T from 0.5 to 1.9 GeV/c and in six centrality classes (0-5%, 5-10%, 10-20%, 20-30%, 30-40%, 40-60%).Negative and positive same-charged hadron pairs were measured separately.In each case, the correlation function C 2 (q) was calculated.Although no significant difference was observed between the positively and the negatively charged pairs, the two cases were treated separately to identify minor discrepancies.While the Bose-Einstein peak for these correlation functions occurs for q < 200 MeV/c, a structure was also observed in the C 2 (q) distributions at higher q values.This higher q-structure can be attributed to a number of possible sources, including energy and momentum conservation, resonance decays, bulk flow phenomena [19], and minijets [19].We addressed this by fitting the following empirically determined functional form [19,53,55] to the long-range background: where N, α 1 , α 2 , R 1 , R 2 are fit parameters.An example fit to the experimental data using Eq. ( 11) is shown in Fig. 1.The correlation function C 2 (q) was then divided by the long range background function BG(q) resulting in the double-ratio correlation function The large-q background is close to constant on the scale of the Bose-Einstein peak, i.e., over a range of few times 10 MeV/c.Consequently, this correction is not as important as was found for proton-proton collisions [53], where the much smaller source radii lead to much wider correlation functions.BG(q) fit (q>0.2GeV/c) Figure 1: An example of a long range background fit to the correlation function C 2 (q) of negatively charged hadron pairs with 1.00 < K T < 1.05 GeV/c in the 10-20% centrality bin.The Bose-Einstein peak is below q = 0.2 GeV/c, therefore the 0.2 < q < 8.0 GeV/c region was used for the long range background fit.

L évy fits to correlation functions
The double-ratio correlation function DR(q) removes most effects other than the quantumstatistical correlation and the final state interactions.This function was fitted with an extension of Eq. ( 8) that allows for a possible residual linear background: Here R is the Lévy scale parameter, α is the Lévy stability index, and λ is the correlation strength.The empirical normalization parameter N and linear background parameter ϵ are needed to remove any remaining background effects in the double-ratio correlation function.
The fitting was performed using the MINUIT2 package [56,57] by standard χ 2 minimization.The asymmetric statistical uncertainties were calculated using the MINOS package [56,57].A fit is deemed statistically acceptable if the confidence level calculated from the χ 2 value and the number of degrees of freedom is above 0.1%.The convergence of each fit together with a positive definite covariance matrix was also required.Furthermore, the stability of the fit parameters versus the binning of DR(q) was tested, and no modification in the values of the parameters was found.Figure 2 shows the results of a typical fit.
Upper and lower limits in q were set for the fitted region.For large q values, the correlation no longer exhibits any features beyond those resulting from the residual long-range background.
For very small q values, the data are strongly affected by the finite momentum resolution and pair reconstruction efficiency of the detectors, as well as the details of the pair selection.These effects were investigated with detailed Monte Carlo simulations based on HYDJET (version 1.8, tune "Drum" [58]) events, where the detector response was simulated using GEANT4 [59], and a considerable decrease of pair reconstruction efficiency was found at low q values, below approximately 50 MeV/c, as shown in Fig. 3.Moreover, it has previously been shown in Refs.[40,60] that the sharp decrease of the correlation function for very low q values does not depend on the exact source shape, but is rather determined by the Gamow factor that corresponds to the spatial integral of the relative wave function squared multiplied by a point-like source [53].This very low-q behavior of the fit function is found to be unaffected when calculated with various source and final-state interaction assumptions [60].For the present analysis,  the lower and upper fit limits were selected independently for each centrality and K T class, with the lower limit ranging between 0.024 and 0.067 GeV/c and the upper limit ranging between 0.12 and 0.55 GeV/c.The limits were selected based on the confidence level and the convergence of the fit and based on the property of the Bose-Einstein peak, which moves towards higher q values for more peripheral collisions and as K T increases.

Systematic uncertainties
The sources of systematic uncertainties considered in the analysis include the criteria for the event selection, the single-track selection, the pair-track selection, and the limits set on the fitted q range.Furthermore, the centrality calibration is also varied to estimate the systematic uncertainty resulting from this source.For each K T and centrality class, the double-ratio correlation function DR(q) was obtained and fitted with the nominal analysis parameters.Variations in the fit parameters were then determined by changing, individually, each of the analysis parameters to both its lowest and highest plausible values, and then fitting DR(q) again using these values.The total systematic uncertainties were obtained by adding in quadrature all of the individual contributions for positive and negative modifications in the final fit parameters separately.Table 1 lists the analysis parameters and the range through which they were varied.These analysis parameters were separated into three categories.The "Track and event selection" category included the criteria on the z position of the vertex, the centrality edges, the selections of p T , δp T , |η|, N pixel-hit , |d xy /σ(d xy )|, and |d z /σ(d z )|, respectively.The "Pair selection" category included the (∆η, ∆ϕ) pair selection.The "Fit limits" category included the q min lower, and the q max upper fit limits.The asymmetric systematic uncertainties for each of the fit parameters are summarized in Table 2. Uncertainties for different K T ranges and the two charge signs are averaged in each centrality range.For every fit parameter, the dominant systematic uncertainty results from the modification of the fit limits.The uncertainty corresponding to the pair selection criteria is highly asymmetric for the three most central cases, it is zero in one of the directions for all three parameters.∆η min = 0.011 ∆ϕ min = 0.016 q min lower fit limit q 0 min (K T , cent) q 0 min − 0.004 q 0 min + 0.004 q max upper fit limit q 0 max (K T , cent) 0.85q 0

Results
Three physical parameters (R, α, λ) were determined for each centrality and K T class by fitting the theoretical formula to the measured correlation functions.The results are presented as a function of m T .Systematic uncertainties are separated into correlated and uncorrelated pointto-point parts.Correlated uncertainties associate the same relative uncertainty for a given parameter for all centrality and K T classes, whereas point-to-point values can vary between classes.This separation is necessary when investigating the m T or centrality dependence of the fit parameters, where only the uncorrelated point-to-point uncertainty is used in determining the χ 2 values for the fits.The correlated part of the systematic uncertainty for each parameter value (R, λ, or α) is taken as the average uncertainty for the parameter.The point-to-point uncertainty for a parameter value is then the difference between the full systematic and the correlated uncertainties.

The L évy scale parameter R
Figure 4 shows the obtained R values as a function of m T .All values are found between 1.6 and 5.8 fm.The decreasing trend in m T that is expected from hydrodynamics is observed.A clear centrality dependence is also visible.The decreasing value of R as collisions become more peripheral supports the geometric interpretation of this scale parameter.
To further analyze the m T dependence of R, 1/R 2 was plotted versus m T as shown in Fig. 5.This plot is motivated by hydrodynamic predictions [11,12] that suggest a linear dependence when a Gaussian source is assumed.Here it is investigated whether a similar linear behavior exists for a Lévy source.Linear fits are shown in Fig. 5, with the fit parameters tabulated in Table 3. Statistical and point-to-point systematic uncertainties were added in quadrature for the fits and used for the determination of the statistical uncertainties of the fit parameters.The confidence levels were statistically acceptable everywhere (i.e., have confidence level above 0.1%, similarly to Refs.[21,52]), suggesting that hydrodynamic scaling of 1/R 2 holds for a Lévy source as well.
The slope A and the intercept B of the linear fits are also related to hydrodynamic models.The slope A has connection to the Hubble constant (H) of the quark-gluon plasma [12,61], with where T f is the freeze-out temperature.With a Gaussian source assumption, the intercept B is related to the size of the source at freeze-out (R f ), with Figure 6 shows that the slope parameter A has a strong centrality dependence, since the average number of nucleons participating in the collision (⟨N part ⟩) is closely related to centrality.The mapping from a centrality class to ⟨N part ⟩ is given in Ref. [62].By assuming a constant T f ≈ 156 MeV [63], the value of the Hubble constant falls between 0.11-0.18c/fm, from the most central to the most peripheral collisions.This indicates that the centrality of the collision affects the velocity of the expansion, resulting in larger velocities in more peripheral collisions.The values of the Hubble constant obtained here are close to the value of 0.17 c/fm measured in high-multiplicity proton-proton collisions at √ s = 13 TeV [53], and are larger than the estimated value of 0.07 c/fm measured in 0-30% centrality AuAu collisions at √ s NN = 200 GeV [21].The intercept parameter B is found to be negative in all cases (as shown in Fig. 6), which prevents a determination of the freeze-out size using Eq. ( 15).In the hydrodynamic models assuming a Gaussian shape, B is positive.The negative value found here might be related to the Lévy source assumed in the present analysis.Hydrodynamic solutions where B is negative are possible [64], although the implications for a Lévy source are not clear.
Figure 7 shows the dependence of R on ⟨N part ⟩ 1/3 for eight representative samples of the 24 m T classes considered in this analysis, where the source volume is expected to scale with ⟨N part ⟩.The results are fitted with a linear function for each m T class.The fits are shown in Fig. 7, with the fit parameters and the corresponding confidence levels tabulated in Table 4.The results are again consistent with a geometrical interpretation of the Lévy scale parameter R.  Table 3: Fit parameters and the corresponding confidence levels (CL) of the linear fits to 1/R 2 versus m T , for positively (left) and negatively (right) charged hadron pairs.h    4. Table 4: Fit parameters and the corresponding confidence levels (CL) of the linear fits to R versus ⟨N part ⟩ 1/3 , for positively (left) and negatively (right) charged hadron pairs.h

The L évy stability index α
Figure 8 shows the values of α as a function of m T for both negatively and positively charged hadron pairs.Within uncertainties, all values are found to fall in the range of 1.6-2.0,with little m T dependence within a given centrality class.However, there is a clear centrality dependence, with the results tending to a Gaussian shape (α = 2) for the most central events.The observation that the correlation functions for some centralities can be described with Lévy functions having indices statistically inconsistent with a value of 2 means that a Gaussian assumption for the source shape is invalid in these cases.These Lévy distribution fits can be interpreted as suggesting the presence of anomalous diffusion resulting from expansion in the hadron gas.
To more clearly show the centrality dependence, as well as the significance of the deviation from 2.0, Fig. 9 shows ⟨α⟩ (the brackets indicating an average over m T ) as a function of ⟨N part ⟩ for both positively and negatively charged hadron pairs.A very similar linear dependence is observed for both charge signs.
The PHENIX Collaboration at RHIC reported a mean value for α of 1.207 for pions pairs with |η| < 0.35 and 228 < m T < 871 MeV/c 2 in 0-30% centrality AuAu collisions at √ s NN = 200 GeV [21].Very little m T dependence was observed at RHIC, a result similar to that reported here.However, our value for ⟨α⟩ in the same centrality range is 45-60% larger.This increase in the Lévy stability index may be connected to the larger energy densities achieved at the LHC, which would be expected to result in less anomalous diffusion.An increasing energy density could also explain the increase of ⟨α⟩ with ⟨N part ⟩.

The correlation strength λ
The measured λ values are shown in Fig. 10 as a function of m T , for both negatively and positively charged hadron pairs.The values decrease with increasing transverse mass and as the  The value of λ is reduced compared with the case of using identified pions only, because the sample contains particles other than pions (mostly kaons and protons).The strength of the correlation is proportional to the number of identical pion pairs, and pairs of nonidentical particles decrease the observed correlation.To see how the lack of particle identification affects the value of λ, a new parameter, λ * , is introduced by re-scaling λ with the square of the pion fraction: The centrality and K T dependent pion fraction was measured by the ALICE Collaboration in the range |η| < 0.5, as reported in Ref. [51].The calculated λ * values as a function of m T are shown in Fig. 10.In each of the centrality classes and m T bins, λ * is smaller than 1, which indicates that there exist additional effects lowering the correlation strength beyond the lack of particle identification discussed above.This result of λ * < 1 can, in fact, be understood on the basis of the core-halo model.In this interpretation, λ * describes the square of the fraction of pions produced in the core, as detailed in Eqs. ( 4) and (5).Furthermore, λ * is found to lack any clear m T dependence for a given centrality class.This is consistent with the observations at RHIC [21], and can be explained in terms of the core-halo model by the negligible momentum dependence of the core fraction in the investigated range.The decrease in λ * with centrality suggests a relatively smaller halo contribution for peripheral collisions.A geometric interpretation of the R parameter is suggested by its dependence on the average number of participating nucleons in the collision.Assuming a pion mass for the charged particles, a linear dependence of 1/R 2 on the transverse mass (m T ) is observed, consistent with a hydrodynamic scaling behavior, even for the case of Lévy sources.Based on the observed linear behavior, it is estimated that the Hubble constant of the quark-gluon plasma created in 5.02 TeV PbPb collisions increases from 0.11 c/fm to 0.18 c/fm when moving from most central to most peripheral collisions.The intercept of the 1/R 2 versus m T linear fits is negative in all cases, requiring further studies for its interpretation.The α parameter is found to have little, if any, m T dependence and to range between 1.6-2.0,increasing with centrality.The α values found in this paper are approximately 45-60% larger than those reported for gold-gold collisions at √ s NN = 200 GeV at RHIC.This increase, while not fully understood, may result from the greater energy densities achieved at the LHC.As a function of m T , a strong and decreasing trend is observed for the λ parameter, which can be explained by the lack of particle identification.After rescaling the λ values to account for the fraction of pions among the charged hadrons, a nearly constant trend of the pion-fraction-corrected λ * with m T is observed.The λ * values are found to be smaller than unity, which can be interpreted on the basis of the core-halo model as a nonnegligible halo contribution.Furthermore, λ * is found to decrease as the collisions become more central.Altogether, these results imply that the hadron emitting source in √ s NN = 5.02 TeV PbPb collisions can be described by Lévy distributions.This allows for a new and precise characterization of this source in high-energy heavy ion collisions.

Figure 2 :Figure 3 :
Figure 2: An example fit to the double-ratio correlation function DR(q) of negatively charged hadron pairs with 1.30 < K T < 1.35 GeV/c in the 20-30% centrality bin.The error bars show the statistical uncertainties.The fitted function is shown in blue, while the red overlay indicates the range used for the fit.The size of the Coulomb correction is indicated in magenta.The lower panel shows the deviation of the fit from the data in each bin in units of the standard deviation in that bin.

Figure 4 :
Figure 4: The Lévy scale parameter R versus the transverse mass m T in different centrality classes, for negatively (left) and positively (right) charged hadron pairs.The error bars show the statistical uncertainties, while the boxes indicate the point-to-point systematic uncertainties.These boxes are slightly shifted along the horizontal axes for better visibility.The correlated systematic uncertainty is also indicated.

Figure 5 :
Figure 5: The 1/R 2 distribution vs. transverse mass m T in different centrality classes, for negatively (left) and positively (right) charged hadron pairs.The error bars show the statistical uncertainties, while the boxes indicate the point-to-point systematic uncertainties.These boxes are slightly shifted along the horizontal axes for better visibility.The correlated systematic uncertainty is also indicated.A linear fit to the data is shown for each centrality bin.The fit parameters are tabulated in Table3.

Figure 6 :
Figure 6: The slope A (left) and the intercept B (right) linear-fit parameters versus ⟨N part ⟩, for positively and negatively charged hadron pairs.The error bars show the statistical uncertainties.The correlated systematic uncertainty is also indicated.The points are slightly shifted along the horizontal axes for better visibility.

Figure 7 :
Figure 7: The Lévy scale parameter R versus ⟨N part ⟩ 1/3 in different m T classes, for negatively (left) and positively (right) charged hadron pairs.The error bars show the statistical uncertainties, while the boxes indicate the point-to-point systematic uncertainties.The correlated systematic uncertainty is also indicated.A linear fit to the data is shown for each m T class.The fit parameters are tabulated in Table4.

Figure 8 :
Figure 8: The Lévy stability index α versus the transverse mass m T in different centrality classes, for negatively (left) and positively (right) charged hadron pairs.The error bars show the statistical uncertainties, while the boxes indicate the point-to-point systematic uncertainties.These boxes are slightly shifted along the horizontal axes for better visibility.The correlated systematic uncertainty is also indicated.

Figure 9 :
Figure 9: The average Lévy stability index ⟨α⟩ versus ⟨N part ⟩, for both positively and negatively charged hadron pairs.The error bars show the statistical uncertainties.The correlated systematic uncertainty is also indicated.The points are slightly shifted along the horizontal axes for better visibility.collisions become more central.

Figure 10 :
Figure 10: The correlation strength λ (upper panel) and λ * , which is re-scaled with the square of the pion fraction (lower panel), versus the transverse mass m T in different centrality classes, for negatively (left) and positively (right) charged hadron pairs.The error bars show the statistical uncertainties, while the boxes indicate the point-to-point systematic uncertainties.These boxes are slightly shifted along the horizontal axes for better visibility.The correlated systematic uncertainty is also indicated.integrated luminosity of 0.607 nb −1 , at a center-of-mass energy per nucleon pair of √ s NN = 5.02 TeV, recorded by the CMS experiment at the LHC.The correlation functions found in different centrality and average transverse momentum classes are analyzed in terms of Lévy sources, including Coulomb effects.The values of the Lévy scale parameter R, the Lévy stability index α, and the correlation strength λ are determined.
No. Z191100007219010; the Ministry of Education, Youth and Sports (MEYS) of the Czech Republic; the Hellenic Foundation for Research and Innovation (HFRI), Project Number 2288 (Greece); the Deutsche Forschungsgemeinschaft (DFG), under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306, and under project number 400140256 -GRK2497; the Hungarian Academy of Sciences, the New National Excellence Program -ÚNKP, the NK-FIH research grants K 124845, K 124850, K 128713, K 128786, K 129058, K 131991, K 133046, K 138136, K 143460, K 143477, 2020-2.2.1-ED-2021-00181, and TKP2021-NKTA-64 (Hungary); the Council of Science and Industrial Research, India; the Latvian Council of Science; the Ministry of Education and Science, project no.2022/WK/14, and the National Science Center, contracts Opus 2021/41/B/ST2/01369 and 2021/43/B/ST2/01552 (Poland); the Fundac ¸ão para a Ciência e a Tecnologia, grant CEECIND/01334/2018 (Portugal); the National Priorities Research Program by Qatar National Research Fund; MCIN/AEI/10.13039/501100011033,ERDF "a way of making Europe", and the Programa Estatal de Fomento de la Investigaci ón Científica y Técnica de Excelencia María de Maeztu, grant MDM-2017-0765 and Programa Severo Ochoa del Principado de Asturias (Spain); the Chulalongkorn Academic into Its 2nd Century Project Advancement Project, and the National Science, Research and Innovation Fund via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, grant B05F650021 (Thailand); the Kavli Foundation; the Nvidia Corporation; the Su-perMicro Corporation; the Welch Foundation, contract C-1845; and the Weston Havens Foundation (USA).

Table 1 :
The sources of the systematic uncertainties with their values in the default, lower, and upper settings.The meaning of the analysis parameters are given in Sections 4.1 and 4.2.

Table 2 :
The relative effect of the different types of systematic sources in each centrality class for the R (upper Table), α (middle Table)and λ (lower Table) parameters (in percentage).The values were averaged over K T and the two charge signs; the upwards (downwards) arrow, ↑ (↓) represents positive (negative) uncertainty in the value of the final fit parameters.