Two-pion femtoscopic correlations in Be+Be collisions at sNN=16.84\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\text {NN}}} = 16.84$$\end{document} GeV measured by the NA61/SHINE at CERN

This paper reports measurements of two-pion femtoscopic correlations in Be+Be collisions at a beam momentum of 150AGeV/c\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\,\hbox {GeV}\!/\!c$$\end{document} (energy available in the center-of-mass system for nucleon pair sNN=16.84\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\text {NN}}} = 16.84$$\end{document} GeV) by the NA61/SHINE experiment at the CERN SPS accelerator. The obtained momentum space correlation functions can be well described by a Lévy distributed source model. The transverse mass dependence of the Lévy source parameters is presented, and their possible theoretical interpretations are discussed. The results show that the Lévy exponent α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} is approximately constant as a function of mT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\text {T}}$$\end{document} , and far from both the Gaussian case of α=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha = 2$$\end{document} or the conjectured value at the critical endpoint, α=0.5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha = 0.5$$\end{document}. The radius scale parameter R shows a slight decrease in mT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\text {T}}$$\end{document} , which can be explained as a signature of transverse flow. Finally, an approximately constant trend of the intercept parameter λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} as a function of mT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\text {T}}$$\end{document} was observed, similar to previous NA44 S + Pb results (obtained with a Gaussian approximation, but unlike RHIC results).


Introduction
This paper reports measurements of quantum-statistical femtoscopic correlation functions for identified, like-sign charged pion pairs produced in central Be+Be collisions at 150A GeV/c beam momentum.
The method of quantum-statistical (Bose-Einstein) correlations was first applied in astrophysical intensity correlation measurements by R. Hanbury Brown and R. Q. Twiss (HBT) [1] in order to determine the apparent angular diameter of stellar objects.Later, a similar quantum-statistical method was applied in momentum correlation measurements for proton-antiproton collisions [2,3] to obtain information on the size parameters of particle emission sources in high-energy particle collisions.Since then, quantumstatistical (HBT) correlation measurements have become a standard tool for experimental characterization of the probability density in particle the emission process, i.e. the source function, which sheds light on the spatio-temporal structure of particle emission.This experimental method also largely contributed to the understanding of the hydro-dynamical nature of the produced strongly interacting matter.In fact, the pair momentum dependence of Gaussian shaped source radii [4,5] can be well explained by a hydro-dynamical expansion.The shape of the particle emitting source was furthermore suggested to be affected by the nature of the quark-hadron transition [6].Hence, exploring HBT correlations is of utmost importance in the quest for understanding the nature of the matter created in relativistic heavy-ion collisions.
The results presented in this paper were obtained by the NA61/SHINE [7] experiment at the CERN Super Proton Synchrotron accelerator.This measurement is part of the NA61/SHINE strong interaction program investigating the properties of the onset of deconfinement and searching for the possible existence of the critical point of strongly interacting matter.This goal is pursued by the NA61/SHINE Collaboration through a beam energy scan with various nucleus-nucleus collisions.This strategy allows to systematically investigate properties of the phase diagram of strongly interacting matter [8].
Within the framework of statistical models, the data on particle production suggest that with increasing collision energy, the temperature increases and the baryon chemical potential of strongly interacting matter decreases at freeze-out, whereas by increasing the nuclear mass number of the colliding nuclei, the temperature decreases [9,10].As a result of the NA61/SHINE research program, a large set of collision data on p+p, p+Pb, Be+Be, Ar+Sc, Xe+La, and Pb+Pb collision systems has already been recorded.An upgrade of the NA61/SHINE experiment was completed in 2022, and further high-statistics data on Pb+Pb collisions will be collected in the near future [11].A basic reference has already been established with p+p, Be+Be, and Ar+Sc interactions on particle spectra and multiplicities [12][13][14][15][16][17].The present paper provides results on Bose-Einstein correlations of identified, like-sign pion pairs in 0-20% centrality selected 7 Be+ 9 Be collisions.The data were recorded in 2011, 2012, and 2013 using a secondary 7 Be beam produced by fragmentation of the primary Pb beam from the CERN SPS [18].
The paper is organized as follows.Section 2 recalls the fundamental theory behind the technique of Bose-Einstein correlations in order to fix notations.In Section 3, the NA61/SHINE detector is described.In Section 4, the details of the analysis procedure are discussed.In Section 5, the obtained experimental results are presented.The paper closes with Section 6, summarizing the results and conclusions.

