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

TL;DR

This study enhances astrochemical models to simulate radiation-induced chemistry in interstellar ices, successfully reproducing experimental results and providing insights into bulk ice reactivity under ion bombardment.

Contribution

We developed and tested a new method within astrochemical models to accurately simulate radiation chemistry effects in bulk interstellar ices, aligning with laboratory experiments.

Findings

Successfully reproduced O$_3$ abundance in irradiated O$_2$ ice.

Accurately modeled H$_2$O$_2$ formation and temperature dependence.

Highlighted the importance of rapid, non-diffusive reactions of radiolysis products.

Abstract

Many astrochemical models today explicitly consider the species that comprise the bulk of interstellar dust grain ice-mantles separately from those in the top few monolayers. Bombardment of these ices by ionizing radiation - whether in the form of cosmic rays, stellar winds, or radionuclide emission - represents an astrochemically viable means of driving a rich chemistry even in the bulk of the ice-mantle, now supported by a large body of work in laboratory astrophysics. In this study, using an existing rate equation-based astrochemical code modified to include a method of considering radiation chemistry recently developed by us, we attempted to simulate two such studies in which (a) pure O$_2$ ice at 5 K and, (b) pure H$_2$O ice at 16 K and 77 K, were bombarded by keV H$^+$ ions. Our aims are twofold: (1) to test the capability of our newly developed method to replicate the results of ice-irradiation experiments, and (2) to determine in such a well-constrained system how bulk chemistry is best handled using the same gas-grain codes that are used to model the interstellar medium (ISM). We find that our modified astrochemical model is able to reproduce both the abundance of O$_3$ in the 5 K pure O$_2$ ice, as well as both the abundance of H$_2$O$_2$ in the 16 K water ice and the previously noted decrease of hydrogen peroxide at higher temperatures. However, these results require the assumption that radicals and other reactive species produced via radiolysis react quickly and non-diffusively with neighbors in the ice.

Figures

Tables

Table 1. Table 1: Physical parameters used in simulations of experiments.

Parameter	Value	Source
\ceO2
Temperature	5 K	Baragiola et al. (1999)
Density	$5.7\times {10}^{22}$ cm-3	Horl & Kohlbeck (1982)
\ceH+ Energy	100 keV	Baragiola et al. (1999)
$S_e$	$4.6\times {10}^{-14}$ eV cm2	Baragiola et al. (1999)
\ceH2O
Temperature	16 K, 77 K	Gomis et al. (2004b)
Density	$3.1\times {10}^{22}$ cm-3	Martoňák et al. (2005)
\ceH+ Energy	200 keV	Gomis et al. (2004b)
$S_e$	$3.3\times {10}^{-14}$ eV cm2	Uehara et al. (2000)

Table 2. Table 2: Desorption energies for relevant species.

Species	$E_D$ (K)	Source
\ceO	1 660aaUsed in models (a) and (b)	He et al. (2015)
	1 040bbUsed in model (c)	Minissale et al. (2013)
\ceO2	930	Jing et al. (2012)
\ceO3	1 833	Jing et al. (2012)
\ceH	450	Garrod & Herbst (2006)
(Thrower et al., 2011; Marchione et al., 2016) \ceH2	430	Garrod & Herbst (2006)
\ceOH	2 850	Garrod & Herbst (2006)
\ceH2O	5 700	Garrod & Herbst (2006)
\ceHO2	4 510	$E_D^O+E_D^{\mathrm{OH}}$
\ceH2O2	5 700	$E_D^{\mathrm{OH}}+E_D^{\mathrm{OH}}$

Table 3. Table A1: Radiolysis reactions, branching fractions, and yields ( G-values) used in this work

