Particle Transport - Scott Crowe

The probabilities associated with radiation particle interactions in a medium can be expressed as differential cross sections which are probability density functions. Monte Carlo methods can be applied to simulate a chain of interactions experienced by a particle in a medium – this is ‘particle transport’. Described below are examples of how sampling techniques are applied to calculate free path lengths and the results of some interactions important in the consideration of radiotherapy dosimetry. The three major photon interactions responsible for dose deposition: the photoelectric effect, Compton scattering and pair production; are examined and the simulation of electrons summarised.

Rayleigh scattering, though common at the energy values typical of radiotherapy treatments, does not significantly affect the distribution of dose to the patient, since no energy is deposited in the tissue. It has been suggested that Rayleigh scattering need only be modelled when dealing with low energy (<1 MeV) simulations (Rogers et al, 2009). Flourescent photons (or associated Auger electrons) produced from atomic relaxation events do not affect the distribution of dose to the patient due to their low energies.

Differential cross sections are calculated as a function of atomic number and incident particle energy: when dealing with compound materials, there are two approaches available to determine what cross section (σ) data to sample. The first is to randomly select a Z according to the atomic composition (Agostinelli et al, 2003):

$Prob(Z_i, E) = \frac{\eta_{ati}\ \sigma(Z_i, E)}{\sum_i [\eta_{ati}\cdot \sigma_i(E)]}$

where η_ati is the number of atoms per unit volume of the i^th element of the material (Agostinelli et al, 2003). Alternatively, the reference differential cross sections can be combined, and a fit parameterised for the effective-Z.

Particles continue to experience interactions until they reach a threshold energy level where they are terminated (chosen to avoid the simulation of interactions insignificant to the dose distribution). The selection of this threshold energy is most important for electrons, which are otherwise captured only when they reach binding energies (usually below 1 keV).

Mean Free Path Length

If a particle is generated at a point in a three dimensional space, with a prescribed motion and energy, the distance the particle travels before an interaction is simulated is called the mean free path length, or step. The following probability density function (Morin, 1988) describes the probability that a collision will occur between at distance l for a particle travelling through a homogeneous medium:

$f(x) = \mu_T e^{- \mu_T \ell}$

where μ_T is the total linear attenuation coefficient (which is also the inverse of the mean free path) of the medium. The CDF is the integral of this

$F(x) = 1 - e^{- \mu_T \ell}$

The selection of a random number in the range [0,1] for F(x) can be mapped to an x value (Metcalfe et al, 2007): this is inverse sampling; the inverse of this function can be expressed as (Raeside, 1976)

$\ell = - \frac{1}{\mu_T} ln(r)$

where l is the sampled path length, the distance the particle travels before an interaction occurs, and r is a random number in the interval (0,1].

The total linear attenuation coefficient, μ_T, represents the probability per unit path length for some interaction (Midgley, 2004), and is proportional to the total cross section σ_T:

$\mu_T = \frac{\rho \sigma_T}{u A}$

where ρ is the density of the medium, u is the atomic mass unit and A is the atomic mass of the element. The total cross section is the sum of the cross-sections for each individual process. For a photon, this could be

$\sigma_T = \sigma_{\mathrm{rayleigh}} + \sigma_{\mathrm{pe}} + \sigma_{\mathrm{compton}} + \sigma_{\mathrm{pair}} + \sigma_{\mathrm{triplet}}$

where σ_rayleigh, σ_pe, σ_compton, σ_pair and σ_triplet are the cross-sections for Rayleigh scattering, the photo-electric effect, Compton scattering, pair and triplet production respectively. Once the free path length is determined, the type of interaction (i.e. compton scattering or pair production) can be selected using the composition method.

Photoelectric Effect

Photoelectrons are atomic electrons emitted when energy is absorbed from an incident photon. Cross sections for the photoelectric effect are stored as a function of atomic number, and may be parameterized as (Agostinelli et al, 2003):

