Electromagnetic interaction models for Monte Carlo simulation of protons and alpha particles

Electromagnetic interactions of protons and alpha particles are modeled in a form that is suitable for Monte Carlo simulation of the transport of charged particles. The differential cross section (DCS) for elastic collisions with neutral atoms is expressed as the product of the DCS for collisions with the bare nucleus and a correction factor that accounts for the screening of the nuclear charge by the atomic electrons. The screening factor is obtained as the ratio of the DCS for scattering of the projectile by an atom with a point nucleus and the parameterized Dirac-Hartree-Fock-Slater (DHFS) electron density, calculated from the eikonal approximation, and the Rutherford DCS for collisions with the bare point nucleus. Inelastic collisions, which cause electronic excitations of the material, are described by means of the plane-wave Born approximation, with an empirical simple model of the generalized oscillator strength (GOS) that combines several extended oscillators with resonance energies and strengths determined from the atomic configurations and from the empirical mean excitation energy of the material. The contributions from inner subshells are renormalized to agree with realistic ionization cross sections calculated numerically from the DHFS self-consistent model of atoms by means of the plane-wave Born approximation. The resulting DCS allows analytical random sampling of individual hard inelastic interactions.


Introduction
Monte Carlo simulation of the transport of fast charged particles in matter is difficult because of the large number of interactions undergone by these particles before being brought to rest [1; 2].This difficulty can be solved by using two alternative strategies: 1) conventional condensed simulation, or class-I simulation, which consists of splitting each particle trajectory into a number of steps of definite length and making use of multiple scattering theories [3; 4; 5; 6] for describing the cumulative effect of the multiple interactions that occur along each step, and 2) mixed, or class-II simulation, where hard interactions involving energy transfers or angular deflections larger than predefined cutoff values are simulated individually, and soft interactions are described collectively by means of a multiple-scattering approach [7; 8; 9].Class-II schemes are superior because hard interactions are treated exactly by random sampling from the corresponding restricted differential cross sections (DCSs), although they require knowledge of the various DCSs and accurate sampling methods for hard interactions must be implemented in the simulation code.In the present article we describe realistic DCSs for elastic and inelastic electromagnetic interactions of protons and alpha particles with matter, together with algorithms for the restricted random sampling of hard interactions.The proposed simulation strategies are applicable to other charged particles heavier than the electron.
For the sake of generality, the theoretical interaction models are formulated for the general case of projectile particles with mass M 1 , assumed to be larger than the electron mass m e , and charge Z 1 e, where e denotes the elementary charge.The considered interactions are elastic collisions with atoms (i.e., interactions that do not cause excitations of the material) and inelastic interactions, which result in electronic excitations of the medium.These interactions are essentially electromagnetic and can be described quite reliably from first-principles calculations or from appropriate models.
A simulation program transports particles in the laboratory (L) frame, where the material is at rest and the projectile moves with kinetic energy E before the interaction.In order to cover the range of kinetic energies of interest in applications, we shall use relativistic collision kinematics.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 where c is the speed of light in vacuum and M 1 is the projectile rest mass, M 1 = m p = 1836.15m e for protons, m a = 7294.30m e for alphas. ( The rest energy of the electron is m e c 2 = 511.00keV.The total energy of the projectile is We recall the general relations where is the speed of the particle in units of c and is the particle's total energy in units of its rest energy.The present article describes the essential physics involved in the calculation of the DCS and general aspects of the sampling algorithms; details and specific formulas are given in a document available as supplementary material.

Elastic collisions
Let us consider elastic collisions of the projectile with neutral atoms.These collisions involve a certain transfer of kinetic energy to the target atom, which manifests as the recoil of the latter.The recoil of the target atom is easily accounted for by sampling the collisions in the center-of-mass (CM) frame, which moves relative to the L frame with velocity where M A is the mass of the atom.That is, A neutral atom of the element of atomic number Z consists of the atomic nucleus and Z bound electrons in their ground state.The atomic nucleus is a system of Z protons and N neutrons, bound together by the nuclear forces.The total number of nucleons, A ≡ Z + N , is called the mass number.The atomic mass of the isotope A Z is estimated by means a mass formula [10] (see the supplementary document) that approximates the experimental atomic masses of naturally occurring isotopes [11] with a relative accuracy better than about 10 −4 , which is sufficient for the present purposes.
The calculated cross sections for each element are obtained as an average over those of the naturally occurring isotopes, weighted by their respective natural abundances [11].Consistently, in the simulations we consider that the mass of a target atom is the average atomic mass of the element [12] where A w is the molar mass of the element, and u = m( 12 C)/12 is the atomic mass unit.This simplification permits reducing the required information for each element (and projectile kind) to a single cross section table, irrespective of the number of isotopes of that element.
In the CM frame the linear momenta of the projectile and the atom before the collision are, respectively, p ′ i = p ′ 0 and p ′ Ai = −p ′ 0 , with Notice that linear momenta 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 ′ Af = 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 differential cross section (DCS) per unit solid angle, dσ/dΩ, in the CM frame.
We follow the approach described by Salvat and Quesada [13] (see also Ref. [14]), i.e., we assume that the interaction potential in the CM frame is central, since this is a prerequisite for applying the partial-wave expansion method to compute the DCS in the CM frame.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.

Interaction potential
The interaction potential between a charged projectile and the target atom is expressed as where r is the distance between the projectile and the center of mass of the atom, V nuc (r) is the interaction energy of the projectile and the bare atomic nucleus, and Φ(r) is the screening function, which accounts for the shielding of the nuclear charge by the atomic electrons.If the nucleus is represented as a point structureless charged particle, the nuclear potential reduces to the Coulomb potential where Z 1 e the projectile charge (Z 1 = 1 for protons, = 2 for alphas).To facilitate calculations, we use approximate screening functions having the analytical form with the parameters given by [15] for elements with atomic numbers Z = 1 to 92, which were determined by fitting the self-consistent Dirac-Hartree-Fock-Slater (DHFS) atomic potential of neutral free atoms.Parameters for heavy elements with Z = 93 − 99 obtained from the same kind of fit were added more recently.The advantage of using the representation (13) of the screening function is that a good part of the calculation of the DCS for atoms with point nuclei can be performed analytically [14].It is worth noticing that the screened atomic potential vanishes for radial distances r much larger than the "atomic radius", where a 0 = ℏ 2 /(m e e 2 ) = 5.292 × 10 −9 cm is the Bohr radius.
The interaction energy of the projectile with a bare nucleus of the isotope A Z having atomic number Z and mass number A can be described by a phenomenological complex optical-model potential where the first term is a real potential that reduces to the Coulomb potential at large radii, and the second term, iW nuc (r), is an absorptive (negative) imaginary potential which accounts for the loss of projectile particles from the elastic channel caused by inelastic interactions with the target nucleus.Except for the Coulomb tail, the nuclear potential is of finite-range, it vanishes when the distance r from the projectile to the nucleus is larger than about twice the "nuclear radius", Parameterizations of optical-model potentials have been proposed by various authors.In the calculations for protons (and neutrons) we use the parameterization of the nuclear global optical-model potential given by Koning and Delaroche [16], which is valid for projectiles with kinetic energies E between 1 keV and about 200 MeV and nuclei with 24 ≤ A ≤ 209.Owing to the lack of more accurate approximations, because the potential values vary smoothly with A, Z and E, we use those parameters for all isotopes with A > 6 and for energies up to 300 MeV, for higher energies the potential parameters at E = 300 MeV are employed.For protons having E < 35 MeV colliding with target isotopes of mass number A such that 6 < A < 24 (Z < 12), we use the optical-model potential of Watson et al. [17], which is applicable to energies from 10 MeV to 50 MeV; for projectile protons with energies higher than 35 MeV, the potential of Koning and Delaroche is adopted because it yields DCSs in better agreement with available experimental information.For alpha particles, the adopted parameterization of the nuclear potential is the one proposed by Su and Han [18], which is valid for nuclides with 20 ≤ A ≤ 209 and projectiles with kinetic energies up to 386 MeV, although we use it for any nucleus.For alphas with higher energies, we use the parameter values at E = 386 MeV.
In principle, given the interaction potential, the collision DCS can be calculated by the method of partial waves [19].As pointed out by Salvat and Quesada [13], in the energy range of interest for transport calculations, the de Broglie wavelength, λ dB = h/p ′ 0 , of the projectile is much smaller than the atomic radius R at and, consequently, the numerical solution of the radial wave equation to determine the phase-shifts and the DCS is very difficult.In addition, the partial-wave series converge extremely slowly, requiring the calculation of a large number (≳ 100, 000) of phase-shifts.Since approximate calculation methods are available for the case of screened Coulomb potentials (i.e., corresponding to atoms with a point nucleus), we first calculate the DCS for elastic collisions with bare nuclei and introduce the effect of electronic screening as a correction factor to the nuclear DCS.

Elastic collisions with bare nuclei
The scattering of nucleons and alpha particles by nuclei can be described by using the partial-wave expansion method in the CM frame.The underlying physical picture is that of a stationary process represented by a distorted plane wave, i.e., by an exact solution of the time-independent relativistic Schrödinger equation for the potential V nuc (r), with the relativistic reduced mass which asymptotically behaves as a plane wave with an outgoing spherical wave.Owing to the assumed spherical symmetry of the target nucleus, the angular distribution of scattered projectiles is axially symmetric about the direction of incidence, i.e., independent of the azimuthal scattering angle in both the CM and L frames.
In the case of scattering of spin-unpolarized protons (and neutrons), the optical-model potential contains spin-orbit terms, and the wave function is a two-component spinor.The DCS per unit solid angle in CM is [19] where the functions f (θ) and g(θ) are, respectively, the direct and spin-flip scattering amplitudes.They are evaluated from their partial-wave expansions, where P ℓ (cos θ ′ ) and P 1 ℓ (cos θ ′ ) are Legendre polynomials and associated Legendre functions of the first kind [20], respectively, and are the S-matrix elements.The quantities δ ℓa , with a = sign[2(j − ℓ)], are the phase-shifts, which depend on the total and orbital angular momenta of the projectile, j and ℓ, respectively.Inelastic interactions with the nucleus cause a loss of projectile particles from the elastic channel.The reaction cross section, σ react , (i.e., the total cross section for inelastic interactions) is given by The quantities T ℓa = 1−|S ℓa | 2 , the so-called transmission coefficients, measure the fraction of flux that is lost from each partial wave.Since alpha particles have zero spin, the wave function of these particles is a scalar.The DCS for elastic collisions of alpha particles with bare nuclei in the CM frame is given by with the scattering amplitude where [19] S ℓ = exp(2iδ ℓ ).
The reaction cross section for inelastic interactions of alpha particles with the nucleus is The phase shifts δ ℓa and δ ℓ are calculated by using the Fortran subroutine package radial of Salvat and Fernández-Varea [19], which implements a robust power series solution method that effectively avoids truncation errors and yields highly accurate radial functions and phase shifts.The calculations for protons and alpha particles with kinetic energies up to about 1 GeV are doable because their de Broglie wavelengths are comparable to the range of the potential (excluding the Coulomb tail, which determines the kind of "external" radial function), ∼ R nuc .It is worth noticing that global opticalmodel potentials were adjusted to yield reaction cross sections in agreement with measurements and, as a consequence, the calculated values of the reaction cross section and of the DCS are equally reliable.
It is well known that optical-model potentials are not very reliable for light target nuclei.For collisions of protons with light isotopes having A ≤ 6 we use the empirical parameterization of the nuclear DCS described by Galyuzov and Kozov [21], which approximates the available experimental data in an energy range wider than the one needed for transport calculations.For these light isotopes, the reaction cross section is estimated from the empirical formula given by Prael and Chadwick [22].

Electronic screening
Let us consider elastic collisions of the projectile and a target atom of the element of atomic number Z, assuming that the atomic nucleus can be regarded as a point particle.The corresponding interaction potential takes the form of a screened Coulomb potential, where we have introduced the analytical screening function (13).The DCS can then be calculated from the wave equation [14] The DCS for collisions of charged particles with a bare point nucleus, described by the unscreened Coulomb potential V C (r), Eq. ( 12), can be obtained from the exact solution of the wave equation (28) [23] for spinless particles.It is given by the relativistic Rutherford formula, where is the momentum transfer.As indicated above, the smallness of the proton wavelength makes the partial-wave calculation of the DCS for scattering by the screened Coulomb potential unfeasible.A practical approach adopted in Refs.[24; 13] is to use DCSs calculated with the eikonal approximation [25; 26; 27], in which the phase of the scattered wave is obtained from a semi-classical approximation to the scattering wave function under the assumption of small angular deflections of the projectile.
The DCS for scattering by a screened Coulomb potential resulting from the eikonal approximation is [14] The function is the eikonal scattering amplitude at the polar scattering angle θ for a particle of mass µ r and momentum p ′ 0 = ℏk.J 0 (x) is the Bessel function of the first kind and zeroth order, and χ(b) is the eikonal phase for projectiles incident with impact parameter b.For the analytical potential (27), the eikonal phase takes the form [28; 14] where K 0 (x) is the modified Bessel function of the second kind and zeroth order.The eikonal scattering amplitude can thus be evaluated by means of a single quadrature.Because the effect of screening decreases when the scattering angle increases (i.e., when the classical impact parameter b decreases), the DCS calculated from the eikonal approximation, Eq. ( 31), tends to the Rutherford DCS at large angles.
Although the eikonal approximation is expected to be valid for scattering angles up to about (kR at ) −1 [25], numerical calculations indicate that the approximation yields fairly accurate DCSs, practically coincident with those obtained from classical-trajectory calculations up to much larger angles, of the order of For still larger angles the calculation loses validity and presents numerical instabilities.Following Salvat [24], the DCS for angles larger than θ eik is approximated by the expression with the coefficients A, B and C obtained by matching the calculated numerical values of the eikonal DCS and its first and second derivatives at θ = θ eik .
The ratio of the calculated DCS to the Rutherford DCS, measures the effect of screening; it approximates unity at large angles (see Ref. [13]).

Elastic-scattering database
Considering that 1) the effect of screening is limited to small angles (large impact parameters), and 2) the DCS for scattering by the bare finite nucleus differs from the Rutherford DCS only at large angles (small impact parameters), it follows that screening and nuclear effects do not interfere.Hence, the CM DCS for collisions of protons and alphas with neutral atoms can be evaluated as [13] dσ The total elastic cross section is finite and given by For simulation purposes, it is convenient to consider the DCS as a function of the angular deflection of the projectile, measured by the quantity which takes values between 0 (forward scattering) and 1 (backward scattering).Notice that and We can also write where p(µ) is the normalized probability density function of µ in a single collision.
A Fortran program named panel has been written to calculate differential and integrated cross sections for elastic collisions of protons and alphas (and neutrons) with neutral atoms.This program computes cross sections for elastic collisions of a projectile particle with a given isotope Z A for the kinetic energies of the projectile specified by the user.Alternatively, it can produce a complete database of DCSs and integrated cross sections for collisions of projectiles of a given kind, with laboratory kinetic energies covering the range from 100 keV to 1 GeV for each element from hydrogen (Z = 1) to einsteinium (Z = 99).As indicated above, the atomic DCSs in the database are obtained as the average over naturally occurring isotopes of each element.
The database grid of energies is logarithmic, with 35 points per decade.For each energy the program calculates the DCS in CM, Eq. ( 37), for a grid of 1000 polar angles θ.In order to reduce the size of the database, and also to improve the accuracy of interpolation in energy, the DCS is tabulated as a function of the variable c 2 times the square of the momentum transfer in CM.The original table is "cleaned", by removing points in regions where the DCS varies smoothly, to define a reduced grid that allows accurate natural cubic spline interpolation in t.The DCS interpolated in this way is estimated to be accurate to four or more digits.For each projectile energy, the database includes the values of the total elastic cross section, Eq. ( 41), the reaction cross section obtained from Eq. ( 22) or ( 26), the first transport cross section (or momentum transfer cross section), and the second transport cross section where ⟨µ n ⟩ denotes the n-th moment of the angular deflection in a single collision.The values of these integrated cross sections serve to assess the accuracy of the DCS interpolation scheme adopted in the simulation.We recall that the total elastic cross section and the reaction cross section have the same values in the CM and L frames. Figure 1 compares results from the empirical formulas of Galyuzov and Kozov [21] with experimental data from various authors, which have been taken from the Experimental Nuclear Reaction Data (EXFOR) Database of the IAEA [29] (https://www-nds.iaea.org/exfor/).The displayed theoretical curves were obtained by assuming that the projectile and the target atom are indistinguishable, i.e., the plotted DCS describes collisions where the projectile is deflected at an angle θ together with collisions in which the recoiling target atom moves in directions with polar angle θ (or, equivalently, where the projectile emerges in directions with polar angle π − θ).Notice that, as both the projectile and the recoiling target are followed by the simulation program, the DCSs in the database are those for the scattered projectile only, which are defined for θ between 0 and π.
As indicated above, collisions of protons with nuclei of light isotopes are described by means of the optical-model potential of Watson et al. [17] for protons with kinetic energies up to 35 MeV.For higher energies, the potential of Koning and Delaroche is adopted [16].The change of model potential at 35 MeV is motivated by the comparison of results from both potentials with experimental data, as illustrated in Fig. 2.
The global potential of Koning and Delaroche [16] is expected to give a quite reliable description of elastic collisions of protons with isotopes having A > 24 (which correspond to natural elements with Z > 11).This is illustrated in Fig. 3 for collisions of protons with atoms of the isotope 208 Pb. Figure 4 compares DCSs of alpha particles with nickel atoms, 62 Ni, with the nuclear DCS calculated from the optical-model potential of Su and Han [18], which is expected to provide quite realistic DCSs for collisions of alphas with any target atom with A ≥ 20.It is worth noticing that more    reliable theoretical cross sections could be obtained by using local opticalmodel potentials (specific of each isotope) rather than the global potential models adopted here.A partial justification of the present approach for transport simulations is that collisions of charged particles much heavier than the electron are preferentially at small angles, where the DCS is mostly determined by the screened Coulomb potential of the nucleus; the details of the nuclear potential affect the DCS only for collisions with intermediate and large scattering angles, which occur with very small probabilities.

Simulation of elastic collisions
Let us assume that the projectile is moving with kinetic energy E in a compound medium whose molecules consist of n i atoms of the element with atomic number Z i (i = 1, . . ., N ).The molecular elastic DCS is obtained from the additivity approximation, i.e., as the sum of DCSs of the various atoms in a molecule, where dσ el (Z i )/dµ denotes the DCS for collisions with the element of atomic number Z i .The total elastic molecular cross sections are expressed similarly, and the ratios p i = σ el (Z i )/σ el define the probabilities of colliding with the various atoms of the molecule.In accordance with the additivity approximation, we disregard aggregation effects and, consequently, the atoms in the molecule are assumed to react as if they were free and at rest.We consider the detailed simulation of elastic collisions of the projectile with an atom of the element of atomic number Z.The kinematics of these collisions is completely determined by the polar scattering angle θ in CM.In the CM frame, after an elastic collision the magnitudes of the linear momenta of the projectile and the target atom are the same as before the collision, and the scattering angles θ, ϕ determine the directions of motion of the two particles.As mentioned above, the final kinetic energy E f and the polar scattering angle θ 1 of the projectile in the L frame are obtained by applying a Lorentz boost with velocity −v CM .The final energy of the projectile in L is     [18].Other details as in Fig. 3.
with the energy loss W given by where is the maximum energy loss in a collision, which occurs when θ = π.The polar angle θ 1 of the final direction of the projectile in L is given by with and where is the speed of the scattered projectile in CM.Notice that the azimuthal angle of the projectile direction in L is the same as in the CM frame.After the collision, in the L frame the target atom recoils with kinetic energy E A = W and direction in the scattering plane with the polar angle θ A given by cos In class II simulations [8; 9] it is necessary to consider the contribution of soft elastic collisions to the elastic transport cross sections and to the stopping cross section.The required quantities are determined by the angular DCS in the L frame and by the energy-loss DCS associated to elastic collisions.The angular DCS is expressed in terms of the scattering angles in the L frame by making use of the inverse of the relation (51), If τ is less than, or equal to unity only the plus sign before the square root has to be considered.For τ > 1, there are two values of the CM deflection θ, given by Eq. ( 56), for each value of θ 1 , which correspond to different final energies of the projectile in L. The DCS in the L frame is given by where the last factor is the DCS in the CM frame.From the relation (56), we obtain (a derivation of this expression is given in the supplementary document) If τ < 1 only the plus sign is valid and the scattering angle θ 1 varies from 0 to π.When τ ≥ 1, the DCS in L vanishes for angles θ 1 larger than for angles θ 1 < θ 1,max , Eq. ( 56) yields two values of θ in (0, π), the expression on the right-hand side of Eq. (58) must then be evaluated for these two angles (with the corresponding plus or minus sign in the numerator), and the resulting values added up to give the DCS in L.
The energy-loss DCS is and the so-called nuclear stopping cross section is given by where σ el,1 is the first transport cross section in the CM frame, Eq. ( 44).The simulation of elastic collisions is performed by using the same strategy as in the penelope and penh codes [24; 8].Mean free paths and other energy-dependent quantities are obtained by log-log linear interpolation of tables, prepared at the start of the simulation run, with a logarithmic grid of 200 laboratory kinetic energies E i that covers the interval of interest.The angular distribution of scattered projectiles in CM, is tabulated at the same grid energies.
The CM scattering angle θ of a projectile with laboratory energy E in the interval (E i , E i+1 ] is sampled from the distribution with which is obtained from the tabulated distributions by linear interpolation in ln E. The sampling is performed by using the composition method: 1) select the value of the index k = i or i+1, with respective point probabilities π i and π i+1 , and 2) sample µ from the distribution p(E k , µ).
With this interpolation by weight method, µ is generated by sampling from only the distributions at the grid energies E i .This sampling is performed by the inverse transform method by using the RITA (rational interpolation with aliasing) algorithm [30; 8].The required sampling tables are prepared by the program at the start of the simulation run.

