Analytic saddlepoint approximation for ionization energy loss distributions

TL;DR

This paper introduces a saddlepoint approximation for ionization energy loss distributions applicable across a range of relativistic velocities, providing a computationally efficient method that closely matches detailed models.

Contribution

A novel saddlepoint approximation method for ionization energy loss distributions that is accurate and computationally efficient across various relativistic velocities.

Findings

Accurately reproduces Vavilov's most probable energy loss within 1%.

Matches the full width at half maximum within 10%.

Performs better as Vavilov's kappa approaches 1.

Abstract

We present a saddlepoint approximation for ionization energy loss distributions, valid at arbitrary relativistic velocities of the incident particle 0<v/c<1, provided that ionizing collisions are still the dominant energy loss mechanism. We derive a closed form solution closely related to Moyal's distribution. This distribution is intended for use in simulations with relatively low computational overhead. The approximation generally reproduces the Vavilov most probable energy loss and full width at half maximum to better than 1% and 10%, respectively, with significantly better agreement as Vavilov's kappa approaches 1.

Equations

$${\left(\frac{d\sigma }{d\epsilon }\right)}_R=\frac{2\pi z^2{e}^4}{m_e{v}^2{\epsilon }^2},$$

$${\left(\frac{d\sigma }{d\epsilon }\right)}_M={\left(\frac{d\sigma }{d\epsilon }\right)}_R[1-{\beta }^2\epsilon /T_m],$$

$$T_m=2m_e\frac{(T+M{)}^2-M^2}{M^2+2(T+M)m_e+m_e^2}.$$

$$-⟨\frac{dE}{dx}⟩=\frac{KZ\rho }{A{\beta }^2}\left[\ln \frac{\sqrt{2m_e{\beta }^2{\gamma }^2{T}_m}}{I}-{\beta }^2\right].$$

$$N=$$

$$=$$

$$x_0=\frac{1}{{\sigma }_T}\frac{A}{Z{N}_A\rho }.$$

$$-⟨\frac{dE}{dx}⟩=\frac{KZ\rho }{2A{\beta }^2}\left[\ln \left(\frac{T_m}{I^{\prime }}\right)-{\beta }^2\right].$$

$$I^{\prime }=I{e}^{{\beta }^2}\frac{I}{2m_e{\beta }^2{\gamma }^2}.$$

$${\sigma }^2(x)=$$

$$=$$

$$\approx$$

$$F(\Delta ,N)=\frac{1}{2\pi i}{\int }_{s_0-i\infty }^{s_0+i\infty }ds\exp [s\Delta +NM(s)],$$

$$M(s)={\int }_1^{{\epsilon }_m}d\epsilon \frac{d\sigma }{d\epsilon }{e}^{-s\epsilon }-1.$$

$$\Delta =-N{M}^{\prime }(s).$$

$$F(\Delta ,N)\approx \frac{\exp \{N[M(s)-s{M}^{\prime }(s)]\}}{\sqrt{2\pi N{M}^{\prime \prime }(s)}}.$$

$$M(s)=R(s)+\frac{{\beta }^2}{{\epsilon }_m}{R}^{\prime }(s).$$

$$R(s)=$$

$$=$$

$$R^{\prime }(s)=-{\int }_1^{{\epsilon }_m}d\epsilon \frac{e^{-s\epsilon }}{\epsilon },$$

$$R^{\prime \prime }(s)={\int }_1^{{\epsilon }_m}d\epsilon e^{-s{\epsilon }_m}=\frac{e^{-s}-e^{-s{\epsilon }_m}}{s}.$$

$$R^{(3)}(s)=-\frac{1}{s}\left[R^{\prime \prime }(s)+e^{-s}-{\epsilon }_m{e}^{-s{\epsilon }_m}\right],$$

$$R^{(4)}(s)=-\frac{1}{s}\left[2R^{(3)}(s)-e^{-s}+{\epsilon }_m^2{e}^{-s{\epsilon }_m}\right].$$

$$R^{\prime \prime }(s)\approx \frac{1}{s}-1+\mathcal{O}(s),$$

$$R^{\prime }(s)=\ln (s)+c.$$

$$\Delta =-N\ln (s)-Nc,$$

$$s=e^{-\Delta /N+c},$$

$$F(\Delta ,N)\propto \sqrt{s}{e}^{-Ns}.$$

$${\Delta }_p=N\ln (2N)-Nc$$

$$\frac{dF}{d\omega }=\frac{1}{\sqrt{2\pi }}\exp \left\{-\frac{1}{2}[\omega +\exp (-\omega )]\right\},$$

