Ionization and excitation cross sections for the interaction of HZE particles in liquid water and application to Monte Carlo simulation of radiation tracks

Relativistic heavy ions of high charge (Z) and energy (E) (HZE) in galactic cosmic rays (GCR) are important contributors to space radiation risk because they cannot be shielded completely and their relative biological effectiveness is very high. To understand these risks, Monte Carlo track structure simulations by radiation transport codes are widely used in radiation biology to provide information on energy deposition and production of radiolytic species that damage cellular structures. In this paper, we show that relativistic corrections can be applied to existing semi-empirical cross section models for the ionization and excitation of water molecules by ions to extend the validity of their energy range up to ∼104 MeV amu−1. Similarly, an effective charge value correction is applied for Z>2 ions. Simulations of HZE tracks have been performed by a new C++ Monte Carlo transport code, named RITRACKS, that uses these cross sections to calculate the stopping power, radial dose, XY-plane projections of track segments and radial distributions of primary radiolytic species (H•, •OH, H2, H2O2 and e−aq) at ∼10−12 s. These new data will be useful to understand results from experiments performed at ion accelerators by discriminating the role of the so-called core and penumbra of the tracks.


Introduction
Astronauts are exposed to galactic cosmic rays (GCR), comprising ions of which the nuclear charge extends from hydrogen to uranium. The energy spectrum of the GCR peaks near 1 GeV per nucleon, and these particles are so penetrating that shielding can only partly reduce the doses absorbed by the astronauts. Predictions of the nature and magnitude of risks posed by exposure to such radiation in space are subject to many uncertainties [1]. The large ionization power of high charge (Z ) and energy (E) (HZE) ions with Z > 2 makes them the major contributors to the risk, even if their fluence is lower than protons. These HZE ions are usually studied up to Z = 26 (iron), because heavier ions are infrequent. As space experiments are difficult, expensive and poorly reproducible, ground experiments performed in facilities like the NASA Space Radiation Laboratory (NSRL) at the Brookhaven National Laboratory (Upton, NY, USA) are key to get experimental data of great relevance for reducing uncertainty on risk assessment. Such experiments have shown that the relative biological effectiveness of HZE ions can be as high as 30 for chromosomal exchanges in interphase by energetic iron ions [2] or for tumors in mice (reviewed in [1]).
Monte Carlo track structure simulations are widely used in radiation biology to provide detailed information on the energy deposition and on the production of radiolytic species that damage DNA and other cellular structures [3]. These simulations are often limited to light ions, but there is interest to extend these simulations to HZE particles for the reasons discussed above. To perform these simulations, total and differential ionization and excitation cross sections are needed. While considerable experimental data are available for water vapor, measurements for the liquid phase are very difficult and are largely lacking [4,5].
In general, the inelastic interaction cross section for charged particles with condensed matter are calculated within the framework of the dielectric response theory [6]. This theory takes into account medium polarization effects like the Fermi density effect. In this paper, we show that semi-empirical models developed previously can be used as an alternative to obtain ionization and excitation cross sections for ions in the relativistic energy range by introducing simple relativistic corrections to extend the energy range up to 10 GeV amu −1 . The cross sections are also extended to higher Z ions by scaling the cross section by a factor Z * 2 , where Z * < Z is an effective charge that takes into account the fact that heavy ions at low energies have electrons attached. Simulations of HZE tracks have been performed by a new C + + program, named RITRACKS, that uses these cross sections to calculate stopping power, radial dose, X Y -plane projections of track segments and radial distribution of primary radiolytic species (H • , • OH, H 2 , H 2 O 2 and e − aq ) at ∼10 −12 s.

