Coulomb Final State Interaction in Heavy Ion Collisions for Lévy Sources

Abstract

Investigation of momentum space correlations of particles produced in high energy reactions requires taking final state interactions into account, a crucial point of any such analysis. Coulomb interaction between charged particles is the most important such effect. In small systems like those created in $e^{+}{e}^{-}$ or $p+p$ collisions, the so-called Gamow factor (valid for a point-like particle source) gives an acceptable description of the Coulomb interaction. However, in larger systems such as central or mid-central heavy ion collisions, more involved approaches are needed. In this paper we investigate the Coulomb final state interaction for Lévy-type source functions that were recently shown to be of much interest for a refined description of the space-time picture of particle production in heavy-ion collisions.

1. Introduction

Coulomb repulsion is the most important final state interaction that has to be considered in Bose-Einstein correlation measurements in high-energy physics. In $e^{+}{e}^{-}$ or $p+p$ collisions, where the particle emitting source is much smaller than the wavelength corresponding to the relative momentum of the particle pair, the well-known Gamow factor (essentially the value of the Coulomb interacting pair wave function at the origin) can be used to “correct” for the Coulomb effect. However, for an extended source, the Gamow factor overestimates the correction. A more advanced approach is to take the source-averaged Coulomb wave function (instead of its value at the origin, which may be valid then for a point-like source), see e.g., Refs. [1,2]. In these papers a method (aptly referred to as the Bowler–Sinyukov method) is also described that is widely used to take the effect of long-lived resonances into account.

Traditionally one assumes simple source function shapes (such as exponential, Gaussian ones) for calculating the source averaged Coulomb wave function (as e.g., in the papers referred above); we may mention that more general sources are also considered e.g., in Ref. [3]. Recently, an even more general type of source functions, namely Lévy sources [4,5] has gotten much interest. Lévy-type source functions simplify to Cauchy as well as to Gaussian ones in special cases, and allow for a more refined treatment of the space–time picture of the particle emission. Moreover, a certain parameter of a Lévy distribution (the so-called Lévy exponent) also may carry information about the order of the phase transition between deconfined and hadronic matter [5,6].

Our objective in this paper is to tackle the effect of Coulomb interaction for the case of Lévy-type sources. (Because of the slow, power-law-like decay of Lévy-type sources at large distances, many previously developed methods are unsuitable for them.) We strive for analytical approximate methods that are well suited for use in the actual treatment of experimental Bose–Einstein correlation functions.

2. Coulomb Effect in Bose–Einstein Correlations: Basic Concepts

In this section we briefly review some notions and well-known formulas pertaining to the work presented hereafter.

In a statistical physical (specifically, hydrodynamical) description of particle production in high-energy collisions, a basic ingredient is the (one-particle) source function (Wigner function), denoted here by $S(x,p)$. Its physical meaning is essentially that the probability of the production of a particle in the infinitesimal phase-space neighborhood of momentum $\mathbf{p}$ and point $\mathbf{r}$ is proportional to $S(\mathbf{r},\mathbf{p})d^3\mathbf{r}{d}^3\mathbf{p}$. Thus it is natural that the one-particle momentum distribution function $N_1\left(\mathbf{p}\right)$ can be expressed as

$$N_1\left(\mathbf{p}\right)=\int d^3\mathbf{r}S(\mathbf{r},\mathbf{p}),\mathrm{with}\mathrm{the}\mathrm{normalization}\int d^3\mathbf{p}{N}_1\left(\mathbf{p}\right)=1.$$

(1)

For a slight convenience we chose the normalization condition of $S(\mathbf{r},\mathbf{p})$, so that $N_1\left(\mathbf{p}\right)$ is now considered to be the probability distribution of the momentum of the produced particles1.

According to a simple quantum mechanical treatment of Bose–Einstein correlation effects, the two-particle momentum distribution function $N_2(\mathbf{q},\mathbf{K})$ can be expressed [7] with the source distribution function $S(\mathbf{r},\mathbf{p})$ as an integral over the two-particle final state wave function,