Inelastic collisions
Let us now consider the description and simulation of inelastic collisions of charged particles, i.e., interactions of the projectile that result in electronic excitations of the material.The most probable effect of inelastic collisions is the excitation of weakly bound (valence or conduction) electrons of the material, which can be described by means of the relativistic plane-wave Born approximation (PWBA) [31; 32].Notice that the wave functions of weakly bound electrons are strongly affected by the state of aggregation of the material and, hence, a realistic description of the response of the material requires the use of empirical information.The interaction model described here accounts for the dependence on the microscopic structure of the material by using the empirical value of the mean excitation energy I [33], which determines the stopping power for high-energy projectiles.
Formally, the adopted model is analogous to the one employed in penelope for inelastic collisions of electrons and positrons, which is slightly modified to yield a finite stopping power for slow projectiles.We disregard the fact that the mass of the target is finite and, consequently, inelastic collisions are described in the laboratory frame, where the stopping material is at rest.For the sake of generality, we consider a molecular medium, with Z M electrons in a molecule.Its electronic structure is described as a number of bound electron subshells, each with f k electrons and binding (ionization) energy U k , which essentially retain their atomic properties, and, in the case of conducting materials, a set of f cb nearly free electrons in the conduction band, with U cb = 0.By construction, Individual inelastic collisions of a projectile (mass M 1 and charge Z 1 e) with kinetic energy E and linear momentum p are conveniently characterized by the energy loss of the projectile, W = E − E f , and the magnitude q of the momentum transfer q ≡ p − p f , where E f and p f are, respectively, the kinetic energy and the linear momentum of the projectile after the interaction.Notice that (cp and To simplify the form of the DCS, it is customary to introduce the so-called recoil energy, Q, which is defined as the kinetic energy of an electron with momentum equal to the momentum transfer [31], in other words, where θ = arccos(p• pf ) is the polar scattering angle.Equivalently, The doubly-differential cross section (DDCS), differential in W and Q, can be expressed as (see, e.g., [31; 8]) and where df (Q, W )/dW is the generalized oscillator strength (GOS), which completely characterizes the response of the material.The first term in expression (68) describes excitations caused by the instantaneous Coulomb interaction; the second term accounts for excitations induced by the transverse interaction (exchange of virtual photons).We should mention that the transverse contribution in Eq. ( 68) results from the approximation of neglecting the differences between longitudinal and transverse GOSs (see, e.g., [31; 34; 35]).These differences are negligible for small Q, which dominate in transverse interactions, as well as for large Q.
For a given energy loss W , the allowed values of the recoil energy lie in the interval (Q − , Q + ), with endpoints given by Eq. (66) with cos θ = +1 and −1, respectively.In other words, When W ≪ E, the lowest allowed recoil energy can be calculated from the approximate relation [36] Q Note that the curves which, when W ≪ E, reduces to It follows that, for given values of E and Q [< Q + (0)], the only kinematically allowed values of the energy loss are those in the interval 0 < W < W m (Q).The energy-loss DCS is defined by The probability distribution function (PDF) of the energy loss in a single inelastic collision is given by where is the total cross section for inelastic interactions.It is convenient to introduce the quantities where ⟨W n ⟩ denotes the n-th moment of the energy loss in a single collision (notice that σ in and σ in are known as the stopping cross section and the energy-straggling cross section, respectively.
The mean free path λ in for inelastic collisions is where N is the number of molecules per unit volume.The electronic stopping power S in and the energy straggling parameter Ω 2 in are defined by and respectively.The stopping power gives the average energy loss per unit path length.The physical meaning of the straggling parameter is less direct; the product Ω 2 in (E) ds gives the variance of the energy distribution of charged projectiles that start moving with energy E after traveling a (small) distance ds within the medium.

The generalized oscillator strength model
Although realistic GOSs may be available for simple systems, given either by analytical formulas (hydrogenic approximation [32] and electron gas [37]) or by numerical tables (obtained, e.g., from DHFS calculations for atoms [34; 35]), they are not suited for general-purpose Monte Carlo simulations, mostly because of the strong correlations between the variables W and Q.To account for these correlations, we should sample the two quantities from their joint PDF, i.e., from the DDCS, a process that requires massive memory storage and accurate interpolations.
Here we use a model of the GOS, adapted from the penelope code [38; 8], that reproduces the most conspicuous features of the GOS, satisfies relevant sum rules, and provides exact analytical formulas for sampling W and Q in individual interactions.Excitations of electrons in a subshell k with binding energy U k are described as a single "oscillator" or one-electron GOS, F k (Q, W ), defined as where with The quantity b (> 0) is a free parameter; a comparison with calculated subshell ionization cross sections by means of the PWBA with the DHFS potential [35] (see Fig. 7 below) indicates that a value b ∼ 4 is adequate.Notice that The first term in expression (82) represents low-Q (distant) interactions, which are described as a single resonance at the energy W k .The second term corresponds to large-Q (close) interactions, in which the target electrons react as if they were free and at rest (W = Q); close interactions are allowed only for energy transfers W larger than U k .It is worth noticing that in the case of conductors the model can be used for describing the GOS of the conduction band (with U cb = 0), and the resulting stopping power only vanishes at E = 0. Figure 5 displays schematically the model GOSs for inner subshells and for the conduction band.The molecular GOS is the sum of contributions for the various electron shells of the atoms in a molecule, where f k is the number of electrons in the k subshell.For bound shells, the resonance energy is defined as where is the plasma energy of a free electron gas with the electron density of the medium, and a is an adjustable parameter, the so-called Sternheimer factor.The term 2f k Ω 2 p /(3Z M ) in expression (86) accounts for the Lorentz-Lorenz correction (the resonance energies in a condensed medium are larger than those of isolated atoms or molecules).In the case of conductors, excitations of the conduction band are represented by a single oscillator with oscillator strength f cb equal to the number of free electrons per molecule, null binding energy (U cb = 0), and resonance energy Note that W cb is the plasmon excitation energy of a free-electron gas with the electron density of the conduction band.When a material is qualified as a conductor, f cb is set equal to the average lowest negative valence of the elements present (f cb = 0 for insulators).For free-electron-like materials, such as metallic aluminum, the value (88) is close to the energy of volume plasmons.
The GOS model (85) satisfies the Bethe sum rule, for all Q.In the limit Q → 0 the GOS reduces to the optical oscillator strength (OOS), which characterizes the optical properties of the medium, and determines the density effect correction to the stopping power of highenergy particles.Indeed, the OOS resulting from our GOS model, with the resonance energies (86), coincides with the OOS assumed by Sternheimer et al. [40; 41] in their calculations of the density effect correction.The Sternheimer factor a is fixed by requiring that the GOS model leads to the empirical value of the mean excitation energy I of the material [42], Thus, the GOS is completely determined by the mean excitation energy I, which is the only free parameter of the model.By default the simulation code uses I values from the ICRU Report 37 [42].Typical values of the Sternheimer factor range between about 2 and 3.The requirements ( 89) and (91) ensure that the stopping power of high-energy particles coincides with the values given by the Bethe formula [43].

Differential and integrated cross sections
The GOS completely characterizes the response of individual molecules to inelastic interactions with the projectile (within the PWBA).The molecular DDCS can be expressed as where d 2 σ k /(dQ dW ) is the DDCS for excitations of a single electron described by the one-electron GOS F k (Q, W ). Hereafter the summation over oscillators k includes a term corresponding to the conduction band, with oscillator strength f cb , resonance energy W cb , and ionization energy equal to zero.The DDCS for collisions with an oscillator is conveniently split into contributions from close collisions and from distant (resonant) longitudinal and transverse interactions, The DDCSs for close collisions and for distant longitudinal interactions are, respectively, and The quantity W ridge is the maximum energy loss in collisions of the projectile with free electrons at rest, which is given by Notice that, when M 1 = m e , W ridge = E.For projectiles heavier than the electron (M 1 ≫ m e ) with kinetic energies much less than their rest energy M 1 c 2 , R ∼ 1 and The response of molecules in a dense medium is modified by the dielectric polarization of the material, which modifies the distant transverse interactions and causes a reduction of the stopping power known as the density-effect correction.The DDCS for distant transverse interactions is approximated as (98) where δ F is the density-effect correction to the stopping power.It is worth mentioning that this approximate DDCS results from 1) neglecting the angular deflection of the projectile in distant transverse interactions, which is generally very small, and 2) requiring that it gives the exact contribution of the distant transverse interactions to the stopping power for high-energy projectiles, in accordance with the corrected Bethe formula for the stopping power [33].
The quantity δ F is calculated as [44; 8] where L is a real-valued function of β 2 defined as the positive root of the equation The function F(L) decreases monotonically with L, and hence, the root L(β 2 ) exists only when 1 − β 2 < F(0); otherwise δ F = 0.In the high-energy limit (β → 1), the L value resulting from Eq. ( 100) is large (L ≫ W k ) and can be approximated as L 2 = Ω 2 p /(1 − β 2 ).Then, using the Bethe sum rule (89) and the relation (91), we obtain The energy-loss DCS for collisions with the k-th oscillator, can also be split into contributions from close, distant longitudinal, and distant transverse interactions, where and These energy-loss DCSs, as well as the one-electron cross sections integrated over an arbitrary interval (W 1 , W 2 ), can be evaluated analytically (see the supplementary document).Evidently, the molecular integrated cross sections for inelastic collisions are σ Figure 6 compares the electronic stopping powers of aluminum, silver, and gold for protons and alpha particles calculated from the present GOS model with realistic values obtained by means of the program sbethe of Salvat and Andreo [43], which uses a corrected Bethe formula.This comparison illustrates our claim that the stopping power obtained from the GOS model effectively tends to the realistic value for high-energy projectiles.

