Coulomb and strong interactions in the final state of HBT correlations for L\'evy type source functions

We present detailed calculations about the expected shape of two-pion Bose-Einstein (or HBT) correlations in high energy heavy ion collisions that include a realistic treatment of final state Coulomb interaction as well as strong interactions (dominated by s-wave scattering). We assume L\'evy type source functions, a generalization that goes beyond the Gaussian approximation. Various recent experimental results found the use of such source types necessary to properly describe the shape of the measured correlation functions. We find that strong final state interactions may play an important role in the shape of the two-pion correlation functions, especially if one considers source parameters beyond the Gaussian HBT radii. Precise experimental determination of these source parameters (such as L\'evy stability exponent, correlation strength, etc.) seems to require the inclusion of the treatment of strong interaction not just for heavier particles (e.g. protons, lambdas) but also in case of two-pion measurements.


I. INTRODUCTION
Heavy ion physics strives to understand the properties of strongly interacting matter produced in high energy nuclear collisions.One of the key observables suited for the experimental investigation of the space-time geometry of such collision events is the femtoscopic correlation of produced particles (called Bose-Einstein correlations in case of identical bosons).Since the discovery of quantum statistical correlations of pions produced in high energy reactions [1,2], more and more experimental data led to a refined understanding of the connection between such correlations and the actual source dynamics, as well as an increased expectation on phenomenological models to reproduce the observations.In conjunction with the discovery of the strongly interacting Quark-Gluon Plasma (sQGP) by the experiments at the Relativistic Heavy Ion Collider [3][4][5][6] a renewed interest arose in the investigation of femtoscopic correlations.For a review of such measurements and connected phenomenological studies, see e.g.Refs.[7,8].
In heavy ion physics, for many years the usual assumption for the source shape was Gaussian.This was corroborated by phenomenological studies such as hydrodynamical model calculations (see e.g.Refs.[9,10]).Recent results showed that to achieve a statistically acceptable description of the measured correlation functions, one must go beyond this simple picture.The application of the source imaging technique discussed in Ref. [11] to correlation functions measured in high energy heavy ion collisions led to one of the first signs of non-Gaussian behavior in such reactions [12]; it was found that the two-pion source function indeed exhibits a power-law behavior.Since then a lot of experimental as well as theoretical work has been done in this direction.Recent results by the PHENIX experiment [13] showed that by utilizing Lévy type sources one can provide an acceptable description of the measured correlations.These type of source functions are expected to emerge from a scenario called anomalous diffusion [14], but there are other possible competing explanations such as jet fragmentation [15] or critical behavior [16].
When one tries to extract information about the source through the analysis of femtoscopic correlations, it is of utmost importance to properly take into account final state interactions (FSI).The shapes of the experimentally measured correlation functions are significantly affected by these interactions (such as Coulomb repulsion and also strong interaction), and taking them into account in the theoretical framework is crucial.The effect of the Coulomb interaction and the methods to properly include it in the description of the correlation functions have been widely studied before, for details see e.g.Refs.[17][18][19].However, final state strong interaction between like-sign pions is generally thought to have a small effect [20], so in most experimental analyses it is neglected.In this paper we present a detailed calculation of the shape of two-pion HBT correlation functions with the assumption of Lévy stable source functions taking into account Coulomb and strong final state interactions.
The structure of the paper is as follows: in Section II.we discuss the basic definitions and properties of the femtoscopic correlations with special emphasis on the choice of the source function.In Section III.we investigate the effect of final state interactions on the pair wave function, and subsequently on the correlation function.In Section IV. we present results of a numerical calculation of the correlation function and investigate the differences between using only Coulomb or both Coulomb and strong interactions.Finally, in Section V. we conclude and summarize our findings.

II. FEMTOSCOPIC CORRELATIONS
In this section we discuss the basic definitions and properties of femtoscopic correlations, with special emphasis on the shape of the source function.