Bose-Einstein correlations 2.1. Bose-Einstein correlation functions
The two-particle Bose-Einstein correlations are defined in terms of the single-and the two-particle invariant momentum distributions N 1 and N 2 as: where p 1 and p 2 are the momenta of the individual particles.If the S (x, p) (source function) denotes the probability density of particle creation at space-time point x and momentum p, the momentum distribution of emitted particles can be expressed via this source function as [19]: where Ψ p (x) and Ψ p 1 ,p 2 (x 1 , x 2 ) are the single-and two-particle wave functions.In the case of the singleparticle wave function, |Ψ p (x)| 2 = 1 holds, whereas for the two-particle wave function, taking into account the Bose-Einstein symmetrisation, one has [20]: where X = x 1 − x 2 is the relative coordinate and Q = p 1 − p 2 is the relative momentum of the pair.In the QX term, a division by ℏ is suppressed, and throughout this paper, we will utilize the ℏ = 1 convention.Substituting the above equation into Eq.( 1), one infers where S is the Fourier transform of S in its first variable.If relative momentum Q is much smaller compared to the average momentum of the pair K = (p 1 + p 2 )/2, then Eq. ( 5) can be expressed as: Alternatively, with a simplified notation where the K-dependence is suppressed and a normalized source is assumed, one may write This choice is motivated by the so-called smoothness approximation [21].The dependence on relative momentum Q is stronger than on the average momentum of the pair K, hence Q is considered as the more important variable of the correlation function, and the other variable is mostly suppressed in the notation.When the correlation function is parameterized based on an ansatz for the source function, its parameters can depend on K.In order to explore the transverse dynamics of the source, the average transverse momentum of the pair, where K x and K y are transverse components of K. Furthermore, motivated by hydro-dynamical considerations, the dependence on the transverse mass m T = m 2 c 4 + K 2 T c 2 is often studied, where m is the particle mass.The m T -dependence of the source parameters, such as its width, the so-called HBT scale, or radius R, was crucial in understanding the transverse expansion dynamics of the strongly interacting matter [22,23].One of the main goals of HBT correlation measurements is to estimate the size (or, rather, the correlation length) of the hadron emitting source.

Core-halo model
Equation (6) implies that the correlation function takes the value of 2 at zero relative momentum, or equivalently, if the notation C 2 (Q → 0) = 1 + λ were used, then λ = 1 would follow.However, the intercept parameter λ is often smaller than one in experimentally measured correlation functions.The widely accepted explanation for this phenomenon is the core-halo model [24,25], namely that some correlated particles are produced in decays of long-lived resonances, creating a spatially extended component of the source, their momentum difference being unresolvable by the detector.The core-halo model treats these as belonging to the halo component of the source, while the primordial particles and the decay products of short-lived resonances represent the core.While the latter has a size of a few fm, the former may extend to thousands of femtometers, due to long-lived resonances.One can then break up the source S into S core and to S halo as follows: In experimental measurements, the accessible range of Q values is not smaller than a few MeV/c, due to the finite two-track resolution of the tracking detectors.Because of the large radius of the halo, in the accessible Q-range S halo (Q, K) ≈ 0 thus, S (Q, K) ≈ S core (Q, K).Given that the Fourier-transform of each of the source components at Q = 0 equals to the number (N) of particles in that component, follows, and therefore one obtains with for the experimentally resolvable Q-range.Although the core-halo model provides a natural explanation for the phenomenon C 2 (Q → 0) < 2, i.e. λ < 1 in experimental data, it is important to note that λ 1 can also be explained by other effects, such as coherent pion production [20,26] or background from improperly reconstructed particles.It is evident, however, that measuring λ is an important tool in understanding particle creation in relativistic heavy-ion collisions.Other effects, such as Coulomb ones, are discussed more in detail in Sec.2.4 and strong interactions are negligible [27].