Integrated angular cross sections
Inelastic collisions cause small deflections of the projectile and contribute to the directional spreading of particle beams when they penetrate matter.
For simulation purposes, it is convenient to describe angular deflections by means of the variable µ, Eq. ( 39), instead of the polar scattering angle θ.The recoil energy Q, Eq. ( 66), can then be expressed as In distant interactions with the k-th oscillator, W = W k and the magnitude p f,k of the linear momentum of the projectile after the collision, is fixed, which implies that µ is a function of Q only.In close collisions (110) The total angular cross section, the first transport cross section, and the second transport cross section for inelastic collisions with the k-th oscillator are defined, respectively, as and where dσ in /dµ is the DCS, differential in the deflection µ.Naturally, both the differential and the integrated angular cross sections per molecule are the sums of contributions from the various oscillators, The contribution of close collisions with the k-th oscillator to the integrated angular cross sections can be calculated in terms of the energy-loss DCS, while that of distant longitudinal interactions is conveniently calculated in terms of the DCS differential in the recoil energy, Distant transverse interactions do not contribute to the transport cross sections because the projectile is not deflected in those interactions.In the simulation program, the integrals in Eqs.(111) are calculated numerically (details of this calculation are given in the supplementary document).

Near-threshold distant interactions
The details of the oscillator GOS model have been tailored to allow exact random sampling of the energy loss W and the recoil energy Q.In addition, the model can be used for describing interactions with both bound electrons and conduction electrons.An exact sampling algorithm, which keeps the correlations between Q and W embodied in the GOS model, is described in the supplementary document.
Each inelastic interaction with the k-th oscillator causes the release of a secondary electron with kinetic energy E s = W − U k in the direction of the momentum transfer, defined by the polar angle θ r given by Eq. (70).
In the case of excitations of a bound subshell, the energy loss distribution associated with distant interactions is described as a single resonance (delta function), while the actual distribution is continuous for energy losses above the ionization threshold.As a consequence, energy loss spectra simulated from the present GOS model will show unphysical narrow peaks at energy losses that are multiples of the resonance energies.To get rid of this kind of artifact, we spread the resonance line by sampling the energy loss in distant interactions from the continuous triangular distribution in the interval from That is, we consider the distribution which gives the correct average value, ⟨W ⟩ = W k (see Fig. 5).Since energy losses larger than W m (Q c ) are forbidden, the value of W d should be smaller than W m (Q c ).When this is not the case, we modify the resonance energy W k , and replace it with the value That is, the quantity W k is replaced with this modified value in all formulas pertaining to the distant excitations of bound subshells.Also, to prevent an anomalous increase of the ionization cross section of bound subshells for projectiles with kinetic energy near the threshold, we multiply the DCS for distant excitations by the factor which reduces to unity when Thus, the maximum allowed energy loss in distant excitations of bound subshells, Eq. (114), is given by which never exceeds W m (Q c ).The energy loss in distant excitations is sampled from the pdf (115) by using the sampling formula where ξ is a random number uniformly distributed in (0,1); this formula results from the inverse transform method [8].The spread distribution and the lowenergy modification of the resonance energy are applied only to bound electron subshells.The energy spectrum of distant interactions with conduction-band electrons is not altered, i.e., the energy loss in these excitations equals W cb independently of the energy of the projectile.

