Simulations of Ice Chemistry in Cometary Nuclei

As mentioned above, molecules are allowed to diffuse thermally between the different layers, through a process of bulk diffusion. Garrod (2013) formulated rates for this process in the case of a single bulk-ice phase topped by one monolayer of surface ice. Diffusion from the bulk ice into the surface layer for species i was dependent on a swapping rate coefficient, k_swap(i), for bulk-ice diffusion of a single molecule from one site to another adjacent site, determined by a diffusion barrier E_swap(i), such that

$$\begin{eqnarray}&&{k}_{\mathrm{swap}}(i)={\nu }_{0}(i)\exp [-{E}_{\mathrm{swap}}(i)/T],\end{eqnarray} \tag{ 1 }$$

where ν₀(i) is the characteristic vibrational frequency of a harmonic oscillator and T is the ice temperature. Surface and bulk-ice diffusion barriers, as well as desorption energies, are given by Garrod (2013) for various species used in this model. Note that this treatment assumes that the ice is both amorphous and chemically diverse, rather than having a pure water and/or crystalline structure, which could inhibit bulk diffusion.

In the present model, the rate of diffusion of species i into the surface from the layer beneath (of whatever thickness) is given by

$$\begin{eqnarray}&&{R}_{\mathrm{swap},{\rm{m}}}(i)={N}_{m}(i)\ \displaystyle \frac{{N}_{S}}{{N}_{M}}\ {k}_{\mathrm{swap}}(i)/6,\end{eqnarray} \tag{ 2 }$$

where N_m(i) is the mantle population of species i, and N_M and N_S are the total mantle and surface populations over all species (thus, in the comet model, N_S/N_M = 1/3). The factor of 6 corresponds to the six competing directions available for diffusion assuming a simple cubic structure.

The (swapping) diffusion of a mantle molecule into the surface is in each case balanced by the diffusion of a surface molecule into the layer beneath; the choice of molecule to undergo this exchange is determined proportionately to the fractional abundances of each species in the surface layer.

Thus, the overall net rate of diffusion of molecules of all kinds into and out of each layer is zero, while the net rate for individual chemical species may not be. Following Garrod (2013), the complementary rates for the above-mentioned exchange of (in this case) material from the surface layer into the layer beneath are given by

$$\begin{eqnarray}&&{R}_{{\rm{swap}},{\rm{s}}}(i)=\displaystyle \frac{{N}_{s}(i)}{{N}_{S}}\cdot \displaystyle \sum _{{all}\,j}{R}_{{\rm{swap}},{\rm{m}}}(j).\end{eqnarray} \tag{ 3 }$$

A similar treatment is extended to all layers in the comet model; however, diffusion initiated in each subsurface layer can lead to exchange with either the layer above or the layer below. Diffusion in each direction has its own version of Equation (2) for each species in each layer, as well as complementary rates for the exchange partners. In general, the rates of interlayer diffusion for all layers and chemical species will be mostly negligible at very low temperatures (<20 K). No diffusion into or out of the reservoir (via the lowest ice layer) is allowed.

Chemical reactions are allowed to occur within each ice layer, with the reactants meeting nominally through the same site-to-site diffusion process as outlined above and described by Equation (1). Rates of reaction between arbitrary species A and B are therefore determined by the expression

$$\begin{eqnarray}&&{R}_{\mathrm{AB}}={\kappa }_{\mathrm{AB}}{N}_{m}(A)\ {N}_{m}(B)\ [{k}_{\mathrm{swap}}(A)+{k}_{\mathrm{swap}}(B)]/{N}_{M \mbox{-} D},\end{eqnarray} \tag{ 4 }$$

where κ_AB represents the reaction efficiency in cases where an activation energy barrier exists, corresponding to the competition between reaction of the two species and their diffusion away from each other. Note that the denominator of Equation (4) contains not N_M, but ${N}_{M-D}={N}_{M}(1-{X}_{\mathrm{dust}})$, the total ice abundance in the mantle layer in question with dust excluded. This substitution is made to account for the fact that diffusive species in the ice are unable to penetrate the space occupied by the dust.

A similar expression to Equation (4) pertains to reactions in the surface layer, although here a lower diffusion barrier is used, and all surface diffusion sites are considered available; thus, the denominator contains simply N_S. Following Garrod (2013), the surface and bulk diffusion barriers are set to fractions 0.35 and 0.7, respectively, of the desorption energy, E_des(i). Many of the radical–radical and neutral H-abstraction reactions included in the chemical network are listed by Garrod et al. (2008) and Garrod (2013).

The above treatment for bulk-ice diffusion and the related treatment for diffusion-mediated reactions do not explicitly consider the presence of porosity in the ice (nor in the grains suspended within them). No infrared signatures of the OH dangling bond feature have so far been detected in the interstellar medium (Keane et al. 2001), suggesting that interstellar ices are not very porous; however, cometary ices are thought to be so, based on density determinations (e.g., Jorda et al. 2016; see also Sections 4, 5.2, and 5.3). This porosity may result from irradiation and the subsequent loss of volatiles (e.g., Johnson 1991). But it is unclear whether comets under cold-storage conditions should initially have significant ice porosity—nor, crucially, whether such pores would show significant connectivity throughout the overall ice structure on scales greater than, say, the characteristic size of a dust grain particle. The presence of a well-connected network of pores or larger cavities spanning different depth scales into the ice could allow a more rapid transport of material, via thermal diffusive hopping of individual molecules on those surfaces, which would be significantly faster than bulk processes. Also, even the presence of enclosed, unconnected pores would provide a surface upon which chemical reactions could occur more rapidly than within a bulk ice. In this initial cometary chemistry model, the mechanisms used in the interstellar treatments are retained. An explicit treatment of porosity and diffusion will be a focus for future work.

In addition to thermal desorption from the surface, which is governed by an expression similar to Equation (1), nonthermal desorption processes are also included for the surface layer (no desorption is allowed from deeper layers). Reactive desorption is included in the model according to the treatment of Garrod et al. (2007). Single-product surface reactions are thus assumed to result in desorption of the product in a fraction of cases no greater than 1% (i.e., a_RRK = 0.01).

UV photodesorption is included for each surface species, i, with rate coefficients set by the expression

$$\begin{eqnarray}&&{k}_{\mathrm{PD}}(i)=Y(i)\,\sigma \,{F}_{\mathrm{UV}}/2,\end{eqnarray} \tag{ 5 }$$

where σ is a generic geometric cross section equal to 10⁻¹⁵ cm²; F_UV is the flux of UV photons, set to 10⁸ cm⁻² s⁻¹ to represent the standard interstellar radiation field (with the factor of 2 corresponding to exposure to photons from one side only); and Y(i) is the photodesorption yield per photon, as determined by Öberg et al. (2009b). Where measured yields are unavailable (as in most cases), a value of 10⁻³ is assumed, in line with the order of magnitude of the measured values. In models where the UV field strength is altered (for the purposes of photodissociation), the value of F_UV above is also adjusted accordingly.

To approximate the composition of cometary ices at the beginning of the life of a comet, representative abundances are assumed based on the dominant, simple components of interstellar ices. Water is assumed to be the most abundant ice species, followed by CO and CO₂ (20% with respect to water), methanol (5%), and methane, formaldehyde, and ammonia (1% each). No other chemical species are initially assumed, although an inert dust component of unspecified composition is included in the model. Each of these abundances is assumed to be uniform throughout all layers of the comet at the beginning of the model. Since all chemical species in the initial ices are highly stable, the gradual evolution of chemistry within the ices in the model is driven by the dissociation of stable molecules into radicals or atoms that may react to form other species.

Sulfur-bearing species are not included in the present models, as no sulfur-bearing precursor is included in the initial abundances. Although H₂S has been detected in various comets with abundances in the range of 0.1%–1% H₂O (see Mumma & Charnley 2011), sulfur chemistry in dense and star-forming regions is still a matter of debate, due to the lack of a clear carrier for the majority of the sulfur budget. For reasons of simplicity, sulfur is therefore omitted from the model, which is nevertheless in keeping with the lack of detection of interstellar solid-phase H₂S. Future models will consider the inclusion of sulfur chemistry, by the adoption of more varied initial ice abundances including a range of sulfur-bearing species.

Unlike interstellar ices, which are simply formed on the surfaces of dust grains, the cometary ices are expected to contain dust grain particles themselves, interspersed with the ice. The cometary ice model must therefore include a component corresponding to dust in each layer. As well as taking up physical space, the dust will also act to extinguish the external UV field that could otherwise induce photodissociation of molecules in the ice.

