Hadron and Ion Ionisation

Method

The class G4hIonisation provides the continuous energy loss due to ionisation and simulates the ‘discrete’ part of the ionisation, that is, \(\delta\)-rays produced by charged hadrons. The class G4ionIonisation is intended for the simulation of energy loss by positive ions with change greater than unit. Inside these classes the following models are used:

G4BetheBlochModel, valid for protons with \(T\) > 2 MeV
G4BraggModel,valid for protons with \(T\) < 2 MeV
G4BraggIonModel, valid for protons with \(T\) < 2 MeV
G4ICRU73QOModel, valid for anti-protons with \(T\) < 2 MeV

The scaling relation (33) is a basic conception for the description of ionisation of heavy charged particles. It is used both in energy loss calculation and in determination of the validity range of models. Namely the \(T_p\) = 2 MeV limit for protons is scaled for a particle with mass \(M_i\) by the ratio of the particle mass to the proton mass \(T_i = T_p M_p/M_i\).

For all ionisation models the value of the maximum energy transferable to a free electron \(T_{max}\) is given by the following relation [ea06]:

(168)\[T_{max} =\frac{2m_ec^2(\gamma^2 -1)}{1+2\gamma (m_e/M)+(m_e/M)^2 },\]

where \(m_e\) is the electron mass and \(M\) is the mass of the incident particle. The method of calculation of the continuous energy loss and the total cross-section are explained below.

Continuous Energy Loss

The integration of (29) leads to the Bethe-Bloch restricted energy loss (\(T < T_{cut}\)) formula [ea06], which is modified taking into account various corrections [Ahl80]:

(169)\[\frac{dE}{dx} = 2 \pi r_e^2 mc^2 n_{el} \frac{z^2}{\beta^2} \left [\ln \left(\frac{2mc^2 \beta^2 \gamma^2 T_{up}} {I^2} \right) - \beta^2 \left( 1 + \frac{T_{up}}{T_{max}} \right) - \delta - \frac{2C_e}{Z} + S + F\right ]\]

where

\[\begin{split}r_e &= \mbox{classical electron radius } = e^2/(4 \pi \epsilon_0 mc^2 ) \\ mc^2 &= \mbox{mass-energy of the electron} \\ n_{el} &= \mbox{electron density in the material} \\ I &= \mbox{mean excitation energy in the material} \\ Z &= \mbox{atomic number of the material} \\ z &= \mbox{charge of the hadron in units of the electron change} \\ \gamma &= E/mc^2 \\ \beta^2 &= 1-(1/\gamma^2) \\ T_{up} &= \min(T_{cut},T_{max}) \\ \delta &= \mbox{density effect function} \\ C_e &= \mbox{shell correction function} \\ S &= \mbox{spin term} = 0 for s=0, \left(\frac{0.5 T_{up}}{E}\right)^2 for s=1/2 \\ E &= \mbox{primary energy} \\ F &= \mbox{high order corrections}\end{split}\]

For spin large that 1/2 the same S term is used in the current model. In a single element the electron density is

\[n_{el} = Z \: n_{at} = Z \: \frac{\mathcal{N}_{av} \rho}{A}\]

(\(\mathcal{N}_{av}\): Avogadro number, \(\rho\): density of the material, \(A\): mass of a mole). In a compound material

\[n_{el} = \sum_i Z_i \: n_{ati} = \sum_i Z_i \: \frac{\mathcal{N}_{av} w_i \rho}{A_i} .\]

\(w_i\) is the proportion by mass of the \(i^{th}\) element, with molar mass \(A_i\).

The mean excitation energy \(I\) for all elements is tabulated according to the NIST recommended values for Geant4 NIST materials, for other materials ICRU recommended values [eal84] are used.

Shell Correction

\(2C_e/Z\) is the so-called shell correction term which accounts for the fact of interaction of atomic electrons with atomic nucleus. This term more visible at low energies and for heavy atoms. The classical expression for the term [eal93] is used

(170)\[C = \sum{C_{\nu}(\theta_{\nu},\eta_{\nu})}, \;\; \nu=K,L,M,..., \; \theta=\frac{J_{\nu}}{\epsilon_{\nu}}, \;\; \eta_{\nu}=\frac{\beta^2}{\alpha^2 Z^2_{\nu}},\]