$${\Delta }_p=\xi \left(\ln \frac{\xi }{I^{\prime }}+\ln 2-c\right),$$

$$\xi =x\frac{KZ\rho }{2A{\beta }^2}=x\frac{2\pi e^4}{m_e{\beta }^2}\frac{\rho N_{A}Z}{A}.$$

Full text

††journal: Nuclear Instruments and Methods in Physics Research Section B

1 Introduction

The motivation for this study is to provide a fast approximate method to generate energy loss distributions for heavy charged particles passing through thick, dense materials. The primary goal is to improve the fidelity of our simulations for proton radiography [1] by adding straggling distributions as a function of energy loss, without adding computational time penalty. While there is a large body of work on this subject, much of which would be suitable for our application in one form or another, the main advantage of this study is a fast, compact solution which is suitable over a large range of energy losses for arbitrary values of $v/c$ for the incident particle.

Our starting point is the approximation applied by Moyal [2] to find a closed form approximation to the Landau distribution [3]. It will be further developed to obtain a more general solution applicable to higher incident energies and larger energy losses like Vavilov’s solutions [4]. This further development of the saddlepoint approximation has been pursued elsewhere [5], but our generalization leads to a simpled closed form solution valid for arbitrary relativisitic velocities.

2 Energy loss formulae

The single collision energy loss spectrum for large energy losses comes directly from the Rutherford [6] and Mott [7] differential cross sections in the nonrelativistic and relativistic cases, respectively. Converting the Rutherford cross section to energy loss in cgs units (e.g., [8]), one finds

[TABLE]

where $z$ is the charge of the incident ionizing particle, $v$ is its velocity, $\varepsilon$ is the energy lost in the collision and $m_{e}$ is the mass of the electron. Since we are primarily concerned with protons, the remainder of this article will assume $z=1$. The same equation can be applied to energy loss in collisions with the nuclei in the stopping material with the substitution $m_{e}\rightarrow M(Z,A)$, showing that the energy loss in these interactions is orders of magnitude smaller.

The same conversion using the Mott cross section for relativistic incident particles gives

[TABLE]

where $\beta=v/c$ is the velocity of the incident particle in units of the speed of light and $T_{m}$ is the maximum kinematically allowed energy loss in a single collision. Further corrections dependent on the spin of the incident particle are neglected since they are proportional to $\varepsilon/E$ or $m_{e}/M$, where $M$ is the incident particle’s velocity.

The maximum kinematically allowed energy loss in a single collision, written in terms of the kinetic energy $T=E-M$, is

[TABLE]

This reduces to $T_{m}\approx 2m_{e}\beta^{2}\gamma^{2}$ for $(T+M)m_{e}/M^{2}\ll 1$. This value agrees with the particle data group (PDG) [9], but it uses kinetic energy instead of $\gamma$ and $\beta$ for reference.

The average energy loss is given by the Bethe-Bloch equation:

[TABLE]

The mean excitation energy $I$ accounts for the average effect of all the bound electron states. The constant $K=4\pi N_{A}e^{4}/m_{e}=0.307075$ MeV cm2/mol. For compound materials, a suitable average can be made for $\langle Z/A\rangle$ and $\langle 1/I\rangle$, but we will assume an elemental material for simplicity. We neglect the density correction $\delta$ [10] in the following sections, but it can be included by replacing the second term in the square brackets with $\beta^{2}\rightarrow\beta^{2}+\delta/2$.

3 Energy loss moments

The moments of the total energy loss cross section determine the straggling distributions. When $T_{m}\approx 2m_{e}\beta^{2}\gamma^{2}$ is a good approximation, the energy loss given by Equation 4 is equal to twice the energy loss given by $(d\sigma/d\varepsilon)_{M}$ with the minimum energy loss set to $I$. We will calculate an adjusted value $I^{\prime}$ to account for the difference between the assumed and true energy loss cross sections.

Given the assumed form of the cross section $(d\sigma/d\varepsilon)_{M}$, the mean number of collisions $N$ after passage through a thickness $x$ is

[TABLE]

The last line defines the effective total cross section, $\sigma_{T}$; it is given by $K/2N_{A}\beta^{2}I^{\prime}$, if the term of order $I^{\prime}/T_{m}$ is neglected. Setting $N=1$ and solving for $x$ gives a value for the effective average distance between collisions,

[TABLE]

The mean energy loss will be taken as a given, determined by the Bethe-Bloch equation with any necessary corrections. In terms of $I^{\prime}$ and the assumed Mott form of the cross section it is

