On Simulating the Proton-irradiation of O₂ and H₂O Ices Using Astrochemical-type Models, with Implications for Bulk Reactivity

The upcoming launch of the James Webb Space Telescope (JWST) promises to usher in a new astrochemical "ice age" by allowing for the collection of an unprecedented amount of data regarding species frozen onto dust grains in the interstellar medium (ISM). Astrochemical models will surely figure prominently in both interpreting these data as well as in informing future observing proposals. This opportunity, however, requires models that are up to the task and unfortunately significant uncertainties remain regarding the chemistry of dust grains and dust grain ice mantles (Cuppen et al. 2017)—a situation that has given it the dubious reputation of "the last refuge of the scoundrel" (Charnley et al. 1992).

Currently, it is a common practice in modern astrochemical codes to distinguish between species in the top few monolayers of dust grain ice mantles (the selvedge), and those underneath, in the interior (the bulk), though there exists considerable disagreement as to whether, or to what degree, such bulk species are chemically active (Cuppen et al. 2017). The hypothesis that chemistry within the ice mantle is less significant than that in the selvedge rests on the following three main lines of reasoning:

• (a)  
  species in the bulk are more tightly bound;
• (b)  
  a richer chemistry is possible in the selvedge, given the constant exchange of species with the surrounding gas; and
• (c)  
  photoprocesses, which further drive chemical reactions, are more efficient in the topmost monolayers.

An implicit assumption behind (a) is that reactions in the bulk occur largely via thermal diffusion, as on surfaces, and thus reaction rates will be low or negligible for species with high binding energies in cold environments. However, some previous experimental studies strongly suggest that reactions in the bulk can occur quickly among neighboring species without the need for diffusion (Baragiola et al. 1999; Bennett & Kaiser 2005; Abplanalp et al. 2016; Ghesquière et al. 2018).

Similarly, just as selvedge chemistry is stimulated via photoprocesses and interactions with the gas—as noted in (b) and (c)—so too have experiments shown that bombardment by energetic ions, which penetrate solids much more efficiently than photons (Spinks & Woods 1990; Gerakines et al. 2001), causes physicochemical changes that drive a rich chemistry within the ice at low temperatures. This radiation chemistry is quite relevant to the ISM, since ices there are continually bombarded by energetic particles of one kind or another, though typically in the form of cosmic rays (Indriolo & McCall 2013; Rothard et al. 2017), stellar winds (Hudson & Moore 2001; Madey et al. 2002), or radionuclide emission (Cleeves et al. 2013).

Thus, the processing of dust grain ice mantles by ionizing radiation is a very real phenomenon and the fact that such interactions can lead to the formation of astrochemically interesting, even prebiotic, species makes radiation chemistry very promising from a modeling standpoint (Kaiser & Roessler 1997; Holtom et al. 2005; Hudson et al. 2008; Abplanalp et al. 2016). Until recently, though, the variety and complexity of the underlying microscopic processes have hindered attempts to add radiation-chemical reactions to astrochemical codes (Abplanalp et al. 2016; Shingledecker et al. 2018).

These processes result from collisions between an incoming particle (the primary ion) and species in the target material. Following the formalism of Bohr (1913), it is customary to divide such collisions into two categories, namely those in which energy is transferred to target nuclei (elastic) or electrons (inelastic). One can therefore express the energy loss per unit path length of the primary ion—called the stopping power—as follows:

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

where n is the density of the target (cm⁻³) and S_n and S_e are, respectively, the so-called nuclear/elastic and electronic/inelastic stopping cross sections (eV-cm²). Inelastic collisions, in turn, can further be divided into those in which target species are ionized or excited. Such ionizations result in the formation of so-called "secondary electrons," which typically have energies under 50 eV and can themselves ionize or excite target species, thereby further propagating the physicochemical effects initiated by the primary ion (Johnson 1990; Spinks & Woods 1990).

Recently, in Shingledecker & Herbst (2018, hereafter, SH), we introduced a method for treating radiation-chemical processes that is simple enough to include in existing astrochemical models, while simultaneously preserving salient features inferred from experimental data. This method was later added to the astrochemical-type Nautilus v1.1 code (Ruaud et al. 2016) and used by us to determine what effects, if any, these new processes would have on cold core models (Shingledecker et al. 2018). Preliminary results from that study indicated that cosmic-ray-induced reactions could indeed result in enhanced abundances for a number of species, including methyl formate and ethanol, and allowed us to include novel reactions such as insertions (Bergner et al. 2017; Bergantini et al. 2018). However, given the uncertainties sometimes involved in comparing the results of astrochemical models with observations of interstellar environments—particularly in cases where grain-surface chemistry is involved—questions lingered regarding how well astrochemical models were actually simulating irradiation chemistry.

To that end, we have carried out simulations of two different ice irradiation experiments using an existing astrochemical-type model (Vasyunin et al. 2017), modified to include the SH method for treating radiation chemistry. Specifically, we have modeled the bombardment of a 5 K pure O₂ ice by 100 keV protons as described in Baragiola et al. (1999), as well as that of a pure water ice at both 16 and 77 K by 200 keV protons, as reported in Gomis et al. (2004a). Doing so has afforded us the opportunity not only to compare how well our simulations replicate the experimental data, but also to directly gauge how well our models perform more generally with handling reactions in bulk ice over a range of temperatures. By thus testing our approaches to ice chemistry in well-constrained experimental systems, confidence can be increased in the results of models of much less well-constrained interstellar environments.

