Lévy-stable two-pion Bose-Einstein correlations in $\sqrt{s_{{}_{NN}}}=200$ GeV Au$+$Au collisions

This analysis uses the beam-beam counters (BBC) for event characterization. Its two arms (“North” and “South”) are located at $\pm 144$ cm along the beam axis ($z$ axis) from the center of PHENIX, corresponding to the 3.0$<|\eta|<$3.9 pseudorapidity interval. Each arm of the BBC comprises 64 quartz Čerenkov counters, covering $2\pi$ in azimuth. They provide minimum-bias (MB) triggering; the MB trigger condition requires at least two hits in coincidence in both BBC arms, thus capturing $92\pm 3$% of the total Au$+$Au inelastic cross section Adare et al. (2013). The charge sum in both BBC arms is used for event centrality determination. The BBCs also measure the average hit time in the north and south arm photomultipliers (PMTs), thus providing collision vertex position measurements along the $z$ direction (from the hit time difference) as well as initial timing information for the collision. With an intrinsic timing resolution of $\approx$40 ps, the $z$-vertex resolution is $\approx$0.5 cm and $\approx$1.5 cm in central and peripheral Au$+$Au collisions, respectively.

PHENIX has two central arm spectrometers (“east” and “west”), each covering $|\eta|<0.35$ in pseudorapidity and $\Delta\varphi=\pi/2$ in azimuth, as seen in Fig. 1. In each central arm, charged particle tracks are reconstructed using hit information from the drift chamber (DC), the first layer of pad chambers (PC1) and the collision $z$-vertex position measured by the BBC Adcox et al. (2003b).

The DCs are located at a radial distance of 202–246 cm from the beam axis. They provide trajectory measurement in the transverse plane, with an angular resolution of $\approx$1 mrad. The PC1s are multiwire proportional chambers with pad readout, located immediately behind the DCs. They provide track position measurement both in the $\varphi$ and in the $z$ direction, with a $z$-resolution of $\approx$1.7 mm.

The PHENIX central arm spectrometer magnet generates a magnetic field approximately parallel to the beam line. It comprises two pairs of independently operable concentric coils, an inner and an outer coil pair, located at radial distances of $\approx$60 cm and $\approx$180 cm, respectively. The DCs are positioned so that they are in the reduced field region. Charged-particle-momentum determination is enabled by the measurement of the bending of the track in the magnetic field. The transverse momentum $p_{T}$ is determined by the bending angle measured by the DC, while the polar angle of the momentum is determined by the $z$ coordinate measured by PC1 and the $z$-vertex coordinate from the BBC. Reconstructed tracks are then projected to the outer detectors used for track verification and timing measurement.

Because at not too low $p_{T}$ the momentum resolution is governed mainly by the angular resolution of the DC, high bending fields are desirable. Thus usually the two coil pairs are operated with currents flowing in the same direction (this is called “$++$” or “$--$” mode), to achieve the designed maximum total field integral of $\int B\cdot dl\approx$1.1 T m (this is the relevant quantity for the bending, and in turn for the momentum measurement).

In 2010, the Hadron Blind Detector (HBD), a specialized Čerenkov counter located around the nominal collision point for the measurement of dielectron pairs, was installed Anderson et al. (2011). The operation of the HBD required a field-free region around the collision point, which was achieved by running the inner and outer coils in the opposite directions (in “$+-$” or “$-+$” modes). This reduced the field integral to $\approx$40% of its maximum value. However, the present analysis deals with low and intermediate $p_{T}$ hadrons (up to $p_{T}\approx$0.85 GeV$/c$), so high $p_{T}$ momentum resolution is not crucial. (The momentum resolution for $p_{T}$ in the dataset used is estimated to be $\delta p_{T}/p_{T}\approx 1.3\%\oplus 1.2\%\times p_{T}[\mathrm{GeV}/c]$ Adare et al. (2012). The $p_{z}$ momentum resolution has, in addition, a component stemming from the BBC $z$-vertex resolution.) Moreover, the reduced magnetic field had a beneficial side effect for the present analysis. Namely, the low momentum acceptance of this dataset is extended to lower values of transverse momentum, enabling a relatively clean identified pion sample down to $p_{T}\approx$0.2 GeV$/c$. This would have been much harder, if not impossible, with the normal $++$ or $--$ field setting, because of too large bending angles and residual bending outside of the DC nominal radius, which is not taken into account in the standard PHENIX track projection algorithm.

In the present analysis, we identify charged pions by their time of flight from the collision point to the outer detectors. We use the lead-scintillator electromagnetic calorimeter (PbSc) as well as the high resolution time-of-flight detectors (TOF east and TOF west) Aizawa et al. (2003).