Ionization of inner subshells and re-normalization
The GOS model given by Eq. (85) provides a quite realistic description of the correlations between the energy loss and the scattering angle in inelastic collisions of charged particles.However, the subshell total cross section obtained from that GOS model may differ appreciably from results of experiments and of more accurate calculations.Inaccuracies in the total cross section for ionization of inner electron subshells become apparent when we consider the emission of x rays induced by impact of charged particles: the number of x rays emitted is proportional to the ionization cross section of the active subshell.
To provide a more accurate description of the emission of x rays and Auger electrons, we have calculated a complete database of cross sections for ionization of inner subshells (K shell, L, M, and N subshells with binding energy larger than 50 eV) of all the elements from hydrogen to einsteinium (Z = 1 to 99), by impact of protons and alpha particles with energies up to 10 GeV.The calculations were based on the relativistic PWBA, as formulated by Bote and Salvat [34] (see also [35]) using longitudinal and transverse GOSs computed with the DHFS potential.Following Chen [45] and Chen and Crasemann [46], we adopted the perturbed-stationary-state approximation of Brandt and Lapicki [47], which improves the PWBA by accounting 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.Details of these calculations are described by Salvat [39].Chen and Crasemann [46] performed similar calculations using the non-relativistic PWBA, also with GOSs obtained from the DHFS potential, and published tables of cross sections for ionization by protons with energies up to 5 MeV.Our results agree closely with theirs, but extend to much higher energies.In addition, to approximately account for the density effect, we reduce the cross sections in the database by a factor equal to the ratio of the cross sections obtained from the GOS model with and without the density effect correction, δ F .Hereafter, the ionization cross section of our calculated database, with this density-effect correction factor, will be referred to as "reference" ionization cross sections.
In our simulation program, the total cross section, σ in , is decomposed into contributions from inner and outer electron subshells, where the first summation is over inner subshells (i.e., K to N7 subshells with binding energies U i greater than the cut-off energy E cut = 50 eV); the second summation is over outer subshells (i.e., those with U j < E cut or with principal quantum number larger than 4). Figure 7 compares the reference ionization cross sections of the inner shells of the cobalt atom (Z = 27) with the predictions of our GOS model for solid cobalt.The various curves correspond to the indicated subshells; notice that σ in,i tends to increase when the binding energy of the active subshell decreases.As the total cross section and the stopping cross section are dominated by contributions from outer subshells with relatively small binding energies, the total cross sections of inner subshells may be modified, up to a certain extent, and those of the outer subshells may be re-normalized so that the input stopping power remains unaltered.The simulation program assumes that hard inelastic collisions with inner subshells ionize the target atom, and the relaxation of the resulting vacancies is simulated by the penelope routines by using the transition probabilities given in the Evaluated Atomic Data Library of Perkins et al. [48].To get the correct number of emitted x rays, the total cross section of each inner shell, f i σ in,i (E), is replaced with the reference cross section σ (ref) in,i (E), without altering the details of the PDF of the energy-loss and scattering angle.That is, the "oscillator strength" f i of the i-th inner shell is replaced with when σ in,i (E) > 0. It is worth noticing that because of the neglect of the motion of atomic electrons in close collisions, the GOS model gives effective ionization thresholds that are higher than those of the reference cross sections.That is, we may have σ in,i (E) = 0 but σ in,i (E) ̸ = 0, in which case the projectile particles can ionize the inner shell at energies lower than the corresponding ionization threshold; under these circumstances, the energy transfer is set equal to the binding energy of the subshell, W = U i , and the projectile's trajectory is not deflected.Of course, this procedure implies increasing the inner-subshell contribution to the stopping power in the (small) quantity The program reads a table of the stopping power, S in (E), from the input material-data file, which is considered to be the actual stopping power of the material.By default, this table is calculated from the GOS model (85) as described above.In order to avoid altering the input stopping power, the total cross sections of outer subshells, f j σ j (E), are multiplied by an energy-dependent scaling factor, N (E), the same for all outer subshells, given by where σ in,j (E) is the one-electron stopping cross section for excitations of the j-th outer subshell, Eq. (106).Formally, this modification is equivalent to replacing the oscillator strengths f j of the outer subshells with f ′ j = N (E) f j .As already mentioned, by default the input stopping power is calculated from the PWBA with the GOS model (85).However, the PWBA with the density-effect correction is valid only for projectiles with relatively high energies.Departures from the PWBA give rise to the Lindhard-Sørensen and Barkas corrections to the Bethe formula [43].To account for these departures, the user may edit the input material-data file and replace the stopping power table with more reliable values.As reference stopping powers one may use those generated by the program sbethe of Salvat and Andreo [43], which are consistent with the recommendations and values given in the ICRU Report 49 [33].

