Femtoscopic correlation measurement with symmetric L\'evy-type source at NA61/SHINE

Measuring quantum-statistical, femtoscopic (including final state interactions) momentum correlations with final state interactions in high-energy nucleus-nucleus collisions reveal the space-time structure of the particle-emitting source created. In this paper, we report NA61/SHINE measurements of femtoscopic correlations of identified pion pairs and describe said correlations based on symmetric L\'evy-type sources in Ar+Sc collisions at 150A GeV/c. We investigate the transverse mass dependence of the L\'evy-type source parameters and discuss their possible interpretations.


Introduction
The NA61/SHINE is a fixed target experiment using a large acceptance hadron spectrometer located in the North Area H2 beam line of the CERN Super Proton Synchrotron accelerator [1]. Its main goals include the investigation and mapping of the phase diagram of strongly interacting matter, as well as measuring cross sections of processes relevant for cosmic rays and neutrino physics. In this paper, we are focusing on mapping the QCD phase diagram. In order to accomplish this, NA61/SHINE performs measurements of different collision systems at multiple energies. The experiment provides excellent tracking down to p T = 0 GeV/c. This performance is achieved by using four large Time Projection Chambers (TPC's), which cover the full forward hemisphere. The experiment also features a modular calorimeter, called the Projectile Spectator Detector. It is located on the beam axis, after the TPC's, and measures the forward energy which determines the collision centrality of the events. A setup of the NA61/SHINE detector system is shown in Fig. 1.
The search for the critical endpoint (CEP) and investigation of the QCD phase diagram requires analysis at different temperatures and baryon-chemical potentials. To study, we need to map the phase diagram using different system sizes at various energies. NA61/SHINE investigations cover several beam momenta (13A, 20A, 30A, 40A, 75A and 150A GeV/c) and collision systems (p+p,p+Pb,Be+Be,Ar+Sc,Xe+La,Pb+Pb). In this paper, we describe the femtoscopic correlations of identical pions emitted from central Ar+Sc collisions at beam momentum of 150A GeV/c. This field is often called femtoscopy as it reveals the femtometer scale structure of particle production.

Femtoscopy with Lévy shaped sources
The method of quantum-statistical (Bose-Einstein) correlations is based on the work of R. Hanbury Brown and R. Q. Twiss (HBT) [2], who applied it first in astrophysical intensity correlation measurements. The method was developed to determine the apparent angular diameter of stellar objects. Shortly afterwards, a similar quantum-statistical method was applied in momentum correlation measurements for proton-antiproton collisions [3,4] by Goldhaber and collaborators. Their objective was to understand pion-pion correlations and gain information on the radius, R, of the interaction volume in high-energy particle collisions. The key relationship for measuring Bose-Einstein correlations shows that the spatial arXiv:2306.08696v2 [nucl-ex] 28 Jun 2023 momentum correlations, C(q), are related to the properties of the particle emitting source, S(x), that describes the probability density of particle creation for a relative coordinate x as: whereS(q) is the Fourier transform of S(x), and q is the relative momentum of the particle pair (with the dependence on the average momentum, K, of the pair suppressed and described in more detail in [5]). The usual assumption for the shape of the source based on the central limit theorem, is a Gaussian. However, such Gaussian shaped sources lead to Gaussian correlation functions. A more general assumption is the Lévy distribution [6,7]. It exhibits a power-law tail and includes a Gaussian limit, as well. Correlation functions based on this approach have been shown to describe data from different experiments, such as LEP [8], RHIC [9], and LHC [10,11] quite well. Several phenomena could explain the appearance of Lévy shaped sources. The non-Gaussianity of the source could be attributed to critical fluctuations and the emergence of spatial correlations on a large scale, which may indicate the existence of similar sources with power-law tails [12]. Further reasons include the fractal structure of QCD jets [13]. In this paper, the measured femtoscopic correlation (including final state interaction) with spherically symmetric Lévy distributions is defined as: where R is the Lévy scale parameter and α defined as the Lévy stability index. In addition, ⃗ ζ is the three-dimensional integration variable with dimensions of MeV/c and ⃗ r is the vector of spatial coordinates. There are two special cases where the distribution can be expressed analytically. One such case is, the already mentioned, Gaussian distribution for α = 2. Besides this, the α = 1 case leads to a Cauchy distribution. An important difference between Lévy distributions and Gaussians is the presence of a power-law tail ∼ r −(d−2+α) in case of α < 2, where d represents the number of spatial dimensions. With the assumption of Lévy sources, the femtoscopic correlation functions can be expressed in the following way: C 2 (q) is a stretched type of exponential, where the λ intercept parameter is defined as: At vanishing relative momentum, the correlation function has a value of 1 + λ. This value is not accessible in the measurements and extrapolation from the region, when two tracks are experimentally resolved, is needed. However, it is commonly observed that the intercept parameter λ is less than 1. The core-halo model, explained in Refs. [14,15], can provide some insights into this parameter. The model assumes that the source S is made up of two parts, the core and the halo (S core and S halo ), respectively. The core contains pions created directly from hadronic freezeout or from extremely short lived (strongly decaying) resonances. The halo consists of pions created from longer-lived resonances and the general background. It may extend to thousands of femtometers, while core part has a size of around a few femtometers. In this picture, the λ parameter turns out to be connected to the ratio of the core and the halo as: Then, one can modify the correlation function to take the effect of the halo into account, by utilizing the Bowler-Sinyukov method [16,17] as: The halo part contributes at very small values of relative momenta, q. Therefore, it does not affect the source radii of the core part [18]. It is well known that critical points are characterized by critical exponents. One, in particular, is related to spatial correlations by the exponent denoted as η. The appearance of the parameter can be explained by the second-order phase transition at the CEP, where fluctuations will appear at all scales causing the spatial correlation function to exhibit a power-law tail with an exponent of −(d − 2 + η), with d denoting the dimension. The Lévy exponent, α, is the exponent in the case of, previously defined, Lévy distributed sources and will also exhibit a power-law tail, with an exponent of −d + 2 − α [19]. Hence, α was suggested to be directly related to or being explicitly equal to, the critical exponent, η [12], in absence of other phenomena affecting the source shape. This is the basis of the idea connecting α and η. However, the Lévy-shape of the source can be attributed to several different factors besides critical phenomena, including QCD jets, anomalous diffusion, critical phenomena, and others [6,7,13,20,21]. Hence, while a non-monotonic behavior of α is expected near the critical point, a detailed understanding of the collision energy and system size dependence of α is needed to draw conclusions about the critical point.
It has been suggested that the universality class of QCD is the same as that of the 3D Ising model [22,23]. The value of η in the 3D Ising model is 0.03631 ± 0.00003 [24]. An alternative solution is to use the universality class of the 3D Ising model with a random external field, which yields an η value of 0.50 ± 0.05 [25]. The statements mentioned suggest that α would decrease to 0.50, or below, at the vicinity of the CEP. To confirm this, measurements of α are needed in different collision systems at various energies.
In this analysis, we are dealing with like-charged particles that are influenced by Coulomb repulsion. The final state Coulomb effect has been neglected in the previously defined correlation function. Thus, the correlation function in Eq. (6) that lacks this effect will be denoted as C 0 2 (q) from now on. The correction necessitated by this effect can be done by simply taking the ratio of C Coul 2 and C 0 2 (q): where C Coul 2 (q) is the interference of solutions of the two-particle Schrödinger equation; with a Coulomb-potential [26,27]. The numerator in Eq. (7) cannot be calculated analytically and requires a large numerical effort to estimate.
An approximate formula for C Coul 2 (q) was obtained in Ref. [10] for the case of Cauchyshaped sources. However, a more precise treatment is required due to our assumption of Lévy-shaped sources. We are utilizing a new method in our analysis for estimating the effect of Coulomb repulsion. The treatment includes the numerical calculation presented in Refs. [27,28], the parametrization of its results, and, finally; the parametrization of the dependence of the physical parameters R, λ, and α. Thus Eq. (6) is modified as: where N is introduced as normalization parameter and K Coulomb (q) denotes the Coulomb correction.
It is important to highlight that the Coulomb correction is calculated in the paircenter-of-mass (PCMS) system, while the measurement is often done in the longitudinally co-moving system (LCMS). The assumption of Coulomb correction in the one-dimensional HBT in LCMS picture is that the shape of the source is spherical, i.e. R out = R side = R long = R ≡ R LCMS . The shape of the source, however, is spherical in the LCMS and not in the PCMS. Therefore, an approximate one-dimensional PCMS size parameter is needed. A study was done, where an average PCMS radius of was calculated [29], with β T = K T m T .

Measurement details
For the measurements this paper is based on, we analyzed Ar+Sc collisions at 150A GeV/c beam momentum in the 0-10% most central events. The available data set contains around 2.7 million events, which was reduced to around 700 000 events in the analysis. The following paragraph describes the various event, track and pair selection performed on the data. As mentioned, we have selected 0-10% of the most central events by measuring the energy contained in the projectile remnants, with the Projectile Spectator Detector (PSD). We have also selected events where no off-time beam particle was detected. These are particles which come within the drift time of the chambers which are not emptied out. Furthermore, all the events chosen have between ±10 cm as that is the maximal distance between the main vertex z position and the center of the scandium target. Particle identification was handled by using dE/dx energy deposit in the TPC gas. The tracks were extrapolated to the interaction plane and matched with the distance against the interaction point. If the distance was ≤ 4 cm in the horizontal plane and ≤ 2 cm in the vertical plane, the track was kept. Moreover, track splitting was handled by selecting the ratio of the total number of reconstructed points on the track to the potential number of points to be between 0.5 and 1.0. Finally, to counter track merging, we have used a selection on the momentum space distance between the two tracks. It uses non-standard momentum coordinates, s x = p x /p xz , s y = p y /p xz and ρ = 1/p xz .
We then analyzed the combination of negative pion pairs and positive pion pairs. It is important to note that this analysis was done with a one-dimensional relative momentum variable q, calculated in LCMS. These pairs were sorted into eight K T (average transverse momentum of the pair) bins in the range of 0-450 MeV/c. In each momentum bin, the relative momentum distribution of coincident pion pairs were obtained. Let us call this the actual pair momentum difference distribution, A(q). A(q) contains quantum-statistical correlations, as well as many other residual effects related to kinematics and acceptance. The effects can be removed by constructing a combinatorial background pair distribution, denoted as B(q), which is measured in the same K T or m T intervals as the A(q) distribution. The method we use involves randomly selecting the same number of particles as the multiplicity of the actual event. The selected particles are from other events of similar  parameters and each is selected from a different event. Let us call the pair momentum difference distribution made from this method as, the background distribution, B(q). This, by construction, enables us to have an uncorrelated pool of events. Then the correlation function is calculated as in a [q 1 , q 2 ] range where quantum-statistical correlations are not expected. Upon adding the contribution from the background, the previously defined Eq. (8) is modified as follows: with N being a normalization parameter responsible for the proper normalization of the A(q)/B(q) ratio, ε describing the linearity of background, and K Coulomb (q) being the Coulomb correction. We then use this formula to our measured C 2 (q) as shown in Fig. 2.
We then use this formula to our measured C 2 (q) as shown in Fig. 2. To determine the fitting range where the effects of detector resolution do not play a significant role, we have used EPOS simulation [30] with GEANT3 for particle propagation [31]. Note that in the low-q region, the fit does not describe the data. This can be explained, according to Monte Carlo simulations of the detector response, by the limited resolution of pairs with small relative momentum.

Results
The three physical parameters (α, R, and λ) were measured in eight bins of pair transverse momentum, K T . The three mentioned parameters were obtained through fitting the measured correlation functions with the formula shown in Eq. (11). The results were investigated regarding their transverse momentum dependence. In the following, we report on the transverse mass dependence of α, λ, and R; where transverse mass is expressed as m T = m 2 π c 4 + K 2 T c 2 , with m π being the pion mass. As explained above, the shape of the source is often assumed to be Gaussian. The Lévy stability exponent, α, can be used to extract the shape of the tail of the source. Our results, shown in Fig. 3, yield values for α between 1.5 and 2.0, which imply a source closer to the Gaussian shape than the one in Be+Be collisions [5], but are still significantly lower than the α = 2 (Gaussian) case. The observed α parameter is also significantly higher than the conjectured value at the critical point (α = 0.5). Altogether, these results suggest that measured correlation functions align with the assumption of a Lévy source, indicating that it is more advantageous over the Gaussian assumption. Further studies are ongoing 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.
As a second parameter, let us look at the Lévy scale parameter, R, visible in Fig. 4. It determines the length of homogeneity of the pion emitting source. The parameter R depends on the transverse-mass as 1/ √ m T . This can be derived using simple, hydrodynamical predictions for Gaussian sources [15,32]: where u 2 T is the average, transverse expansion and T 0 is the hadronisation temperature. In our case, rather surprisingly, despite the non-Gaussian nature of the source, this formula works in describing the measured femtoscopic radii. More precisely, as mentioned above, observing an R ∼ 1/ √ m T is particularly interesting as this type of m T dependence should rise in case of Gaussian sources (α = 2) [33]. It is not entirely clear why this happens, the indicated m T dependence could form in the QGP or at a later stage. This phenomenon was also observed at RHIC [9] and in simulations at RHIC and LHC energies [20,21]. The final parameter being investigated is the intercept (also known as the correlation strength) parameter, λ, defined in Eq. (5). The dependence of λ on m T is shown in Fig. 5. One may observe a slight dependence on m T , but this can still can be considered constant in the investigated range. When compared to measurements from RHIC Au+Au collisions [9,34,35] and from SPS Pb+Pb interactions [36,37], an interesting phenomenon is observed. At the energies of SPS, there is no visible "hole" at lower m T values, but at RHIC energies, the "hole" appears at m T values of around a few hundred MeV. This "hole" was interpreted in Refs. [9] and [38] to be a sign of in-medium mass modification. The results presented in Fig. 5, at the given statistical precision, do not indicate the presence of such a low-m T hole. This trend might imply that this phenomenon can be turned off at SPS. Furthermore, it can be highlighted that the λ values we obtained are significantly below unity. A possible answer can be given by the halo part of the core-halo model. It may indicate that a significant fraction of pions are the decay products of long-lived resonances.

Conclusions
In the report above, we discussed the NA61/SHINE measurement of one-dimensional, identified, two-pion, femtoscopic correlation functions; in the 0-10% most central Ar+Sc collisions at 150A GeV/c. We discussed the transverse mass dependencies of the Lévy source parameters. Results on the Lévy scale parameter, α, showed a significant deviation from Gaussian sources and are not in the vicinity of the conjectured value at the critical point. The Lévy scale parameter, R, shows a visible decrease with m T . The correlation strength parameter, λ, does not show any significant m T dependence, but, maps different patterns at RHIC and similar trends at SPS energies. With these results at hand, we plan to measure Bose-Einstein correlations in larger systems, as well as at smaller energies, to continue mapping the phase diagram of the strongly interacting matter.