The PbSc is a sampling calorimeter located approximately 5.1 m radial distance from the beam axis. It covers $|\eta|<$0.35 in both arms, and in terms of $\varphi$, it covers all $\pi/2$ acceptance of the west arm, and $\pi/4$ (i.e. half) of the east arm, as seen on Fig. 1. It is a finely segmented detector, consisting of 15,552 individual channels (“towers”). After careful tower-by-tower and energy dependent calibration, a timing resolution of $\approx$ 400–600 ps (depending on deposited energy, incident angle, individual channel electronics imperfections, etc.) was achieved for pions. The part of the east arm acceptance not covered by the PbSc is covered by the lead-glass (PbGl) calorimeter, which has a much worse timing resolution for hadrons and thus was not used for the present analysis.

The TOF east detector is also located at approximately a 5.1 m from the beam axis, and covers much of the PbGl acceptance in the east arm. It is made of 960 plastic scintillator slats, with 2 PMTs attached to each side of them. After calibration, the timing resolution was found to be $\approx$150 ps. Adare et al. (2016). The TOF west detector takes advantage of the multigap resistive plate chamber (MRPC) technology. It has two separate panels, each covering $\Delta\varphi\approx\pi/16$ in the west arm, at around 4.8 m radial distance from the beam pipe. Each panel comprises 64 MRPCs and has 256 individual copper readout strips. After calibration, a timing resolution of $\approx$90 ps was achieved.

The MB-triggered data sample used in this analysis comprises $\approx 7.3\times 10^{9}$ $\sqrt{s_{{}_{NN}}}=200$ GeV Au$+$Au events recorded by PHENIX during the 2010 running period. This sample is reduced to $\approx 2.2\times 10^{9}$ events when we apply a 0%–30% centrality selection. The event $z$-vertex position was constrained between $\pm 30$ cm in order to have an efficient BBC response as well as to avoid scattering in the central magnet steel.

We selected tracks of good quality, i.e. those where the DC and PC1 information was unambiguously matched. To reduce in-flight decays as well as random associations between tracks and hits in the PbSc/TOF detectors, a track matching cut of $2\sigma$ was applied for the difference between the projected track position and the closest hit position in these detectors, in both the $\varphi$ and $z$ directions. As part of the systematic uncertainty investigation, we studied the dependence of the final results on these selection criteria.

For the present analysis, a clean sample of identified pions was necessary. Charged pion identification was performed with the help of time-of-flight information ($t$) from the PbSc/TOF detectors and the BBC, as well as using path length information ($L$) from the track model and the momentum value $p$ measured by the DC/PC1. The reconstructed squared mass $m^{2}$ of a track is then

and pions were selected by applying a $2\sigma$ cut in the $m^{2}$ distribution of the PbSc and the TOF detectors. For the $p_{T}$ range of interest in this analysis, the contamination in the pion sample caused by misidentified kaons or protons is negligible. A more important contamination in the pion sample comes from the random association of tracks and hits in the PbSc or the TOF detectors at low momentum, reaching $\approx$2%–3% for the TOF detectors, and as high as 8%–10% for the PbSc at or below $p_{T}\approx 0.2$ GeV$/c$. This background quickly diminishes for even slightly higher $p_{T}$ (at $p_{T}\approx 0.25$ GeV$/c$), as inferred from the observed $m^{2}$ distributions. However, even at low $p_{T}$ this is a gross overestimation of the contamination. Most of the tracks are pions, even those for which the track projection algorithm didn’t find the proper hit because of the residual bending at low momentum. The systematic uncertainty stemming from mis-identified particles is mapped out by varying the mentioned standard $2\sigma$ cut on the $m^{2}$ spectrum of pions, as detailed in Section V. In this analysis, we apply a $p_{T}>0.16$ GeV$/c$ selection, including all identified pions above this threshold into our sample.

In general, the two-particle correlation function $C_{2}(p_{1},p_{2})$ is defined as

where $N_{1}(p_{1})$, $N_{1}(p_{2})$ and $N_{2}(p_{1},p_{2})$ are the one- and two-particle invariant momentum distributions at four-momenta $p_{1}$ and $p_{2}$, and the superscript “spm” denotes that here the correlation function is written as a function of the single particle momenta.

There can be many causes of correlated particle production, such as collective flow, jets, resonance decays, conservation laws. In heavy ion collisions, the main cause of like-sign pion pairs correlation at small relative momentum is the quantum-statistical Bose-Einstein or HBT correlation stemming from the indistinguishability (and thus the symmetrical pair wave-function) of two identical bosons. This source of correlations grows with the mean number of pairs at small relative momentum, which is approximately proportional to the mean multiplicity squared. Other possible sources of correlations (for example pion pair production from resonance decays) increase only linearly with the mean multiplicity. Hence, for the large multiplicity heavy ion collisions, Bose-Einstein correlations dominate the correlation function at small relative momenta.