A. Basic definitions
The general definition of the two-particle correlation function as a function of the single particle four-momenta is the following: 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.The pair momentum distribution can be calculated from the S(x, p) source distribution and the Ψ p1,p2 (x 1 , x 2 ) symmetrized pair wave function: (2) Using the pair source D(r, K), defined as equation ( 1) can be reinterpreted as This way, instead of the single-particle variables p 1 , p 2 , x 1 , x 2 one can use the following pair variables: the pair separation four-vector r, the pair center of mass fourvector ρ, the relative momentum k = (p 1 −p 2 )/2, and the average momentum K = (p 1 + p 2 )/2.Since the Lorentzproduct of the k and K four-vectors are zero, one may transform the k dependent correlation function to depend on the three-vector component k only.Furthermore, if the energy of the particles contributing to the correlation function are similar, then K is approximately on shell, so the correlation function can be measured as a function of k and K.At this point it is also useful to introduce the corehalo picture, in which the particle emitting source has two components: a hydrodynamically behaving fireballlike core which contains particles created directly from the freeze-out (or from decays of short-lived resonances), and a surrounding halo which contains particles that are the decay products of long-lived resonances (such as η, η ′ , K 0 S , ω).This picture is particularly important for pions, but the general structure of the model may be relevant for other mesons as well.If one assumes that the single-particle source has two components (S = S core + S halo ), it follows that the pair source D will have three -a core-core, a core-halo, and a halo-halo component: Experimentally however, only the core-core part is relevant, the width of the Fourier transform of the other two is below the minimal resolvable momentum difference.
Introducing the correlation strength parameter λ and coupling the core-halo model with the Bowler-Sinyukov procedure the correlation function can be written as More details about the core-halo model and the importance of the λ correlation strength parameter can be found e.g. in Ref. [13] .
To calculate the shape of the C 2 (k, K) two-particle correlation function, one needs an assumption on the shape of the pair source D (c,c) (r, K), and a proper description of the effect of final state interactions enclosed in the Ψ  we combine the previous calculations to derive the shape of the correlation function.

B. Lévy-stable source functions
Stable distributions are of utmost importance when studying the limiting distributions of random variables based on a sum of elementary processes.It is well known, that in case of one dimensional random variables, the stable distributions can be given through the following formula: where the characteristic function is given as: where Φ = tan( πα 2 ), α ̸ = 1, − 2 π log |q|, α = 1.
In our case, the symmetric, centered (β = 0, µ = 0) stable distributions may play a role of the source distribution, if that results from a statistical process.In multiple dimensions, the situation is far less clear.It is however known that the following distribution in N dimensions is stable [21]: from which in case of spherical symmetry (R ij = R 2 δ ij ), we obtain The two main parameters of such distributions are the index of stability, α, and the scale parameter, R. In case of α < 2 the distribution exhibits a power-law behavior, while the α = 2 case corresponds to the Gaussian distribution.The most important property of this distribution is that any moment greater than α is not defined and it retains the same α under convolution of random variables.From the latter it is apparent that if the single particle source S core (r) is a Lévy-stable distribution, then the pair-source D (c,c) (r) also has a Lévy shape with the same index of stability α: ) An illustration of the shape of such distributions can be seen on Fig. 1.The average momentum dependence appears through the two parameters of D (c,c) (r): The dependence of the Lévy source parameters on the pair average momentum K is non-trivial, and is often the subject of the experimental investigations.

III. FINAL STATE INTERACTIONS
To make the paper as self-contained as possible, in this section we review the methodology of the calculation of a correlation function that includes the effect of the final state Coulomb and strong interactions.In doing so, we closely follow along the lines of Ref. [18].

