Motivated by the specific needs of protontherapy and proton-induced x-ray emission, we have recently extended the penelope code system [1] to introduce the simulation of interactions and transport of protons [2, 3]. The resulting code, named penh, performs class-II Monte Carlo simulation of coupled electron-photon-proton transport in material structures consisting of homogeneous bodies limited by quadric surfaces. In class II simulation schemes, hard interactions (that is, interactions involving scattering angles or energy losses larger than preselected cutoffs) are simulated by random sampling from their restricted differential cross sections (DCSs) while the cumulative effect of soft (sub-cutoff) interactions is described by means of a multiple-scattering approach. Class II simulation has distinct advantages in front of conventional condensed simulation (see, e.g., Ref. [1]).

The penh code accounts for elastic and inelastic collisions of protons with atoms, as well as for nuclear reactions induced by proton impact. Elastic collisions of protons with nuclei are described by means of numerical DCSs obtained as the product of the DCS for scattering by the bare nucleus, which was computed by the method of partial-waves with the global optical-model potential of Koning and Delaroche [4], and a screening factor accounting for the effect of the screening of the Coulomb field of the nucleus by the atomic electrons, which was calculated from the eikonal approximation for the Dirac-Hartree-Fock-Slater potential of neutral atoms [2, 3, 5]. Proton-induced electronic excitations are simulated by means of DCSs obtained from the Born approximation with the Sternheimer-Liljequist model of the generalized oscillator strength (GOS) [1], which is modified by rescaling the relative contributions of the various electron subshells to give cross sections for the impact ionization of inner subshells equal to unpublished results obtained from accurate atomic GOSs [2]. Nuclear reactions induced by proton impact are simulated by using information from data files in ENDF-6 format [6], which provide cross sections, multiplicities, and angle-energy distributions of all reaction products: light ejectiles (up to alphas), gammas, as well as recoiling heavy residuals.

Although electromagnetic interactions are faithfully described for protons in the energy range from 100 keV up to 1 GeV, the unavailability of nuclear reaction data for energies higher than about 200 MeV, limit the possible applications of the code. Simulations with penh are also limited by the fact that the code can only follow electrons/positrons, photons, and protons. To lessen the effect of this limitation, light charged ejectiles other than protons (deuterons, tritons, ³He, and alphas) are tracked as weighted protons [3]. Neutrons and heavy products are not followed and their kinetic energies were assumed to be deposited at the reaction site. Comparison with results from other codes shows that, as one could predict, the calculated dose distributions from proton beams are generally accurate in the volume swept by the protons, but are generally too low outside that volume, mostly because of the neglect of neutron transport. In this respect, it is worth recalling that neutrons practically do not interact with electrons and, as a consequence, they have large mean free paths. In the energy range covered by the proton simulation code, nominally from 100 keV to 1 GeV, the mean free path of neutrons in water ranges from about 5 cm at 100 keV to ∼ 66 cm at 1 GeV.

The aim of the present article is twofold. First, we present the theory and computational aspects of the calculation of the DCSs for elastic collisions of neutrons and protons, and of the cross sections for ionization of inner electron subshells of atoms by proton impact. We also describe the structure and contents of the associated numerical databases, which completely define the interaction models used in penh for nucleons. Second, we present a simplified algorithm for the simulation of the effect of neutron transport in Monte Carlo calculations of dose distributions from proton beams, which only uses the information provided in the calculated elastic-scattering database, that is, the total and DCS for elastic collisions and the reaction cross section, which is one of the basic parameters used to set the nuclear optical-model potential. The main simplification is that neutron-induced nuclear reactions are described as purely absorptive, that is, when the transported neutron induces a reaction, the simulation of the neutron is discontinued and a fraction of its kinetic energy is assumed to be locally absorbed. It will be shown that, in spite of its crudeness, this transport algorithm accounts for the contribution of neutrons to the spatial dose distribution fairly accurately. Although random histories of neutrons with protontherapy energies (say, from 100 keV to ∼ 300 MeV) may be simulated detailedly (i.e., interaction by interaction) we adopt a class II tracking scheme, which is analogous to the ones used in penelope-penh for electrons/positrons and protons, because it allows speeding up simulations of high-energy neutrons.