Tracking algorithm
The interaction models described above permit the formulation of a class-II tracking scheme [8; 9] with a fixed energy-loss cutoff W cc , which is set by the user, and an energy-dependent cutoff deflection µ c for elastic collisions that is defied internally by the program in terms of two user-defined simulation parameters, C 1 and C 2 .Particle trajectories are generated by using the random-hinge method [7; 8], which operates similarly to detailed simulations, i.e., the transported particle is moved in straight "jumps", and the energy and direction of movement change only through discrete events (hard interactions and hinges).Here we sketch the simulation algorithms briefly, additional details can be found in the manual of the code system penelope and in the article by Asai et al. [9].

Elastic collisions
In our simulation code the cutoff deflection µ c , which separates hard and soft elastic collisions, is determined by two energy-independent user parameters, C 1 and C 2 , which typically should be given small values, between 0 and 0.2.These two parameters are used to fix the mean free path between hard elastic events (i.e., the average step length between consecutive hard elastic collisions), which is defined as where λ el,1 = [N σ el,1 ] −1 is the first transport mean free path, see Eq. (44), and is the CSDA range calculated from the input electronic stopping power.The identity then fixes the cutoff µ c as a function of the energy E of the projectile, which may be different for the various atoms in a molecule.The recipe (123) forces limits the average fractional energy loss along the step.An increase of C 1 or C 2 leads to increased values of both the mean free path between hard events, λ el , and the cutoff deflection, µ c , in certain energy ranges [8].Of course, an increase of λ (h) el implies a reduction in the number of hard events along a particle track with an accompanying reduction of the simulation time.
The angular deflection effect of the soft interactions that occur between each consecutive pair of hard interactions is determined by the transport cross sections of orders ℓ = 0 and 1 of the soft interactions in the L frame.The contributions from elastic collisions are where µ 1 is the angular deflection in the L frame.It is important to notice that soft inelastic collisions also cause a small deflection of the projectile.The scattering effect of these interactions is accounted for by considering their contributions to the soft transport cross sections, where is the sum of contributions of all oscillators restricted to energy losses less than W cc .The combined (elastic plus inelastic) soft scattering process is then described by the transport mean free paths of orders ℓ = 1 and 2. Assuming that the energy loss is small, the first and second moments of the angular deflection after a path length s, under the sole action of soft elastic and soft inelastic interactions, are [6; 8] In practical simulations the angular deflection µ s after a path length s is sampled from an artificial distribution, P (µ s ), which is required to have the same moments, of orders n = 1 and 2 as the real distribution, Eqs. ( 131), but is otherwise arbitrary [8; 9].