Ionization cross sections by ions
Rudd [7] has proposed a semi-empirical equation for the differential ionization cross section of the molecular orbitals 1 b 1 , 3 a 1 , 1 b 2 , 2 a 1 and 1 a 1 of liquid water by protons: Here, i is the index of the orbital, I i is the binding energy of the electron in the target, w = W/I i , W is the energy of the secondary electron, E p is the energy of the incident proton, T = E p (m/M p ) is the energy of an electron of mass m which would have the same velocity as a proton of mass M p of the same energy and v = √ T /I i is a scaled velocity of the incident particle. Furthermore, S i = 4πa 2 0 N i (R/I i ) 2 , a 0 = 5.3 × 10 −11 m is the Bohr radius, R = 13.6 eV is the Rydberg energy, N i is the number of electrons of the orbital and w c i = 4v 2 − 2v − R/4I i is the scaled cutoff energy. The functions F 1 (v) and F 2 (v) are defined as: where The ten parameters A 1 , B 1 ,..E 1 , A 2 , B 2 ..D 2 and α have been determined to be the best fit of the experimental data and are shown in table 1. The total ionization cross section is found by numerically integrating Here, E max is the maximum energy transfer in a single collision. In classical mechanics, the maximum energy transferred to a resting electron of mass m by an ion of mass M and energy E is given by [8]: The semi-empirical equations for the ionization cross sections are valid in the classical energy range. To extend this formulation to relativistic energy range, v 2 is replaced with where c is the speed of light. Using Taylor's series, one can verify that this expression reduces to v 2 = T /I i when T mc 2 . Similarly, E max is replaced by its equivalent relativistic expression [8] An ion of energy per nucleon τ has the same velocity of a proton of a given energy [6]; thus, the relativistic quantities β and γ are conviently expressed as a function of τ : Except at extreme relativistic energy, γ m/M 1 and E max reduces to The integration (8) has been performed using the parameters for liquid water and is shown in figure 1. The results are compared with the calculations of Dingfelder et al [9]. For energies < 10 MeV, Dingfelder et al used Rudd's equation, so the calculations give exactly the same results. For energies >10 MeV, v 2 reduces to a constant value and the cross section becomes constant. Dingfelder et al [6] used cross sections based on the dielectric response theory. There is a slight difference with our results, but also with Rudd's cross section at 10 MeV.  [11] for proton energy <10 MeV and Dingfelder et al [6] for proton energy >10 MeV are also shown for sake of comparison. The data from each molecular orbital from Dingfelder et al [6] are not shown for clarity.
For track structure simulation, we need to generate secondary electrons that are distributed energetically as dσ/dW . It can be written as This distribution has the form where k is a constant at least equal to 1, g(x) is a normalized distribution that can be generated by inversion and the so-called rejection function ψ(x) is [0, 1] valued. Random numbers that follows such a distribution can be generated by a rejection method [10]. The function g(x) can be generated from a random number U unifomly distributed between 0 and 1 by the direct inversion method: where The normalization constant c is given by the integration of g(w) between 0 and w max . The rejection function is The algorithm described previously has been used 10 6 times to sample the energy loss for 0.1 MeV protons, for all the molecular ionization orbitals. The results are recorded in a normalized histogram, and shown in figure 2. The dots are the results from the sampling, whereas the lines are the analytical predictions given by equation (16).

Excitation cross sections by ions
The excitation cross section is roughly one order of magnitude less than the ionization cross section, except at low energy (<50-100 keV). The relative contributions of each discrete electronic excitation level so-calledÃ 1 B 1 ,B 1 A 1 , Ryd A + B, Ryd C + D (Rydberg series), diffuse 10.1 21.3 bands and plasmon excitation (collective excitation) is currently a matter of debate [11]. The expression of these excitation levels can be calculated by using the semi-empirical formulation shown in Dingfelder et al [9] or in Kutcher and Green [12]. In our program, for reasons that are discussed in Cobut et al [13], only the excitation levelsÃ 1 B 1 ,B 1 A 1 and plasmon excitation are considered. The differential cross section can be written: where Here, W is the energy loss by excitation, ρ(W ) is the differential cross section for charged particle scattering on free electrons at rest, T = E(m/M p ), a 0 is the Bohr radius, ε 0 is vaccuum permittivity, m and e are the mass and charge of the electron; M p , E, v and z are the mass, energy, velocity and charge of the incident proton. Note that the definition of W and v are different from the previous section.
The two excitation levels are modeled by a Gaussian function: The plasmon is described by the function The parameters f 0,i , α i and W 0,i are shown in table 2.
The excitation cross section is calculated by integrating (21): The integration limits are set empirically to 0.01 and 100 eV. It is necessary to have W > 0 to avoid the divergence of the integral. However, the functions f i are peaked and these integration The excitation cross sections for protons are shown in figure 3. The results of our calculations are compared with those from Dingfelder [9]. Even though the models are very different, the calculations of the cross sections yield similar results for the energy range 0.1-100 MeV. At high energy, the empirical formula of Dingfelder is non-relativistic, because the cross section continues to decrease over 100 MeV. The energy loss due to excitation is much smaller than the one for ionization. For example, the mean excitation energy of electrons in the energy range 10-10 7 eV varies slowly from ∼8-14 eV in water vapour [14]. This is in agreement with our values, except for the energy loss by plasmon excitation which is in the range 20-30 eV.
The differential cross section (21) for excitation levelsÃ 1 B 1 andB 1 A 1 is the product of a Gaussian function ψ(W ) with a function g(W ): Here, since dσ/dW has the form kg(W )ψ(W ), it is also possible to use a rejection algorithm to generate energy loss that follows this distribution. The Gaussian functions are generated by the Box-Muller method [10]. The rejection function is equation (29). Similarly, to sample the energy loss in a plasmon excitation event, a rejection method is also used. Random numbers distributed as f 0,pl are generated by an inversion method [10]. The sampled values W s are given by: where u 1 , u 2 and C are: Figure 5. Calculation of the LET as a function of the energy for 1 H + , 4 He 2+ , 12 C 6+ , 20 Ne 10+ and 56 Fe 26+ using the cross sections shown above and including charge effects and relativistic corrections.
The rejection function is also given by equation (29). These two algorithms have been used 10 6 times to generate W s values for excitation levelsÃ 1 B 1 ,B 1 A 1 and plasmon excitation by 1 MeV protons, and the values of W s are stored in a normalized histogram. The results are shown in figure 4. The dots are the results from the algorithm sampling, whereas the lines are the analytical predictions.