Process	$f_{\mathrm{br}}$	$G$-valueaaGiven in molecules/100 eV	Type
\ceH2O
$H_{2}O\rightsquigarrow {\mathrm{OH}}^{\ast }+H^{\ast }$	0.9	$2.57\times {10}^{-1}$	I
$H_{2}O\rightsquigarrow O^{\ast }+H_2^{\ast }$	0.1	$2.57\times {10}^{-1}$	I
$H_{2}O\rightsquigarrow \mathrm{OH}+H$	0.9	$1.58\times {10}^{-1}$	II
$H_{2}O\rightsquigarrow O+H_2$	0.1	$1.58\times {10}^{-1}$	II
$H_{2}O\rightsquigarrow H_2{O}^{\ast }$	1.0	$1.58\times {10}^{-1}$	III
\ceO2
$O_2\rightsquigarrow O^{\ast }+O^{\ast }$	1.0	$6.81\times {10}^0$	I
$O_2\rightsquigarrow O+O$	1.0	$2.91\times {10}^0$	II
$O_2\rightsquigarrow O_2^{\ast }$	1.0	$2.91\times {10}^0$	III
\ceO3
$O_3\rightsquigarrow O_2^{\ast }+O^{\ast }$	1.0	$5.57\times {10}^1$	I
$O_3\rightsquigarrow O_2+O$	1.0	$1.80\times {10}^1$	II
$O_3\rightsquigarrow O_3^{\ast }$	1.0	$1.80\times {10}^1$	III
\ceHO2
${\mathrm{HO}}_2\rightsquigarrow O^{\ast }+{\mathrm{OH}}^{\ast }$	0.5	$3.70\times {10}^0$	I
${\mathrm{HO}}_2\rightsquigarrow H^{\ast }+O_2^{\ast }$	0.5	$3.70\times {10}^0$	I
${\mathrm{HO}}_2\rightsquigarrow O+\mathrm{OH}$	0.5	$3.71\times {10}^0$	II
${\mathrm{HO}}_2\rightsquigarrow H+O_2$	0.5	$3.71\times {10}^0$	II
${\mathrm{HO}}_2\rightsquigarrow {\mathrm{HO}}_2^{\ast }$	1.0	$3.71\times {10}^0$	III
\ceH2O2
$H_2{O}_2\rightsquigarrow {\mathrm{OH}}^{\ast }+{\mathrm{OH}}^{\ast }$	0.5	$5.10\times {10}^1$	I
$H_2{O}_2\rightsquigarrow H^{\ast }+{\mathrm{HO}}_2^{\ast }$	0.5	$5.10\times {10}^1$	I
$H_2{O}_2\rightsquigarrow \mathrm{OH}+\mathrm{OH}$	0.5	$4.13\times {10}^1$	II
$H_2{O}_2\rightsquigarrow H+{\mathrm{HO}}_2$	0.5	$4.13\times {10}^1$	II
$H_2{O}_2\rightsquigarrow H_2{O}_2^{\ast }$	1.0	$4.13\times {10}^1$	III
\ceO
$O\rightsquigarrow O^{\ast }$	1.0	$3.70\times {10}^0$	I
$O\rightsquigarrow O^{\ast }$	1.0	$1.93\times {10}^0$	III
\ceH2
$H_2\rightsquigarrow H^{\ast }+H^{\ast }$	1.0	$3.70\times {10}^1$	I
$H_2\rightsquigarrow H+H$	1.0	$1.02\times {10}^1$	II
$H_2\rightsquigarrow H_2^{\ast }$	1.0	$1.02\times {10}^1$	III
\ceOH
$\mathrm{OH}\rightsquigarrow O^{\ast }+H^{\ast }$	1.0	$3.70\times {10}^0$	I
$\mathrm{OH}\rightsquigarrow O+H$	1.0	$5.66\times {10}^0$	II
$\mathrm{OH}\rightsquigarrow {\mathrm{OH}}^{\ast }$	1.0	$5.66\times {10}^0$	III

Table 4. Table B1: Solid-phase reactions used in water simulations