[TABLE]

This equation defines $I^{\prime}$ for a given value of $\langle dE/dx\rangle$. Then it is possible to get a precise value for the average energy loss while correctly capturing the high energy loss part of the cross section that is most relevant to straggling distributions with energy losses $\Delta>T_{m}$. The solution for $I^{\prime}$ in terms of Equation 4 is

[TABLE]

For high energy incident particles, the last factor on the right side will be very small and $I^{\prime}$ will be smaller than any physical minimum energy loss. The straggling distributions derived from these assumptions will be too narrow for small energy losses, when compared with realistic cross sections including minimum energy losses and resonant enhancement of the collision cross section with bound electrons. Inclusion of the density correction $\delta$ changes the argument in the exponential to $\beta^{2}+\delta/2$.

The asymptotic width of the straggling distribution in a central limit theorem approximation is

[TABLE]

which agrees with Gaussian approximations presented elsewhere [11, 12]. However, the Gaussian approximation is generally not satisfactory for moderate energy losses satisfying the assumption that $\beta$ and $d\sigma/d\varepsilon$ do not change.

4 Saddlepoint approximation

The objective is to find energy loss distributions for massive, energetic particles passing through thick objects. For the sake of simple, analytic forms the treatment will be approximate. The initial state of interest is a monochromatic beam of particles incident on a finite thickness of material.

We begin with a version of of the Landau-Bothe equation [3, 13, 14] written in the same fashion as Moyal [2],

[TABLE]

in which $\Delta$ is the energy loss, $N$ is the mean number of collisions, and $R(s)$ is a moment generating function given by

[TABLE]

This cross section is normalized so that $F(\Delta,N)$ is written in terms of expected number of collisions instead of the equivalent thickness. The integration limits are based on factorization of the minimum energy loss in a collision, $I^{\prime}$, such that $\varepsilon_{m}=T_{m}/I^{\prime}$.

Equation 14 is a standard Laplace inversion formula and the moment generating function in the exponential encapsulates the form of the cross section and fluctuations in the number of collisions. For given energy loss $\Delta$ and expected number of collisions $N$, the integral for $F(N,\Delta)$ has a saddle point given by

[TABLE]

Taylor expansion of the argument to second order gives a Gaussian integral with amplitude

[TABLE]

The saddle point condition in Equation 16 defines the energy loss as a function of $s$ and $F(\Delta,N)$ gives the amplitude.

The saddle point amplitude requires the moment generating function and its derivatives. The integral in Equation 15 can be written in terms of either the Rutherford or Mott forms of the cross section. If $M(s)$ is used for the relativistic Mott form and $R(s)$ is the moment generating function for the Rutherford cross section, then we have

[TABLE]

Integration by parts gives

[TABLE]

The first derivative is related to the exponential integral,

[TABLE]

while the second derivative takes a simple closed form:

[TABLE]

The third derivative is

[TABLE]

and for completeness the fourth derivative is

[TABLE]

These derivatives are necessary to solve for the most probable energy loss using the relativistic form of the cross section.

5 Moyal’s solution

Moyal gave a closed form solution based on the same limits originally studied by Landau. The relevant limits are $s\ll 1\ll s\varepsilon_{m}$, or equivalently, $N\gg 1$ and $\Delta\ll\varepsilon_{m}$. In this limit, the leading order contribution to the second derivative is

[TABLE]

which we integrate to find

[TABLE]

The precise integration constant is the Euler-Mascheroni constant $c=0.5772\ldots$ The energy loss in terms of $s$ is

[TABLE]

the value of $s$ in terms of the energy loss is

[TABLE]

and the leading dependence of $NR(s)-NsR^{\prime}(s)$ is given by $-Ns$. Putting it all together, the result is

[TABLE]

Setting the derivative equal to zero, the maximum value of $F(\Delta,N)$ is at $s=1/2N$ and the most probable energy loss is

[TABLE]

Using a reduced variable $\omega=(\Delta-\Delta_{p})/N$ to write the solution in terms of the most probable energy loss, the result is Moyal’s distribution,

[TABLE]

where the normalization constant follows from a Gaussian integral upon the substitution $u=\exp(-\omega/2)$.

If we rewrite the variables, the result for the most probable energy loss is

[TABLE]

where $\xi$ has the same definition originally given by Landau

[TABLE]