The paper is organized as follows. In Section 2 we discuss details of both the model and chemical network used in this work. In Section 3 our major results are described and their significance is discussed. Finally, our conclusions and suggestions for future work are presented in Section 4.

The basis of this work is the MONACO code, described in Vasyunin et al. (2017), which uses a rate-equation approach for reactions in the gas phase and the so-called "modified rate-equation" method for reactions on and in dust grain ice mantles (Garrod 2008; Garrod et al. 2009). As in Vasyunin & Herbst (2013) we assume the selvedge, a term for the altered layers near the top of the mantle, is comprised of the top four monolayers of the ice (Vasyunin & Herbst 2013).

We have modified the code to account for irradiation-chemical reactions on and in ice using the SH method, the foundation of which is the assumption that one of the following outcomes is possible upon collision between an energetic particle, the primary ion or a secondary electron, and some target species, A:

$$\begin{eqnarray}&&{\rm{A}}\rightsquigarrow {{\rm{A}}}^{+}+{{\rm{e}}}^{-},\end{eqnarray} \tag{ P1 }$$

$$\begin{eqnarray}&&{\rm{A}}\rightsquigarrow {{\rm{A}}}^{+}+{{\rm{e}}}^{-}\to {{\rm{A}}}^{* }\to {{\rm{B}}}^{* }+{{\rm{C}}}^{* },\end{eqnarray} \tag{ P2 }$$

$$\begin{eqnarray}&&{\rm{A}}\rightsquigarrow {{\rm{A}}}^{* }\to {\rm{B}}+{\rm{C}},\end{eqnarray} \tag{ P3 }$$

$$\begin{eqnarray}&&{\rm{A}}\rightsquigarrow {{\rm{A}}}^{* },\end{eqnarray} \tag{ P4 }$$

where the curly arrow indicates a collision. In processes (P1) and (P2), A is ionized, though with fast charge recombination leading to the formation of suprathermal dissociation products, designated with an asterisk, assumed to occur in (P2). Processes (P3) and (P4) represent collisions in which a suprathermal product is electronically excited directly into some higher state, which in the case of (P3) is assumed to lead to the formation of thermal products. The efficiencies of (P1)–(P4) are characterized by radiochemical yields, called G-values, which give the number of molecules created (or destroyed) per 100 eV transferred from bombarding particles to the system (Dewhurst et al. 1952; Spinks & Woods 1990). A full list of the radiolysis processes used here is given in Table 3 of Appendix A. The yields given there were arrived at through a two-step process where initial values were estimated using the SH method and subsequently adjusted to improve agreement with the experimental data. It should be noted that following Shingledecker et al. (2018), we assume the rate of process (P1) is zero, i.e., that all charged species produced in the bulk via ionization are quickly neutralized.

First-order rate coefficients for each radiolysis process are calculated using

$$\begin{eqnarray}&&{k}_{\mathrm{rad}}=\left(\displaystyle \frac{G}{100\,\mathrm{eV}}\right){S}_{{\rm{e}}}\phi ,\end{eqnarray} \tag{ 2 }$$

with ϕ being the radiation flux and S_e being the previously described inelastic/electronic stopping cross section. In this work, we assume the flux to be 1.0 × 10¹¹ H⁺ cm⁻² s⁻¹, and use the stopping cross sections given in Table 1.

We assume that, in the bulk, suprathermal species produced via radiolysis rapidly undergo either (a) reaction with a neighbor, or (b) are quenched by the surrounding material. Rate coefficients for the reaction between A and B—where either one could be suprathermal—are calculated using the expression

$$\begin{eqnarray}&&{k}_{\mathrm{fast}}={f}_{\mathrm{br}}\left[\displaystyle \frac{{\nu }_{0}^{{\rm{A}}}+{\nu }_{0}^{{\rm{B}}}}{{N}_{\mathrm{bulk}}}\right]\exp \left(-\displaystyle \frac{{E}_{\mathrm{act}}^{\mathrm{AB}}}{{T}_{\mathrm{ice}}}\right),\end{eqnarray} \tag{ 3 }$$

where f_br is the branching fraction, N_bulk is the total number of bulk species in the simulated ice, and ${E}_{\mathrm{act}}^{\mathrm{AB}}$ is the activation energy for reaction. Absent other data, we assume here that suprathermal species react barrierlessly. The characteristic vibrational frequency, ${\nu }_{0}^{{\rm{A}}}$, is given by

$$\begin{eqnarray}&&{\nu }_{0}^{{\rm{A}}}=\sqrt{\displaystyle \frac{2{n}_{{\rm{s}}}{E}_{{\rm{D}}}^{{\rm{A}}}}{{\pi }^{2}{m}_{{\rm{A}}}}},\end{eqnarray} \tag{ 4 }$$

with ${E}_{{\rm{D}}}^{{\rm{A}}}$ and m_A being the desorption energy and mass of A, respectively, and n_s is the phsysisorption site density—here 6.55 × 10¹⁴ cm⁻² (Vasyunin et al. 2017). This characteristic vibrational frequency is also used as the first-order rate coefficient for the quenching of the suprathermal species by the ice. It should be noted that, particularly in water ice, these electronic excitations, also called excitons, may diffuse from the interior of the ice to the selvedge, where they can drive the desorption of species into the gas (Thrower et al. 2011; Marchione et al. 2016). Thus, our assumption that these suprathermal species are rapidly quenched should be seen as a first-order approximation of their real behavior.

