A novel method for calculating Bose–Einstein correlation functions with Coulomb final-state interaction

Measurements of Bose–Einstein correlations played a crucial role in the discovery and the subsequent detailed exploration of the Quark-Gluon-Plasma (QGP) created in high-energy collisions of heavy nuclei. Such measurements gave rise to femtoscopy, a flourishing sub-field of high-energy physics, and there are important new directions to explore and discoveries to be made in the near future. One of these important current topics is the precise investigation of the shape of the correlation functions utilizing Lévy-stable sources. In this paper, we utilize Lebesgue’s dominated convergence theorem and Fubini’s theorem to present a novel method of calculating the shape of the two-particle correlation functions, including the Coulomb final-state interaction. This method relies on an exact calculation of a large part of the necessary integrals of the Coulomb wave function and can be utilized to calculate the correlation function for any source function with an easily accessible Fourier transform. In this way, it is eminently applicable to Lévy-stable source functions. In this particular case, we find that the new method is more robust and allows the investigation of a wider parameter range than the previously utilized techniques. We present an easily applicable software package that is ready to use in experimental studies.


Introduction
High energy heavy-ion physics aims at recreating the conditions characteristic to the first few microseconds of our Universe and studying the matter that is present under such conditions: the strongly interacting Quark Gluon Plasma (QGP) [1][2][3][4].Experiments observe the particles produced in heavy ion collisions in order to draw conclusions about the properties of this matter.An important class of observables is that of momentum correlations of identical particles.In the case of bosonic particles (e.g., pions or kaons), these are called Bose-Einstein-correlations because they arise as a result of the quantum statistical properties of the particles.They are also called HBT correlations after R. Hanbury Brown and R. Q. Twiss [5], who first discovered a related effect in astronomy.These correlations are connected to the space-time geometry of the particle emitting source [6]; thus, they provide a vital tool to study the matter produced in these collisions.These types of measurements have been extensively studied for more than 60 years now [7], especially since the era of the Relativistic Heavy Ion Collider (RHIC) started around the early 2000s [8,9].
The usual way to experimentally study the source function is to assume a source shape (e.g., a Gaussian distribution) characterized by a handful of parameters, calculate the correlation function arising from such a source, then test this assumption on the measured momentum correlation function.If the calculated correlation shape fits the data well (in terms of statistical acceptability), one can extract the source parameters and investigate their dependence on average momentum, centrality, collision energy, etc.One of the most essential steps in such an analysis is thus the calculation of the correlation function for a given source function.
As detailed below, if one neglects the final state interactions of the produced particles, the correlation function is given by the Fourier transform of the source function, a relatively simple calculational step in most cases.However, in the case of charged particles (such as charged pions, the most frequent target of measurements), the correlation function is given by a more complicated integral transform of the source function.In the process of fitting the source parameters to measurements, a frequently used method has been to calculate the Coulomb integral only for a pre-defined set of the parameter values and derive a "Coulomb correction factor" from it.In this way, however, there is an inherent distortion in the fitting procedure since the Coulomb effect is calculated at some arbitrary parameter values instead of the supposedly final, fitted ones.As experimental measurements and phenomenological descriptions become more and more refined, it becomes absolutely necessary that in the fitting procedure, the Coulomb interaction is considered in a self-consistent way.This necessitates the calculation of the mentioned Coulomb transform of the assumed source function for arbitrary source parameters during the fitting procedure.
The simplest assumption for the source function is that of a Gaussian distribution, and in the past decades, many aspects of Gaussian correlation measurements were explored [10][11][12].However, as the experimental resolution and the accumulated data increased throughout the years, an expectation arose to better describe the correlation function's shape [13].One possibility is presented by a spherical harmonics expansion [14]; another is to use a functional form for the source that goes beyond the Gaussian approximation.To that end, Lévy-stable distributions were utilized [15][16][17].In this way, a statistically acceptable successful description of the data was achieved.The utilization of more and more general source function assumptions, such as the mentioned Lévy distributions, is thus another pressing motivation for developing precise calculational tools to treat the Coulomb interaction.
In this paper, we present a new method for the calculation of the Bose-Einstein correlation function with the Coulomb interaction taken into account that is applicable to a wide variety of source functions.The new technique relies on the Fourier transform of the source function and is superior to already known methods in several ways: computationally, it is much cheaper and more stable than existing methods, and conceptually it highlights the actual effect of the Coulomb interaction on the correlation function in a straightforward way.In this paper, we present the method only for spherically symmetric source functions; however, the generalization to the non-spherical case is also being worked out.Another plausible future generalization is taking more complicated final state interaction schemes, such as strong interaction, into account, along the lines of Ref. [18].
The structure of this paper is as follows.In Section 2, we detail the basic definitions and known methodology pertaining to Bose-Einstein correlations and the effect of the Coulomb interaction, paying particular attention to the role of Lévy distributions as source functions.In Section 3, we present the new method.It has some delicate mathematical intricacies that are interesting on their own (so apart from the practical benefits of our new method, its derivation might interest the purely mathematically inclined reader as well).Some mathematical details are highlighted in Section 3, and some are left to the Appendix.In Section 4, we verify our new method for the case of Lévy-stable source functions.We find that the technique is significantly helpful and augments experimental measurements by making the correlation function fitting procedure easier.We also present a software package ready to use for such measurements [19].Section 5 summarizes our findings, pointing to possible future applications and generalizations.