Experimentally the method of the measurement is the so called event-mixing. To discuss that in this subsection, let us denote any experimental choice for the measure of the two-pion relative momentum by $q$, defining our particular choice later in subsection III.4. In the present subsection we discuss only those properties of the two-pion Bose-Einstein correlation functions that are generally valid, independently of the particular experimental choice of $q$ for the measure of the relative momentum of the pion pair.

Let us define $A(q,K)$ as the actual $q$ distribution of pion pairs for a given average four-momentum $K$, where both members of the pair stem from the same event. Note also that our choice for $K$ is detailed later in Section III.4. This $A(q,K)$ distribution will contain effects which have to be excluded from the Bose-Einstein correlation function (such as resonance decay effects, kinematics, acceptance effects etc.). For this purpose, one defines a background distribution with pairs of pions from different events. Let us denote this background distribution with $B(q,K)$. A usual method is to construct the background distribution by keeping an event pool of a predefined size, and correlating each pion of the investigated event with all same charged pions of the background pool. However, in this case, multiple particle pairs will come from the same event pair. In this analysis we use the method described in Achard et al. (2011) that eliminates any possible residual correlation of this type as well. For each “actual” event, we form a “mixed” event by choosing pions (of the same number as in the actual event for each charge) from other randomly selected events within the background pool (that has to be larger than the maximal multiplicity of pions of a given charge), under the condition that no two tracks may originate from the same event. After this procedure, each “mixed” event comprises pions originating from different events. The background distribution is then created from the (same charge) pairs of this mixed event. It must also be noted that in order for the background event to exhibit the same kinematics and acceptance effects, one has to build the background event from the same event class (i.e. from events of similar centrality and of similar $z$ coordinate of the collision vertex). We used 3% wide centrality and 2 cm wide $z$-vertex bins to achieve that goal.

If we now take the ratio of the actual and the background distributions, we get the prenormalized correlation function as

where the integral is performed over a range where the correlation function is not supposed to exhibit quantum statistical features. Let us note that the method described above is applied to pairs belonging to a given range of average momenta, and in that case $K$ denotes the mean of these average momenta in the given range. Furthermore, in the mixing technique described above, the number of actual and background pairs is the same – aside from the effect of two-track cuts, which is outlined in the next subsection.

When forming pairs to construct the aforementioned actual $A(q)$ and background $B(q)$ pair distributions, one has to take into account detector inefficiencies and peculiarities of the track reconstruction algorithm which sometimes doubles or splits one track into two (creating so-called ghost tracks). It is also possible that two different tracks are not well distinguished when they approach one another too closely. To remove these possible track splitting and track merging effects, we studied track separation distributions in each detector involved, in each of the transverse momentum bins used in this analysis. Then we applied the following cuts in the $\Delta\varphi-\Delta z$ plane (in units of radians and cm, respectively) of pairs of hits in the given detector, associated with track pairs:

We applied these two-track cuts to both the actual and the background sample.

In addition to these cuts, if we found multiple tracks that are associated with hits in the same tower of the PbSc, slat of the TOF east, or strip of the TOF west detector, we removed all but one of them. This ensured that we do not take into account any ghost tracks that would have remained in the sample after the above mentioned pair cuts.

Our analysis method is somewhat different from those of earlier measurements of Bose-Einstein correlations in heavy ion collisions, in particular with respect to the kinematic variables and the application of Lévy-stable distributions. Thus we proceed carefully here and provide a thorough and detailed description of the concepts and procedures that we applied in the determination of the proper kinematic variables and the shape analysis of the Bose-Einstein correlation functions.

The correlation function, as defined in Eq. (2), depends on single particle and pair momentum distributions. These can be calculated in the Wigner function formalism, assuming chaotic particle emission, from the single particle and pair wave functions, as detailed in Refs. Yano and Koonin (1978); Makhlin and Sinyukov (1988); Pratt et al. (1990); Csörgő (2002). For the pair momentum distribution, neglecting dynamical two-particle correlations, one obtains the Yano-Koonin formula Yano and Koonin (1978)

by means of the phase-space density of the particle-emitting source $S(x,p)$, sometimes referred to as “source distribution” or simply as “source”, and $\Psi^{(s)}_{p_{1},p_{2}}(x_{1},x_{2})$, the symmetrized pair wave function. Neglecting final state Coulomb and strong interactions, as well as possible higher order wave-function symmetrization effects on the level of two-particle correlation functions, the pair wave-function is a properly symmetrized plane wave, i.e. in this case,

This approximation in turn leads to the expression of the pure quantum-statistical correlation function ($C^{(0)}_{2}$) as Yano and Koonin (1978); Makhlin and Sinyukov (1988); Pratt et al. (1990); Csörgő (2002)

where complex conjugation is denoted by ^∗, the $(0)$ index signals that the Coulomb effect is not taken into account, the superscript “spm” denotes that the correlation function is written as a function of the single particle momenta, and from now on