In the model presented here, dust is treated as a distinct, inert chemical species that is present in every ice layer. It may not diffuse within the ice, and it is not allowed to undergo thermal desorption from the ice surface, but it may nevertheless be transferred between layers in the same way as other species, to compensate for gains and losses primarily caused by desorption of volatiles from the surface. In this way, it is possible for the dust to become more concentrated in the surface layers, as volatiles desorb while the dust is left behind.

The dust component is assigned a uniform initial abundance throughout the comet, which can be expressed as a fraction of each layer of the ice/dust mixture (i.e., a volume fraction). Because the majority of the cometary ice is composed of water and thus depends on the availability of interstellar oxygen, the determination of the dust fraction in this model is based on an estimate of the mass ratio of dust grain material to interstellar oxygen content. Assuming a canonical Milky Way gas-to-dust ratio (for atomic hydrogen) of 100:1, the volume ratio of dust to water ice is given by

$$\begin{eqnarray}&&\displaystyle \frac{{V}_{d}}{{V}_{w}}=\displaystyle \frac{0.01}{X({\rm{O}})\,{m}_{w}\,{F}_{w}}\displaystyle \frac{{\rho }_{w}}{{\rho }_{d}},\end{eqnarray} \tag{ 6 }$$

where X(O) is the elemental abundance of oxygen with respect to hydrogen nuclei, m_w is the mass ratio of water to atomic hydrogen (=18), F_w is the fraction of oxygen atoms locked up in water ice, and ρ_d and ρ_w are the mass densities of dust grain material and of water ice, respectively, the latter of which was assigned in Section 2 a value of 0.9 g cm⁻³. Based on the initial ice abundances provided in Section 2.2, if the total oxygen budget is present in the form of H₂O, CO, CO₂, H₂CO, or CH₃OH ice, then the amount locked up specifically in water is F_w = 100/166. Following past interstellar chemistry models (e.g., Hasegawa et al. 1992), a dust density of ρ_d = 3 g cm⁻³ is assumed. The remaining quantity required to evaluate Equation (6) is the overall interstellar oxygen abundance. Two values have commonly been assumed in past models: X(O) = 3.2 × 10⁻⁴ (Garrod et al. 2008), and the so-called low metal abundance value of X(O) = 1.76 × 10⁻⁴ used by Graedel et al. (1982). These result in respective volume ratios of V_d/V_w = 0.865 and 1.57, which are equivalent to ratios of dust to total ice content of 0.584 and 1.06 (based on the ice composition in Section 2.2), or ratios of dust to total-ice oxygen of 0.521 and 0.947. For the models presented here, an intermediate value of V_d/V_w = 1.11 is assumed, which corresponds to a ratio of dust to total ice of 111/148 = 0.75. Based on the assumed initial abundances of all ice species, dust therefore comprises ∼42.9% by volume of each layer of cometary material (dust+ice) at the beginning of the model. As with other species, the dust fraction in the reservoir beneath the lowest ice layer is held constant over time.

The precise initial value used here is not definitive, as it depends on the several assumptions above; the dust content could plausibly diverge by several tens of percent from the value assumed here. However, the absence or otherwise of porosity in the ices will not affect the above calculations.

Adopting a canonical interstellar dust grain radius of 0.1 μm, the grains are around 300 times larger than the individual molecules that make up the remainder of each layer. In the model presented here, the transfer of material from layer to layer, as required to counteract any net losses or gains in each layer, is treated in the same way for grain material as for molecules. However, this could in practice lead to a layer becoming entirely composed of grain material, filling even every molecule-sized space. The rigid structure of the grains evidently precludes this, so a maximum fractional population for grains is imposed on each layer. This maximum is assumed to be equal to the well-known maximum space-filling factor for identically sized spheres, i.e., $\pi /3\sqrt{2}\simeq 74 \%$ (e.g., Hales 2005). The dust component of each layer, which begins at a volume fractional abundance of ∼42.9%, may therefore concentrate by a factor of up to ∼1.73 in this model before reaching the geometrical maximum. While the ∼74% volume limit is invariant with the choice of grain size, the inclusion of a range of dust grain sizes could allow a higher threshold (see also Section 2.4 regarding the size distribution). The solution of the close-packing problem for multiple sphere sizes is not trivial; however, computational determinations of upper bounds for binary sphere packing indicate maxima around 86% (De Laat et al. 2014).

The model assumes that the simulated comet is exposed to a UV radiation field capable of inducing photodissociation within the ice, and which is similar in strength to the standard interstellar radiation field. Interstellar gas-phase photodissociation rates are used as a basis for the photodissociation of molecules in the ice layers (similarly to Garrod 2013, and earlier models).

In the interstellar model, photodissociation rates are defined by two parameters, a and c; these correspond, respectively, to the photodissociation rate coefficient (s⁻¹) produced by the unattenuated interstellar radiation field (of strength 1 Draine unit) and to a dimensionless modifier to the visual extinction (A_V) that accounts for the wavelength dependence of absorption by dust (with typical values in the range of 1–3), giving the expression

$$\begin{eqnarray}&&{k}_{\mathrm{PD}}(j)=a(j)\,\exp [-c(j)\,{A}_{{\rm{V}}}],\end{eqnarray} \tag{ 7 }$$

where k_PD(j) is the rate coefficient of photodissociation process j for some specific molecule.

In the present models, the surface layer of comet ice is exposed to the maximum field, producing photodissociation at the maximum rates. Photodissociation in the deeper layers is modified according to the amount of dust, as well as the amount of each individual molecule, in the layers above, both of which act to reduce the flux of dissociating photons. The attenuation of the UV field is determined independently for each photodissociation process. It is implicitly assumed that all photodissociation occurs through line absorption at specific frequencies; however, in general, the presence of dust mixed into the ice has a much larger effect on the attenuation of all photodissociation rates than do the individual absorption processes that are modeled implicitly here.

The attenuation of the radiation field associated with a specific photodissociation rate is calculated repetitively through each layer of ice, beginning in the outer surface layer and working inward. For layers that are thicker than one monolayer (i.e., all layers beneath the surface layer), the attenuation by each individual monolayer within that layer is considered. A mean photodissociation rate—averaged over all depths within that layer—is then calculated, which is applied (in the chemical calculations) equally to all material within that full layer.

To obtain the total fractional attenuation due to absorption in a single monolayer (${f}_{\mathrm{ML}}$) for a specific photodissociation process j, three fixed quantities are employed: the fractional attenuation per monolayer of dust for UV/visual-band photons, f_dust; the wavelength-dependent modifier of the dust absorption for the specific photodissociation process, c(j); and the fractional attenuation per monolayer due to absorption by the species being dissociated, f_PD(j). These are modified by the fractional abundance of dust (by volume) in that particular monolayer, X_dust, and by the fractional abundance, X(i), of the molecule being photodissociated, i, each abundance taking a value from 0 to 1, and whose precise values are time-dependent variables calculated in the code. Thus,

$$\begin{eqnarray}&&{f}_{\mathrm{ML}}(j)={f}_{\mathrm{PD}}(j)\,X(i)+{f}_{\mathrm{dust}}\,{X}_{\mathrm{dust}}\,c(j),\end{eqnarray} \tag{ 8 }$$

where f_ML(j) is the total fractional loss in intensity due to transmission through an individual monolayer, and $(1-{f}_{\mathrm{ML}}(j))$ is thus the fraction of the intensity transmitted by that monolayer. For an ice layer of total thickness N in monolayers, the fraction of intensity transmitted after passage through the full layer is therefore ${(1-{f}_{\mathrm{ML}}(j))}^{N}$. The mean rate of photodissociation applied to material within a particular layer is a factor ${{\rm{\Sigma }}}_{n=1}^{N}{(1-{f}_{\mathrm{ML}}(j))}^{n}/N$ multiplied by the total rate after attenuation by any layers above the layer in question. The value of f_PD(j) is assumed for all species to be 0.007, following the determination by Andersson & van Dishoeck (2008) for the attenuation of the photodissociation rate of water, per water monolayer.

