The role of electron capture and energy exchange of positively charged particles passing through matter

The conventional treatment of the Bethe-Bloch equation for protons accounts for electron capture at the end of the projectile track by the small Barkas correction. This is only a possible way for protons, whereas for light and heavier charged nuclei the exchange of energy and charge along the track has to be accounted for by regarding the projectile charge q as a function of the residual energy. This leads to a significant modification of the Bethe-Bloch equation, otherwise the range in a medium is incorrectly determined. The LET in the Bragg peak domain and distal end is significantly influenced by the electron capture. A rather significant result is that in the domain of the Bragg peak the superiority of carbon ions is reduced compared to protons.


Introduction
The application of the Bethe-Bloch equation ( Barkas shell I π π π π ρ ρ ρ ρ β β β β E I is the atomic ionization energy, weighted over all possible transition probabilities of atomic/molecular shells, and q denotes the charge number of the projectile (e.g. proton). The meaning of the correction terms a shell , a Barkas , a 0 and a Bloch are explained in literature, see Bethe (1930), Bloch (1933), Bethe et al (1953), Bethe (1953), ICRU49 (1993), andBoon (1998). Since the Bloch correction a Bloch will be introduced in equation (12), we present, for completeness, the remaining correction terms according to The parameter C, referring to shell corrections, is determined by different models (ICRU49 and references therein). A unique parameterization of C depending on Z, A N , and E I does not exist. It is therefore recommended to select C according to proper domains of validity. It must be noted that several models have been proposed to account for shell transitions. Therefore, the recommendations of ICRU49 have been applied in this work.
The function F ARB in equation (4)  represents a correction of BBE due to the electron capture of the positively charged protons at lower energies in the domain of the Bragg peak and behind leading to a slightly increased range R csda , whereas the negatively charged anti-protons cannot capture electrons from the environmental electrons. Therefore their range is slightly smaller. With regard to protons this kind of correction works, i.e. the charge q 2 = 1 is assumed along the total proton track, whereas for charged ions such as He or C 6 is appears to be insufficient to keep the nuclear charge constant along the total track and to restrict the electron capture only to the small Barkas correction (Barkas et al 1963). This means that all positively charged projectile particles stand in permanent exchange of energy E and charge q with environment, and, as a consequence, q 2 is a function of the actual residual energy, i.e. q 2 = q 2 (E), and only for E = E 0 (initial energy) q 2 = q 0 2 is valid. A correct modification of BBE by accounting for q 2 (E) makes the Barkas correction superfluous.
A further critical aspect of BBE, which leads to a modification by accounting for q 2 (E) is the range R csda of the electronic stopping power. Thus a naïve application of BBE would lead to the conclusion that a carbon ion would require the initial energy per nucleon E 0 (carbon ion) = 3 x E 0 (proton), since the square of the carbon charge amounts to 36 and the nuclear mass unit is 12 x nuclear mass unit of the proton. However, the ratio is not 3 to obtain the same range R csda , but about 25/12. The Monte-Carlo code GEANT4 assumes an average charge q Average = 5.06 for the simulations of the carbon tracks. This is, however, not satisfactory, since electron capture is a dynamical process. Therefore the range of carbon ions has been subjected to many studies (Betz 1972, Dingfelder et al 1998, Hollmark et al 2004, Hubert et al 1990, Kanai et al 1993, Kusano et al 2007, Martini 2007, Matveev et al 2006, Sigmund et al 2006, Sigmund 2006, Sihver et al 1998, Yarlagadda et al 1978, Zhang and Newhauser 2009, Ziegler et al 1988 due to the increasing importance of carbon ions in radiotherapy. It is also possible to substitute the electron mass m by the reduced mass m ⇒ µ. However, this leads for protons to a rather small correction (i.e., less than 0.1 % for protons). For complex systems E I and some other contributions like a shell and a Barkas can only be approximately calculated by simple quantummechanical models (e.g., harmonic oscillator); the latter terms are often omitted and E I is treated as a fitting parameter, but different values are proposed and used (ICRU49). The restriction to the logarithmic term leads to severe problems, if either v → 0 or 2m v 2 /E I → 1. It should be added that a correct treatment of the electron capture removes the singularity of positively charged ions, since q 2 (E) → 0, if the residual energy E assumes zero.

The integration BBE for protons
In the following, we consider at first the integration of BBE for protons, i.e. we consider the Barkas correction in the conventional way. In previous publications (Ulmer 2007, Ulmer andMatsinos 2010) we have presented an analytical integration of BBE. BBE is the physical base in the transport of protons and electrons.
In order to obtain the integration of BBE, we start with the logarithmic term and perform the substitutions: With the help of substitution (and without any correction terms), BBE leads to the integration: The boundary conditions of the integral are: The general solution is given by the Euler exponential integral function Ei(ξ) with P.V. = principal value: Some details of Ei(ξ) and its power expansions can be found in Abramowitz and Stegun (1970). The critical case ξ = 0 results from E critical = ME I /4m (for water with E I = 75.1 eV, the critical energy E critical amounts to 34.474 keV; for Pb with E I ≈ 800 eV to about 0.4 MeV). Since the logarithmic term derived by Bethe implies the Born approximation, valid only if the transferred energy E transfer >> the energy of shell transitions, the above corrections, exempting the Bloch correction, play a significant role in the environment of the Bragg peak, and the terms a 0 and a shell remove the singularity. With respect to numerical integrations (Monte Carlo), we note that, in the environment of E = E critical , the logarithmic term may become crucial (leading to overflows); rigorous cutoffs circumvent the problem. Therefore, the shell correction is an important feature for low proton energies. In similar fashion, we can take account of the Barkas correction. Since this correction is also important for low proton energies, it is difficult to make a quantitative distinction to the shell correction, and different models exist in the literature implying overall errors up to 2 % (ICRU49). Using the definitions/suggestions of the correction terms according to ICRU49 and the substitutions we obtain:  (9) does not exist, but it can be evaluated via a procedure valid for integral operators (Feynman 1962), which reads for commutative operators: A'+B' is equated to the complete denominator on the right-hand side of equation .The small Barkas correction and the Bloch correction a Bloch (see equations 12,13) are identified with A' and the other (more important) terms with B': We should already point out here that the expansion (10) is also applicable, if the electron capture is accounted for. The integration of equation with the help of relation (10) leads to standard tasks (i.e. to a series of usual exponential functions). In the following, we add the Bloch correction to the denominator of equation. In order to use the procedure, we define now the non-relativistic energy E nr by: E nr = 0.5·Mv 2 and write the relativistic energy expression E rel (the rest energy Mc 2 is omitted) in terms of an expansion:  (12) The integration of equation is carried out with the boundary conditions. Since these conditions are defined by logarithmic values, which have to be inserted to an exponential function series, the result yields a power expansion for R CSDA in terms of E 0 : The coefficients α n are determined by the integration procedure and only depend on the parameters of the BBE. For applications to therapeutic protons, i.e., E 0 < 300 MeV, a restriction to N = 4 provides excellent results ( Figure 1). For water, we have to take E I = 75.1 eV, Z/A N = 10/18, ρ = 1 g/cm 3 ; Formula becomes: The values of the parameters of Formulas with restriction to N = 4 are displayed in Tables 1 and 2.   (17), if E 0 is in MeV, E I in eV and R CSDA in cm.
a 1 a 2 a 3 a 4 6.94656·10 -3 8.13116·10 -4 -1.21068·10 -6 1.053·10 -9 The determination of A N and Z is not a problem in case of atoms or molecules, where weight factors can be introduced according to the Bragg rule; for tissue heterogeneities, it is already a difficult task. Much more difficult is the accurate determination of E I , which results from transition probabilities of all atomic/molecular states to the continuum (δ-electrons). Thus, according to the report ICRU49 of stopping powers of protons in different media, there are sometimes different values of E I proposed (e.g., for Pb: E I = 820 eV and E I = 779 eV). If we use the average (i.e., E I = 800.5 eV), the above formula provides a mean standard deviation of 0.27 % referred to stopping-power data in ICRU49, whereas for E I = 820 eV or E I = 779 eV we obtain 0.35 % -0.4 %. If we apply the above formula to data of other elements listed in ICRU49, the mean standard deviations also amount to about 0.2 % -0.4 %. . A Gompertz-function is defined by: By inserting the integration boundaries u = 2·ln·4m·E 0 /(M·E I ), i.e., E = E 0 and u → ∞ (E = 0), the integration leads to a sequence of exponential functions; the power expansion is replaced by: For therapeutic protons, the restriction to N = 2 provides the same accuracy ( Figure 2) as formula (16); the parameters are given in Table 3 (a 1 is the same as in Table 2). In the following, we will verify that the latter formula provides some advantages with respect to the inversion E 0 = E 0 (R CSDA ).  (16)) and two exponential functions (equation (17a)).

The Inversion problem: calculation of E 0 (R CSDA ) and E(z)
Above formulas can also be used for the calculation of the residual distance R CSDA -z, relating to the residual energy E(z); we have only to perform the substitutions R CSDA → R CSDA -z and E 0 → E(z) in these formulas. In various problems, the determination of E 0 or E(z) as a function of R CSDA or R CSDA -z is an essential task. The power expansion implies again a corresponding series E 0 = E 0 (R CSDA ) in terms of The coefficients c k are calculated by a recursive procedure; we have given the first three terms in formula (18). Due to the small value of a 1 = 6.8469·10 -4 , this series is ill-posed, since there is no possibility to break off the expansion; it is divergent and the signs of the coefficients c k are alternating, see Abramowitz and Stegun (1970). The inversion procedure of equation leads to the formula: The inverse formula of equation (17a) For therapeutic protons, a very high precision is obtained by the restriction to N = 5 (Table 4 and Figure   4). Formula (20) is also suggested by a plot S(R CSDA ) = E 0 (R CSDA )/R CSDA according to equation (17a). This plot is shown in Figure 3 and gives rise for an expansion of S(R CSDA ) in terms of exponential functions. This plot is obtained by an interchange of the plot E 0 versus R CSDA and a calculation according to Relation.  One way to obtain the inversion Formula is to find S(R CSDA ) by a sum of exponential functions with the help of a fitting procedure. Thus it turned out that the restriction to five exponential functions is absolutely sufficient and yields a very high accuracy. A more rigorous way (mathematically) has been described in the LR of Ulmer (2007). The residual energy E(z), appearing in equation (20), is the desired analytical base for all calculations of stopping power and comparisons with GEANT4. The stopping power is determined by dE(z)/dz and yields the following expression: The aforementioned restriction to N = 5 is certainly extended to equation which can be considered as a representation of the BBE in terms of the residual energy E(z). Due to the low-energy corrections (a 0 , a shell , a Barkas ) the energy-transfer function dE(z)/dz remains finite for all z (i.e., 0 ≤ z ≤ R CSDA ). This is, for instance, not true for the corresponding results according to Formulas at z = R CSDA . The calculation of E(z) and dE/dz according to equations, referred to as LET, is presented in Figure 5. The figure shows that, within the framework of CSDA, the LET of protons is rather small, except at the distal end of the proton track. A change from the interacting reference medium water to any other medium can be carried out by the calculation of R CSDA , where the substitutions have to be performed and used in equations: It is also possible to apply Formulas in a stepwise manner (e.g., voxels of CT). This procedure will not be dE/dz discussed here, since it requires a correspondence between (Z·ρ/A N ) Medium and information provided by CT, see Schneider et al (1996). With regard to heterogeneous media with only CT data as basis information the application of BBE is a more difficult task.

Qualitative properties of the electron transfer described by BBE and electron capture
According to BBE the energy spectrum of produced by carbon ions should be the same as that produced by protons, and the only difference between protons and carbon ions should be the intensity of the released collision electrons, i.e. the amplification factor should be 36 for carbon ions. It is well-known that this property is not valid for the following reasons: The average ionization energy for carbon ions turned out to be E I = 80 eV instead of E I = 75 eV for protons (ICRU49, Paul 2007). The paper of Paul (2007) is based on investigations of some other authors (Bichsel et al 2000, Dingfelder et al 1998, Kraft 200, Krämer et al 2000, Sigmund 1997)The second reason is the electron capture of the carbon ion. Thus a carbon ion can capture a free electron, which has been excited immediately before. Figure 7 shows this effect. However, only electrons with a slow relative velocity to the carbon ion can account for this process (v relative about 0). Since the transition time of the capture electron to a lower atomic state of the carbon ion is less than 10 -10 sec with a simultaneous emission of light (UV or visible), it is possible that the captured electrons goes lost again, and only a stripping effect occurs for a short time. If the C 6+ ions has been finally transferred to a stable C 5+ ion, the identical process can be repeated until at the end track a neutral carbon atom is obtained having only a thermal energy. In the environment of the Bragg peak the effective charge of the carbon ion is about the same that of a proton, namely +e 0 . Since the electron capture can only occur for electrons of which the relative velocity is slow, the upper energy limit of the energy exchange E ex is the Fermi edge E F , which is for an electron gas not higher than the thermal energy k B T. If the charge of carbon ion amounts to +6·e 0 and, at least, > +e 0 , the environmental atomic electrons suffer lowering of the energy levels due to the Coulomb interaction, which leads to an increase of E I . Therefore the stated value of E I = 80 eV represents an average value produced the fast carbon ion starting with +6·e 0 and ending with an uncharged, neutral carbon atom.

Application Fermi-Dirac statistics to electron capture
In the following it is the task to obtain a quantum-statistical description of electron capture and stripping of electrons, i.e. those electrons which reduce the effective charge of the carbon ion for a short time and go lost before a transition to a stable atomic state of carbon can occur. For this purpose we consider the quantum statistical energy exchange E ex between projectile particle such as proton, He ion or carbon ion.
The related mathematical procedure can be used to describe processes like energy straggling, lateral scatter and energy/charge exchange between projectile ion and released electrons below the Fermi edge E F . However, before we can account for the latter problem we have to consider the related mathematical tools.
In general, if H represents the Hamiltonian (either non-relativistic or relativistic) and f(H) an operator functions, then for continuous operators H the connection holds: At first we apply this relation in the non-relativistic case to derive the Gaussian convolution for the description of energy straggling. If the stopping power S(z) = dE(z)/dz of protons is calculated by BBE or by phenomenological equations (13, 22) based on classical energy dissipation, then the energy fluctuations are usually accounted for by: Qualitative figure of projectile interaction of a charged particle BBE:

Electron capture by carbon ion
Carbon ion electron This kernel may either be established by non-relativistic transport theory (Boltzmann equation) or, as we prefer here, by a quantum statistical derivation. Let ϕ be a distribution function and Φ a source function, mutually connected by the operator F H (operator notation of a canonical ensemble): It must be noted that the operator equation (26) was formally introduced (Ulmer 2010) to obtain a Gaussian convolution as Green's function and to derive the inverse convolution. F H may formally be expanded in the same fashion as the usual exponential function exp(ξ); ξ may either be a real or complex number. This expansion is referred to as Lie series of an operator function. Only in the thermal limit (equilibrium), can we write E ex = k B T, where k B is the Boltzmann constant and T is the temperature. This equation can be solved by the spectral theorem provided by the discipline 'functional analysis': π π π σ σ σ σ π π π π γ γ γ γ π π π π γ γ γ γ π π π π σ σ σ σ π π π π z u z u K The kernel K according to equation (27) may either be established by non-relativistic transport theory (Boltzmann equation) or, as we prefer here, by a quantum statistical derivation. It is a noteworthy result (Ulmer 2007) that a quantum stochastic partition function leads to a Gaussian kernel as a Green's function, which results from a Boltzmann distribution function and a non-relativistic exchange Hamiltonian H. An operator formulation of a canonical ensemble is obtained by the following way: let ϕ be a distribution (or output/image) function and Φ a source function, which are mutually connected by the operator. In a 3D version, linear combinations of K(σ, u -x) and the inverse kernel K -1 are also used in scatter problems of photons (Ulmer 2010). As an example, we consider the Schrödinger equation of a free electron transferring energy from the projectile to the environment and obeying a Boltzmann distribution function f(H) =exp(-H/E ex ): The above relation provides: In the case of thermal equilibrium, we can replace the exchange energy E ex by k B T.

Dirac equation and Fermi-Dirac statistics
With regard to our task the Dirac equation to describe the particle motion is an adequate starting-point: ) 30 ( 4 2 2 2 2 0 1  From the view-point of the many-particle-problem Fermi-Dirac statistics is adequate mean by the notation of operator functions: By that, the above expression assumes the shape: Since 1/cosh(ξ) = sech(ξ) holds, we are able to carry out the following expansion, which is given by the Euler numbers E l (see e.g. Abramowitz and Stegun 1970) . Convergence is only established for ξ ≤ π/2. Therefore we have derived a modified expansion which provides convergence for arbitrary arguments of ξ (Ulmer and Matsinos 2010): The spectral theorem of functional analysis provides: By performing all integrations we obtain the distribution functions in the energy space (equation (39)) and position space (equation (39a)): However, we should like to point out that according to the preceding section this determination is only valid for protons and cannot be applied to heavy ions without a change of the parameters.
The transition to the integration (continuum approach up to second order) provides: An essential result is that we are able to modify the previous formula between initial energy E 0 and the range R csda : Please note that the parameters have slightly to be modified α = 0.0069465598; β = 0.0008132157; γ = -0.00000121069; δ = 0.000000001051.
If N =1 and q eff = 0.995 the above formula is valid for protons. However, it turns out that the determination of the effective charge q eff depends on the initial energy E 0 .This can be verified by Figure 7.

Results
In the following we present results of calculations for protons, He ions and carbon ions; the initial energy amounts to 400 MeV/nucleon. This appears to be a reasonable restriction with regard to therapeutic conditions. Thus Figure 8 shows that at the end of the projectile track all charged ions nearly behave in the same manner.  The succeeding figure 10 presents the decrease of the actual charge of carbon ions in dependence of the initial energy E 0 /nucleon. Thus we can conclude that for residual energies E < 50 MeV/nucleon the behavior of the carbon ions does not depend on the initial energy E 0 . With regard to the therapeutic efficacy the behavior of the LET in the environment of the Bragg peak is very significant. For a comparison, we first regard a previous result (Ulmer and Matsinos 2010) referring to the LET of protons. According to Figure 11 the stopping power of protons at the end track depends significantly on the initial energy E 0 and on the beam-line (energy spectrum at the impinging plane). The electron capture of the proton at the end track is ignored.
However, the previous figure 9 clearly shows that with regard to protons the electron capture only becomes more and more significant, when the actual proton energy is smaller than E = 2 MeV. The electron capture of protons at the end track would make the LET of protons zero independent of the initial energy. The succeeding figure 13 presents E(z) and S(z) = dE(z)/dz of protons and S(z) of carbon ions with taking account for electron capture. The initial proton energy amounts to 270 MeV, whereas the initial carbon ion energy is 400 MeV/nucleon. Most significant is the height of the Bragg peak, which is resulting from the electron capture only a factor 1.7 higher than that of protons. In both cases the csda approach is assumed.
Since protons are much more influenced by energy straggling and scatter, their peak height are reduced again, whereas for carbon ions scatter and energy straggling do not play a very significant role due to the mass factor 12.  With regard to the decrease of fluence of primary carbon ions we have derived some modifications of the corresponding decrease curves for protons. However, it appears not to be appropriate to go into further details. A further aspect is the use of the code GEANT4. Since this Monte-Carlo code represents an open programming package, some suitable additional reaction channels have been introduced.

Discussion
The main purpose of this communication was the derivation of a systematic theory of electron capture of charged particles and the role for the LET. There are purely empirical trials to include charge capture in Monte-Carlo codes. However, it appears that a profound basis for the calculation of q 2 (E), E(z), S(z) and R csda (E 0 ) depending besides the initial energy E 0 also on the nuclear mass number N is required to account for further influences of Bragg curves such as the density of the medium and its nuclear mass/charge A N and Z. The unmodified use of BBE leads to wrong results and the Barkas correction, which does not affect the factor q 2 of BBE, only works for protons or antiprotons, whereas for projectile particles like He or carbon ions this correction cannot be considered as small. The presented theory includes the Barkas effect without any correction model.