$\sigma_{\mathrm{pe}} \left(Z,E\right)= \frac{a\left(Z, E_\gamma\right)}{E_\gamma} + \frac{b\left(Z, E_\gamma\right)}{E_\gamma^2} + \frac{c\left(Z, E_\gamma\right)}{E_\gamma^3} + \frac{d\left(Z, E_\gamma\right)}{E_\gamma^4}$

where Z is the atomic number, E_γ is the energy of the incident photon and the functions a, b, c and d are fit to the reference data (for instance, of Storm and Israel (1970)) in several energy intervals, separated by photoabsorption edges (Agostinelli et al, 2003). Providing the incident photon energy exceeds the binding energy of the atomic shell, a photoelectron is emitted with the energy:

$E_{e^-} = E_\gamma - E_{\mathrm{shell}}$

where E_shell is the binding energy (for example, the k-edge energy). The properties of the emitted photoelectron must also be selected, for example, the polar angle (θ) could be determined by sampling the Sauter-Gavrila distribution (Kawrakow and Rogers, 2003):

$\cos \theta = \frac{1}{\beta} \left[1 - {\left( \sqrt{{\left(\kappa + \frac{1}{1 + \beta} \right)}^2 + 4 \beta \gamma^2 \left( \kappa + \gamma^2 \right) r } - \kappa \right)}^{-1} \right]$

where β is the electrons velocity in terms of c, the speed of light, and r is a uniformly distributed random number,

$\kappa = \frac{\gamma}{2} \left( \gamma - 1\right) \left(\gamma - 2\right) \mathrm{\ and\ } \gamma = \frac{1}{\sqrt{1 - \beta^2}}.$

Where the incident photon energy does not exceed the binding energy for the electron, the energy would be deposited at that point. If relaxation events are not modelled (they aren’t required for accurate radiotherapy simulations) any energy difference can also be deposited at that point.

Compton Scattering

Compton scattering occurs when a photon imparts energy to an orbital electron. This is the principal mechanism for the deposition of dose by a radiotherapy beam. The cross-section for the compton interaction can be parameterised as (Agostinelli et al, 2003)

$\sigma_{\mathrm{compton}}(Z, E_\gamma) = \left[ P_1(Z) \frac{\log(1 + 2X)}{X} + \frac{P_2(Z) + P_3(Z)X + P_4(Z)X^2}{1 + aX + bX^2 + cX^3} \right]$

where Z is the atomic number, E_γ is the energy of the incident photon, X = ln(E_γ / m_e c²) and P_i(Z) = Z(d_i + e_iZ + f_iZ²). The remaining coefficients (a to f) are obtained by fitting the function to reference data. The resultant polar angle (θ) of a compton interaction can be simulated knowing the Klein-Nishina cross-section (Klein and Nishina, 1929; Kawrakow and Rogers, 2003):

$\frac{d\sigma}{d \cos \theta} = \pi r_e^2 Z {\left( \frac{E_\gamma^\prime}{E_\gamma} \right)}^2 \left[ \frac{E_\gamma^\prime}{E_\gamma} + \frac{E_\gamma}{E_\gamma^\prime} - \sin^2 \theta \right]$

where r_e is the classical radius of an electron and E’_γ is the resultant photon energy, which can be related to the scattering angle (Agostinelli et al, 2003):

$E_\gamma^\prime = E_\gamma \frac{m_e c^2}{m_e c^2 + E_\gamma (1 - \cos \theta)}$

Rejection sampling is used to obtain a θ value using the probability density function P₁ (cos θ) (Kawrakow and Rogers, 2003):