Bose-Einstein correlation functions and the Coulomb interaction
The definition of the observable two-particle correlation function in a heavy-ion collision experiment is as follows: where p 1 and p 2 are the particle momentum (three-)vectors, N 1 (p) ≡ E dn d 3 p is the single-particle invariant momentum distribution (meaning the number of particles produced within a given invariant volume element in momentum space, and E is the particle energy corresponding to momentum p), and is the two-particle invariant momentum distribution, i.e., the number of particle pairs produced with momenta p 1 and p 2 (and so with energies E 1 and E 2 ).
Before proceeding, let us introduce some appropriate notations and variables.Three-vectors are denoted by boldface letters, four-vectors by standard typeset letters; e.g., p 1 is the three-momentum, p 1 ≡ (E 1 , p 1 ) is the four-momentum of particle 1. Lorentz products of four-vectors are denoted by a dot; for the metric tensor we use the g µν = diag(1, −1, −1, −1) convention.The scalar product (Descartes inner product) of three-vectors is denoted by writing the vectors next to each other.Below we deal with functions of two sets of space-time coordinates, x 1 ≡ (t 1 , r 1 ) and x 2 ≡ (t 2 , r 2 ); instead of these coordinates we can use the center-of-mass coordinate X ≡ (T, R) and the relative coordinates x ≡ (t, r) defined as (2) Similarly, besides the momenta of the particles, p 1 and p 2 , we use the average momentum four-vector K, as well as the relative momentum four-vector Q, and use the three-vector slices of these, K and Q: Since we are dealing with identical particles, their masses are equal, m 1 = m 2 ≡ m.With respect to the momentum variables, we usually suppress the ℏ reduced Planck constant, meaning that momentum vectors are also understood as wave number vectors if the dimensions of the quantities imply so; for example, we write a plane wave-like wave function of the r variable like e iQr but measure Q in MeV/c like a momentum variable.The C 2 (p 1 , p 2 ) correlation function introduced in Eq. ( 1) above equals unity if the particle production is uncorrelated.In high-energy nuclear collisions, there are a number of processes that lead to C 2 ̸ = 1, i.e., correlated particle production (collective flow, jets, resonance decays, etc.).For identical bosons (such as charged pions), the main source of correlation at low momentum difference is what is called Bose-Einstein correlation, a quantum mechanical correlation stemming from the indistinguishability of the particles.For the simplest theoretical treatment of this, one can introduce the source function S(x, p) ≡ S(t, r, p) as the probability of a particle produced at space-time point x ≡ (t, r) with momentum p.With this, the invariant momentum distribution is readily expressed as where d 4 x stands for dt d 3 r, i.e., the integration is taken over the whole space-time; however, usually S has a finite support or at least a rapid decrease in space and in time.Note that for the space-time coordinates, we have used t and r as independent variables.On the other hand, for the momentum variable, we only write out p explicitly, with E omitted, since for a given type of particle with mass m, E is already determined by p as E = m 2 +p 2 .Alternatively, one can argue that the Ψ p (x) single-particle wave function describing a particle with momentum p is a plane wave, and thus The essence of the description of the Bose-Einstein correlation lies in the statement that the analog expression of N 2 (p 1 , p 2 ) that takes into account particle production at two distinct space-time points x 1 and x 2 must be written up utilizing the two-particle wave-function Ψ (2) and that for bosons, Ψ p1,p2 (x 1 , x 2 ) is symmetric in the x 1 and x 2 variables.If the final state interactions of the particles are neglected, the wave function is a symmetrized plane wave, denoted by Ψ (2,0) p1,p2 : For the modulus square of this wave function, one obtains and from this, after some simplifications, one gets the following expression for the correlation function (where the 0 superscript denotes that final state interactions are neglected): with S(Q, p) being the Fourier transform of the source function: Typically the correlation function is written up as a function of the average and relative momentum variables K and Q, respectively, instead of p 1 , p 2 .The reason for this is that it turns out that the dependence on Q is much stronger than on K, so one almost always makes the so-called smoothness approximation [21] that in Eq. ( 9), p 1 ≈ p 2 ≈ K in the second argument of S. One then has and one can conceptually separate the Q and K dependence by assuming that the source function S(x, K) has a given functional form (e.g., a Gaussian) as a function of x, and the parameters of this functional form may depend on K.By measuring the correlation function as a function of Q at different K values, one can reconstruct the K-dependence of the source parameters.
We can express the Ψ (2,0) two-particle wave function in terms of the center-of-mass system variables and the relative variables as where In this form, we see that the relative motion and the center-of-mass motion disentangle.In the case of nonvanishing final state interactions, the center-of-mass motion still retains its plane-wave e −2iQX form.However, the form of the wave function pertaining to the relative motion changes.In this case, the two-particle wave function (now with the 0 superscript omitted, owing to the presence of final state interaction) is Here the ψ Q,K (x) part of the wave function describes the relative motion.It is symmetric to the x ↔ −x exchange (corresponding to the x 1 ↔ x 2 change), and its concrete form depends on the interaction.We see that only this function determines the object of interest, the modulus of the two-particle wave function: the plane wave corresponding to the free motion of the center of mass cancels.
With this notation, we may rewrite the expression of the two-particle momentum distribution and the (utilizing the smoothness approximation, p 1 ≈ p 2 ≈ K, just as above) as So with a new definition, we have where we introduced the D(x, K) relative coordinate distribution or spatial correlation function as Note that while S(x, p) is not necessarily an even function of x, this D(x, K) function is even in x.Taking this into account, and that from Eq. ( 13), we have |ψ , we can write the correlation function in the case of no final state interaction as We see that owing to the presence of the denominator in Eqs. ( 17) and ( 19), the overall normalization of D (the integral over space-time, at any given K) cancels from the measurable C 2 correlation function.To simplify our formulas, in the following, we thus assume that D is normalized to unity at any given K.
For actual calculations, it would be desirable to have the time variable (in an appropriate reference frame) canceled from the production described by S(x, p).The freeze-out duration ∆τ is in principle not zero (indeed, a prime goal of HBT studies is to measure ∆τ , and thus the order of phase transition).However, as multiple calculations show (see, e.g., Ref. [22] for an early exposition), the effect of ∆τ being non-zero can essentially be factored into the geometric (space-like) radii of the source function (with the time component of the relative momentum, q 0 being expressed with the spatial components along the way).So in the following, we write our formulas as if D(x, K) was non-zero only if t=0, highlighting this by omitting the time-dependence in the notation of our quantities, and remember that this is understood in the sense that the effect of finite time particle emission was treated by incorporating its effects into the geometrical radii of the source function.We thus have specifically, for the case when final state interactions are neglected, We also appropriately change our notation at this point: since we no longer need four-vectors, in the case of a three-vector (denoted by bold letter), the same letter in standard typeset indicates the magnitude of the vector, e.g., Q ≡ |Q|, r ≡ |r|.
In the case of the Coulomb interaction, the ψ(x) wave function (a solution to the Schrödinger equation satisfying appropriate boundary conditions) has a known form in the non-relativistic case.This implies that one has to use the center-of-mass frame of the particle pair (PCMS frame), where their motion can be approximated to be non-relativistic.In the following thus, we take Q to be understood in this PCMS frame and use the k = Q/2 notation, and (as customary for such treatment of the Coulomb interaction) assume that the simultaneity condition expounded above holds in PCMS.So, in this case, we have where we dropped the K arguments for simplicity, but remember that the source function depends on K through the K-dependence of its assumed parameters, and use the appropriately symmetrized ψ k (r) Coulomb interacting out state in the following form, well known in quantum mechanical scattering theory (see, e.g., Ref. [23] for details): The η quantity appearing here is the so-called Sommerfeld parameter: with α ≡ q 2 e 4πε0 1 ℏc ≈ 1 137 being the fine structure constant.The quantity N is the normalization factor of the wave function, of whose square modulus is the so-called Gamow factor : where Γ(z) is the Gamma function, and its simple properties (such as zΓ(z) = Γ(z+1) and Γ(z)Γ(1−z) = π sin(πz) ) lead to the expression written up for |N | 2 .Finally, M (a, b, z) is the confluent hypergeometric function, Besides the (modulus square of the) wave function, an assumption is needed about the functional form of the D(r) pair distribution in order to calculate the correlation function as in Eq. ( 22).As mentioned in Section 1, recent measurements utilized Lévy-stable distributions as an assumption for D(r).The appearance of Lévystable sources in high energy heavy-ion collisions was expounded in [24].Since then, several mechanisms were proposed as the possible cause for this particular source shape, such as jet fragmentation [25], critical phenomena [26], directional averaging and non-sphericality [27], event averaging [27], resonance decays and hadronic rescattering [28,29].Out of these, each may result in Lévy-shaped sources at different collision energies and systems: jet fragmentation in e + e − or pp collisions, critical phenomena at energies lower than the top RHIC or LHC energies, while hadronic rescattering may occur generally.In [29], it was shown that in realistic simulations of the hadronic medium, Lévy sources appear on an event-by-event basis at various stages of the evolution.
The symmetric Lévy distribution is a generalization of the Gaussian function as it has an additional parameter α, called the Lévy exponent.In the spherically symmetric case, it has the following form: By the known property of Fourier transforms, one sees that this function is normalized to unity: In the α = 2 case we recover a Gaussian distribution with radius R; while for α = 1, a Cauchy distribution is obtained: For α ̸ = 2 and α ̸ = 1, one cannot express L(α, R, r) with a combination of simple functions.However, its asymptotic expression is known: for every α ̸ = 2, it decreases like a power of r: as exemplified in the above α=1 case as well.
Lévy distributions generalize Gaussian distributions also in the sense that they retain the stability property of Gaussian: the convolution of two Lévy distributions (with the same α parameter) is again such a Lévy function, as can be directly seen by inserting the above definition (a Fourier integral) and performing the integrals.Specifically, From another point of view, with a similar calculation, we obtain that if the S(r, p) single-particle source function is a Lévy-type function, then so is the D(r, p) two-particle source function: We could even allow p-dependent parameters as well as an arbitrary r 0 (p) shift in S(r, p); this does not change the shape of D(r, p), the quantity that appears in the expression of the correlation function.
The stability property written up in Eq. ( 31) is the main reason why Lévy distributions are expected to appear in many different circumstances (as limiting distributions, just as Gaussians).In the scenarios mentioned above, they also naturally arise as particle production source functions in heavy-ion collisions, the main object of interest in this paper.Indeed, Lévy distributions have already been successfully applied to describe such sources [15], after earlier measurements indicated that a Gaussian description falls short of such data [13].However, as explained below, it turns out that the actual calculation of the C 2 correlation function for Lévy-like D(r, p) is cumbersome even numerically; there is thus a natural need for the development of such calculational methods.This becomes even more pressing if one abandons spherical symmetry or if one allows further generalizations of the source function beyond the Lévy assumption.