HZE tracks simulations
Within the first-order plane wave Born approximation, the interaction cross section for bare heavy ions of velocity v are obtained from proton interaction cross sections by scaling the differential cross section for a proton with the same velocity v by the square of the charge Z of the ion [6], A heavy ion of mass M with kinetic energy E ion = (M/M p )τ has the same velocity of a proton (mass M p ) with kinetic energy τ . This is true for both non-relativistic and relativistic ions. It means that the empirical equations of the ionization and excitation cross sections can be used for an ion of energy per nucleon τ , but multiplied by a scaling factor Z 2 . Heavy ions at small energies have electrons attached and thus have a reduced charge Z * < Z . If Z * is defined to give the correct observed stopping power, it is not equal to the mean charge per particle of a beam leaving an absorber [15]. An equation from Booth and Grant [16] is commonly used for the effective charge: This correction factor becomes important (>5%) for ions with β < 0.0217Z 2/3 , which is ∼18 MeV amu −1 for Fe ions and ∼0.22 MeV for protons. Monte Carlo simulations of HZE tracks have been performed with a new C + + transport code, named RITRACKS, using methods described in Cobut et al [13] but using the relativistic cross sections and the effective charge value [17]. Briefly, the distance between two energy loss events is found by sampling an exponential distribution with parameter λ = [N σ (E)] −1 , where N is the number of molecules per unit volume and σ (E) is the sum of the cross sections for an ion of energy E. A type of interaction is then sampled according to its energy-dependent probability, which is the ratio of the interaction cross section to the total cross section at this energy. The sampling algorithms described above are used to calculate the energy loss of the ion and of the energy of the secondary electrons. The great majority of electron ejected during the interaction of an heavy charged particle with matter have low energy (<100 eV) [5] and can be followed subsequently by using a Monte Carlo transport code like the one described in Cobut et al [13]. However, the maximum energy transfer (equation (12)) for a 1 GeV amu −1 ion is 3.3 MeV, and it increases to 136 MeV for a 10 GeV amu −1 ion. Although relatively infrequent, the transport of these high energy electrons needs to be described adequately because their penetration range defines the geometry of the track [17]. This will be the subject of a forthcoming paper.
The stopping power (LET) is calculated by simulating several histories of short track segments (10 µm) and keeping the energy constant along the track, but otherwise using the same Monte Carlo transport methods described above. The calculations of stopping power (figure 5) for 1 H + , 4 He 2+ , 12 C 6+ , 20 Ne 10+ and 56 Fe 26+ are in good agreement with the calculations derived from other codes [16], ICRU values [18] and CRC Handbook values [19], for all ions and the energy range considered (∼10 −1 -10 4 MeV amu −1 ). It should be noted that other corrections probably have to be added at ultrarelativistic energies (>10 4 MeV amu −1 ), but they are not considered in this model.
In figure 7, the radial distribution profile of the primary radiolytic species (H • , • OH, H 2 , H 2 O 2 and e − aq ) at ∼10 −12 s of these tracks is shown. These distributions are consistent with similar work performed with 70 keV µm −1 ions [3], which has demonstrated that the extent of penumbra increases for higher Z ions at the same LET.
Another important quantity is the radial dose distribution, which is of interest for microdosimetric purposes and which is also useful to validate the current work with previous calculations. It is calculated by recording the energy deposited and the radial distance from the Y -axis of the energy loss events. The energy deposited in each differential volume element is then divided by its volume and converted to Gy. The calculated radial dose is shown in figure 8 for (a) 1 H + (1 MeV amu −1 , LET ∼ 33 keV µm −1 ), (b) 20 Ne 10+ (377 MeV amu −1 , LET ∼ 31 keV µm −1 ) and (c) 4 He 2+ (0.55 MeV amu −1 ) and 56 Fe 26+ (1 GeV amu −1 ) (LET ∼ 150 keV µm −1 for both ions). In all cases, the results from our calculations are in good agreement with previous amorphous tracks calculations [20] and experimental data [21].
From figures 6-8, the extent of the penumbra increases with increasing charge and energy per nucleon for the same LET value, because the maximum energy transferred from the ion to an electron in an ionization increases with the energy per nucleon.
Monte Carlo simulation of liquid water radiolysis has greatly contributed to our understanding of ionizing radiation by giving detailed information on the energy deposition and spatial distribution of the radiolytic species as a function of different parameters such as the ion type and the LET of the HZE particle. This work shows that ions of the same LET can have completely different radial dose and radiolytic species distribution profiles, and this should provide further insight into the biological effects of HZE nuclei by discriminating the role of the so-called core and penumbra of the tracks. These studies are of particular significance for improving the understanding of radiobiology experiments with heavy ions, radiation protection issues and quantification of risks associated with the effects of GCR to which astronauts are exposed during prolonged space exploration missions.