As an example of how the above equations are utilized in our code, consider the time-dependent abundance of B*. Assuming B* is produced via the radiolysis of A as in (P2) and destroyed both via quenching and reaction with some other bulk species, X, then the rate of change in the number of B* as a function of time can be described by

$$\begin{eqnarray}&&\displaystyle \frac{d\,{N}_{{{\rm{B}}}^{* }}}{{dt}}={k}_{\mathrm{rad}}^{{\rm{P}}2}{N}_{{\rm{A}}}-{\nu }_{0}^{{\rm{B}}}{N}_{{{\rm{B}}}^{* }}-{k}_{\mathrm{fast}}{N}_{X}{N}_{{{\rm{B}}}^{* }},\end{eqnarray} \tag{ 5 }$$

where the reaction between B* and X is, as with all suprathermal species here, assumed to occur barrierlessly.

In order to gauge how well astrochemical models are able to reproduce bulk chemistry, we have run the following three sets of simulations for each experiment considered here using the parameters listed in Table 1:

• (a)  
  a set of models in which thermal radicals and atomic oxygen were assumed to react non-diffusively with neighboring species in the ice with rate coefficients calculated using Equation (3);
• (b)  
  a set of models in which all bulk rate coefficients were calculated using the standard diffusive formula (Hasegawa et al. 1992) along with diffusion barriers obtained from the desorption energies listed in Table 2 and no tunneling under diffusion barriers for any species;
• (c)  
  a set of models similar to (b) but with tunneling through diffusion barriers allowed for H, H₂, and O.

Table 2.  Desorption Energies for Relevant Species

Species	ED (K)	Source
O	1660a	He et al. (2015)
	1040b	Minissale et al. (2013)
O2	930	Jing et al. (2012)
O3	1833	Jing et al. (2012)
H	450	Garrod & Herbst (2006)
(Thrower et al. 2011; Marchione et al. 2016) H2	430	Garrod & Herbst (2006)
OH	2850	Garrod & Herbst (2006)
H2O	5700	Garrod & Herbst (2006)
HO2	4510	${E}_{{\rm{D}}}^{{\rm{O}}}+{E}_{{\rm{D}}}^{\mathrm{OH}}$
H2O2	5700	${E}_{{\rm{D}}}^{\mathrm{OH}}+{E}_{{\rm{D}}}^{\mathrm{OH}}$

Notes. Diffusion barriers in the selvedge were assumed to be 0.5 × E_D, while those in the bulk were assumed to be 0.7 × E_D.

^aUsed in models (a) and (b). ^bUsed in model (c).

Download table as:  ASCIITypeset image

For the third set of simulations, we treat tunneling under diffusion barriers by H and H₂ using the standard formalism in Hasegawa et al. (1992) with a barrier width of 1.0 Å. For atomic oxygen, following results from Minissale et al. (2013), we use a desorption energy of 1040 K and a barrier width of 0.7 Å. In all simulations, H and H₂ were assumed to tunnel through reaction barriers (Tielens & Hagen 1982). The competition mechanism for systems with chemical activation energies (Chang et al. 2007; Herbst & Millar 2008) is used here in all models only for reactions on the surface.

For the chemical network, we used the reactions listed in Table 4 of Appendix B, as well as the neutral–neutral oxygen reactions listed in Table 7 of Shingledecker et al. (2017). In cases where we have been unable to find values for activation energies, we have assumed a barrier of 10⁴ K, except for radical–radical reactions, which were assumed to be barrierless. For every thermal reaction of the type given in Table 4, we included variants involving one suprathermal reactant. To illustrate this point, for every reaction of the form ${\rm{A}}\ +\ {\rm{B}}\ \to \,\mathrm{products}$, we include the following suprathermal variants:

$$\begin{eqnarray}&&{{\rm{A}}}^{* }\ +\ {\rm{B}}\ \to \,\mathrm{products},\end{eqnarray} \tag{ R1 }$$

$$\begin{eqnarray}&&{\rm{A}}\ +\ {{\rm{B}}}^{* }\ \to \,\mathrm{products}.\end{eqnarray} \tag{ R2 }$$

A number of reactions listed in Table 4 have barriers that are quite high for astrochemical networks, but are more common in the combustion-chemical literature from which most are taken. Even though the rates of such reactions involving thermal species will indeed be negligible in our models, the variants involving suprathermal species should be efficient even at the temperatures considered here. For example, even a 40,000 K barrier corresponds to only ∼3.4 eV, a realistic energy for an electronically excited suprathermal species.

In order to examine how well astrochemical models using the SH method are capable of reproducing experimental data, as well as to determine the accuracy of bulk chemistry in such codes, we have simulated the ice-irradiation systems described in two previous radiation-chemical studies discussed above; the results of these simulations are described below.

3.1. H⁺ Bombardment of Pure O₂ Ice