stands for the difference of the four-momenta of particles 1 and 2 ($q_{0}$ denotes energy difference, i.e. the zeroth component of the relative four-momentum $q$) and $\widetilde{S}(q,p)$ denotes the Fourier transform of the source

For source distributions and typical kinematic domains encountered in heavy ion collisions, the dependence of $\widetilde{S}(q,p)$ as defined in Eq. (12) is much smoother Lisa et al. (2005) in the original $p$ momentum variable than in the relative momentum $q$, coming from the Fourier transform. Hence, it is customary to apply the $p_{1}\approx p_{2}\approx K$ approximation in Eq. (10), where

is the average four-momentum of the pair ($K_{0}$ denotes the average energy of the pair, i.e. the zeroth component of the average four-momentum $K$). With this,

The validity of these approximations was reviewed in Refs. Csörgő (2002); Wiedemann and Heinz (1999) and for typically exponential single particle spectra the approximation was found to be within 5% of the more detailed and substantiated calculations.

If the above approximations are justified, the two-particle Bose-Einstein correlation function is unity plus a positive definite function of the relative momentum $q$. In the $\sqrt{s_{{}_{NN}}}=200$ GeV 0%–30% centrality Au$+$Au data reported in this analysis, we found that Eq. (14) is consistent with the data; we did not observe the nonpositive definite, oscillatory behavior that was observed in $e^{+}e^{-}$ collisions at LEP Achard et al. (2011), and in $p$$+$$p$ collisions at the LHC Khachatryan et al. (2011); Astalos (2015). Note that in $e^{+}e^{-}$ collisions at LEP and in $p$$+$$p$ collisions at the LHC the smoothness approximation indicated above is not valid, but the Yano-Koonin formula of Eq. (8) still holds Achard et al. (2011); Khachatryan et al. (2011).

In general, as described above, the correlation function depends on four-momenta $p_{1}$ and $p_{2}$ or, equivalently, on $q$ and $K$. However, the Lorentz product of $q$ and $K$ is zero, i.e. $qK=q_{0}K_{0}-{\boldsymbol{q}}{\boldsymbol{K}}=0$. Here $\boldsymbol{q}$ and $\boldsymbol{K}$ are defined as three-vector components of $q$ and $K$ as

This in turn implies

Based on this relation, one may transform the $q$-dependent correlation function to depend on $\boldsymbol{q}$ instead. If the particles contributing to the correlation function are similar in energy, then $K$ is approximately on-shell; thus the correlation function can be measured as a function of $\boldsymbol{K}$ and $\boldsymbol{q}$.

As the dependence on $\boldsymbol{K}$ in heavy ion reactions is typically smoother than on $\boldsymbol{q}$, one may think of $\boldsymbol{q}$ as the “main” kinematic variable. Then one may assume a parameterization of the $\boldsymbol{q}$ dependence, and explore the dependence of the parameters on $\boldsymbol{K}$. Close to midrapidity, instead of $\boldsymbol{K}$, the dependence on

or, alternatively, on the transverse mass

may be investigated, with $m$ being the particle (e.g. pion) mass. Note that the average four-momentum $K$ is not on mass-shell, but $m_{T}$ would be the transverse mass of a particle with momentum $K$. Furthermore, $m_{T}$ also corresponds to the average transverse mass of the particle pair, $M_{T}=0.5(m_{T,1}+m_{T,2})$ in the limit of vanishing relative momentum $|{\bf q}|\rightarrow 0$. As earlier results were frequently given in terms of $K_{T}$, which is a unique function of $m_{T}$ of Eq. (18), we decided to use $m_{T}$ instead of $M_{T}$ to characterize the transverse momentum of a pair of identical pions.

Let us also note that Eq. (14) can be reinterpreted if we introduce the pair distribution as

where $r$ is the pair separation four-vector and $\rho$ is the four-vector of the center of mass of the pair. Then the correlation function can be expressed as

where $\widetilde{D}$ is defined with the Fourier transformation as

Thus the two-particle Bose-Einstein correlation function is connected to the pion pair distribution $D(r,K)$, so this is the quantity that can be reconstructed from two-particle correlation data directly. Different source distributions that keep $D(r,K)$ invariant yield equivalent results from the point of view of two-particle Bose-Einstein correlation measurements.

At any fixed value of the average pair momentum $K$, the correlation function $C_{2}(q,K)$ can be measured as a function of various decompositions of the components of the relative momentum $\boldsymbol{q}$. The Bertsch-Pratt (BP) or side-out-longitudinal decomposition Pratt (1986); Bertsch et al. (1988) is frequently used. Here

