Event-by-Event Investigation of the Two-Particle Source Function in Heavy-Ion Collisions with EPOS

Exploring the shape of the pair-source function for particles such as pions or kaons has been an important goal of heavy-ion physics, and substantial effort has been made in order to understand the underlying physics behind the experimental observations of non-Gaussian behavior. In experiments, since no direct measurement of the source function is possible, quantum-statistical momentum correlations are utilized to gain information about the space-time geometry of the particle emitting source. Event generators, such as EPOS, however, provide direct access to the freeze-out coordinates of final state particles, and thus the source function can be constructed and investigated. The EPOS model is a sophisticated hybrid model where the initial stage evolution of the system is governed by Parton-Based Gribov–Regge theory, and subsequently a hydrodynamic evolution is utilized, followed by hadronization and hadron dynamics. EPOS has already proven to be successful in describing several different experimental observations for systems characterized by baryon chemical potential close to zero, but so far the source shape has not been explored in detail. In this paper we discuss an event-by-event analysis of the two-particle source function in sNN = 200 GeV Au+Au collisions generated by the EPOS model. We find that when utilizing all stages of the model, Lévy-shaped distributions (unlike Gaussian distributions) provide a good description of the source shape in the individual events. Hence it is clear that it is not the event averaging that creates the non-Gaussian features in the pair distributions. Based on this observation, we determine Lévy-parameters of the source as a function of event centrality and particle momentum.


Introduction
A long-standing goal of high-energy nuclear physics has been to understand and explore the space-time geometry of the particle emitting source created in heavy-ion collisions [1]. One main observable that is of great interest is the two-particle source function, sometimes also called spatial correlation function or pair-separation distribution. Detailed investigation of this quantity is important for a multitude of reasons (connected to hydrodynamic expansion, critical behavior, etc.), however, it is not something that is easy to reconstruct experimentally. There is a whole sub-field of high-energy nuclear-and particle-physics called femtoscopy, which deals with such measurements of lengths and time intervals on the femtometer (fm) scale [2]. Since it was shown by G. Goldhaber et al. that intensity correlations of identical pions can be used to gain information about the pair-source function [3,4], femtoscopy has propelled to the forefront of investigations, and today it is still one of the most extensively studied field of high-energy physics. Besides the ample experimental studies, phenomenological investigations also placed great emphasis on describing the shape of the source function. Hydrodynamical model calculations suggest [5][6][7][8][9][10] that the source-shape is Gaussian, and this was adopted by many measurements as well [11,12].
Source imaging studies [13,14] on the other hand suggest that the two-particle source function of pions has a long-range component, obeying a power-law behavior. It was also shown recently by various experimental measurements, that a generalization of the Gaussian source shape, the Lévy distribution can provide a much more suitable description of the observed sources [15,16]. These kind of source shapes arise in many different scenarios [17] such as anomalous diffusion [18], jet fragmentation [19], critical behavior [20], or resonance decays. It has been shown that even averaging over many events with different sizes can contribute to the appearance of a power-law component [21,22]. In order to have a better understanding of the underlying processes behind the experimental results, more effort is needed from the phenomenology side. Among the important tools for such investigations are the event-generators that encompass different theoretical and phenomenological methods to model nuclear reactions. One of such event generators is the EPOS model [23]-the Energy conserving quantum mechanical multiple scattering approach, based on Partons (parton ladders), Off-shell remnants, and Splitting of parton ladders. In this paper we present a detailed event-by-event analysis of the two-pion source distribution in √ s NN = 200 GeV Au+Au collisions generated by EPOS. The event-by-event nature of our analysis helps in deciding if the role of event averaging is crucial in the apparent non-Gaussian but Lévy nature of the observed sources. The paper is structured as follows. In Section 2 we introduce the EPOS model and its several stages of evolution. In Section 3 we discuss the basic definitions and properties of the pair source distribution. In Section 4 we discuss the details of the event-by-event analysis and the applied methods. In Section 5 we present the results of the analysis and discuss our findings. Finally, in Section 6 we summarize and conclude.

The EPOS Model
EPOS, Energy conserving quantum mechanical multiple scattering approach, based on Partons (parton ladders), Off-shell remnants, and Saturation of parton ladders, is a phenomenological model based on Monte Carlo techniques. It opens up the possibility of investigating various phenomena and observables such as particle production, momentum distributions or flow correlations, providing a better understanding of the evolution of the system created in elementary (proton-proton) collisions and also during complex reactions involving heavy-ions. The theoretical framework included in the model provides a coherent description of the space-time expansion of matter based on a precise spectrum of studies of both elementary processes such as electron-positron annihilation or lepton-nucleon scattering and more compound collisions of protons or nuclei. The model was designed to describe processes appearing in collisions at µ B ≈ 0, at very high (top RHIC or LHC) energies and for various systems, such as Au+Au, Pb+Pb or p+p.
In this Section, all of these stages are described.

