Kaon femtoscopy with L\'evy-stable sources from $\sqrt{s_{_{\rm NN}}} = 200$ GeV $\mathrm{Au}+\mathrm{Au}$ collisions at RHIC

Femtoscopy has the capacity to probe the space-time geometry of the particle-emitting source in heavy-ion collisions. In particular, femtoscopy of like-sign kaon-pairs may shed light on the origin of non-Gaussianity of the spatial emission probability density. The momentum-correlations between like-sign kaon-pairs are, hence, measured in data recorded by the STAR experiment, from $\sqrt{s_{_{\rm NN}}} = 200$ GeV $\mathrm{Au}+\mathrm{Au}$ collisions at RHIC, BNL. Preliminary results hint at the, possible, existence of non-Gaussian, L\'evy-stable sources; and signal the, likely, presence of an anomalous diffusion process; for the identically-charged kaon-pairs so produced. More statistically significant studies, at lower centre-of-mass energies, may contribute to the search for the critical end point of QCD as well.


Introduction
Following the discovery of the quark-gluon plasma, one of the main thrusts of high-energy nuclear-physics has been the understanding and exploration of the space-time geometry of the particle-emitting source created in heavy-ion collisions [1]. The quantity mainly investigated to this end is the two-particle source function, sometimes called the spatial correlation function (CF) or the pair-source distribution. Even though this quantity is not easy to reconstruct experimentally, detailed studies of its behaviour are merited by a multitude of reasons; including its connections to hydrodynamic expansion [2,3], critical phenomena [4], light nuclei formation [5], etc. Phenomenological studies and experimental analyses, both, emphasise the importance of describing the shape of the source function. Earlier, hydrodynamical, studies undertaken seemed to suggest a Gaussian source-shape [3,6]. Multiple measurements were also conducted with the Gaussian assumption [27,7]. But, recent source-imaging studies suggest that the two-particle source function for pions has a long-range component, obeying a power-law behaviour [8,9,10,28,11,12].
Femtoscopy [29] is the sub-field of high-energy heavy-ions physics that allows the investigation of the space-time structure of femtometre-scale processes encountered in high-energy nuclear-and particle-physics experiments. Femtoscopic correlations in heavy-ion collisions are currently understood to be caused partly by Bose-Einstein statistics [30,31,32,33]. Alternatively, they are called Hanbury-Brown-Twiss (HBT) correlations in recognition of pioneering works by Hanbury-Brown and Twiss [34,35,36] on intensity interferometry in the field of observational astronomy to extract the apparent angular sizes of stars from correlations between the signals of two detectors. Additionally, correlations can arise out of final-state interactions, like electromagnetic interactions and strong-interactions, undergone by the investigated particles. These correlations between pairs of identical bosons can, hence, be used to explore the properties of the matter created in heavy-ion collisions and to map the geometry of the particle-emitting source [1].

Correlations
Femtoscopy works on the principle that the momentum-correlation function of a pair of particles is related to the probability density of particle creation at a space-time point X, for a particle with four-momentum P . This probability density, S(X, P ), is also called the sourcefunction. Defining N 1 (P ) -given by multiplying the particle-creation probability density with n , the average number of particles -as the invariant-momentum distribution and N 2 (P 1 , P 2 ) -given by multiplying the pair-creation probability density with n(n − 1) , the average number of pairs -as the pair-momentum distribution, the two-particle correlation function can be written as [37]: where with ψ P1,P2 (X 1 , X 2 ) being the two-particle wave-function that simplifies to: in the interaction-free case for bosons, when only the quantum-statistical effects are taken into account. Thus, the correlation function can be re-defined as [13]: where denote the pair-momentum difference, the average momentum and the Fourier-transform (FT) of the source, respectively, assuming that Q K holds for the kinematic range under investigation. The measurement of the correlation functions is done with respect to Q; over a range of well-defined K-values; and then the properties of the correlation functions are analysed as functions of the average-K, for each of those ranges.
A significant fraction of the particles created in a heavy-ion collision is secondary, i.e. they come from decays. Whereas the primordial particles are created directly from the hydrodynamic expansion of the collision-volume, the decay-particles arise from interactions that take place much later. Hence, the source can be assumed to consist of the following two components [14]: 1. a core, S C (X, K) -consisting of primordial particles created by the hydrodynamicallyexpanding, strongly-interacting quark-gluon plasma; along with the decays of resonances with half-lives less than a few fm/c; and 2. a halo, S H (X, K) -consisting of the products created by the decay of long-lived resonances, including but not limited to η, η , K 0 S , ω, making it possible to decompose S(X, K) as [14]: As explained in detail in Ref. [14], the Fourier-transformed emission function of the halo vanishes for resolvable relative momenta, i.e. the Q-values that lie in the experimentally achievable region. Hence, the halo does not contribute to the measured correlation function, which in turn is determined by the core component. Hence, the measured correlation function, when extrapolated to Q = 0, does not take a value of 2 as expected from Eqn. (4), but takes the form: This "intercept parameter" of the correlation function, also called the correlation-strength, λ(K), may depend on the pair-momentum K and can be understood on the basis of the corehalo model by re-writing the correlation function using Eqns. (4), (6) and (7): with N C (K) = S C (X, K)d 4 X and N H (K) = S H (X, K)d 4 X being, respectively, the core's and the halo's contributions and with λ(K) being: Realising that X ≡ X( r, t), the spatial correlation function: can be used to re-write C(Q, K) as [13]:

