A practical solution of the Bethe equation in the energy range applicable to radiotherapy and radionuclide production

While the dose deposition of charged hadrons has received much attention over the last decades starting in 1930 with the publication of the Bethe equation, there are still practical obstacles in implementing it in fields like radiotherapy and isotope production on cyclotrons. This is especially true if the target material consists of non-homogeneous materials, either consisting of a mixture of different elements or experiencing phase changes during irradiation. While Monte-Carlo methods have had great success in describing these more difficult target materials, they come at a computational cost, especially if the problem is time-dependent. This can greatly hinder optimal advancement in therapy and isotope targetry. Here, a regular perturbation method is used to solve the Bethe equation in the limit of small relativistic effects. Particular focus is given to incident energy level relevant to radionuclide production and radiotherapy applications, i.e. 10–200 MeV. We present a series solution for the range and dose distribution in terms of elementary functions, as opposed to special functions which will aid in uptake by practitioners.

where n is the electron density of the target material; e is the electron charge; m e is the electron mass; I is the mean excitation potential of the target material; z is the multiple of the electron charge of the projectile; β ≡ u c / is the ratio of the projectile velocity to the speed of light in vacuum, c; and ε 0 is the vacuum permittivity. The Bethe equation is a result of a first-order Born approximation with the projectile being faster than the bound electrons in the target atoms.
The energy of a proton beam E is related to the mass m p and velocity of a proton via with representative values for common materials given in Table 1. The reduction given by Eqs. 2-5 was advanced by Grimes et al. 13 , which we use directly in this work. Using the continuous-slowing-down-approximation (CSDA), Eq. 4 has been integrated from the projectiles' initial velocity to its resting state using a number of different techniques in order to determine its range R, defined as the root of = u R ( ) 0 8, [14][15][16] . Of note is a recent solution of the Bethe equation advanced by Grimes et al. 13 who show that in the non-relativistic limit ( where E 1 is the exponential integral; u o is the initial velocity of the particles; b is dimensionless and related to the ratio of the incident energy to the mean excitation potential; and the big-O notation is used to describe the error term in the approximation. It indicates that the absolute value of the error is approximately ε 2 . We also note that the grouping on the left-hand side of the equation is dimensionless indicating that the characteristic length  c for this problem is ≡  u A / c o 4 , as R has units of length. The third dimensionless group that governs this problem is ε, the ratio of the particle velocity to the speed of light. This solution indicates that the (dimensionless) range  R/ c may be uniquely expressed as a function of b and ε, and to lowest order is independent of relativistic effects. More www.nature.com/scientificreports www.nature.com/scientificreports/ importantly, this solution methodology may be extended and used to estimate the dose distribution. However, to do so one requires the use of an inverse exponential integral function, a non-standard function, which Grimes et al. 13 indicate that no closed-form expression exists which hinders the utility of this expression. Overcoming this difficulty motivates this work.
In this work, we present a solution of the Bethe equation using a regular perturbation method. Using a solution methodology different to that performed by Grimes et al. 13 , we are able to create a reasonable solution which avoids the use of the exponential integral. We present our solution in the small relativistic limit and then extend the number of terms to cover the energy range up to ~200 MeV for protons. We find that with this novel solution the complexity of using the Bethe equation has been reduced dramatically, without significant loss of precision. Because of this, we anticipate uptake of these expressions by the radioisotope and radiotherapy communities 17 for optimization purposes, for example target design.