In this system, the only radical used is atomic oxygen, which can be thermal or suprathermal. We first simulated the irradiation of a pure O₂ ice at 5 K by 100 keV H⁺, following the study by Baragiola et al. (1999), which was previously modeled using the much more detailed Monte Carlo code, CIRIS (Shingledecker et al. 2017). An irradiated O₂ ice is, in many respects, an ideal system for testing simulations of radiation chemistry, given the limited number of possible neutral species, i.e., O, O₂, and O₃. However, in the Monte Carlo code even this seemingly simple system required a network of ∼50 reactions or processes and involved a level of detail beyond the practical capabilities of rate equation-based astrochemical codes, such as explicitly calculating the tracks of incident ions and secondary electrons or following each ice species on a hop-by-hop basis.

Despite these limitations though, as shown in Figure 1(a), even our simplified approach is capable of reasonably reproducing the concentration of O₃ as a function of fluence. Here, fluence—the product of the H⁺ flux and irradiation time—represents the number of ions that have bombarded the ice per cm². As can be seen in Figure 1, the best results were obtained in model (a), where thermal atomic oxygen is assumed to react non-diffusively in the bulk, yielding excellent agreement between calculated and experimental abundances in the steady stage regime. An analogous assumption was made in our previous Monte Carlo simulations using CIRIS, with similarly good results (Shingledecker et al. 2017). The slight underprediction of O₃ at fluences under ∼10¹³ H⁺ cm⁻¹ is likely due to additional processes—such as dissociative electron attachment (DEA)—which are also driven by energetic ion bombardment but not currently considered in this simplified treatment of radiolysis (Arumainayagam et al. 2010). Also noticeable in Figure 1(a) is the fact that atomic oxygen abundances remained low throughout the simulation, despite their continuous production from the dissociation of O₂, and to a lesser degree, O₃.

Zoom In Zoom Out Reset image size

Comparison with Figure 1(b) reveals that model (b), using standard Langmuir–Hinshelwood rate coefficients, is the least accurate. There, though the modeled O₃ abundance is only a factor of a few lower than the experimental values, an unrealistically large abundance of atomic oxygen is predicted, contrary to what has been suggested by previous studies (Gerakines et al. 1996; Baragiola et al. 1999; Bennett & Kaiser 2005). This buildup is caused by the very slow rates of thermal diffusion reactions at 5 K, given the high atomic oxygen desorption energy of 1660 K (He et al. 2015).

In order to determine the possible effects of quantum tunneling by atomic oxygen through diffusion barriers, as suggested by Minissale et al. (2013), we ran further simulations using their best-fit results of a barrier width of 0.7 Å and a desorption energy of 1040 K (assuming E_b = 0.5E_D). Tunneling for H and H₂ is treated using the method of Hasegawa et al. (1992), in which a barrier width of 1 Å is assumed. Comparison of the results of this model, shown in Figures 1(c), with (b) reveals that though agreement between the calculated and experimental ozone abundances is improved, the abundance of atomic oxygen remains too high.

The results of these three models demonstrate that the best agreement with experimental data—both the measured O₃ abundance as well as the inferred low O abundance—is obtained when atomic oxygen is assumed to react quickly in the bulk. We note, however, that the predicted atomic oxygen abundances in model (a) may in fact be somewhat too low, based on the observation by Bennett & Kaiser (2005) of an IR feature associated with an [O₃...O] complex, though no estimate of the atomic oxygen abundance (or an upper limit) was derived. In such irradiated ices, it is likely that some fraction of the oxygen atoms cannot react quickly, either due to trapping by O₃ or steric effects, in which case their abundance would indeed be greater than that predicted in our models. These effects could be accounted for in the method used here by decreasing the value of the pre-exponential frequency, ν₀, thereby increasing the abundance of O. Here, we use the characteristic vibrational frequency given in Equation (4) (which typically has a value on the order of 10¹² s⁻¹) as a rough approximation, though in fact the value of this parameter can be as low as 10⁻³ s⁻¹ (Theulé et al. 2013), and rate coefficients obtained using (4) should probably be interpreted as upper limits.

3.2. H⁺ Bombardment of Pure H₂O Ice

We next simulated the irradiation of a pure H₂O ice by 200 keV H⁺ at both 16 and 77 K using the three different bulk chemistry schemes described in Section 2. Here, in addition to atomic oxygen, we further assume that all radicals in our network—H, OH, and HO₂—react quickly in model (a). This system, though substantially more complex than O₂, is of greater astronomical interest, given the ubiquity of water ice in planetary bodies (Carlson et al. 2009; Altwegg et al. 2015) as well as interstellar dust grain ice mantles (Gibb et al. 2004).

Shown in Figure 2 are both our calculated H₂O₂ abundances vs proton fluence as well as the approximate steady-state experimental values from Gomis et al. (2004a) relative to H₂O of ∼1.0% at 16 K and 0.25% at 77 K. As with the pure O₂ ice, model (a) once again provides the best agreement with experimental data and notably is the only one of the three to predict a drop in hydrogen peroxide abundance at higher temperatures. Moore & Hudson (2000), who first detected this trend, speculated that it was due to the reaction

$$\begin{eqnarray}&&\mathrm{OH}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+{\mathrm{HO}}_{2}.\end{eqnarray} \tag{ R3 }$$

Zoom In Zoom Out Reset image size