Lévy shaped sources and the QCD critical endpoint
When the source size (i.e., the HBT scale parameter R) or the correlation strength (i.e., the intercept parameter λ) of the Bose-Einstein correlation is to be measured, a full three-dimensional source reconstruction can be performed if the available statistics allow it.While in our analysis it cannot be proven experimentally that the source is fully spherical, studies of 2D correlation functions in the transverse and longitudinal momentum difference variable show no need to go beyond the 1D approximation.Nevertheless, for a more complete study with increased statistics, utilization of spherical harmonics decomposition [28] or performing further studies in three dimensions may be adequate.Alternatively, a parametric ansatz for the source shape may be used, and its derived correlation function is fitted to the data in order to determine its shape parameters.Quite naturally, Gaussian sources lead to Gaussian correlation functions.
In the present analysis, a more general ansatz is used, i.e. that of Lévy shaped sources [29,30], exhibiting possible power-law tails and also incorporating the Gaussian limit.Correlation functions based on this ansatz have been shown to describe LEP [31], RHIC [32], and LHC [33,34] data as well.
The spherically symmetric Lévy distribution is defined as where the parameters of this distribution are α and R, the Lévy stability index and Lévy scale parameter, respectively, while r is the vector of spatial coordinates and the vector ζ represents the integration variable.
In the case of α = 2, one recovers the Gaussian distribution, while α = 1 is equivalent to the Cauchy distribution, and for α < 2, the Lévy distribution exhibits a power-law tail.Hence determining the parameter α by a fit to experimental data yields a way to estimate the deviation of the source from a Gaussian or a Cauchy shape.
Ideally, the correlation function C 2 is investigated as a function of momentum difference in the entire three dimensions, but in case of statistically insufficient data samples, or spherically symmetric sources it is advantageous to measure the correlation functions versus a single-dimensional momentum variable.
A natural choice may be the invariant momentum difference, equivalent to the magnitude of the threemomentum difference in the pair comoving (i.e.pair-center-of-mass) system (PCMS).Another possible choice is the magnitude of the three-momentum difference in the longitudinally comoving system (LCMS): where the coordinate system is set up such that z is the direction of the beam, also sometimes called the longitudinal direction; and the transverse plane coordinates are x and y, which can be chosen arbitrarily.The momentum difference in this direction can be expressed in the LCMS as: where E 1 and E 2 are the energies of the respective particles.The LCMS can be advantageous since hadron emission turns out to be approximately spherically symmetric in this frame at RHIC energies [32] and at GSI, HADES as well [35].We note in passing that preliminary investigation of the full threedimensional correlation function indicate that, indeed this is a natural variable for parametrizing Bose-Einstein correlations, also at SPS energies for Be+Be collisions.
Assuming a three-dimensional spherically symmetric Lévy shaped source function and the core-halo model, the corresponding parametric form of the two-particle Bose-Einstein correlation function becomes Its three parameters λ, R, and α implicitly depend on the average transverse momentum K T , or, alternatively, on the transverse mass m T .
The shape parameter α carries information on the nature of the quark-hadron transition.Namely, lattice QCD calculations [36][37][38] and other theoretical expectations show two important regions of the baryochemical potential (µ B ) axis of the phase diagram of strongly interacting matter.A phase transition of analytical or "cross-over" type is expected at low µ B values and a first-order phase transition at high values of µ B .Therefore, a critical endpoint of the phase transition line is expected, where a second-order phase transition takes place.The mapping of the phase diagram of strongly interacting matter, and determining the position of the critical endpoint is one of the main goals of high-energy heavy-ion physics experiments, such as CBM [39], HADES [40], PHENIX [41], STAR [42], and NA60+ [43].
The Lévy stability index α is related to the spatial critical exponent η [44], since at the critical endpoint, fluctuations appear at all scales and spatial correlations will exhibit a power-law tail of the form ∼ r −1−η , and Lévy distributed sources also exhibit a power-law tail ∼ r −1−α (in three dimensions).
Theoretical expectations suggest that the universality class of QCD to be the same as that of the 3D Ising model [45,46].The value of the exponent η around the critical point in the 3D Ising model is 0.03631 ± 0.00003 [47], and with a random external field it is seen to be 0.50 ± 0.05 [48].This argument suggests that close to the critical endpoint (CEP) of the phase transition line of strongly interacting matter, α should also decrease to values near or even possibly below 0.5.While finite size effects and dynamics may modify this simple picture, measuring the Lévy stability index α is still expected to provide a signature of the critical point of the phase diagram of strongly interacting matter.