Model Derivation
Our methodology is inspired by Grimes et al. 13 who demonstrate that the problem can be approximated by a series solution. Here, we present an alternative methodology in creating the series and begin our solution by scaling Eq. 4 with o c and then switch the roles of the dependent and independent variables to avoid divergence of the velocity gradient near the Bragg peak. This yields If we seek a solution of the form and reduce the complexity of the notation by transforming the dependent variable using In an unmotivated trick, we do so so again and find that With this we postulate that the solution of the problem may be expressed as a series in 1/z o as the magnitude of z o is large. When this integration procedure is repeated, we find that ; hence < − < z z 1 1/ 2/ 0 0 2 , implying a convergent series. This solution may be expressed in a more convenient form by using Eq. 11. Indeed, we find that when truncated after three terms, representing an error of ~1% at 100 MeV. We extend the leading-order solution by repeating the procedure given above and solve Eqs. 13-15 using the same procedure to find that We demonstrate the usefulness of this solution in Fig. 1 where excellent agreement is found when compared to literature values 7 and a numerical solution of the Bethe equation using Romberg integration by applying Richardson extrapolation on a trapezoidal rule for a select number of materials -a gas, liquid, a solid and a complex material; a larger selection of material is shown in Figs. 2, 3 and 4. These calculations were performed using the parameters given in Table 1. At ε O( ) 2 the solution is fairly reasonable to 50 MeV; at ε O( ) 8 the expression yields reasonable estimates to ~200 MeV. The immediate benefit of Eq. 20 is that we are able to estimate the stopping distance without use of the exponential integral. Increasing the number of terms does not improve the fit below 1 MeV as the Barkas and Bloch corrections are not included. Interestingly, we find that the empirical expression 18,19 p 0 behaves similarly to our solution. Here, α is a material dependent constant and p depends on the proton velocity which typical ranges from 1.5 to 2. As can be seen in Fig. 5, in the region < < b 10 1 0 1 3 , R varies only weakly implying that over small energy intervals, R may be taken as constant, yielding ∝ R E 0 2 . This is also evident if we solve the Bethe equation without including the logarithmic term. estimating dose distribution. To estimate the dose distribution u(x), we proceed in an analogous manner to that given above. Here, we examine the solution to the leading order equation, i.e. Eq. 12 and then report the solution to the higher order terms. As Eq. 12 is separable, we are in a position to integrate this expression by parts, repeatedly, to generate a series in 1/z. We abandon the initial condition ( = x z ( ) 0 o ) and in its place, use the stopping condition defining the range,( . By doing so, this dramatically aids in the evaluation of the integral and allows for the development of the series. To highlight this procedure, we rewrite Eq. 12 as x R z z 2 2 which can be integrated to give  www.nature.com/scientificreports www.nature.com/scientificreports/   when truncated after three terms. When this process is repeated and extended to the next term we find (2(ln( )) ln( ) 1) 8(ln( )) (9(ln( )) 3ln( ) 2) 36(ln( )) Finally, we apply the solution methodology to the higher order terms as (120(ln( )) 2ln( ) 1) 512(ln( )) 31(25(ln( )) 5ln( ) 2) 2400(ln( )) (28) A representative comparison is given in Fig. 6 where the numerical and approximate solution are compared to SRIM 8 of the dose distribution of 100 MeV protons impinging in four different materials. The data points are average energy results of protons intersecting the individual depth and the error bars are the corresponding standard deviations of the individual energy spreads. Again, very good agreement is observed between the analytical and our approximate solution and with SRIM simulations.

Robustness of the Solution.
In this work we present a series solution to approximate the solution of the Bethe equation. We have tested the error in this approximation by comparison to a numerical solution, tabulated data, and SRIM calculations. We find excellent agreement over a defined proton energy range applicable to radiotherapy and radionuclide production. One immediate benefit of this solution is that estimate water-equivalent-thickness for radiotherapy applications can be improved, over the empirical Bragg-Kleemann law 19 .
Care must be given in interpretation of our solution. In essence, we advance a simpler solution to the Bethe equation. Although the Bethe equation is routinely used, there are limitations in its use. In practice, there are a number of physical mechanisms which cause deviations from the range predicted by Eq. 6. As is evident in our results, the first type of deviations occurs at low incident energy, where the first-Born approximation reaches its limit. Here we find corrections to the form of the Bethe equation by Barkas or Bloch 6 . In this energy regime the Barkas, Bloch and Shell correction need to be applied. In addition, at high projectile energies, relativistic and density effects and radiative losses need to be taken into account 6 .
In the second category, we find that the distance travelled by the particle along its initial direction is slightly reduced by small deflections caused by scattering. The impact of this is minimal for low Z materials, with the www.nature.com/scientificreports www.nature.com/scientificreports/ 'detour factor' of the projected range greater than 99.87% for protons of energy >100 MeV in water 6 or for metals where the stopping distance is small. Outside of this range, the detour factor becomes substantial, especially for gases 16 where Schlyer et al. 20 shows that the standard deviation of the probability distribution in scattering angle is related in a complex manner to the local energy of the beam.
In the third category, we find a number of works highlighting the effect of a two-way coupling between the beam and the absorbing material on the irradiation profile. In perhaps the most widely cited series of work, Heselius and his co-workers 21,22 demonstrate experimentally symmetry-breaking in the beam profile when density of the gas target is stratified due to localized heating. Indeed, this problem is relatively unexplored on the theoretical front with most authors de-coupling the temperature effects [23][24][25][26] . In perhaps the most comprehensive study, O'Brien and her co-workers 25,27,28 couples Monte-Carlo simulation of the beam profile with Reynolds-Averaged Navier-Stokes (RANS) 29 modelling of the transport equations to estimate the temperature, density and beam profile. Interestingly, with these latter two categories, these deviations are not evident in our result in this range, even with gas.

conclusion
We present a simple and easy to implement approximate solution of the Bethe equation with good agreement in the proton energy range of 10 MeV to approximately 200 MeV, applicable for the energy range of proton radiotherapy and medical isotope production. Although tabulated values are available for uniform materials, this solution methodology dramatically reduces the effort to estimate the range in or the dose distribution for mixed media, i.e. composite bodies or multi-layer materials, or estimating water equivalent thickness. Within solid target design, it is now possible to combine computational fluid dynamics codes with our method to estimate internal heating; instead of coupling these with Monte-Carlo packages. Our approach can also easily be coupled with fluid materials for reduced-order modeling of the interaction caused by irradiation and consequent fluid-dynamic effects. This in turn simplifies and accelerates the solution finding and is therefore a potentially convenient tool to advance radiotherapy and medical isotope-production targetry.