The peak energy loss grows like $x\ln x$ for small thicknesses; clearly, this behavior must change as $\Delta_{p}$ approaches the average energy loss. The full width at half maximum (FWHM) of Moyal’s distribution is $3.6\xi$, compared to a FWHM of $4\xi$ for Landau’s distribution. This reflects the limitation of the saddle point approximation and the first order treatment of $R^{\prime}(s)$, but this is a worst case scenario in comparison to the Vavilov distribution, as we will show later. Figure 1 is a good illustration of the difference between these two distributions. Although the saddlepoint distribution has been calculated using the generalized form of the next section, it is effectively identical to the Moyal result in this case.

6 Generalized form

A more general form is necessary for higher energy ionizing particles and larger energy losses. We seek a compact closed form which is applicable in such a scenario. Our starting point is an assumption analogous to Equation 26, which satisfies three derivatives with respect to $s$, about the most probable energy loss, or equivalently the most probable value of $s$.

The most probable value of $s$ is found by setting the derivative of Equation 17 to zero. The result is

[TABLE]

In the non-relativistic limit $\beta^{2}=0$, $s_{p}$ varies from $1/2N$ in the Landau-Moyal limit to $1/4N$ in the large energy loss limit. The general result for the minimum value of $s_{p}$ is

[TABLE]

with a limiting minimum value of $s_{p}=1/3N$ as $\beta^{2}\rightarrow 1$.

Our generalization of Equation 26 is

[TABLE]

in which $s_{p}$ is the most probable value of $s$, $m_{1}=M^{\prime}(s_{p})$, $m_{2}=M^{\prime\prime}(s_{p})$, $r=m_{2}/m_{3}$ and $m_{3}=|M^{(3)}(s_{p})|$. This form ensures the correct value of $s_{p}$ and the first three derivatives at $s_{p}$. It also reduces to Moyal’s appproximation for $R^{\prime}(s)$ in the Landau-Moyal limit. When we invert the approximate form for $M^{\prime}(s)$, the first two derivatives in terms of the energy loss are

[TABLE]

$\Delta_{p}$ is the most probable energy loss and this value of $\omega$ also reproduces Moyal’s reduced energy loss variable in the small energy loss limit.

Now we wish to find a distribution analogous to Equation 31 with more general applicability. The more complicated cross section does not allow the same simple reduction of the exponent in the saddlepoint amplitude. Instead we integrate the approximate form for $M^{\prime}(s)$ to obtain $M(s)$,

[TABLE]

which has used the substitution $u=1+(s-s_{p})/r=e^{-\omega}$. It is fine to neglect a potential integration constant since we will renormalize the distribution obtained from the saddlepoint approximation in any case. A simplified result for the exponent is

[TABLE]

The factor of $1/\sqrt{M^{\prime\prime}}$ cancels the last term since it takes the form $e^{-\omega/2}$. The final renormalized distribution is

[TABLE]

where we have defined

[TABLE]

If we integrate the distribution in terms of $te^{-\omega}$, it takes the integral form of the gamma function. This gives the normalization constant. The value of $t$ is 1/2 in the small energy loss limit, leading to Moyal’s result. Larger energy losses lead to larger values of $t$. For values of $t\gg 1$, the Gaussian approximation is recovered if one expands the argument to second order: $\omega+e^{-\omega}\approx 1+\omega^{2}/2$.

7 Efficient computations

The primary difficulty with the outlined approach is to determine $s_{p}$ and $R^{\prime}(s_{p})$. To this end we give simple methods to caclulate these values. The condition for the most probable value of $s$ is

[TABLE]

In this equation we substitute $s_{1}=\textrm{min}(s_{p})$ from Equation 35 and $s_{2}=\textrm{max}(s_{p})=1/2N$ to find $y_{1}(s_{1})=\delta s_{1}$ and $y_{2}(s_{2})=\delta s_{2}$. This is a monotonic function and it is linear in $s$ to an excellent approximation over the range in question. Whichever value of $y$ is closer to zero, we iterate to the solution by setting $s(2)=s(1)-\delta s$. This procedure rapidly converges to the precise value of $s_{p}$ and $\Delta_{p}$.

The higher derivatives have simple closed forms. The precise value of the first derivative of the moment generating function is given by

[TABLE]

where $E_{1}$ is the exponential integral (see, for example, Abramowitz and Stegun [16]). In the Landau-Moyal limit only the second term is necessary. It is useful to have quick solutions for this integral without special math libraries. For values such that $s\varepsilon_{m}>1$, a convenient solution is

[TABLE]

which has truncated terms beginning with $s^{2}/4$ and $-\exp(-s\varepsilon_{m})/(s\varepsilon_{m})^{2}$. A convenient formula for small values of $s\varepsilon_{m}$ is