where \(\alpha\) is the fine structure constant, \(\beta\) is the hadron velocity, \(J_{\nu}\) is the ionisation energy of the shell \(\nu\), \(\epsilon_{\nu}\) is Bohr ionisation energy of the shell \(\nu\), \(Z_{\nu}\) is the effective charge of the shell \(\nu\). First terms \(C_K\) and \(C_L\) can be analytically computed in using an assumption non-relativistic hydrogenic wave functions [Wal52][Wal56]. The results [Kha68] of tabulation of these computations in the interval of parameters \(\eta_\nu\) = 0.005–10 and \(\theta_{\nu}\) = 0.25–0.95 are used directly. For higher values of \(\eta_{\nu}\) the parameterization [Kha68] is applied:

\[C_{\nu} = \frac{K_1}{\eta} + \frac{K_2}{\eta^2} + \frac{K_3}{\eta^3},\]

where coefficients \(K_i\) provide smooth shape of the function. The effective nuclear charge for the \(L\)-shell can be reproduced as \(Z_L = Z - d\), \(d\) is a parameter shown in Table 20.

Table 20 Effective nuclear charge for the \(L\)-shell [eal93].
\(Z\)	3	4	5	6	7	8	9	\(>\)9
\(d\)	1.72	2.09	2.48	2.82	3.16	3.53	3.84	4.15

For outer shells the calculations are not available, so \(L\)-shell parameterization is used and the following scaling relation [eal93][Bic92] is applied:

(171)\[C_{\nu} = V_{\nu}C_L(\theta_L,H_{\nu}\eta_L), \;\; V_{\nu}=\frac{n_{\nu}}{n_L}, \;\; H_{\nu}=\frac{J_{\nu}}{J_L},\]

where \(V_{\nu}\) is a vertical scaling factor proportional to number of electrons at the shell \(n_{\nu}\). The contribution of the shell correction term is about 10% for protons at \(T\) = 2 MeV.

Density Correction

\(\delta\) is a correction term which takes into account the reduction in energy loss due to the so-called density effect. This becomes important at high energies because media have a tendency to become polarized as the incident particle velocity increases. As a consequence, the atoms in a medium can no longer be considered as isolated. To correct for this effect the formulation of Sternheimer [SP71] is used:

\(x\) is a kinetic variable of the particle : \(x = \log_{10}(\gamma \beta) = \ln(\gamma^{2} \beta^{2})/4.606 \), and \(\delta(x)\) is defined by

(172)\[\begin{split}\begin{array}{rll} \mbox{for} & x < x_0 : & \delta(x) = 0 \\ \mbox{for} & x \in [x_0,\ x_1] : & \delta(x) = 4.606 x - C + a(x_1 - x)^m \\ \mbox{for} & x > x_1 : & \delta(x) = 4.606 x - C \end{array}\end{split}\]

where the matter-dependent constants are calculated as follows:

(173)\[\begin{split}h\nu_p &= \mbox{plasma energy of the medium } = \sqrt{4\pi n_{el} r_e^3} mc^2/\alpha = \sqrt{4\pi n_{el} r_e} \hbar c \\ C &= 1 + 2 \ln (I/h\nu_p) \\ x_a &= C/4.606 \\ a &= 4.606(x_a - x_0)/(x_1 - x_0)^m \\ m &= 3 .\end{split}\]

For condensed media