We find that, in agreement with Moore & Hudson (2000), reaction (R3) is indeed behind the drop in H₂O₂ abundance at ∼80 K, because even though it is assumed that radicals such as OH react quickly, (R3) only becomes competitive at 77 K due to the 755.3 K barrier (Ginovska et al. 2007). In models (b) and (c) the effect of reaction (R3) is reduced by the slow diffusion rates of OH and H₂O₂. As can be seen in Figure 3, without the assumption that radicals react quickly in the bulk, OH in particular becomes quite abundant in the ice, even at 77 K.

Zoom In Zoom Out Reset image size

The increase in hydrogen peroxide abundance at 77 K in models (b) and (c) is due to the following reactions:

$$\begin{eqnarray}&&\mathrm{OH}+\mathrm{OH}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2},\end{eqnarray} \tag{ R4 }$$

$$\begin{eqnarray}&&{\rm{H}}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2},\end{eqnarray} \tag{ R5 }$$

where here, the higher rate of H₂O₂ formation is due to the increased mobility of the reactants. Enabling quantum tunneling through diffusion barriers, as in model (c), further speeds up the rate of reaction (R5), as well as the formation of the precursor species, HO₂, via

$$\begin{eqnarray}&&{\rm{H}}+{{\rm{O}}}_{2}\longrightarrow {\mathrm{HO}}_{2},\end{eqnarray} \tag{ R6 }$$

thereby contributing to the even higher H₂O₂ abundance at 77 K in model (c) compared with model (b).

We can gain further insights into how closely our simulations are replicating the experiment by examining G(H₂O₂), the radiolytic hydrogen peroxide formation yield. Unfortunately, we cannot simply compare the G-values given in Table 3 with experimentally determined ones directly, since our values are more representative of the immediate creation (or destruction) of target species, i.e., the efficiencies of each of the microscopic radiolytic processes given in (P1)–(P4), than the single effective experimental value, which is sensitive to the temperature-dependent chemistry of the system (Spinks & Woods 1990). However, following the method used in Moore & Hudson (2000), we can estimate what the experimental G-value might be from the slope of a linear fit to the abundance curves in Figure 4 over the pre-steady-state regime—corresponding to doses of ca. 0–10 eV, where dose is the product of the fluence and S_e. From this, we calculate the yield of H₂O₂ at 16 K to be 0.1 molecules/100 eV—exactly what was mentioned by Moore & Hudson (2000) as the yield in pure H₂O.

Zoom In Zoom Out Reset image size

Shown in Figure 5 are the abundances of H₂, O₂, and O₃ versus fluence. In model (a), the abundance of O₂ at 16 K is kept low because of destruction via (6) to form HO₂ but increases at 77 K due, in part, to more efficient formation via

$$\begin{eqnarray}&&{\mathrm{HO}}_{2}+{{\rm{O}}}_{3}\longrightarrow {{\rm{O}}}_{2}+{{\rm{O}}}_{2}+\mathrm{OH},\end{eqnarray} \tag{ R7 }$$

which has a small barrier of 490 K (Burkholder et al. 2014). Similarly, the increase in molecular oxygen abundance at 77 K in model (c) is further driven by

$$\begin{eqnarray}&&\mathrm{OH}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+{{\rm{O}}}_{2},\end{eqnarray} \tag{ R8 }$$

where the abundances of OH and HO₂ are enhanced relative to model (b) because of the effects of* quantum tunneling through diffusion barriers by H, H₂, and O—as shown in Figure 3. The decreased abundance of these radicals at 77 K in model (b), combined with their destruction with atomic oxygen, leads to the drop in [O₂] in Figure 5(b).

Zoom In Zoom Out Reset image size

Unfortunately, comparison of our O₂ results with experimental data, as with H₂O₂, is complicated by the fact that homonuclear diatomic molecules are IR-inactive. Thus, their abundances cannot be measured using standard techniques, such as Fourier Transform infrared spectroscopy. Nevertheless, it is well known that H₂ and O₂ form during water ice radiolysis based on analysis of both sputtering products as well as post-irradiation temperature-programmed desorption of the ice via mass spectrometry (Johnson & Quickenden 1997; Teolis et al. 2017). In principle, though, such measurements should be possible using Raman techniques (Rothard et al. 2017) and would be of great value, in part by enabling us to further refine both our radiochemical yields and chemical network. Interest in constraining O₂ abundances in irradiated water was recently renewed following its detection in the coma of comet 67P/C-G by Bieler et al. (2015). As can be seen in Figure 5, the maximum abundance of O₂ with respect to water achieved here is ∼0.1% in model (a) at 77 K, a value that increased only negligibly at still higher temperatures. Thus, our models predict that the radiolysis of pure H₂O ice is not the dominant mechanism behind the ∼3.8% O₂ abundances relative to water measured by Rosetta (Bieler et al. 2015). In that study, moreover, no evidence for ozone was found, though an upper limit of 1 × 10⁻⁴% relative to water was established. Interestingly, only in model (a) are the ozone abundances predicted to remain below this limit, even at 77 K.