For a single, canonically sized grain suspended in ice, a geometric cross section of ∼3.14 × 10⁻¹⁰ cm² is presented. To calculate (as required by this model) the average attenuation of the UV/visual-band photon flux per monolayer of dust, f_dust, an evaluation of the number of monolayers of grain material provided by the entire grain in the given cross section is required. This is easily calculated as the height of the canonical grain in monolayers, ≃622. However, a solid, spherical grain does not occupy all of the space in the cylinder traced out by its cross section; the sphere of the dust grain occupies exactly 2/3 of this space, meaning that the average fractional abundance of dust over each of the 622 monolayers of its height is X_dust = 2/3.

In this model, the attenuation per monolayer is assigned a value f_dust = 0.0055, which results in the transmitted UV/visual-band intensity being reduced by a factor ${(1-{f}_{\mathrm{dust}}\times 2/3)}^{622}\simeq 0.1$ over a full grain height (ignoring any contribution from f_PD(j)), corresponding to an absorption efficiency of ∼0.9 over the whole grain. The latter value is approximately equal to the mean of the typical absorption coefficients, Q_abs, for carbonaceous and silicaceous interstellar grains of radius 0.1 μm at a wavelength of 5500 Å, which take values typically around 0.1 and 1.7, respectively (e.g., Krügel 2008). Note that the simple one-dimensional treatment presented here considers only absorption and not scattering, while the close proximity between dust grains in the comet ice will also influence the optical transmission. The assumed value of f_dust = 0.0055, which is applied to the photodissociation rates via Equation (8), is thus quite imprecise.

The adoption of a single grain size in the above calculations is itself a substantial simplification; choosing an appropriate representative size is therefore important. In the interests of consistency with interstellar treatments, the adoption of a grain radius of 0.1 μm, as above, makes sense for this initial chemical modeling study.

The dust in observed comets in fact displays a distribution of sizes, from submicron to macroscopic scales; as summarized by Agarwal et al. (2007), much of it is composed of aggregates of submicron-sized subunits. The subunits are understood to have mass densities of a few g cm⁻³ (consistent with the value adopted in Section 2.3), while the aggregates would be porous, with densities in the range of 0.1–1 g cm⁻³. Agarwal et al. analyzed the dust environment of comet 67P, based on a selection of remote observational data. Their models suggest that the total dust mass is concentrated in grains of size 10 μm or larger, with one model showing a distinct peak around this value. Recent work by Mannel et al. (2019) presents imaging and analysis of dust particles obtained directly from 67P, studied using Rosetta's MIDAS Atomic Force Microscope. This work reinforces the idea that the dust is composed of aggregates of small subunits. Those authors find a mean size of ∼0.1 μm for the subunits. Experimental work by Price et al. (2010), which aimed to reproduce the impacts on detectors in the NASA Stardust spacecraft of small dust particles from comet 81P/Wild 2, also derived comparable sizes.

The fact that the larger grains observed in real comets appear to be low-density, loose, and/or porous aggregates of smaller subunits of 0.1 μm therefore seems consistent with the simple UV penetration model based around absorption by grains exclusively of that size. The adoption of a more comprehensive grain-size distribution in the chemical models, with implications for the UV treatment and the packing density of grains, is left for future studies.

The above treatment assumes that all UV impinging on the comet is interstellar in origin, rather than solar. Although the relative contributions of interstellar and solar UV at a particular distance from the Sun will vary with wavelength, a comparison at a representative wavelength is instructive. Curdt et al. (2001) presented solar spectra in the 670–1609 Å range (i.e., ∼18.5–7.7 eV); the majority of the wavelength-integrated dissociation cross sections for molecules of interest tend to fall within this range (see, e.g., Heays et al. 2017). Interstellar UV is dominated by emission from stars hotter than the Sun; thus, in general, solar UV will be most competitive with interstellar UV at the longer-wavelength end of the range. The Curdt et al. data provide a (quiet-Sun) solar radiance of ∼0.6 W sr⁻¹ m⁻² Å⁻¹ at 1600 Å. At an assumed inner radius for the Oort Cloud of 1000 au, this corresponds to a field strength of ∼3.3 × 10³ photons s⁻¹ cm⁻² Å⁻¹. The interstellar UV field of Draine (1978) gives a field strength of ∼1.9 × 10⁵ photons s⁻¹ cm⁻² Å⁻¹ at this wavelength, roughly 60 times greater than the solar value. Thus, to a first approximation, solar UV may be ignored, especially if the comet is significantly more distant than 1000 au.

Interstellar UV photon flux drops sharply at energies above the ionization potential of atomic hydrogen (13.6 eV). It is therefore possible that there is some small contribution of solar UV to the photodissociation rates that is not otherwise captured in the interstellar rates. The UV output of the early Sun was also likely higher than at present (e.g., Zahnle & Walker 1982), perhaps sufficient to become competitive with interstellar values. A more detailed treatment of the combined interstellar and solar photodissociation rates is left for future work.

The penetration treatment used for cosmic-ray-induced dissociation is presented in Section 4, along with the results of the models using that treatment.

To complete the treatment of reactions within and upon the cometary ice, a further mechanism is included in this model beyond those typically considered in interstellar ice chemistry simulations. The calculated reaction rates described in Section 2.1 are purely diffusive; thus, reaction depends on the mobility of the reactants and does not consider their origins, other than the implicit assumption that the reactants are not immediately in contact at the start of an individual diffusive reaction process. However, the production of reactive species as the result of photodissociation in the ice may lead to the spontaneous appearance of a reactive atom or radical in close proximity to another that is already present in the ice, leading immediately to reaction without any mediating (thermal) diffusion process. This possibility is especially important in the case of very low ice temperatures (i.e., 5–20 K), in which the thermal diffusion of atoms and/or radicals may be prohibitively slow.

To incorporate into the model the immediate reaction of photodissociation fragments with other reactants present in the ice, the total rate of production for each chemical species, as a result of all photodissociation processes that can produce it, is calculated (which is done separately within each ice layer). This total production rate, R_prod(i), is then multiplied by the fraction of the ice composed of each other species with which the photofragment may react, providing a rate of immediate reaction between the two species. Thus, for some reaction j, occurring between chemical species A and B, the total rate (s⁻¹) at which it occurs as a direct and immediate result of the production of either reactant by photodissociation is given by

$$\begin{eqnarray}&&{R}_{\mathrm{reac}}(j)={\kappa }_{\mathrm{AB}}\ \left({R}_{\mathrm{prod}}(A)\ \displaystyle \frac{{N}_{m}(B)}{{N}_{{\rm{M}}-{\rm{D}}}}+{R}_{\mathrm{prod}}(B)\ \displaystyle \frac{{N}_{m}(A)}{{N}_{{\rm{M}}-{\rm{D}}}}\right),\end{eqnarray} \tag{ 9 }$$

where, again, R_prod represents exclusively the production rate caused by photodissociation of precursor species, with other quantities defined as in Section 2.1. A similar expression is used for the same process occurring in the upper surface layer.

As an example, a hydrogen atom may be produced by the photodissociation of a selection of molecules in the ice, including H₂O, CH₃OH, and NH₃. The sum of these production rates is multiplied by the fractional abundance of, for example, CH₃, to give the rate for immediate production of CH₄. A further term corresponding to photoproduced CH₃ reacting with existing H would complete the prescription given in Equation (9). Because the fractional abundance of CH₃ and all other species ranges between 0 and 1, the resultant reaction rate will in practice always be less than the overall production rate of photofragments. It should be noted that in this treatment, as with the calculation of bulk-ice reaction rates, the determination of fractional abundances within the layers beneath the surface ignores the presence of dust because the dust and ice are distinct phases that cannot be fully mixed.

This new treatment is considered for all existing reactions in the network for which one or both reactants may be produced via photodissociation of other species. For reactions involving activation energy barriers, the usual barrier treatment is used, i.e., the fastest of thermal or quantum tunneling through the activation barrier, in competition with thermal diffusion of the reactants. Thus, the photoproducts involved are implicitly assumed not to have excess thermal energy for overcoming activation barriers, nor to be in an electronically or vibrationally excited state that might otherwise alter the barrier.

This initial version of the cometary ice/dust model is intended to represent (stable) conditions in the Oort Cloud. The main physical parameters are the choice of ice temperature and UV field strength. Seven temperature values are used in the models presented here: 5, 10, 20, 30, 40, 50, and 60 K. While the two lowest values are more representative of the typical values in the Oort Cloud, the higher temperatures are chosen to test how sensitive the solid-phase chemistry is to elevated temperatures, which might be achieved periodically through stochastic astrophysical events such as supernovae or the close passage of hot stars. It is also assumed here that the same temperature permeates the entire comet ice at every layer. Propagation of temperature structure within comets has been modeled in detail by other authors (see, e.g., Prialnik et al. 2008; Guilbert-Lepoutre et al. 2016); the effects of temperature structure within the chemical models are left for later work.

