Cosmic-Ray Tracks in Astrophysical Ices: Modeling with the Geant4-DNA Monte Carlo Toolkit

During the bombardment of some target material (such as a dust-grain ice mantle) by an energetic primary ion (such as a cosmic ray), the primary ion will collide with the atoms and molecules that comprise the material. Following Bohr (1913), it is customary to divide these collisions into two main categories, namely, inelastic (electronic) and elastic (nuclear) components. For a particle of some energy, E, moving through some material of mass density, ρ, one can thus describe the energy lost per unit path length (${dE}/{dx}$) as

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dx}}=\rho ({S}_{{\rm{e}}}+{S}_{{\rm{n}}}),\end{eqnarray} \tag{ 1 }$$

where S_e and S_n are, respectively, the electronic and nuclear loss functions, with units of cm² eV g⁻¹. The nuclear component S_n is substantially smaller than S_e at energies relevant to cosmic rays and is implicated in changes to the physical structure of the target through, e.g., sputtering and the formation of lattice defect sites, more so than changes to the composition of the target through chemical reactions. Thus, we ignore the elastic component of the primary ion energy loss in this work and instead focus on the contribution of the inelastic component (Sigmund 1969; Johnson 1991).

In general, inelastic collisions between the primary ion and atoms in the material result in the ionization and excitation of target species (Spinks & Woods 1990; Shingledecker & Herbst 2018). Ionizing collisions result in the formation of secondary electrons, which have a broad energy spectrum, but the average energies of which do not exceed around 50–70 eV, depending on the target material (Spinks & Woods 1990). Along the trajectories of these secondary electrons, further ionizations and excitations occur. Collectively, the trajectories of the primary ion and all secondary electrons in the target are referred to as the track. Most secondary electrons are stopped near the site of their formation, and thus, the track can approximately be pictured as a cylinder characterized by some radius, r_cyl, which we will refer to as the track "core."

Within this cylindrical region surrounding the path of the primary ion, the short-lived excited (suprathermal) species drive a rich variety of reactions at even very low temperatures ($\lt 10$ K) that can result in the formation of COMs and even prebiotic molecules (Holtom et al. 2005; Lafosse et al. 2006; Hudson et al. 2008). In Abplanalp et al. (2016), it was shown for the first time that reactions involving these suprathermal species are critical for reproducing the chemistry of cold cores, with later investigations showing that their inclusion in astrochemical models results in significant enhancements of the abundance of COMs such as methyl formate under TMC-1 conditions (Shingledecker et al. 2018).

The energy deposited in the track core also results in a sudden, sharp increase in the temperature of this region (Leger et al. 1985; Bringa & Johnson 2004; Ivlev et al. 2015). This rise in temperature further stimulates chemical changes in the target and, in particular, drives reactions with energy barriers that otherwise could not occur at the equilibrium temperature of the ice mantle. Thus, taken together, the combination of suprathermal reactions and thermal chemistry in the hot track core represent two promising mechanisms that can help explain the observations of COMs in cold prestellar cores.

At the interface between the track core and the surrounding vacuum, i.e., the point at which the primary ion enters the target, a hot spot forms in the material with an area of $\sim \pi {r}_{\mathrm{cyl}}^{2}$. In this process, known as impulsive spot heating, the increased temperature of the ice surface within this area significantly increases the rate of thermal desorption. As shown by Leger et al. (1985), assuming ${r}_{\mathrm{cyl}}\lt {r}_{\mathrm{grain}}$, where r_grain is the radius of the grain, this process is independent of the actual grain size.

Conversely, the subsequent process of whole grain heating is not independent of grain size and occurs as a result of the distribution of the heat deposited in the core throughout the rest of the ice mantle and underlying dust grain. The timescale of this heating has an a² dependence on the grain radius, a, such that for interstellar grains with average radii of a ≈10⁻⁵ cm, it occurs on the order of nanoseconds (Leger et al. 1985). Fast, exothermic radical–radical recombinations triggered by such heating could result in the catastrophic loss of the ice mantle through a so-called grain explosion, first noted by Greenberg & Yencha (1973).

A separate mechanism that could likewise trigger grain explosions is impulsive spot heating. This process was studied in detail by Ivlev et al. (2015), who showed that, depending on (a) the value of r_cyl, and thus, the volume of the core, as well as (b) the amount of energy deposited therein, a dramatic bifurcation in the fate of ice mantles can occur, characterized by the value of a dimensionless parameter referred to as λ by Ivlev and coworkers. For λ greater than some critical value ${\lambda }_{\mathrm{CR}}=9.94$, the bombardment by an energetic ion will similarly result in the sudden loss of the ice mantle via grain explosion.