Final state Coulomb effect
The final state phenomena, such as the electromagnetic interactions between charged hadrons, were neglected in the above considerations.Namely, the quantum-statistical correlation functions discussed so far were obtained with the plane-wave assumption for the wave function.In the following, these will be denoted by C 0 (q).Let's suppose that the final-state electromagnetic interactions are included in the correlation function.Then the correlation function has to be calculated not via the interference of plane waves, but rather via the interference of solutions of the two-particle Schrödinger equation having a Coulomb-potential, describing the final state electromagnetic interactions.The ratio of these two correlation functions is called the Coulomb correction [27,49]: The numerator in Eq. ( 15) cannot be calculated analytically and is quite tedious to estimate numerically.
To simplify experimental analysis, in Ref. [33], an approximate formula was obtained and utilized subsequently for the case of Cauchy-shaped sources (α = 1).However, the Coulomb correction may also depend on the Lévy stability index α, hence a more precise treatment is required.To this end, a numerical calculation was performed in Refs.[49,50] and the results were parameterized, then the dependence on R, λ and α has been parameterized as well.In this analysis, we utilize the results obtained in Refs.[49,50] for estimating the Coulomb effect.
In order to take into account the effect of the halo mentioned in Sec.2.2, the Bowler-Sinyukov method [51,52] is utilized.The halo part only contributes at very small values of relative momenta, and hence it does not affect the source radii of the core component [53].This justifies the mentioned Bowler-Sinyukov method, in which the fit ansatz function is as follows: Here K Coulomb (q) is the Coulomb correction, and a normalisation parameter N was also introduced.
In addition, another effect has to be taken care of, which is related to the fact that the Coulomb correction is calculated in the PCMS while the measurement is done in the LCMS.When measuring one-dimensional HBT correlations in the LCMS, the assumption is that the source is of spherical shape, meaning However, the source is only spherical in the LCMS, hence an approximate one dimensional PCMS size parameter needs to be estimated.
This was done in Ref. [54], where an average PCMS radius of was obtained, where β T = K T m T .Furthermore, the following fact has to be taken into account: the momentum difference in the Coulomb correction is expressed in the PCMS, as a function of q PCMS = q inv (the invariant four-momentum difference equals the three-momentum difference in the PCMS).Since the reconstruction of q inv for a given pair (knowing only q) is not possible, the measurement should be performed as a function of both q as well as q inv .The estimation performed in Ref. [54] showed that a simple, approximate relation of the two may be given as q inv ≈ 1 − β 2 T /3 • q.Implementing both of the above mentioned effects results in the following formula for the Coulomb correction expressed in terms of q and R LCMS , based on the 3D calculation in PCMS: where now the dependence on R is indicated explicitly.This modified Coulomb correction is then used in Eq. ( 16).Note that K Coulomb depends also on the Lévy-index α, but the mentioned PCMS-LCMS transformation leaves this parameter unaffected, hence we suppressed this from the function arguments above.Furthermore, it should also be underlined that the effect of modifying the Coulomb correction based on the PCMS-LCMS difference discussed above is small, particularly negligible, compared to the listed systematic uncertainty sources of Section 4.6.

The NA61/SHINE detector
The NA61/SHINE fixed-target experiment uses a large acceptance hadron spectrometer located in the North Area H2 beam line of the CERN Super Proton Synchrotron accelerator [7].The main goals of the experiment include the investigation of the phase diagram of strongly interacting matter.A schematic of the layout employed during the Be+Be data taking is shown in Fig. 1.

Detectors
The key components of the experiment for the detection of particles produced in the collisions are the five large-volume Time Projection Chambers (TPCs) for tracking.The two most upstream ones are the Vertex TPCs (VTPCs), residing in the two superconducting bending magnets.The magnets have 9 T•m maximum combined bending power.Downstream of the VTPCs, the two Main TPCs (MTPCs) are located symmetrically to the beam line in order to extend the tracking lever arm and to perform particle identification by measuring their ionisation energy loss in the TPC gas.One smaller TPC is located in the gap between VTPCs, and is called Gap-TPC (GTPC; denoted GAP TPC in Fig. 1).The VTPCs and GTPC are operated with an Ar(90):CO 2 (10) gas mixture and the MTPCs with an Ar(95):CO 2 (5) mixture.The further downstream Time-of-Flight (ToF) detectors are not used in the present analysis.
The Projectile Spectator Detector (PSD) at the end of the setup is a segmented forward hadron calorimeter, centered on the nominal deflected beam trajectory.It measures the energy contained in the projectile spectators which is used for event centrality characterization.Central collisions are selected by lower values of this very forward energy.Although for smaller systems, such as Be+Be collisions, the very forward energy is not expected to be tightly correlated with the actual impact parameter of the collision, the terms central and centrality are still adopted following the convention widely used in heavy-ion physics.
The beam line instrumentation is schematically shown in the inset of Fig. 1.A set of scintillation counters as well as beam position detectors (BPDs) [7] upstream of the target provide timing reference, selection, identification and precise measurement of the position and direction of individual beam particles.