Three UV field strengths are tested in the models, based around the standard interstellar Draine field, G₀: low (G = 0.1 G₀), medium (G = G₀), and high (G = 10 G₀).

One model is run for each value of temperature and UV field strength, beginning from the same set of chemical initial conditions. Each model is run for a period of 5 Gyr, approximately the age of the solar system. Further models are run for the 5, 10, and 20 K cases only, in which the action of cosmic-ray impacts is approximated, in tandem with the UV-related effects; see Section 4.

The CPU time required to run each model varies, dependent on the degree of interaction between different ice layers (which affects the stiffness of the differential equations); however, all are on the order of 1 week, using a single thread per model run.

Figures 1 and 2 show the 5 K model results at t = 10⁶ yr and 5 × 10⁹ yr, respectively. It may be seen that the most drastic chemical changes occur in the upper 1 μm of the ice, with the initial ice constituents dropping significantly over time, except in the very surface layer. In Figure 1, very similar behavior is seen between the low- and high-UV cases. Conversely, by the end of the model runs (Figure 2), the abundance of dust has reached its maximum value, to a depth of approximately 10 μm for the low-UV case, 100 μm for the medium-UV case, and 1 mm for the high-UV case. The increase in dust fraction produces a commensurate drop in the main ice components, which may be more easily distinguished from the photodestruction of these species as the dust enrichment extends to greater depths over time. The concentration of dust in the outer layers is driven by the loss of material through desorption from the surface layer. At 5 K, this is caused exclusively by photodesorption, with water loss making up the largest individual component. The total loss rate from the surface at the end of the 5 K low-UV model run is 6.3 × 10⁻⁵ ML yr⁻¹, of which water makes up ∼55%. Other such rates, along with the total amount of ice lost in each model and the depth to which the dust reaches peak concentration, are shown in Table 1.

The photodissociation of the initial ice components results in the production of new molecules from the radical/atomic fragments. The abundance behavior of these products is similar across the various UV flux models. The large abundance of H₂, which here is the dominant molecule in the ice to a depth of a fraction of a micron, is mainly the product of the photodissociation of water (H₂O + hν → OH + H) followed by immediate reaction between the newly dissociated H atom and a contiguous H atom already present in the ice. The influence of thermal bulk-ice diffusion on reaction rates is negligible at 5 K compared to this immediate reaction process, even for H atoms. A significant amount of H and O (around 10% and 8%, respectively, for the medium-UV model) is also present in the upper subsurface layers, and these abundances are mostly constant to a depth of ∼1 μm. In the surface layer itself, the abundances of hydrides that are initially present in the ice remain relatively large throughout the model runs, due to the lower barriers to thermal diffusion, allowing photofragments to find each other and react, with hydrogen dominating the diffusion and reaction process.

O₂ is very abundant in the upper layers, as seen in the middle panels of the figures. Its production is mostly the result of the liberation and immediate reaction of oxygen atoms via the dissociation of CO₂ and OH (produced mainly from water), although, once formed, O₂ may also be photodissociated. The abundance of O₂ produced through this immediate photodissociation/reaction mechanism is around 1% of the total dust/ice composition, although in the region where O₂ is abundant, water itself is strongly depleted, so that the O₂/H₂O ratio is approximately unity. H₂O₂, N₂, and a range of hydrocarbons are also produced through similar mechanisms, via the destruction of the initial ice components. Oxygen atoms produced by photodissociation may also react with O₂, producing ozone (O₃, not shown in figures).

At relatively early times (10⁶ yr), product molecules are seen to reach a peak at their greatest depth (around 1 μm), before strongly dropping off in the deeper layers. However, at the end time of the models, the product abundance profiles are generally flatter down to a depth of 1 μm. The deep-peaking effect at early times is caused by the maintenance of high abundances of parent species deeper into the ice, where photodissociation rates are lower and therefore chemical timescales are longer. Later in the models, the abundances of parent species are reduced, and the chemistry reaches stability, balancing between photodestruction of parent and product species and the formation/reformation of both from the resulting radicals.

The right panels show the abundances of a selection of COMs, either detected in cometary environments, observed in the gas phase toward star-forming sources, or considered in interstellar chemical models. These species show the same early-peaking effect as the simpler products. The majority of complex species plotted show abundances of around 10⁻⁸–10⁻⁴ to a depth of 1 μm. The most abundant of these is formic acid (HCOOH), which is formed through the reaction HCO + OH → HCOOH. The HCO radical is produced by reactions of liberated H atoms with CO already present in the ice; again, no thermal diffusion of H is required in this case, due to the production of H atoms in the presence of the abundant CO molecule. Photodissociation of water then produces OH that reacts with the HCO. The large abundances of both CO and H₂O therefore result in significant formic acid production.

Many of the other COMs plotted in the figures are products of CH₃O or CH₂OH radicals. In interstellar ices, these radicals are usually thought to be formed either through H atom addition to formaldehyde (H₂CO) on dust grain surfaces or via the abstraction of a hydrogen atom from methanol (CH₃OH) on grains. In the 5 K comet models, the reactions CH₃ + O → CH₃O and CH₂ + OH → CH₂OH are the main formation routes. The relatively low abundances of H₂CO and CH₃OH (as compared with CO), combined with the barriers involved in either H addition or abstraction, make the interstellar routes less favored at 5 K. Reactions of the CH₃O and CH₂OH radicals with HCO tend to be somewhat more efficient than reactions with CH₃, due to the rapid H + CO formation routes for HCO.

Abundances of complex organics are similar between different UV field models, especially toward the end of the model runs. This is because, under the conditions at 5 K in which diffusion between layers is generally minimal, the chemistry of those molecules reaches a balance between formation and destruction that is dependent almost exclusively on photodissociation (in the upper layers where COMs are more abundant); photodissociation of smaller molecules produces COMs, while those products are also destroyed by photodissociation. Both formation and destruction therefore scale in the same proportion as the UV field.

Hydrocarbons and many of the complex organics show a strong tendency to increase in the outer surface layer, due to the greater ability of, e.g., methyl radicals to diffuse on the surface versus within the bulk ice. The gradual concentration of C₂ in the surface layer (due to desorption of more volatile species), following its formation in deeper layers, also results in higher abundances of its hydrogenation products such as C₂H₂ (see middle panels).

Glycine (NH₂CH₂COOH) is also plotted, although it barely reaches fractional abundances above 10⁻¹⁰. The reactions by which it may form in the current network all involve the addition of radicals and are described in more detail by Garrod (2013). Here, production occurs mainly through the reaction NH₂ + CH₂COOH → NH₂CH₂COOH. The amino radical is formed through direct dissociation of ammonia, as well as through the repetitive addition of H to atomic N (ultimately sourced from NH₃). The radical CH₂COOH is a dissociation product of acetic acid, whose formation requires the production of the CH₃CO radical. The latter is formed by the barrier-mediated reaction of CH₃ with CO.

Results for the 10 K models are plotted in Figures 3 and 4. By the end of the model runs, the concentration of dust is similar to that of the 5 K models, as most of the molecular material is still not thermally mobile within the bulk ice. Thus, water photodesorption still dominates the loss of surface material at rates determined purely by the UV field strength. However, aside from the effects of the concentration of dust in the upper layers, the abundance of water now remains close to its maximum value throughout the models. This is due to the thermal diffusion of H atoms in the bulk ice, which react with O and OH to re-form water. The mobility of atomic H means that an individual H atom has a short lifetime in the ice, disfavoring the reaction H + H → H₂, which depends on high fractional abundances of H. Hence, the H₂ abundance in the bulk ice drops to low values at temperatures of 10 K and above. The H₂ abundance produced increases with the UV field strength, and in the high-UV plot it can be seen that at 10 K this molecule is mobile enough in the ice to have thermally diffused into the deepest ice layers by the end of the model.

At 10 K, O₂ production in the ice is somewhat lessened, versus the 5 K models, due to H atoms reacting more frequently with atomic O and reducing its abundance in the ice. Production of O₂ still requires immediate reaction between photodissociation products. Hydrocarbon production is more pronounced at 10 K, with molecules such as C₂H₂, C₂H₆, and C₃H₈ reaching abundances in the upper micron that are 2–3 orders of magnitude greater at the end time of the model, versus the 5 K results. C₂ and C₃ also show more modest increases.