A new method for the treatment of the Coulomb effect
As written up in Eq. ( 22) above, for calculating the correlation function with the Coulomb interaction taken into account, one has to integrate the modulus of the Coulomb interacting wave function ψ k (r) weighted with the D(r) two-particle source function.This is deemed to be feasible only numerically, owing to the complicated form of ψ k (r), Eq. ( 23).Nevertheless, as seen above in Eq. ( 21), in case of vanishing Coulomb force, |ψ k (r)| is a plane wave, and our desired integral reduces to a Fourier transformation of D(r).
There is an important, if not the only interesting, class of possible source function assumptions that have the property of being defined and best calculable as a Fourier transform of some simple function.A Gaussian source function obviously falls into this category; however, a preeminent motivating example is when the D(r) two-particle source function is a Lévy distribution, whose Fourier transform, e −|Rq| α is an easily computable fast decreasing function.In the following, we thus assume that We assume that f (q) is an integrable function (over the whole q space); this is a natural assumption if we want to interpret this Fourier transform as a regular (Lebesgue-)integral.We also assume D(r) to be an integrable function (over the whole r space; this is naturally fulfilled whenever S(r) is integrable, which is necessary to interpret S(r) as the function whose integral, according to Eq. ( 4), gives the single particle invariant momentum distribution).These imply that both f (q) and D(r) are bounded and continuous, as well as that the (inverse) Fourier transform of D(q) is also a regular integral: and thus the normalization condition for D(r) means a simple condition on f : In particular, the interaction-free correlation function has f as the main component; from Eqs. ( 21) and (34), In the Coulomb interacting case, our goal is to calculate the correlation function using Eq. ( 22); put together with Eq. ( 33), we have where ψ k (r) is the Coulomb wave function written up in Eq. ( 23).The straightforward way to proceed used, e.g., in Refs.[18,30] in the case of Lévy distributions, is then to calculate D(r) from f (q) as a Fourier integral, then perform the integral over r, in which |ψ k (r)| 2 enters.Numerically, this is a daunting task in some cases.E.g., for the Lévy distribution, the result of the Fourier transform (the function D(r) itself) is only slowly decreasing (with a power-law-like behavior for large r), while |ψ k (r)| 2 is an oscillating function; asymptotically a plane wave (up to logarithmic corrections).There is also a conceptual awkwardness in this methodology: we take a Fourier transform of f (q) and then subject the resulting D(r) function to an "almost inverse Fourier transform" (i.e., the r-integral with the "almost plane wave" |ψ k (r)| 2 as a kernel) to finally arrive at C 2 (Q).In the interaction-free case, |ψ k (r)| 2 is really a plane wave, and the result, Eq. ( 36) indeed shows that these two Fourier transformations cancel each other.One would thus very much prefer a calculational scheme where this "back and forth" transformation is not needed in its full numerical complexity.(For example, calculating C 2 (Q) at high Q values, where owing to the decrease of f (q) it is increasingly closer to unity, still requires significant computational power because at higher k, the oscillations of ψ k (r) become faster and faster.) The natural idea to resolve these problems is that we would like to perform the Fourier transform of the interacting |ψ k (r)| 2 function and "act" with this resulting integral kernel (a function of k and q) on the f (q) function.In other words, we would like to interchange the order of integrals in Eq. (37) so that we could write Regrettably, however, this is not possible in such a simple form: since ψ k (r) is asymptotically a plane wave, its modulus is definitely not an integrable function whose Fourier transform could simply be calculated by an integral.In the following, we work around this problem by carefully treating integrability and the interchange of integrals and limits, a rarely found exercise in physics-motivated standard calculations.The resulting formulas, however, are worth such careful investigation.In doing so, we utilize some fundamental theorems about (Lebesgue) integrability. 1 We outline the main steps of our calculations below but leave some mathematical and calculational details to the Appendix.
It turns out that we can proceed by inserting an exponential "regularization", e −λr into the integrand of Eq. ( 37), with a positive real λ > 0 parameter, which at the end goes to 0. With finite λ we can interchange the order of integrals and finally arrive at = lim (Remember that Q is the variable of the resulting observable C 2 correlation function, Q = 2k, while q is a mere auxiliary integration variable in this calculation.)In Eq. ( 39), every step is justified mathematically, in particular, the exchange of the limit and the integral in Step 1 (by virtue of Lebesgue's theorem), and the exchange of the order of integrals in Step 2 (by virtue of Fubini's theorem); see Appendix A.1 for details.Nevertheless, the point is that in the resulting final formula, the λ → 0 limit cannot be exchanged with the q-integral.The way to proceed is to calculate the r-integral for finite λ values, and then "simplify" the λ → 0 limit to arrive at a final formula where there are no explicit limits that would have to be evaluated numerically (which would be very challenging if not impossible).
In this paper, we proceed with pair source functions that are spherically symmetric; this also implies spherical symmetry of f (q), which we highlight everywhere with an s in the subscript.In the spherically symmetric case, we can perform the solid angle integral in the Fourier integral, arriving at the following expression of D(r): Following the exact same steps as in Eq. (39), with this expression for the spherically symmetric D(r) we have Note that while ψ k (r) depends on the direction of k = Q/2, the result of this integral indeed only depends on the magnitude Q, after performing the r-integral.
As a next step, we substitute ψ k (r) from Eq. ( 23).The modulus square results in four terms (from each pairings of the different M (a, b, z) functions).With an appropriate r → −r change of variable in two of these, we are left with only two different terms: The integrals of Eqs. ( 43) and (44) are calculated in Appendix A.2; the utilized method is similar to that of Nordsieck [31].In our case, the result is where we make use of the (ordinary) hypergeometric function: This power series is valid only for |z| < 1, however, the domain of 2 F 1 (a, b, c, z) can be extended to any z complex number arguments by analytic continuation; see e.g.Ref.
[32].The next step is then to try to evaluate C 2 (Q) as in Eq. ( 42) from these expressions of D 1λs and D 2λs .For this, we need to take the λ → 0 limit in a careful way.It turns out that as λ approaches 0, the functional forms of D 1λs and D 2λs become ill-behaved; the desired limit (of the result of the q-integral) cannot simply be expressed as some limiting q-function multiplied by f s (q) and then integrated.Instead, in the final expression, one has to take a specific linear functional of f s (q). 2 With leaving some intermediate steps to Appendix A.3, in the following Section, we begin by writing up the main result of our calculation for C 2 (Q).45) and ( 46), a careful utilization of Lebesgue's theorem (detailed in Appendix A.3) allows one to perform the remaining λ → 0 limit in Eq. (42).We thus arrive at the main result of our new method for calculating the Coulomb interacting C 2 (Q), which can be summarized as follows:

Results and discussion
where the A 1s and A 2s terms are the following functionals of f s (q): If η → 0 (either by letting k → ∞, or by formally turning off the Coulomb force), A 1s and A 2s take on finite values (because both the η denominator and the imaginary values of the denoted quantities go to 0).So for η → 0, in Eq. ( 48) the contributions of A 1s and A 2s vanish (where they are again multiplied by η); we wrote the η factors in the way we did to highlight this feature.Note also that in Eq. ( 48), the effect of Coulomb interaction enters in a straightforward, "traceable" way.Omitting |N | 2 as well as the A 1s , A 2s terms, and recalling that Q = 2k, we simply have Eq.(36), the result for the interaction-free case.Including the |N | 2 factor (but not A 1s and A 2s ) is the so-called Gamow correction, the 2 The situation is similar in the well-known much simpler case when one approximates the Dirac delta with ever narrower and higher peaks with area of unity.Consider the following identity: Based on this, one can say that the narrower and higher Lorentz curves "approximate the δ(x−x 0 ) delta function"; the meaning of this statement is essentially the same as the displayed identity.From a practical point of view, the benefit of this formula is that instead of numerically performing some λ → 0 limit of the result of the left-hand side integral, the right-hand side provides a much simpler calculational statement (by requiring just to evaluate the f function at a given x 0 point); however, this simpler statement is no longer a true integral transformation acting on f .simplest approximate treatment of the Coulomb interaction: it treats the source as point-like for the Coulomb integration (but not for the calculation of the interaction-free correlation function).The A 1s and A 2s terms can thus be thought of as the effect of the source being not point-like; the parts of the correlation function that is neglected by the Gamow prescription.
In our normalization, f s (0) = 1; however, we retained f s (0) in Eq. ( 49) to highlight that the fraction containing the f s function is a well-defined function even around q = 0 (provided that f s is continuously differentiable).Similarly, the fraction containing f s in the expression of A 2s , Eq. ( 50), is also a well-defined function for continuously differentiable f s , even around q = 2k.In fact, we can loosen the requirement of f s being continuously differentiable.The more general assumption under which the above formulas are derived is explained (along with other details) in Appendix A.3 around Eqs. ( 82) and ( 83).Here we just remark that this more general assumption is satisfied already if f (q) is everywhere continuously differentiable, and also in case of Lévy distributions (for which f (q) is not differentiable at the origin if α ≤ 1).
Power functions of complex variables are understood in the strictly single-valued function sense: for z, w ∈ C, z w := exp(w Ln z), where exp is an entire function, its inverse however, Ln z ≡ log(|z|) + iargz, has a branch cut along the R − 0 negative real line (owing to the jump of arg z, the phase of z, from π to −π, when crossing this line in a counterclockwise direction).This branch cut along R − 0 is inherited by the power function z w as a function of z, whenever w / ∈ Z.So when taking a power of a quantity that is on this branch cut, we have to specify which side of the cut we are on.This is the reason for the +i0 term in Eq. ( 50).Also when one defines the 2 F 1 (a, b, c, z) hypergeometric function for |z| ≥ 1, that is, outside the convergence radius of the power series in Eq. ( 47), one encounters a branching point at z = 1, and the usual convention places the branch cut on the [1, ∞] ⊂ R + real line.The −i0 term in Eq. ( 49) specifies which side of this cut we are on.