Initial Stage of the Evolution
In the theoretical framework of the model the crucial element is the sophisticated treatment of both the hadron-hadron scattering and the initial stage of the collisions at ultrarelativistic energies. It is highly relevant in the understanding of possible parton-hadron phase transitions. In EPOS, a merged approach of the Gribov-Regge Theory (GRT) and the eikonalised parton model is utilised to provide proper treatment of the first interactions happening just after a collision. This approach satisfies conservation laws, and treats the subsequent Pomerons (interactions) equally (as opposed to other multiple interaction approaches, for example, Pythia, where the first interaction is not treated exactly the same way as the others) [24].
The formalism used for the calculation of the cross-sections is the same as the one used for calculating particle production. It is based on the Feynman diagrams of the QCD-inspired effective field theory and provides energy conservation. The nucleons are divided into a certain number of "constituents" carrying the incident momentum fraction. The fractions sum to unity in order to ensure momentum conservation. A nucleon is called a spectator if it is not part of the interaction region of the colliding nuclei. If a nucleon is not a spectator, then its constituents can either be participants taking part in the elementary interactions with constituents from the opposite side, or a remnant, which although part of the interaction region, does not take part in the elementary interactions. This is illustrated in Figure 1. The splitting between participants and remnants shows the momentum sharing between constituents ensuring conservation of the given variable [25].
The particle production is based on the String Model approach [26,27]. The parton ladders are recognized as a quasi-longitudinal color field (elementary flux tubes) and are treated as classical strings [23]. The intermediate gluons introduce the transverse motion into the kinky string evolution. The schematic picture of the flux tube with the transverse kink is shown on Figure 2.

Core-Corona Approach
If the density of the strings is very high, they cannot decay independently. This situation is characteristic for heavy-ion collisions and high-multiplicity pp collisions. Henceforth, in EPOS a dynamical process of the division of the strings segments into core and corona is introduced [23,28,29].
The core-corona division is based on the abilities of a given string segment to leave the "bulk matter". The transverse momentum of the element and the local string density are considered as criteria for the division. If the string segment belongs to a very dense area, it will not escape but will contribute to the core, which will be governed in the next step by the hydrodynamical evolution. When the segment originates from the part of the string close to the kink, it is characterized by high transverse momenta. It escapes the bulk matter and will contribute to the corona. It consequently will show up as a hadron in a jet. There is also a possibility that the string segment is close to the surface of the dense part of the medium, and its momentum is high enough to leave it; then it also becomes a corona particle.

Viscous Hydrodynamical Evolution, Event-by-Event Treatment and EoS
In EPOS a 3D+1 viscous hydrodynamics approach is applied, called vHLLE (viscous relativistic Harten-Lax-van Leer-Einfeld Riemann solver-based algorithm). In the simulations, the separate treatment of individual events is highly important-smooth initial conditions for all events are not applied. The event-by-event approach in hydrodynamical evolution is based on the random flux tube initial conditions [23]. It has a relevant impact on the final observables such as spectra or various harmonics of flow. The viscous hydrodynamics uses Equation of State X3F ("cross-over" and "3 flavor conservation") which is compatible with lattice QCD data from Ref. [30]. It corresponds to µ B = 0 MeV, hence this feature limits the applicability of the model to describe the region of the QCD phase diagram characterized by finite baryon density [23].

Hadronization and Hadronic Cascades
The expanding medium in the processes of hydrodynamical evolution reaching the given freeze-out condition is transformed into the particle spectra. In EPOS 3, the criterion characterizing the hadronization hypersurface is that the energy density equals 0.57 GeV/fm 3 . EPOS 3 furthermore utilizes the Cooper-Frye formula [31] when determining distributions. The final part of the simulation uses a so-called hadronic afterburner, based on UrQMD [32,33]. The hadronic scattering has a significant impact on the final observables [34].