The present article is organized as follows. In Section 2 we describe the calculation of DCSs for elastic collisions of nucleons and protons with atoms (nuclei), the structure of the numerical database for elastic scattering of these particles, and the implementation of a class II algorithm for the simulation of multiple elastic scattering of neutrons. Section 3 deals with the ionization of inner electron subshells of atoms by impact of protons, which is based on total cross sections calculated from the relativistic plane-wave Born approximation with corrections for binding/polarization and Coulomb deflection effects. After a brief presentation of the theory and sample results, the associated numerical database is described. The validation of the penh calculation of dose distributions from proton beams is considered in Section 4. Section 5 illustrates the simulation of proton-induced x-ray emission from solid targets. Finally, in Section 6 we offer a few concluding remarks.

We consider elastic collisions of nucleons with atoms of the element of atomic number Z. In order to cover the range of kinetic energies of interest in proton therapy, up to about 300 MeV, we shall use relativistic collision kinematics. The simulation code transports particles in the laboratory (L) frame, where the material is at rest and the projectile moves with kinetic energy E before the collision. For simplicity, we consider that the z axis of the reference frame is parallel to the linear momentum of the projectile, which is given by p = c − 1 E E + 2 m c 2 z ̂ , (1) where c is the speed of light in vacuum and m is the projectile mass (= m _p for protons, and = m _n for neutrons). The total energy of the projectile is W = E + m c 2 = m 2 c 4 + c 2 p 2 . (2) We recall the general relations p = β γ m c and E = γ − 1 m c 2 (3) where β = v c = E E + 2 m c 2 E + m c 2 (4) is the speed of the particle in units of c and γ = 1 1 − β 2 = E + m c 2 m c 2 (5) is the particle’s total energy in units of its rest energy.

Elastic collisions involve a certain transfer of kinetic energy to the target atom, which is easily accounted for by sampling the collisions in the center-of-mass (CM) frame, which moves relative to the L frame with velocity v C M = β C M c = c 2 p E + m c 2 + M A c 2 , (6) where M _A is the mass of the atom. In the CM frame the linear momenta of the projectile and the atom before the collision are, respectively, p i ′ = p 0 ′ and p A i ′ = − p 0 ′ , with p 0 ′ = M A c 2 m c 2 + M A c 2 2 + 2 M A c 2 E p . (7) Quantities in the CM frame are denoted by primes. After the elastic collision, in CM the projectile moves with momentum p f ′ = p 0 ′ in a direction defined by the polar scattering angle θ′ and the azimuthal scattering angle ϕ′, and the target atom recoils with equal momentum p A f ′ = p 0 ′ in the opposite direction. The final energies and directions of the projectile and the atom in the L frame are obtained by means of a Lorentz boost with velocity − v _CM. Thus, elastic collisions are completely determined by the DCS dσ/dΩ′ in the CM frame.

We follow the approach described in our previous work on proton transport [3], i.e., we assume that the interaction potential in the CM frame is central, because this is a prerequisite for applying the partial-wave expansion method to compute the DCS. Our approach can be qualified as semi-relativistic, because we are using strict relativistic kinematics but we do not account for the breaking of the central symmetry of the interaction when passing from the L to the CM frame.