Utilizing the new formula for Lévy-stable sources
Before delving into the details of the applied computational methods, in Figure 1, we show the main results of our calculations, i.e., the C 2 (Q) correlation functions calculated for Lévy-stable source distributions L(α, R, r), separately for pion-pion and kaon-kaon pairs.Different colors denote different Lévy scale R values, while the dashed and solid lines correspond to the Lévy exponent α values of 0.6 and 2, respectively.Correlation functions with α between these two extremal values are illustrated with a filled area.In the interactionfree case, for any given R, correlation functions take the same value for any α at one specific Q > 0 value: In the case when Coulomb interaction is included, the curves with different α values still intersect each other at approximately the same point (for a given R), in the vicinity of Q ≈ ℏc/R.However, in this case, the C 2 (Q ≈ ℏc/R) values have an R dependence.Moreover, the intersection of all the curves is only approximate (albeit a very good one); higher zooming around the apparent "nodes" would tell the curves apart.

Testing the final numerical integral
As described above in Eq. ( 48), for the correlation function, two numerical integrals have to be performed.For a reasonable f s (q) function, such as the one for Lévy distributions, the integrands there are particularly wellbehaving (non-oscillating, smooth) functions.The integrals over q ∈ R + 0 are best transformed to integrals over the [0, 1] interval, with a new integration variable x introduced separately on the q ∈ [0, 2k] and the q ∈ [2k, ∞] intervals as q = 2kx, and q = 2k x , respectively.The integrands turn out to be well-behaving "tame" functions of this new x variable; they are written up as The integrals can then be evaluated using the rectangle method or the trapezoidal rule, but since the individual function evaluations have a high cost in terms of CPU time and in practical applications such as optimization (fit) procedures, the final result has to be calculated several thousands of times at a given Q, it is beneficial to search for methods requiring fewer evaluations.We find that the Gauss-Kronrod quadrature formula [33] provides an acceptable solution for this.This method works well in our case because the function to be integrated is changing fast near 0 and 1 but very slowly in the middle.Hence a varying bin width has to be applied, and such naturally arises from the Gauss-Kronrod iteration.
We utilized the boost C++ library to perform this integral and find that 15 nodes provide a fast converging result.The number of function evaluations and the integral results are shown in Fig. 2.Here we chose a large Q value, as the accuracy of the original numerical calculations, when one directly calculates the integral with |ψ k (r)| 2 , decreases towards large Q, and hence accuracy can be well tested.The integral result changes on the order of 10 −6 , and even on that scale, it converges when the number of maximum iterations is set to 3 − 4. It is also important to note that the typical order of magnitude of the bin-by-bin statistical uncertainties of a correlation function measured in high-energy experiments is significantly higher than the change in the integral values shown in Fig. 2. If the tolerance parameter of the quadrature is larger than 10 −3 , then the integral result differs slightly more from the plotted results.However, beyond that, if tolerance is at least 10 −3 , increasing it further does not modify the integral result.On the other hand, the number of function evaluations increases fast, especially if the tolerance is very small and the number of iterations is large.Hence an optimal solution (analyzing results similar to Fig. 2 for many R, α, and Q values) is provided by a tolerance value of 10 −3 and a limit of 3 for the maximal number of iterations.With this, one integral requires up to a few hundred function evaluations, suitable for simultaneously fitting many experimental C 2 (Q) datasets.

Comparison with the original numerical integration method
After carefully testing the reliability of the Gauss-Kronrod numerical integral, we proceeded with comparing the new calculation to a previously utilized integral method based on Eq. (37), and used, e.g., in Ref. [18].The previous method requires significantly more computational effort to achieve a reasonable precision at small α or large R values.In that case, the values of the correlation function were pre-calculated for various α and R values and saved in a large lookup table; from this table, an interpolation can be used to get the correlation function for any parameter combinations.In the following, let us call our new approach "wave function Fourier method" and denote the result with WFF C 2 (Q).The previous method with the pre-calculated lookup table is denoted with table C 2 (Q).
To illustrate the difference between the two methods, we plot for various α and R values, as shown in Figure 3.We find that ∆C 2 (Q) has a small dependence on Q, generally its magnitude is larger at larger Q values.It is also evident that the difference between the two methods is largest at the smallest α values, as expected.To better illustrate the dependence of ∆C 2 (Q) on R and α, in Figure 4 we plot the Q-average ⟨∆C 2 (Q)⟩ values.It can be clearly seen that the ∆C 2 (Q) difference between the two methods is smallest at large α and small R values.
The conclusion of this comparison is that, on one hand, we can be reasonably assured that our new methodology works well; as cautious users, we really strive for numerical verification, even if one is convinced by the  flawlessness of the mathematical derivation.On the other hand, our new methodology offers, at a significantly lower numerical cost, a way of calculating the C 2 (Q) correlation function more robust than the previous method, especially at large Q values.For this reason, we find that the new method is ready to be used in experimental analyses of measured C 2 (Q) correlation functions.