Triggers
The schematic of the placement of the beam and trigger detectors is shown in the inset of Fig. 1.These consist of a scintillation counter (S1) recording the presence of the beam particle, a set of veto scintillation counters with hole (V0, V1, V1 p ) used for rejecting beam halo particles, and a threshold Cherenkov charge detector (Z).Trigger signals indicating the passage of valid beam particles are defined by the coincidence T1 = S1 • V0 • V1 • V1 p • Z(Be) for high momentum data taking.
Central collisions were selected through the analysis of the signal from the 16 central modules of the PSD [55].The low-energy part of the deposited energy spectrum was selected to contain 20% of the most central collisions.The interaction trigger condition was thus T2 = T1•PSD for the higher energies.
The data consists of ≈ 2.828 • 10 6 events before event and track selection.

The 7 Be beam and 9 Be target
The NA61/SHINE beam line is designed to handle primary as well as secondary beams.The beam instrumentation was optimized accordingly.In the Be+Be runs, a secondary beam was used, fragmented from a primary Pb beam from the SPS accelerator.A threshold Cherenkov charge tagging detector, called the Z detector, was used in order to identify and select the Z = 4 fragment nuclei.In order to have a low material budget for the Z detector, a thin quartz wafer Cherenkov radiator was used.Additionally, the amplitudes of the signals measured in the three Beam Position Detectors (BPDs, see Fig. 1) were used to improve the Z resolution.A detailed description of the technique for the identification of 7 Be fragments is given in Ref. [18].
The target was a 12 mm thick plate of 9 Be placed approximately 80.0 cm upstream of VTPC-1.The total mass concentrations of impurities in the target were measured at 0.287% [56].No correction was applied for this negligible contamination.

Analysis procedure
Femtoscopic correlations are studied in this paper for pions reconstructed as originating from the primary interaction in the 20% most central 7 Be+ 9 Be collisions selected by the total energy emitted into the forward direction covered by the PSD detector.In the following, we describe the event, track and pair selection procedure, and all the steps required to obtain the measured source parameters.

Event selection
The events considered for analysis had to satisfy the following conditions: (i) there are no off-time beam particles detected within a time window of 4.5 µs around the particle triggering the event (this is the time needed to have good position resolution of the main vertex at all times), (ii) the event has a well-fitted main interaction vertex, (iii) the maximal distance between the main vertex z position and the centre of the beryllium target is between ± 5 cm (vertex z), (iv) the 0-20% most central collisions, based on PSD energy measurement, are accepted.

Track selection
The tracks selected for the analysis had to satisfy the following conditions: (i) the fit of the particle track converged, (ii) the distance between the track extrapolated to the interaction plane and the interaction point (impact parameter) should be smaller or equal to 4 cm in the horizontal (bending -|B x |) plane and 2 cm in the vertical (drift -|B y |) plane1 , (iii) the total number of reconstructed points in all TPCs on the track should be at least 30 (nPoint) and, at the same time, the sum of the number of reconstructed points in VTPC-1 and VTPC-2 should be at least 15 (VTPC points) or the number of reconstructed points in the GTPC should be at least 5 (GTPC points), (iv) the ratio of the total number of reconstructed points on the track to the potential number of points should be between 0.5 and 1.02 (nPointRatio), (v) identified particle's rapidity is in the interval ±2 around midrapidity.