$$N_2({\mathbf{p}}_1,{\mathbf{p}}_2)=\int d^3{\mathbf{r}}_1{d}^3{\mathbf{r}}_{2}S({\mathbf{r}}_1,{\mathbf{p}}_1)S({\mathbf{r}}_2,{\mathbf{p}}_2){\left|{\psi }_{{\mathbf{p}}_1,{\mathbf{p}}_2}^{\left(2\right)}({\mathbf{r}}_1,{\mathbf{r}}_2)\right|}^2.$$

(2)

The two-particle wave function ${\psi }^{\left(2\right)}$ must be symmetric in the space variables (for bosons); this is the main reason for the appearance of quantum statistical (Bose–Einstein) correlations.

With some trivial simplifications, we thus get the correlation function $C_2(\mathbf{q},\mathbf{K})$ as

$$C_2(\mathbf{q},\mathbf{K})=\int d^3\mathbf{r}D(\mathbf{r},{\mathbf{p}}_1,{\mathbf{p}}_2){\left|{\psi }_{\mathbf{q}}^{\left(2\right)}\right|}^2,\mathrm{where}D(\mathbf{r},{\mathbf{p}}_1,{\mathbf{p}}_2)=\int d^3\mathbf{R}S\left(\mathbf{R}+\frac{\mathbf{r}}{2},{\mathbf{p}}_1\right)S\left(\mathbf{R}-\frac{\mathbf{r}}{2},{\mathbf{p}}_2\right).$$

(3)

The momenta of the two particles were denoted by ${\mathbf{p}}_1$ and ${\mathbf{p}}_2$, and we used the combinations of these, the $\mathbf{q}$ relative momentum and the $\mathbf{K}$ average momentum as

$$\mathbf{q}={\mathbf{p}}_1-{\mathbf{p}}_2,\mathbf{K}=\frac{1}{2}\left({\mathbf{p}}_1+{\mathbf{p}}_2\right).$$

(4)

The notation $D(\mathbf{r},{\mathbf{p}}_1,{\mathbf{p}}_2)$ was introduced for the so-called two-particle source function, obtained as indicated, by integrating over the average spatial position $\mathbf{R}$ of the particle pair (with $\mathbf{r}\equiv {\mathbf{r}}_1-{\mathbf{r}}_2$ thus standing for the relative coordinate).

The (symmetrized) two-particle wave function may depend on all momentum and coordinate components; however, its modulus does not depend on the average momentum $\mathbf{K}$ or the average coordinate $\mathbf{R}$ (owing to translational invariance).