Inelastic collisions
As indicated above, the simulation of inelastic collisions is tuned by the cutoff energy transfer W cc set by the user, which separates soft and hard interactions.Hard inelastic interactions with energy-loss higher than W cc are simulated individually from the corresponding restricted DDCS.To simplify the programming, distant interactions with an oscillator are considered to be hard only if U k ≥ W cc , i.e., distant excitations of oscillators with U k < W cc are all soft.This classification avoids the need of splitting the continuous distribution (115).The sampling of hard interactions is performed exactly by using the algorithms described in the supplementary document, modified so as to deliver energy losses larger than W cc .Along each trajectory step (to or from a hard interaction), soft interactions with W < W cc may occur.The cumulative effect of these soft interactions is described by means of a multiple scattering approach determined by the restricted stopping power, A difficulty of class-II algorithms arises from the fact that the energy of the particle decreases along the step between two consecutive hard interactions.
Because the cutoff energy W cc does not change with E, we can assume that, at least for small fractional energy losses, the DCSs for soft energy-loss events vary linearly with E. Under this assumption we can calculate the first moments of the distribution of the energy loss W s of a particle with initial energy E 0 after traveling a path length s under only the influence of soft events [8].The mean and variance of this distribution are, respectively, The energy loss caused by soft events along a trajectory step is sampled from an artificial pdf with parameters obtained from the stopping cross section and the energy-straggling cross section for soft interactions [8].The accumulated angular deflection caused by soft interactions along a step is sampled from an artificial distribution with its first and second moments determined by the first and second transport cross sections restricted to soft interactions.These integral characteristics of soft interactions are readily obtained from the expressions given above with the appropriate limits of the integrals.