Reaction	$E_A$ (K)	$f_{\mathrm{br}}$	Source
\ceH + H -¿ H2	0	1.0	Wakelam et al. (2015)
\ceH + H2 -¿ H2 + H	1 900	1.0	Cohen et al. (1992)
\ceH + O -¿ OH	0	1.0	Atkinson et al. (2004)
\ceH + O2 -¿ HO2	0	1.0	JL
\ceH + O3 -¿ O2 + OH	450	1.0	Tielens & Hagen (1982)
\ceH + OH -¿ H2O	0	1.0	Atkinson et al. (2004)
\ceH + H2O -¿ OH + H2	9 700	1.0	Baulch et al. (1992)
\ceH + HO2 -¿ O + H2O	0	0.0194	Atkinson et al. (2004)
\ceH + HO2 -¿ O2 + H2	0	0.0857	Atkinson et al. (2004)
\ceH + HO2 -¿ H2O2	0	0.0894	Atkinson et al. (2004)
\ceH + H2O2 -¿ H2O + OH	1 400	0.999	Baulch et al. (1992)
\ceH + H2O2 -¿ HO2 + H2	1 900	0.0006	Baulch et al. (1992)
\ceH2 + O -¿ H2O	9 700	1.0	Javoy et al. (2003)
\ceH2 + O2 -¿ HO2 + H	28 000	1.0	Karkach & Osherov (1999)
\ceH2 + O3 -¿ OH + O2	10 000	1.0	See text
\ceH2 + OH -¿ H2O + H	1 800	1.0	Burkholder et al. (2014)
\ceH2 + H2O -¿ H2O + H	10 000	1.0	See text
\ceH2 + HO2 -¿ H2O2 + H	13 000	1.0	Tsang & Hampson (1986)
\ceH2 + H2O2 -¿ H2 + OH + OH	10 000	1.0	See text
\ceO + O2 -¿ O3	180	1.0	Benderskii & Wight (1996)
\ceO + OH -¿ HO2	0	1.0	Atkinson et al. (2004)
\ceO + H2O -¿ H2O2	8 800	1.0	Karkach & Osherov (1999)
\ceO + HO2 -¿ O2 + OH	0	1.0	M84
\ceO + H2O2 -¿ OH + HO2	2 000	1.0	Atkinson et al. (2004)
\ceO2 + OH -¿ HO2 + O	25 000	1.0	Srinivasan et al. (2005)
\ceO2 + H2O -¿ HO2 + OH	37 000	1.0	Mayer & Schieler (1968)
\ceO2 + HO2 -¿ O3 + OH	10 000	1.0	See text
\ceO2 + H2O2 -¿ HO2 + HO2	18 000	1.0	Donaldson & Francisco (2003)
\ceO3 + OH -¿ HO2 + O2	940	1.0	Atkinson et al. (2004)
\ceO3 + H2O -¿ H2O2 + O2	10 000	1.0	See text
\ceO3 + HO2 -¿ O2 + O2 + OH	490	1.0	Burkholder et al. (2014)
\ceO3 + H2O2 -¿ O2 + O2 + H2O	10 000	1.0	See text
\ceOH + OH -¿ H2O2	0	1.0	Atkinson et al. (2004)
\ceOH + H2O -¿ H2O2 + H	40 000	0.5	Lamberts & Kästner (2017)
\ceOH + H2O -¿ H2O + OH	2 100	0.5	Uchimaru et al. (2003)
\ceOH + HO2 -¿ H2O + O2	0	1.0	Schwab et al. (1989)
\ceOH + H2O2 -¿ H2O + HO2	760	1.0	Buszek et al. (2012)
\ceH2O + HO2 -¿ H2O2 + OH	17 000	1.0	Lloyd (1974)
\ceH2O + H2O2 -¿ OH + OH + H2O	10 000	1.0	See text
\ceHO2 + HO2 -¿ H2O2 + O2	0	1.0	See text
\ceH2O2 + H2O2 -¿ H2O + OH + HO2	10 000	1.0	See text

Equations

$$\frac{\mathrm{d}E}{\mathrm{d}t}\approx n(S_{\mathrm{n}}+S_{\mathrm{e}})$$