$P_1 (\cos \theta) = \frac{\frac{d\sigma}{d \cos \theta} \ d \cos \theta}{\int^1_{-1} \frac{d\sigma}{d \cos \theta} \ d \cos \theta} = \frac{{\left( \frac{E_\gamma^\prime}{E_\gamma} \right)}^2 \left[ \frac{E_\gamma^\prime}{E_\gamma} + \frac{E_\gamma}{E_\gamma^\prime} - \sin^2 \theta \right] d\cos\theta}{\int^1_{-1} {\left( \frac{E_\gamma^\prime}{E_\gamma} \right)}^2 \left[ \frac{E_\gamma^\prime}{E_\gamma} + \frac{E_\gamma}{E_\gamma^\prime} - \sin^2 \theta \right] d\cos\theta}$

where the Klein-Nishina cross section has been integrated between the maximum and minimum values for cos θ: -1 and 1. This can be re-written using ε = E’_γ / E_γ (Kawrakow and Rogers, 2003)

$P_1 (\epsilon) = P_1 (\cos\theta) \frac{d\epsilon}{d\cos\theta} = N \left( \frac{1}{\epsilon} + \epsilon - \sin^2\theta \right)$

where N is a normalisation constant. The maximum and minimum ε values can be defined as

$\epsilon_{\mathrm{min}} = \frac{1}{1+2 E_\gamma}, \ \epsilon_{\mathrm{max}} = 1$

Sampling the ratio ε is effectively sampling cos θ and can be done using the rejection method (Kawrakow and Rogers, 2003). First ε_min is calculated knowing the incident E_γ, then the upper rejection boundary g_max is calculated using:

$g_{\mathrm{max}} = \frac{1}{\epsilon_\mathrm{min}} + \epsilon_{\mathrm{min}}$

Two random numbers, r₁ and r₂, are generated and an ε value selected (between ε_min and 1) using

$\epsilon = r_1 + (1 - r_1) \epsilon_\mathrm{min}$

which is substituted into a rejection function g:

$g = \frac{1}{g_\mathrm{max}} \left( \frac{1}{\epsilon} + \epsilon - \sin^2\theta \right)$

where sin²θ is determined using ε. If r₂ ≤ g then ε has been sampled (that is, not rejected), and with it the corresponding polar angle θ.

The recoil electron would receive the difference in energies (E_γ – E’_γ) and is deflected by an angle ψ such that momentum is conserved. It should be noted that binding effects and Doppler broadening are not simulated using the Klein-Nishina formula, and while it would be possible to augment it with an impulse approximation (Rogers et al, 2009), this would be unnecessary for accurate radiotherapy simulations. This is because only a small number of recoil electrons are produced at low energies, and the binding energy can be ignored when it is small compared to the recoil electron energy (Agostinelli et al, 2003).

Pair Production

Pair production occurs when a photon (of an energy greater than or equal to 1.022 MeV) interacts with a nucleus, and results in the production of an electron-positron pair. This interaction is most likely for higher energy simulations. The cross-section can be expressed as (Agostinelli et al, 2003):

$\sigma_{\mathrm{pair}}(Z, E_\gamma) = Z(Z+1) \left[ F_1(X) + F_2(X) Z + \frac{F_3(X)}{Z}\right]$

where E_γ is the incident photon energy, X = ln(E_γ / m_e c²) (where m_e is the electron rest mass) and the functions F_i parameterised to best fit the experimental data for the atomic number Z. The properties of the resultant pair are governed by the law of conservation of momentum: they receive equal momentum and the polar angles must be equal in magnitude.

Triplet production occurs when an atomic electron is ejected during a pair production event (Nelson et al, 1985). Pair and triplet production contributions are usually combined as general pair production (Podgorak, 2005) and may be simulated as regular pair production events (Nelson et al, 1985). This is because the ratio of triplet production to pair production cross sections is of the order 1/Z (Nelson et al, 1985), and is, at most (in the case of hydrogen), only approximately 30% (Podgorak, 2005).

Electron Simulation