A. The pair wave function
Firstly let us introduce the Sommerfeld parameter η that appears frequently during calculations concerning the quantum mechanical Coulomb problem: Here µ is the reduced mass of the particle pair.Note that one often uses the fine structure constant α ≡ in the definition of η, with which it could be written as η = µα ℏck ; nevertheless, we avoid this in this paper because we denote the Lévy index also by α, as indicated in the previous section.
A normalization constant N appears in many contexts in the Coulomb wave function.Its definition is and its modulus square, which is called the Gamow factor, can be calculated with elementary functions (owing to the well known step and reflection properties of the gamma function) as The Schrödinger equation in a repulsive Coulomb potential can be written as For the treatment of the final state interactions, one has to utilize the scattering wave solutions whose asymptotic form is a plane wave plus a spherical wave.Such solutions for the Coulomb potential are well known: Here F(a, b, z) is the (renormalized) confluent hypergeometric function (Kummer's function); its definition and some basic properties are recited in Appendix A. (A well-known property shows that the two forms of each functions introduced here are indeed equal.) The connection between these wave functions is From the asymptotic expression of the confluent hypergeometric function one can verify that the asymptotic form of these wave functions is Here the notation f c (ϑ) stands for the Coulomb scattering amplitude, which is defined as One indeed sees that asymptotically the ψ (+) k (r) and the ψ (−) k (r) wave functions contain a plane wave plus an outgoing or an incoming spherical wave, respectively.(There are logarithmic factors stemming from the long range nature of the Coulomb interaction that distort both of them; these factors do not influence the physical meaning of the wave functions.)The ψ (+) k (r) and the ψ (−) k (r) functions are called in and out scattering states, respectively. 1 The scattering states written up here can be expanded in terms of energy eigenstates which are also angular momentum eigenstates.For given l and m angular momentum quantum numbers, one has two linearly independent angular momentum eigenstate solutions of the ( 16) Schrödinger equation: their angle dependence is that of the Y lm (ϑ, φ) spherical harmonic function, and their radial parts are called regular and singular Coulomb waves, respectively.We denote them here by F k,l (r) and G k,l (r) (as they depend on the k wave number magnitude and the l total angular momentum quantum number but not on the magnetic quantum number m); their expression is where the so-called Tricomi's function, U (a, b, z) is another solution of the confluent hypergeometric equation (see Appendix A for some details).They are chosen for the set of linearly independent solutions because F k,l is finite at the r=0 origin, and their asymptotic form is quite simple and straightforward: for r → ∞ we have where the so-called Coulomb phase shift δ c k,l is defined as 1 It is a known fact that when calculating transition matrix elements, one has to utilize the ψ (−) k (r) state (the out state) for the wave function of the final state; this might seem somewhat counter-intuitive, since this function contains an incoming spherical wave.Similarly, one has to use ψ (+) k (r) for the initial state.See e.g.Ref. [22] for some details.
One can also take a linear combination of these two functions whose asymptotic form contains an additional arbitrary ∆ k,l phase shift.(28) whose asymptotic form duly is The above scattering-like solutions of the Schrödinger equation can be expanded in partial waves as Owing to the short range of strong interaction, we can treat its effect by introducing the ∆ s k,0 s-wave ,,strong" phase shift, and modifying the s-wave component of the exact Coulomb wave function to a s-wave which contains this additional phase shift (see more details in e.g.Ref [23]).This is done by replacing the F k,0 function in the l=0 term in the expansion (30) with the above defined M s k,0 (r) function which contains the additional ∆ s k,0 phase shift: so the wave function incorporating the Coulomb and strong interaction effects, Ψ cs k (r), becomes Substituting the formulas for the respective wave functions encountered here, we get For identical bosonic particles (e.g.pions) one needs the symmetrized two-particle wave function: Finally, one needs to calculate the modulus square of the wave function.The [r → −r] term within the braces in the following expression represents terms similar to the ones that stand before it, just with a mirrored r: B. The two-particle correlation function In this section we combine the previously discussed approaches, and write up the complete functional form of the correlation function by plugging in Equation ( 12) and ( 35) to Equation (6).
where the I (c,c) (k) integral can be written as Substituting Eq. ( 35) into Eq.(37) we get the following expression: where the following integrals were introduced: The last step is to explore the dependence of the strong phase shift ∆ s k,0 on k.Using the notation of Ref. [18] we can relate ∆ s k,0 to the full (Coulomb+strong) scattering amplitude f c (k): The scattering amplitude f c (k) can be expressed as [23] where h(η) is related to the digamma function ψ as The k dependence of f c (k) partly comes from the function K(k), which can be expressed with the δ k,0 phaseshift (where the (2) superscript denotes the I = 2 isospin channel, the only allowed channel in case of identical charged pion pairs): If there would be no Coulomb, only strong interaction, δ k,0 would be identical to the previously introduced ∆ s k,0 strong phase-shift.One can find different parametrizations for δ k,0 in the literature, in the following we mention some of them.A simple parametrization can be found in J. Bijnens et al. [24]: where a (2) 0 is called the scattering length, and r (2) 0 is called the effective range.The latter can also be connected to a b (2) 0 slope parameter as This effective-range parametrization is thought to be useful when the scattering length is much larger than the range of the scattering potential [25], which is not the case for identical pion scattering.Another parametrization [26] better suited for our investigations can be written up with the help of the center-of-mass energy s = 4(m 2 π + k 2 ) as , where ( 49) b(2) The s (2) 0 parameter corresponds to the value of s where the phase shift passes through 90 • .It usually has a negative value, indicating that for the I = 2 channel the phase remains below 90 • .The parametrization can also be extended with higher order terms, the values of the parameters can be found e.g. in Colangelo-Gasser-Leutwyler (CGL) [27]: a where the parameter values are the following: z 2 = 143.5 MeV, B 0 = −79.4,B 1 = −63.0,√ ŝ = 1050 MeV.A comparison of the previously mentioned parametrizations can be seen on Fig. 2. In the k range important for our investigations (k ≲ 100 MeV/c) the different parametrizations give almost identical results, so in the following we utilized the most recent one from Ref. [28].

IV. NUMERICAL RESULTS
In this chapter we present the results of the numerical calculation of C 2 (k).Using numerical integral calculations we created a lookup table for the function defined in Equation (38) for a wide range of values of k, R and α.This lookup table then was used to obtain the value of the function for any k, R and α by interpolation (within the available range).
If we omit the I (3) and I (4) terms from Eq. ( 38), we get back the pure Coulomb part.In the following, we compare the correlation function containing only the Coulomb interaction with the one containing both the Coulomb and the strong interactions, and try to give an estimate on the change in the values of the Lévy source parameters that is caused by the proper treatment of the strong interaction compared to the neglection of it.
From here on, we change the relative momentum variable to Q = 2k to better compare to the notation of published experimental results.

A. Comparison of Coulomb and strong FSI effects
Fig. 3. shows the calculated correlation functions for three different Lévy-scale values at the same index of stability α and same correlation strength λ.It is clearly visible that turning on the strong interaction affects the strength of the correlation functions, however, the effect on the Lévy-scale R and the index of stability α is not so transparent at this point.
To investigate the effect of the strong interaction in more detail, we generated histograms by sampling the calculated functions containing both Coulomb and strong interactions.To make the generated correlation function resemble real data, we randomly scatter the points around the calculated function and assign a relative error proportional to 1/Q (which is a realistic assumption if one considers typical experimental scenarios).We then fit the generated data with the help of the ROOT Minuit2 minimizer framework, with a similar method to what is described in Ref. [13].To check the validity of the fitting method, first we fit the generated histogram with the corresponding functional form to see if we get back the input parameter values.Fig. 4 shows such a fit to the generated data.The fit converged with an acceptable χ 2 /NDF value, the error matrix turned out to be accurate, and for the output parameter we got back within errors the same ones as were given as input.We repeated this test for multiple different input parameter values and found that our fitting method is indeed reliable.
As a next step, we took the same generated data and fitted it with a function containing only the effect of the Coulomb interaction.Fig. 5 shows an example for such a fit on panel (a).The fit converged again, the error matrix again turned out to be accurate.The resulting χ 2 value becomes just slightly higher than before, nevertheless, the fit is still acceptable.Although in this case the function containing only the Coulomb interaction can provide an acceptable fit to the generated data which contains also the strong interaction, the values of the fit parameters differ from the input parameter values.It seems that in this case one underestimates the value of λ from such a fit, and overestimates α.Within this precision, it seems that the value of R is unaffected.
One can also assume that if the data is more precise, meaning that the fluctuation and the statistical uncertainty of the generated points are smaller, the fit will not provide an acceptable χ 2 anymore.To check this, we also generated such C 2 (Q) histrograms, and found that the Coulomb fits converged, but indeed the χ 2 values increase by a considerable amount resulting in statistically unacceptable fits.An example for this can be seen on panel (b) of Fig. 5.One can also observe that on the subplot showing the values of the difference of the fit from the data divided by the uncertainty of the datapoint, a characteristic oscillating structure appears.

B. Quantitative estimation of the strong FSI effect
To give a better estimation on the change in the parameter values when fitting data containing strong interaction with a function containing only the Coulomb effect, we generated and fitted histograms similar to panel (b) of Fig. 5, spanning a wide range in parameter space of λ input = 0.3 − 1.0, R input = 3 fm -9 fm and α input = 1.0 − 2.0.For each fit parameter, we plotted the output versus the input values.The plotted output values represent a weighted average of output values coming from the same input for the given parameter but different inputs for the other two parameters.The results of this investigation can be seen on Fig. 6, panel (a)-(c).
By fitting data containing the Coulomb and strong final state interactions with a functional form describing only the Coulomb part, it seems that the correlation strength λ is underestimated by about 5% on average.The effect on the Lévy-scale parameter R is negligible at small values of it, while at higher values of R (up to about 9 fm) it is also slightly underestimated, by about 1%.The Lévy exponent α is overestimated by about 1-2%.
The estimations given here for the change in parameter values are by no means universal, they also depend on other factors such as numerical precision of the integral calculations, fit limits (Q min dependence), the precision of the generated data (see for example the difference between Fig. 5 (a) and (b)), or the parametrization of the strong phase-shift.The important conclusion from our investigations is that if the data is precise enough (which could be the case for recent measurements at RHIC or LHC), one most likely has to incorporate the strong interaction in the fits to achieve a statistically acceptable description of pion-pion correlation functions.

V. SUMMARY AND CONCLUSIONS
In this paper we presented a detailed calculation of the shape of two-pion HBT correlation functions with the assumption of Lévy stable source functions, and taking into account the Coulomb and strong final state interactions.Strong final state interactions were treated in the s-wave approximation.
A numerical calculation of the correlation function revealed that the strong final state interaction can have a non-negligible effect on the shape of pion-pion correlation function.As a first step towards the more thorough evaluation, we presented a quantitative estimation of the magnitude of this effect.As a general trend, we can ascertain that fits without the strong interaction effect typically underestimate the strength of the correlation, λ, and the Lévy scale R, while overestimate the Lévy exponent α.The magnitudes of these deviations are generally found to be no more than a few percent.
However, typical fits to measured correlation functions can become statistically unacceptable if the strong inter- FIG.5: Numerically generated two-pion correlation histogram incorporating Coulomb and strong final state interactions, fitted with a functional form containing only the Coulomb effect.When the generated data is less precise (a), the fit is statistically acceptable, but the output parameter values differ from the input.The difference is even more pronounced when the generated data is more precise (b), in this example the value of λ decreased by about 4%, the value of R decreased by about 1%, and the value of α increased by about 3%.It is also important to note that in this case the χ 2 /NDF value is not acceptable anymore.action is neglected.If one aims at a high level of precision (feasible in case of precise enough data coming from today's typical heavy ion experiments), one can arrive at refined conclusions about the source function if the small deviations (caused by the strong interaction) are treated properly in the fitting procedure.
As an outlook, we note that there is some room for improvement in the methodology of the numerical calculations presented here.Such improvements might yield so precise predictions that it becomes possible to actually give constraints on like-sign pion strong interactions (i.e.scattering lengths) based on HBT correlation measurements in heavy ion collisions, a topic long thought to be interesting to investigate [29].We look forward to a concrete experimental test of the predictions made here about the shape of the correlation function that gets influenced by strong final state interaction.

( 2 )
k (r) pair wave function.In the following, in Section II B. we discuss the details of Lévy type source functions, and in Section III A. we proceed by discussing the calculation of Ψ

( 2 )
k (r) with the Coulomb and strong final state interactions included.Finally, in Section III B.
can be found in a more recent paper from García-Martín et al. (GM)[28]:

FIG. 3 :
FIG. 3: Two-pion correlation functions calculated for Lévystable sources.Three different Lévy-scale values are compared at the same index of stability α = 1.5 and same correlation strength λ = 1.The functions containing only the Coulomb interaction and the ones including both the Coulomb and strong interactions are shown separately.

FIG. 4 :
FIG.4: Numerically generated two-pion correlation histogram, fitted with the corresponding functional to test the validity of the fitting method.The output parameter values are within errors the same as the input.

FIG. 6 :
FIG. 6: Output versus input values from fits similar to Fig. 5.(b).The correlation strength λ is shown on panel (a), the Lévy scale parameter R is shown on panel (b) and the Lévy exponent α is shown on panel (c).The identity line is shown with a dashed line, while a linear fit is shown with a continuous line.For a given input parameter, the weighted average of the output values are shown with markers, and the standard deviation is shown with a band.