Concluding comments
We have presented DCSs for elastic and inelastic collisions of protons and alpha particles suited for class-II Monte Carlo simulations of the transport of charged particles in matter.The DCS for elastic collisions are calculated from realistic nuclear optical-model potentials by using highly accurate partial-wave methods, and corrected to account for the effect of screening of the nuclear charge by the atomic electrons.Atomic DCSs in the CM frame have been calculated for the elements with atomic numbers 1 to 99; they have been included in an extensive database for protons, alpha particles (and neutrons) with kinetic energies between 100 keV and 1 GeV.
Inelastic collisions are described by means of the PWBA, in order to provide a description of electron binding effects and of the correlations between the energy loss and the deflection angle of the projectile in inelastic events.The proposed GOS model satisfies the Bethe sum rule, and partially incorporates the effect of aggregation by using an empirical value of the mean excitation energy I as a defining parameter.As a consequence our DCSs lead to the correct electronic stopping for high energy projectiles.A simple renormalization of the DCS of inner subshells, to agree with ionization cross sections calculated with the DHFS self-consistent potential, ensures that simulations will generate the correct number of ionizations and the ensuing emission of x rays and Auger electrons.In addition, a further renormalization of the DCSs of outer electron subshells permits incorporating more realistic stopping powers for projectiles with intermediate and low energies.
The proposed interaction models can be used in class-II simulations of charged-particle transport.They permit the formulation of adequate sampling algorithms for hard interactions, i.e., elastic collisions with angular deflections larger than µ c and inelastic collisions with energy loss larger then W cc , with arbitrary cutoffs.An exact sampling algorithm for inelastic collisions is described in the supplementary document.These models and databases have been implemented in a Fortran simulation code named penhan that, in conjunction with penelope [8], simulates the coupled transport of electrons, positrons, photons, protons, and alpha particles in matter.A detailed description of penhan, which is available from the authors under request, will be published elsewhere.