Lévy-distribution
Usually, the shape of the source distribution is assumed to be Gaussian. However, evidence of a non-Gaussian source for correlated pions has been found in multiple studies, necessitating a generalisation of the distribution to a Lévy-stable one [38]: where L is the one-dimensional Lévy-function, R is the Lévy-scale-parameter, λ is the correlationstrength, 0 < α ≤ 2 is the Lévy-exponent and Q is the integration-variable. These parameters are, generally, understood to depend on K. Salient features of the distribution include moments greater than α being undefined and D( r, K) necessarily having a Lévy-shape, with the same α, in cases where S(X, K) is Lévy-shaped. The distribution exhibits a power-law behaviour for α < 2, with α = 1 representing a Cauchy distribution and α = 2 representing a Gaussian distribution. Multiple factors; such as anomalous diffusion, jet-fragmentation and proximity to the critical end point (CEP); can contribute to the appearance of Lévy-stable sources. However, the high-multiplicity, nucleon-on-nucleon nature of the heavy-ion collisions analysed makes it unlikely for jet-fragmentation to be the dominant reason for the appearance of Lévy-sources in this study -as it has been identified as the cause behind Lévy-stable sources in e + e − collisions at LEP [11]. On the other hand, the high centre-of-mass-energy of the collisions explored here rules out the possibility of the system being close to the critical end point [15,4,16]. Interestingly, it is trivial to establish that α is related to one of the critical exponents, η. In case of a second-order phase transition, the η exponent describes the power-law behaviour of the spatial correlation function at the critical end point with an exponent of −(d − 2 + η); where d is the dimension. In a three-dimensional analysis, with d = 3, this exponent would compute to −(1 + η). However, it is established that the three-dimensional Lévy-distribution describes the power-law tail of the spatial correlation function with an exponent of −(1 + α). Thus, comparing the exponents at the critical end point, it can be easily seen that the Lévyexponent, α, is identical to the critical exponent, η; a conjecture explored in Ref. [17]. The second-order QCD phase transition is expected to be in the same universality class as the threedimensional Ising model and in that case, the η-exponent has a value of 0.03631±0.00003 at the critical end point [18]. However, the universality class of the random-field, three-dimensional Ising model may also be of relevance here, where the value of η is 0.50±0.05 [39]. Thus, extracting α at collision energies lower than the ones used in this analysis; while taking into account finite-size and finite-time effects; might yield insightful information about the nature of the quark-hadron phase transition and shed light on the location of the critical end point in the QCD phase diagram [40,41,4,19,42].
As coordinate-space distributions extracted from experimental data show a heavy tail, the limitations of the hydrodynamical approach -assuming idealised freeze-out with a sudden jump in the mean-free-path from zero to infinity -become clear. This requires the, more realistic, approach using hadronic re-scattering where the system cools as it dilutes with an expanding hadron gas, its mean-free-path diverges to infinity in a finite time-interval and rescattering occurs in the presence of a time-dependent mean-free-path. This signals the existence of anomalous diffusion -experimentally observed as power-law-shaped tails in coordinatespace distributions of the source -in the system; as opposed to normal diffusion -with the Gaussian source exhibiting a strongly-decaying tail -caused by the Brownian motion of the particles constituting the system. One-dimensional, normal diffusion's momentum-space diffusion-equation is [16]: where k N is the normal-diffusion constant, q is the momentum, t is the time and W (q, t) is the momentum-space probability distribution. The coordinate-space solution to Eqn. (14) is given by the Gaussian expression: For anomalous diffusion, the coordinate-space diffusion-equation; in terms of the spatial probability distribution W (x, v, t); is the generalised Fokker-Planck equation [16]: Here, η α is the generalised friction constant of dimension is the Riemann-Liouville operator: and L FP is the Fokker-Planck operator: with V (x) being related to the force F (x) by F (x) = − dV (x) dx ; as explained in Refs. [43,20]. The momentum-space solution to Eqn. (16) is given by: This W (q, t) happens to be the FT; or the characteristic function; of Lévy-stable sourcedistributions, with α being the Lévy-exponent from Eqn. (13) and k A being the anomalousdiffusion constant. If a centred, spherically symmetric, Lévy-stable source-distribution is assumed; and all final-state interactions are neglected; the one-dimensional, two-particle correlation function takes the simplified form: with λ as the correlation-strength from Eqns. (10) and (8), R as the Lévy-scale, α as the Lévyexponent and q as the absolute value of the three-momentum-difference in the longitudinally co-moving system (LCMS) [8]: R can be interpreted as the homogeneity-length of the particle-species, while α represents the extent of the anomalous diffusion occurring in the system. The spherical symmetry in q LCMS is ideal for a one-dimensional analysis of a three-dimensional, spherically symmetric system. Subsequent measurements, with higher precision, might necessitate a move towards a full, three-dimensional analysis. Until then, the approximations used in Ref. [44] may be utilised for a preliminary analysis. Momentum-correlations of like-sign kaon-pairs, at √ s NN = 200 GeV, can be utilised to calculate C(q) and to, thus, ascertain the shape of the pair-source distribution. If anomalous diffusion is the sole source of non-Gaussianity, then it is expected that the extent of the anomaly will depend on the total cross-section, or equivalently, on the mean-free-path. Since the meanfree-path is larger for kaons than for pions, the former's diffusion is expected to be more anomalous in the hadron gas. Hence, the α for kaons, α κ , is expected to be smaller [16]. Thus, measuring α κ may shed light on the role of anomalous diffusion in the hadron gas as the origin of the appearance of Lévy-distributions.