$$A⇝A^{+}+e^{-}$$

$$A⇝A^{+}+e^{-}\to A^{\ast }\to B^{\ast }+C^{\ast }$$

$$A⇝A^{\ast }\to B+C$$

$$A⇝A^{\ast }$$

$$k_{\mathrm{rad}}=\left(\frac{G}{100\mathrm{eV}}\right)S_{\mathrm{e}}\varphi ,$$

$$k_{\mathrm{fast}}=f_{\mathrm{br}}\left[\frac{{\nu }_0^A+{\nu }_0^B}{N_{\mathrm{bulk}}}\right]\exp \left(-\frac{E_{\mathrm{act}}^{AB}}{T_{\mathrm{ice}}}\right),$$

$${\nu }_0^A=\sqrt{\frac{2n_{\mathrm{s}}{E}_{\mathrm{D}}^A}{{\pi }^2{m}_A}},$$

$$\frac{d{N}_{B^{\ast }}}{dt}=k_{\mathrm{rad}}^{\text{ref}p2}{N}_A-{\nu }_0^B{N}_{B^{\ast }}-k_{\mathrm{fast}}{N}_X{N}_{B^{\ast }},$$

$${\mathrm{A}}^{\ast }+\mathrm{B}\to products$$

$$\mathrm{A}+{\mathrm{B}}^{\ast }\to products.$$

$$\text{ce}OH+H2O2->H2O+HO2$$

$$\text{ce}OH+OH->H2O2$$

$$\text{ce}H+HO2->H2O2,$$

$$\text{ce}H+O2->HO2,$$

$$\text{ce}HO2+O3->O2+O2+OH.$$

$$\text{ce}OH+HO2->H2O+O2$$

Full text

Our aims are twofold: (1) to test the capability of our newly developed method to replicate the results of ice-irradiation experiments, and (2) to determine in such a well-constrained system how bulk chemistry is best handled using the same gas-grain codes that are used to model the interstellar medium (ISM). We find that our modified astrochemical model is able to reproduce both the abundance of \ceO3 in the 5 K pure \ceO2 ice, as well as both the abundance of \ceH2O2 in the 16 K water ice and the previously noted decrease of hydrogen peroxide at higher temperatures. However, these results require the assumption that radicals and other reactive species produced via radiolysis react quickly and non-diffusively with neighbors in the ice.

astrochemistry - molecular processes - cosmic rays

††journal: ApJ††software: MONACO (Vasyunin et al., 2017)

1 Introduction

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, though, calls for models that are equal to the task and, unfortunately, there remain significant uncertainties 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 in the selvedge rests on the following three main lines of reasoning:

1. (a)

species in the bulk are more tightly bound, 2. (b)

a richer chemistry is possible in the selvedge, given the constant exchange of species with the surrounding gas, and 3. (c)