Figure 4 :
Figure 4: Elastic DCS in the CM frame for collisions of alpha particles with neutral atoms of nickel, 62 Ni.The solid curves represent results from partial-wave calculations with the global optical-model potential of Su and Han [18].Other details as in Fig. 3.

Figure 5 :
Figure 5: Oscillator model for the subshell GOS, represented by the solid lines with thickness proportional to the GOS value.The continuous curve is the maximum allowed energy loss as a function of the recoil energy, W m (Q), Eq. (73) for protons with E = 5 keV.(a) GOS of a bound subshell with U k = 1 keV.For distant interactions the possible recoil energies lie in the interval from Q − to Q c , and the energy loss W varies between U k and W d , Eq. (114).(b) Oscillator-GOS model for excitations of the conduction band of conductors (U cb = 0).

Figure 6 :
Figure 6: Stopping power of inelastic collisions S in /ρ for protons and alpha particles in aluminium, silver (×10) and gold (×100) as a function of the kinetic energy of the projectile.Solid curves are results from the present GOS model.Dashed curves are results from the corrected Bethe formula implemented in the program sbethe [43].

Figure 7 :
Figure 7: Ionization sections of the inner subshells of cobalt atoms by impact of protons and alphas, as functions of the kinetic energy of the projectile.Solid curves represent the reference ionization cross sections obtained from the accurate calculations described in the text.Dashed curves are the predictions from the present GOS model for solid cobalt.

S 0 W 2
of numerical consistency, we also include the stopping due to soft elastic collisions, which accounts for energy transfers W = W max µ to recoiling target nuclei (nuclear stopping) , dσ el (Z, E) dµ dµ ,(136)where both W max , Eq. (50), and µ c , Eq. (125), are specific of each target element.The global stopping power and energy-straggling parameter of soft interactions are S s (E) = S

S
s (E 0 ) s , (138b) where the factors in curly braces account for the global effect of the energy dependence of the soft energy-loss DCS, within the linear approximation.