The slowing down of fast protons in matter is mostly due to inelastic collisions, i.e., collisions causing electronic excitations of the material. Protons also slow down due to elastic collisions with nuclei (the so-called nuclear stopping), an effect that is automatically accounted for by the simulation scheme adopted in penh. The plane-wave Born approximation (PWBA) [15, 16] is suited for describing inelastic collisions of protons with velocities much larger than those of the active target electrons. We consider inelastic collisions of a proton with initial kinetic energy E and momentum p [(cp)² = E (E + 2 mc²)] with atoms of atomic number Z, characterized by the energy transfer W = E − E′ and the momentum transfer q = p − p′, where E′ and p′ are, respectively, the energy and momentum of the projectile after the collision. The corresponding DCS is most conveniently considered as a function of the energy transfer W and the recoil energy Q, defined as the kinetic energy of an electron with linear momentum equal to the momentum transfer [15], Q = c p 2 + c p ′ 2 − 2 c p c p ′ cos ⁡ θ + m e 2 c 4 − m e c 2 . (42) The DCS obtained from the PWBA for ionizing collisions with electrons of a subshell a with binding energy E _a can be expressed as [17]. d 2 σ a d W d Q = 2 π z 0 2 e 4 m e v 2 2 m e c 2 W Q Q + 2 m e c 2 d f a Q , W d W + 2 m e c 2 W Q Q + 2 m e c 2 − W 2 2 β 2 − W 2 Q Q + 2 m e c 2 d g a Q , W d W (43) where z ₀ = + 1 is the proton charge in units of e, and the functions df (Q, W)/dW and dg (Q, W)/dW are, respectively, the longitudinal and transverse generalized oscillator strengths (GOS). Bote and Salvat [17] have calculated these GOS for all electron shells of atoms with Z = 1 to 99 by using an independent-electron model with the Dirac-Hartree-Fock-Slater self-consistent potential of free atoms. The energy-loss DCS is obtained by integration over the kinematically allowed interval of recoil energies, (Q ₋, Q ₊), with endpoints given by Eq. 42 with cos θ = + 1 and −1, d σ a d W = ∫ Q − Q + d 2 σ a d W d Q d Q . (44)

In penelope-penh proton inelastic collisions are described by means of the PWBA with a simplified model of the GOS [1, 2], which is modulated so as to reproduce the cross sections for ionization of inner electron subshells read from the database. The main limitation of the PWBA is due to the neglect of the distortion of the projectile wave functions caused by the field of the target atom. For electrons and positrons, this distortion can be largely accounted for by using the distorted-wave Born approximation (DWBA), in which the projectile states are represented as distorted plane waves (see, e.g. [18], and references therein). Because of the slow convergence of the partial-wave series, this kind of calculation is only possible for projectile electrons and positrons with kinetic energies up to about 30E _a . Bote and Salvat [17] used an optimized computation strategy, which combines the DWBA and the PWBA, to generate a database of electron-impact ionization cross sections for the K shell and the L and M subshells of all the elements from hydrogen to einsteinium (Z = 1–99) and for energies of the projectile from 50 eV up to 1 GeV. The results were found to agree well with available experimental data [19].

Unfortunately, the calculation of cross sections for inelastic collisions from the DWBA is not feasible for charged particles heavier than the electron, because the smallness of the de Broglie wavelength of the projectile renders the calculation of free spherical waves extremely difficult. Chen et al. [20], and Chen and Crasemann [21, 22] went beyond the PWBA by using the perturbed-stationary-state approximation of Brandt and Lapicki [23], which accounts for 1) alterations in the binding of the active electron due to the presence of the projectile near the nucleus of the target atom, and 2) the deflection of the projectile path caused by the Coulomb field of the nucleus. In our PWBA calculations of ionization cross sections for protons, these effects are introduced by means of semi-classical correcting factors, which are described in the following paragraphs.

• Binding effect

In collisions where the projectile proton penetrates deep into the target atom, the presence of the projectile modifies the binding energy of the active electron and, in the case of positively charged projectiles, leads to a reduction of the DCS. For the K shell and L subshells, Brandt and Lapicki [23] performed a first-order perturbation analysis, assuming that the projectile follows a straight trajectory and using hydrogenic wave functions. They obtained an ionization-energy shift of the active target electron given by Δ E a = 2 z 0 E a Z a Θ a g a ξ − h a ξ , (45) where Z _a = Z − δ _a is the effective nuclear charge felt by the electrons in the unperturbed orbitals, with δ _K = 0.3 and δ _Li = 4.15. The quantity Θ _a is the reduced ionization energy, Θ a = 2 n a 2 E a / Z a 2 E h , (46) where E _h = m _e e ⁴/ℏ ² = 27.211 eV is the Hartree energy. The last factor in Eq. 45 is a function of the dimensionless parameter ξ ≡ Z a E h n a E a m e M 2 E E h . (47) In the calculations of ionization of K and L electrons we used the parameterization of the g _a (ξ) and h _a (ξ) functions given by Chen and Crasemann [21].