In this work, we have simulated the bombardment of pure O₂ and H₂O ices by energetic protons using a general rate-equation-based astrochemical code, modified to include radiation-chemical processes using the SH method. These models were carried out with the MONACO program (Vasyunin et al. 2017) and a network consisting of the radiolysis processes listed in Table 3 of Appendix A and the reactions noted in Appendix B. As illustrated in Figures 1 and 2, we were able to qualitatively reproduce both the abundance of O₃ in pure O₂ (Baragiola et al. 1999) and H₂O₂ in pure H₂O (Gomis et al. 2004a) utilizing the radiochemical processes given in Table 3. Thus validated, these processes, along with the reactions given in Table 4, can be added to existing chemical networks in order to better account for physicochemical effects caused by cosmic-ray bombardment of dust grain ice mantles.

Moreover, by simulating well-constrained experiments rather than interstellar environments we have been afforded a unique opportunity to compare the accuracy and physical realism of several approaches to modeling bulk chemistry over a variety of temperatures relevant to the ISM. As reported here, we have found that the standard approaches to bulk chemistry based on thermal diffusion or quantum tunneling through diffusion barriers less accurately reproduced the experimental data—particularly at low temperatures—than our model in which radicals and atomic oxygen were assumed to react quickly with neighboring species in the ice. This finding is in agreement with recent experiments by Ghesquière et al. (2018), who found no evidence for true bulk diffusion.

Regrettably, despite the large body of work in laboratory astrophysics on the irradiation of interstellar ice-analogs, it has not been possible until recently to incorporate many of the results of these experiments into astrochemical codes (Shingledecker et al. 2018). However, our work presented here shows that not only can such models simulate radiation-chemical reactions, they might even be fruitfully used as a replacement for the simpler kinetic models sometimes employed (e.g., Baragiola et al. 1999; Gomis et al. 2004b) in understanding and analyzing experimental data.

In summary, this study represents an attempt to shrink the existing gap between experiments and models, an increasingly urgent task in light of the upcoming launch of JWST. However, there is ample opportunity for even further refinements to our approach considering, for example, the effects of the implantation and subsequent reactions of the bombarding H⁺ ions, the effects of ice heating along the particle track, or the effects caused by the nuclear/elastic component of the stopping, which begins to dominate over the electronic/inelastic component considered here at lower particle energies (Spinks & Woods 1990). In addition to the synthesis of molecules, charged particle bombardment is known to drive the non-thermal desorption of even large molecules such as benzene (Thrower et al. 2011; Marchione et al. 2016). From experiments it is known that, particularly in water ice, excitons migrating to the surface represent one such mechanism that can stimulate this desorption. Given lingering questions about how molecules formed in dust-grain ice mantles are introduced into the surrounding gas in cold environments, future improvements to our approach in this area are warranted. Finally, experiments in which the abundances of multiple species are followed during irradiation would further advance our knowledge of radiation-chemical processes in ices and help to reduce uncertainties in future modeling research. In particular, Raman spectroscopic analysis—where even the behavior of IR-inactive species like O₂ could be monitored—represents a powerful yet perhaps underutilized technique that should be considered in future studies.

C.N.S. gratefully acknowledges the support of the Alexander von Humboldt Foundation. E.H. acknowledges the support of the National Science Foundation through grant AST-1514844. A.V. acknowledges the support of the Russian Science Foundation through grant 18-12-00351.

Software: MONACO (Vasyunin et al. 2017).

The radiolysis processes and radiochemical yields used here are listed in Table 3.

Table 3.  Radiolysis Reactions, Branching Fractions, and Yields (G-values) Used in This Work