While the abundances of product molecules are still restricted to layers in the upper 1 μm, by the end of the 10 K model runs these species have been able to diffuse gradually to somewhat greater depths, or their radical or atomic precursors have done so, allowing the product species to be formed in the deeper layers to a small extent.

At 10 K, the abundances of complex molecules are noticeably higher in the upper layers, versus 5 K, with the fractional abundances of many species reaching around 10⁻⁴ with respect to total ice and dust content. The survival of greater abundances of formaldehyde and methanol in the upper layers of the ice allows production of CH₃O/CH₂OH through addition or abstraction of H to be maintained, resulting in more COM production. Glycine is also produced in greater abundance—again through NH₂ addition to CH₂COOH, although the latter radical is formed more easily in the 10 K models through the abstraction of a hydrogen atom from acetic acid, which itself is formed as a result of H abstraction from acetaldehyde followed by OH addition. Hydroxylamine, NH₂OH, suggested by some to be a possible precursor to glycine, is strongly produced through the addition of NH₂ and OH radicals, as well as by repetitive hydrogenation of NO by H atoms. The final abundance of NH₂OH continues to be substantial for all model temperatures, up to some limiting depth.

At 20 K (see Figures 5 and 6), the effects of thermal mobility within the bulk ice are more noticeable. The region of maximum dust concentration extends much deeper and more uniformly for all UV field strengths, with high dust abundances to depths of 1–10 mm at the end of the models. Thermal diffusion of CO within the ice becomes sufficiently rapid to allow it gradually to reach the surface layer, where it may be photodesorbed; thermal desorption of CO is still less rapid than the corresponding photoprocess by several orders of magnitude. However, for H₂, which is also able to reach the surface via diffusion, thermal desorption is indeed the dominant rate, and the production of H₂ within the bulk ice is still significant enough at 20 K for H₂ desorption to provide the greatest rate of loss of material from the surface (see Table 1).

Product molecules, both simple and complex, achieve relatively low abundances at 20 K, as compared with lower temperatures, because these molecules spread to much greater depths either through direct diffusion of their own or by bulk-ice swapping with highly diffusive species. O₂ and N₂ in particular are able to maintain modest abundances throughout the ice by the end of the high-UV run. H₂O₂ and HCN are also abundant to depths of around 1 m.

Ethane (C₂H₆) becomes the dominant hydrocarbon in the ice at these temperatures, with a fractional abundance around 10⁻⁴ near to the surface layer; this value is nevertheless somewhat less than the ∼10⁻³–10⁻² detected in comets. The C₂H₆ abundance appears fairly insensitive to the UV field strength. The abundances of complex molecules like methyl formate (CH₃OCHO) are significantly lower than at 10 K, due to their spread to greater depths. The production of complex molecules is still occurring mainly in the upper 1 μm of the ice, with the products then gradually diffusing to deeper layers. However, diffusion rates are high enough that production of new molecules through thermal diffusion of radicals in the ice has begun to dominate the chemistry (rather than immediate reaction with contiguous species, following photodissociation). The barrier against bulk diffusion for the CH₃ radical (823 K) in these models is somewhat lower than that for HCO (1120 K); the effect of this can be seen in the abundance ratios of CH₃- versus HCO-related product molecules such as CH₃OCHO and CH₃OCH₃, the latter of which attains abundances several times higher than the former. The achievement of these ratios is directly attributable to the changeover to a thermal-diffusion-dominated bulk-ice chemistry. For such reasons, glycine abundance falls off at 20 K, as the radicals involved in its formation are not very mobile at these temperatures and do not become mobile below around 50–60 K.

It may be noted also that at 20 K, in the high-UV case, the production of CH₃O and CH₂OH in the upper 1 μm is dominated by direct photodissociation of methanol, rather than the addition/abstraction of an H atom to/from H₂CO/CH₃OH. Production of HCO is also dominated by photodissociation of H₂CO, rather than the addition of H to CO, which delivers negligible amounts of HCO, due to the high thermal diffusion rate of atomic H competing with quantum tunneling through the activation energy barrier of the H + CO → HCO reaction.

The 30 K models (Figures 7 and 8) achieve strong dust enrichment to depths greater than 1 m by the end of the runs and achieve maximum dust abundances to depths of 1 cm within 10⁶ yr. Thermal desorption of small atoms and molecules such as H, H₂, O, and O₂ dominate the losses from the surface, but the precise values are strongly dependent on the UV field strength of the model. Table 1 shows that at a sustained temperature of 30 K, after 5 Gyr the comet has lost ∼13 m of ice.

The strong degree of diffusion at these temperatures again results in low abundances of many product species, especially at early times. Some molecules, such as acetylene (C₂H₂), tend to survive only weakly at greater depths and late times, due to destruction by H atoms, which are highly mobile; acetylene is ultimately hydrogenated all the way to ethane (C₂H₆). On the other hand, hydrogen peroxide (H₂O₂) is not very diffusive, due to similar binding and diffusion characteristics to water itself. This means that H₂O₂ has its highest abundances in the upper 1 μm where UV irradiation is strong. The more volatile initial ice constituents, especially CO and CH₄, desorb strongly at 30 K (although they do not dominate the total desorption rate), while CO may also react rapidly with OH to form CO₂ + H; as a result, CO₂ tends toward higher abundances in the shallower layers, over long timescales.

Although the local abundances of many molecules appear somewhat low in the 30 K models, their aggregated abundances through the ice are substantial and are significantly higher for most species than in the 20 K models, as they reach to deeper layers of the cometary material.

At 40 and 50 K, the loss of surface material is consistently dominated by H, H₂, and O desorption (Table 1). At 60 K, CO, CO₂, and H₂CO dominate. The total loss rates from the surface (corresponding to the final values) increase with temperature and UV field strength. However, the overall loss of material over the course of each model is not substantially different for temperatures between 30 and 60 K.

With increasing temperature, the abundances of CO, CH₄ and H₂CO become smaller due to diffusion, reaction, and surface desorption, so that by a time of 1 Myr in the 40 K models, CO and CH₄ abundances are negligible compared with the less volatile initial constituents. At a temperature of 60 K, all three of these species decrease substantially in abundance throughout the ice within the first 1000 yr of model time.

The majority of complex organics are formed most efficiently around 40 K, although ethylene glycol, (CH₂OH)₂, is a notable exception; it is formed through the addition of CH₂OH radicals, which have a bulk diffusion barrier in this model of 3556 K (see Garrod 2013). In the higher-temperature models, the smaller radicals such as CH₃ and HCO become very diffusive and volatile, decreasing their abundances, as they are able to desorb from the nearby surface layer in preference to diffusing into the deeper layers. At 60 K, nevertheless, complex molecule abundances extend through every layer. Higher UV fields at high temperatures still have a noticeable effect on the ability to produce such molecules. At 60 K, glycine fractional abundance is around 10⁻⁸ for the high-UV case, while it is orders of magnitude lower for the low UV field strength model.

It should be noted again that these high-temperature models, with fixed temperature, are run to a final time of 5 Gyr, for consistency with the other models. However, it is unlikely that stochastic heating events could produce these temperatures for such long timescales (see Section 5).

Figures 12 and 13 show results from the cosmic-ray models at times 10⁶ yr and 5 × 10⁹ yr, respectively. The top row of panels correspond to a temperature of 5 K, the middle panels to 10 K, and the bottom panels to 20 K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The low temperatures in the cosmic-ray models mean that the rates of thermal diffusion are fairly small, especially at 5 K; thus, reaction rates are dominated by the immediate reaction of radicals following their production via the dissociation of parent species (whether through UV or cosmic-ray processing). In the upper micron of material, the chemical behavior is similar to the models with only UV processing included; for many molecules a step-like feature is seen, as the UV-induced dissociation drops off. However, at greater depths, various product species take on non-negligible abundances, even after as little as 1 Myr (Figure 12), due to the weaker—but more penetrating—cosmic-ray-induced dissociation. At 5 K, some complex organics reach fractional abundances of around 10⁻¹⁰ to depths on the order of 10 m at this time. In the 10 K model, complex organics have somewhat higher abundances at those depths, while their abundances become as high as 10⁻⁵ to depths of 1 mm. Such behavior is not seen in the UV-only models, either in the 5 K or the 10 K case. The enhancements down to ∼1 mm, and the weaker enhancement to ∼10 m, are the result of the cosmic-ray processing alone. At 20 K, the behavior is more comparable with the UV-only models to depths of around 1 mm, but the enhancements seen to depths >10 m are unique to the cosmic-ray model. At 10⁶ yr, this model produces COM fractional abundances on the order of 10⁻⁸ to depths of around 10 m.