Particle identification
Particle identification (PID) for pions was performed using dE/dx cuts, shown in Fig. 2. Tracks within the two lines are considered pions.The fine-tuning of the dE/dx cut parameters was performed as follows.
(i) a reasonable interval in ln(p/1(GeV/c) was selected where the pion contribution dominates (from -1.6 to 4, 0.2 to 55 GeV/c), (ii) the data was binned in ln(p) into 80 slices, (iii) in each bin, a Gaussian was fitted to the dE/dx data, in order to establish the most probable value of the pion dE/dx peak, (iv) the standard deviation (σ) obtained from these Gaussian fits was used to determine the pion dE/dx response width (found to be between 0.05 and 0.18, depending on momentum), which was the basis of the pion dE/dx selection as shown in Fig. 2 (PID cut).

Pair selection and event mixing
Due to the possible imperfections of the detector and of the tracking algorithm, hits created by a single particle may be reconstructed as two tracks.This is called track splitting and leads to a track pair with small momentum difference (< 20 MeV/c).Furthermore, the hits of two close particles may be reconstructed as a single track: this is called track merging.In a correlation analysis, it is important to minimize the effect of these track reconstruction problems.Track splitting is already largely removed by the track selection cuts (in particular, (iv) of Sec.4.2).The contribution from track merging was estimated by Monte Carlo (MC) simulations using EPOS simulation [57] and GEANT3 for particle propagation [58] and reconstruction, and an appropriate lower limit in the momentum difference was defined, as described below.
The basic quantity of correlation measurements is the pair distribution.From pairs of pions created in the same event, one obtains the so-called actual pair distribution A(q).Such distributions were measured in several intervals of average pair transverse momentum K T or pair transverse mass m T .This pair distribution is influenced by single-particle momentum distributions, kinematic acceptance of the detector, phase-space effects, and other phenomena not connected to quantum-statistics or final state interactions.These can be removed by constructing a combinatorial background pair distribution B(q), measured in the same K T or m T intervals as the A(q) distribution.Calculating this background distribution starts with the event mixing procedure, where an artificial, mixed event is created from particles originating from different events.Subsequently, pairs formed within the mixed event are used to create the background distribution B(q).By construction, this method ensures that no two particles are selected from the same background event, creating an uncorrelated pool of events.
The obtained background distribution B(q) exhibits all the previously mentioned non-quantum-statistical effects (acceptance, momentum distribution, phase-space, etc.), hence dividing A(q) by B(q) leaves us with a ratio which exhibits quantum-statistical and final-state interaction effects as well as the effect of reconstruction inefficiencies (and, in addition, momentum conservation, which is not relevant in the range of the investigated q-range).Thus, the measured correlation function is defined as where [q 1 , q 2 ] is a large-q range where quantum-statistical effects no longer affect the correlation function.The integrals in Eq. ( 19) provide the normalization of the correlation function to unity at high relative momentum.An example for C 2 (q) is shown in Fig. 3 for both data and EPOS simulation.It is  readily apparent that at low q values Bose-Einstein correlation and Coulomb repulsion effects determine C 2 (q) data points.These effects are not present in the simulations, hence the simulated C 2 (q) values are approximately constant.The reconstructed C 2 (q), however, suffers from track merging effects (where the two tracks forming the pairs are close spatially), strongly suppressing C 2 (q) at low q values.Hence the deviation of the simulated and reconstructed correlation function provides a good estimate of the range where inefficiencies are important.This allows to determine the range in q over which fits can be considered reliable.The fit range is then selected for each K T bin, e.g. for K T = 0.20 − 0.35 GeV/c the interval where fit is considered good is q = 0.049 − 0.8 GeV/c.This method then ensures that track merging is not present in the fitting interval.

Estimation of source shape parameters via fitting
The measured correlations were fitted with the formula described in Eq. ( 16).Due to the often modest number of entries in the signal and background distributions [59], Poisson maximum-likelihood fitting was used [60].The corresponding penalty function (χ 2 λ,p ) to minimise is where λ and p denote the fact that we are using a likelihood χ 2 for Poisson distributed histograms, c i references the number of counts, n i is the number of entries in the i th histograming bin obtained from the data, and y i is its corresponding parametric model value to be fitted to the data.Goodness-of-fit was determined using regular χ 2 methods in two ranges: the full range and the Bose-Einstein peak range.Fits were done both for positively and negatively charged pion pairs, as well as their combinations, in four m T intervals.A fit was accepted if the algorithm converged, the covariance matrix was positive definite, and the confidence value corresponding to the χ 2 and NDF was larger than 0.1%.An example fit is shown in Fig. 4. To estimate the statistical uncertainties of the fit parameters, the Minos method was utilised [61] (also called as likelihood based confidence intervals) which by its nature, yields asymmetric statistical errors.16) within the range of 0.049 GeV/c to 0.8 GeV/c, and the black dotted line indicates the extrapolated function outside of the fit range.