\[\begin{split}\begin{array}{ll} I < 100 \: \mbox{eV} & \left \{ \begin{array}{rll} \mbox{for } C \leq 3.681 & x_0 = 0.2 & x_1 = 2 \\ \mbox{for } C > 3.681 & x_0 = 0.326 C - 1.0 & x_1 = 2 \end{array} \right . \\ I \geq 100 \: \mbox{eV} & \left \{ \begin{array}{rll} \mbox{for } C \leq 5.215 & x_0 = 0.2 & x_1 = 3 \\ \mbox{for } C > 5.215 & x_0 = 0.326 C - 1.5 & x_1 = 3 \end{array} \right . \end{array}\end{split}\]

and for gaseous media

\[\begin{split}\begin{array}{rlll} \mbox{for} & C < 10. & x_0 = 1.6 & x_1 = 4 \\ \mbox{for} & C \in [10.0,\ 10.5[ & x_0 = 1.7 & x_1 = 4 \\ \mbox{for} & C \in [10.5,\ 11.0[ & x_0 = 1.8 & x_1 = 4 \\ \mbox{for} & C \in [11.0,\ 11.5[ & x_0 = 1.9 & x_1 = 4 \\ \mbox{for} & C \in [11.5,\ 12.25[ & x_0 = 2. & x_1 = 4 \\ \mbox{for} & C \in [12.25,\ 13.804[ & x_0 = 2. & x_1 = 5 \\ \mbox{for} & C \geq 13.804 & x_0 = 0.326 C -2.5 & x_1 = 5 . \end{array}\end{split}\]

High Order Corrections

High order corrections term to Bethe-Bloch formula (169) can be expressed as

(174)\[F = G - S + 2(z L_1 + z^2 L_2),\]

where G is the Mott correction term, S is the finite size correction term, \(L_1\) is the Barkas correction, \(L_2\) is the Bloch correction. The Mott term [Ahl80] describes the close-collision corrections tend to become more important at large velocities and higher charge of projectile. The Fermi result is used:

\[G = \pi\alpha z\beta.\]

The Barkas correction term describes distant collisions. The parameterization is expressed in the form:

\[L_1 = \frac{1.29 F_A(b/x^{1/2})}{Z^{1/2}x^{3/2}}, \;\; x = \frac{\beta^2}{Z\alpha^2},\]

where \(F_A\) is tabulated function [ARB73], b is scaled minimum impact parameter shown in Table 22 [eal93]. This and other corrections depending on atomic properties are assumed to be additive for mixtures and compounds.

Table 22 Scaled minimum impact parameter b.
\(Z\)	1 (\(H_2\) gas)	1	2	3 - 10	11 - 17	18	19 - 25	26 - 50	> 50
\(d\)	0.6	1.8	0.6	1.8	1.4	1.8	1.4	1.35	1.3

For the Bloch correction term the classical expression [eal93] is following:

\[z^2L_2 = -y^2 \sum^{\infty}_{n=1} \frac{1}{n(n^2 + y^2)}, \;\; y = \frac{z\alpha}{\beta}.\]

The finite size correction term takes into account the space distribution of charge of the projectile particle. For muon it is zero, for hadrons this term become visible at energies above few hundred GeV and the following parameterization [Ahl80] is used:

\[S = \ln(1 + q), \;\; q = \frac{2 m_e T_{max}}{ \varepsilon^2},\]

where \(T_{max}\) is given in relation (168), \(\varepsilon\) is proportional to the inverse effective radius of the projectile (Table 24).

Table 24 The values of the \(\varepsilon\) parameter for different particle types.
mesons, spin = 0 (\(\pi^{\pm}\), \(K^{\pm}\))	0.736 GeV
baryons, spin = 1/2	0.843 GeV
ions	0.843 \(A^{1/3}\) GeV

All these terms break scaling relation (33) if the projectile particle charge differs from \(\pm\)1. To take this circumstance into account in G4ionIonisation process at initialisation time the term \(F\) is ignored for the computation of the \(dE/dx\) table. At run time this term is taken into account by adding to the mean energy loss a value

\[\Delta T' = 2 \pi r_e^2 mc^2 n_{el} \frac{z^2}{\beta^2} F\Delta s,\]

where \(\Delta s\) is the true step length and \(F\) is the high order correction term (174).

Parameterizations at Low Energies

For scaled energies below \(T_{lim}\) = 2 MeV shell correction becomes very large and precision of the Bethe-Bloch formula degrades, so parameterisation of evaluated data for stopping powers at low energies is required. These parameterisations for all atoms is available from ICRU’49 report [eal93]. The proton parametrisation is used in G4BraggModel, which is included by default in the process G4hIonisation. The alpha particle parameterisation is used in the G4BraggIonModel, which is included by default in the process G4ionIonisation. To provide a smooth transition between low-energy and high-energy models the modified energy loss expression is used for high energy

\[S(T) = S_H (T) + (S_L(T_{lim}) - S_H(T_{lim}))\frac{T_{lim}}{T}, \;\; T > T_{lim},\]

where \(S\) is smoothed stopping power, \(S_H\) is stopping power from formula (169) and \(S_L\) is the low-energy parameterisation.

The precision of Bethe-Bloch formula for \(T\)>10 MeV is within 2%, below the precision degrades and at 1 keV only 20% may be guaranteed. In the energy interval 1–10 MeV the quality of description of the stopping power varied from atom to atom. To provide more stable and precise parameterisation the data from the NIST databases are included inside the standard package. These data are provided for 316 predefined materials (98 elemental and 180 compounds). Note that 278 are “real” NIST materials taken from [NISa][NISb][SBS84] and the remainder are based on their chemical formulae (16 HEP Materials, 3 Space Science Materials and 19 Biomedical Materials). The data from the PSTAR database are included into G4BraggModel. The data from the ASTAR database are included into G4BraggIonModel. So, if Geant4 material is defined as a NIST material, than NIST data are used for low-energy parameterisation of stopping power. If material is not from the NIST database, then the ICRU’49 parameterisation is used. It is suggested to refer to the class G4NistMaterialBuilder to determine the correct nomenclature and composition for each material.

Nuclear Stopping

Nuclear stopping due to elastic ion-ion scattering since Geant4 v9.3 can be simulated with the continuous process G4NuclearStopping. By default this correction is active and the ICRU’49 parameterisation [eal93] is used, which is implemented in the model class G4ICRU49NuclearStoppingModel.

Total Cross Section per Atom

For \(T \gg I \) the differential cross section can be written as

(175)\[\frac{d\sigma }{dT} = 2\pi r_e^2 mc^2 Z \frac{z_p^2}{\beta^2} \frac{1}{T^2} \left[ 1 - \beta^2 \frac{T}{T_{max}} + s\frac{T^2}{2E^2} \right]\]

[ea06], where s = 0 for spinless particles and s = 1 for others. The correction for spin 1/2 is exact and it is not for other values of spin. In described models there is an internal limit \(T_{cut} \geq I\). Integrating from \(T_{cut}\) to \(T_{max}\) gives the total cross section per atom :

(176)\[\sigma (Z,E,T_{cut}) = \frac {2\pi r_e^2 Z z_p^2}{\beta^2} mc^2 \times \left[ \left( \frac{1}{T_{cut}} - \frac{1}{T_{max}} \right) - \frac{\beta^2}{T_{max}} \ln \frac{T_{max}}{T_{cut}} + s\frac{T_{max} - T_{cut}}{2E^2} \right]\]

In a given material the mean free path is:

\[\begin{array}{lll} \lambda = (n_{at} \cdot \sigma)^{-1} & \mbox{or} & \lambda = \left( \sum_i n_{ati} \cdot \sigma_i \right)^{-1} \end{array}\]

The mean free path is tabulated during initialization as a function of the material and of the energy for all kinds of charged particles.

Simulating Delta-ray Production

A short overview of the sampling method is given in Monte Carlo Methods. Apart from the normalization, the cross section (175) can be factorized:

\[\frac{d\sigma}{dT}=C f(T) g(T) \;\; \mbox{with} \;\; T \in [T_{cut}, \ T_{max}]\]

where

\[\begin{split}f(T) &= \frac{1}{T^2} \\ g(T) &= 1 - \beta^2 \frac{T}{T_{max}} + s\frac{T^2}{2E^2} .\end{split}\]

The energy \(T\) is chosen by

sampling \(T\) from \(f(T)\)
calculating the rejection function \(g(T)\) and accepting the sampled \(T\) with a probability of \(g(T)\).

After the successful sampling of the energy, the direction of the scattered electron is generated with respect to the direction of the incident particle. The azimuthal angle \(\phi\) is generated isotropically. The polar angle \(\theta\) is calculated from energy-momentum conservation. This information is used to calculate the energy and momentum of both scattered particles and to transform them into the global coordinate system.

Ion Effective Charge

As ions penetrate matter they exchange electrons with the medium. In the implementation of G4ionIonisation the effective charge approach is used [ZBL85]. A state of equilibrium between the ion and the medium is assumed, so that the ion’s effective charge can be calculated as a function of its kinetic energy in a given material. Before and after each step the dynamic charge of the ion is recalculated and saved in G4DynamicParticle, where it can be used not only for energy loss calculations but also for the sampling of transportation in an electromagnetic field.

The ion effective charge is expressed via the ion charge \(z_i\) and the fractional effective charge of ion \(\gamma_i\):

(177)\[z_{eff} = \gamma_i z_i.\]

For helium ions fractional effective charge is parameterized for all elements

(178)\[\begin{split}(\gamma_{He})^2 & = \left (1-\exp\left [-\sum_{j=0}^5{C_jQ^j}\right ]\right) \left ( 1 + \frac{ 7 + 0.05 Z }{1000} \exp( -(7.6-Q)^2 ) \right )^2, \\ Q & = \max ( 0, \ln T),\end{split}\]

where the coefficients \(C_j\) are the same for all elements, and the helium ion kinetic energy \(T\) is in keV/amu.

The following expression is used for heavy ions [BK82]:

(179)\[\gamma_i = \left ( q + \frac{1-q}{2} \left (\frac{v_0}{v_F} \right )^2 \ln {\left ( 1 + \Lambda^2 \right )} \right ) \left ( 1 + \frac{(0.18+0.0015Z)\exp(-(7.6-Q)^2)}{Z_i^2} \right ),\]

where \(q\) is the fractional average charge of the ion, \(v_0\) is the Bohr velocity, \(v_F\) is the Fermi velocity of the electrons in the target medium, and \(\Lambda\) is the term taking into account the screening effect:

(180)\[\Lambda = 10 \frac{v_F}{v_0} \frac{(1-q)^{2/3}}{Z_i^{1/3}(6+q)}.\]

The Fermi velocity of the medium is of the same order as the Bohr velocity, and its exact value depends on the detailed electronic structure of the medium. The expression for the fractional average charge of the ion is the following:

(181)\[q = [1 -\exp(0.803y^{0.3}-1.3167y^{0.6}-0.38157y-0.008983y^2)],\]

where \(y\) is a parameter that depends on the ion velocity \(v_i\)

(182)\[y = \frac{v_i}{v_0Z^{2/3}} \left ( 1 +\frac {v_F^2}{5v_i^2} \right ).\]

The parametrisation of the effective charge of the ion applied if the kinetic energy is below limit value

(183)\[T < 10 z_i \frac{M_i}{M_p}\;\mbox{MeV},\]

where \(M_i\) is the ion mass and \(M_p\) is the proton mass.

Bibliography

Ahl80(1,2,3): Steven P. Ahlen. Theoretical and experimental aspects of the energy loss of relativistic heavily ionizing particles. Reviews of Modern Physics, 52(1):121–173, jan 1980. URL: https://doi.org/10.1103/RevModPhys.52.121, doi:10.1103/revmodphys.52.121.
ARB73(1,2): J. C. Ashley, R. H. Ritchie, and Werner Brandt. Z13-dependent stopping power and range contributions. Physical Review A, 8(5):2402–2408, nov 1973. URL: https://doi.org/10.1103/PhysRevA.8.2402, doi:10.1103/physreva.8.2402.
Bic92: Hans Bichsel. Stopping power and ranges of fast ions in heavy elements. Physical Review A, 46(9):5761–5773, nov 1992. URL: https://doi.org/10.1103/PhysRevA.46.5761, doi:10.1103/physreva.46.5761.
BK82: W. Brandt and M. Kitagawa. Phys. Rev. B, 25:5631, 1982.
ea06(1,2,3): W-M Yao et al. Review of particle physics. Journal of Physics G: Nuclear and Particle Physics, 33(1):1–1232, jul 2006. URL: https://doi.org/10.1088/0954-3899/33/1/001, doi:10.1088/0954-3899/33/1/001.
eal84: M.J. Berger et al. Icru report 37. Journal of the International Commission on Radiation Units and Measurements, os19(2):, dec 1984. URL: https://doi.org/10.1093/jicru/os19.2.Report37, doi:10.1093/jicru/os19.2.report37.
eal93(1,2,3,4,5,6,7,8,9): M.J. Berger et al. Report 49. Journal of the International Commission on Radiation Units and Measurements, os25(2):NP–NP, may 1993. ICRU Report 49. URL: https://doi.org/10.1093/jicru/os25.2.Report49, doi:10.1093/jicru/os25.2.report49.
Kha68(1,2): Govind S. Khandelwal. Shell corrections for k- and l-electrons. Nuclear Physics A, 116(1):97–111, jul 1968. URL: https://doi.org/10.1016/0375-9474(68)90485-5, doi:10.1016/0375-9474(68)90485-5.
NISa: NIST. Nist database of atomic weight and isotopic composition. https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses. [Online; accessed 11-december-2017].
NISb: NIST. Nist material composition database. https://physics.nist.gov/cgi-bin/Star/compos.pl. [Online; accessed 11-december-2017].
SP71: R.M. Sternheimer and R.F. Peierls. General expression for the density effect for the ionization loss of charged particles. Physical Review B, 3(11):3681–3692, jun 1971. URL: https://doi.org/10.1103/PhysRevB.3.3681, doi:10.1103/physrevb.3.3681.
SBS84: RM Sternheimer, MJ Berger, and Stephen M Seltzer. Density effect for the ionization loss of charged particles in various substances. Atomic Data and Nuclear Data Tables, 30(2):261–271, 1984.
Wal52: M. C. Walske. The stopping power ofK-electrons. Physical Review, 88(6):1283–1289, dec 1952. URL: https://doi.org/10.1103/PhysRev.88.1283, doi:10.1103/physrev.88.1283.
Wal56: M. C. Walske. Stopping power ofL-electrons. Physical Review, 101(3):940–944, feb 1956. URL: https://doi.org/10.1103/PhysRev.101.940, doi:10.1103/physrev.101.940.
ZBL85: J.F. Ziegler, J.P. Biersack, and U. Littmark. The Stopping and Ranges of Ions in Solids. Pergamon Press, vol.1 edition, 1985.