Abundances of H₂, O₂, and H₂O₂ are also noticeably enhanced at a time of 10⁶ yr in the 5 and 10 K models with cosmic-ray processing. In the 10 and 20 K models, the abundances of hydrocarbons are enhanced at greater depths also, although the hydrocarbon abundances at depths less than 1 μm are suppressed versus the UV-only models.

At a time of 5 Gyr in the 5 K models, H₂O₂ is very abundant to depths of around 10 m or more, while O₂ maintains a fractional abundance greater than 10⁻³ to a depth of several meters and greater than 10⁻⁴ to a depth of 10 m. H₂ is abundant throughout the ice and maintains an abundance greater than 10⁻⁴ to several tens of meters. In the 10 K models, O₂ and H₂O₂ are somewhat lower in abundance but retain the same general behavior, while H₂ becomes a very minor component of the ice, versus the 5 K model. All three species achieve only very small abundances at 20 K. Meanwhile, HCN and the saturated hydrocarbons C₂H₆ and C₃H₈ retain a consistent abundance around 10⁻⁴–10⁻⁵ at all three temperatures at the end of the cosmic-ray models. N₂ maintains an abundance of at least 10⁻⁵ to a depth of a few meters in the 5 and 10 K models but is lower in the 20 K model.

The production of O₂ at large depths into the ice, where only cosmic-ray processing is important, occurs through the dissociation of O₂H (which is formed via the reaction between O and OH) and by direct addition of oxygen atoms. OH radicals are produced by the dissociation primarily of water, but also H₂O₂, while atomic oxygen is formed by dissociation of CO₂, OH, CO, and several other species, including O₂ itself.

At 5 Gyr, many complex molecules are quite abundant throughout the 5 and 10 K ices, maintaining fractional abundances as high as 10⁻⁸ down to around 10 m. Values as high as 10⁻⁴ are achieved to these depths for many species. Glycine attains a maximum abundance around 10⁻⁸, to a depth of a few meters. Local abundances of the most complex organics become significantly lower in the 20 K model, caused by the gradual thermal mixing of molecules into the deepest layers over the full 5 Gyr timescale; however, as a result, substantial abundances of these molecules are maintained to the very greatest depths sampled in the models. By the end time of the cosmic-ray models, the larger hydrocarbon species and many other complex organics tend to achieve lower abundances in the upper ice layers than they do at depths greater than 1 μm; the combined UV- and cosmic-ray-induced destruction of these species is rapid in the upper layers.

Tables 2–5 show aggregated abundances to a depth of 15 m, presented as a fraction of water in that same depth, for the 60 K high-UV (UV-only) model and for each of the cosmic-ray-processing models at various times. The last two columns in the tables show abundances of various molecules for comet Hale–Bopp (C/1995 O1), determined via remote observations, and for Churyumov–Gerasimenko (67P), determined by the COSAC and ROSINA mass spectrometers.

The 15 m depth over which the aggregated abundances are calculated is appropriate to demonstrate the overall effect on abundances of, in particular, cosmic-ray processing (Tables 3–5). However, beyond this depth, the resolution of the models is insufficient to provide accurate aggregated abundances, as the depth at which cosmic-ray processing falls off (and thus molecular abundances fall sharply) lies somewhere within the next deepest layer, which has a thickness of ∼30 m, making the integration quite imprecise. Future models will adopt higher depth resolution around the threshold region.

In order to be able to compare directly between the abundances produced by the models and those obtained from observations of active-phase comets, consideration must be given as to how much material has been lost by those bodies since their entry into the inner solar system. The loss of surface material during a comet's active phase would mean that at least some fraction of the high molecular abundances residing in the upper layers (based on the present models) would likely be lost. Typical dynamical lifetimes for comets are a few hundred thousand years, so in the case of JFCs, for example, there may have been thousands of orbits (stripping away thousands of meters of the surface) prior to observation. The present models are therefore most applicable to relatively dynamically new comets (those that have only entered the inner solar system a few tens of times at most), or dynamically new comets (those that are entering the inner solar system for the first time).

It is useful, then, to consider the depth to which the modeled abundances must extend in order to be directly comparable with the Rosetta values, any further chemical processing during the active phase notwithstanding. The erosion rate of comet 67P has been calculated to be approximately 1.0 m per orbit (±50%) by Bertaux (2015), while its orbital period is determined to be 6.44 yr, providing an average rate of 0.16 m yr⁻¹. Assuming that no significant erosion occurred before its 1840 encounter with Jupiter that brought it into a closer orbit (although its earlier orbital activity is uncertain; Królikowska 2003), this corresponds to ∼13–41 m of material lost during the comet's active phase. For the solid-phase chemistry simulated in the model to be important to the recently observed chemical abundances of 67P, such chemical processing should therefore extend to such a depth, or greater. The 15 m aggregated abundances calculated from the models thus satisfy at least the lower limit of this range. Substantial abundances of molecules of interest could extend beyond 15 m; however, as mentioned above, the present models lack the depth resolution to determine this.

The other comet with which the modeled abundances are compared is Hale–Bopp; prior to 1997, it is estimated to have made its last perihelion in 2215 BC, which has also been speculated to have been its first (Marsden 1997). Based on submillimeter measurements of particle mass loss during its most recent solar approach, Jewitt & Matthews (1999) calculated a mean value of around 10 m of material lost per orbit. This would suggest, assuming only one previous perihelion with similar mass-loss characteristics, that only around 10 m of material should have been lost prior to the 1997 perihelion, and that at least 5 m of native material within the upper 15 m or so would have survived to be released during that well-observed event. Direct comparison of the aggregated values presented in Tables 2–5 with the observed molecular abundances determined for Hale–Bopp therefore seems justified. It should be noted, however, that the penetration depths of cosmic rays used in the models correspond to the density of material in 67P in particular, which could differ from that of Hale–Bopp.

Although some UV-only models show significant abundance enhancements of COMs to a depth of around 10 m or more, for most temperatures they cannot reproduce values detected in comets 67P or Hale–Bopp (or not to appropriate depth; see below) and are thus not presented in the tables. The 60 K UV-only model is perhaps the most successful such model at reproducing the observed values and is thus included to demonstrate the maximum likely contribution from UV surface processing. Unlike the other UV-only models, the 60 K models maintain substantial abundances of COMs throughout all layers as a result of efficient diffusion. At this temperature, the aggregated abundances begin to approach the measured values for 67P or Hale–Bopp in some cases but do not actually reach them, aside from formamide (NH₂CHO), with formic acid (HCOOH) and methyl formate (CH₃OCHO) falling only a little short. In general, it appears that none of the UV-only models are capable of reproducing either set of observed values at appropriate depths. However, for depths less than the 15 m used in the tables, fractional abundances of complex organics produced over the 5 Gyr lifetime are comparable to those observed in the gas phase in star-forming cores (with the presumption that many of these molecules may have formed in dust grain ice mantles). This in itself would suggest that UV processing of cometary ices at moderately high temperatures could allow a comparable degree of COM production in post-formation cometary material as might be expected to be incorporated from the interstellar or proto-nebular stage (however, see also Section 5.3).

In contrast to the UV-only models, those including cosmic-ray processing of the ice sustain COM production to great depths, regardless of temperature. Crucially, the comparison of aggregated abundances with measured values for Hale–Bopp and 67P is much better for these models. Indeed, in some cases, certain abundances are a little higher than the observed values by the end of the models, in particular formic acid (HCOOH). However, it is only in the 5 and 10 K models of cosmic-ray processing that a good match is provided for both the simple and complex species; diffusion into the deepest layer of the ice at 20 K diminishes abundances throughout. Weaker bulk diffusion effects, as discussed in Section 5.3, would raise the temperature at which the models could reproduce observational abundances.

At 5 and 10 K, COMs including formic acid, acetaldehyde (CH₃CHO), methyl formate (CH₃OCHO), and formamide (NH₂CHO) achieve aggregated abundances that are in reasonably good agreement with abundances observed for Hale–Bopp, to within a factor of 2–3 above or below the observations, while the abundance of dimethyl ether (CH₃OCH₃) is also in line with its observed upper limit. Ethylene glycol, (CH₂OH)₂, is underproduced with respect to Hale–Bopp and 67P abundances, although its only formation mechanism included in the network—the addition of two CH₂OH radicals—also does a poor job in interstellar models. An alternative production mechanism, such as the barrier-mediated hydrogenation of glycolaldehyde by H atoms, may be required to explain its abundance. Glycolaldehyde (CH₂(OH)CHO), uniquely among the COMs, comes close to the detected abundance in 67P. Other complex organics broadly fall around 1–2 orders of magnitude below their detected abundances in 67P.