Systematic uncertainties
In the analysis, one has to consider that the parameters obtained from fits depend on several experimental choices and cuts, such as the PID cut, the width of bins, or the fitting range.These dependencies are the dominating contributors to the systematic uncertainties.In order to estimate these, the fits were performed with the loose and tight event and track selection criteria, and also with slightly varied fit intervals.The standard set of cut values together with the alternative values for systematic error estimation are shown in Table 1.The systematic uncertainty calculation was performed for positively and negatively, like-sign charged pairs summed together.
The combined systematic uncertainties were obtained as follows.Let P denote the fit parameter vector (α, λ, R).Denote by P j n (i) the corresponding estimated parameter vector obtained for the i-th m T bin (i = 0, . . ., 3), with the n-th selection criterion (n = 0, . . ., 8) listed in Table 1 set to the j-th setting ( j = 0, 1, 2 meaning the standard, tight and loose values).The downward (δP − ) and upward (δP + ) systematic uncertainty of the parameter vector P was estimated as follows: Table 1: The standard cuts used to obtain the final results, as well as the loose and tight cuts applied for estimation of systematic uncertainties. δP where P 0 (i) is the parameter vector in i-th m T bin with standard cut ( j = 0), J + n and J − n are the array of j values with which P j n (i) > P 0 (i), and P j n (i) < P 0 (i) occurs respectively, and N j+ n and N j− n denote their corresponding multiplicity.

Results
The three physical parameters (α, λ and R) were measured in four bins of pair transverse momentum K T or pair transverse mass m T .The parameters were obtained via fitting a parametric Lévy ansatz on the source via the formula Eq. ( 16) to the measured correlation functions.
The transverse mass dependence of the intercept parameter λ is shown in Fig. 5.One may observe, within uncertainties, λ(m T ) ≈ const. in the available m T range.An interesting phenomenon becomes apparent when compared to measurements at RHIC energy Au+Au collisions [32,62,63] and at SPS energy S+Pb and Pb+Pb collisions [64,65].At the SPS energies, there is no visible decrease of λ at lower m T values, but at RHIC energies, a "hole" appears at m T values around 2-300 MeV.This "hole" was interpreted in Refs.[32,66] to be a sign of in-medium mass modification of the η ′ meson.The NA61/SHINE results do not indicate the presence of a low-m T hole.Furthermore, it is important to note that the obtained values for λ are smaller than unity, which in the framework of the core-halo model may indicate that a significant fraction of pions are the decay products of long-lived resonances.
The measured values of the radial scale parameter R of the Lévy shaped source function, determining the homogeneity length of the pion emitting source in the LCMS, are shown in Fig. 6  m T .Interestingly, the resulting R parameter values are similar to those measured in p+p collisions at the CMS [33,67].We also observe a slight decrease of R with increasing m T .This can be explained by the presence of radial flow, based on simple hydrodynamical models [68,69] where one obtains a 1/R 2 ∝ m T type of transverse mass dependence: This function was fitted to the R values measured in each m T bin, as shown in Fig. 6, resulting in a good fit quality (χ 2 /NDF = 1.7/2, corresponding to a confidence level of 44%).The obtained fit parameters are: A = 4.5 ± 2.9 (stat.)fm and B = 0.12 ± 0.23 (stat.)GeV, comparable to those measured in p+p collisions at CMS [67] (although the large uncertainties prevent a quantitative comparison).
The Lévy stability exponent α describes the shape of the tail of the source distribution.The NA61/SHINE results, shown in Fig. 7, yield values for α between 0.9 and 1.5, and are significantly lower than the Gaussian (α = 2) case, and also significantly higher than the conjectured critical endpoint value (α = 0.5).The obtained α values are in a similar range as the ones obtained in Au+Au collisions at RHIC energies [32].The shape of the pion emitting source is apparently independent of m T , within uncertainties.Therefore one can calculate a simple average of the four α values via a constant fit, shown in Fig. 7, and resulting in a good fit quality (χ 2 /NDF = 6.0/2, corresponding to a confidence level of 11%).This results in an average value of α = 1.07 ± 0.06 (stat.), which describes a source shape close to a Cauchy distribution (where α = 1).Further studies are foreseen at NA61/SHINE using different collision energies and system sizes in order to map the evolution of the Lévy stability index α as a function of collision energy and system size.Furthermore α ≪ 2 is interesting in particular as we observe a R ∼ 1/ √ m T , visible in Fig. 6.It is not entirely clear why it is the case (the indicated m T dependence could form in the QGP or at a later stage), as one would expect this to show only at α = 2 [70]; this phenomenon was observed also at RHIC [32].