Assuming ${\mathbf{p}}_1\approx {\mathbf{p}}_2$ (in the sense that the pair wave function changes much more rapidly in terms of their difference, $\mathbf{q}$, as does the product of the source functions, we get

$$C_2(\mathbf{q},\mathbf{K})\approx \frac{\int d^3\mathbf{r}D(\mathbf{r},\mathbf{K}){\left|{\psi }_{\mathbf{q}}^{\left(2\right)}\left(\mathbf{r}\right)\right|}^2}{\int d^3\mathbf{r}D(\mathbf{r},\mathbf{K})},$$

(5)

where we introduced $D(\mathbf{r},\mathbf{K})\equiv D(\mathbf{r},\mathbf{K},\mathbf{K})\approx D(\mathbf{r},{\mathbf{p}}_1,{\mathbf{p}}_2)$. In the special case of no final state interactions (i.e., when the ${\psi }^{\left(2\right)}$ wave function is a symmetrized plane wave), we get the well known relation

$$|{\psi }_{\mathrm{free}}^{\left(2\right)}{|}^2=1+cos\left(\mathbf{q}\mathbf{r}\right)\Rightarrow C_2^{\left(0\right)}(\mathbf{q},\mathbf{K})\approx 1+\frac{|\tilde{S}{(\mathbf{q},\mathbf{K})|}^2}{|\tilde{S}{(\mathbf{0},\mathbf{K})|}^2},with\tilde{S}(\mathbf{q},\mathbf{K})=\int d^3\mathbf{r}S(\mathbf{r},\mathbf{K})e^{i\mathbf{q}\mathbf{r}},$$

(6)

thus $\tilde{S}$ being the Fourier transform of the source function. In this formula the $\left(0\right)$ superscript denotes the neglection of final state interactions.

Returning to the general, interacting case, if one assumes (according to the core-halo model, see Ref. [8]) that a certain fraction of the particle production (denoted by $\sqrt{\lambda }$) happens in a narrow, few fm diameter region (“core”), and the rest from the decay of long-lived resonances (the contribution of which comes from a much wider region), then one can write the source as

$$S(\mathbf{r},\mathbf{p})=\sqrt{\lambda }{S}_c(\mathbf{r},\mathbf{p})+(1-\sqrt{\lambda })S_h^{R_h}(\mathbf{r},\mathbf{p}),$$

(7)

with a normalization that respects the requirement that $\lambda$ determines the relative weight of the two components:

$$\int d\mathbf{r}S(\mathbf{r},\mathbf{p})=\int d\mathbf{r}{S}_c(\mathbf{r},\mathbf{p})=\int d\mathbf{r}{S}_h(\mathbf{r},\mathbf{p})=1.$$

(8)

Here the indices c and h stand for core and halo, respectively. The $R_h$ “radius” parameter (the characteristic size of the halo part) will be assumed to be much higher than the experimentally resolvable distance, $r_{\max }\approx \hslash /Q_{\min }$, where $Q_{\min }\approx 1-2$ MeV, the minimal mometum difference that can be resolved experimentally.

One can also introduce the core–core, core–halo and halo–halo two-particle source functions as

$$\begin{array}{c}\hfill D(\mathbf{r},\mathbf{K})=\lambda D_{cc}(\mathbf{r},\mathbf{K})+2\sqrt{\lambda }(1-\sqrt{\lambda })D_{ch}(\mathbf{r},\mathbf{K})+{(1-\sqrt{\lambda })}^2{D}_{hh}(\mathbf{r},\mathbf{K}),\end{array}$$

(9)

where the following obvious definitions were used:

$$D_{AB}(\mathbf{r},\mathbf{K})\equiv \int d^3\mathbf{R}{S}_A\left(\mathbf{R}+\frac{\mathbf{r}}{2},\mathbf{K}\right)S_B\left(\mathbf{R}-\frac{\mathbf{r}}{2},\mathbf{K}\right),\mathrm{for}A,B=c\mathrm{or}h.$$

(10)

With yet another notation we can write the terms of $D(\mathbf{r},\mathbf{K})$ as

$$\begin{array}{c}\hfill D(\mathbf{r},\mathbf{K})=\lambda D_{cc}(\mathbf{r},\mathbf{K})+(1-\lambda )D_{\left(h\right)}(\mathbf{r},\mathbf{K}),\mathrm{with}{D}_{\left(h\right)}=\frac{2\sqrt{\lambda }(1-\sqrt{\lambda })D_{ch}+{(1-\sqrt{\lambda })}^2{D}_{hh}}{1-\lambda },\end{array}$$

(11)

where thus the $D_{\left(h\right)}$ term contains all the halo contributions, and $D_{cc}$ is just the core-core component. Perhaps it is useful to explicitly state the (evident) normalization conditions of all these two-particle functions:

$$\int D(\mathbf{r},\mathbf{K})d^3\mathbf{r}=\int D_{cc}(\mathbf{r},\mathbf{K})d^3\mathbf{r}=\int D_{ch}(\mathbf{r},\mathbf{K})d^3\mathbf{r}=\int D_{hh}(\mathbf{r},\mathbf{K})d^3\mathbf{r}=\int D_{\left(h\right)}(\mathbf{r},\mathbf{K})d^3\mathbf{r}=1.$$

(12)

Using these definitions, the correlation function can be expressed as

$$C_2(\mathbf{q},\mathbf{K})\approx \lambda \int d^3\mathbf{r}{D}_{cc}(\mathbf{r},\mathbf{K})|{\psi }_{\mathbf{q}}^{\left(2\right)}{\left(\mathbf{r}\right)|}^2+(1-\lambda )\int d^3\mathbf{r}{D}_{\left(h\right)}(\mathbf{r},\mathbf{K}){\left|{\psi }_{\mathbf{q}}^{\left(2\right)}\left(\mathbf{r}\right)\right|}^2.$$

(13)

By taking the $R_h\to \infty$ limit in the second term2, one arrives at the well-known Bowler–Sinyukov formula [1,2] as

$$C_2(\mathbf{q},\mathbf{K})=1-\lambda +\lambda \int d^3\mathbf{r}{D}_{cc}(\mathbf{r},\mathbf{K}){\left|{\psi }_{\mathbf{q}}^{\left(2\right)}\left(\mathbf{r}\right)\right|}^2.$$

(14)

Specifically, in the free case (with plane-wave wave functions) one arrives at the

$$C_2^{\left(0\right)}(\mathbf{q},\mathbf{K})=1+\lambda \frac{|{\tilde{S}}_c{(\mathbf{q},\mathbf{K})|}^2}{|{\tilde{S}}_c{(\mathbf{0},\mathbf{K})|}^2},$$

(15)

formula (including the normalization term, which is unity in this paper). The experimental observation is that—although the free correlation function defined in Equation (6) takes the value of 2 at 0 relative momentum: $C_2^{\left(0\right)}(\mathbf{0},\mathbf{K})=2$—the measured value is $1+\lambda$. The core-halo model thus naturally explains this fact in terms of the finite momentum resolution of any experiment. In the core-halo model the intercept of the real, measurable correlation function at $\mathbf{q}=0$ thus tells the fraction of pions coming from the core. In the Coulomb interacting (realistic) case, the interpretation of $\lambda$ as any intercept parameter is not so simple, however. The Bowler–Sinyukov method, Equation (15) gives a means to take the core-halo model into account when treating the Coulomb effect.

To investigate the $\lambda$ parameter (which, as it is directly connected to the proportion of resonance decay particles, may have interesting physical consequences, see e.g., Refs. [5,9]) one needs a firm grasp on the effect of the final state interactions in Bose–Einstein correlation functions. For the most important such effect, the Coulomb effect, the ${\psi }_{\mathbf{q}}^{\left(2\right)}\left(\mathbf{r}\right)$ wave function (the two-body scattering solution of the Schrödinger equation with Coulomb repulsion) is well known in the center-of-mass system of the outgoing particles (the so-called PCMS system). Its expression is

$${\psi }_{\mathbf{q}}^{\left(2\right)}\left(\mathbf{r}\right)=\frac{1}{\sqrt{2}}\frac{\Gamma (1+i\eta )}{e^{\pi \eta /2}}\left\{e^{i\mathbf{k}\mathbf{r}}F\left(-i\eta ,1,i(kr-\mathbf{k}\mathbf{r})\right)+[\mathbf{r}\leftrightarrow -\mathbf{r}]\right\},\mathrm{where}\mathbf{k}=\frac{\mathbf{q}}{2}.$$

(16)

Here $F(·,·,·)$ is the confluent hypergeometric function, $\Gamma (·)$ is the Gamma function, and

$$\begin{array}{c}\hfill \eta =\frac{q_e^2}{4\pi {\epsilon }_0}\frac{1}{\hslash c}\frac{m_{\pi }{c}^2}{qc}={\alpha }_{\mathrm{EM}}\frac{m_{\pi }c}{q}\end{array}$$

(17)

is the Sommerfeld parameter, with $q_e^2/\left(4\pi {\epsilon }_0\right)$ being the Coulomb-constant, ${\alpha }_{\mathrm{EM}}$ the fine-stucture constant of the electromagnetic interaction, and $m_{\pi }$ the pion mass (as from now on, we restrict this analysis to pion pairs).

For a given source function $S(\mathbf{r},\mathbf{K})$, the ratio of the (measurable) correlation function $C_2\left(\mathbf{q}\right)$ and the $C_2^{\left(0\right)}\left(\mathbf{q}\right)$ function is usually called the Coulomb correction3, $K\left(\mathbf{q}\right)$:

$$K\left(\mathbf{q}\right)=\frac{C_2\left(\mathbf{q}\right)}{C_2^{\left(0\right)}\left(\mathbf{q}\right)}\Rightarrow C_2^{\left(0\right)}\left(\mathbf{q}\right)=C_2\left(\mathbf{q}\right)·\frac{1}{K\left(\mathbf{q}\right)}.$$

(18)

If one focuses on the simple property of the $C_2^{\left(0\right)}\left(\mathbf{q}\right)$ function as being the Fourier transform of the source, then one might want to recover $C_2^{\left(0\right)}\left(\mathbf{q}\right)$ from the measured $C_2\left(\mathbf{q}\right)$: for this, one uses the Coulomb correction factor. Indeed, many assumptions have been used to estimate the $K\left(\mathbf{q}\right)$ factor: the simplest case is the so-called Gamow factor that treats the source as a point-like one when calculating $K\left(\mathbf{q}\right)$:

$$S\left(\mathbf{r}\right)={\delta }^{\left(3\right)}\left(\mathbf{r}\right)\Rightarrow K\left(\mathbf{q}\right)=K_{\mathrm{Gamow}}\left(q\right)={\left|{\psi }_{\mathbf{q}}^{\left(2\right)}\left(0\right)\right|}^2=\frac{2\pi \eta }{e^{2\pi \eta }-1}.$$

(19)

A method that suits the scope of heavy-ion collisions a little more would be to pre-calculate $K\left(\mathbf{q}\right)$ for a single specific given assumption for $S\left(\mathbf{r}\right)$, then apply this correction (with the Bowler–Sinyukov method) and find the $S\left(\mathbf{r}\right)$ from a fit to the Fourier transform of the recovered $C^{\left(0\right)}\left(\mathbf{q}\right)$. However, it is clear that this process should be done iteratively: after the first “round” of such fits, one would have to re-calculate the Coulomb correction. When this iteration converges, one in principle arrives at the proper $S\left(\mathbf{r}\right)$.

3. Numerical Table for the Coulomb Correction for Lévy Source

Recent studies have shown that the assumption of a Lévy-type of source function is well suited for the description of two-particle Bose–Einstein correlation functions. The details of the validity of the Lévy-shape assumption is exhaustively expounded in Refs. [4,5]. The (spherically symmetric) Lévy distribution utilized here has two parameters, scale parameter (radius) R and Lévy index $\alpha$, and is expressed as

$$\begin{array}{c}\hfill \mathcal{L}(\alpha ,R,\mathbf{r}):=\int \frac{d^3\mathbf{q}}{{\left(2\pi \right)}^3}{e}^{i\mathbf{q}\mathbf{r}}exp\left(-\frac{1}{2}{\left|{\mathbf{q}}^2{R}^2\right|}^{\alpha /2}\right).\end{array}$$

(20)

In the $\alpha =2$ case one gets a Gaussian distribution, in the $\alpha =1$ case the Cauchy distribution is recovered. For other $\alpha$ values, no simple analytic expression exists for the result of this Fourier transform-like integral. As a remark, we note that the concept of this symmetric Lévy distribution can be generalized without much effort to the non-spherically symmetric case by replacing $R^2$ with a symmetric $3\times 3$ matrix $R_{kl}^2$.

In order to apply Lévy-type sources in a self-consistent way, the Coulomb integral defined in Equation (14) has to be calculated. This cannot be carried out in a straightforward analytic manner. In the following we demonstrate two approaches that can be employed to handle the Coulomb final state effect in the presence of a Lévy source.

The integral in Equation (14) cannot be evaluate analytically for a Lévy source so it has to be calculated numerically. For experimental purposes, the results can be loaded to a binary file as a lookup table and can be used in the fitting procedure (thus circumventing the need for an iterative process for the Coulomb correction). Interpolation also should be applied since the correlation function only can be filled into the lookup table for discrete values of the parameters. This interpolation, however, could cause numerical fluctuations in the ${\chi }^2$ landscape and could mislead the fit algorithm, so an iterative procedure should be applied in the following manner:

• Fit with the function defined in Equation (14) ⇒${\alpha }_0,R_0,{\lambda }_0$,
• Fit with $C_2^{\left(0\right)}(\lambda ,R,\alpha ;Q)\frac{C_2({\lambda }_0,R_0,{\alpha }_0;Q)}{C_2^{\left(0\right)}({\lambda }_0,R_0,{\alpha }_0;Q)}$⇒${\lambda }_1,R_1,{\alpha }_1$,
• Repeat while ${\lambda }_1,R_1,{\alpha }_1$ and ${\lambda }_0,R_0,{\alpha }_0$ differ less then 1%.

In this manner, the fit parameters $\lambda \left(K\right),R\left(K\right),\alpha \left(K\right)$ can be obtained, similarly to Ref. [5].

4. Parametrization of the Coulomb Correction for Lévy Source

In this section, let us review a different approach, where based on the numerical table mentioned above, a parametrization can be formulated. In other words, one can get the Coulomb correction values from the table and parametrize its R and $\alpha$ dependences. This approach was encouraged by the successful parametrization of the $\alpha =1$ case (the Cauchy case) done by the CMS collaboration (see Ref. [10], Equation (5) for details). This can be considered as our starting point for the more general, Lévy case (for arbitrary $\alpha$). The expression used by CMS for the Cauchy distribution, $\alpha =1$ was

$$K{\left(\mathbf{q}\right)}_{\mathrm{Cauchy}}=K_{\mathrm{Gamow}}\left(\mathbf{q}\right)\times \left(1+\frac{{\alpha }_{\mathrm{EM}}\pi m_{\pi }R}{1.26\hslash c+qR}\right),\mathrm{where}{\alpha }_{\mathrm{EM}}=\frac{q_e^2}{4\pi {\epsilon }_0}\frac{1}{\hslash c}\approx \frac{1}{137}.$$

(21)

Generally, this is a correction of the Gamow correction. This simple formula has the advantage of having only 1 numerical constant parameter (the 1.26 in the denominator). However, it assumes $\alpha =1$, and we look for a generalization for arbitrary Lévy $\alpha$ values.

A more general correction for the Gamow correction which is able to describe the Coulomb correction for a Lévy source has to fulfill the following requirements:

• It should follow not only the R, but the $\alpha$ dependence.
• In $\alpha =1$ case, it should reduce to Equation (21).

To fulfill these, we replace R with $R/\alpha$ to introduce the $\alpha$-dependence and take higher order terms in $\frac{qR}{\alpha \hslash c}$ into consideration. Our trial formula is then assumed to be

$$\begin{array}{cc}\hfill K_{\mathrm{Lévy}}(q,\alpha ,R)& =K_{\mathrm{Gamow}}\left(\mathbf{q}\right)\times K_{\mathrm{mod}}\left(\mathbf{q}\right),\mathrm{with}\hfill \\ \hfill K_{\mathrm{mod}}\left(\mathbf{q}\right)& =1+\frac{A(\alpha ,R)\frac{{\alpha }_{\mathrm{EM}}\pi m_{\pi }R}{\alpha \hslash c}}{1+B(\alpha ,R)\frac{qR}{\alpha \hslash c}+C(\alpha ,R){\left(\frac{qR}{\alpha \hslash c}\right)}^2+D(\alpha ,R){\left(\frac{qR}{\alpha \hslash c}\right)}^4},\hfill \end{array}$$

(22)

and the task is to find a suitable choice for the $A(\alpha ,R)$, $B(\alpha ,R)$, $C(\alpha ,R)$, $D(\alpha ,R)$ functions that yield an acceptable approximation of the results of the numerical integration (contained in our lookup table). The assumed form seems to be sufficient since it simplifies to Equation (21) if $\alpha =1$ and $C=D=0$, and could follow the observed weak $\alpha$ dependence of the Coulomb integral (see Figure 1).

We fit the above (22) formula to the numerically calculated results for $\alpha$ parameter values between 0.8 and 1.7 and R parameter values between 3 fm and 12 fm, where the ranges were motivated by the results of Pioneering High Energy Nuclear Interaction eXperiment (PHENIX) [5]. With this we obtained the A, B, C, D values as a function of the given $\alpha$ and R parameters. As a next step, we also parametrized these dependencies empirically, and found that the following expressions give satisfactory agreement with the lookup table:

$$\begin{array}{cc}\hfill A(\alpha ,R)& ={(a_A\alpha +a_B)}^2+{(a_{C}R+a_D)}^2+a_E{(\alpha R+1)}^2\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill B(\alpha ,R)& =\frac{1+b_A{R}^{b_B}-{\alpha }^{b_C}}{{\alpha }^{2}R({\alpha }^{b_D}+b_E{R}^{b_F})}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill C(\alpha ,R)& =\frac{c_A+{\alpha }^{c_B}+c_C{R}^{c_D}}{c_E}{\left(\frac{\alpha }{R}\right)}^{c_F}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill D(\alpha ,R)& =d_A+\frac{R^{d_B}+d_C{\alpha }^{d_F}}{R^{d_D}{\alpha }^{d_E}}.\hfill \end{array}$$

(26)

The parameters in these functions turn out best to have the values as follows:

$$\begin{array}{cccccc}{a}_A=0.36060,\hfill & a_B=-0.54508,\hfill & a_C=0.03475,\hfill & a_D=-1.30389,\hfill & a_E=0.00378,\hfill & \\ b_A=2.04017,\hfill & b_B=0.55972,\hfill & b_C=2.47224,\hfill & b_D=-1.26815,\hfill & b_E=-0.11767,\hfill & b_F=0.52738,\\ c_A=-1.00015,\hfill & c_B=0.00012,\hfill & c_C=0.00008,\hfill & c_D=0.26986,\hfill & c_E=0.00003,\hfill & c_F=1.75202,\\ d_A=0.00263,\hfill & d_B=-0.13124,\hfill & d_C=-0.83149,\hfill & d_D=1.57528,\hfill & d_E=0.27568,\hfill & d_F=0.04937.\end{array}$$

This parametrization describes the R and $\alpha$ dependence of the Coulomb integral in a range where the Coulomb correction deviates from 1 by more than a factor of ∼10${}^{-4}$–10${}^{-5}$. We find that this region is 0 GeV/c$<q<$ 0.2 GeV/c. As an example, for $R=6$ fm and with different $\alpha$ values, we plotted the results of the parametrization on Figure 1.

It turns out that the functional form specified above does yield a satisfactory fit at lower values of q, below 0.1–0.2 GeV/c. However, at higher values, the fit that is acceptable at low q, inevitably starts to deviate from the desired values, i.e., cannot be used to extrapolate beyond the fitted q range. The intermediate q region above and around 0.1 GeV/c can instead be described with an exponential-type function parametrized based on intermediate q fits to the numerical table, with the following functional form:

$$E\left(q\right)=1+A(\alpha ,R)exp\left\{-B(\alpha ,R)q\right\},$$

(27)

where the $A(\alpha ,R)$ and $B(\alpha ,R)$ functions have a form as

$$\begin{array}{cc}\hfill A(\alpha ,R)& =A_a+A_b\alpha +A_{c}R+A_d\alpha R+A_e{R}^2+A_f{\left(\alpha R\right)}^2,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill B(\alpha ,R)& =B_a+B_b\alpha +B_{c}R+B_d\alpha R+B_e{R}^2+B_f{\left(\alpha R\right)}^2.\hfill \end{array}$$

(29)

The parameters were chosen based on a fit to numerically calculated Coulomb correction values, and the optimal case was found to be represented by these parameter values:

$$\begin{array}{cccccc}{A}_a=0.20737,\hfill & A_b=-0.00999,\hfill & A_c=-0.02671,\hfill & A_d=-0.00373,\hfill & A_e=0.00119,\hfill & A_f=0.00016,\\ B_a=25.80500,\hfill & B_b=4.01674,\hfill & B_c=0.00873,\hfill & B_d=-0.25606,\hfill & B_e=0.01077,\hfill & B_f=-0.00270.\end{array}$$

The exponential damping factor of Equation (27) is “joined” to the proper parametrization valid for the interesting q range by a Wood–Saxon-type of cut-off function:

$$F\left(q\right)=\frac{1}{1+exp\left(\frac{q-q_0}{D_q}\right)},$$

(30)

where $q_0=0.08$ GeV and $D_q=0.02$ GeV. We investigated different cut-off functions, such as $1/(1+{(q/q_0)}^n)$, but found that the results are rather independent from this choice.

Putting all of the above together, our final parametrization, valid for $\alpha =$ 0.8–1.7 and $R=$ 3–12 fm values, is thus

$$K{(q,\alpha ,R)}^{-1}=F\left(q\right)\times K_{\mathrm{Gamow}}^{-1}\left(q\right)\times K_{\mathrm{mod}}^{-1}(q;\alpha ,R)+(1-F\left(q\right))\times E\left(q\right)$$

(31)

and the Coulomb corrected correlation function which could be fitted to data, can be written in the form of

$$C_2(q;\alpha ,R)=\left[1-\lambda +K(q;\alpha ,R)\lambda \left(1+exp\left[{\left|qR\right|}^{\alpha }\right]\right)\right]·\left(\mathrm{assumed}\mathrm{background}\right).$$

(32)

We used this formula to reproduce PHENIX results Ref. [5] 4; this can be seen on Figure 2. The two fits are compatible with each other. For an example code calculating the formula of (31), please see Ref. [11]. Example curves resulting from the above (32) formula (with the background being unity) are shown in Figure 3. These clearly show how R changes the scale, and $\alpha$ changes the shape of the correlation functions. Parameter $\lambda$ provides an overall normalization to the distance of these curves from unity, as described by Equation (32).

We investigated the parametrization by means of its relative deviation from the lookup table. The results can be seen in Figure 4. In the case when $\alpha =1.2$ with different R values, we present a two-dimensional histogram of the relative differences in in Figure 5. The maximum of these relative differences is around 0.05%.

5. Conclusions

We investigated the Coulomb correction of Bose–Einstein correlations in high energy heavy ion reactions under the assumption of Lévy source functions. We outlined two equivalent methods that are suited for an experimental analysis. One of them is a numerical lookup table, another one is a parametrization obtained from the former. We investigated the accuracy of the methods and found that a not very complicated ad-hoc parametrization, in the well defined parameter range of $R=$ 3–12 fm and $\alpha =$ 0.8–1.7, provides an experimentally acceptable description of the results of the numerical integration that is required for the handling of the Coulomb effect. Our parametrization can thus be used effectively in Bose–Einstein correlation analyses that assume Lévy-type source functions.