The success in reproducing the abundances of COMs in Hale–Bopp in particular may give credence to the idea that its chemistry is indicative of a young and relatively uneroded surface, whose outer layers from the cold-storage phase survived, allowing their complex molecular content to be observed in its coma.

Interestingly, the 5 K model produces aggregate abundances of O₂ that are comparable to the observed value for 67P. Again, the agreement is not so good at 10 K and very poor at 20 K. However, the results would indicate that cosmic-ray processing of water over several billion years may be sufficient to produce the observed fractional abundances of this molecule in comet 67P; fractional abundances of at least 10⁻³ with respect to water are maintained to a depth of around 10 m, although at the depth resolution used in this model the precise depth is rather uncertain. Resolution on the order of 1 m or so at depths of 1 m and beyond may be necessary to investigate fully the cosmic-ray processing of cometary ices.

The production of significant amounts of O₂ to large depths in the cosmic-ray processing models is dependent on several factors, including the high porosity/low density of the cometary material assumed here (based on estimates for comet 67P), which determines how deeply the cosmic rays may penetrate. The precise origin of this porosity, and whether it was present at the earliest times in the evolution of 67P, or whether it was itself the result of some kind of processing, will determine whether the density estimate is appropriate for use in these models.

It should also be noted, however, that the rate and/or efficiency of water dissociation (and of other molecules) at various depths by galactic cosmic rays is not well constrained and may be a more important parameter. The estimate of the water dissociation rate used in the present models (discussed in Section 4) is around 50 times lower than the value implied by the calculations of SH18, for the case where the cosmic-ray flux is unchanged by penetration into the material (i.e., the zero-depth rate).

Bieler et al. (2015) measured an average O₂ abundance with respect to water in comet 67P of 3.80% ± 0.85%, with individual measurements ranging from 1% to 10%. If O₂ can be formed and maintained during the cold-storage phase at an abundance of around 1% with respect to water, to a depth on the order of 10 m, this may help to explain the mysteriously high abundance of this molecule in 67P. If O₂ is indeed formed this way, within the comet itself during cold storage, its production at large depths would be dependent on cosmic rays of energies >10 GeV; such high-energy particles are likely to be of similar flux both in the Oort Cloud and within the heliosphere; thus, all comet families would be expected to be similarly affected (prior to any surface ablation). The 5 K models presented here produce appropriate O₂ abundances to a depth of around 10 m, which is shy of the depth required if 67P has lost 27 m of ice in its recent history, but close to the lower limit on this estimate (13 m). However, an increase in the estimated cosmic-ray-induced dissociation rate of water by even a factor of a few would push the efficient production of O₂ to a depth comfortably consistent with the required range. An increase by a factor of 50, to be consistent with SH18, would be expected to give values well in agreement with observations. However, the relatively low depth resolution of the present models in this depth range means that the precise depth of O₂ production is poorly constrained.

The failure of the cosmic-ray-processing model at 20 K to provide O₂ in significant abundance in the same way as the 5 K model is largely due to the thermal mobility of the O₂ molecule, while at temperatures much greater than 20 K, O₂ might be expected to leave the ice entirely. Mousis et al. (2016) suggest that the retention of O₂ to high temperatures in 67P itself is due to its clathration with surrounding H₂O molecules when it is formed via radiolysis, although those authors suggest that the formation process occurs in dust grain ice mantles in the proto-solar nebula, prior to incorporation into a comet. Nevertheless, if clathration is indeed the mechanism that allows O₂ to remain trapped up to the high temperatures of an active-phase comet, high-energy cosmic-ray processing of the comet itself would also appear to be a plausible mechanism for its formation, even for a JFC such as 67P, based on the present models. It is also likely that the deepest layers of a comet would retain a very low temperature (e.g., Guilbert-Lepoutre et al. 2016), even during the active phase, which would minimize thermal mobility over long periods. As discussed below, the back-diffusion effect may also reduce somewhat the net mixing rate, although this would be a relatively small effect if O₂ retains a large abundance. The formation of a dust-rich layer on the comet surface could also help to retain any mobile O₂.

Bieler et al. (2015) also measured hydrogen peroxide abundance in 67P, obtaining a value H₂O₂/O₂ ≃ 6 × 10⁻⁴. The comet chemical models actually produce H₂O₂ more abundantly than O₂ in the deeper ice layers. Due to their related formation mechanisms, and the possibility for O₂ to be directly hydrogenated to H₂O₂, it is possible that the production of hydrogen peroxide is too efficient in the models. O₂ might therefore be expected to attain higher abundances than the models suggest, at the expense of H₂O₂. The depth attained by O₂ in the ice would also be greater as a result, more similar to that achieved by H₂O₂ in the models. The sum of these two molecules down to the aggregated depth of 15 m is certainly sufficient to agree adequately with the observed O₂ abundance, based on the 5 K model. The 10 K model also is not far from the observed value.

Rubin et al. (2015) detected N₂ emanating from 67P, although, since the ratio with H₂O is uncertain, its abundance is not included for comparison in the tables in this paper. Rubin et al. measured an N₂:CO ratio of $5.7(\pm 0.66)\times {10}^{-3}$. The cosmic-ray-processing models presented here produce aggregated abundance ratios for N₂:CO of 2.1 × 10⁻⁵ and 1.1 × 10⁻⁴ for the 5 and 10 K models, respectively, while the ratio for the 20 K model falls short by a larger margin. All of the cosmic-ray-processing model values thus fall short of the measurements; also, the depth reached by N₂ in the 5 and 10 K models is a little less than is achieved for O₂. The precise values of the modeled ratios for these, as for all other molecules, will depend also on the initial abundances assumed for them and/or their parent molecules, i.e., water, ammonia, CO, and CO₂. It is also possible that N₂ and other small molecules may be formed in the gas phase, so the solid-phase model abundances may not bear a direct comparison.

The achievement of significant COM abundances to depths of tens of meters in the models at 60 K depends on the long-range mobility of either the complex organics themselves or their radical or atomic precursors, over periods on the order of a billion years. The models assume that bulk diffusion is reasonably rapid, and there is some laboratory evidence for such behavior for certain molecules mixed with water (e.g., CO₂; Fayolle et al. 2011). Laboratory experiments of this kind typically involve large quantities of the diffusing species in relatively thin ices (on the order of 1 μm, often much less), meaning that the effects of back-diffusion, i.e., random walk, on the diffusion timescale are small. However, the extension to ice many meters in thickness, with diffusers of very low fractional abundance, could well produce different behavior. The interstellar gas-grain chemistry model of Garrod et al. (2017) applied a back-diffusion treatment for transfer between ice layers that was based on the parameterization of the results from a simple three-dimensional Monte Carlo diffusion model. Here, the number of bulk diffusion events, N_move, required for a single diffuser (placed at an arbitrary starting position) to reach the surface of an ice of finite thickness, N_thick monolayers, followed the expression ${N}_{\mathrm{move}}=2{({N}_{\mathrm{thick}}+1/2)}^{2}$. For transfer of low-abundance species between the deepest (i.e., thickest) layers in the comet model, this would suggest diffusion rates many orders of magnitude lower than assumed in the present models. However, in the other extreme, diffusers that make up a fraction of the total ice on the order of unity will show a negligible net back-diffusion effect, regardless of ice thickness. Back-diffusion is therefore highly sensitive to the abundance of the diffuser and requires a more careful calculation, as per Garrod et al. (2017), when dealing with very thick ices. For this reason, the Garrod et al. expression for low-abundance diffusers (given above) was not used in the comet models. Furthermore, its blanket application in fact would violate the necessary equality of diffusion rates across an arbitrary boundary for high-abundance species, due to the mismatch in layer thicknesses employed in these models.

It is also likely that the degree (or even the possibility) of bulk diffusion will be dependent on the specific composition of the ice. The degree of diffusion allowed to nonwater molecules in a highly water-dominated ice may be negligible, due to the buildup of crystalline structure. Ices composed of more balanced mixtures of water, CO, CO₂, and other commonly detected solid-phase molecules may make possible the degree of COM diffusion required by the UV-only models to produce significant COM abundances at large depths into the ice. The degree of mixing exhibited in these models for various chemical species—complex or otherwise—into the deeper layers in particular may therefore represent something of an upper limit.