Summary
In summary, we presented a novel method for calculating Bose-Einstein correlation functions with the Coulomb effect incorporated.Our method can be applied to any particle emitting source function assumptions that are expressed as a Fourier transform.In this way, our method starts directly from the interaction-free correlation function, which is the Fourier transform of the source function.The integrals necessary for our calculational scheme previously have eluded exact calculations, often resulting in experimental analyses using a fixed source size when correcting for the Coulomb effect.With precise data, however, this is no longer a viable method; our approach allows for simple and more exact handling of the Coulomb final state interactions.We demonstrated that for Lévy-stable sources (including the Gaussian limiting case), the results are close to those obtained with numerical integrals previously.On the other hand, at extreme values of the Lévy source parameters (small Lévy exponent α, and large Lévy scale R), the new method is more precise and reliable.Furthermore, the new method allows for a direct and fast fitting of correlation functions obtained in experiments.We also published a software package [19] that is easily applicable and ready to use in experimental analyses.
The natural next step is to extend the methodology to the case when the assumption of spherical symmetry is rid of.In that case, one would have a final formula that is a linear functional of the Fourier transform of the pair source function, D(r), denoted by f (q) in this paper.It turns out that the necessary integrals (similar to but slightly more complicated than Eqs.( 43) and (44) in this paper) are also readily calculable, and then the necessary limit when the regularizing parameter is removed (denoted as λ → 0 in this paper) is also manageable.The resulting formulas are, however (owing to the inherently more complex nature of the nonspherically symmetric case) somewhat more complicated, thus their analysis and implementation go beyond the scope of this paper.

Acknowledgements
This research was supported by the NKFIH grants K-136138 and TKP2021-NKTA-64.

A Assorted mathematical and calculational details
In the Appendices, we present the calculations that were left out of the body of the paper.We rely heavily on complex analysis as well as on some intricacies of (Lebesgue) integrability; the latter of which we give a brief summary here.(A more detailed discussion can be found in any of the standard textbooks on mathematical analysis, such as Ref. [34].)

A.1 Discussion of Eq. (39)
The focal point of our method of circumventing the need for a back-and-forth integral transformation of the Fourier transform of the source function (denoted by f (q) in Eq. ( 33) and from that on) is Eq.(39), the mathematical justification of whose steps is thus essential.The main ingredients are conditions of integrability, Lebesgue's dominated convergence theorem (or simply the Lebesgue theorem), and Fubini's theorem.A summary of these theorems from a practical point of view is: • Concerning Lebesgue integrability, the key is the modulus of a function.A (real-valued, or complex-valued, or any normed space valued) function F is integrable if and only if |F | is integrable.Also, if there is a G function for which |G(x)| ≥ |F (x)| for almost all x, and G is integrable, then F is also integrable.(Here, we denoted the integration variable by x; it can be of any type of real integration variable.)Conversely, if |G(x)| ≥ |F (x)|, and F is not integrable, then neither is G.
• Let F λ (x) be integrable functions, with λ as a parameter (either a continuous one or a discrete index), and for almost all x let the pointwise limiting function F (x) := lim λ F λ (x) exist (for any reasonable type of limit; for example, λ → 0, or λ ≡ n → ∞).Lebesgue's theorem states that if there is a (λ-independent) G function for which for almost any x, |G(x)| ≥ |F λ (x)|, for any λ, then the integrals of the F λ functions (the F λ values) converge, the limiting function F is integrable, and the limit of the integrals is equal to the integral of the limiting function, lim λ F λ = F , i.e., the limit and the integral are interchangeable.
The key is the existence of the "dominant" G function: in the case when there is no such function, none of the statements of the theorem are necessarily true.A simple (counter)example is the approximation of the Dirac delta with functions whose pointwise limit is everywhere zero; in this well-known case, one maybe does not even recognize this peculiarity.However, for less explored cases (like ours in this paper), one has to be careful when interchanging limits and integrals, and then this theorem comes in handy in many cases.
• Fubini's theorem concerns multiple integrals, iterated integrals, and the justification of interchange of integrals.The main point is if an F (x, y) function is such that its modulus is integrable in one order as a repeated integral; i.e., if the dx dy |F (x, y)| integral exists, then F itself is integrable in both orders, and these integrals coincide: dx dy F (x, y) = dy dx F (x, y), i.e., the integrals are interchangeable.
The transformation in Eq. ( 39) is again written up, with the intermediate steps slightly more detailed, as = lim In Step 1 we inserted the e −λr regularizing factor; its limit being lim λ→0 e −λr = 1.Because D(r) is integrable, so is it multiplied by |ψ k (r)| 2 (which is a bounded function), and thus |ψ k (r)| 2 • |D(r)| is a good dominant function (that is, integrable and greater or equal than the integrand, independently of λ).So by virtue of Lebesgue's theorem, we could interchange the λ → 0 limit and the r-integral in Step 2. In Step 3, we inserted the Fourier integral expression of D(r), careful (yet) about the order of integrals.However, since |e iqr | = 1, the double integrand here has the modulus that is a product of |f (q)|, an integrable q-function (as assumed), and an integrable r-function, e −λr |ψ k (r)| 2 .So the repeated integral of the modulus exists, from which it follows (by Fubini's theorem) that we can interchange the original integrals as well; this is Step 4 here.
In the resulting right-hand side, we cannot perform the λ → 0 limit in the integrands; that would mean the exchange of the original r-and q-integrals, which is not possible, as stated around Eq. (38).Instead, we have to calculate the integrals here, and then perform the λ → 0 limit.