From the preceding discussion, it should hopefully be clear that knowledge of r_cyl is required for accurate considerations of both the chemistry and interfacial dynamics, which occur as a result of the bombardment of a dust-grain ice mantle by a cosmic ray. To date, previous astrochemical works dealing with these topics have typically relied on the estimation from Leger et al. (1985) of ${r}_{\mathrm{cyl}}\approx 5\,\mathrm{nm}$, independent of primary ion velocity (Shen et al. 2004; Ivlev et al. 2015; Kalvāns 2016). Conversely, Bringa & Johnson (2000) proposed that r_cyl should increase with ${dE}/{dx}$ due to fast ($\lt 1$ ps) transport of excitation energy over a few lattice spacings. Based on fits to experimental data, they proposed an expression for r_cyl that has a linear dependence on ${dE}/{dx}$ and an overall density dependence of ${\rho }^{1/3}$ (Bringa & Johnson 2000).

In principle, however, the radius of the track core is qualitatively described by the stopping range of an electron with the average energy of the ejected secondary electrons. Assuming the continuous slowing down approximation, where all particles with the same energy are assumed to travel the same average distance (Johnson 1990), this range for secondary electrons with an average energy W_av is given by

$$\begin{eqnarray}&&{r}_{\mathrm{cyl}}\approx {\int }_{0}^{{W}_{\mathrm{av}}}{\left(\displaystyle \frac{{dE}}{{dx}}\right)}_{\mathrm{electron}}^{-1}{dE}=\displaystyle \frac{1}{\rho }{\int }_{0}^{{W}_{\mathrm{av}}}\left(\displaystyle \frac{1}{{S}_{{\rm{e}},\mathrm{electron}}}\right){dE}.\end{eqnarray} \tag{ 2 }$$

Here, ${\left({dE}/{dx}\right)}_{\mathrm{electron}}$ is the energy deposited per unit path length for electrons, and similarly, ${S}_{{\rm{e}},\mathrm{electron}}$ is the electronic loss function, also for electrons. In principle, one could use Equation (2) to estimate the track-core radius; however, such an approach would require accurate analytical expressions for the electronic stopping losses, the derivation of which are beyond the scope of the current work. Nevertheless, Equation (2) is still useful, as it allows us to qualitatively understand the results obtained using the Monte Carlo methods utilized here, for example, by comparing the dependence of our results on material density with the $1/\rho$ dependence one would expect from Equation (2).

Thus, we are presented with three conflicting predictions as to the value of r_cyl. In this work, we seek to resolve this confusion and establish more explicitly the value of this critical parameter based on a leading-edge Monte Carlo code, designed to yield astrochemically relevant values for amorphous solid water over a range of cosmic-ray proton energies. The rest of this work is arranged as follows: in Section 2, we provide details of our model and computational approach; in Section 3, we present the results of calculations and discuss their astrophysical significance; and finally, in Section 4, we summarize our conclusions.

For this work, we have employed the Geant4 v10.6 Monte Carlo simulation toolkit (Agostinelli et al. 2003; Apostolakis et al. 2009; Allison et al. 2016), which was initially designed to simulate systems relevant for high-energy physics. However, the flexibility of the code allows for its application to problems in a wide variety of fields, including medical physics and astrophysics. The toolkit was later extended by the Geant4-DNA project to simulate microdosimetry through the addition of additional physical processes, such as excitation and charge exchange, that were not included in the original Geant4 code and allow for the accurate modeling of collisional events down to energies of a few eV (Incerti et al. 2010a, 2010b; Bernal et al. 2015; Incerti et al. 2018). The original motivation for this extension, as indicated by the addition of "DNA" to the name of the toolkit, was to investigate the effects of ionizing radiation on biological systems, including especially DNA and RNA damage.

For this study, we utilized the G4EmDNAPhysics_option2 physics list. A full description of the processes and valid energy ranges for particles considered in Geant4-DNA can be found in Incerti et al. (2018) and includes elastic electron scattering, shell ionization cross sections (five shells), excitation cross sections (five levels), full secondary electron cascade generation from individual shells (using shell-specific differential ionization cross sections), vibrational excitation, and molecular attachment. In the context of Geant4-DNA, and indeed, of all similar MC codes, the track is defined as the collection of the above-mentioned interaction "points," which occur at a given set of x, y, and z coordinates in our simulated volume. Among other restrictions of quantum origin, the spatial extent of the "point" cannot be smaller than the dimensions of the target molecule which, in the case of water, is governed by the ∼0.3 nm diameter of the molecule.