Measurement
The data used for this analysis is obtained from RHIC's gold-on-gold collisions at 200 GeV centre-of-mass energy per nucleon pair; performed in 2016; as measured by the solenoidal tracker at RHIC (STAR) experiment. Different particle species, depending on their mass and charge, produce different shapes when their ionisation-energy-loss, dE/dx, is plotted as a function of momentum. These shapes can help distinguish the particle species from each other and help isolate the kaons, as observed in Fig. 1. For this investigation, the analysis processes about 410 million events in the 0-30% centrality range. They undergo strict track-and pair-selection criteria that, following Ref. [45], including the identification of kaons via energyloss measured by the time-projection chamber (TPC); pair-selection based on the fraction-ofmerged-hits (FMH) and the splitting-level (SL); and limitations on the track's momentum, rapidity and distance-of-closest-approach (DCA).
The one-dimensional, like-sign, two-kaon, femtoscopic correlation functions are then experimentally constructed, using the event-mixing technique [21]. The actual pair-distribution, A(q), is formed from kaon-pairs, with members of the pair belonging to the same event. This distribution is affected by various effects, such as kinematics and acceptance. To correct for these effects, a background-distribution, B(q), is constructed with the pairs, where the members originate from separate events. In this analysis, the method for event-mixing described in Refs. [8,11] is used to correct for residual correlations. For a set of event-classes based on the centrality and the location of the collision-vertex, a background event pool is established. Then, for each real event, a mixed event of the same event-class is created from this pool, making sure that each particle in this mixed event belongs to a different real event. Subsequently, pairs within the mixed event contribute to the formation the aforementioned background-distribution, B(q). Finally, the pre-normalised correlation function is calculated as: for three different ranges of transverse-mass, m T , defined as m T = m 2 + (K T /c) 2 ; with m as the kaon-mass and K T as the average transverse-momentum of the pair. The normalisation integrals are performed over a range where the correlation function is not supposed to exhibit quantum-statistical features. It is to be noted that, the method outlined here is applied to pairs belonging to a given range of average momenta. Furthermore, in the event-mixing technique described above, the number of actual and background pairs is the same, aside from the effect of the pair-selection criteria mentioned earlier.
With the momentum correlations, obtained from experimental data, measured and the empirical values of the correlation function calculated; preparations are made for fitting the Lévy-function, detailed in Eqn. (20), to C(q). However, the assumptions behind Eqn. (20) make it unsuitable for a direct fit to experimentally-obtained data. As mentioned in Sec. 1, final-state interactions have considerable effects on the momentum correlations between likesign kaon-pairs and hence, they need to be imbibed into the analysis as corrections to Eqn. (20) -to make it suitable and, physically, relevant as a fit for the correlation function obtained above. Multiple factors can contribute to the final-state modifications of the momentumcorrelations, but leading amongst them are Coulomb interactions; since a gas of charged hadrons can never be, entirely, devoid of Coulomb repulsion.
The final-state Coulomb-interactions are incorporated into the CF by using the Bowler-Sinyukov formula, that includes a correction-term for Coulomb-repulsion, as [46,22]: where N is a normalisation factor, ε is responsible for a linear, residual, long-range background and K is the Coulomb correction [13]: with D( r) being the, simplified, spatial pair-distribution and ψ Coul q ( r) being the solution to the two-particle Schrödinger-equation, in the presence of a Coulomb potential. In this study, K(q; α, R) for kaons is calculated numerically employing the procedure used in Refs. [13,28,23]. The inclusion of other final-state contributions, like the strong interaction, can resolve the possible underestimation -regarding R and λ -and overestimation -regarding α -of the Lévy-parameters [24]; but, the statistical significance of such precise corrections turns out to be negligible in the context of the current measurement.