with $q_{\mathrm{\rm long}}$ pointing in the beam direction, $q_{\mathrm{\rm out}}$ in the direction of the average transverse momentum $(K_{x},K_{y})$, and the “side” direction orthogonal to these two directions. The transformation to the BP variables corresponds to a rotation in the transverse plane, depending on the direction of the average momentum. For the BP decomposition, it is particularly favorable to use the longitudinal co-moving system (LCMS) of the pair, where the average momentum is perpendicular to the beam axis. Here the BP decomposition of the average momentum is simply $\boldsymbol{K}_{\mathrm{BP}}\equiv(K_{T},0,0)$, as $K_{T}=K_{\mathrm{\rm out}}$, and the temporal information of the source is coupled to the out component of the Bose-Einstein correlation function Wiedemann and Heinz (1999); Csörgő (2002).

However, the Bertsch-Pratt variables require three-dimensional Bose-Einstein correlation measurements, so a detailed shape analysis in terms of them can suffer from a lack of statistical precision. For example, it is very difficult to identify any non-Gaussian structure in a three-dimensional analysis of correlation functions. For this reason, sometimes the two-particle correlation function is measured as a function of a one-dimensional momentum variable Achard et al. (2011); F. Siklér, for the CMS Collaboration (2014). The Lorentz invariant relative momentum, corresponding to the Lorentz length of $q^{\mu}$, is defined as

In the LCMS, using the Bertsch-Pratt variables $q_{\mathrm{inv}}$ is expressed as

where $\beta_{t}=2K_{T}/(E_{1}+E_{2})$ is the “average transverse speed” of the pair.

Let us introduce also the rest frame of the pair, here referred to as pair center-of-mass system (PCMS), and define the relative three-momentum in this system as $\boldsymbol{q}_{\mathrm{PCMS}}$. Then the variable $q_{\mathrm{inv}}$ can be expressed as

Equation (24) shows that $q_{\mathrm{inv}}$ can be very small at moderate $K_{T}$, even for not very small $q_{\mathrm{\rm out}}$ values. It is also well known that the Bertsch-Pratt radii ($R_{\mathrm{\rm out}}$, $R_{\mathrm{\rm side}}$, $R_{\mathrm{\rm long}}$) are of similar magnitude in $\sqrt{s_{{}_{NN}}}=200$ GeV Au$+$Au reactions at RHIC, so the Bose-Einstein correlation functions are nearly spherically symmetric in the LCMS frame Adler et al. (2004); Adams et al. (2005b); Afanasiev et al. (2008, 2009). This also implies that the correlation function boosted to the PCMS frame is definitely not spherically symmetric (especially for intermediate or high $K_{T}$, i.e. for $\beta_{t}$ values approaching 1). The conclusion is that $q_{\mathrm{inv}}$ is not a proper one-dimensional variable of Bose-Einstein correlations of pions in $\sqrt{s_{{}_{NN}}}=200$ GeV Au$+$Au collisions.

We look for a novel one-dimensional variable whose small value is only possible in the case when $q_{\mathrm{\rm out}}$, $q_{\mathrm{\rm side}}$, $q_{\mathrm{\rm long}}$ are all small. Hence, we introduce LCMS three-momentum difference $\boldsymbol{q}_{\mathrm{LCMS}}$ This quantity is invariant for Lorenz boosts in the beam direction. For the sake of simplicity, we hereafter define

which can be expressed with the lab-system components of the individual particle momenta as

Because the correlation functions are approximately spherically symmetric in the LCMS, the measured correlation functions are approximately independent of the orientation of $\boldsymbol{q}_{\mathrm{LCMS}}$.

We thus conclude that $Q$ can be introduced in a reasonable manner as the proper one-dimensional variable of the Bose-Einstein correlations in $\sqrt{s_{{}_{NN}}}=200$ GeV Au$+$Au collisions.

In order to perform a detailed shape analysis in the LCMS, we thus measured them as univariate functions of $Q$ (for $K_{T}$ values in various ranges). Thus this one-dimensional analysis in the LCMS in terms of $Q$ can be viewed as an approximation to a three-dimensional analysis with the approximation that the three HBT radii are equal.

In principle, a more complete picture of the source geometry can be obtained by a three-dimensional Lévy analysis, utilizing Eqs. (49)-(52) of Ref. Csörgő et al. (2004). Given that the details of these studies go beyond the scope of the current manuscript, let us make only some general remarks here. If the source is a symmetric three-dimensional Gaussian, then in a one-dimensional analysis (in our $Q$ variable, measured in the LCMS), one would obtain $\alpha=2$ for the Lévy shape parameter. If the source is an asymmetric 3D Gaussian, then non-Gaussian 1D correlation functions would be obtained, but also strong deviations from the Lévy shape could be observed. We investigated this using the method of Lévy expansion of the correlation functions Novák et al. (2016) for each $m_{T}$ bin, and found no first order deviations from the Lévy shape. However, an $m_{T}$ averaged correlation function shows deviations from the pure Lévy shape, which may be attributed to the $m_{T}$ dependence of $\alpha$. These observations suggest that the observed Lévy shapes do not originate from an asymmetric three-dimensional Gaussian source.