photo-processes, 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 that 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 quickly occur among neighboring species, without the need of 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 photo-processes 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 (Gerakines et al., 2001; Spinks & Woods, 1990), 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 (Madey et al., 2002; Hudson & Moore, 2001), 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 (Hudson et al., 2008; Holtom et al., 2005; Kaiser & Roessler, 1997; 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:

[TABLE]

where $n$ is the density of the target (cm*-3*) and $S_{\mathrm{n}}$ and $S_{\mathrm{e}}$ are, respectively, the so-called nuclear/elastic and electronic/inelastic stopping cross sections (eV-cm2). 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 (Spinks & Woods, 1990; Johnson, 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, if any, effects 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 \ceO2 ice by 100 keV protons as described in Baragiola et al. (1999) as well as that of a pure water ice at both 16 K and 77 K by 200 keV protons, as reported in Gomis et al. (2004a). Doing so has not only afforded us the opportunity 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 the 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 rest of the paper is as follows: in §2 we discuss details of both the model and chemical network used in this work. In §3 our major results are described and their significance is discussed. Finally, our conclusions and suggestions for future work are presented in §4.

2 Model and Network

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$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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 to the system from bombarding particles (Dewhurst et al., 1952; Spinks & Woods, 1990). A full list of the radiolysis processes used here is given in Table A1 of Appendix A. The yields given there were arrived at through a two-step process where first, initial values were estimated usnig 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

[TABLE]

with $\phi$ being the radiation flux and $S_{\mathrm{e}}$ the previously described inelastic/electronic stopping cross section. In this work, we assume the flux to be $1.0\times 10^{11}$ \ceH+ $\mathrm{cm}^{-2}$ $\mathrm{s^{-1}}$, 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

[TABLE]

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

[TABLE]

with $E^{A}_{\mathrm{D}}$ and $m_{A}$ being the desorption energy and mass of $A$, respectively, and $n_{\mathrm{s}}$ the phsysisorption site density - here $6.55\times 10^{14}$ cm*-2* (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$, the rate of change in the number of $B^{*}$ as a function of time can be described by

[TABLE]

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:

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 Eq. (3), 2. (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, and 3. (c)

a set of models similar to (b) but with tunneling through diffusion barriers allowed for H, \ceH2, and O.

For the third set of simulations, we treat tunneling under diffusion barriers by H and \ceH2 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 \ceH2 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 B1 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^{4}$ K, except for radical-radical reactions, which were assumed to be barrierless. For every thermal reaction of the type given in Table B1, we included variants involving one suprathermal reactant. To illustrate this point, for every reaction of the form $\mathrm{A+B\rightarrow}\;products$, we include the following suprathermal variants

[TABLE]

[TABLE]

A number of reactions listed in Table B1 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 $\sim 3.4$ eV, a realistic energy for an electronically excited suprathermal species.

3 Results and Discussion

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 ice-irradiation systems described in two previous radiation-chemical studies discussed above - the results of which are described below.

3.1 H+ Bombardment of Pure O2 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 \ceO2 ice at 5 K by 100 keV \ceH+, following the study by Baragiola et al. (1999), which has previously been successfully modeled by us using the much more detailed Monte Carlo code, CIRIS (Shingledecker et al., 2017). An irradiated \ceO2 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, \ceO2, and \ceO3. However, in the Monte Carlo code even this seemingly simple system required a network of $\sim 50$ reactions or processes and involved a level of detail beyond the practical capabilities of rate equations-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 Fig. 1a, even our simplified approach is capable of reasonably reproducing the concentration of \ceO3 as a function of fluence. Here, fluence - the product of the \ceH+ flux and irradiation time - represents the number of ions that have bombarded the ice per cm2. As can be seen in Fig. 1, the best results were obtained in model (a), where 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 under-prediction of \ceO3 at fluences under $\sim 10^{13}$ \ceH+ cm*-1* 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 Fig. 1a is the fact that atomic oxygen abundances remained low throughout the simulation, despite its continuous production from the dissociation of \ceO2 and, to a lesser degree, \ceO3.

Comparison with Fig. 1b reveals that model (b), using standard Langmuir-Hinshelwood rate coefficients, is the least accurate. There, though the modeled \ceO3 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_{\mathrm{b}}=0.5E_{\mathrm{D}}$). Tunneling for H and \ceH2 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 Fig. 1c, with 1b 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 best agreement with experimental data - both the measured \ceO3 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 [\ceO3…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 \ceO3 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, $\nu_{0}$, thereby increasing the abundance of O. Here, we use the characteristic vibrational frequency given in Eq. (4) (which typically has a value on the order of $10^{12}$ s*-1*) as a rough approximation, though in fact, the value of this parameter can be as low as $10^{-3}$ s*-1* (Theulé et al., 2013), and rate coefficients obtained using (4) should probably be interpreted as upper limits.

3.2 H+ Bombardment of Pure H2O Ice

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

Shown in Fig. 2 are both our calculated \ceH2O2 abundances vs proton fluence as well as the approximate steady-state bf experimental values from Gomis et al. (2004a) relative to \ceH2O of $\sim 1.0$% at 16 K and $0.25$% at 77 K. As with the pure \ceO2 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

[TABLE]

We find that, in agreement with Moore & Hudson (2000), reaction (R3) is indeed behind the drop in \ceH2O2 abundance at $\sim 80$ K, since, 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 \ceH2O2. As can be seen in Fig. 3, without the assumption that radicals react quickly in the bulk, OH, in particular, becomes quite abundant in the ice, even at 77 K.

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

[TABLE]

[TABLE]

where here, the higher rate of \ceH2O2 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, \ceHO2, via

[TABLE]

thereby contributing to the even higher \ceH2O2 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$(\ceH2O2), the radiolytic hydrogen peroxide formation yield. Unfortunately, we cannot simply compare the $G$-values given in Table A1 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 Fig. 5 over the pre-steady-state regime - corresponding to doses of ca. 0-10 eV, where dose is the product of the fluence and $S_{\mathrm{e}}$. From this, we calculate the yield of H2O2 at 16 K to be $0.1$ molecules/100 eV - exactly what was mentioned by Moore & Hudson (2000) as the yield in pure \ceH2O.

Shown in Fig. 6 are the abundances of \ceH2, \ceO2, and \ceO3 versus fluence. In model (a), the abundance of \ceO2 at 16 K is kept low because of destruction via (R6) to form \ceHO2 but increases at 77 K due, in part, to more efficient formation via

[TABLE]

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

[TABLE]

where the abundances of OH and \ceHO2 are enhanced relative to model (b) because of the effects of* quantum tunneling through diffusion barriers by H, \ceH2, and O - as shown in Fig. 3. The decreased abundance of these radicals at 77 K in model (b), combined with destruction with atomic oxygen, leads to the drop in [\ceO2] in Fig. 6b.

Unfortunately, comparison of our \ceO2 results with experimental data, as with \ceH2O2, is complicated by the fact that homonuclear diatomic molecules, lacking permanent dipoles, are IR inactive. Thus, their abundances cannot be measured using standard techniques, such as Fourier Transform Infrared Spectroscoppy (FTIR). Nevertheless, it is well known that \ceH2 and \ceO2 form during water ice radiolysis based on analysis of both sputtering products as well as post-irradiation temperature-programmed desorption (TPD) 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 \ceO2 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 Fig. 6, the maximum abundance of \ceO2 with respect to water achieved here is $\sim 0.1$% in model (a) at 77 K, a value which increased only negligibly at still higher temperatures. Thus, our models predict that the radiolysis of pure \ceH2O ice is not the dominant mechanism behind the $\sim 3.8$ % \ceO2 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\times 10^{-4}$ % relative to water was established. Interestingly, only in model (a) are the ozone abundances predicted to remain below this limit, even at 77 K.

4 Conclusions & Outlook

In this work, we have simulated the bombardment of pure \ceO2 and \ceH2O 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 A1 of Appendix A and the reactions noted in Appendix B. As illustrated in Figs. 1 and 2, we were able to qualitatively reproduce both the abundance of \ceO3 in pure \ceO2 (Baragiola et al., 1999) and \ceH2O2 in pure \ceH2O (Gomis et al., 2004a) utilizing the radiochemical processes given in Table A1. Thus validated, these processes, along with the reactions given in Table B1, 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 did more poorly at reproducing 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-analogues, 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. Gomis et al. (2004b); Baragiola et al. (1999)) 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 by, for example, considering the effects of the implantation and subsequent reactions of the bombarding \ceH+ ions, of ice heating along the particle track, or of 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 well 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 \ceO2 could be monitored - represents a powerful, yet perhaps underutilised, technique that should be considered in future studies.

CNS gratefully acknowledges the support of the Alexander von Humboldt Foundation. EH acknowledges the support of the National Science Foundation through grant AST-1514844. AV acknowledges the support of the Russian Science Foundation through grant 18-12-00351.

Appendix A Radiolysis Processes

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

\startlongtable

Appendix B Chemical Network

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

\startlongtable