[TABLE]

In this limit it will usually be acceptable to neglect $s^{n}$.

8 Comparisons with Vavilov distribution

Now we demonstrate the performance of our approximate distribution relative to the well-known Vavilov distribution. The results for the Vavilov distribution were generated using the ROOT package from CERN [17]; in particular, we have used the VavilovAccurate implementation of the programs by Schorr [18].

The saddlepoint approximation reproduces the peak energy loss and distribution in the vicinity of the peak energy loss. The approximation is worse for smaller $N$ and farther from the peak of the distribution. To characterize the differences, we have compared the saddlepoint distributions to the Vavilov distribution as a function of Vavilov’s $\kappa$, which is related to Landau’s $\xi$ by

[TABLE]

The comparisons were made over the range $\kappa\in[0.001,10]$ for the values $\beta^{2}=\left\{0.1,0.5,0.9926\right\}$, corresponding to proton kinetic energies of 50.8 MeV, 389 MeV and 10 GeV. Over this range, the most probable energy loss is reproduced with discrepancies less than 1% in all cases. The overall worst performance of the saddlepoint approximation is at small values of the energy loss or $\kappa$. Figure 1 shows the energy loss of 10 GeV protons passing through about 2 mm of Pb; in this case the difference between the most probable energy losses is 0.15%.

The FWHM illustrates the limitations of the saddlepoint method. In the Moyal-Landau limit, the Landau distribution has FWHM$=4.02\xi$ and the Moyal distribution has FWHM=$3.59\xi$, so the Moyal distribution is 10.7% too narrow by this metric. It is worth noting that the Landau distribution is generally deficient in this limit itself: straggling distributions in this limit are wider due to enhanced cross sections at energy losses comparable with electron binding energies. In all cases, the comparison between the Vavilov and saddlepoint FWHM improves with increasing energy loss. This is shown for three values of $\beta^{2}$ in Figure 2. The difference in the FWHM at small energy losses is characteristic of the Landau and Moyal distributions. The disagreement decreases with further energy loss. At $\kappa=1$ the discrepancy has decreased to $<3\%$ for all three cases; the largest difference is 2.4% for $\beta^{2}=0.1$ at $\kappa=1$.

For completeness, we include two more comparisons. The mean energy losses in these comparisons were obtained using the PDG values for the average excitation energy and no other corrections were included. Since we are interested in proton radiography with kinetic energy 800 MeV at the Los Alamos Neutron Science Center, we include protons at this energy ($\beta^{2}=0.7086$) passing through 1 cm of Cu in Figure 3. The saddlepoint approximation with $t=1.4$ performs well at and above the mean energy loss in this case, which is a result of using Equation 36 to get the moments correct at the peak energy loss. The difference between the peak energy losses is about 40 keV with a peak energy loss of 14 MeV, or less than 0.3%. The FWHM is about 5% smaller for the saddlepoint distribution.

Finally, to show a case with a relatively large energy loss where the approximation performs very well, Figure 4 shows distributions for protons with $\beta^{2}=0.1$ and Vavilov $\kappa=10$. The distributions are indistinguishable on a linear scale. Even for $\kappa=10$, the distribution is not symmetric. The saddlepoint approximation is much better than a symmetric Gaussian. An additional advantage is that saddlepoint approximation has a simpler functional form than a skewed Gaussian.

9 Summary

We have given a closed form saddlepoint approximation to Vavilov’s distributions for the ionization energy loss of energetic charged particles. This generally applicable distribution in Equation 44 take the form of a simple extension to Moyal’s closed form approximation for the Landau distribution. The limit where the distribution takes the Moyal form is where it performs worst in comparison with the Landau-Vavilov forms. Our approximate distribution correctly reproduces the most probable energy loss to better than 1% over orders of magnitude in energy loss. The width of the approximate distribution is about 10% narrower than the Landau distribution in the thin limit and ultimately obtains the same form as the Vavilov distribution. In practice, the approximate distribution is nearly identical to the Vavilov result for $\kappa\geq 1$.

We have described how to compute our distribution. No special functions are required, only exponentials and logarithms. The simple closed form of distribution has advantages. The normalization constant only requires the calculation of $t^{t}/\Gamma(t)$. The cumulative distribution function resulting from Equation 44 can be expressed in terms of common special functions such as the incomplete gamma function and the exponential integral. We find the performance of this approximate method more than adequate for our application to proton radiography. We hope it might find broader application.

Acknowledgement

This work was supported by the Advanced Radiography Science Campaign (C3) at Los Alamos National Laboratory.