Process	fbr	G-valuea	Type
H2O
${{\rm{H}}}_{2}{\rm{O}}\ \rightsquigarrow \ {\mathrm{OH}}^{* }\ +\ {{\rm{H}}}^{* }$	0.9	2.57 × 10−1	I
${{\rm{H}}}_{2}{\rm{O}}\ \rightsquigarrow \ {{\rm{O}}}^{* }\ +\ {{\rm{H}}}_{2}^{* }$	0.1	2.57 × 10−1	I
${{\rm{H}}}_{2}{\rm{O}}\ \rightsquigarrow \ \mathrm{OH}\ +\ {\rm{H}}$	0.9	1.58 × 10−1	II
${{\rm{H}}}_{2}{\rm{O}}\ \rightsquigarrow \ {\rm{O}}\ +\ {{\rm{H}}}_{2}$	0.1	1.58 × 10−1	II
${{\rm{H}}}_{2}{\rm{O}}\ \rightsquigarrow \ {{\rm{H}}}_{2}{{\rm{O}}}^{* }$	1.0	1.58 × 10−1	III
O2
${{\rm{O}}}_{2}\ \rightsquigarrow \ {{\rm{O}}}^{* }\ +\ {{\rm{O}}}^{* }$	1.0	6.81 × 100	I
${{\rm{O}}}_{2}\ \rightsquigarrow \ {\rm{O}}\ +\ {\rm{O}}$	1.0	2.91 × 100	II
${{\rm{O}}}_{2}\ \rightsquigarrow \ {{\rm{O}}}_{2}^{* }$	1.0	2.91 × 100	III
O3
${{\rm{O}}}_{3}\ \rightsquigarrow \ {{\rm{O}}}_{2}^{* }\ +\ {{\rm{O}}}^{* }$	1.0	5.57 × 101	I
${{\rm{O}}}_{3}\ \rightsquigarrow \ {{\rm{O}}}_{2}\ +\ {\rm{O}}$	1.0	1.80 × 101	II
${{\rm{O}}}_{3}\ \rightsquigarrow \ {{\rm{O}}}_{3}^{* }$	1.0	1.80 × 101	III
HO2
${\mathrm{HO}}_{2}\ \rightsquigarrow \ {{\rm{O}}}^{* }\ +\ {\mathrm{OH}}^{* }$	0.5	3.70 × 100	I
${\mathrm{HO}}_{2}\ \rightsquigarrow \ {{\rm{H}}}^{* }\ +\ {{\rm{O}}}_{2}^{* }$	0.5	3.70 × 100	I
${\mathrm{HO}}_{2}\ \rightsquigarrow \ {\rm{O}}\ +\ \mathrm{OH}$	0.5	3.71 × 100	II
${\mathrm{HO}}_{2}\ \rightsquigarrow \ {\rm{H}}\ +\ {{\rm{O}}}_{2}$	0.5	3.71 × 100	II
${\mathrm{HO}}_{2}\ \rightsquigarrow \ {\mathrm{HO}}_{2}^{* }$	1.0	3.71 × 100	III
H2O2
${{\rm{H}}}_{2}{{\rm{O}}}_{2}\ \rightsquigarrow \ {\mathrm{OH}}^{* }\ +\ {\mathrm{OH}}^{* }$	0.5	5.10 × 101	I
${{\rm{H}}}_{2}{{\rm{O}}}_{2}\ \rightsquigarrow \ {{\rm{H}}}^{* }\ +\ {\mathrm{HO}}_{2}^{* }$	0.5	5.10 × 101	I
${{\rm{H}}}_{2}{{\rm{O}}}_{2}\ \rightsquigarrow \ \mathrm{OH}\ +\ \mathrm{OH}$	0.5	4.13 × 101	II
${{\rm{H}}}_{2}{{\rm{O}}}_{2}\ \rightsquigarrow \ {\rm{H}}\ +\ {\mathrm{HO}}_{2}$	0.5	4.13 × 101	II
${{\rm{H}}}_{2}{{\rm{O}}}_{2}\ \rightsquigarrow \ {{\rm{H}}}_{2}{{\rm{O}}}_{2}^{* }$	1.0	4.13 × 101	III
O
${\rm{O}}\ \rightsquigarrow \ {{\rm{O}}}^{* }$	1.0	3.70 × 100	I
${\rm{O}}\ \rightsquigarrow \ {{\rm{O}}}^{* }$	1.0	1.93 × 100	III
H2
${{\rm{H}}}_{2}\ \rightsquigarrow \ {{\rm{H}}}^{* }\ +\ {{\rm{H}}}^{* }$	1.0	3.70 × 101	I
${{\rm{H}}}_{2}\ \rightsquigarrow \ {\rm{H}}\ +\ {\rm{H}}$	1.0	1.02 × 101	II
${{\rm{H}}}_{2}\ \rightsquigarrow \ {{\rm{H}}}_{2}^{* }$	1.0	1.02 × 101	III
OH
$\mathrm{OH}\ \rightsquigarrow \ {{\rm{O}}}^{* }\ +\ {{\rm{H}}}^{* }$	1.0	3.70 × 100	I
$\mathrm{OH}\ \rightsquigarrow \ {\rm{O}}\ +\ {\rm{H}}$	1.0	5.66 × 100	II
$\mathrm{OH}\ \rightsquigarrow \ {\mathrm{OH}}^{* }$	1.0	5.66 × 100	III

Notes. Following Shingledecker et al. (2018) we use Type I to indicate process (P2), Type II to indicate process (P3), and Type III to indicate process (P4).

^aGiven in molecules/100 eV.

Download table as:  ASCIITypeset image

The chemical network used in this work consists of the reactions listed in Table 4, as well as the neutral–neutral reactions in Table 7 of Shingledecker et al. (2017).

Table 4.  Solid-phase Reactions Used in Water Simulations