If the final-state-strong and Coulomb interactions can be neglected, then Eq. (14) implies that the correlation function takes the value $2$ at vanishing relative momentum, $C^{(0)}_{2}(Q=0,K)=2$. However, experimentally the two-track resolution (corresponding to a minimum value of $Q_{\mathrm{\rm min}}$ of at least 6–8 MeV, depending on track momentum) prevents the measurement correlation functions at $Q=0$. So the correlation function is measured at nonzero relative momenta and then extrapolated to $Q=0$. This extrapolated value in general can be different from the exact value at $Q=0$, and this can be quantified by defining

where $\lambda$ may depend on average momentum $K$.

In our analysis we measure the $C_{2}$ correlation functions as a ratio of actual and background distributions $A$ and $B$, and we have carefully checked in our dataset that $\lim_{Q\rightarrow 0}A(Q,K_{T})=0$ and $\lim_{Q\rightarrow 0}B(Q,K_{T})=0$ in every transverse momentum range, indicating that the split tracks have been removed from our data sample. The two-track resolution, embodied into the values of two-track cuts as seen in Section III.3, corresponds to a maximum spatial resolution of $R_{\mathrm{\rm max}}\approx\hbar/Q_{\mathrm{\rm min}}\approx 25-30$ fm. In our analysis, source details on spatial scales larger or equal to $R_{\rm max}$ cannot be experimentally resolved.

This (perhaps with different $R_{\mathrm{\rm max}}$ values) is a general feature of any similar experiment, and it leads to the core-halo picture of Bose-Einstein correlations in high energy heavy ion reactions Bolz et al. (1993); Csörgő et al. (1996). The core-halo picture treats the particle emitting source as a composite one, corresponding to particle emission from a hydrodynamically behaving fireball-type core, surrounded by a halo of long-lived resonances. Such a picture is particularly relevant for pion production. Several long-lived resonances with decay widths of $\Gamma\ll Q_{\mathrm{\rm min}}$ (like the $\eta$, $\eta^{\prime}$, $K^{0}_{S}$ mesons, and, depending on the experimental two-track resolution, maybe the $\omega$ meson) decay to pions that contribute to the halo region. The general structure of the core-halo model may hold not only for pion production but for the production of other mesons as well.

In short, $\lim_{Q\rightarrow 0}C_{2}(Q,K)=1+\lambda(K)$ is in general different from the exact value of $C_{2}(Q=0,K)$ which (independently of $K$) is 2 for a thermal, fully chaotic particle source. In most data sets, $\lambda<1$ holds, see again the overview papers in Refs. Boal et al. (1990); Weiner (2000); Wiedemann and Heinz (1999); Csörgő (2002); Lisa et al. (2005); Tannenbaum (2006); Lisa and Pratt (2008); A. Kisiel, for the ALICE Collaboration (2011); Heinz and Snellings (2013); Adamczyk et al. (2015).

In the core-halo picture, for thermal particle emission, the intercept $\lambda$, the extrapolation of the measured resolvable part of the correlation function to zero relative momentum, is the square of the fraction of pions coming from the core, defined as

because both pions have to come from the core if they are to contribute to the resolvable correlation function. This requires a physical assumption, that the phase-space density of the pion emitting source is made up of two components, i.e.

each component having a Fourier transform defined as

where we again used the four-vector variables $q=p_{1}-p_{2}$ and $K=(p_{1}+p_{2})/2$. Then each component has a space-time integral corresponding to the contribution of the given component to the momentum distribution. We then may define

Here the first equation in Eq. (34) and Eq. (35) represents our physical assumption about the phase-space density of the core and the halo, while the second equation in Eq. (34) and Eq. (35) indicates a mathematical identity about the Fourier transform. Taking these and Eq. (31) into account, we obtain

For the experimentally resolvable $q$ values, this system of physical assumptions yields the approximation

thus the correlation function ($C^{(0)}_{2}(q,K)$) shown in Eq. (14) can be expressed as

Hence, in the core-halo picture, at any given momentum

holds; see Ref. Csörgő et al. (1996) for details. Thus parameter $\lambda$ can be interpreted as the squared fraction of pions from the core with respect to the total number of pions with a given average momentum $K$. The $q$ dependent part in Eq. (38), i.e. the shape of the Bose-Einstein correlation function is connected to the core, $S_{\rm core}$. This source component is the one that may correspond to the perfect fluid, the hydrodynamically evolving central part of the fireball created in high energy heavy ion collisions.