the Two-Particle Source Function
In this section we discuss the basics definitions and properties of the two-particle source function. The pair source distribution is defined as the auto-correlation of the single particle phase-space density S(x, p): where instead of the single particle variables x and p the pair-variables appear -the pair center of mass four-vector ρ = (x 1 + x 2 )/2, the pair separation four-vector r = x 1 − x 2 , and the average momentum K = (p 1 + p 2 )/2. The D(r, K) distribution is the quantity that can be reconstructed indirectly from femtoscopic momentum correlation measurements, and experiments usually investigate the source-parameters that describe the shape of this distribution, see details, for example, in [15]. It was recently shown by different experiments that for pions this pair-source exhibits a power-law behavior, and can be described with a Lévy-stable distribution [15,16]. In case of spherical symmetry, the symmetric, centered stable distribution is defined as where the temporal dimension is removed from the dependence, made possible by the mass-shell condition, as detailed in Ref.
[15]. The two important parameters that describe such a distribution are the Lévy-scale parameter R and the Lévy-exponent α. One can think of the latter as the parameter that is responsible for "how far" the distribution is from the Gaussian. In the α = 2 case L(r; α, R) is identical to a Gaussian distribution, while in case of α < 2 it exhibits a power-law behavior. An illustration of the shape of such distributions can be seen on Figure 3. Since this distribution retains the same α exponent under convolution of random variables, if the single-particle source densities have a Lévy-shape then it follows that the two-particle source will also have such a shape, only the scale-parameter will be different. This can be summarized via the three-dimensional single-particle source S(x) and the pair-source D(r) (where we now suppress the momentum dependence, which is in turn contained in the parameters of the distribution):

Analysis Method
For the analysis presented below we used √ s NN = 200 GeV Au+Au events generated by EPOS359. Using like-sign pion pairs, we measured the one-dimensional pair-source distribution in the longitudinal co-moving system (LCMS). The LCMS pair-separation vector can be expressed in lab-frame single-particle coordinates as Using this variable, one can construct the spatio-temporal distance distribution D(r LCMS , t). After angle-and time-integration we obtain the one-dimensional distance distribution as: where we now suppress the K dependence indicated in Equation (1). Note that the dependence on the lab-frame time-coordinate disappears after the time integral of Equation (5), since we only keep the dependence on r LCMS , the final variable. In fact when analyzing the EPOS output, we only calculate the number of pairs in a given r LCMS bin, hence dependence on all other coordinates is naturally integrated out. When selecting pions we used the single-particle rapidity and transverse momentum requirements of |η| < 1 and 0.2 GeV/c < p T < 1.0 GeV/c. For each individual event we constructed the D(r LCMS ) distribution for 5 different average transverse momentum k T classes in equal bins ranging from 0.2 to 0.4 GeV/c. Note that k T = 0.5 K 2 x + K 2 y is the transverse component of K used in Equation (1). We chose this k T region to be around the peak of the pair k T distribution to have adequate statistics (number of pairs) in the individual k T bins. To investigate centrality dependence as well, we separated the measurements to the centrality classes of 0-5%, 5-10%, 10-20%, 20-30%. In total we used 63,000 EPOS events, 10,500 for the first two centrality classes and 21,000 for the rest.
As mentioned before, EPOS has different stages of evolution including hydrodynamic expansion and hadronic rescattering. In order to identify the effect of the different stages on the shape of the pair-source distribution, as well as the contribution from the resonance decay products, we separated our investigation to four different cases as follows: where primordial pions include pions coming from the thermal medium, that is, primordial pions are those that are not decay products. For each single event we fitted a Lévy distribution to the constructed D(r LCMS ) distribution within the range of 2 fm to 100 fm.
The fit was considered good if the confidence level calculated from the χ 2 and NDF values was greater than 0.1%. An example of such fits for the four different cases can be seen on   In case (a) we found that the events exhibit sharp cutoff features and mostly can be fitted well with a Gaussian. In case (b) the inclusion of the decay product pions results in power-law like structures, appearing at different regions in r LCMS . The shape of the events can be very different depending on the number (and origin) of decay pions in the sample. Due to the increased fluctuations and different event shapes the event-by-event extraction of the source parameters could not be done in the previous two cases. The fit settings would have to be fine-tuned for each event separately for this to work, which makes it impossible to do for thousands of events. In case (c) and (d) however, distinct non-Gaussian structures (power-law tails) are present in all events, shapes can be described by Lévy distributions in a statistically acceptable manner, and the extraction of the event-by-event source parameters is feasible. Furthermore, their distribution in the event sample can also be determined.
An example for such a distribution can be seen on Figure 5. In this example the distribution was reconstructed from fitting 21,000 events, out of which 18,460 fits were successful for case (c) and 18,768 for case (d) according to the confidence level criteria. Note that the distribution of source parameters was quite the same for non-acceptable fits as well; however, those do not necessarily represent the acquired source distributions, hence they were omitted from the further calculations. As mentioned before, we measured these twodimensional R vs. α distributions for four different centrality classes, and five different k T regions. From these we can extract the mean and standard deviation values, and investigate their centrality and k T dependence. There are multiple ways to determine the mean and standard deviation parameters-on one hand, one could do normal distribution fits to the obtained 2D histograms, and on the other hand, one can simply calculate the first and second momenta of the distributions. We checked both and since the results were quite similar, for the sake of simplicity we chose the latter one. Let us reiterate the point here that we analyzed individual EPOS events, and determined the R and α parameters of the pair source distribution in those individual events. We saw Lévy-shaped distributions when we included hadronic scattering, with slightly different α parameter values when decay pions were also included (besides primordial pions).