In the interstellar medium, water ice forms on dust grains through the adsorption and subsequent reaction of, e.g., O, H, OH, and O₂ (Cuppen et al. 2010; Ioppolo et al. 2010). This ice exists mostly as amorphous solid water (ASW). The properties of this glassy metastable material depend on the physical conditions under which it formed, combined with the effects resulting from any changes of these conditions and of any subsequent processing. The two main types of ASW of astrophysical relevance are low-density amorphous (LDA) and high-density amorphous (HDA) ices, which have densities of 0.94 and 1.1 g cm⁻³, respectively (Narten et al. 1976; Jenniskens & Blake 1994). In the ISM, the bombardment of LDA by energetic particles similar to cosmic rays has been found to result in its compactification, leading possibly to HDA (Palumbo 2005; Mitterdorfer et al. 2014).

By default, the Geant4-DNA simulations include data only for liquid water; however, given the structural similarity of it with the glassy ASW, we here approximate ASW by scaling the density of the material in our model. Because the mean free path of electrons is dependent on the density of the material, we here perform calculations at densities relevant to both LDA and HDA. In the context of MC transport of energetic charged particles, several studies have examined the energy-loss function properties of solid water (amorphous and hexagonal ice) and found them to be very similar to those of liquid water (Emfietzoglou et al. 2007a, 2007b, 2008). Therefore, we expect the energy loss of charged particles in ASW and liquid water to be quite similar, and as such, the inclusion of ASW-specific interaction cross sections will not appreciably change the results or conclusions described in the following sections. Perhaps the main uncertainty in the present work comes from the lack of rigorous corrections to the first Born approximation for inelastic electron scattering below about 100 eV. Geant4-DNA has already implemented such corrections in some of its physics models, including the one used in this work, but they are mostly phenomenological. It is possible that these uncertainties may influence very low-energy electron transport at the few nanometer scale; however, an investigation in this matter is beyond the scope of our study.

In our simulations, we represent the ice as a cube with edges 1 μm in length. This ice is then bombarded with protons—the major constituent of cosmic rays—with energies between 100 keV and 100 MeV, which cover both the peak and subsequent falloff of the electronic stopping power, as depicted in Figure 1. Incident primary ions are assumed to collide with the ice normal to the surface in the center of the topmost side. Because the stopping length of protons in the energy range considered here is larger than the 1 μm thickness of our ice, they are able to pass completely through, at which point they are considered to have left the system and are not followed any further. For each model run, the simulation begins with the first collision of the primary ion and ends when all secondary electrons reach energies of around 7.4 eV, below which Geant4-DNA does not currently simulate them.

Zoom In Zoom Out Reset image size

Finally, in order to aid the calculation of track-core widths, we have disabled elastic scattering processes for the incident protons⁸ as well as neutral hydrogen,⁹ because the incident H⁺ can capture an electron from the target material. As one can see from Figure 1, because the nuclear component of the stopping power is more than ∼2 orders of magnitude less than the electronic component, and moreover, because it does not affect the width of the track core, this assumption should not hinder the accuracy of the resulting calculated values of r_cyl.

The results of our simulations for LDA ice are listed in Table 1. There, in addition to listing the total radius counting all collisions along the track (${r}_{\mathrm{tot}}^{\mathrm{LDA}}$), we further list radii in which only ionization or excitation collisions of the total track were counted. These radii are plotted as a function of initial particle energy in the left panel of Figure 3. In addition, in Figure 6 of Appendix A, we show the 1D and 2D representations of the total track, as well as several subtracks considered here. To better illustrate our method of obtaining the values given in Table 1, we show in Figure 4 a representative example of the cumulative distribution plots for a 0.1 MeV H⁺, as well as contours indicating the track-core radii.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From Figure 3, one can see that the calculated radii for ionization, r_ion, are the smallest among those shown. This result is expected as secondary electrons only have energy to undergo ∼a few ionization collisions before falling below the ionization threshold of the material. Moreover, ionization cross sections are generally on the order of ∼1 order of magnitude larger than excitation cross sections (Johnson 1991), and therefore, the former are expected to occur with a similarly higher probability than the latter (Shingledecker et al. 2017). For secondary energies above the ionization threshold of the material, ionizing collisions will dominate and, thus, will occur within a shorter distance to the point of formation than for excitation collisions.

Because we follow the average energy deposited per primary ion rather than the total energy deposition, the fact that more energy is lost in ionizing collisions (being equal to the ionization energy of the material) than excitation collisions ($\lt 10$ eV) means that r_energy closely follows r_ion, though at a slightly larger value due to the contribution of energy deposited during excitation collisions.