Reaction	EA (K)	fbr	Source
${\rm{H}}+{\rm{H}}\longrightarrow {{\rm{H}}}_{2}$	0	1.0	Wakelam et al. (2015)
${\rm{H}}+{{\rm{H}}}_{2}\longrightarrow {{\rm{H}}}_{2}+{\rm{H}}$	1900	1.0	Cohen et al. (1992)
${\rm{H}}+{\rm{O}}\longrightarrow \mathrm{OH}$	0	1.0	Atkinson et al. (2004)
${\rm{H}}+{{\rm{O}}}_{2}\longrightarrow {\mathrm{HO}}_{2}$	0	1.0	JL
${\rm{H}}+{{\rm{O}}}_{3}\longrightarrow {{\rm{O}}}_{2}+\mathrm{OH}$	450	1.0	Tielens & Hagen (1982)
${\rm{H}}+\mathrm{OH}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}$	0	1.0	Atkinson et al. (2004)
${\rm{H}}+{{\rm{H}}}_{2}{\rm{O}}\longrightarrow \mathrm{OH}+{{\rm{H}}}_{2}$	9700	1.0	Baulch et al. (1992)
${\rm{H}}+{\mathrm{HO}}_{2}\longrightarrow {\rm{O}}+{{\rm{H}}}_{2}{\rm{O}}$	0	0.0194	Atkinson et al. (2004)
${\rm{H}}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{O}}}_{2}+{{\rm{H}}}_{2}$	0	0.0857	Atkinson et al. (2004)
${\rm{H}}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}$	0	0.0894	Atkinson et al. (2004)
${\rm{H}}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+\mathrm{OH}$	1400	0.999	Baulch et al. (1992)
${\rm{H}}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {\mathrm{HO}}_{2}+{{\rm{H}}}_{2}$	1900	0.0006	Baulch et al. (1992)
${{\rm{H}}}_{2}+{\rm{O}}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}$	9700	1.0	Javoy et al. (2003)
${{\rm{H}}}_{2}+{{\rm{O}}}_{2}\longrightarrow {\mathrm{HO}}_{2}+{\rm{H}}$	28,000	1.0	Karkach & Osherov (1999)
${{\rm{H}}}_{2}+{{\rm{O}}}_{3}\longrightarrow \mathrm{OH}+{{\rm{O}}}_{2}$	10,000	1.0	See text
${{\rm{H}}}_{2}+\mathrm{OH}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+{\rm{H}}$	1800	1.0	Burkholder et al. (2014)
${{\rm{H}}}_{2}+{{\rm{H}}}_{2}{\rm{O}}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+{\rm{H}}$	10,000	1.0	See text
${{\rm{H}}}_{2}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}+{\rm{H}}$	13,000	1.0	Tsang & Hampson (1986)
${{\rm{H}}}_{2}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {{\rm{H}}}_{2}+\mathrm{OH}+\mathrm{OH}$	10,000	1.0	See text
${\rm{O}}+{{\rm{O}}}_{2}\longrightarrow {{\rm{O}}}_{3}$	180	1.0	Benderskii & Wight (1996)
${\rm{O}}+\mathrm{OH}\longrightarrow {\mathrm{HO}}_{2}$	0	1.0	Atkinson et al. (2004)
${\rm{O}}+{{\rm{H}}}_{2}{\rm{O}}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}$	8800	1.0	Karkach & Osherov (1999)
${\rm{O}}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{O}}}_{2}+\mathrm{OH}$	0	1.0	M84
${\rm{O}}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow \mathrm{OH}+{\mathrm{HO}}_{2}$	2000	1.0	Atkinson et al. (2004)
${{\rm{O}}}_{2}+\mathrm{OH}\longrightarrow {\mathrm{HO}}_{2}+{\rm{O}}$	25,000	1.0	Srinivasan et al. (2005)
${{\rm{O}}}_{2}+{{\rm{H}}}_{2}{\rm{O}}\longrightarrow {\mathrm{HO}}_{2}+\mathrm{OH}$	37,000	1.0	Mayer & Schieler (1968)
${{\rm{O}}}_{2}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{O}}}_{3}+\mathrm{OH}$	10,000	1.0	See text
${{\rm{O}}}_{2}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {\mathrm{HO}}_{2}+{\mathrm{HO}}_{2}$	18,000	1.0	Donaldson & Francisco (2003)
${{\rm{O}}}_{3}+\mathrm{OH}\longrightarrow {\mathrm{HO}}_{2}+{{\rm{O}}}_{2}$	940	1.0	Atkinson et al. (2004)
${{\rm{O}}}_{3}+{{\rm{H}}}_{2}{\rm{O}}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}+{{\rm{O}}}_{2}$	10,000	1.0	See text
${{\rm{O}}}_{3}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{O}}}_{2}+{{\rm{O}}}_{2}+\mathrm{OH}$	490	1.0	Burkholder et al. (2014)
${{\rm{O}}}_{3}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {{\rm{O}}}_{2}+{{\rm{O}}}_{2}+{{\rm{H}}}_{2}{\rm{O}}$	10,000	1.0	See text
$\mathrm{OH}+\mathrm{OH}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}$	0	1.0	Atkinson et al. (2004)
$\mathrm{OH}+{{\rm{H}}}_{2}{\rm{O}}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}+{\rm{H}}$	40,000	0.5	Lamberts & Kästner (2017)
$\mathrm{OH}+{{\rm{H}}}_{2}{\rm{O}}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+\mathrm{OH}$	2100	0.5	Uchimaru et al. (2003)
$\mathrm{OH}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+{{\rm{O}}}_{2}$	0	1.0	Schwab et al. (1989)
$\mathrm{OH}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+{\mathrm{HO}}_{2}$	760	1.0	Buszek et al. (2012)
${{\rm{H}}}_{2}{\rm{O}}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}+\mathrm{OH}$	17,000	1.0	Lloyd (1974)
${{\rm{H}}}_{2}{\rm{O}}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow \mathrm{OH}+\mathrm{OH}+{{\rm{H}}}_{2}{\rm{O}}$	10,000	1.0	See text
${\mathrm{HO}}_{2}+{\mathrm{HO}}_{2}\longrightarrow {{\rm{H}}}_{2}{{\rm{O}}}_{2}+{{\rm{O}}}_{2}$	0	1.0	See text
${{\rm{H}}}_{2}{{\rm{O}}}_{2}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}\longrightarrow {{\rm{H}}}_{2}{\rm{O}}+\mathrm{OH}+{\mathrm{HO}}_{2}$	10,000	1.0	See text

Note. In the absence of other experimental or theoretical data, we here assume E_A = 0 K for radical–radical reactions and 10⁴ K for all others.

Download table as:  ASCIITypeset image