Results and Discussion
As described above, we performed fits to individual EPOS events from the final stage of EPOS (CORE+CORONA+UrQMD), and investigated averages of the resulting R and α parameters. We repeated this exercise for various centralities (0-5%, 5-10%, 10-20% and 20-30%) and k T regions (five equal bins from 0.2 to 0.4 GeV/c). We analyzed two cases separately: first the case of using only primordial pions, and then a case where both primordial and decay pions were included in the sample. Results for the mean R and α values versus m T = m 2 + k 2 T are shown in Figure 6.  Table A1 for both cases.
One can observe a clear decreasing trend in R with both m T and centrality, for both the case of primordial pions, as well as the case where decays were also included. This trend with m T is similar to the observed R 2 Gauss ∝ m T trend observed universally across collision centrality, particle type, colliding energy, and colliding system size [11,35], even though it is based on Gaussian source radii. The decrease of R with increasing centrality shows the relation of the Lévy-scale to the initial fireball size. One can also observe that R is only weakly affected by the inclusion of decay pions; the values are slightly higher in the latter case.
The Lévy-stability index α shows less prominent centrality dependence, although a small decrease for more peripheral events is visible in case we include decay pions as well. This feature, that is, the centrality dependence of α, was not yet investigated in experimental publications (with the exception of preliminary data in Ref. [36]). One may also observe a weak decrease with m T similarly to what was observed in Ref. [15]. Furthermore it is clearly visible that when decay pions are also included, the α parameter decreases. This is expected as decay pions produce an even stronger tail, creating a smaller α value.
Concrete values of R can be compared to measured values from Ref. [15]. There R values of 7-8 fm were measured for the 0-30% centrality class and the m T = 0.25-0.45 GeV/c 2 window. Our calculation yields similar values in the 10-20% and 20-30% centrality class. In Ref.
[15], however, α values around 1.2 with a weak m T dependent decrease were found. The trend in m T is similar in our analysis, but the magnitude of the α values is somewhat different. The reasons for this could be multifold, they can range from the unavoidable event averaging present in the experiment to initial or final state effects not present in our simulations. The exploration of this difference is beyond the scope of present paper.
Note that the filled bands on both plots indicate the standard deviation of the R and α distributions over the investigated event sample (let us remind the reader that we investigate and fit pair source distributions in individual events). The statistical uncertainty of these data points is basically negligible, due to the large number of events this average was performed over. This means that the trends observed in Figure 6 and discussed above are true features of the EPOS event sample investigated in this paper.

Summary and Conclusions
We investigated a sample of EPOS events individually, and using identical pion pairs we reconstructed the pair source function in every individual event. In the case when only primordial pions were analyzed before hadronic scattering, a Gaussian shape was observed. However, when decay pions were also included, already power-law like structures appeared, and after the inclusion of hadronic scatterings (via UrQMD), Lévyshaped pair distributions arose in the individual events. It is hence clear that it is not the event averaging that creates the non-Gaussian features in the pair distributions (and the arising correlation functions in femtoscopical measurements).
Subsequently, using the final stage of EPOS events (CORE+CORONA+UrQMD), we analyzed the event sample mean of the event-by-event Lévy-scale R and Lévy-index α values for the case of using only primordial pions and the case of including both primordial and decay pions. We observed clear trends as a function of m T and centrality for both cases. These observations show that in a realistic hydrodynamics-based simulation deviations from the Gaussian source shape appear when one includes hadronic scattering and decays. The values and trends of R are compatible with experimentally measured values, although we did not perform a detailed data comparison here. The weak decrease of α with m T is similar to what was observed experimentally, but the values of α are somewhat larger than measured values.
In the future we plan to utilize similar techniques to explore the dependence of these results on particle species as well. Further investigations might also include expanding the analysis to multiple dimensions, different collision energies, as well as reconstructing femtoscopical correlation functions.
Finally let us note that one of the important conclusions of our analysis is that the Lévy-shaped source assumption provides an acceptable description of the pion pair-source in EPOS3.

Appendix A
The data from Figure 6 is listed in Table A1.