Summary and conclusions
Measurement of two-particle femtoscopic correlations in 150A GeV/c Be+Be collisions with the NA61/SHINE detector system was presented.The correlation functions were measured in several bins of pair transverse mass m T , and their fundamental shape parameters were extracted via fitting a Lévy shaped ansatz for the particle source function.Correction for the final state Coulomb interaction was performed.The m T -dependence of the shape parameters λ, R and α were studied.
The results show that the Lévy exponent α is approximately constant as a function of m T , and far from both the Gaussian case of α = 2 or the conjectured value at the critical endpoint, α = 0.5.The radius scale parameter R shows a slight decrease in m T , which can be explained as a signature of transverse flow.Finally, an approximately constant trend of the intercept parameter λ as a function of m T was observed, clearly different from measurement results at RHIC, but similar to previous NA44 measurements based on a Gaussian approximation.The NA61/SHINE experimental program plans further measurements at  different energies and system sizes of these Lévy shape parameters.This will complete a systematic study of the energy and system size dependence of the source shape parameters.

B. Systematic uncertainties
Tables 2-4 contain the physical parameter values from Lévy fits and their systematic uncertainty contributions from various sources.Values in the table indicate the change of the given fit parameter when the indicated analysis setting was modified, as compared to its default value.In these tables, an uncertainty value of 0.00% means that the given setting modified the fit parameter value only in the opposite direction.

Figure 1 :
Figure 1: Schematic of the NA61/SHINE detector setup, used during the Be+Be data taking.

Figure 3 :
Figure3: The ratio of A(q) and B(q) as a function of q.Left: measured NA61/SHINE data.Right: EPOS simulation and EPOS reconstructed data (GEANT3 + NA61/SHINE rec.chain).

Figure 4 :
Figure 4: Example fit with Bose-Einstein correlation function at K T = 0.20 − 0.35 GeV/c for the sum (π+ + π + ) + (π − + π − ).Blue points with error bars represent the data, the green dash-dotted line shows the fitted function with Coulomb correction given by Eq. (16) within the range of 0.049 GeV/c to 0.8 GeV/c, and the black dotted line indicates the extrapolated function outside of the fit range.

Figure 5 :
Figure 5: The intercept parameter λ, for 0-20% central Be+Be at 150A GeV/c, as a function of m T .Boxes denote systematic uncertainties, bars represent statistical ones.

Figure 6 :
Figure 6: The radial scale parameter R, for 0-20% central Be+Be at 150A GeV/c, as a function of m T .The fit to R(m T ) with Eq. (23) is shown with a solid line.Boxes denote systematic uncertainties, bars represent statistical ones.

Figure 7 :
Figure 7: The Lévy stability index α, for 0-20% central Be+Be at 150A GeV/c, as a function of m T .Special cases corresponding to a Gaussian (α = 2) or a Cauchy (α = 1) source are shown, as well as α = 0.5, the conjectured value corresponding to the critical endpoint, while the constant α fit is shown with a solid line.Boxes denote systematic uncertainties, bars represent statistical ones.

Figures 8 -
Figures 8-10 show the individual correlation function measurements and the corresponding fits in each of the K T intervals, separately.

Figure 8 :
Figure 8: Measured Bose-Einstein correlation function and the corresponding fit, similarly to Fig. 4, but for K T = 0.00 − 0.10 GeV/c.Fit range in q in this case was 0.028 GeV/c to 0.8 GeV/c.

Figure 9 :
Figure 9: Measured Bose-Einstein correlation function and the corresponding fit, similarly to Fig. 4, but for K T = 0.10 − 0.20 GeV/c.Fit range in q in this case was 0.035 GeV/c to 0.8 GeV/c.

Figure 10 :
Figure 10: Measured Bose-Einstein correlation function and the corresponding fit, similarly to Fig. 4, but for K T = 0.35 − 0.60 GeV/c.Fit range in q in this case was 0.049 GeV/c to 0.8 GeV/c.
Figure2: The dE/dx measurement for negatively charged (left panel) and positively charged (right panel) particles versus natural logarithm of momentum p (in laboratory frame).The two lines represent the interpolated selection boundaries based on Gaussian fits to dE/dx distributions at several momenta, and are 3 standard deviations from the pion mean for negatively charged particles, while for positively charged particles, the lower line is 1.5 standard deviations distance from the mean, to remove the more significant kaon and proton contributions present in this case.