For M and outer shells, Chen et al. [20] considered that the effective ionization energy is the one of the “united” atom (i.e., of the atom with atomic number z ₀ + Z). We have adopted a similar approach, which avoids the need of considering ionization energies of other atomic species. Expressing the ionization energies of the unperturbed states as (screened hydrogenic levels) E a = Z − δ a ′ 2 2 n a 2 E h , (48a) and noting that the screening constant δ a ′ is nearly the same for neighbouring elements, we can approximate the effective ionization energy in the form E a ′ = z 0 + Z − δ a ′ 2 2 n a 2 E h . (48b)

This gives the following ionization-energy shift Δ E a = E a ′ − E a = z 0 2 + 2 z 0 Z − δ a ′ 2 n a 2 E h . (49)

• Coulomb deflection

For projectiles with small speeds, the PWBA, and the equivalent straight-trajectory semi-classical approximation [see, e.g., Ref. [24], and references therein], overestimate the ionization cross sections because they neglect the effect of the Coulomb field of the nucleus on the trajectory of the projectile. In the semi-classical treatment, the energy-loss DCS for a projectile following a classical hyperbolic orbit in the Coulomb potential of the bare target nucleus can be obtained by multiplying the energy-loss DCS, calculated by assuming that the projectile follows a straight trajectory, by a correction factor. This Coulomb-deflection factor was approximated as [23]. F a C o u l E ; W = 1 − 1 3 x 1 / 3 + 5 3 x 2 / 3 exp − 2 π x , (50) where x = Z 0 Z M 3 8 m e M r e d 2 E h E 3 W . (51)

The ionization cross section, including the binding and Coulomb-deflection corrections, is given by σ a i o n = ∫ E a + Δ E a W max F a C o u l E ; W d σ a c o n t d W d W . (52)

We have generated a database of proton-impact ionization cross sections for the K shell and the L, M and N subshells of the elements with Z = 1 to 99, by numerical integration of the energy-loss DCS obtained from the longitudinal and transverse GOSs calculated by Bote and Salvat [17]. Figure 1 displays these ionization cross sections for electron subshells of the titanium and gold atoms. For comparison purposes, the plots include also values calculated by Chen and Crasemann [21, 22] for the DHFS potential but only with the longitudinal GOS. It is worth mentioning that our database includes relativistic effects in a more consistent way, and it covers the energy range up to 10¹¹ eV for all elements from hydrogen (Z = 1) to einsteinium (Z = 99).

In the simulation code, the ionization cross sections are multiplied by an energy-dependent factor (deduced from the Sternheimer-Liljequist GOS model) that accounts for the reduction of the cross section caused by the density-effect correction [3].

The databases presented in the previous Sections for elastic collisions and impact ionization of protons provide realistic models for electromagnetic interactions of protons, and the class II simulation scheme avoids the use of multiple scattering approximations and their associated uncertainties.

As mentioned above, the original penh code [3] gave a spatial dose distribution that was nearly correct near the proton beam axis, but the simulated dose was too small far from the beam axis. This was first noted by Verbeek et al. [25] and latter confirmed by the authors through comparison with results from the fluka and Geant4 codes, which include neutron transport. The main cause of this underestimation of the distant dose was attributed to the neglect of neutron transport, and this motivated the inclusion of the present simple approach for neutrons. Preliminary simulations with the parameter F N A B S = 1 gave depth-dose distributions from pencil beams of protons that were closer to, but slightly exceeded those from Geant4. The discrepancy can be readily corrected by using a smaller value of F N A B S .