While it has been assumed here that the diffusion within the ices occurs through bulk diffusion, cometary ices are expected to be somewhat porous; indeed, Jorda et al. (2016) determined the density of comet 67P to be ∼0.53 g cm⁻³, suggesting significant porosity. This value was used to determine the cosmic-ray penetration in the models presented here, based on a pure water ice. A porous ice structure would provide an internal surface that could allow relatively rapid diffusion compared with bulk processes. This would have the local effect of increasing the rates of reactions that are moderated by diffusion, as well as making mixing between layers more rapid, assuming that the porous structures are well connected and not enclosed. Again, surface random walk would likely have a significant influence on the mixing rates of low-abundance species through such a mechanism. Willis & Garrod (2017) studied the influence of two-dimensional random walk on the reaction rates of diffusive species on interstellar dust grain surfaces, taking into consideration the abundance of the diffuser with respect to the available number of surface sites. While the effect was small in the case of discrete microscopic dust grains, on the order of a factor of a few, when extended to comet scales, the random walk effect could arguably slow down diffusion by a factor of around 10–100 for a molecule of fractional abundance 10⁻⁸ (depending on the available surface area and its geometry).

The chemical and/or cosmic-ray processing of the ice could plausibly result in the gradual production of porous structures over time, even starting from a relatively nonporous initial ice. More detailed models that explicitly include pores/cavities, as well as chemistry upon their internal surfaces, would be necessary to confirm such effects.

Cosmic-ray processing also produces COMs to large depths in these models. Glycine is not formed in great abundance (unlike other COMs); while this is technically consistent with the detected values for 67P, which range from 0 to 0.0025 with respect to water, the models do not reach anywhere close to the maximum detected value. Altwegg et al. (2016) suggest that the variation they detect is caused by changing conditions on the comet rather than uncertainty in the values. Furthermore, their results apparently suggest that the detected glycine is associated with its release from liberated dust, rather than ice, which is itself heated within the coma. Thus, the ratio of glycine with respect to water may not be a meaningful reflection of the overall ice composition. Altwegg et al. (2016) also compare the abundance of glycine with that of methylamine (CH₃NH₂) and of ethylamine (C₂H₅NH₂), finding ratios of 1.0 ± 0.5 and 0.3 ± 0.2, respectively. All of these molecules are included in the present models, although ethylamine has only a sparse network in comparison with the other two.

In the cosmic-ray-processing models, glycine is at least an order of magnitude or so less abundant than methylamine at all times and depths. The production of glycine relies mostly on the addition reaction NH₂ + CH₂COOH → NH₂CH₂COOH, as mentioned in Section 3.1, rather than on the reaction NH₂CH₂ + COOH → NH₂CH₂COOH, which might be expected to follow more closely the abundance behavior of methylamine. However, the methylamine itself is formed by the addition of NH₂ and CH₃ radicals in the ice, so methylamine production should nevertheless be expected to show similar trends. This is indeed the case, as seen in Figure 14, in which the abundances of glycine, methylamine, and ethylamine are plotted for the cosmic-ray-processing models. Although methylamine abundances are uniformly lower than those of glycine, the same general behavior is seen for both, especially after the full 5 Gyr time span of the 5 and 10 K models. The models thus show some consistency with the detected behavior. Of course, if some or all of either molecule was formed in the ices prior to incorporation into the comet, then this analysis carries less weight. It may be seen that the abundance of ethylamine tracks very closely with that of glycine at 5 and 10 K, in good agreement with detected behavior. However, ethylamine in the present network does not have a formation mechanism directly involving NH₂; thus, the precise agreement may simply be fortuitous, although, like methylamine, the general agreement in its abundance trend may be more reliable.

Zoom In Zoom Out Reset image size

Based on the models presented here, if the true glycine abundance in the nucleus of 67P is close to the detected upper limit with respect to water, then its production is likely not the result of UV or cosmic-ray processing during cold storage.

Setting aside that the present models do not consider a native composition at the start of the cold-storage phase that includes COMs, the broad observational implication of the UV and cosmic-ray processing investigated here is that new (and relatively dynamically young) comets in particular would show an enhanced complex molecular content versus older comets. Beyond a depth of perhaps several tens of meters, no effect is expected (although the influence of thermal diffusion is difficult to quantify, as discussed above). Molecules that are direct products of radicals formed from the most abundant ice components should be especially abundant in new comets, according to the models. These would include species already observed in comets (Hale–Bopp in particular) such as CH₃CHO, HCOOH, NH₂CHO, CH₃OCHO, and CH₂(OH)CHO, as well as related molecules including ethanol (C₂H₅OH) and dimethyl ether (CH₃OCH₃). Hydroxylamine (NH₂OH), although not yet detected in the interstellar medium, nor apparently in any comets, is predicted to have significant abundances also, as well as glyoxal (OHCCHO).

The observed abundances of nitrile species (e.g., HCN and CH₃CN) and the cyanopolyyne HC₃N are generally not well reproduced by the models. This may be a function of the simple starting ice composition. The inclusion of, say, HCN or HNCO in the initial ices may have a strong impact on the later products, although the more complex molecules may also have been present themselves in the comets without requiring further processing. Given the observed, substantial abundances of H₂S and some other sulfur-bearing species in comets, it might be expected that the large abundances of oxygen-bearing molecules like methyl formate could be accompanied by similar sulfuretted structures.

The various comparisons made here between the model abundances and the observed values must be tempered by the fact that the active phase itself may promote some degree of thermal chemistry within the ices, relying on the presence of unreacted radicals in the ice, that could become more mobile at elevated temperatures and thus more reactive. Other thermally activated reactions (including those between stable species) could be enhanced also. The production of yet more complex species than the typically observed "hot-core" type molecules, which includes the example of glycine as discussed above, may be more dependent on the elevated temperatures achieved during the active phase. Thus, the ideal comparison of models with observations would include the simulation of the time-dependent physical conditions associated with entry into the inner solar system.

The chemical compositions of volatiles ejected from the two largest fragments of comet Schwassmann–Wachmann 3 (73P) were observed by Dello Russo et al. (2007). A JFC with a likely origin in the Kuiper Belt, 73P split into multiple fragments during its 1995 apparition, exposing material from deep within the comet nucleus. The observations by Dello Russo et al., taken during the comet's 2006 return, show that the two fragments are of very similar composition, indicating a lack of chemical differentiation within the nucleus. It is unclear whether the depth scale of processing found in the chemical models presented here, of order 10 m, would be sufficient to differentiate material originating from the two fragments of 73P, which have sizes of hundreds of meters. The molecules observed by those authors are also fairly small compared with many of those investigated here. The modeled abundances of native ice species with respect to water, aggregated to a depth of 15 m as shown in Tables 2–5, vary the most for methanol (CH₃OH), which is depleted by a factor of around 4 at the end of the 5 and 10 K cosmic-ray-processing models. Methanol varies by a factor of around 1.3 between 73P-B and 73P-C. If any chemically processed material lay only in the upper 10% of the nuclear material in those fragments, with the remaining material being otherwise identical, then the model results would seem compatible with those observations.

Gibb et al. (2007), on the other hand, found that the Oort Cloud comet C/2001 A2 (LINEAR), which fragmented into two pieces, indeed showed evidence of nuclear heterogeneity, based on variability in the abundances of H₂CO and CH₄. These two molecules also are found to vary from their initial abundances in the 5 K cosmic-ray models, increasing their abundances by factors of 1.7 and 2.2, respectively. Observations of the dynamically new comet C/2012 S1 (ISON) (Dello Russo et al. 2016; DiSanti et al. 2016), which broke up upon its close approach to the Sun, showed strong variations in the abundances of some molecules (notably NH₃, H₂CO, HCN, and C₂H₂) at heliocentric distances less than 0.5 au, versus those outside this radius. This may also suggest heterogeneity in the nuclear ices, although the precise origins of the detected molecules are uncertain.

Comets C/2001 A2 and C/2012 S1 would each appear to be a better dynamical match for the assumed physical and chemical conditions of the models presented here than 73P, but none of those cometary molecular observations would appear to rule out strongly the suggested scenario of chemical processing in the Oort Cloud that is investigated in these simulations.