Electrons passing through a medium constantly undergo scattering interactions, of which many are ‘soft’, having little effect on the dose distribution. These include elastic collisions (such as Rutherford Scattering) which are zero-energy loss events (Kawrakow and Rogers, 2003) and inelastic collisions where the energy loss is below a defined threshold. The direct simulation of all of these interactions would be inefficient (Foote and Smyth, 1995) so in the interest of minimising simulation time, ‘soft’ interactions are often ‘condensed’ using multiple scattering theories (Agostinelli et al, 2003; Kawrakow and Rogers, 2003) Condensed history approaches approximate the energy loss and displacement of an electron experiencing multiple interactions over a step (Foote and Smyth, 1995). The step length is typically dependent on the presence of voxel boundaries in the simulation geometry (such that energy loss can be deposited into single dose bins) and user-specified energy loss limits (set to regulate step length). These steps occur between discrete non-soft, or ‘catastrophic’, interactions.

‘Catastrophic’ interactions: Bremsstrahlung production, ionisation (via Moller or Bhabha scattering) and annihilation; are not simulated using multiple scattering theories. The same basic approach used for photon interactions is applied here: for a given energy and atomic number, a general cross-section model is fitted to reference data. The subsequent energy and direction of the incident and produced particles is sampled using probability density functions.

Bremsstrahlung radiation is produced when the velocity of an electron is altered by an atomic nucleus: the difference in energy before and after this deceleration is emitted as a photon. This process is arguably the most computationally intensive: the cross-section equation is complex, the energy state of the emitted photons can require 36 linear parameters for parameterization (Agostinelli et al, 2003) and the angular distribution of the emitted photons requires the evaluation of natural log functions (Kawrakow and Rogers, 2003).

There are two ionisation (or `catastrophic’ inelastic collision) interactions: Moller (electron-electron) and Bhabha (electron-positron) scattering. These can be used when binding effects are ignored (Kawrakow and Rogers, 2003), which is reasonable at high energies. Once the cross section has been sampled to determine the kinetic energy of the initially-at-rest electron the polar scattering angles can be calculated using that value and the kinetic energy of the incident particle.

Annihilation occurs when a positron (produced by a pair production event) interacts with an atomic electron, producing two photons (with properties satisfying the laws of conservation of momentum and energy). Almost all positrons produced in the patient will be terminated in this way. Available cross section data, such as Heitler’s (1954), assumes the production of two 511 keV photons, with equal and opposite momentum.

References

Agostinelli S et al, 2003. GEANT4 – a simulation toolkit, Nuclear Instruments and Methods in Physics Research A, 506(3): 250-303.
Foote B J and Smyth V G, 1995. The modelling of electron multiple scattering in EGS4/PRESTA and its effect on ionisation chamber response, Nuclear Instruments and Methods in Physics Research B, 100(1): 22-30.
Heitler W, 1954. The Quantum Theory of Radiation.
Kawrakow I and Rogers D W O, 2003. The EGSnrc Code System.
Klein O and Nishina Y, 1929. ber die Streuung von Strahlung durch freie Elektronen nach der neuen relativistischen Quantendynamik von Dirac, Zeitshcrift fur Physik, 52: 853.
Podgorak E B, 2005. Radiation Physics for Medical Physicists.
Raeside D E, 1976. Monte Carlo principles and applications, Physics in Medicine and Biology, 21: 181-197.
Rogers D W O, Walters B and Kawrakow I, 2009. BEAMnrc Users Manual.
Storm E and Israel H, 1970. Photon cross sections from 1 keV to 100 MeV for elements Z = 1 to 100, Nuclear Data Tables, A7: 565-681.
Metcalfe P, Kron T and Hoban P, 2007. The Physics of Radiotherapy X-Rays and Electrons.
Midgley S M, 2004. A parameterization scheme for the x-ray linear attenuation coefficient and energy absorption coefficient, Physics in Medicine and Biology, 49: 307-325.
Morin R L, 1988. Monte Carlo Simulation in the Radiological Sciences.
Nelson W R, Hirayama H and Rogers D W O, 1985. The EGS4 Code System.