A.2 Calculation of D 1λs and D 2λs
In this Appendix, we derive the Eqs.( 45) and (46), the results of the integrals D 1λs and D 2λs , defined in Eqs. ( 43) and (44).We utilize Nordsieck's method [31], who applied a similar technique to simplify a very similar integral that occurs in the theory of bremsstrahlung and pair creation.We write up the definitions (43)-(44) again, with a slight change of notation, to express the sin(qr) function with exponentials: where (±) 2λs (q) = 1 q d 3 r e −λr r e ±iqr M 1+iη, 1, −i(kr−kr) M 1−iη, 1, i(kr+kr) . (56) We use the following complex contour integral representation of the confluent hypergeometric function: where the path is any closed curve that encircles the real line segment [0, 1] once counterclockwise on the complex t plane; this segment (a branch cut) is the only set on the t plane where the integrand is not analytic.
Inserting twice this expression (with the integration variables denoted by t and u, whose paths are such as just specified) into Eqs.( 55) and (56) for the two confluent hypergeometric functions in each, we get We would like to exchange the order of the integral over r and the contour integrals because for the r-integral, we could then use the following auxiliary formula, valid for any β ∈ C complex number and B ≡ (B x , B y , B z ) three-vector with any complex components: Here B 2 = B 2 x +B 2 y +B 2 z even for complex B x , B y , and B z (in particular, without complex conjugation), and the real part of the B vector is The condition written up in Eq. ( 60) is necessary and sufficient for the integral to exist.This is because a function is (Lebesgue) integrable if and only if its modulus is integrable, and 1 r e −βr+Br = 1 r e −Re β•r+(ReB)r , which is integrable if and only if the multiplier of r in the exponent, Re β, is strictly bigger than the length of the vector there, |Re B|.For real β, B x , B y , B z , the stated result is elementary, and because both the integral itself and the stated result are analytic functions of β, B x , B y , B z (provided the integral exists), the result is valid for all the allowed complex values, as stated.
If we can exchange the r-integral and the contour integrals in Eqs. ( 58) and ( 59), the we would get According to Fubini's theorem, the interchange is justified if the modulus of the integrand is integrable in either order (provided that one has written up the contour integrals in a parametrized way, with integrals taken over line segment of the parameter).In our case, the moduli of the complex powers of t and u do not cause concern: since u and t run at a positive distance from the branch cut, these factors are bounded and thus do not spoil integrability.The factor of the modulus of the integrand that indeed has a role is that of the r-dependent part.
In the case of D 1λs , Eq. ( 61), the modulus of the part of the integrand that is of interest is e −λr±iqr−i(t−u)(kr+kr) = e −[λ−kIm(t−u)]r+Im(t−u)kr , which is integrable over r if and only if for the following holds (by virtue of the condition stated in Eq. (60): We arrived at the simpler condition by sorting through all the different possibilities for Im u and Im t.Similarly, in the case of D 2λs , Eq. (62), we have e −λr±iqr−i(t−u)kr+i(t+u)kr = e −[λ−kIm(t−u)]r−Im(t+u)kr , which is integrable over r if and only if If these conditions hold for all possible u and t, then the results of these r-integrals are continuous bounded functions of u and t, so their contour integrals exist.So the condition for the integrals over r and u, t to be interchangeable is that the inequalities stated in Eqs. ( 63) and (64) hold for any u and t on their integration paths.Thus if we require these additional constraints on the u-and t-paths, we can interchange the integrals, and utilizing Eq. (60), we get the following from Eqs. ( 61) and (62): where the paths of the u and t variables run as specified: they both encircle the [0, 1] branch cut and obey at all points the conditions (63) in case of D 1λs , and (64) in case of D 2λs , respectively.(For any finite λ and k, it is indeed possible to specify the integration paths in this way.)The case of D 2λs is simpler because the integrand factorizes in u and a t; we have written Eq.(66) in a way that highlights this.The next step is to perform the u-integral at a fixed t.As a function of u, the integrand in both cases is analytic everywhere except for the branch cut along u ∈ [0, 1] which is encircled by the path, as well as for a simple pole, denoted in case of D 1λs and D 2λs by u 1s (t) and u 2s (t), respectively: We see that Im(t−u 1s ) = λ 2k and Im(u 2s ) = − λ 2k : this means that (for any t that is in its allowed domain) u 1s and u 2s do not satisfy the conditions in Eqs. ( 63) and (64).From this, one concludes that u 1s and u 2s lies outside of the integration contour on the u plane: were it otherwise, a narrower integration path chosen for u could cross u 1s and u 2s , but this is impossible because then u 1s and u 2s could not violate the conditions.
One can then expand the u-integration contour to infinity: the integrand decreases rapidly enough (as ∼ 1 u 2 in both cases), so the non-vanishing contribution comes from the simple poles u 1s and u 2s , so from the residue theorem we have (with an extra minus sign from the negative sense of the paths encircling the poles) where in case of D 1λs , u 1s (t) is as in Eq. (67), while in case of D 2λs we inserted u 2s from Eq. (67) right away. 3n Eq. ( 69), we can use this same method for the t-integral.The integrand has a simple pole in t at t 2s = iλ 2k ± q 2k , and the same reasoning as for u 2s above yields that because Im t 2s = λ 2k , t 2s does not satisfy the condition in Eq. ( 64), so it lies outside of the integration contour.Thus an expansion of the contour to infinity yields (again because of the ∼ 1 t 2 decrease) just the contribution from the pole at t 2s .From this, we get This leads to the result for D 2λs , stated in Eq. (46) in the main text, after some simplification.
In the case of D 1λs , the remaining integrand in Eq. (68) does not have a simple pole as a function of t (besides the branch cut on [0, 1]), but another branch cut on a straight line segment between t = ± q 2k +i λ 2k and t = 1± q 2k +i λ 2k (i.e., where u 1s (t) = 0 and u 1s (t) = 1).This branch cut, as its imaginary part, λ 2k is in violation of condition (63), lies outside the t integration contour.In this case, we can transform the integral in Eq. ( 68) Finally, we thus arrive at the following expression, where we can utilize the (73) integral representation: and from this we get the result given for D 1λs in the body of the paper, Eq. ( 45), by substituting N , β, B, and x, and finally combining together D A.3 Derivation of the final formula, Eq. ( 48) In this Appendix, we derive Eq. ( 48) along the expression of the quantities denoted by A 1s and A 2s in the body of the paper, Eqs. ( 49) and (50).The starting point is Eq. ( 42), the expression of C 2 (Q) with D 1λs and D 2λs as a λ → 0 limit.We write this up again but first with the integral denoted as a three-dimensional q-integral, and perform the first crucial step as we "separate" the f s (0) and f s (2k) values.This idea comes from observing that the values of f s (q) at q = 0 and q = 2k play a special role in the interaction-free case as C (0) 2 (Q) = f s (0)+f s (2k), moreover, f s (0) and f s (2k) enter by virtue of the "plain" and the "cross" terms from the modulus square of the wave function, just as D 1λs and D 2λs in the interacting case.We thus write Now for λ>0, the integrals of D 1λs and D 2λs , denoted temporarily by I 1 and I 2 , do exist: I 1 := d 3 q D 1λs (q) = 4π ∞ 0 dq q 2 D 1λs (q), I 2 := d 3 q D 2λs (q) = 4π which can be directly seen from their expressions, Eqs.(45) and (46): q 2 D 1λs (q) and q 2 D 2λs (q) are continuous and bounded on every compact [0, q max ] interval, and they decrease as ∼ 1 q 2 , owing to the taking of the imaginary part in them; recall that 2 F 1 (a, b, c, x) = 1 + ab c x + . . .as x → 0. So the I 1 and I 2 integrals are given by the value of the Fourier transforms of D 1λs and D 2λs at r = 0 (where r is the variable of their Fourier transforms). 5But we actually calculated D 1λs and D 2λs as Fourier transforms, so we just need to evaluate the integrands of the defining (43)-(44) integrals at r = 0, from which, knowing that M (a, b, x=0) = 1, we have I 1 = I 2 = 8π 3 , so q 2 D 1λs(q) 2π 2 f s (q)−f s (0) + lim λ→0 ∞ 0 dq q 2 D 2λs (q) 2π 2 f s (q)−f s (2k) .(79 Our statement is now that these remaining limits result in A 1s and A 2s as written up in Eqs. ( 49) and (50).Substituting D 1λs and D 2λs from Eqs. ( 45) and (46), and using the convenient F + (x) ≡ 2 F 1 (iη, 1+iη, 1, x) notation, we thus have to prove that Im (q+2k) iη (q−2k+i0) iη .(81) Note that the pointwise λ → 0 limit of the integrands on the left-hand sides are just the ones on the right-hand sides.We want to apply Lebesgue's theorem (see Appendix A.1 above) so that we can interchange the limit and the integral.To this end, we need a "dominant" function, i.e., an integrable function whose modulus is greater or equal to that of the integrand, independently of λ, for any λ > 0. We do not explicitly write up this function, however, we give the various estimations that are needed to construct this function.
1.The f s function is bounded everywhere, f s (q) ≤ K. Also, the main point of introducing the subtraction of f s (0) and f s (2k) values as in Eqs. ( 80) and ( 81) is that for the fractions appearing from these can be dominated by a function that is integrable around q = 0 and q = 2k if f s is continuously differentiable everywhere: in this case for a finite q max value, |f (q) − f (0)| q , and |f (q) − f (2k)| q − 2k are bounded on [0, q max ], and thus integrable here.Had we not have subtracted f s (0) and f s (2k), the fractions fs(q) q and fs(q) q−2k would have not be integrable around q = 0 and q = 2k.
In the case of our preeminent practical example of source functions, the Lévy distribution, f s (q) is not necessarily differentiable around q = 0.However, for our estimation to find the dominant function here, it is enough if the fractions written up in Eq. ( 82) are not bounded but can be dominated by an integrable function.This is readily satisfied also in the case of Lévy distributions, where the following estimation is true instead: For Lévy distributions, for any q and q ′ , |f (q) − f (q ′ )| ≤ K ′ |q − q ′ | α ; this is enough.(83) 2. For the various factors appearing, the following estimations can be made.For real X and Y , the pure imaginary power X iY is bounded by e −πY from below and by e πY from above.The F + (x) hypergeometric function is bounded on the whole complex x plane. 6  3. Owing to the continuous differentiability of F + (x) around x = 0, there exists a C ′ constant for which |F + (x) − 1| ≤ C ′ |x| for small enough x.Also, for pure imaginary powers, the following estimation can be made: if Y ∈ R, 0 ≤ R < 1 and |z| < R, then (1+z From these, with some effort, one can construct the dominant function; this is best done by separately treating the q ∈ [0, Q max ] and the q ∈ [Q max , ∞] intervals, where Q max is an arbitrary value that is at a safe distance above 2k; say: Q max = 4k.In this way, the exchange of the integral and the limit in Eqs.(80) and (81) becomes justified, which leads to the desired results for A 1s and A 2s .

4. 1
The new formula for C 2 (Q)Starting from the expression of C 2 (Q) from Eq. (42), and substituting D 1λs and D 2λs from Eqs. (

Figure 1 :
Figure 1: Example correlation functions for pions (left) and kaons (right), plotted for four different R and two α values.At a given R value, the shape of the correlation function with increasing α values from α = 0.6 to α = 2 goes smoothly through the shaded region.

Figure 2 :
Figure 2: The result of the integral for R = 6 fm and α = 1.4,at Q = 0.2 GeV/c (left) and the number of required function calls (right).Note that the integral result is represented as(1 − C 2 (Q)) • 10 4, so the differences are on the order of 10 −6 .Both plots are shown as a function of the tolerance and the number of maximal iterations allowed in the Gauss-Kronrod integral.

Figure 3 :Figure 4 :
Figure3: Difference between the correlation function calculated with a numerical integral method described in Ref.[18] ( table C 2 (Q)) and the correlation function calculated with the wave-function Fourier method described in the current paper ( WFF C 2 (Q)).∆C 2 (Q) is plotted for 6 different α values and two R values, for pions (left) and kaons (right) separately.At a given α value, ∆C 2 (Q) goes smoothly through the shaded region when increasing R values from R = 3 fm to R = 12 fm.