If we assume that the source ($S$) is a sum of the core and the halo components as shown in Eq. (31), then it follows that the pair distribution ($D$) shown in Eq. (19), is a sum of the three components,

where subscript ‘c’ denotes the core and ‘h’ denotes the halo. It can be easily shown that the core-core component denoted by $\mathrm{(c,c)}$ is resolvable, but the core-halo or $\mathrm{(c,h)}$ type of pion pairs or the halo-halo or $\mathrm{(h,h)}$ components are unresolvable (i.e. the width of their Fourier transform is below the minimal resolvable momentum difference). With this compared to Eq. (20), the correlation function of Eq. (38) can be re-expressed as

In summary, $\lim_{q\rightarrow 0}C_{2}(q,K)\neq 2$ is an experimental finding, and so it is customary to introduce $\lambda$ as an experimental parameter, defined as $\lim_{q\rightarrow 0}C_{2}(q,K)-1$, and measured by extrapolating the correlation function to zero relative momentum. The core-halo model is then an interpretation of the value $\lambda$. It also relates the relative momentum dependent, resolvable part of the Bose-Einstein correlation function to $S_{\rm core}$, the core component of particle emission in high energy heavy ion collisions. From this interpretation it is particularly clear that while long-lived resonance effects dominate the variances of the source, they lead to a peak in the unresolvable part of the Bose-Einstein correlation function, with measurable effects only on $\lambda$. Particle emission from the hydrodynamically expanding fireball however, i.e. the core component of the source, is observable from the $q$-dependent shape analysis of the Bose-Einstein correlation functions.

Thus one of the motivations for measuring the $\lambda$ parameter is that it carries indirect information on the decays of long-lived resonances to the observable pion spectra. Of particular interest is the contribution of the $\eta^{\prime}$ meson to the low momentum pion yield. It is expected Kapusta et al. (1996) that in the case of chiral $U_{A}(1)$ symmetry restoration in heavy-ion collisions, the in-medium mass of the $\eta^{\prime}$ meson (the ninth pseudoscalar meson, a would-be Goldstone boson) is decreased, thus its production cross section is heavily enhanced at low momentum. This (because the decay chain of the $\eta^{\prime}$ meson produces many charged pions) implies that at low transverse momentum, the $\lambda$ parameter decreases Vance et al. (1998). A recent study Csörgő et al. (2010) of existing $\lambda(m_{T})$ measurements (presented in greater detail in Ref Vértesi et al. (2011)) reported an indirect observation of a mass drop of the $\eta^{\prime}$ meson in $\sqrt{s_{{}_{NN}}}=200$ GeV Au$+$Au collisions at RHIC.

However, many of the earlier $\lambda(m_{T})$ measurements were made with the assumption that the shape of the correlation function is a Gaussian one. Given the fact that the detailed analysis presented below indicates that the Gaussian approximation is a statistically unfavored assumption, we attempt here a precise shape analysis of the correlation functions. This is required for a precise measurement of the intercept parameter $\lambda$, as its value depends on the shape of the correlation function through the extrapolation of the measured correlation function to vanishing relative momentum.

Let us note here that the modification of the observable intercept parameter $\lambda$ from unity can result from various reasons besides the core-halo model, for example coherence in the pion production Weiner (2000); Csörgő (2002). If a fraction of pions are created in a coherent manner, then two- and three-particle Bose-Einstein correlation functions at zero relative momentum are simply related to the fraction of coherently produced pions and to the fraction of pions coming from the core Csörgő (2002). Thus a simultaneous measurement of $\lambda$ in two- and three-pion correlation functions offers the possibility of separating the component of a possibly coherent pion production, in addition to the resonance decay contribution. Such a simultaneous analysis of second, third and higher order correlations was recently reported at the LHC Adam et al. (2016). Also, more exotic quantum statistical effects like squeezed coherent states may modify the values of the intercept parameter (however, in the present analysis we have no compelling reason to consider this possibility). Hence, one of the goals of the paper is to measure $\lambda(m_{T})$ precisely, without any physical assumption about the mechanism of the pion production.

In the following, we utilize a generalization of the usual Gaussian shape of the Bose-Einstein correlations, namely we analyze our data using Lévy-stable source distributions. We have carefully tested that this source model is in agreement with our data in all the transverse momentum regions studied. All the Lévy fits were statistically acceptable, as discussed in Section VI. We note that using the method of Lévy expansion of the correlation functions Novák et al. (2016), we investigated deviations from the Lévy shape. We have found that the coefficient of the first correction term is within uncertainties consistent with zero. Hence, we restrict the presentation of our results to the analysis of the correlation functions in terms of Lévy-stable source distributions.