One can also see that the total radius of the track core, r_tot, is nearly identical to the radius of excitation collisions, r_exc, from 0.1 to 100 MeV. Over this energy range, these radii display values of $\sim 9.0\pm 1.0\,\mathrm{nm}$. The close association between the total and excitation radii is again not surprising, because, once falling below the ionization threshold, excitation collisions of various kinds, e.g., electronic and vibrational, are the dominant inelastic processes until the electron falls below the excitation energy threshold of the material.

These radii characterize the behavior of secondary electrons in the material and are thus sensitive to their energies upon formation. As noted by Rudd et al. (1992), for an incident proton with energy E_ion, momentum p, and mass m_p, the maximum secondary electron energy, W_max, which can be produced by that ion is given by

$$\begin{eqnarray}&&{W}_{\max }(p)=2{m}_{e}{v}_{0}^{2}=4\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{{\rm{p}}}}\right){E}_{\mathrm{ion}},\end{eqnarray} \tag{ 3 }$$

where m_e and v₀ are, respectively, the electron mass and incident ion velocity. Thus, the maximum secondary electron energy capable of being produced by a 0.1 MeV proton is 218 eV, while for a 100 MeV proton, the value is 218 keV. However, the probability of producing such a secondary electron with energy W_max is small, with the majority of secondary electrons having a broad energy spectrum but with average energies of not more than 50–70 eV, depending on the target material. The secondary electron energies indicated by our radii do not show the linear dependence on E_ion implied by Equation (3). In fact, our results are in agreement with semiempirical estimations of secondary electron energies by Rudd (1988), who found that the average secondary electron energy, W_av, increased until it saturates at around incident proton energies of a few tenths to a few MeV, again depending on the target (see Figure 10 of Rudd 1988). The results presented in Figure 10 of Rudd (1988), along with the expression for r_cyl given in Equation (2), together provide a very reasonable description of the energy dependence of the radii shown in Figure 3, where the track-core radii will likewise follow changes in W_av. This ability to qualitatively describe the findings presented in Figure 3 in light of Equation (2) and the work of Rudd (1988) is thus a reassuring indicator of the reliability of our approach.

We note that the relationship between our derived r_cyl values and incident proton energies, though in agreement with the expected average secondary electron energies calculated by Rudd (1988), is not what one would expect based on the formula derived by Bringa & Johnson (2000), which is linearly dependent on the stopping power $({dE}/{dx})$. As shown in Figure 1, the electronic stopping power varies by around three orders of magnitude over the 0.1–100 MeV range considered here; however, as shown in Figure 3, our track-core radii vary by only a factor of a few over that range.

Shown in Figure 5 are results for HDA ice. Calculated radii for HDA ice are also listed in Table 1 and plotted as a function of energy in the right panel of Figure 3. The 1D and 2D representations of the total track, as well as several subtracks considered here, are given in Figure 7 of Appendix B. As with the LDA ice, we see a similar trend where the smallest predicted radii are for ionization, with values that level off at $\sim 5.5\pm 0.5\,\mathrm{nm}$. Likewise, following roughly the same trend but at slightly larger values, the energy deposition radii have values of ∼4.2 nm. Also similar to the LDA ice, the total track-core radii closely follow the radii of excitation collisions, again with an average radius of $\sim 8.0\pm 0.5\,\mathrm{nm}$, which is smaller than the LDA value.

Equipped now with results for both LDA and HDA ices, we next turn our attention to a consideration of the dependence of the track radii on the material density. Recall from Equation (2) in Section 1 that, in principle, the track radius is determined by the stopping range of an electron with the average energy of the ejected secondary electrons, resulting in a density dependence of ${r}_{\mathrm{cyl}}\propto 1/\rho$. Using the densities we employed here for LDA and HDAs ice of 0.94 and 1.1 g cm⁻³, we predict a ratio of ${r}_{\mathrm{cyl}}^{\mathrm{HDA}}/{r}_{\mathrm{cyl}}^{\mathrm{LDA}}=0.85$. Taking the ratio of the r_tot values at 6 MeV, we similarly obtain of ${r}_{\mathrm{tot}}^{\mathrm{HDA}}/{r}_{\mathrm{tot}}^{\mathrm{LDA}}=0.85$, thereby nicely recovering the expected density dependence.

These results are somewhat at odds with the findings of Bringa & Johnson (2000), whose expression for track-core radii has an overall density dependence of ${r}_{\mathrm{cyl}}\propto {n}^{1/3}$. Similarly, the peak radii we obtain in our simulations of 9.9 nm and 8.4 nm for LDA and HDA ices, respectively, are roughly double the constant 5 nm value assumed in Leger et al. (1985) and widely used in astrophysical models. However, as with the LDA results, the track-core radii we obtain for HDA ice qualitatively follow the behavior we would expect from Equation (2) and Figure 10 of Rudd (1988).

Zoom In Zoom Out Reset image size