Results & Discussion
As illustrated by Fig. 2, the Coulomb-corrected Lévy-distribution function agrees with the measured C(q) over the entire q LCMS -range. The femtoscopic peak [32,33] and the Coulombhole [22] are both observed, as expected. The values of the normalisation factor, N , and the linear-background factor, ε, are observed to be close to 1 and 0, respectively.
The systematic uncertainties are obtained by combining the uncertainties arising from variations in the event-and pair-selection criteria, denoted by ∆ cuts , mentioned above and those arising out of variations to the range of the fit, denoted by ∆ fits . At this preliminary stage, systematic uncertainties arising out of variations in the track-selection criteria are not included. Thus, the final systematic-uncertainties, ∆ total , are obtained as: Fig. 3 shows the kaon-homogeneity length, R; otherwise known as the Lévy-scale; as a function of m T . It is observed to exhibit large, systematic uncertainties; a very weak dependence on m T ; and a possible decrease w.r.t. it, as found from previous studies [7,25,3,2,6,12]. Note, however, that the hydrodynamical studies predicting the decrease of the Lévy-scale, as a function of m T , are based on the Gaussian-source assumption [3,2]. Hence, more investigations on this topic are in order. The extracted values of the Lévy-scale, in this charged-kaon analysis, are also found to be similar to PHENIX's like-sign pion results [8], with R π ∼ 5-7 fm for the m Trange of 600-700 MeV/c 2 . A more detailed comparison of the m T -dependence of Lévy-scales of different particle-species could shed light on the origin of the appearance of Lévy-stable sources, given that; from calculations based on hydrodynamics; a species-independent m T -scaling was predicted in Ref. [26].    Figure 5: α, as a function of m T , for 0-30% centrality. The hollow, blue squares denote positively-charged kaon-pairs and the solid, blue circles denote negatively-charged kaon-pairs; along with their error bars. The systematic uncertainties are shown as hollow (K + K + ) and shaded (K − K − ) rectangles.
The intercept of the correlation function -the correlation strength, λ -is shown in Fig.  4. Values extracted from the fits show that it's close to unity; as is to be expected from the small fraction of decay-kaons present in the system.
The extent of the anomalous diffusion might be gleaned from the Lévy-exponent, α, as seen in Fig. 5. It also illustrates the values corresponding to the Gaussian and Cauchy distributions, with dashed and dotted blue lines, respectively. The Lévy-exponent is observed to have values between those two extremes, thereby indicating a power-law behaviour and the, likely, existence of anomalous diffusion. The extracted values of α ∼ 1.0-1.5 for kaons are, again, similar to PHENIX's pion-results; with α π ∼ 1.2 in the same transverse-mass range. α κ not being smaller than α π hints at the existence of other factors, on top of anomalous diffusion, contributing to the appearance of non-Gaussian source-shapes. However, the current statistics prevent the drawing of more definitive conclusions.

Summary & Outlook
Preliminary analysis of data collected by STAR, from RHIC's 2016 BES √ s NN = 200 GeV Au + Au collisions, suggests a non-Gaussian, Lévy-stable source-shape for pairs of the identically-charged kaons produced in the collisions. The Lévy-stability-exponent α κ , is observed to be comparable to that of like-sign pion-pairs obtained from PHENIX. However, anomalous diffusion may not be solely responsible for the heavy tails observed in the source distributions, as suggested by α κ being comparable to α π . It is to be noted that, a complete systematic uncertainty analysis; which is currently ongoing; is required to achieve definitive conclusions about any and all claims made. Seeing as Lévy-stable sources can arise in strongly-interacting systems due to their proximity to the QCD critical end point, at higher chemical potentials, similar studies at lower beam energies would likely strengthen the search for the QCD critical end point.