Past measurements of two-pion Bose-Einstein correlation functions in Au$+$Au collisions that went beyond the Gaussian approximation show that the precise shape of Bose-Einstein correlations is indeed not Gaussian Afanasiev et al. (2008); Adler et al. (2007). The shape exhibits a power-law-like long-range component. In expanding systems, a generalized form of the central limit theorem and investigation of generalized random walk (also called anomalous diffusion) suggests the appearance of Lévy distributions as source functions Metzler et al. (1999); Csörgő et al. (2004). The one-dimensional, symmetric Lévy distribution is the generalization of the Gaussian distribution defined by the Fourier transform

Here $R$ is called the Lévy length scale parameter, and $\alpha$ is called the Lévy index of stability. In the $\alpha=2$ case we recover a Gaussian form; in the $\alpha=1$ case, we have a Cauchy distribution. For $\alpha<2$, the Lévy distributions have a power-law-like tail, $\mathcal{L}(\alpha,R,\boldsymbol{r})\propto(r/R)^{-(3+\alpha)}$ for $r/R\to\infty$ (with $r\equiv|\boldsymbol{r}|$). Equivalently, for the angle-averaged Lévy distribution one gets

Thus Lévy distributions for $\alpha<2$ have an infinite second moment or root-mean-square (RMS) radius. However, even in this case, the scale parameter $R$ provides a measure of the characteristic size of the system. In particular, the integral of the Lévy distribution is finite and proportional to $R^{3}$. Note also that if the core part of the source ($S_{\mathrm{core}}$) has a Lévy shape, then the core-core pair distribution ($D_{\mathrm{(c,c)}}$) also has a Lévy shape, due to the fact that the autocorrelation of two identical Lévy distributions is also a Lévy distribution with the same index of stability $\alpha$,

Thus the Lévy-type source distributions offer a more general description of the shape of the correlation function than a Gaussian would do. They provide a better handle on the $\lambda$ intercept parameter as well. The Gaussian limit corresponds to the special $\alpha=2$ case, so one can experimentally check how far given data are from the Gaussian limit. We illustrate the shape of Lévy-type source distributions ($S_{\mathrm{core}}=\mathcal{L}(\alpha,R,\boldsymbol{r})$) with various $\alpha$ values in Fig. 2.

There is yet another motivation for Lévy distributions. Namely, the exponent $\alpha$ of the Lévy distribution (that determines the power-law-like behavior of the distribution at large distances) is related to the critical exponent $\eta$ of a system at a second order phase transition Csörgő (2008). This exponent characterizes the power-law structure of the spatial correlation at the critical point. If an order parameter $\phi$ is introduced, its correlation function (in three dimensions, as a function of distance $r$) will be

As noted above in Eq. (43), the Lévy source distribution has the same limiting behavior, thus in this case, $\eta=\alpha$. According to lattice quantum chromodynamics (QCD) Aoki et al. (2006); Bhattacharya et al. (2014); Soltz et al. (2015) the quark-hadron transition is analytic (cross-over) at vanishing baryochemical potential $\mu_{B}=0$, and is expected to be a first order phase transition at high values of $\mu_{B}$. There may be a critical endpoint (CEP) at certain intermediate values of $\mu_{B}$, where one has a second order phase transition, with a specific value of the $\eta$ exponent. This value is $0.03631(3)$ in the 3D Ising model El-Showk et al. (2014), and $0.50\pm 0.05$ in the random field 3D Ising model Rieger (1995). Given that the second order QCD phase transition is expected to be in the same universality class as the 3D Ising model Halasz et al. (1998); Stephanov et al. (1998), the QCD critical point may be signaled by Lévy sources with a specific $\alpha$ exponent. To locate and characterize the CEP is one of the most pressing present day challenges of experimental heavy-ion physics. It is thus desirable to measure $\alpha$ for various colliding systems and collision energies, to map various parts of the $(\mu_{B},T)$ plane, in a quest to find the location of the CEP of the quark-hadron transition. We present below the first determination of the Lévy index of stability in $\sqrt{s_{{}_{NN}}}=200$ GeV Au$+$Au collisions.

Using the plane-wave approximation, and assuming a spherically symmetric, three dimensional Lévy-type source and using the core-halo model, the shape of the two-particle correlation function turns out to have the simple form of

with $Q$ being the independent variable as introduced in Eq. (26), and with three fit parameters, which may depend on average momentum $K$. The scale parameter $R$, the strength (intercept) $\lambda$ and the Lévy index $\alpha$ (note that the fitting procedure is detailed in Section VI.1). However, one cannot fit the above functional form to the measured correlation functions before properly taking the final state Coulomb repulsion of the identically charged pions into account.

In the treatment of this effect, we follow the general lines of the Sinyukov-Bowler method Sinyukov et al. (1998); Bowler (1991). Coupling this with the core-halo picture, one has to average the modulus squared of the final state pair wave-function over the “core-core” spatial pair distribution $D_{\mathrm{(c,c)}}(\boldsymbol{r},K)$, obtaining