Figure 2 compares depth-dose distributions of proton beams incident normally on a water phantom with energies of 100 and 200 MeV simulated with the penh and Geant4 codes. The red histograms were produced by assuming that neutrons resulting from proton-induced nuclear reactions are absorbed at the reaction site. The blue histograms resulted from assuming that neutrons are transported as described in Section 2.4 with F N A B S = 1 . It is seen that the transport and absorption of neutrons increases the depth-dose beyond the Bragg peak by nearly a factor of 10, slightly exceeding the depth-dose generated by the Geant4 code (green histograms). Reducing the value of F N A B S to 0.85 and to 0.8 for 100 and 200 MeV protons gives depth-dose distributions nearly coincident with those from Geant4. The latter are seen to agree closely with the penh result, except in the deep fall of the Bragg peak, probably because of differences in the proton transport physics of the two codes. Another relevant effect of neutron transport is a slight decrease of the depth-dose at shallow depths, which is visible in the plot for the 200 MeV beam; this effect is expected to increase in importance with the energy of incident protons because of the higher energies of the released neutrons.

As penh uses fairly reliable cross sections for proton impact ionization, together with a careful modeling of electromagnetic interactions, it may be used to simulate proton-induced x-ray emission (PIXE) spectra [26]. The code describes the relaxation of ionized atoms by means of the penelope subroutines, which use transition probabilities from the Evaluated Atomic Data Library of Perkins et al. [27] and empirical values of the x-ray energies [1]. The simulation of PIXE spectra by penh is analogous to that of electron induced x-ray emission by the penelope code, which has been proved to be effective for quantification in electron-probe microanalysis [28]. It is worth noticing that penh accounts for both the attenuation and the generation of x-ray fluorescence within the target. It also follows the bremsstrahlung emitted by secondary electrons, which produces a smooth background of in the simulated PIXE spectra. Primary protons also contribute to the background through the emission of atomic bremsstrahlung (see [29] and references therein), however, this mechanism is not accounted for in penh. In its present form, our simulation code can be useful to determine the influence of the composition and local geometry of the irradiated target on the line intensities of the emitted x rays.

Figure 3 displays results from two simulations of x-ray spectra emitted from a silver target bombarded with 3 MeV protons impinging in a direction at an angle of 45° from the target surface. In those simulations neutrons were not followed, and the variance reduction techniques of interaction forcing and emission splitting of bremsstrahlung photons and x rays [1] were applied. The upper plot shows the energy spectrum of all photons that emerge from the target. The continuous background corresponds to bremsstrahlung emitted by secondary electrons. Since the maximum energy of secondary electrons released in inelastic collisions of 3 MeV photons is about 6.5 keV, the spectral background ends at this energy. A weak and noisy background component that extends to higher energies in the simulated spectrum (not shown in the upper plot of Figure 3) corresponds to gamma rays released in proton-induced nuclear reactions. The lower plot of Figure 3 shows the energy spectrum of photons that emerge in directions forming angles less than 30° with the normal to the target, restricted to the energy interval of L lines. The most prominent x-ray lines can be readily identified from their numerical values in the penelope database.

Simulated photon spectra correspond to an ideal detector with unit efficiency. They are output in the form of histograms with a fixed bin width, which implies a resolution in energy of the order of the bin width. To get results with an appearance closer to measurement data, the simulated spectrum should be convolved with the response function of the detector [28]. Comparison of spectra so obtained with measured spectra may help to identify features not evident from the experimental result.

In its present form, the penh code provides a consistent description of electromagnetic interactions of electrons, positrons and protons with matter for projectiles with energies up to 1 GeV. The use of nuclear data from ENDF-formatted files allows accounting for proton-induced nuclear reactions, and the release of gammas and secondary particles resulting from these interactions, in the energy range covered by available libraries, which usually extends up to about 200 MeV. For protons with somewhat higher energies, the code can extrapolate the nuclear data, with the risk of distorting the results.

Since the interaction models implemented in penh lose validity at low energies, the code should not be used for electrons, positrons and photons with E ≲ 1 keV, and for protons and neutrons with E < 100 keV.

The method adopted for tracking neutrons is intended only to correct for the effect of neutron transport on the simulated dose distributions, of interest mostly in proton therapy. Processes where neutrons may have more relevance cannot be dealt with penh.